Hello again!

Well, 'today' turned into 'tomorrow' but here I am again with good news.
I managed to successfully implement my work into viewweight.

User can input an arbitrary amount of points (in g/cm3) to generate a
voronoi mesh and program will return the total mass of the material in kg.

Obviously, now maybe 20% of the code inside viewweight is useless or wrong
for the current implementation, but the main operation works (getting the
mass of an object assuming only one material, who's density distribution is
defined through the point mesh). Anything related to the density file right
now is unused (we might wanna use it to store the points?). Some of the
printf happening in the end is now wrong since we are not using
per-material densities nor are we using a density table.

Anyway, what changes should be carried out from now on still would need to
be discussed.

Code attached!!!

Mario Meissner.

On 14 October 2017 at 08:12, Mario Meissner <mr.rash....@gmail.com> wrote:

> Hi Sean,
>
> Thank you for the update! I will come back with more changes today.
>
> Cheers!
>
> On Oct 13, 2017 21:02, "Christopher Sean Morrison" <brl...@mac.com> wrote:
>
>> Hey Mario,
>>
>> That’s looking good.  You don’t have anything to worry about.  The past
>> month has been solid effort.
>>
>> I talked with Adam yesterday and gave him an update and ETA on the report.
>>
>> Cheers!
>> Sean
>>
>> On Oct 13, 2017, at 11:42 AM, Mario Meissner <mr.rash....@gmail.com>
>> wrote:
>>
>> Hello again!
>>
>> Attached code is cleaner and more thoroughly tested.
>> This is what I will implement into rtweight, which shouldn't be a too
>> difficult task.
>> I hope to be able to come back to you with rtweight changes tomorrow.
>>
>> Adam, the coordinator for SOCIS, contacted me again urging us to finish
>> the report for me as soon as possible.
>> I haven't heard much feedback about my progress in my last month of
>> coding, which left me a bit uncertain. I hope we reached a decent milestone
>> by now. We should get back to Adam as soon as possible with the report.
>>
>> Greetings,
>> Mario.
>>
>> On 11 October 2017 at 09:52, Mario Meissner <mr.rash....@gmail.com>
>> wrote:
>>
>>> Hello everyone.
>>>
>>> I received an email from the SOCIS organizer saying that the coding
>>> period is finishing really soon. However I haven't received any news from
>>> your side. What is gonna happen now? Can we (and will we if so) extend it
>>> any further? I would be able to continue, but in an intermittent basis.
>>> This week I will integrate my code into rtweight.
>>>
>>> Mario
>>>
>>> On Oct 9, 2017 11:05, "Mario Meissner" <mr.rash....@gmail.com> wrote:
>>>
>>>> Hi people!
>>>>
>>>> I proudly present my latest work, attached here. It features a clean
>>>> algorithm to compute the average density of a segment that goes through a
>>>> cloud of density regions arranged in a Voronoi Tesselation. For now
>>>> confirmed to work for one and two point cases, but I suppose only minor
>>>> bugs will hinder it from properly handling any arbitrary amount of points.
>>>>
>>>> Note that it still requires some refactoring and cleaning since there
>>>> are some remains of the old code laying around.
>>>>
>>>> The last step would be to integrate this into rtweight, which shouldn't
>>>> be too much of a problem.
>>>>
>>>> As always, feedback welcome.
>>>>
>>>> Thank you,
>>>> Mario.
>>>>
>>>>
>>>>
>>>> On 6 October 2017 at 19:41, Mario Meissner <mr.rash....@gmail.com>
>>>> wrote:
>>>>
>>>>> Hi Daniel,
>>>>> thank you! In the end the problem was that I didn't run cmake after
>>>>> svn up.
>>>>>
>>>>> Attached goes today's work. I coded the function 'region_boundary()'
>>>>> which is an upgrade to the now unused intersection(). It does the same
>>>>> thing but now also checks if the intersection point lays inside or outside
>>>>> the segment, AND if the intersection is a valid region boundary or not (by
>>>>> checking distances to all other points, something I was doing outside of
>>>>> the intersection call until now). All in all a much cleaner function that
>>>>> will allow stages 1 and 2 to be much simpler and possibly run in one 
>>>>> single
>>>>> loop.
>>>>>
>>>>> On Sunday I will probably finish implementing my new algorithm (which
>>>>> is now trivial since the hard work was this function).
>>>>>
>>>>> Note that current code is only meant to be a testing platform for the
>>>>> new function. On line 386 I included a snippet that calls my new function
>>>>> and prints its results, then I break out with a return. Feel free to try
>>>>> out some new intersections. Currently I tested with three points, first 
>>>>> one
>>>>> ignored, and the second and third ones acting as 'previous' and 'current'
>>>>> respectively.
>>>>>
>>>>> Cheers!
>>>>>
>>>>> On 6 October 2017 at 17:41, Daniel Roßberg <danielmrossb...@gmail.com>
>>>>>  wrote:
>>>>>
>>>>>> spectrum.c was recently removed.  So, update your checked out sources
>>>>>> and check your librt CMakeaLists.txt that it matches with the one on the
>>>>>> head of the reposity (you've probably made changes to it).
>>>>>>
>>>>>> Regards,
>>>>>>     Daniel
>>>>>>
>>>>>> Am 06.10.2017 5:27 nachm. schrieb "Mario Meissner" <
>>>>>> mr.rash....@gmail.com>:
>>>>>>
>>>>>>> I am getting compilation errors again, on 70355.
>>>>>>>
>>>>>>> fatal error C1083: Cannot open source file:
>>>>>>> 'C:\Users\MEISSNERBLANCOJOHANN\Desktop\brlcad\trunk2\src\librt\spectrum.c':
>>>>>>> No such file or directory
>>>>>>>
>>>>>>> And indeed no such file exists (I cannot find it anywhere).
>>>>>>>
>>>>>>> Hints? :/
>>>>>>>
>>>>>>> On 4 October 2017 at 23:14, Mario Meissner <mr.rash....@gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> So I ended up designing a new algorithm and would love to get some
>>>>>>>> feedback before I implement it.
>>>>>>>> Here it goes:
>>>>>>>>
>>>>>>>> * Sort the points from 'left to right' (distance of projection to
>>>>>>>> inhit)
>>>>>>>> * Look for the closest point to inhit.
>>>>>>>> * Add this point to the list of contributions.
>>>>>>>> * Discard every point 'to the left' of this point.
>>>>>>>> * From the next point onwards, do, for each point:
>>>>>>>>     * Take the next point and the last one we added to
>>>>>>>> contributions
>>>>>>>>     * Check if the intersection of their middle-plain and the
>>>>>>>> segment lays inside the segment and
>>>>>>>>
>>>>>>>> has no closer point than the two we are using..
>>>>>>>>  * If so, this point's region is crossed by the segment, so we add
>>>>>>>> it to the contribution list.
>>>>>>>>  * If not, then this point's region is not crossed.
>>>>>>>> * Look for the closest point to the outhit.
>>>>>>>> * If this point is the same as the last point we added, finish.
>>>>>>>> * If not, add it and finish.
>>>>>>>>
>>>>>>>> Attached go three pics that are a sample of the example situations
>>>>>>>> I used to develop the algorithm.
>>>>>>>> Ideally, I would need someone to find a situation where this
>>>>>>>> algorithm fails.
>>>>>>>>
>>>>>>>> Thank you in advance,
>>>>>>>> Mario.
>>>>>>>>
>>>>>>>> On 3 October 2017 at 04:14, Christopher Sean Morrison <
>>>>>>>> brl...@mac.com> wrote:
>>>>>>>>
>>>>>>>>>
>>>>>>>>> Hi Mario,
>>>>>>>>>
>>>>>>>>> As always, thank you for all the updates.  Here are some responses
>>>>>>>>> to questions you posted last week.
>>>>>>>>>
>>>>>>>>> I'm passing two arguments to my segment_density, which are inhit
>>>>>>>>>> and outhit points.
>>>>>>>>>> I now realize that these are pointers and could change
>>>>>>>>>> unexpectedly within the call if we run multiple threads.
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> The in/out hit point pointers shouldn’t be changing within the
>>>>>>>>> call.  I suspect something else is going on...
>>>>>>>>>
>>>>>>>>>  How can I make this safe? Should I BU_GET safe heap memory to
>>>>>>>>>> store the points for the call, or pass the coordinates one by one as 
>>>>>>>>>> actual
>>>>>>>>>> numbers?
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> In general, the calling code should pass structures to functions
>>>>>>>>> (which can be on the heap or on the stack) and those functions should 
>>>>>>>>> just
>>>>>>>>> fill them in or read from them (but not both, separate functions for
>>>>>>>>> reading and writing).
>>>>>>>>>
>>>>>>>>> For debugging purposes I would like to set only one thread for
>>>>>>>>> rtweight, how can I do that? The fact that I get 4 or more threads 
>>>>>>>>> running
>>>>>>>>> my segment_density at the same time makes it difficult to track down 
>>>>>>>>> issues.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Check out rtweight’s manual page (run "brlman rtweight” outside
>>>>>>>>> mged), it’s the -P# option.
>>>>>>>>>
>>>>>>>>> Cheers!
>>>>>>>>> Sean
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> ------------------------------------------------------------
>>>>>>>>> ------------------
>>>>>>>>> Check out the vibrant tech community on one of the world's most
>>>>>>>>> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>>>>> _______________________________________________
>>>>>>>>> BRL-CAD Developer mailing list
>>>>>>>>> brlcad-devel@lists.sourceforge.net
>>>>>>>>> https://lists.sourceforge.net/lists/listinfo/brlcad-devel
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>> ------------------------------------------------------------
>>>>>>> ------------------
>>>>>>> Check out the vibrant tech community on one of the world's most
>>>>>>> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>>> _______________________________________________
>>>>>>> BRL-CAD Developer mailing list
>>>>>>> brlcad-devel@lists.sourceforge.net
>>>>>>> https://lists.sourceforge.net/lists/listinfo/brlcad-devel
>>>>>>>
>>>>>>>
>>>>>> ------------------------------------------------------------
>>>>>> ------------------
>>>>>> Check out the vibrant tech community on one of the world's most
>>>>>> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>>>>>> _______________________________________________
>>>>>> BRL-CAD Developer mailing list
>>>>>> brlcad-devel@lists.sourceforge.net
>>>>>> https://lists.sourceforge.net/lists/listinfo/brlcad-devel
>>>>>>
>>>>>>
>>>>>
>>>>
>> <rtexample.diff>--------------------------------------------
>> ----------------------------------
>> Check out the vibrant tech community on one of the world's most
>> engaging tech sites, Slashdot.org <http://slashdot.org/>! http://sd
>> m.link/slashdot_______________________________________________
>> BRL-CAD Developer mailing list
>> brlcad-devel@lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/brlcad-devel
>>
>>
>>
>> ------------------------------------------------------------
>> ------------------
>> Check out the vibrant tech community on one of the world's most
>> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
>> _______________________________________________
>> BRL-CAD Developer mailing list
>> brlcad-devel@lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/brlcad-devel
>>
>>
Index: src/rt/viewweight.c
===================================================================
--- src/rt/viewweight.c (リビジョン 70365)
+++ src/rt/viewweight.c (作業コピー)
@@ -48,7 +48,28 @@
 #include "./rtuif.h"
 #include "./ext.h"
 
+fastf_t MARGIN = 0.0001;
+struct density_point {
+       struct bu_list l;
+       point_t point;
+       fastf_t density;
+};
 
+struct contribution {
+       struct density_point * density_point;
+       fastf_t inhit_dist;
+       point_t left_boundary;
+       point_t right_boundary;
+};
+
+struct boundary_t {
+       int is_valid;
+       point_t boundary_point;
+};
+
+struct density_point * user_plist;
+const int MAX_POINTS = 99;
+
 extern struct resource resource[];
 
 /* Viewing module specific "set" variables */
@@ -56,9 +77,292 @@
     {"",       0, (char *)0,   0,      BU_STRUCTPARSE_FUNC_NULL, NULL, NULL}
 };
 
