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.
2017-08-11 23:28 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
> 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!!
>
> 2017-08-11 23:26 GMT+02:00 Mario Meissner <mr.rash....@gmail.com>:
>
>>
>> ---------- Forwarded message ----------
>> From: MarioPC <mr.rash....@gmail.com>
>> Date: 2017-08-10 20:33 GMT+02:00
>> Subject: RE: [brlcad-devel] Example of heterogeneous density +
>> masscalculation
>> To: Mario Meissner <mr.rash....@gmail.com>
>>
>>
>> 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.
>>
>>
>>
>> Mario.
>>
>>
>>
>> *De: *Mario Meissner <mr.rash....@gmail.com>
>> *Enviado: *miércoles, 9 de agosto de 2017 23:04
>> *Para: *BRL-CAD Developer Mailing List
>> <brlcad-devel@lists.sourceforge.net>
>> *Asunto: *Re: [brlcad-devel] Example of heterogeneous density +
>> masscalculation
>>
>>
>>
>> Hi Sean!
>>
>> Thank you a lot for the details. The visualization of Voronoi points was
>> really helpful. I think I understand your vision and will work on
>> implementing the first steps as soon as the stuff I mention below is done.
>>
>> I'm running into some trouble with my code. I need to call
>> readDensityValues() (user input through stdin) but I don't know exactly
>> where I should place the call.
>>
>> Placing it into view_init, view_2init, application_init or view_setup
>> breaks it: buffer doesn't behave correctly. My exit command isn't
>> recognized anymore, and when I Crtl-C, it runs the while loop two more
>> times (printing usage) before shutting off. It feels like there are many
>> threads running at once and interfering with each other. Might this be the
>> case?
>>
>> I'm afraid I cannot deliver functioning code today, but if we solve this
>> problem I will be able to do so tomorrow. For now here goes the not-working
>> state.
>>
>>
>>
>> Next steps:
>> -> store the points into a file (instead of sharing a memory pointer)
>>
>> -> make the point input procedure a separate program, or ask the user at
>> the start of rtweight execution if he wants to create/override the point
>> list.
>>
>> Maybe this can solve the problem mentioned in the above paragraphs.
>> How do you think we should do it?
>>
>>
>>
>> I also want to mention that the during the next three weeks I will have a
>> little bit more restricted available time (as I mentioned in my
>> application). I will steadily advance, but it might feel a little bit
>> slower than up to now.
>>
>>
>>
>> PD. What you mean with 'patch' is the result of running svn diff into a
>> file? If so, here it goes!
>>
>> Mario.
>>
>>
>>
>>
>>
>> 2017-08-09 7:58 GMT+02:00 Christopher Sean Morrison <brl...@mac.com>:
>>
>>
>>
>> Fair enough. My next goal will be to make our initial minimal feature
>> working. In this case this means having a way to let the user define points
>> that will be used for query_density(), instead of hardcoding them. I guess
>> for now I will use the .density file as I have no idea how the whole
>> "material objects" thing should work, yet.
>>
>>
>>
>> If you have a specific idea/plan in mind, I’d love to hear it. However,
>> I can’t really envision a good way to extend the tabular .density file
>> format with the information you need without introducing a problem.
>>
>>
>>
>> Two options that come to mind would be to modify rtweight to take a
>> secondary input file (temporary) or implement a material object stub that
>> has what you need (more complicated, but in the right direction) [1]. For
>> time, I think the second option is the way to go given the time remaining
>> and that’ll keep you directly working with the rtweight sources.
>>
>>
>>
>> [1] http://brlcad.org/wiki/Material_and_Shader_Objects
>>
>>
>>
>> I guess it's up to you to decide what I should focus on. I could work on
>> the basics of this new approach so that future work on it can then build on
>> it and eventually finish it, but that would mean quite a change to what our
>> initial goal for the socis period was. I'm totally fine with either way.
>>
>>
>>
>> If we knowingly plan to extend the current materialID system to later
>> point to an object, then you’ll just have to make sure whatever you read in
>> from file fits that mental model. Directly modifying rtweight to take a
>> second specification file will more easily allow for testing N-point cases,
>> so we’ll just have to be careful it can easily translate into an object.
>>
>>
>>
>> I think this [checking if there are overlapping points] can also be
>> handled by changing the way vectors are being defined (or perhaps what they
>> mean). Instead of having a concept of origins and detachment of point
>> densities, let density_vector actually refer to a density_point and let the
>> other struct fields encode the nature of the contribution (effectively
>> combining your struct combination data into one concept).
>>
>>
>>
>> I'm afraid I'm not following you here. Could you detail this out a bit
>> more? The referenced density_point would be the origin of the vector? (if
>> not, what is it?).
>>
>>
>>
>> No, at least not necessarily. Density points would be simply a density
>> specification at a specific point in space. How that value is interpreted
>> or used nearby would be a separate concept (and handled in code in a
>> separate struct). In fact, it becomes possible to define reasonable
>> behavior without any vector information.
>>
>>
>>
>> Consider two side by side boxes (with a gap between them) that are in the
>> same region. Under the current system, there’s a materialID=3 which maps
>> to a density of 7.75 g/cm^3. If I specified a density point somewhere in
>> one of the boxes as being 8.05 g/cm^3, I would intuitively expect that to
>> be the density of that box (and only that box). No falloff, no smooth
>> transition to 7.75, no override of the universe to 8.05.
>>
>>
>>
>> How to achieve that in code is going to be really tricky for concave
>> shapes, but I think it’s the simplest and most natural behavior a user
>> would expect. (keep the simple easy)
>>
>>
>>
>> The 'nature' of the contribution would be something like rate: linear,
>> quadratic, etc.?
>>
>>
>>
>> It’s only once we add a second point or define a transition vector (which
>> implicitly entails two points) that we have to care about this, but
>> essentially yes.
>>
>>
>>
>> And how does this change help us check if there are overlapping points?
>>
>>
>>
>> So again consider defining points but not defining vectors. If we have
>> one box and define two density points in the box, that is the density at
>> those points. Defining the same point as having two separate densities
>> would be a definition error. Separate them by any meaningful distance and
>> we’ve effectively defined a density “binds" that maps well to a Voronoi
>> tessellation of space, like this:
>>
>>
>>
>> http://alexbeutel.com/webgl/voronoi.html
>>
>>
>>
>> With that conceptualization, adding in vectors from one point to another
>> merely changes the space definition to being a continuous smooth transition
>> between pairs of points. Make sense? That becomes a fully generalize
>> solution that should capture nearly any material definition while keeping
>> the simple easy, and complex possible.
>>
>>
>>
>> From the rest of the answers in your emails I think I got the idea that
>> you want to let the user specify both points and vectors as if they were
>> two means of describing density. Then a vector should be something like a
>> superset of the point, where we can specify some more information about it
>> to
>>
>> describe the continuous transition between the two points the vector
>> spans.
>>
>>
>>
>> I think you got it. More specifically, vectors augment/extend the
>> limited info of a density point, describing how it transitions out from a
>> point.
>>
>>
>>
>> Don't we need to reference two density_points then?
>>
>>
>>
>> Either as two points (origin pt1 to destination pt2 with an implicit dir)
>> or as one point and a direction vector, not unitized (so adding pt1+dir
>> gives pt2). Either should work fine.
>>
>>
>>
>> Doesn't this method depend too much on the direction from which the rays
>> come? For example, ptA could be really far from 'IN' if you shoot from
>> above but really close to it if you shoot from the side. If ptA was the
>> heaviest point, then shooting from above will turn a huge part of the
>> material into really heavy, and shooting from the side should weight much
>> less. I'm not sure if I'm talking nonsense here, I'll check with an actual
>> example to see it working, but this was my first thought when I saw it.
>>
>>
>>
>> Hopefully the Voronoi webgl demo helps clear that up. There was only
>> weighting when we were assuming linear interpolation or vectors to/from all
>> points. Using the Voronoi / halfspace approach is direction invariant.
>>
>>
>>
>> For now I'll make my code fully 'complete', as you mention. Then I will
>> get rid of my 'origins' and make projections onto the shotline (but this
>> will require some more discussions as I'm not 100% sure how this would all
>> work and fit together, especially points with no vectors [the question
>> above]). At some point in the future we will need to take the distance to
>> the shotline into consideration to make a 'fair' average of contributions,
>> as you suggested.
>>
>>
>>
>> It’s also okay if you need to move back to rtexample to prove the concept
>> first. It’s when the feature is introduced into a production facility that
>> it should be done “complete”, one bit at a time.
>>
>>
>>
>> For what it’s worth, projecting onto the shotline is probably only going
>> to work on convex geometry. For concave shapes, it’s a bit more complex to
>> know where a given density point applies. It’s certainly resolvable with
>> the neighboring ray information (because they represent a
>> quasi-voxelization of space).
>>
>>
>>
>> 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,438 @@
#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;
+ 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 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) {
+ /* Segment vector (for projections) */
+ vect_t segment_vect;
+ VSUB2(segment_vect, outhit, inhit);
+ fastf_t segment_mag = MAGNITUDE(segment_vect);
- /* 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 );
+ /* -- Compute projection list -- */
- /* entry hit point, so we type less */
- hitp = pp->pt_inhit;
+ /* Loop prep */
+ struct projection * new_projection;
+ /*struct density_point * projections;
+ BU_GET(projections, struct density_point);
+ BU_LIST_INIT(&(projections->l));*/
+ vect_t dpoint_vect;
+ vect_t proj_vect;
+ vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+ struct density_point * point_iter;
+ int array_pointer = 0;
+ struct projection ** proj_array;
+ bu_calloc(MAX_POINTS, sizeof(struct projection *), proj_array);
+
- /* 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);
+ /* Loop through points and project them onto segment */
+ for (BU_LIST_FOR(point_iter, density_point, &(user_plist->l))) {
+ BU_GET(new_projection, struct density_point);
+ VSUB2(dpoint_vect, point_iter->point, inhit);
+ new_projection->original = point_iter;
+ VPROJECT(dpoint_vect, segment_vect, new_projection->point,
paral_projection);
+ /* FIXME: Check that projection actually falls INTO the segment
*/
+ /*BU_LIST_PUSH(&(projections->l), new_projection);*/
+ VSUB2(proj_vect, new_projection->point, inhit);
+ new_projection->distance = MAGNITUDE(proj_vect);
+ proj_array[array_pointer] = new_projection;
+ array_pointer++;
+ }
- /* primitive we encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+ /* Sort projection list based on distance to inhit */
+ /* I wonder if this will work just like this? */
+ qsort(proj_array, array_pointer, sizeof(struct projection *),
compare_dpoints);
- /* 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);
+ /* -- Go through projections and calculate (with precision) average
density -- */
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
+ /* Loop preparation */
+ struct density_point * proj_iter;
+ fastf_t acc_avg_density = 0;
+ vect_t temp_vect;
+ point_t * prev_pt_pointer = inhit;
- /* 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);
+ /* First projection point, we take everything from inhit to this point
+ as having density of this point */
+ //Dequeue first point
- /* 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);
+ /* Loop for the rest. For each point, take half the value from last
point,
+ and half the value from current point */
+ struct density_point * previous_point;
+ for (int i = 0; i < array_pointer; i++) {
+ VSUB2(temp_vect, proj_iter->point, *prev_pt_pointer);
+ acc_avg_density += proj_iter->density * MAGNITUDE(temp_vect) /
segment_mag;
+ prev_pt_pointer = &(proj_iter->point);
+ }
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ /* Last portion, from the last point to the outhit point has density
value
+ of this last point */
+ //DOSTUFF
- /* 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;
+ }
+ /* Else (for now) there is one point, so homogeneous density */
+ else {
+ /* This code is really silly and will dissappear once the above
loops work */
+ struct density_point * point;
+ BU_LIST_POP(density_point, &(user_plist->l), point);
+ return point->density;
+ }
+}
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+/*
+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) {
- /* 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);
+ /* 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? */
+ }
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ /* Local variables */
- /* 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.
- */
+ /* I count the number of times I */
+ int count = 0;
+ struct density_point *entry;
+ int correct, reading = 1;
+ char input[30];
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ /* 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;
+ int c;
+ /* Inner loop. While user inputs points in correct format we
loop here */
+ while (correct) {
+
+ bu_fgets(input, 30, stdin);
+ //gets(input);
+
+ //while ((c = getchar()) != '\n' && c != EOF) {}
//Clean buffer
+
+ /* If user typed exit we need to exit both loops */
+ if (!strcmp(input, "exit\n")) {
+ bu_log("Exiting...");
+ 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, pp->pt_outhit);
+
+ /* Get the length */
+ length = pp->pt_outhit->hit_dist - pp->pt_inhit->hit_dist;
+ length = length / 10; //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
