Hi people!
I proudly present my latest work, attached here. It features a clean
algorithm to compute the average density of a segment that goes through a
cloud of density regions arranged in a Voronoi Tesselation. For now
confirmed to work for one and two point cases, but I suppose only minor
bugs will hinder it from properly handling any arbitrary amount of points.
Note that it still requires some refactoring and cleaning since there are
some remains of the old code laying around.
The last step would be to integrate this into rtweight, which shouldn't be
too much of a problem.
As always, feedback welcome.
Thank you,
Mario.
On 6 October 2017 at 19:41, Mario Meissner <mr.rash....@gmail.com> wrote:
> Hi Daniel,
> thank you! In the end the problem was that I didn't run cmake after svn
> up.
>
> Attached goes today's work. I coded the function 'region_boundary()' which
> is an upgrade to the now unused intersection(). It does the same thing but
> now also checks if the intersection point lays inside or outside the
> segment, AND if the intersection is a valid region boundary or not (by
> checking distances to all other points, something I was doing outside of
> the intersection call until now). All in all a much cleaner function that
> will allow stages 1 and 2 to be much simpler and possibly run in one single
> loop.
>
> On Sunday I will probably finish implementing my new algorithm (which is
> now trivial since the hard work was this function).
>
> Note that current code is only meant to be a testing platform for the new
> function. On line 386 I included a snippet that calls my new function and
> prints its results, then I break out with a return. Feel free to try out
> some new intersections. Currently I tested with three points, first one
> ignored, and the second and third ones acting as 'previous' and 'current'
> respectively.
>
> Cheers!
>
> On 6 October 2017 at 17:41, Daniel Roßberg <danielmrossb...@gmail.com>
> wrote:
>
>> spectrum.c was recently removed. So, update your checked out sources and
>> check your librt CMakeaLists.txt that it matches with the one on the head
>> of the reposity (you've probably made changes to it).
>>
>> Regards,
>> Daniel
>>
>> Am 06.10.2017 5:27 nachm. schrieb "Mario Meissner" <mr.rash....@gmail.com
>> >:
>>
>>> I am getting compilation errors again, on 70355.
>>>
>>> fatal error C1083: Cannot open source file:
>>> 'C:\Users\MEISSNERBLANCOJOHANN\Desktop\brlcad\trunk2\src\librt\spectrum.c':
>>> No such file or directory
>>>
>>> And indeed no such file exists (I cannot find it anywhere).
>>>
>>> Hints? :/
>>>
>>> On 4 October 2017 at 23:14, Mario Meissner <mr.rash....@gmail.com>
>>> wrote:
>>>
>>>> So I ended up designing a new algorithm and would love to get some
>>>> feedback before I implement it.
>>>> Here it goes:
>>>>
>>>> * Sort the points from 'left to right' (distance of projection to inhit)
>>>> * Look for the closest point to inhit.
>>>> * Add this point to the list of contributions.
>>>> * Discard every point 'to the left' of this point.
>>>> * From the next point onwards, do, for each point:
>>>> * Take the next point and the last one we added to contributions
>>>> * Check if the intersection of their middle-plain and the segment
>>>> lays inside the segment and
>>>>
>>>> has no closer point than the two we are using..
>>>> * If so, this point's region is crossed by the segment, so we add it
>>>> to the contribution list.
>>>> * If not, then this point's region is not crossed.
>>>> * Look for the closest point to the outhit.
>>>> * If this point is the same as the last point we added, finish.
>>>> * If not, add it and finish.
>>>>
>>>> Attached go three pics that are a sample of the example situations I
>>>> used to develop the algorithm.
>>>> Ideally, I would need someone to find a situation where this algorithm
>>>> fails.
>>>>
>>>> Thank you in advance,
>>>> Mario.
>>>>
>>>> On 3 October 2017 at 04:14, Christopher Sean Morrison <brl...@mac.com>
>>>> wrote:
>>>>
>>>>>
>>>>> Hi Mario,
>>>>>
>>>>> As always, thank you for all the updates. Here are some responses to
>>>>> questions you posted last week.
>>>>>
>>>>> I'm passing two arguments to my segment_density, which are inhit and
>>>>>> outhit points.
>>>>>> I now realize that these are pointers and could change unexpectedly
>>>>>> within the call if we run multiple threads.
>>>>>>
>>>>>
>>>>> The in/out hit point pointers shouldn’t be changing within the call.
>>>>> I suspect something else is going on...
>>>>>
>>>>> How can I make this safe? Should I BU_GET safe heap memory to store
>>>>>> the points for the call, or pass the coordinates one by one as actual
>>>>>> numbers?
>>>>>>
>>>>>
>>>>> In general, the calling code should pass structures to functions
>>>>> (which can be on the heap or on the stack) and those functions should just
>>>>> fill them in or read from them (but not both, separate functions for
>>>>> reading and writing).
>>>>>
>>>>> For debugging purposes I would like to set only one thread for
>>>>> rtweight, how can I do that? The fact that I get 4 or more threads running
>>>>> my segment_density at the same time makes it difficult to track down
>>>>> issues.
>>>>>
>>>>>
>>>>> Check out rtweight’s manual page (run "brlman rtweight” outside mged),
>>>>> it’s the -P# option.
>>>>>
>>>>> 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
>>>>>
>>>>>
>>>>
>>>
>>> ------------------------------------------------------------
>>> ------------------
>>> 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
>>>
>>>
>> ------------------------------------------------------------
>> ------------------
>> 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 (リビジョン 70355)
+++ src/rt/rtexample.c (作業コピー)
@@ -1,65 +1,65 @@
/* R T E X A M P L E . C
- * BRL-CAD
- *
- * Copyright (c) 2004-2016 United States Government as represented by
- * the U.S. Army Research Laboratory.
- *
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * version 2.1 as published by the Free Software Foundation.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with this file; see the file named COPYING for more
- * information.
- */
+* BRL-CAD
+*
+* Copyright (c) 2004-2016 United States Government as represented by
+* the U.S. Army Research Laboratory.
+*
+* This program is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License
+* version 2.1 as published by the Free Software Foundation.
+*
+* This program is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+* Lesser General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public
+* License along with this file; see the file named COPYING for more
+* information.
+*/
/** @file rt/rtexample.c
- *
- * This is a heavily commented example of a program that uses librt to
- * shoot a single ray at some geometry in a .g database.
- *
- * The primary BRL-CAD ray-tracing API consists of calls to the
- * function rt_shootray(). This function takes a single argument, the
- * address of a data structure called the "application" structure.
- * This data structure contains crucial information, such as the
- * origin and direction of the ray to be traced, what to do if the ray
- * hits geometry, or what to do if it misses everything. While the
- * application struct is a large and somewhat complex looking, (it is
- * defined in the raytrace.h header) there are really very few items
- * in it which the application programmer must know about. These are:
- *
- * a_rt_i The "raytrace instance" obtained via rt_dirbuild()
- * a_ray The ray origin and direction to be shot
- * a_hit A callback function for when the ray encounters geometry
- * a_miss A callback function for when the ray misses everything
- *
- * Most of the work an application performs will be done in the "hit"
- * routine. This user-supplied routine gets called deep inside the
- * raytracing library via the rt_shootray() function. It is provided
- * with 3 parameters:
- *
- * ap Pointer to the application structure passed to
rt_shootray()
- * PartHeadp List of ray "partitions" which represent geometry hit
- * segp List of ray "segments" that comprised partitions
- *
- * Most applications can ignore the last parameter. The PartHeadp
- * parameter is a linked-list of "partition" structures (defined in
- * the raytrace.h header). It is the job of the "hit" routine to
- * process these ray/object intersections to do the work of the
- * application.
- *
- * This file is part of the default compile in source distributions of
- * BRL-CAD and is usually installed or provided via binary and source
- * distributions. To compile this example from a binary install:
- *
- * cc -I/usr/brlcad/include/brlcad -L/usr/brlcad/lib -o rtexample rtexample.c
-lbu -lrt -lm
- *
- * Jump to the START HERE section below for main().
- */
+*
+* This is a heavily commented example of a program that uses librt to
+* shoot a single ray at some geometry in a .g database.
+*
+* The primary BRL-CAD ray-tracing API consists of calls to the
+* function rt_shootray(). This function takes a single argument, the
+* address of a data structure called the "application" structure.
+* This data structure contains crucial information, such as the
+* origin and direction of the ray to be traced, what to do if the ray
+* hits geometry, or what to do if it misses everything. While the
+* application struct is a large and somewhat complex looking, (it is
+* defined in the raytrace.h header) there are really very few items
+* in it which the application programmer must know about. These are:
+*
+* a_rt_i The "raytrace instance" obtained via rt_dirbuild()
+* a_ray The ray origin and direction to be shot
+* a_hit A callback function for when the ray encounters geometry
+* a_miss A callback function for when the ray misses everything
+*
+* Most of the work an application performs will be done in the "hit"
+* routine. This user-supplied routine gets called deep inside the
+* raytracing library via the rt_shootray() function. It is provided
+* with 3 parameters:
+*
+* ap Pointer to the application structure passed to rt_shootray()
+* PartHeadp List of ray "partitions" which represent geometry hit
+* segp List of ray "segments" that comprised partitions
+*
+* Most applications can ignore the last parameter. The PartHeadp
+* parameter is a linked-list of "partition" structures (defined in
+* the raytrace.h header). It is the job of the "hit" routine to
+* process these ray/object intersections to do the work of the
+* application.
+*
+* This file is part of the default compile in source distributions of
+* BRL-CAD and is usually installed or provided via binary and source
+* distributions. To compile this example from a binary install:
+*
+* cc -I/usr/brlcad/include/brlcad -L/usr/brlcad/lib -o rtexample rtexample.c
-lbu -lrt -lm
+*
+* Jump to the START HERE section below for main().
+*/
#include "common.h"
@@ -67,264 +67,479 @@
#include <math.h>
#include <string.h>
#include <stdio.h>
+#include <string.h>
#include "vmath.h" /* vector math macros */
#include "raytrace.h" /* librt interface definitions */
-/**
- * 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.
- */
+fastf_t MARGIN = 0.0001;
+struct density_point {
+ struct bu_list l;
+ point_t point;
+ fastf_t density;
+};
+
+struct contribution {
+ struct density_point * density_point;
+ //point_t projection_point;
+ //vect_t projection_vect;
+ fastf_t inhit_dist;
+ //fastf_t prev_dist;
+ //fastf_t projection_dist;
+ point_t left_boundary;
+ point_t right_boundary;
+};
+
+struct boundary_t {
+ int is_valid;
+ point_t boundary_point;
+};
+
+struct density_point * user_plist;
+const int MAX_POINTS = 99;
+
+int compare_dpoints(const void *, const void *);
+fastf_t absolute(fastf_t value) {
+ if (value < 0) {
+ return -value;
+ }
+ else {
+ return value;
+ }
+}
+
+/* Check if c is between a and b */
+int isBetween(point_t c, point_t a, point_t b) {
+ vect_t c_minus_a, b_minus_a;
+ VSUB2(c_minus_a, c, a);
+ VSUB2(b_minus_a, b, a);
+ fastf_t dot_product = VDOT(c_minus_a, b_minus_a);
+ fastf_t squared_length_ba = MAGNITUDE(b_minus_a)*MAGNITUDE(b_minus_a);
+ return dot_product > 0 && dot_product < squared_length_ba;
+}
+
+/*
+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 head 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;
+ 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);
- /* primitive we encountered on entry */
- stp = pp->pt_inseg->seg_stp;
+ return count;
+}
- /* 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);
+/* Receives two points: current and previous, and a segment, and finds out
+if the segment will cross the 'current point's region or not.
+Uses global userpoint list variable */
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
+void region_boundary(struct boundary_t * boundary, int current_index, int
previous_index,
+ struct contribution ** points_array, int array_size, point_t
segment_start, point_t segment_end) {
- /* 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);
+ point_t current, previous;
+ VMOVE(current, points_array[current_index]->density_point->point);
+ VMOVE(previous, points_array[previous_index]->density_point->point);
+ int valid_intersection, valid_boundary;
+ point_t plane_point, boundary_point;
+ vect_t plane_normal;
+ fastf_t epsilon = 1e-6;
+ fastf_t fac;
+ vect_t u, w;
+ VSUB2(plane_normal, current, previous);
+ VSCALE(plane_normal, plane_normal, 0.5);
+ VADD2(plane_point, previous, plane_normal);
+ VSUB2(u, segment_end, segment_start);
+ fastf_t dot = VDOT(plane_normal, u);
+ VSUB2(w, segment_start, plane_point);
+ fac = -VDOT(plane_normal, w) / dot;
+ VSCALE(u, u, fac);
+ VADD2(boundary_point, segment_start, u);
+ /* Check if intersection is valid */
+ valid_intersection = abs(dot) > epsilon &&
+ isBetween(boundary_point, segment_start, segment_end);
- /* 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);
+ /* Break out if we don't intersect */
+ if (!valid_intersection) {
+ boundary->is_valid = 0;
+ return;
+ }
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ valid_boundary = 1;
+ int index;
+ vect_t vector;
+ point_t point;
+ VSUB2(vector, boundary_point, current);
+ /* This is the distance between the intersection and our points */
+ fastf_t boundary_distance = MAGNITUDE(vector);
+ fastf_t point_distance;
+ /* Now check distances to see if boundary is valid */
+ for (index = 0; index < array_size; index ++) {
+ if (index == current_index || index == previous_index) continue;
+ VMOVE(point, points_array[index]->density_point->point);
+ VSUB2(vector, point, boundary_point);
+ point_distance = MAGNITUDE(vector);
+ if (point_distance < boundary_distance) {
+ valid_boundary = 0;
+ break;
+ }
+ }
- /* 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 (valid_boundary) {
+ boundary->is_valid = 1;
+ VMOVE(boundary->boundary_point, boundary_point);
+ }
+ else {
+ boundary->is_valid = 0;
+ }
+}
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+int search_closest_point(point_t point, struct contribution ** point_list, int
array_size) {
+ int pos_min;
+ fastf_t dist;
+ fastf_t min_dist = FLT_MAX;
+ vect_t vect;
+ point_t candidate;
+ for (int i = 0; i < array_size; i++) {
+ VMOVE(candidate, point_list[i]->density_point->point);
+ VSUB2(vect, candidate, point);
+ dist = MAGNITUDE(vect);
+ if (dist < min_dist) {
+ pos_min = i;
+ min_dist = dist;
+ }
+ }
+ return pos_min;
+}
- /* 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);
- }
+/* Receives a segment and returns average density
+Uses memory pointer to user point list */
+fastf_t segment_density(point_t inhit, point_t outhit) {
- /* 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.
- */
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ /* If no points, mass is 0 */
+ if (BU_LIST_IS_EMPTY(&(user_plist->l))) {
+ bu_log("No points were provided, thus object has no mass.\n");
+ return 0;
+ }
+
+ /* If there is only one point we skip the paraphernalia */
+ if (bu_list_len(&user_plist->l) == 1) {
+ struct density_point * dp;
+ dp = (struct density_point *)user_plist->l.forw;
+ return dp->density;
+ }
+
+ vect_t full_segment_vector;
+ VSUB2(full_segment_vector, outhit, inhit);
+ fastf_t full_segment_length = MAGNITUDE(full_segment_vector);
+
+ /* Projections Loop prep */
+ struct contribution * new_contribution;
+ struct density_point * dpoint_iter;
+ vect_t dpoint_vect;
+ vect_t prev_to_current_vect;
+ vect_t dpoint_to_projection_vect;
+ vect_t temp;
+ vect_t projection_vector;
+ vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+ /* Bu_calloc has given me some problems, so temporary hardcoded max */
+ /* FIXME: I should be allocating space dynamically */
+
+
+ /* Array variables */
+ /* Array containing the points that will contribute (i.e. we will cross
the region) */
+ struct contribution * contributions[99];
+ /* Array containing all points */
+ struct contribution * userpoints[99];
+ int userpoints_size = 0;
+ int contributions_size = 0;
+ /* Index of the userpoint we are currently working with */
+ int userpoints_index = 0;
+ /* Index of the next position where we will store a contribution */
+ int contributions_index = 0;
+
+ /* Copy all points to the array (and convert them to contribution type)
*/
+ for (BU_LIST_FOR(dpoint_iter, density_point, &(user_plist->l))) {
+ BU_GET(new_contribution, struct contribution);
+ /* Create userpoint vector */
+ VSUB2(dpoint_vect, dpoint_iter->point, inhit);
+ /* Compute projected distance to inhit */
+ VPROJECT(dpoint_vect, full_segment_vector,
+ projection_vector, paral_projection);
+ new_contribution->inhit_dist = MAGNITUDE(projection_vector);
+ new_contribution->density_point = dpoint_iter;
+ userpoints[userpoints_index] = new_contribution;
+ userpoints_size++;
+ userpoints_index++;
+ }
+
+ /* Sort all points from 'left to right', i.e. distance to inhit */
+ qsort(userpoints, userpoints_size, sizeof(struct contribution *),
compare_dpoints);
+
+ /* Reset the indeces */
+ userpoints_index = 0;
+ contributions_index = 0;
+
+ /* Search the closest point to inhit since it's gonna be the first
region we cross */
+ userpoints_index = search_closest_point(inhit, userpoints,
userpoints_size);
+ contributions[contributions_index] = userpoints[userpoints_index];
+ VMOVE(contributions[contributions_index]->left_boundary, inhit);
+ contributions_size++;
+ contributions_index++;
+
+ /* Position inside userpoints where the last point we added to
contributions is */
+ int last_added_userpoint = userpoints_index;
+ userpoints_index++;
+ /* Now loop through the rest in order, checking when we cross the next
region */
+ struct boundary_t * boundary;
+ BU_GET(boundary, struct boundary_t);
+ for (userpoints_index; userpoints_index < userpoints_size;
userpoints_index++) {
+ region_boundary(boundary, userpoints_index,
last_added_userpoint,
+ userpoints, userpoints_size, inhit, outhit);
+ if (boundary->is_valid) {
+ /* If valid, add current userpoint to the contributions
list */
+ contributions[contributions_index] =
userpoints[userpoints_index];
+
VMOVE(contributions[contributions_index]->left_boundary,
boundary->boundary_point);
+ VMOVE(contributions[contributions_index -
1]->right_boundary, boundary->boundary_point);
+ contributions_index++;
+ contributions_size++;
+ last_added_userpoint = userpoints_index;
+ }
+ }
+ BU_FREE(boundary, struct boundary_t);
+
+ /* Finally the last right_boundary is outhit */
+ VMOVE(contributions[contributions_size - 1]->right_boundary, outhit);
+
+ /* STAGE THREE: COMPUTE MASS ACCORDING TO CONTRIBUTION DISTRIBUTION */
+ vect_t contribution_segment;
+ fastf_t average_density, contribution_segment_len;
+ average_density = 0.0;
+
+
+ /* Loop through contributions and compute average density */
+ for (int i = 0; i<contributions_size; i++) {
+ VSUB2(contribution_segment, contributions[i]->left_boundary,
contributions[i]->right_boundary);
+ contribution_segment_len = MAGNITUDE(contribution_segment);
+ average_density += contributions[i]->density_point->density *
contribution_segment_len / full_segment_length;
+ }
+ return average_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 fastf_t RAY_AREA = 0.2 * 0.2; //0.2cm height, 0.2cm 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 cm, volume crossed: %lf \n",
length, volume);
+ bu_log(" Average density through the segment: %lf \n",
avg_density);
+ 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));
+ 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);
+ }
+ //END DEBUG
- /* optional parameter to rt_dirbuild() that can be used to capture
- * a title if the geometry database has one set.
- */
- char title[1024] = {0};
- /* 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]);
- }
+ /* 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]);
+ }
- /* 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]);
- }
+ 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]);
+ }
- /* Display the geometry database title obtained during
- * rt_dirbuild if a title is set.
- */
- if (title[0]) {
- bu_log("Title:\n%s\n", title);
- }
+ if (title[0]) {
+ bu_log("Title:\n%s\n", title);
+ }
- /* 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++;
- }
+ while (argc > 2) {
+ if (rt_gettree(rtip, argv[2]) < 0)
+ bu_log("Loading the geometry for [%s] FAILED\n",
argv[2]);
+ argc--;
+ argv++;
+ }
- /* 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);
+ rt_prep_parallel(rtip, 1);
+ RT_APPLICATION_INIT(&ap);
+ ap.a_rt_i = rtip;
+ ap.a_onehit = 0;
+ VSET(ap.a_ray.r_pt, -10000, 500, 500);
+ VSET(ap.a_ray.r_dir, 1.0, 0.0, 0.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);
- /* initialize all values in application structure to zero */
- RT_APPLICATION_INIT(&ap);
+ return 0;
+}
- /* your application uses the raytrace instance containing the
- * geometry we loaded. this describes what we're shooting at.
- */
- ap.a_rt_i = rtip;
- /* 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;
+int compare_dpoints(const void *a, const void *b) {
+ struct contribution *x = *(struct contribution **)a;
+ struct contribution *y = *(struct contribution **)b;
+ return x->inhit_dist - y->inhit_dist;
+}
- /* 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);
- /* Simple debug printing */
- VPRINT("Pnt", ap.a_ray.r_pt);
- VPRINT("Dir", ap.a_ray.r_dir);
- /* This is what callback to perform on a hit. */
- ap.a_hit = hit;
-
- /* 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.
- */
-
- return 0;
-}
-
/*
- * Local Variables:
- * mode: C
- * tab-width: 8
- * indent-tabs-mode: t
- * c-file-style: "stroustrup"
- * End:
- * ex: shiftwidth=4 tabstop=8
- */
+* Local Variables:
+* mode: C
+* tab-width: 8
+* indent-tabs-mode: t
+* c-file-style: "stroustrup"
+* End:
+* ex: shiftwidth=4 tabstop=8
+*/
------------------------------------------------------------------------------
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