-
 const char title[] = "RT Weight";
 
+/*
+Reads density points from stdin and stores them into density_point structs
+Receives the head of a bu_list of density_point structs
+and appends all points generated during the reading session to the list
+Returns the number of successfully read points
+*/
+int
+read_density_points(struct density_point * dplist_head)
+{
+
+       /* If we are given no proper bu_lists ... */
+       if (!BU_LIST_IS_INITIALIZED(&(dplist_head->l))) {
+               bu_log("No valid heads of bu_list has been provided.\n");
+               /* Maybe also check if bu_list is of correct type? */
+       }
+
+       /* Local variables */
+       int count = 0;
+       struct density_point * entry;
+       int correct;
+       int reading = 1;
+       char input[30];
+       struct density_point * entries[99]; //Hardcoded max for now
+
+                                                                               
/* Outer loop. We come back here whenever user inputs wrong format
+                                                                               
and until he types exit */
+       while (reading) {
+               bu_log("Please insert density points in the format:");
+               bu_log(" X Y Z density_value(in g / cm3). \n");
+               bu_log("Type 'exit' to finish. \n");
+               correct = 1;
+
+               /* Inner loop. While user inputs points in correct format we 
loop here */
+               while (correct) {
+                       bu_fgets(input, 30, stdin);
+                       /* If user typed continue we need to exit both loops */
+                       if (!strcmp(input, "exit\r")) {
+                               bu_log("Finished inserting points.\n");
+                               correct = 0;
+                               reading = 0;
+                       }
+
+                       /* User didn't type continue */
+                       else {
+
+                               /* Set up new entry */
+                               BU_GET(entry, struct density_point);
+
+                               /* If the input is in correct format we load it 
and continue */
+                               if (sscanf(input, "%lf %lf %lf %lf", 
&entry->point[0],
+                                       &entry->point[1], &entry->point[2], 
&entry->density) == 4) {
+                                       BU_LIST_PUSH(&(dplist_head->l), 
&(entry->l));
+                                       entries[count] = entry;
+                                       count++;
+                               }
+
+                               /* Otherwise we dont load it and jump back to 
outer loop
+                               to print usage again */
+                               else {
+                                       correct = 0;
+                               }
+                       }
+               } /* Close inner loop */
+       } /* Close outer loop */
+
+
+
+       return count;
+}
+
+/* Check if c is between a and b */
+int isBetween(point_t c, point_t a, point_t b) {
+       vect_t c_minus_a, b_minus_a;
+       VSUB2(c_minus_a, c, a);
+       VSUB2(b_minus_a, b, a);
+       fastf_t dot_product = VDOT(c_minus_a, b_minus_a);
+       fastf_t squared_length_ba = MAGNITUDE(b_minus_a)*MAGNITUDE(b_minus_a);
+       return dot_product > 0 && dot_product < squared_length_ba;
+}
+
+int compare_dpoints(const void *a, const void *b) {
+       struct contribution *x = *(struct contribution **)a;
+       struct contribution *y = *(struct contribution **)b;
+       return x->inhit_dist - y->inhit_dist;
+}
+
+/* Receives two points: current and previous, and a segment, and finds out
+if the segment will cross the 'current point's  region or not.
+Uses global userpoint list variable */
+void region_boundary(struct boundary_t * boundary, int current_index, int 
previous_index,
+       struct contribution ** points_array, int array_size, point_t 
segment_start, point_t segment_end) {
+
+       point_t current, previous;
+       VMOVE(current, points_array[current_index]->density_point->point);
+       VMOVE(previous, points_array[previous_index]->density_point->point);
+       int valid_intersection, valid_boundary;
+       point_t plane_point, boundary_point;
+       vect_t plane_normal;
+       fastf_t epsilon = 1e-6;
+       fastf_t fac;
+       vect_t u, w;
+       VSUB2(plane_normal, current, previous);
+       VSCALE(plane_normal, plane_normal, 0.5);
+       VADD2(plane_point, previous, plane_normal);
+       VSUB2(u, segment_end, segment_start);
+       fastf_t dot = VDOT(plane_normal, u);
+       VSUB2(w, segment_start, plane_point);
+       fac = -VDOT(plane_normal, w) / dot;
+       VSCALE(u, u, fac);
+       VADD2(boundary_point, segment_start, u);
+       /* Check if intersection is valid */
+       valid_intersection = abs(dot) > epsilon &&
+               isBetween(boundary_point, segment_start, segment_end);
+
+       /* Break out if we don't intersect */
+       if (!valid_intersection) {
+               boundary->is_valid = 0;
+               return;
+       }
+
+       valid_boundary = 1;
+       int index;
+       vect_t vector;
+       point_t point;
+       VSUB2(vector, boundary_point, current);
+       /* This is the distance between the intersection and our points */
+       fastf_t boundary_distance = MAGNITUDE(vector);
+       fastf_t point_distance;
+       /* Now check distances to see if boundary is valid */
+       for (index = 0; index < array_size; index++) {
+               if (index == current_index || index == previous_index) continue;
+               VMOVE(point, points_array[index]->density_point->point);
+               VSUB2(vector, point, boundary_point);
+               point_distance = MAGNITUDE(vector);
+               if (point_distance < boundary_distance) {
+                       valid_boundary = 0;
+                       break;
+               }
+       }
+
+       if (valid_boundary) {
+               boundary->is_valid = 1;
+               VMOVE(boundary->boundary_point, boundary_point);
+       }
+       else {
+               boundary->is_valid = 0;
+       }
+}
+
+int search_closest_point(point_t point, struct contribution ** point_list, int 
array_size) {
+       int pos_min;
+       fastf_t dist;
+       fastf_t min_dist = FLT_MAX;
+       vect_t vect;
+       point_t candidate;
+       for (int i = 0; i < array_size; i++) {
+               VMOVE(candidate, point_list[i]->density_point->point);
+               VSUB2(vect, candidate, point);
+               dist = MAGNITUDE(vect);
+               if (dist < min_dist) {
+                       pos_min = i;
+                       min_dist = dist;
+               }
+       }
+       return pos_min;
+}
+
+
+
+/* Receives a segment and returns average density
+Uses memory pointer to user point list */
+fastf_t segment_density(point_t inhit, point_t outhit) {
+
+
+       /* If no points, mass is 0 */
+       if (BU_LIST_IS_EMPTY(&(user_plist->l))) {
+               bu_log("No points were provided, thus object has no mass.\n");
+               return 0;
+       }
+
+       /* If there is only one point we skip the paraphernalia */
+       if (bu_list_len(&user_plist->l) == 1) {
+               struct density_point * dp;
+               dp = (struct density_point *)user_plist->l.forw;
+               return dp->density;
+       }
+
+       vect_t full_segment_vector;
+       VSUB2(full_segment_vector, outhit, inhit);
+       fastf_t full_segment_length = MAGNITUDE(full_segment_vector);
+
+       /* Projections Loop prep */
+       struct contribution * new_contribution;
+       struct density_point * dpoint_iter;
+       vect_t dpoint_vect;
+       vect_t projection_vector;
+       vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+                                                        /* Bu_calloc has given 
me some problems, so temporary hardcoded max */
+                                                        /* FIXME: I should be 
allocating space dynamically */
+
+
+                                                        /* Array variables */
+                                                        /* Array containing 
the points that will contribute (i.e. we will cross the region) */
+       struct contribution * contributions[99];
+       /* Array containing all points */
+       struct contribution * userpoints[99];
+       int userpoints_size = 0;
+       int contributions_size = 0;
+       /* Index of the userpoint we are currently working with */
+       int userpoints_index = 0;
+       /* Index of the next position where we will store a contribution */
+       int contributions_index = 0;
+
+       /* Copy all points to the array (and convert them to contribution type) 
*/
+       for (BU_LIST_FOR(dpoint_iter, density_point, &(user_plist->l))) {
+               BU_GET(new_contribution, struct contribution);
+               /* Create userpoint vector */
+               VSUB2(dpoint_vect, dpoint_iter->point, inhit);
+               /* Compute projected distance to inhit */
+               VPROJECT(dpoint_vect, full_segment_vector,
+                       projection_vector, paral_projection);
+               new_contribution->inhit_dist = MAGNITUDE(projection_vector);
+               new_contribution->density_point = dpoint_iter;
+               userpoints[userpoints_index] = new_contribution;
+               userpoints_size++;
+               userpoints_index++;
+       }
+
+       /* STAGES ONE AND TWO - CREATE LIST OF FINAL CONTRIBUTION POINTS */
+
+       /* Sort all points from 'left to right', i.e. distance to inhit */
+       qsort(userpoints, userpoints_size, sizeof(struct contribution *), 
compare_dpoints);
+
+       /* Reset the indeces */
+       userpoints_index = 0;
+       contributions_index = 0;
+
+       /* Search the closest point to inhit since it's gonna be the first 
region we cross */
+       userpoints_index = search_closest_point(inhit, userpoints, 
userpoints_size);
+       contributions[contributions_index] = userpoints[userpoints_index];
+       VMOVE(contributions[contributions_index]->left_boundary, inhit);
+       contributions_size++;
+       contributions_index++;
+
+       /* Position inside userpoints where the last point we added to 
contributions is */
+       int last_added_userpoint = userpoints_index;
+       userpoints_index++;
+       /* Now loop through the rest in order, checking when we cross the next 
region */
+       struct boundary_t * boundary;
+       BU_GET(boundary, struct boundary_t);
+       for (userpoints_index; userpoints_index < userpoints_size; 
userpoints_index++) {
+               region_boundary(boundary, userpoints_index, 
last_added_userpoint,
+                       userpoints, userpoints_size, inhit, outhit);
+               if (boundary->is_valid) {
+                       /* If valid, add current userpoint to the contributions 
list */
+                       contributions[contributions_index] = 
userpoints[userpoints_index];
+                       
VMOVE(contributions[contributions_index]->left_boundary, 
boundary->boundary_point);
+                       VMOVE(contributions[contributions_index - 
1]->right_boundary, boundary->boundary_point);
+                       contributions_index++;
+                       contributions_size++;
+                       last_added_userpoint = userpoints_index;
+               }
+       }
+       BU_FREE(boundary, struct boundary_t);
+
+       /* Finally the last right_boundary is outhit */
+       VMOVE(contributions[contributions_size - 1]->right_boundary, outhit);
+
+       /* STAGE THREE: COMPUTE MASS ACCORDING TO CONTRIBUTION DISTRIBUTION */
+       vect_t contribution_segment;
+       fastf_t average_density, contribution_segment_len;
+       average_density = 0.0;
+
+
+       /* Loop through contributions and compute average density */
+       for (int i = 0; i<contributions_size; i++) {
+               VSUB2(contribution_segment, contributions[i]->left_boundary, 
contributions[i]->right_boundary);
+               contribution_segment_len = MAGNITUDE(contribution_segment);
+               average_density += contributions[i]->density_point->density * 
contribution_segment_len / full_segment_length;
+       }
+       return average_density;
+
+}
+
 void
 usage(const char *argv0)
 {
@@ -100,7 +404,6 @@
     fastf_t volume;
 };
 
