Sean!
Implemented the math and did some bug hunting. I successfully ran an
N-Point Voronoi example (without vectors for now).
I shoot a ray from above at my 1000mm box. I set density points at heights
600, 400 and 200 with density values of 20, 10 and 5 respectively. Expected
result is 13.5 average density throughout the segment from height 1000 to
0.
Attached goes today's code and output.
Notice that I have left my debugging logs on since I think they will be
very useful for you to see how code behaves, and for me now that I'm gonna
start including vectors.
I desperately need some feedback now, hope you'll have some time to review
this.
https://puu.sh/xaUwS/c2f30fc275.png
Mario.
2017-08-14 19:44 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
> The sorting function is now correctly implemented and works as expected
> with projection structs, sorting them by distance to inhit. Some debug code
> is included that prints the unsorted and then the sorted array of
> projections.
> Also much of the 'basement' for Voronoi is now ready, so now only the
> overlaying math required to compute the average density is missing.
> I think it can be done by tomorrow. I however am not sure at all if this
> is what you meant when we talked about this model, I'm working a bit
> towards unknown territory at the moment.
> Attached goes diff.
> Mario.
>
> 2017-08-11 23:46 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>
>> Hello Sean.
>> Things got a bit messy on the previous thread, I apologize for that.
>> I somehow managed to send emails only to myself for the last three days.
>>
>> I decided to start a new fresh thread and include the text of the last
>> three emails here, as well as the newest code I have. It is a new iteration
>> over email3 with some of today's code commented out so that it compiles and
>> performs the 'complete' functionality I mentioned on my email2. You can get
>> a feel of what I pretend to do next by having a look at these commented out
>> lines of code as well.
>>
>> You can ignore the old thread and just reply to this one email instead.
>> Sorry again for the mess.
>>
>> EMAIL 1
>>
>> Okay,
>>
>> So since stuff doesn’t seem to work very well on viewweight yet, and
>> since we need to heavily modify existing code anyways, I decided to move
>> back to rtexample.
>>
>> If we want to use Voronoi tesselation as our base model and expand it
>> with transition vectors, many we have to rethink how we query data as well.
>> It makes no sense anymore to just query the density at a specific point
>> since we cannot just interpolate. Somewhere in between the two points, a
>> sudden change of density might occur. Interpolation would asume some kind
>> of transition, which there might not be (the value we obtain by
>> interpolation might be way off the real value). So, instead of asking for
>> the density at a specific point, lets give our function the whole segment
>> so that it can correctly return the average density, or even the mass. This
>> way the function knows where the segment starts and ends, can compute the
>> projections and apply Voronoi tesselation to determine the different
>> densities. Then with that info we compute the average density and return it.
>>
>>
>>
>> With rtexample I won’t have the problem of sharing a pointer to my
>> point/vector list since I have a main function where I can do all that. My
>> readDensityPoints() function will also work as expected here.
>>
>>
>>
>> I’ll start working on this, maybe with just 0 and 1 points for now, and
>> hopefully make my first ‘complete’ piece of code.
>>
>>
>> EMAIL 2
>>
>> Sean,
>>
>> I’ve started working on Voronoi model for non-continuous density
>> evaluation using points. For now only 0 and 1 point works, but one complete
>> functionality is there, which is what you requested. User inputs 0 or 1
>> points and then program returns mass of what the ray saw when it crossed
>> geometry. Feedback on the code I started writing to evaluate Voronoi points
>> is appreciated and welcomed. Note that its by no means complete and
>> probably full of mistakes. However knowing that it is the right direction
>> is quite helpful.
>>
>>
>>
>> I think for tomorrow I can finish n-point Voronoi evaluation (especially
>> if I got feedback).
>>
>> I’m also pretty confused about why my user input function doesnt work on
>> viewweight. I’d love to get that working to go back to it.
>>
>> Until then I’ll work on functions to save my structs to a file and load
>> them again.
>>
>>
>>
>> Hereby goes diff rtexample.
>>
>> Cheers!!
>>
>>
>> EMAIL 3
>> Hi!
>> I've tried to clean up the code a bit more, and continued with the
>> implementation of Voronoi points. I hit a wall though:
>> I need to order my projections from the inhit to the outhit points. I was
>> all happy using my bu_list to store them but now I realize that this is a
>> really inconvenient option to manage them if I want to sort them
>> afterwards. I thus switched to an array so that I can use stdlib qsort().
>> Is this a good decision? WIP.
>>
>> As for storing the data into a file... how should I do it so that it has
>> the shape or structure of a 'material density object'?
>> Unfortunately today's code is not yet working. I send it anyways in the
>> hopes of receiving much appreciated feedback.
>>
>> Cheers!!
>> Mario.
>>
>
>
Index: src/rt/rtexample.c
===================================================================
--- src/rt/rtexample.c (revisi¾n: 70059)
+++ src/rt/rtexample.c (copia de trabajo)
@@ -67,258 +67,466 @@
#include <math.h>
#include <string.h>
#include <stdio.h>
+#include <string.h>
#include "vmath.h" /* vector math macros */
#include "raytrace.h" /* librt interface definitions */
+struct density_point {
+ struct bu_list l;
+ point_t point;
+ fastf_t density;
-/**
- * rt_shootray() was told to call this on a hit.
- *
- * This callback routine utilizes the application structure which
- * describes the current state of the raytrace.
- *
- * This callback routine is provided a circular linked list of
- * partitions, each one describing one in and out segment of one
- * region for each region encountered.
- *
- * The 'segs' segment list is unused in this example.
- */
-int
-hit(struct application *ap, struct partition *PartHeadp, struct seg
*UNUSED(segs))
-{
- /* iterating over partitions, this will keep track of the current
- * partition we're working on.
- */
- struct partition *pp;
+};
- /* will serve as a pointer for the entry and exit hitpoints */
- struct hit *hitp;
+struct density_vector {
+ struct bu_list l;
+ vect_t vector;
+ point_t origin;
+ fastf_t factor;
+};
- /* will serve as a pointer to the solid primitive we hit */
- struct soltab *stp;
+struct projection {
+ struct density_point * original;
+ point_t point;
+ vect_t vect_from_inhit;
+ fastf_t distance;
+};
- /* will contain surface curvature information at the entry */
- struct curvature cur = RT_CURVATURE_INIT_ZERO;
+struct contribution {
+ struct bu_list l;
+ fastf_t contrib;
+ struct density_vector * dvp;
+};
- /* will contain our hit point coordinate */
- point_t pt;
+struct density_point * user_plist;
+const int MAX_POINTS = 99;
- /* will contain normal vector where ray enters geometry */
- vect_t inormal;
+int compare_dpoints(const void *, const void *);
+int readDensityPoints(struct density_point *);
- /* will contain normal vector where ray exits geometry */
- vect_t onormal;
+/* Receives a segment and returns average density
+ Uses memory pointer to user point list */
+fastf_t segment_density(point_t inhit, point_t outhit) {
- /* iterate over each partition until we get back to the head.
- * each partition corresponds to a specific homogeneous region of
- * material.
- */
- for (pp=PartHeadp->pt_forw; pp != PartHeadp; pp = pp->pt_forw) {
+ /* DEBUG */
+ VPRINT("inhit", inhit);
+ VPRINT("outhit", outhit);
- /* print the name of the region we hit as well as the name of
- * the primitives encountered on entry and exit.
- */
- bu_log("\n--- Hit region %s (in %s, out %s)\n",
- pp->pt_regionp->reg_name,
- pp->pt_inseg->seg_stp->st_name,
- pp->pt_outseg->seg_stp->st_name );
+ /* Segment vector (for projections) */
+ vect_t segment_vect;
+ VSUB2(segment_vect, outhit, inhit);
+ fastf_t segment_len = MAGNITUDE(segment_vect);
+
+ /* Projections Loop prep */
+ struct projection * new_projection;
+ struct projection * projections;
+ BU_GET(projections, struct projection);
+ /*BU_LIST_INIT(&(projections->l)); we will use array? */
+ vect_t dpoint_vect;
+ vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+ struct density_point * point_iter;
+ int array_pointer = 0;
+ /* 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 */
+ for (BU_LIST_FOR(point_iter, density_point, &(user_plist->l))) {
+ BU_GET(new_projection, struct projection);
+ /* Create userpoint vector */
+ VSUB2(dpoint_vect, point_iter->point, inhit);
+ new_projection->original = point_iter;
+ /* Project vector onto segment */
+ VPROJECT(dpoint_vect, segment_vect,
+ new_projection->vect_from_inhit, paral_projection);
+ /* FIXME: Check that projection actually falls INTO the segment
*/
+ /*BU_LIST_PUSH(&(projections->l), new_projection); we will use
array?*/
+ new_projection->distance =
MAGNITUDE(new_projection->vect_from_inhit);
+ VADD2(new_projection->point, new_projection->vect_from_inhit,
inhit);
+ proj_array[array_pointer] = new_projection;
+ array_pointer++;
+ /* DEBUG */
+ VPRINT("projection vector", new_projection->vect_from_inhit);
+ }
- /* entry hit point, so we type less */
- hitp = pp->pt_inhit;
+ /* 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]->distance);
+ }
+ qsort(proj_array, array_pointer, sizeof(struct projection *),
compare_dpoints);
- /* construct the actual (entry) hit-point from the ray and the
- * distance to the intersection point (i.e., the 't' value).
- */
- VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+ /* DEBUG */
+ puts("Ordered projection distances:");
+ for (int i = 0; i < array_pointer; i++) {
+ bu_log("Point at distance %lf \n", proj_array[i]->distance);
+ }
- /* primitive we encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+ /* -- Go through projections and calculate (with precision) average
density -- */
- /* compute the normal vector at the entry point, flipping the
- * normal if necessary.
- */
- RT_HIT_NORMAL(inormal, hitp, stp, &(ap->a_ray), pp->pt_inflip);
+ /* Loop preparation */
+ struct projection * proj_iter;
+ fastf_t acc_avg_density = 0.0;
+ vect_t temp_vect;
+ struct projection * prev_proj;
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
+ /* First projection point, we take everything from inhit to this point
+ as having density of this point */
+ proj_iter = proj_array[0];
+ acc_avg_density += proj_iter->original->density * proj_iter->distance /
segment_len;
+ /* DEBUG */
+ bu_log("Proportion %f, added %lf to accumulator.\n",
+ proj_iter->distance / segment_len, proj_iter->original->density
*
+ proj_iter->distance / segment_len);
+ prev_proj = proj_iter;
+ /* Loop for the rest. For each point, take half the value from last
point,
+ and half the value from current point */
+ for (int i = 1; i < array_pointer; i++) {
- /* This next macro fills in the curvature information which
- * consists on a principle direction vector, and the inverse
- * radii of curvature along that direction and perpendicular
- * to it. Positive curvature bends toward the outward
- * pointing normal.
- */
- RT_CURVATURE(&cur, hitp, pp->pt_inflip, stp);
+ /* DEBUG */
+ bu_log("Entered projection loop (2 or more points) \n");
- /* print the entry curvature information */
- VPRINT("PDir", cur.crv_pdir);
- bu_log(" c1=%g\n", cur.crv_c1);
- bu_log(" c2=%g\n", cur.crv_c2);
+ proj_iter = proj_array[i];
+ VSUB2(temp_vect, proj_iter->vect_from_inhit,
prev_proj->vect_from_inhit);
+ acc_avg_density += prev_proj->original->density
+ * MAGNITUDE(temp_vect) / (2 * segment_len);
+ acc_avg_density += proj_iter->original->density
+ * MAGNITUDE(temp_vect) / (2 * segment_len);
+ prev_proj = proj_iter;
+ }
+
+ /* Last portion, from the last point to the outhit point has density
value
+ of this last point */
+ VSUB2(temp_vect, segment_vect, proj_iter->vect_from_inhit);
+ acc_avg_density += proj_iter->original->density * MAGNITUDE(temp_vect)
/ segment_len;
+ /* DEBUG */
+ bu_log("Proportion %f, added %lf to accumulator.\n",
+ MAGNITUDE(temp_vect) / segment_len,
proj_iter->original->density *
+ MAGNITUDE(temp_vect) / segment_len);
+ /* 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;
+ }
+
+ else {
+ bu_log("Avg density we saw through segment: %lf \n",
acc_avg_density);
+ return acc_avg_density;
+ }
+}
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+/*
+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 readDensityPoints(struct density_point * dplist_head) {
- /* construct the actual (exit) hit-point from the ray and the
- * distance to the intersection point (i.e., the 't' value).
- */
- VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+ /* If we are given no bu_list... */
+ if (!BU_LIST_IS_INITIALIZED(&(dplist_head->l))) {
+ bu_log("No valid head of bu_list has been provided.\n");
+ /* Maybe also check if bu_list is of type density_points? */
+ }
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+ /* Local variables */
- /* compute the normal vector at the exit point, flipping the
- * normal if necessary.
- */
- RT_HIT_NORMAL(onormal, hitp, stp, &(ap->a_ray), pp->pt_outflip);
+ /* I count the number of times I */
+ int count = 0;
+ struct density_point *entry;
+ int correct, reading = 1;
+ char input[30];
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ /* 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: 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) {
- /* A more complicated application would probably fill in a new
- * local application structure and describe, for example, a
- * reflected or refracted ray, and then call rt_shootray() for
- * those rays.
- */
+ bu_fgets(input, 30, stdin);
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ /* If user typed exit we need to exit both loops */
+ if (!strcmp(input, "exit\n")) {
+ bu_log("Exiting...\n");
+ correct = 0;
+ reading = 0;
+ }
+
+
+
+ /* User didn't type exit */
+ 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));
+ 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;
+} /* Close function */
+
+
+int
+hit(struct application *ap, struct partition *PartHeadp, struct seg
*UNUSED(segs))
+{
+
+ struct partition *pp;
+ struct hit *hitp;
+ struct soltab *stp;
+ struct curvature cur = RT_CURVATURE_INIT_ZERO;
+ point_t pt;
+ vect_t inormal;
+ vect_t onormal;
+
+ /* Area that the ray covers, in mm^2 */
+ const int RAY_AREA = 4; //2mm height, 2mm width
+
+ for (pp = PartHeadp->pt_forw; pp != PartHeadp; pp = pp->pt_forw) {
+ bu_log("\n--- Hit region %s (in %s, out %s)\n",
+ pp->pt_regionp->reg_name,
+ pp->pt_inseg->seg_stp->st_name,
+ pp->pt_outseg->seg_stp->st_name);
+
+ hitp = pp->pt_inhit;
+
+ VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+
+ stp = pp->pt_inseg->seg_stp;
+ RT_HIT_NORMAL(inormal, hitp, stp, &(ap->a_ray), pp->pt_inflip);
+
+ rt_pr_hit(" In", hitp);
+ VPRINT(" Ipoint", pt);
+ VPRINT(" Inormal", inormal);
+
+ hitp = pp->pt_outhit;
+
+ VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+
+ stp = pp->pt_outseg->seg_stp;
+ RT_HIT_NORMAL(onormal, hitp, stp, &(ap->a_ray), pp->pt_outflip);
+
+ rt_pr_hit(" Out", hitp);
+ VPRINT(" Opoint", pt);
+ VPRINT(" Onormal", onormal);
+
+ /* Now we compute the total mass this element saw while
crossing the region */
+ fastf_t mass, volume, length, avg_density;
+
+ /*Here we query the avg density of a segment */
+ avg_density = segment_density(pp->pt_inhit->hit_point,
pp->pt_outhit->hit_point);
+
+ /* Get the length */
+ length = pp->pt_outhit->hit_dist - pp->pt_inhit->hit_dist;
+ length = length / 10.0; //mm to cm
+
+ /*We also need to know the volume */
+ volume = length * RAY_AREA;
+
+ mass = volume * avg_density;
+ bu_log(" Segment length: %lf, volume crossed: %lf \n", length,
volume);
+ bu_log(" Mass the ray saw through this region: %lf g \n",
mass);
+ }
+ return 1;
}
/**
- * This is a callback routine that is invoked for every ray that
- * entirely misses hitting any geometry. This function is invoked by
- * rt_shootray() if the ray encounters nothing.
- */
+* This is a callback routine that is invoked for every ray that
+* entirely misses hitting any geometry. This function is invoked by
+* rt_shootray() if the ray encounters nothing.
+*/
int
miss(struct application *UNUSED(ap))
{
- bu_log("missed\n");
- return 0;
+ bu_log("missed\n");
+ return 0;
}
-
-/**
- * START HERE
- *
- * This is where it all begins.
- */
int
main(int argc, char **argv)
{
- /* Every application needs one of these. The "application"
- * structure carries information about how the ray-casting should
- * be performed. Defined in the raytrace.h header.
- */
- struct application ap;
+ /* Prepare user point list and let user fill it */
+ BU_GET(user_plist, struct density_point);
+ BU_LIST_INIT(&(user_plist->l));
+ readDensityPoints(user_plist);
- /* The "raytrace instance" structure contains definitions for
- * librt which are specific to the particular model being
- * processed. One copy exists for each model. Defined in
- * the raytrace.h header and is returned by rt_dirbuild().
- */
- static struct rt_i *rtip;
+ /* Prep raytracing and launch one ray */
+ struct application ap;
+ static struct rt_i *rtip;
+ char title[1024] = { 0 };
+ if (argc < 3) {
+ bu_exit(1, "Usage: %s model.g objects...\n", argv[0]);
+ }
- /* optional parameter to rt_dirbuild() that can be used to capture
- * a title if the geometry database has one set.
- */
- char title[1024] = {0};
+ rtip = rt_dirbuild(argv[1], title, sizeof(title));
+ if (rtip == RTI_NULL) {
+ bu_exit(2, "Building the database directory for [%s] FAILED\n",
argv[1]);
+ }
- /* Check for command-line arguments. Make sure we have at least a
- * geometry file and one geometry object on the command line.
- */
- if (argc < 3) {
- bu_exit(1, "Usage: %s model.g objects...\n", argv[0]);
- }
+ if (title[0]) {
+ bu_log("Title:\n%s\n", title);
+ }
- /* Load the specified geometry database (i.e., a ".g" file).
- * rt_dirbuild() returns an "instance" pointer which describes the
- * database to be raytraced. It also gives you back the title
- * string if you provide a buffer. This builds a directory of the
- * geometry (i.e., a table of contents) in the file.
- */
- rtip = rt_dirbuild(argv[1], title, sizeof(title));
- if (rtip == RTI_NULL) {
- bu_exit(2, "Building the database directory for [%s] FAILED\n",
argv[1]);
- }
+ while (argc > 2) {
+ if (rt_gettree(rtip, argv[2]) < 0)
+ bu_log("Loading the geometry for [%s] FAILED\n",
argv[2]);
+ argc--;
+ argv++;
+ }
- /* Display the geometry database title obtained during
- * rt_dirbuild if a title is set.
- */
- if (title[0]) {
- bu_log("Title:\n%s\n", title);
- }
+ rt_prep_parallel(rtip, 1);
+ RT_APPLICATION_INIT(&ap);
+ ap.a_rt_i = rtip;
+ ap.a_onehit = 0;
+ VSET(ap.a_ray.r_pt, 0.0, 0.0, 10000.0);
+ VSET(ap.a_ray.r_dir, 0.0, 0.0, -1.0);
+ VPRINT("Pnt", ap.a_ray.r_pt);
+ VPRINT("Dir", ap.a_ray.r_dir);
+ ap.a_hit = hit;
+ ap.a_miss = miss;
+ (void)rt_shootray(&ap);
- /* Walk the geometry trees. Here you identify any objects in the
- * database that you want included in the ray trace by iterating
- * of the object names that were specified on the command-line.
- */
- while (argc > 2) {
- if (rt_gettree(rtip, argv[2]) < 0)
- bu_log("Loading the geometry for [%s] FAILED\n", argv[2]);
- argc--;
- argv++;
- }
+ return 0;
+}
- /* This next call gets the database ready for ray tracing. This
- * causes some values to be precomputed, sets up space
- * partitioning, computes bounding volumes, etc.
- */
- rt_prep_parallel(rtip, 1);
+/* Returns the density value associated to the indicated point in space,
+for the material in question.
+For now the function receives an initial set of points but we will need
+to talk about how this function gets its initial information to work with.
+*/
+fastf_t old_query_density(point_t query_point) {
- /* initialize all values in application structure to zero */
- RT_APPLICATION_INIT(&ap);
+ /* Density point list */
+ struct density_point * dplist;
+ BU_GET(dplist, struct density_point);
+ BU_LIST_INIT(&(dplist->l));
- /* your application uses the raytrace instance containing the
- * geometry we loaded. this describes what we're shooting at.
- */
- ap.a_rt_i = rtip;
+ /* When no points provided we use this as fallback value */
+ fastf_t DEFAULT_DENSITY = 8.0;
- /* stop at the first point of intersection or shoot all the way
- * through (defaults to 0 to shoot all the way through).
- */
- ap.a_onehit = 0;
+ /* Set up density_vector list */
+ struct density_vector * dvlist;
+ BU_GET(dvlist, struct density_vector);
+ BU_LIST_INIT(&(dvlist->l));
- /* Set the ray start point and direction rt_shootray() uses these
- * two to determine what ray to fire. In this case we simply
- * shoot down the z axis toward the origin from 10 meters away.
- *
- * It's worth nothing that librt assumes units of millimeters.
- * All geometry is stored as millimeters regardless of the units
- * set during editing. There are libbu routines for performing
- * unit conversions if desired.
- */
- VSET(ap.a_ray.r_pt, 0.0, 0.0, 10000.0);
- VSET(ap.a_ray.r_dir, 0.0, 0.0, -1.0);
+ /* This point is the origin for all density_vectors,
+ and holds origin density value of material */
+ struct density_point * origin;
+ BU_GET(origin, struct density_point);
- /* Simple debug printing */
- VPRINT("Pnt", ap.a_ray.r_pt);
- VPRINT("Dir", ap.a_ray.r_dir);
+ /* If point list is empty, no points, we assume homogeneous default
density.
+ I set the origin point to default in this case. */
+ if (BU_LIST_IS_EMPTY(&(dplist->l))) {
+ origin->density = DEFAULT_DENSITY;
+ VSETALL(origin->point, 0.0);
+ /* dvlist will remain empty */
+ }
+ else {
+ /* If there is at least one point, take first point and make it
origin. */
+ BU_LIST_POP(density_point, &(dplist->l), origin);
- /* This is what callback to perform on a hit. */
- ap.a_hit = hit;
+ /* While there are still points, draw a vector from origin to
point
+ and store it in list. */
+ struct density_vector * new_vect;
+ struct density_point * point_iter;
+ while (BU_LIST_WHILE(point_iter, density_point, &(dplist->l))) {
+ BU_GET(new_vect, struct density_vector);
+ VSUB2(new_vect->vector, point_iter->point,
origin->point);
+ new_vect->factor = point_iter->density /
origin->density;
+ BU_LIST_PUSH(&(dvlist->l), new_vect);
+ BU_LIST_DEQUEUE(&(point_iter->l));
+ }
- /* This is what callback to perform on a miss. */
- ap.a_miss = miss;
+ }
- /* Shoot the ray. */
- (void)rt_shootray(&ap);
- /* A real application would probably set up another ray and fire
- * again or do something a lot more complex in the callbacks.
- */
+ /* We draw the vector corresponding to the queried point */
+ vect_t query_vector;
+ VSUB2(query_vector, query_point, origin->point);
- return 0;
+ /* After this, we will always have at least one point, origin.
+ If no more points are present, we use this point's value for
+ homogeneous density. */
+ if (BU_LIST_IS_EMPTY(&(dvlist->l))) {
+ return origin->density;
+ }
+ else {
+
+ /* We now dequeue our vector list and fill in the array of
projection values */
+
+ vect_t paral_projection; //I need this so VPROJECT doesn't
complain?
+ vect_t orth_projection; //Orthogonally projected query_vector
onto density_vectors
+ /* Set up
contribution lists */
+ struct contribution * contrib_list;
+ BU_GET(contrib_list, struct contribution);
+ BU_LIST_INIT(&(contrib_list->l));
+ /* Variables for the loop */
+ struct density_vector * vect_iter;
+ struct contribution * new_contrib;
+ fastf_t contrib, scalar_proj, proportion;
+
+ /* Loop through vectors and build contribution list */
+ while (BU_LIST_WHILE(vect_iter, density_vector, &(dvlist->l))) {
+ VPROJECT(query_vector, vect_iter->vector,
orth_projection, paral_projection);
+ scalar_proj = MAGNITUDE(orth_projection);
+ proportion = scalar_proj / MAGNITUDE(vect_iter->vector);
+ BU_GET(new_contrib, struct contribution);
+ /* In 'scalar_proj' proportion, our contribution value
will be the density of dv,
+ which is original * factor. The 'rest' is original
density */
+ contrib = origin->density * vect_iter->factor *
proportion +
+ origin->density * (1 - proportion);
+ new_contrib->contrib = contrib;
+ new_contrib->dvp = vect_iter;
+ BU_LIST_PUSH(&(contrib_list->l), new_contrib);
+ BU_LIST_DEQUEUE(&(vect_iter->l));
+ }
+
+
+ /* Now go through contributions and do work on it
+ For now take average */
+ struct contribution * contrib_iter;
+ fastf_t acc_contrib = 0.0; //Accumulated contributions (total)
+ int numContribs = bu_list_len(&(contrib_list->l)); //Amount of
entries
+ while (BU_LIST_WHILE(contrib_iter, contribution,
&(contrib_list->l))) {
+ acc_contrib += contrib_iter->contrib;
+ BU_LIST_DEQUEUE(&(contrib_iter->l));
+ };
+
+ /*DEBUG
+ VPRINT("query_vector", query_vector);
+ VPRINT("orth_proj", vector_orth_proj);
+ bu_log("scalar_proj is %lf \n", scalar_proj);
+ bu_log("proportion is %lf \n", proportion);
+ */
+
+ return acc_contrib / numContribs;
+
+ }
+
}
+int compare_dpoints(const void *a, const void *b) {
+ struct projection *x = *(struct projection **)a;
+ struct projection *y = *(struct projection **)b;
+ return x->distance - y->distance;
+}
+
/*
* 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
