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,465 @@
#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 contribution {
+ struct bu_list l;
+ fastf_t contrib;
+ struct density_vector * dvp;
+};
- /* will contain surface curvature information at the entry */
- struct curvature cur = RT_CURVATURE_INIT_ZERO;
+struct density_point * user_plist;
- /* will contain our hit point coordinate */
- point_t pt;
- /* 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;
- /* 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) {
+/* Receives a segment and returns average density
+ Uses memory pointer to user point list */
+fastf_t segment_density(point_t inhit, point_t outhit) {
+ /* 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 */
+ /* I dont use this yet but I will */
+ struct density_point * projections;
+ struct density_point * new_projection;
+ BU_GET(projections, struct density_point);
+ BU_LIST_INIT(&(projections->l));
+ vect_t temp_vect;
+ vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+ struct density_point * point_iter;
+ for (BU_LIST_FOR(point_iter, density_point, &(user_plist->l))) {
+ BU_GET(new_projection, struct density_point);
+ VSUB2(temp_vect, point_iter->point, inhit);
+ new_projection->density = point_iter->density;
+ VPROJECT(temp_vect, segment_vect, new_projection->point,
paral_projection);
+ BU_LIST_PUSH(&(projections->l), new_projection);
+ }
- /* entry hit point, so we type less */
- hitp = pp->pt_inhit;
+ /* This is WIP loop to go through projections and accumulate final
density */
+ /* For now however we will asume 0 or 1 points, so this is commented
out */
+ /*
+ struct density_point * proj_iter;
+ fastf_t acc_avg_density = 0;
+ vect_t temp_vect;
+ point_t * prev_pt_pointer = inhit;
+ struct density_point * previous_point;
+ for (BU_LIST_FOR(proj_iter, density_point, &(projections->l))) {
+ 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);
+ }
+ */
- /* 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);
+ /* 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 encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+}
- /* 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);
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
+/* 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) {
+
+ /* Density point list */
+ struct density_point * dplist;
+ BU_GET(dplist, struct density_point);
+ BU_LIST_INIT(&(dplist->l));
- /* 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);
+ /* When no points provided we use this as fallback value */
+ fastf_t DEFAULT_DENSITY = 8.0;
+
+ /* Set up density_vector list */
+ struct density_vector * dvlist;
+ BU_GET(dvlist, struct density_vector);
+ BU_LIST_INIT(&(dvlist->l));
- /* 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);
+ /* 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);
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ /* 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);
- /* 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);
+ /* 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));
+ }
+
+ }
+
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+ /* We draw the vector corresponding to the queried point */
+ vect_t query_vector;
+ VSUB2(query_vector, query_point, origin->point);
- /* 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);
+ /* 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 */
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ 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;
- /* 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.
- */
+ /* 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);
+ */
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ return acc_contrib / numContribs;
+
+ }
+
}
-/**
- * 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;
-}
+/*
+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) {
+ /* 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? */
+ }
-/**
- * START HERE
- *
- * This is where it all begins.
- */
+ /* Local variables */
+
+ /* I count the number of times I */
+ int count = 0;
+ struct density_point *entry;
+ int correct, reading = 1;
+ char input[30];
+
+ /* 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 */
+
+
+
+
+ /**
+ * 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
-main(int argc, char **argv)
+hit(struct application *ap, struct partition *PartHeadp, struct seg
*UNUSED(segs))
{
- /* 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;
+ /* iterating over partitions, this will keep track of the current
+ * partition we're working on.
+ */
+ struct partition *pp;
- /* 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;
+ /* will serve as a pointer for the entry and exit hitpoints */
+ struct hit *hitp;
- /* optional parameter to rt_dirbuild() that can be used to capture
- * a title if the geometry database has one set.
- */
- char title[1024] = {0};
+ /* will serve as a pointer to the solid primitive we hit */
+ struct soltab *stp;
- /* 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]);
- }
+ /* will contain surface curvature information at the entry */
+ struct curvature cur = RT_CURVATURE_INIT_ZERO;
- /* 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]);
- }
+ /* will contain our hit point coordinate */
+ point_t pt;
- /* Display the geometry database title obtained during
- * rt_dirbuild if a title is set.
- */
- if (title[0]) {
- bu_log("Title:\n%s\n", title);
- }
+ /* will contain normal vector where ray enters geometry */
+ vect_t inormal;
- /* 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++;
- }
+ /* will contain normal vector where ray exits geometry */
+ vect_t onormal;
- /* 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);
+ /* Area that the ray covers, in mm^2 */
+ const int RAY_AREA = 4; //2mm height, 2mm width
- /* initialize all values in application structure to zero */
- RT_APPLICATION_INIT(&ap);
+ /* 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) {
- /* your application uses the raytrace instance containing the
- * geometry we loaded. this describes what we're shooting at.
- */
- ap.a_rt_i = rtip;
+ /* 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);
- /* 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;
+ /* entry hit point, so we type less */
+ hitp = pp->pt_inhit;
- /* 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);
+ /* 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);
- /* Simple debug printing */
- VPRINT("Pnt", ap.a_ray.r_pt);
- VPRINT("Dir", ap.a_ray.r_dir);
+ /* primitive we encountered on entry */
+ stp = pp->pt_inseg->seg_stp;
- /* This is what callback to perform on a hit. */
- ap.a_hit = hit;
+ /* 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);
- /* This is what callback to perform on a miss. */
- ap.a_miss = miss;
+ /* print the entry hit point info */
+ rt_pr_hit(" In", hitp);
+ VPRINT(" Ipoint", pt);
+ VPRINT(" Inormal", inormal);
- /* 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.
- */
+ /* exit point, so we type less */
+ hitp = pp->pt_outhit;
- return 0;
+ /* 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);
+
+ /* primitive we exited from */
+ stp = pp->pt_outseg->seg_stp;
+
+ /* 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);
+
+ /* print the exit hit point info */
+ 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;
+
+ /* FIXME: 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
+
+ /*FIXME: 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.
+*/
+int
+miss(struct application *UNUSED(ap))
+{
+ bu_log("missed\n");
+ return 0;
+}
+
+int
+main(int argc, char **argv)
+{
+ /* 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);
+
+ /* 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]);
+ }
+
+ 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]);
+ }
+
+ if (title[0]) {
+ bu_log("Title:\n%s\n", title);
+ }
+
+ while (argc > 2) {
+ if (rt_gettree(rtip, argv[2]) < 0)
+ bu_log("Loading the geometry for [%s] FAILED\n",
argv[2]);
+ argc--;
+ argv++;
+ }
+
+ 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);
+
+ return 0;
+}
+
/*
* 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