-
 extern int rpt_overlap;        /* report region verbosely */
 extern fastf_t cell_width;             /* model space grid cell width */
 extern fastf_t cell_height;            /* model space grid cell height */
@@ -144,9 +447,9 @@
 
        /* if we don't have a valid material density to work with, use a 
default material */
        if (density[reg->reg_gmater] < 0) {
-           bu_log("Material type %d used, but has no density file entry.\n", 
reg->reg_gmater);
-           bu_log("  (region %s)\n", reg->reg_name);
-           bu_log("  Mass is undefined.\n");
+           //bu_log("Material type %d used, but has no density file entry.\n", 
reg->reg_gmater);
+           //bu_log("  (region %s)\n", reg->reg_name);
+           //bu_log("  Mass is undefined.\n");
            bu_semaphore_acquire(BU_SEM_SYSCALL);
            reg->reg_gmater = 0;
            bu_semaphore_release(BU_SEM_SYSCALL);
@@ -157,18 +460,22 @@
            const fastf_t partition_volume = depth * cell_height * cell_width / 
(hypersample + 1);
 
            /* convert reg_los percentage to factor */
-           const fastf_t los_factor = (fastf_t)reg->reg_los * 0.01;
-
+               /*FIXME: reg->reg_los is always 0 for whatever reason */
+           //const fastf_t los_factor = (fastf_t)reg->reg_los * 0.01;
+               const fastf_t los_factor = 1;
            dp->volume = partition_volume;
            VBLEND2(dp->centroid, 0.5, ihitp->hit_point, 0.5, ohitp->hit_point);
