Implemented the mentioned changes and fixed projections not being skipped
correctly when they should. Unless my next tests fail, I think this version
is quite stable.
Code attached.
In this test:
https://puu.sh/xn8Hh/0d5c617dfe.png
I added an extra projection point (the third) that should fall onto Z 500
but from far away, obviously its region will not be crossed by our ray, so
the projection should be skipped. So it does.
It also correctly yields the expected result (although in this case the
boundaries do happen to be midpoints) which tells me the code is behaving
at least correctly for this case.
More thorough tests with non-mid-point boundaries will follow.
Mario.
On 29 August 2017 at 18:05, Mario Meissner <mr.rash....@gmail.com> wrote:
> Made a mistake in the last line, should be:
> PR1B = (- DP1PR1^2 + DP2PR2^2 + PR2PR1^2) / (2 * PR2PR1)
>
> Sorry for the confusing nomenclature.
> I also attach the paper version with easier nomenclature.
> I think this will work so I'm gonna start coding it.
> Mario.
>
> On 29 August 2017 at 17:19, Mario Meissner <mr.rash....@gmail.com> wrote:
>
>> Note:
>> Abbreviations:
>> Density Point: DP
>> Projection: PR
>> Boundary: B
>>
>> So I would do the following:
>> Once we have the projections, we can draw two triangles:
>> -DP1, PR1, B
>> -DP2, PR2, B
>>
>> We want to obtain, for example, the distance from projection1 to
>> boundary, so that we can correctly know the percentage of density of
>> region1 that there is between projection1 and projection2.
>>
>> We know that both distances DP1B and DP2B have to be the same, per
>> definition of voronoi tesselation.
>>
>> If we use Pitagoras, we can say that:
>> DP1B^2 = DP1PR1^2 + PR1B^2
>> DP2B^2 = DP2PR2^2 + PR2B^2
>>
>> Knowing DP1B = DP2B we can say that:
>> DP1PR1^2 + PR1B^2 = DP2PR2^2 + PR2B^2
>>
>> We know two of these values and don't know the other two. But the two
>> unknowns have a relation: PR1B + PR2B = PR1-PR2, which we know.
>> We could say that PR2B = PR2PR1 - PR1B
>>
>> So lets substitute:
>> DP1PR1^2 + PR1B^2 = DP2PR2^2 + (PR2PR1 - PR1B)^2
>>
>> Unfold the parenthesis:
>> DP1PR1^2 + PR1B^2 = DP2PR2^2 +
>> PR2PR1^2 + PR1B^2 - 2 * PR2PR1 * PR1B
>>
>> Now we clear the variable we want to compute:
>> PR1B^2 - PR1B^2 + 2 * PR2PR1 * PR1B = - DP1PR1^2 + DP2PR2^2 +
>> PR2PR1^2
>> PR1B^2 = (- DP1PR1^2 + DP2PR2^2 + PR2PR1^2) / (2 * PR2PR1)
>>
>> If I didn't do anything wrong (I'll later double check on paper as well
>> to be sure...), this should be it. We can calculate where the boundary
>> happens between all the projection points we obtained. I had assumed its
>> the middle point by previous discussion we had but now we saw that it can
>> greatly vary.
>>
>> I think it's easier this way (it's a quite simple O(1) operation) than to
>> compute the actual mesh and then evaluate against it. It probably would me
>> more elegant and maybe more efficient but it seems much more complex to me.
>>
>> Let me know what you think, I'll start working on it.
>>
>> Mario.
>>
>>
>> On 29 August 2017 at 15:44, Mario Meissner <mr.rash....@gmail.com> wrote:
>>
>>> Hi Sean!
>>> Back from my trip now.
>>>
>>> The reason why I worked on projections is because you told me to do so:
>>>
>>> Not the distance between points, but distance to their transition
>>>> vectors. From the prior e-mail, the way this should fully generalize can
>>>> be handled by projecting points onto the ray, and projecting any transition
>>>> vectors onto the ray. Then for a given ray, it can fully handle just about
>>>> any transition type.
>>>>
>>>> To see how this will work, it will likely be necessary to redo
>>>> rtexample and get 2-points working without any vectors. You project each
>>>> point onto the ray.
>>>>
>>>> ptA ptB
>>>> | |
>>>> IN v v OUT
>>>> ray o--ptA'----ptB'---o
>>>>
>>>> Density from IN to ptA’ is density(ptA).
>>>> Density from ptB’ to OUT is density(ptB).
>>>> Density from ptA’ to ptB’ is density(ptA) until the midpoint between
>>>> ptA’ and ptB’. Thus the section’s average density would be density(ptA)/2
>>>> + density(ptB)/2But
>>>>
>>>
>>> But I now realize as well that taking the mid-point between projections
>>> does not give us the actual boundary between the polygons.
>>> We could work with the mesh like you mention, or we could do a slight
>>> modification to the current code to compute the actual boundaries while
>>> walking through the points (probably our already available projections,
>>> we'll see). It's probably a simple trigonometry problem, and I'll try to
>>> come up with the operation we need to do to obtain the boundaries in a
>>> segment.
>>> In the attached picture we can see that we need the two blue lines to be
>>> equal in length. Aligning some equations to fit that property may easily
>>> give us the distance from each projected point where the boundary lays.
>>> Unless I'm overseeing something important I think this way would be
>>> sufficient and really easily implementable. The example I'll send this
>>> afternoon may prove me wrong, in which case I'll consider the mesh.
>>> I think I can implement this and leave the code clean before leaving.
>>> Thank you for your feedback!
>>> Mario.
>>>
>>> On 20 August 2017 at 00:03, Christopher Sean Morrison <brl...@mac.com>
>>> wrote:
>>>
>>>>
>>>> > On Aug 17, 2017, at 6:28 AM, Mario Meissner <mr.rash....@gmail.com>
>>>> wrote:
>>>> >
>>>> > So I think now that n-points are working for convex regions, I think
>>>> the next steps would be:
>>>> > • Integrate vectors into existing point system.
>>>> > • Let the user input vectors through the existing input
>>>> interface.
>>>> > • Consider one of the vectors and make it work.
>>>> > • Consider all vectors. N-point and N-vectors working.
>>>> > • Check that all edge cases work correctly (no points,
>>>> no vectors, many of everything, etc.) and clean up code.
>>>> > • Modify implementation so that it properly works in all
>>>> situations (i.e. concave shapes).
>>>>
>>>> Again, I don’t think you should introduce a notion of vectors. For
>>>> now, points will be adequate and present plenty of challenges.
>>>>
>>>> Consider shooting a ray along the X-axis (pnt -1000 0 0; dir 1,0,0)
>>>> with a density point at 0,10,0 and 10,-5,0. Say you enter geometry at
>>>> 1,0,0 and exit at 9,0,0
>>>>
>>>> density
>>>> point A
>>>> o (0,10,0) ___(9,0,0) ray exits
>>>> | /
>>>> | v
>>>> o-> o—x========x-o
>>>> ray (1,0,0) |
>>>> ray o (10,-5,0)
>>>> enters B density point
>>>>
>>>>
>>>>
>>>> With just those two points (A & B), the voronoi density field splits
>>>> halfway between A and B, which has a midpoint above the x-axis and runs
>>>> diagonally across the ray path.
>>>>
>>>> If we simply project, it’ll be the wrong contribution. That’s where I
>>>> was suggesting in the prior e-mail to actually construct the voronoi mesh
>>>> and you’ll get that actual edge between A and B, and calculating where the
>>>> ray intersects that edge becomes trivial.
>>>>
>>>> > • Is it a good idea to compute the distance to the
>>>> previous point and store it into the structure in a loop before starting
>>>> with density evaluation? I think it would clean up the last part of the
>>>> code considerably, but also uses more memory. I like the idea so I probably
>>>> will do so.
>>>>
>>>> We are in no way constrained by memory.
>>>>
>>>> > • Are we really good to go using an array? I decided to do so
>>>> because I needed the sort function, but there might be alternatives I'm not
>>>> aware of.
>>>>
>>>> Shouldn’t need to sort with a voronoi mesh, but nothing wrong with
>>>> using plain arrays + size variables.
>>>>
>>>> > • As mentioned, everything in my code is relative to the inhit
>>>> point. Is this a good approach? For example, projection struct has a vector
>>>> that goes from inhit to the projection point, but I don't store the origin
>>>> of the vector anywhere, so someone using the vector that does not know this
>>>> fact may not know what it's origin is. My thought is that these variables
>>>> will not leave my environment so a simple comment explaining that
>>>> everything is relative should be enough.
>>>>
>>>> If you eliminate vectors and projections, I think this one becomes
>>>> moot. That said, definitely comment the code.
>>>>
>>>> > I will assume that what I propose doing is good unless I hear
>>>> feedback telling me otherwise.
>>>>
>>>> Your e-mail that followed this was golden, as you essentially
>>>> demonstrated that straight-up projections weren’t going to work for even a
>>>> simple 3 point case without changes. That implies you went in the right
>>>> direction, learned something, and now need to adjust accordingly. ;)
>>>>
>>>> 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
>>>>
>>>
>>>
>>
>
Index: src/rt/rtexample.c
===================================================================
--- src/rt/rtexample.c (revisi¾n: 70059)
+++ src/rt/rtexample.c (copia de trabajo)
@@ -67,258 +67,552 @@
#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.
- */
+struct density_vector {
+ struct bu_list l;
+ vect_t vector;
+ point_t origin;
+ fastf_t factor;
+};
+
+/* One possible way of defining vectors is doing it over already defined
points.
+Here we can specify what kind of transition happens between two points. */
+struct transition {
+ struct bu_list l;
+ struct density_point * origin;
+ struct density_point * termination;
+ fastf_t coefficient;
+};
+
+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 contribution {
+ struct bu_list l;
+ fastf_t contrib;
+ struct density_vector * dvp;
+};
+
+struct density_point * user_plist;
+const int MAX_POINTS = 99;
+
+int compare_dpoints(const void *, const void *);
+
+/*
+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
-hit(struct application *ap, struct partition *PartHeadp, struct seg
*UNUSED(segs))
+read_density_points(struct density_point * dplist_head)
{
- /* 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;
+ /* 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? */
+ }
- /* will serve as a pointer to the solid primitive we hit */
- struct soltab *stp;
+ /* Local variables */
+ int count = 0;
+ struct density_point * entry;
+ struct transition * trans;
+ short int origin;
+ short int termination;
+ fastf_t coeff;
+ int correct;
+ int reading = 1;
+ char input[30];
+ struct density_point * entries[99]; //Hardcoded max for now
- /* will contain surface curvature information at the entry */
- struct curvature cur = RT_CURVATURE_INIT_ZERO;
+ /* 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;
- /* will contain our hit point coordinate */
- point_t pt;
+ /* 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;
+ }
- /* will contain normal vector where ray enters geometry */
- vect_t inormal;
+ /* User didn't type continue */
+ else {
- /* will contain normal vector where ray exits geometry */
- vect_t onormal;
+ /* Set up new entry */
+ BU_GET(entry, struct density_point);
- /* 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) {
+ /* 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++;
+ }
- /* 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 );
+ /* Otherwise we dont load it and jump back to
outer loop
+ to print usage again */
+ else {
+ correct = 0;
+ }
+ }
+ } /* Close inner loop */
+ } /* Close outer loop */
- /* entry hit point, so we type less */
- hitp = pp->pt_inhit;
+
- /* 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);
+ 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) {
- /* primitive we encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+ /* 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;
+ }
- /* 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);
+ /* DEBUG */
+ VPRINT("inhit", inhit);
+ VPRINT("outhit", outhit);
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
+ /* 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_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_int, density_point,
&(user_plist->l))) {
+ VSUB2(temp, new_projection->point,
dpoint_iter_int->point);
+ /* If distance to some user point is closer than origin
then no projection */
+ if (dpoint_iter_ext != dpoint_iter_int &&
+ MAGNITUDE(temp) < new_projection->origin_dist) {
+ /* FIXME: Deallocate space we used for
new_projection */
+ is_valid = 0;
+ bu_log("Skipped projection. \n");
+ break;
+ }
+ }
+
+ if (is_valid) {
+
+ /* FIXME: Check that projection actually falls INTO the
segment */
+ 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);
+ }
+
+ }
- /* 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);
+ /* 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);
+ }
- /* 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);
+ qsort(proj_array, array_pointer, sizeof(struct projection *),
compare_dpoints);
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ /* 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);
+ }
- /* 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);
+ /* 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 */
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+ /* 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;
- /* 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);
+ /* 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);
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ /* 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;
+ /* 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);
- /* 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.
- */
+ 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;
+ }
+
+ 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;
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ /* 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;
}
+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));
+ struct transition * trans_list;
+ BU_GET(trans_list, struct transition);
+ BU_LIST_INIT(&(trans_list->l));
+ read_density_points(user_plist, trans_list);
- /* 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;
+ //DEBUG
+ struct density_point * user_plist_iter;
+ struct transition * trans_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);
+ }
+ for (BU_LIST_FOR(trans_iter, transition, &(trans_list->l))) {
+ VPRINT("transition origin: ", trans_iter->origin->point);
+ VPRINT("transition termination: ",
trans_iter->termination->point);
+ printf("With coeff: %lf", trans_iter->coefficient);
- /* optional parameter to rt_dirbuild() that can be used to capture
- * a title if the geometry database has one set.
- */
- char title[1024] = {0};
+ }
+ //END DEBUG
- /* 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]);
- }
- /* 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]);
- }
+ /* 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]);
+ }
- /* Display the geometry database title obtained during
- * rt_dirbuild if a title is set.
- */
- if (title[0]) {
- bu_log("Title:\n%s\n", title);
- }
+ 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]);
+ }
- /* 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++;
- }
+ if (title[0]) {
+ bu_log("Title:\n%s\n", title);
+ }
- /* 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);
+ while (argc > 2) {
+ if (rt_gettree(rtip, argv[2]) < 0)
+ bu_log("Loading the geometry for [%s] FAILED\n",
argv[2]);
+ argc--;
+ argv++;
+ }
- /* initialize all values in application structure to zero */
- RT_APPLICATION_INIT(&ap);
+ 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);
- /* your application uses the raytrace instance containing the
- * geometry we loaded. this describes what we're shooting at.
- */
- ap.a_rt_i = rtip;
+ return 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;
+/* 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) {
- /* 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);
+ /* Density point list */
+ struct density_point * dplist;
+ BU_GET(dplist, struct density_point);
+ BU_LIST_INIT(&(dplist->l));
- /* Simple debug printing */
- VPRINT("Pnt", ap.a_ray.r_pt);
- VPRINT("Dir", ap.a_ray.r_dir);
+ /* When no points provided we use this as fallback value */
+ fastf_t DEFAULT_DENSITY = 8.0;
- /* This is what callback to perform on a hit. */
- ap.a_hit = hit;
+ /* Set up density_vector list */
+ struct density_vector * dvlist;
+ BU_GET(dvlist, struct density_vector);
+ BU_LIST_INIT(&(dvlist->l));
- /* This is what callback to perform on a miss. */
- ap.a_miss = miss;
+ /* 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);
- /* Shoot the ray. */
- (void)rt_shootray(&ap);
+ /* 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);
- /* A real application would probably set up another ray and fire
- * again or do something a lot more complex in the callbacks.
- */
+ /* 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));
+ }
- return 0;
+ }
+
+
+ /* We draw the vector corresponding to the queried point */
+ vect_t query_vector;
+ VSUB2(query_vector, query_point, origin->point);
+
+ /* 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->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
