Peopleee!
Here goes today's work. I'm proud to announce the example I modeled a few
days ago now works correctly!!
Code still contains my debug prints and lots of random 'puts' to know where
and why stuff crashes. I feel I will need it in the future so for now it
remains.
I probably should keep looking for flaws by setting up other examples, and
while doing so I should also clean it up a bit since the code grew over 600
lines now, quite difficult to read and understand.
Screenshot:
https://puu.sh/xuIXk/9626ce2fb6.png
Mario.
On 8 September 2017 at 19:56, Mario Meissner <mr.rash....@gmail.com> wrote:
> Hi Daniel!
>
> Thank you a lot for your reply. I was focused on generating the bisector
> to intersect two lines, but intersecting the actual plane is much simpler.
> I made a little function to intersect a plane with a line (partially taken
> from this stackoverflow post
> <https://stackoverflow.com/questions/5666222/3d-line-plane-intersection>,
> so credits to user ideasman42
> <https://stackoverflow.com/users/432509/ideasman42>). I am now finishing
> of the 'patch' to my code that is supposed to now correctly skip
> unnecessary points and create a list with those who's regions we will
> actually cross. Code probably tomorrow.
>
> Mario.
>
> On 8 September 2017 at 09:35, Daniel Roßberg <danielmrossb...@gmail.com>
> wrote:
>
>> Hi Mario,
>>
>> In 3D space, rectangular to a vector you have a plane and the vector is a
>> normal vector of this plane. In this plane are all vectors which are
>> perpendicular to your first vector, i.e. the vectors which products with
>> your first vector are 0. This gives you a first equation with two
>> variables.
>>
>> Next, you want to intersect the plane with a ray. I.e., one vector from
>> your plan has to be a multiple of the ray which gives you another equation.
>>
>> However, there should be examples in the BRL-CAD source code for this.
>> E.g., the ray-trace of the half primitive is very similar to your problem.
>> There, the intersection of a plane, given by a normal vector and a point in
>> it, and a ray, given by a start point and a direction, will be computed.
>>
>> I hope this helps.
>>
>>
>> Regards,
>> Daniel
>>
>>
>>
>> I'm facing a wall. I need to find a perpendicular vector with the
>> condition that it will intersect a third vector.
>> Basically, I have the ray and a vector drawn from my current point to the
>> previous point.
>> I want the perpendicular bisector of that vector to cross my ray, so that
>> I can obtain the intersection between the bisector and the ray.
>>
>> For that I need two things:
>> 1) Obtain the perpendicular vector. in 3D there are infinite options, but
>> only one will cross the ray. How do I get the one I need?
>> 2) Obtain the intersection point.
>>
>> I browsed through vmath.h but I didn't find anything that could help me.
>> Do you know any macro that could be useful here?
>> Otherwise I will need to make some of my own for this! If so, any ideas
>> on how to approach it?
>>
>> Mario.
>>
>>
>> ------------------------------------------------------------
>> ------------------
>> 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: 70248)
+++ src/rt/rtexample.c (copia de trabajo)
@@ -67,258 +67,622 @@
#include <math.h>
#include <string.h>
#include <stdio.h>
+#include <string.h>
#include "vmath.h" /* vector math macros */
#include "raytrace.h" /* librt interface definitions */
+fastf_t MARGIN = 0.0001;
+struct density_point {
+ struct bu_list l;
+ point_t point;
+ fastf_t density;
+};
+struct density_vector {
+ struct bu_list l;
+ vect_t vector;
+ point_t origin;
+ fastf_t factor;
+};
-/**
- * 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.
- */
+/* One possible way of defining vectors is doing it over already defined
points.
+Here we can specify what kind of transition happens between two points. */
+struct transition {
+ struct bu_list l;
+ struct density_point * origin;
+ struct density_point * termination;
+ fastf_t coefficient;
+};
+
+struct projection {
+ //If we store the density here again we might not need the original
+ //fastf_t density;
+ struct density_point * original;
+ point_t point;
+ vect_t vect;
+ fastf_t inhit_dist;
+ fastf_t prev_dist;
+ fastf_t origin_dist;
+};
+
+struct contribution {
+ struct bu_list l;
+ fastf_t contrib;
+ struct density_vector * dvp;
+};
+
+struct density_point * user_plist;
+const int MAX_POINTS = 99;
+
+int compare_dpoints(const void *, const void *);
+
+/*
+Reads density points from stdin and stores them into density_point structs
+Receives the head of a bu_list of density_point structs
+and appends all points generated during the reading session to the list
+Returns the number of successfully read points
+*/
int
-hit(struct application *ap, struct partition *PartHeadp, struct seg
*UNUSED(segs))
+read_density_points(struct density_point * dplist_head)
{
- /* iterating over partitions, this will keep track of the current
- * partition we're working on.
- */
- struct partition *pp;
- /* will serve as a pointer for the entry and exit hitpoints */
- struct hit *hitp;
+ /* If we are given no proper bu_lists ... */
+ if (!BU_LIST_IS_INITIALIZED(&(dplist_head->l))) {
+ bu_log("No valid heads of bu_list has been provided.\n");
+ /* Maybe also check if bu_list is of correct type? */
+ }
- /* will serve as a pointer to the solid primitive we hit */
- struct soltab *stp;
+ /* Local variables */
+ int count = 0;
+ struct density_point * entry;
+ struct transition * trans;
+ short int origin;
+ short int termination;
+ fastf_t coeff;
+ int correct;
+ int reading = 1;
+ char input[30];
+ struct density_point * entries[99]; //Hardcoded max for now
- /* will contain surface curvature information at the entry */
- struct curvature cur = RT_CURVATURE_INIT_ZERO;
+ /* Outer loop. We come back here whenever user inputs wrong format
+ and until he types exit */
+ while (reading) {
+ bu_log("Please insert density points in the format:");
+ bu_log(" X Y Z density_value(in g / cm3). \n");
+ bu_log("Type 'exit' to finish. \n");
+ correct = 1;
- /* will contain our hit point coordinate */
- point_t pt;
+ /* Inner loop. While user inputs points in correct format we
loop here */
+ while (correct) {
+ bu_fgets(input, 30, stdin);
+ /* If user typed continue we need to exit both loops */
+ if (!strcmp(input, "exit\n")) {
+ bu_log("Finished inserting points.\n");
+ correct = 0;
+ reading = 0;
+ }
- /* will contain normal vector where ray enters geometry */
- vect_t inormal;
+ /* User didn't type continue */
+ else {
- /* will contain normal vector where ray exits geometry */
- vect_t onormal;
+ /* Set up new entry */
+ BU_GET(entry, struct density_point);
- /* iterate over each partition until we get back to the head.
- * each partition corresponds to a specific homogeneous region of
- * material.
- */
- for (pp=PartHeadp->pt_forw; pp != PartHeadp; pp = pp->pt_forw) {
+ /* If the input is in correct format we load it
and continue */
+ if (sscanf(input, "%lf %lf %lf %lf",
&entry->point[0],
+ &entry->point[1], &entry->point[2],
&entry->density) == 4) {
+ BU_LIST_PUSH(&(dplist_head->l),
&(entry->l));
+ entries[count] = entry;
+ count++;
+ }
- /* print the name of the region we hit as well as the name of
- * the primitives encountered on entry and exit.
- */
- bu_log("\n--- Hit region %s (in %s, out %s)\n",
- pp->pt_regionp->reg_name,
- pp->pt_inseg->seg_stp->st_name,
- pp->pt_outseg->seg_stp->st_name );
+ /* Otherwise we dont load it and jump back to
outer loop
+ to print usage again */
+ else {
+ correct = 0;
+ }
+ }
+ } /* Close inner loop */
+ } /* Close outer loop */
- /* entry hit point, so we type less */
- hitp = pp->pt_inhit;
+
- /* construct the actual (entry) hit-point from the ray and the
- * distance to the intersection point (i.e., the 't' value).
- */
- VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+ return count;
+}
- /* primitive we encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+point_t * intersection(point_t p0, point_t p1, point_t p_co, vect_t p_no) {
+ fastf_t epsilon = 1e-6;
+ fastf_t fac;
+ vect_t u, w;
+ point_t intersect;
+ VSUB2(u, p1, p0);
+ fastf_t dot = VDOT(p_no, u);
+ if (abs(dot) > epsilon) {
+ VSUB2(w, p0, p_co);
+ fac = -VDOT(p_no, w) / dot;
+ VSCALE(u, u, fac);
+ VADD2(intersect, p0, u);
+ return &intersect;
+ }
+ else {
+ return NULL;
+ }
+}
- /* 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);
+/* Receives a segment and returns average density
+ Uses memory pointer to user point list */
+fastf_t segment_density(point_t inhit, point_t outhit) {
- /* 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 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;
+ }
- /* 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);
+ /* DEBUG */
+ VPRINT("inhit", inhit);
+ VPRINT("outhit", outhit);
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ /* Segment vector (for projections) */
+ vect_t segment_vect;
+ VSUB2(segment_vect, outhit, inhit);
+ fastf_t total_segment_len = MAGNITUDE(segment_vect);
+
+ /* Projections Loop prep */
+ struct projection * new_projection;
+ struct density_point * dpoint_iter_ext; //External
+ struct density_point * dpoint_iter_mid; //Middle
+ 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];
- /* 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);
+ /* Loop through points and project them onto segment if they need to be
*/
+ for (BU_LIST_FOR(dpoint_iter_ext, density_point, &(user_plist->l))) {
+ is_valid = 1;
+ BU_GET(new_projection, struct projection);
+ /* Create userpoint vector */
+ VSUB2(dpoint_vect, dpoint_iter_ext->point, inhit);
+ new_projection->original = dpoint_iter_ext;
+ /* Project vector onto segment */
+ VPROJECT(dpoint_vect, segment_vect,
+ new_projection->vect, paral_projection);
+ /* Compute absolute point of projection */
+ VADD2(new_projection->point, new_projection->vect, inhit);
+ new_projection->inhit_dist = MAGNITUDE(new_projection->vect);
+ /* Compute distance to original */
+ VSUB2(proj_to_original_vect, dpoint_iter_ext->point,
new_projection->point);
+ new_projection->origin_dist = MAGNITUDE(proj_to_original_vect);
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+ /* Check if we need to project or not based on distance to all
other points */
+ for (BU_LIST_FOR(dpoint_iter_mid, density_point,
&(user_plist->l))) {
+ VSUB2(temp, new_projection->point,
dpoint_iter_mid->point);
+ /* If distance to some user point is shorter than to
the origin then we need to check boundaries */
+ if (dpoint_iter_ext != dpoint_iter_mid &&
+ MAGNITUDE(temp) < new_projection->origin_dist) {
- /* 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 here it means projection has
closer neighbor than original,
+ we now check if intersection falls inside
segment and has closer neighbor.
+ Also it should not be our first point. */
+ if (array_pointer > 0) {
+ puts("array pointer bigger than zero");
+ vect_t curr_to_prev;
+ vect_t plane_normal;
+ struct projection * current =
proj_array[array_pointer];
+ struct projection * previous =
proj_array[array_pointer - 1];
+ VSUB2(curr_to_prev,
current->original->point,
+ previous->original->point);
+ vect_t boundary;
+ point_t * intersectionp;
+ vect_t intersect_to_original;
+ fastf_t distance_to_originals;
+ VSCALE(curr_to_prev, curr_to_prev, 0.5);
+ puts("computing intersection");
+ intersectionp = intersection(inhit,
outhit, proj_array, current->original->point);
+ VSUB2(intersect_to_original,
*intersectionp, current->original->point);
+ /* Look for points closer to
intersection than the two original ones */
+ puts("arrived at inner for");
+ for (BU_LIST_FOR(dpoint_iter_int,
density_point, &(user_plist->l))) {
+ VSUB2(temp,
new_projection->point, dpoint_iter_int->point);
+ if (dpoint_iter_ext !=
dpoint_iter_int &&
+ MAGNITUDE(temp) <
distance_to_originals + MARGIN) {
+ puts("inner loop found
something");
+ /* FIXME: Deallocate
space we used for new_projection */
+ is_valid = 0;
+ bu_log("Skipped
projection. \n");
+ break;
+ }
+ }
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ }
+ }
+ /* We are at the end of one iteration on the outer
check loop.
+ If we detected a nonvalid point we should stop looking.
*/
+ if (!is_valid) break;
+ }
+ if (is_valid) {
+
+ 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);
+ }
+
+
+
+
+ }
- /* 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.
- */
+ /* 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);
+ }
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ qsort(proj_array, array_pointer, sizeof(struct projection *),
compare_dpoints);
+
+ /* 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 * curr_proj;
+ fastf_t acc_avg_density = 0.0;
+ struct projection * next_proj;
+
+ /* WARNING: From now on I refer to 'segment' as each piece of ray
inbetween two
+ points (projection or in/out */
+
+ /* First segment, we take everything as having density of first point */
+ curr_proj = proj_array[0];
+ acc_avg_density += curr_proj->original->density *
+ curr_proj->inhit_dist / total_segment_len;
+
+ /* DEBUG */
+ bu_log("First projection. Proportion %f, added %lf to accumulator.\n",
+ curr_proj->inhit_dist / total_segment_len,
curr_proj->original->density *
+ curr_proj->inhit_dist / total_segment_len);
+
+ /* Distance from current projection point to boundary */
+ fastf_t boundary_dist;
+ /* Distance between current projection point and next one */
+ fastf_t partial_segment_length;
+ vect_t partial_segment;
+ puts("loop start");
+ printf("Array pointer is %d", array_pointer);
+ next_proj = proj_array[0];
+ /* Loop through remaining segments. For each segment, compute the
boundary
+ and take density of each side accordingly, add it to accumulator */
+ for (int i = 0; i < array_pointer - 1; i++) {
+ puts("test zero");
+ curr_proj = proj_array[i];
+ next_proj = proj_array[i + 1];
+ puts("test one");
+ VSUB2(partial_segment, curr_proj->point, next_proj->point);
+ puts("test two");
+ partial_segment_length = MAGNITUDE(partial_segment);
+ boundary_dist = (next_proj->origin_dist *
next_proj->origin_dist +
+ partial_segment_length *
partial_segment_length -
+ curr_proj->origin_dist *
curr_proj->origin_dist) /
+ (2 * partial_segment_length);
+ puts("test three");
+ acc_avg_density += curr_proj->original->density *
+ boundary_dist / total_segment_len;
+ acc_avg_density += next_proj->original->density *
+ (partial_segment_length - boundary_dist) /
total_segment_len;
+ }
+ puts("end of for loop");
+ curr_proj = next_proj;
+ /* Last segment, from the last point to the outhit point has density
value
+ of this last point */
+ vect_t last_to_outhit;
+ VSUB2(last_to_outhit, outhit, curr_proj->point);
+ puts("did vsub");
+ acc_avg_density += curr_proj->original->density *
+ MAGNITUDE(last_to_outhit) / total_segment_len;
+ puts("did acc");
+ /* DEBUG */
+ bu_log("Last projection. Proportion %f, added %lf to accumulator.\n",
+ curr_proj->prev_dist / total_segment_len,
curr_proj->original->density *
+ curr_proj->prev_dist / total_segment_len);
+
+ /* Return final value */
+ bu_log("Avg density we saw through segment: %lf \n", acc_avg_density);
+ return acc_avg_density;
}
+int
+hit(struct application *ap, struct partition *PartHeadp, struct seg
*UNUSED(segs))
+{
+
+ struct partition *pp;
+ struct hit *hitp;
+ struct soltab *stp;
+ struct curvature cur = RT_CURVATURE_INIT_ZERO;
+ point_t pt;
+ vect_t inormal;
+ vect_t onormal;
+ /* Area that the ray covers, in mm^2 */
+ const int RAY_AREA = 4; //2mm height, 2mm width
+
+ for (pp = PartHeadp->pt_forw; pp != PartHeadp; pp = pp->pt_forw) {
+ bu_log("\n--- Hit region %s (in %s, out %s)\n",
+ pp->pt_regionp->reg_name,
+ pp->pt_inseg->seg_stp->st_name,
+ pp->pt_outseg->seg_stp->st_name);
+
+ hitp = pp->pt_inhit;
+
+ VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+
+ stp = pp->pt_inseg->seg_stp;
+ RT_HIT_NORMAL(inormal, hitp, stp, &(ap->a_ray), pp->pt_inflip);
+
+ rt_pr_hit(" In", hitp);
+ VPRINT(" Ipoint", pt);
+ VPRINT(" Inormal", inormal);
+
+ hitp = pp->pt_outhit;
+
+ VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+
+ stp = pp->pt_outseg->seg_stp;
+ RT_HIT_NORMAL(onormal, hitp, stp, &(ap->a_ray), pp->pt_outflip);
+
+ rt_pr_hit(" Out", hitp);
+ VPRINT(" Opoint", pt);
+ VPRINT(" Onormal", onormal);
+
+ /* Now we compute the total mass this element saw while
crossing the region */
+ fastf_t mass, volume, length, avg_density;
+
+ /*Here we query the avg density of a segment */
+ avg_density = segment_density(pp->pt_inhit->hit_point,
pp->pt_outhit->hit_point);
+
+ /* Get the length */
+ length = pp->pt_outhit->hit_dist - pp->pt_inhit->hit_dist;
+ length = length / 10.0; //mm to cm
+
+ /*We also need to know the volume */
+ volume = length * RAY_AREA;
+
+ mass = volume * avg_density;
+ bu_log(" Segment length: %lf, volume crossed: %lf \n", length,
volume);
+ bu_log(" Mass the ray saw through this region: %lf g \n",
mass);
+ }
+ return 1;
+}
+
+
/**
- * This is a callback routine that is invoked for every ray that
- * entirely misses hitting any geometry. This function is invoked by
- * rt_shootray() if the ray encounters nothing.
- */
+* This is a callback routine that is invoked for every ray that
+* entirely misses hitting any geometry. This function is invoked by
+* rt_shootray() if the ray encounters nothing.
+*/
int
miss(struct application *UNUSED(ap))
{
- bu_log("missed\n");
- return 0;
+ bu_log("missed\n");
+ return 0;
}
-
-/**
- * START HERE
- *
- * This is where it all begins.
- */
int
main(int argc, char **argv)
{
- /* Every application needs one of these. The "application"
- * structure carries information about how the ray-casting should
- * be performed. Defined in the raytrace.h header.
- */
- struct application ap;
+ /* Prepare user point list and let user fill it */
+ BU_GET(user_plist, struct density_point);
+ BU_LIST_INIT(&(user_plist->l));
+ struct transition * trans_list;
+ BU_GET(trans_list, struct transition);
+ BU_LIST_INIT(&(trans_list->l));
+ read_density_points(user_plist);
- /* 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
