Hy Sean!
Attached goes code that includes the user input update to support
transitions.
Also contains work on the fix I talked about in my last email. I've tested
the current state of the code but it doesn't seem to skip projections
correctly yet. However my intentions are reflected in code so that you can
see what I want to do. It's either a stupid mistake in the code or in my
tests, or some conceptual error I've overseen.
I won't be available from tomorrow until the 28th, and from the 1st to the
3rd. I'm going on an 8 day trip and quickly after I'm moving to Germany to
start one academic year abroad. I hope to be able to set myself up there
quickly to continue working on this.
Cheers!
Mario.
On 18 August 2017 at 21:44, Mario Meissner <mr.rash....@gmail.com> wrote:
> I found a flaw in my current code.
> I considered all points to be valid and projected them to calculate the
> final density. However, not always do we need all points to be projected.
> In the attached picture we can see that the diagonal red line doesn't use
> the black region so we shouldn't be projecting it.
> I think one way to solve it is to project only those who's distance to the
> original point is closer than to any other point. We will need to loop
> through all points for every point though, which is at least O(n^2).
> I also ended up using absolute points to compute those distances so ignore
> my talk about stuff being relative.
>
> Currently working on the fix so I'll send new code for that tomorrow.
> Just wanted to let you know that I found the flaw and am working on it.
> Also sorry for the avalanche of emails, but I prefer mailing too much than
> too little.
>
> Mario.
>
>
> 2017-08-18 19:02 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>
>> Sean,
>> So I've given the transition concept a deep thought and I think I may
>> have found one way it could work.
>> We define a transition object that links to two points and has a float
>> value stored. This value indicates the rate of change between the two
>> points. How does it indicate it in a way that is useful for us to perform
>> computations on it?
>> Let's consider a transition situation. One point higher (value) than the
>> other, and separated by a certain distance. Lets now take the value of the
>> function right in the middle between those two points. If the transition
>> were linear, then that value should exactly be the mean value of the two
>> points. If it's not, then some sort of quadratic rate of change is
>> happening (let's assume for now that in between the points the values can't
>> grow above or below the higher and lower points respectively, and let's for
>> now only give support to parabolas). The value we store in the structure
>> could somehow represent where the value in the middle of the points sits.
>> For example a value of 1 could represent that it's exactly the mean value,
>> meaning linear transition. A value closer to 0 would mean below the mean
>> value, and above 1 greater. Once we have this new point (and we might want
>> to compute some more points as well using that value) we can safely compute
>> the numerical integral using Newton-Cotes or Simpson. Now, before starting
>> to work on this I would highly recommend to discuss it a bit more.
>>
>> I'm also a bit worried about what might happen if there are
>> non-transition points in between a transition. Should it be an illegal
>> state? It's almost impossible to avoid it, since depending from where you
>> shoot the rays, a projection could always potentially land inside a
>> transition vector.
>>
>> Id also love to see something like the Voronoi diagram you showed me but
>> with support for transitions. I'm still not sure how it would work. If I
>> look at the diagrams and try to imagine how it would look like with one
>> transition on it, I can't come up with anything logical. How does the
>> transition interact with non-transition points? Between non-transition
>> points there is a line indicating when we change from one value to another.
>> What happens between a non-transition point, and a transition point? Is
>> there a line, or not? There should be because there is no transition
>> between that point and the vector, but ... how and where? Is there any
>> example available? I tried looking for it but didn't find anything.
>>
>> I hope you'll find some time to discuss this soon.
>> Mario.
>>
>> 2017-08-18 17:40 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>>
>>> Hello everyone.
>>> Quick question regarding bu_lists.
>>>
>>> When I BU_LIST_PUSH to the tail of the list, and then BU_LIST_FOR
>>> through them, why do I get the last element first? I expected PUSH to put
>>> the element behind all other elements, and then FOR to roll through them
>>> from front to back, showing me the elements in the order I pushed them.
>>>
>>> Do I need to use APPEND and keep a pointer to my last added element to
>>> get the behaviour I need? Why?
>>>
>>> Mario.
>>>
>>> 2017-08-18 13:38 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>>>
>>>> Good morning / afternoon !
>>>>
>>>> I've started considering the inclusion of vectors into the current
>>>> code, and realized that linear vectors may not make much sense in terms of
>>>> new value.
>>>> Yes, with a vector we can define a smooth transition between points
>>>> instead of an abrupt change in the middle of two points. However, the final
>>>> result will be the same. The points meaning abrupt or linear transition,
>>>> after all, is only relevant if we want to know the density at a specific
>>>> point in space, which we don't really use here currently (because I find
>>>> that it's much more comfortable to work with segments directly).
>>>>
>>>> If we only needed to define a density model for the purpose of
>>>> returning the correct value with the segment_density() function, then
>>>> linear vectors are superfluous. We can perfectly just define the two points
>>>> that make up the vector (origin and termination) and compute the result
>>>> using the current code. Results should be exactly the same. Voronoi
>>>> distribution and linear transitions are equivalent in this regard. Now, if
>>>> this model is going to be used by some external agent as well then we may
>>>> need to properly define it.
>>>>
>>>> My question here is, how should I model these vectors if at all,
>>>> considering that they won't be of any use in my current code? I thought one
>>>> way we could implement it is with a struct transition (probably rendering
>>>> obsolete density_vector) that references two points, and contains a
>>>> coefficient value that tells us how 'fast' the transition is (with 1
>>>> meaning linear). This way we include linear transitions in an elegant way
>>>> but give freedom to define any kind of them. Implementing them may be a
>>>> challenge but it seems like our best bet. I'll start working on this. If
>>>> you think this is a bad idea then early feedback is welcome so I don't step
>>>> too deep into the mud.
>>>>
>>>> Mario.
>>>>
>>>> 2017-08-17 12:28 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>>>>
>>>>> 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).
>>>>>
>>>>> I'm also noticing that the current code is getting quite messy so I
>>>>> think its time to make a cleanup before continuing. A few aspects I want
>>>>> to
>>>>> work on:
>>>>>
>>>>> - What data should the projection structure contain?
>>>>> - Do we need the absolute position (point) of the projection? I
>>>>> think no, because in my code I am assuming everything is relative
>>>>> to inhit,
>>>>> so currently the point and its computation are commented out.
>>>>> - 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.
>>>>> - 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.
>>>>> - 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.
>>>>>
>>>>> I will assume that what I propose doing is good unless I hear feedback
>>>>> telling me otherwise.
>>>>> Today I will work on a bit of cleanup and then the first point of my
>>>>> proposed steps.
>>>>> Attached goes code with minor cleanup changes over yesterday. If
>>>>> possible, use this for review.
>>>>> Happy coding!
>>>>> Mario.
>>>>>
>>>>> 2017-08-15 22:28 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>>>>>
>>>>>> 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,583 @@
#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,
+ struct transition * translist_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_LIST_IS_INITIALIZED(&(translist_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 'continue' to proceed defining the transition
between points. \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, "continue\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;
+ reading = 1;
+ /* Outer loop. We come back here whenever user inputs wrong format
+ and until he types exit */
+ while (reading) {
- /* 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);
+ bu_log("Please indicate which of the above points form a
special transition.\n");
+ bu_log("Refer to your previously defined points with numbers
(1,2,3...) \n");
+ bu_log("Use the format X Y C where X and Y are references to
your points");
+ bu_log(" and C is the coefficient value of the transition. \n");
+ bu_log("Type 'exit' to finish. \n");
+ correct = 1;
- /* primitive we encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+ /* Now we read the transition values between defined points. */
+ /* Inner loop. While user inputs transitions in correct format
we loop here */
+ while (correct) {
+ bu_fgets(input, 30, stdin);
+ /* If user typed exit we need to exit both loops */
+ if (!strcmp(input, "exit\n")) {
+ bu_log("Finished user input. \n");
+ correct = 0;
+ reading = 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);
+ /* User didn't type exit */
+ else {
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
+ /* Prepare transition */
+ BU_GET(trans, struct transition);
- /* 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);
+ /* If the input is in correct format we load it
and continue */
+ if (sscanf(input, "%hi %hi %lf", &origin,
&termination, &coeff) == 3) {
+ if (origin <= count && termination <=
count && coeff > 0) {
+ trans->origin = entries[origin];
+ trans->termination =
entries[termination];
+ trans->coefficient = coeff;
+
BU_LIST_PUSH(&(translist_head->l), &(trans->l));
+ }
+ else {
+ bu_log("Invalid point
references or coefficient value have been given.\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);
+ /* Otherwise we dont load it and jump back to
outer loop
+ to print usage again */
+ else {
+ correct = 0;
+ }
+ }
+ } /* Close inner loop */
+ } /* Close outer loop */
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ 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) {
- /* 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 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;
+ }
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+ /* DEBUG */
+ VPRINT("inhit", inhit);
+ VPRINT("outhit", outhit);
- /* 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);
+ /* 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 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);
+ /* 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");
+ continue;
+ }
+ }
+
+ if (is_valid) {
+ 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);
+ /* 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);
+ }
+
+ }
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ /* 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);
+ }
- /* 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.
- */
+ qsort(proj_array, array_pointer, sizeof(struct projection *),
compare_dpoints);
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ /* 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 * proj_iter;
+ fastf_t acc_avg_density = 0.0;
+ struct projection * prev_proj;
+
+ /* 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->inhit_dist
/ segment_len;
+ /* DEBUG */
+ bu_log("First projection. Proportion %f, added %lf to accumulator.\n",
+ proj_iter->inhit_dist / segment_len,
proj_iter->original->density *
+ proj_iter->inhit_dist / segment_len);
+ prev_proj = proj_iter;
+
+ /* Loop for remaining projections. 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++) {
+ proj_iter = proj_array[i];
+ acc_avg_density += prev_proj->original->density
+ * proj_iter->prev_dist / (2 * segment_len);
+ acc_avg_density += proj_iter->original->density
+ * proj_iter->prev_dist / (2 * segment_len);
+ prev_proj = proj_iter;
+ }
+
+ /* Last portion, from the last point to the outhit point has density
value
+ of this last point */
+ acc_avg_density += proj_iter->original->density * proj_iter->prev_dist
/ segment_len;
+
+ /* DEBUG */
+ bu_log("Last projection. Proportion %f, added %lf to accumulator.\n",
+ proj_iter->prev_dist / segment_len,
proj_iter->original->density *
+ proj_iter->prev_dist / 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