-
-           if (density[reg->reg_gmater] >= 0) {
+               
+               //puts("about to call segment_density");
                /* density needs to be converted from g/cm^3 to g/mm^3 */
-               const fastf_t density_factor = density[reg->reg_gmater] * 0.001;
-
+               const fastf_t density_factor = 
segment_density(ihitp->hit_point, ohitp->hit_point) * 0.001;
+               //const fastf_t density_factor = 0.0000007; // wood
+               //bu_log("Read %lf mass from segment_density\n", 
density_factor);
                /* Compute mass in terms of grams */
                dp->weight = partition_volume * los_factor * density_factor;
-           }
+               //bu_log("Partition volume: %lf, los_factor: %lf, 
density_factor %lf \n Stored %lf in dp->weight \n", 
+                       //partition_volume, los_factor, density_factor, 
dp->weight);
+           
        }
     }
     return 1;  /* report hit to main routine */
@@ -292,7 +599,21 @@
     if (minus_o) {
        fclose(outfp);
     }
+       /* Prepare user point list and let user fill it */
+       BU_GET(user_plist, struct density_point);
+       BU_LIST_INIT(&(user_plist->l));
+       read_density_points(user_plist);
 
+       //DEBUG
+       struct density_point * user_plist_iter;
+       puts("Main after user input. These are the loaded points:");
+       for (BU_LIST_FOR(user_plist_iter, density_point, &(user_plist->l))) {
+               VPRINT("density_point", user_plist_iter->point);
+       }
+       //END DEBUG
+
+
+       puts("Finish view_init");
     return 0;          /* no framebuffer needed */
 }
 
@@ -451,7 +772,7 @@
 
            weight *= conversion;
            total_weight += weight;
-
+               bu_log("added %lf to total weight\n", weight);
            BU_ALLOC(ptr, fastf_t);
            *ptr = weight;
            /* FIXME: shouldn't the existing reg_udata be bu_free'd first (see 
previous loop) */
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
BRL-CAD Developer mailing list
brlcad-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/brlcad-devel

Reply via email to