Hi Daniel.
Thank you for your response. Negative values may have been only due to
using the wrong formatting tag while printing values. I have cleaned up
some unused variables and fixed the intersection by using your suggestion.
Mass is 0 in the total_weight variable, as if the code wasn't really adding
up my segment_density values correctly. I might have deleted some important
code related to the storing or reading of the values. I will double check
on that now.
Attached goes current version of my code.
Mario.
On 21 September 2017 at 13:03, Daniel Roßberg <danielmrossb...@gmail.com>
wrote:
> Hi Mario,
>
> I have been trying to track down the reason why some densities are
>> negative. A few questions I have ended up asking myself:
>>
>> - Is it possible for MAGNITUDE() to return a negative number? It would
>> make no sense to me, but I suspect it's happening.
>>
>
> No, MAGNITUDE() is the square root of a (necessarily) positive value.
> You can safely assume it to be not negative.
>
> - In my little function intersection I noticed that I am returning a
>> pointer to a local variable, that could be destroyed in any moment after
>> the call ends. Instead of that, I want to return a safe variable. How do I
>> do that in this case?
>>
>> I don't understand why it doesn't let me return directly the variable
>> point_t. It's a pointer, because it's an array, right? As it didn't let me
>> return the point, I decided to return a pointer to the point instead,
>> however as I mentioned this is not safe, so I have to change it now. Any
>> explanation on why I can't return the point, and what I should do to make
>> it safe?
>>
>
> Well, there is a standard solution for this kind of issue: Putting the
> data on the heap, e.g.
> point_t *intersect = (point_t *) bu_malloc(sizeof(point_t),
> "intersection");
> with the disadvantage that the caller has explicitly to free the memory
> with bu_free().
>
> It's true that C functions can't return arrays (because of historic
> reason). But, returning an array wouldn't help you here as you are using
> the function's return as a flag for intersection too. On the other hand,
> structs can be returned:
> struct intersection_t {
> int intersect;
> point_t intersection_point;
> };
>
> struct intersection_t intersection(point_t p0, point_t p1, point_t
> p_co, vect_t p_no) {
> ....
> should work.
>
>
>> Finally, I think I could need some help regarding the negative values
>> issue, at least once I'm sure the two points above aren't the cause.
>>
>
> Feel free to ask ;)
>
>
> Regards,
> Daniel
>
> ------------------------------------------------------------
> ------------------
> 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 (リビジョン 70318)
+++ src/rt/viewweight.c (作業コピー)
@@ -48,7 +48,32 @@
#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;
+};
+
+struct intersection_t {
+ int intersect;
+ point_t intersection_point;
+};
+
+struct density_point * user_plist;
+const int MAX_POINTS = 99;
+
extern struct resource resource[];
/* Viewing module specific "set" variables */
@@ -59,6 +84,280 @@
const char title[] = "RT Weight";
+int compare_dpoints(const void *, const void *);
+struct intersection_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\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;
+}
+
+/* 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);
+ //bu_log("inhit_dist: %lf \n", new_projection->inhit_dist);
+ /* 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);
+ struct intersection_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);
+ if (intersectionp == NULL) {
+ bu_log("Intersection error\n");
+ return 0;
+ }
+ VSUB2(intersect_to_original,
intersectionp->intersection_point, 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)
{
@@ -100,7 +399,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 +442,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);
@@ -162,13 +460,13 @@
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);
+ //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;
- }
+
}
}
return 1; /* report hit to main routine */
@@ -292,7 +590,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 */
}
@@ -608,7 +920,34 @@
{
}
+struct intersection_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;
+ struct intersection_t * intersect_struct;
+ BU_GET(intersect_struct, struct intersection_t);
+ 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_struct->intersection_point, p0, u);
+ intersect_struct->intersect = 0;
+ return intersect_struct;
+ }
+ else {
+ return NULL;
+ }
+}
+
+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;
+}
+
/*
* Local Variables:
* mode: C
------------------------------------------------------------------------------
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
