I'm working on integrating the code into viewweight. As a first step I thought it would be a good idea to not touch the existing code (and let it read from the .density file) and just override the hit() function to call my segment_density instead of taking the file data. However I'm getting segmentation faults whenever I try to run the code. This happened to me before and I managed to solve it by changing the codification of the .density file, however this time I don't seem to be able to fix it.
Here is result screencap, and I also attach the modified viewweight and the density file. https://puu.sh/xAmiE/b5fa0dfab8.png Mario. On 8 September 2017 at 21:04, Mario Meissner <mr.rash....@gmail.com> wrote: > Peopleee! > > Here goes today's work. I'm proud to announce the example I modeled a few > days ago now works correctly!! > Code still contains my debug prints and lots of random 'puts' to know > where and why stuff crashes. I feel I will need it in the future so for now > it remains. > I probably should keep looking for flaws by setting up other examples, and > while doing so I should also clean it up a bit since the code grew over 600 > lines now, quite difficult to read and understand. > > Screenshot: > https://puu.sh/xuIXk/9626ce2fb6.png > > Mario. > > On 8 September 2017 at 19:56, Mario Meissner <mr.rash....@gmail.com> > wrote: > >> Hi Daniel! >> >> Thank you a lot for your reply. I was focused on generating the bisector >> to intersect two lines, but intersecting the actual plane is much simpler. >> I made a little function to intersect a plane with a line (partially taken >> from this stackoverflow post >> <https://stackoverflow.com/questions/5666222/3d-line-plane-intersection>, >> so credits to user ideasman42 >> <https://stackoverflow.com/users/432509/ideasman42>). I am now finishing >> of the 'patch' to my code that is supposed to now correctly skip >> unnecessary points and create a list with those who's regions we will >> actually cross. Code probably tomorrow. >> >> Mario. >> >> On 8 September 2017 at 09:35, Daniel Roßberg <danielmrossb...@gmail.com> >> wrote: >> >>> Hi Mario, >>> >>> In 3D space, rectangular to a vector you have a plane and the vector is >>> a normal vector of this plane. In this plane are all vectors which are >>> perpendicular to your first vector, i.e. the vectors which products with >>> your first vector are 0. This gives you a first equation with two >>> variables. >>> >>> Next, you want to intersect the plane with a ray. I.e., one vector from >>> your plan has to be a multiple of the ray which gives you another equation. >>> >>> However, there should be examples in the BRL-CAD source code for this. >>> E.g., the ray-trace of the half primitive is very similar to your problem. >>> There, the intersection of a plane, given by a normal vector and a point in >>> it, and a ray, given by a start point and a direction, will be computed. >>> >>> I hope this helps. >>> >>> >>> Regards, >>> Daniel >>> >>> >>> >>> I'm facing a wall. I need to find a perpendicular vector with the >>> condition that it will intersect a third vector. >>> Basically, I have the ray and a vector drawn from my current point to >>> the previous point. >>> I want the perpendicular bisector of that vector to cross my ray, so >>> that I can obtain the intersection between the bisector and the ray. >>> >>> For that I need two things: >>> 1) Obtain the perpendicular vector. in 3D there are infinite options, >>> but only one will cross the ray. How do I get the one I need? >>> 2) Obtain the intersection point. >>> >>> I browsed through vmath.h but I didn't find anything that could help me. >>> Do you know any macro that could be useful here? >>> Otherwise I will need to make some of my own for this! If so, any ideas >>> on how to approach it? >>> >>> Mario. >>> >>> >>> ------------------------------------------------------------ >>> ------------------ >>> 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 (revision 70280)
+++ src/rt/viewweight.c (working copy)
@@ -48,7 +48,30 @@
#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 projection {
+ //If we store the density here again we might not need the original
+ //fastf_t density;
+ struct density_point * original;
+ point_t point;
+ vect_t vect;
+ fastf_t inhit_dist;
+ fastf_t prev_dist;
+ fastf_t origin_dist;
+};
+
+/*FIXME: This is global for now because I need this pointer through all the
+hit calls. How can I make it available for all calls without a global
variable? */
+struct density_point * user_plist;
+
+const int MAX_POINTS = 99;
+
extern struct resource resource[];
/* Viewing module specific "set" variables */
@@ -59,6 +82,273 @@
const char title[] = "RT Weight";
+int compare_dpoints(const void *, const void *);
+point_t * intersection(point_t, point_t, point_t, vect_t);
+
+/*
+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\n")) {
+ 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;
+}
+
+/* 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 user points available it means either no mass or default
density.
+ I will display a message indicating no points are defined so no mass*/
+ if (BU_LIST_IS_EMPTY(&(user_plist->l))) {
+ bu_log("No points available, thus object has no mass.\n");
+ return 0;
+ }
+
+ /* DEBUG */
+ VPRINT("inhit", inhit);
+ VPRINT("outhit", outhit);
+
+ /* Segment vector (for projections) */
+ vect_t segment_vect;
+ VSUB2(segment_vect, outhit, inhit);
+ fastf_t total_segment_len = MAGNITUDE(segment_vect);
+
+ /* Projections Loop prep */
+ struct projection * new_projection;
+ struct density_point * dpoint_iter_ext; //External
+ struct density_point * dpoint_iter_mid; //Middle
+ struct density_point * dpoint_iter_int; //Internal
+ vect_t dpoint_vect;
+ vect_t prev_to_current_vect;
+ vect_t proj_to_original_vect;
+ vect_t temp;
+ vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+ int array_pointer = 0;
+ int is_valid;
+ /* Bu_calloc has given me some problems, so temporary hardcoded max */
+ /* FIXME: I should be allocating space dynamically */
+ struct projection * proj_array[99];
+
+ /* Loop through points and project them onto segment if they need to be
*/
+ for (BU_LIST_FOR(dpoint_iter_ext, density_point, &(user_plist->l))) {
+ is_valid = 1;
+ BU_GET(new_projection, struct projection);
+ /* Create userpoint vector */
+ VSUB2(dpoint_vect, dpoint_iter_ext->point, inhit);
+ new_projection->original = dpoint_iter_ext;
+ /* Project vector onto segment */
+ VPROJECT(dpoint_vect, segment_vect,
+ new_projection->vect, paral_projection);
+ /* Compute absolute point of projection */
+ VADD2(new_projection->point, new_projection->vect, inhit);
+ new_projection->inhit_dist = MAGNITUDE(new_projection->vect);
+ /* Compute distance to original */
+ VSUB2(proj_to_original_vect, dpoint_iter_ext->point,
new_projection->point);
+ new_projection->origin_dist = MAGNITUDE(proj_to_original_vect);
+
+ /* Check if we need to project or not based on distance to all
other points */
+ for (BU_LIST_FOR(dpoint_iter_mid, density_point,
&(user_plist->l))) {
+ VSUB2(temp, new_projection->point,
dpoint_iter_mid->point);
+ /* If distance to some user point is shorter than to
the origin then we need to check boundaries */
+ if (dpoint_iter_ext != dpoint_iter_mid &&
+ MAGNITUDE(temp) < new_projection->origin_dist) {
+
+ /* If we are here it means projection has
closer neighbor than original,
+ we now check if intersection falls inside
segment and has closer neighbor.
+ Also it should not be our first point. */
+ if (array_pointer > 0) {
+ puts("array pointer bigger than zero");
+ vect_t prev_to_curr;
+ point_t plane_point;
+ struct projection * previous =
proj_array[array_pointer - 1];
+ VPRINT("previous",
previous->original->point);
+ VPRINT("current",
new_projection->original->point);
+ puts("read current and previous, now
computing plane");
+ VSUB2(prev_to_curr,
new_projection->original->point, previous->original->point);
+ point_t * intersectionp;
+ vect_t intersect_to_original;
+ fastf_t distance_to_originals;
+ VSCALE(prev_to_curr, prev_to_curr, 0.5);
+ puts("computing intersection");
+ VADD2(plane_point, previous->point,
prev_to_curr);
+ intersectionp = intersection(inhit,
outhit, plane_point, new_projection->original->point);
+ VSUB2(intersect_to_original,
*intersectionp, new_projection->original->point);
+ distance_to_originals =
MAGNITUDE(intersect_to_original);
+ /* Look for points closer to
intersection than the two original ones */
+ puts("arrived at inner for");
+ for (BU_LIST_FOR(dpoint_iter_int,
density_point, &(user_plist->l))) {
+ VSUB2(temp,
new_projection->point, dpoint_iter_int->point);
+ if (dpoint_iter_ext !=
dpoint_iter_int &&
+ MAGNITUDE(temp) <
distance_to_originals + MARGIN) {
+ puts("inner loop found
something");
+ /* FIXME: Deallocate
space we used for new_projection */
+ is_valid = 0;
+ bu_log("Skipped
projection. \n");
+ break;
+ }
+ }
+
+ }
+ }
+ /* We are at the end of one iteration on the outer
check loop.
+ If we detected a nonvalid point we should stop looking.
*/
+ if (!is_valid) break;
+ }
+ if (is_valid) {
+
+ if (array_pointer == 0) {
+ new_projection->prev_dist =
new_projection->inhit_dist;
+ }
+ else {
+ VSUB2(prev_to_current_vect,
new_projection->vect,
+ proj_array[array_pointer - 1]->vect);
+ new_projection->prev_dist =
MAGNITUDE(prev_to_current_vect);
+ }
+ proj_array[array_pointer] = new_projection;
+ array_pointer++;
+ /* DEBUG */
+ VPRINT("projection vector", new_projection->vect);
+ }
+
+
+
+
+ }
+
+ /* Sort projection list based on distance to inhit */
+ /* DEBUG */
+ puts("Unordered Projection distances:");
+ for (int i = 0; i < array_pointer; i++) {
+ bu_log("Point at distance %lf \n", proj_array[i]->inhit_dist);
+ }
+
+ qsort(proj_array, array_pointer, sizeof(struct projection *),
compare_dpoints);
+
+ /* DEBUG */
+ puts("Ordered projection distances:");
+ for (int i = 0; i < array_pointer; i++) {
+ bu_log("Point at distance %lf \n", proj_array[i]->inhit_dist);
+ }
+
+ /* Density Computation Loop preparation */
+ struct projection * curr_proj;
+ fastf_t acc_avg_density = 0.0;
+ struct projection * next_proj;
+
+ /* WARNING: From now on I refer to 'segment' as each piece of ray
inbetween two
+ points (projection or in/out */
+
+ /* First segment, we take everything as having density of first point */
+ curr_proj = proj_array[0];
+ acc_avg_density += curr_proj->original->density *
+ curr_proj->inhit_dist / total_segment_len;
+
+ /* DEBUG */
+ bu_log("First projection. Proportion %f, added %lf to accumulator.\n",
+ curr_proj->inhit_dist / total_segment_len,
curr_proj->original->density *
+ curr_proj->inhit_dist / total_segment_len);
+
+ /* Distance from current projection point to boundary */
+ fastf_t boundary_dist;
+ /* Distance between current projection point and next one */
+ fastf_t partial_segment_length;
+ vect_t partial_segment;
+ puts("loop start");
+ printf("Array pointer is %d \n", array_pointer);
+ next_proj = proj_array[0];
+ /* Loop through remaining segments. For each segment, compute the
boundary
+ and take density of each side accordingly, add it to accumulator */
+ for (int i = 0; i < array_pointer - 1; i++) {
+ curr_proj = proj_array[i];
+ next_proj = proj_array[i + 1];
+ VSUB2(partial_segment, curr_proj->point, next_proj->point);
+ partial_segment_length = MAGNITUDE(partial_segment);
+ boundary_dist = (next_proj->origin_dist *
next_proj->origin_dist +
+ partial_segment_length * partial_segment_length -
+ curr_proj->origin_dist * curr_proj->origin_dist) /
+ (2 * partial_segment_length);
+ acc_avg_density += curr_proj->original->density *
+ boundary_dist / total_segment_len;
+ acc_avg_density += next_proj->original->density *
+ (partial_segment_length - boundary_dist) /
total_segment_len;
+ }
+ puts("end of for loop");
+ curr_proj = next_proj;
+ /* Last segment, from the last point to the outhit point has density
value
+ of this last point */
+ vect_t last_to_outhit;
+ VSUB2(last_to_outhit, outhit, curr_proj->point);
+ acc_avg_density += curr_proj->original->density *
+ MAGNITUDE(last_to_outhit) / total_segment_len;
+ /* DEBUG */
+ bu_log("Last projection. Proportion %f, added %lf to accumulator.\n",
+ curr_proj->prev_dist / total_segment_len,
curr_proj->original->density *
+ curr_proj->prev_dist / total_segment_len);
+
+ /* Return final value */
+ bu_log("Avg density we saw through segment: %lf \n", acc_avg_density);
+ return acc_avg_density;
+}
+
void
usage(const char *argv0)
{
@@ -162,13 +452,19 @@
dp->volume = partition_volume;
VBLEND2(dp->centroid, 0.5, ihitp->hit_point, 0.5, ohitp->hit_point);
+ const fastf_t density_factor =
segment_density(ihitp->hit_point, ohitp->hit_point);
+
+ /* Compute mass in terms of grams */
+ dp->weight = partition_volume * los_factor * density_factor;
+
+ /*
if (density[reg->reg_gmater] >= 0) {
- /* density needs to be converted from g/cm^3 to g/mm^3 */
+ /* density needs to be converted from g/cm^3 to g/mm^3
const fastf_t density_factor = density[reg->reg_gmater] * 0.001;
-
+ */
/* Compute mass in terms of grams */
dp->weight = partition_volume * los_factor * density_factor;
- }
+
}
}
return 1; /* report hit to main routine */
@@ -608,7 +904,32 @@
{
}
+int compare_dpoints(const void *a, const void *b) {
+ struct projection *x = *(struct projection **)a;
+ struct projection *y = *(struct projection **)b;
+ return x->inhit_dist - y->inhit_dist;
+}
+point_t * intersection(point_t p0, point_t p1, point_t p_co, vect_t p_no) {
+ fastf_t epsilon = 1e-6;
+ fastf_t fac;
+ vect_t u, w;
+ point_t intersect;
+ VSUB2(u, p1, p0);
+ fastf_t dot = VDOT(p_no, u);
+ if (abs(dot) > epsilon) {
+ VSUB2(w, p0, p_co);
+ fac = -VDOT(p_no, w) / dot;
+ VSCALE(u, u, fac);
+ VADD2(intersect, p0, u);
+ return &intersect;
+ }
+ else {
+ return NULL;
+ }
+}
+
+
/*
* Local Variables:
* mode: C
.density
Description: Binary data
------------------------------------------------------------------------------ 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
