Attached goes new code with the changes I mentioned. Only tested once on a
simple case.
Should give you an idea of what I meant with Saturdays mail.
Feedback is as always very welcome.
I will now test it more thoroughly and then move on to revamp the first
step too, which is severely flawed still.
I have yet to come up with a more elegant way of checking which 'tiles' our
segment will cross and which not.
Back with more tomorrow.
Mario.
On 23 September 2017 at 20:20, Mario Meissner <mr.rash....@gmail.com> wrote:
> Hello everyone.
>
> I wrote this on my log a few hours ago but I thought it's a good idea to
> go into more detail here since it's major change in direction.
>
> First, I will get rid of the projection structure, which is now just a
> contribution. Its projection is now just something secondary that can be
> used to predict if we need to skip the point or not.
>
> So, my new approach will be clearly separated into three steps:
> -Create a list of contributions. This stage is completely recycled from my
> current code, and for now the methodology to detect points that will
> contribute to our segment will remain the same. It will be the weakest part
> of the implementation once this major revamp is done, and may require some
> work.
> -Compute the boundaries between the contributions. This means computing
> when each contribution starts and ends within the segment. They are stored
> inside the contribution structure.
> -Loop through the contributions and take average density proportional to
> how big each intersection is.
>
> Notice how much simpler the average_density computation is now. Stage two
> is also simple because I have already implemented the intersection function.
>
> As this is a major change it will take me a few days to implement. I've
> spent most of my time today to refactor step 1 to use contributions instead
> of projections. I hope by tomorrow I can send an in-progress version with
> some cases working.
>
> I really think this new version cleans up the mess quite a lot, simplifies
> the concepts, makes it all more readable, and most importantly, will
> support more situations than before, removing several flaws.
>
> Mario.
>
> On 23 September 2017 at 18:01, Mario Meissner <mr.rash....@gmail.com>
> wrote:
>
>> 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.
>>> 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?
>>>
>>
>> I guess I should just do the same I did with intersection, since in the
>> end it's the same issue.
>>
>> 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.
>>
>
>
Index: src/rt/rtexample.c
===================================================================
--- src/rt/rtexample.c (リビジョン 70336)
+++ 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,549 @@
#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 density_vector {
+ struct bu_list l;
+ vect_t vector;
+ point_t origin;
+ fastf_t factor;
+};
+
+/* One possible way of defining vectors is doing it over already defined
points.
+Here we can specify what kind of transition happens between two points. */
+struct transition {
+ struct bu_list l;
+ struct density_point * origin;
+ struct density_point * termination;
+ fastf_t coefficient;
+};
+
+struct 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 intersection_t {
+ int intersect;
+ point_t intersection_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;
+}
+
+void compute_intersection(struct intersection_t * intersection, point_t
line_point0,
+ point_t line_point1, point_t plane_point, vect_t plane_normal) {
+ fastf_t epsilon = 1e-6;
+ fastf_t fac;
+ vect_t u, w;
+
+ VSUB2(u, line_point1, line_point0);
+ fastf_t dot = VDOT(plane_normal, u);
+ if (abs(dot) > epsilon) {
+ intersection->intersect = 1;
+ VSUB2(w, line_point0, plane_point);
+ fac = -VDOT(plane_normal, w) / dot;
+ VSCALE(u, u, fac);
+ VADD2(intersection->intersection_point, line_point0, u);
+ intersection->intersect = 1;
+ }
+ else {
+ intersection->intersect = 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
-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);
+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 (absolute(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;
+ }
+}
- /* print the entry hit point info */
- rt_pr_hit(" In", hitp);
- VPRINT( " Ipoint", pt);
- VPRINT( " Inormal", inormal);
- /* 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);
+/* Receives a segment and returns average density
+Uses memory pointer to user point list */
+fastf_t segment_density(point_t inhit, point_t outhit) {
- /* 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);
+ /* 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;
+ }
- /* exit point, so we type less */
- hitp = pp->pt_outhit;
+ /* DEBUG */
+ //VPRINT("inhit", inhit);
+ //VPRINT("outhit", outhit);
- /* construct the actual (exit) hit-point from the ray and the
- * distance to the intersection point (i.e., the 't' value).
- */
- VJOIN1(pt, ap->a_ray.r_pt, hitp->hit_dist, ap->a_ray.r_dir);
+ /* If 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;
+ }
- /* primitive we exited from */
- stp = pp->pt_outseg->seg_stp;
+ vect_t full_segment_vector;
+ VSUB2(full_segment_vector, outhit, inhit);
+ fastf_t full_segment_length = MAGNITUDE(full_segment_vector);
- /* 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);
+ /* Projections Loop prep */
+ struct contribution * new_contribution;
+ 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 dpoint_to_projection_vect;
+ vect_t temp;
+ vect_t paral_projection; //I need this to VPROJECT doesn't complain?
+ int i = 0;
+ int is_valid;
- /* print the exit hit point info */
- rt_pr_hit(" Out", hitp);
- VPRINT( " Opoint", pt);
- VPRINT( " Onormal", onormal);
- }
+ /* Bu_calloc has given me some problems, so temporary hardcoded max */
+ /* FIXME: I should be allocating space dynamically */
+ struct contribution * contrib_array[99];
- /* 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 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_contribution, struct contribution);
- /* Hit routine callbacks generally return 1 on hit or 0 on miss.
- * This value is returned by rt_shootray().
- */
- return 1;
+ /* TEST */
+ VSET(new_contribution->right_boundary, 0.0, 0.0, 0.0);
+
+ /* Create userpoint vector */
+ VSUB2(dpoint_vect, dpoint_iter_ext->point, inhit);
+ new_contribution->density_point = dpoint_iter_ext;
+ /* Project vector onto segment and compute absolute point of
the projection.
+ This will be necessary to evaluate if this contribution will
contribute or not. */
+ VPROJECT(dpoint_vect, full_segment_vector,
+ new_contribution->projection_vect, paral_projection);
+ VADD2(new_contribution->projection_point,
new_contribution->projection_vect, inhit);
+ new_contribution->inhit_dist =
MAGNITUDE(new_contribution->projection_vect);
+ //bu_log("inhit_dist: %lf \n", new_projection->inhit_dist);
+ /* Compute distance to the projection */
+ VSUB2(dpoint_to_projection_vect,
new_contribution->density_point->point,
+ new_contribution->projection_point);
+ new_contribution->projection_dist =
MAGNITUDE(dpoint_to_projection_vect);
+
+ /* 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_contribution->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_contribution->projection_dist) {
+
+ /* If we are here it means the 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 (i > 0) {
+ //puts("array pointer bigger than
zero");
+ vect_t prev_to_curr;
+ point_t plane_point;
+ struct contribution * previous =
contrib_array[i - 1];
+ //VPRINT("previous",
previous->original->point);
+ //VPRINT("current",
new_projection->original->point);
+ //puts("read current and previous, now
computing plane");
+ VSUB2(prev_to_curr,
new_contribution->density_point->point, previous->density_point->point);
+ struct intersection_t * intersection;
+ BU_GET(intersection, struct
intersection_t);
+ vect_t intersect_to_original;
+ fastf_t distance_to_originals;
+ VSCALE(prev_to_curr, prev_to_curr, 0.5);
+ //puts("computing intersection");
+ VADD2(plane_point,
previous->projection_point, prev_to_curr);
+ compute_intersection(intersection,
inhit, outhit, plane_point, prev_to_curr);
+ if (intersection->intersect == 0) {
+ bu_log("Error: No intersection
happened. \n");
+ return 0;
+ }
+ VSUB2(intersect_to_original,
intersection->intersection_point, new_contribution->density_point->point);
+ BU_FREE(intersection, struct
intersection_t);
+ distance_to_originals =
MAGNITUDE(intersect_to_original);
+ /* 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_contribution->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;
+ }
+ }
+
+ }
+ }
+ /* 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) {
+
+ /* FIXME: Maybe I don't need all this
+ if (i == 0) {
+ new_contribution->prev_dist =
new_contribution->inhit_dist;
+ }
+ else {
+ VSUB2(prev_to_current_vect,
new_contribution->projection_vect,
+ contrib_array[i - 1]->projection);
+ new_contribution->prev_dist =
MAGNITUDE(prev_to_current_vect);
+ }
+ */
+ contrib_array[i] = new_contribution;
+ i++;
+ /* DEBUG */
+ //VPRINT("projection vector", new_projection->vect);
+ }
+
+ }
+
+ /* Sort projection list based on distance to inhit */
+ /* DEBUG */
+ //puts("Unordered Projection distances:");
+ for (int i = 0; i < i; i++) {
+ //bu_log("Point at distance %lf \n",
contrib_array[i]->inhit_dist);
+ }
+
+ qsort(contrib_array, i, sizeof(struct contribution *), compare_dpoints);
+
+ /* DEBUG */
+ //puts("Ordered projection distances:");
+ for (int i = 0; i < i; i++) {
+ //bu_log("Point at distance %lf \n",
contrib_array[i]->inhit_dist);
+ }
+
+ int num_contributions = i;
+
+ /* STAGE TWO: COMPUTE BOUNDARIES */
+
+ /* First contribution starts at inhit */
+ VMOVE(contrib_array[0]->left_boundary, inhit);
+ /* Go through remaining points (second to last) and compute boundary to
it's left neigbor */
+ vect_t curr_to_prev;
+ point_t plane_point;
+ struct intersection_t * intersection;
+ for (i = 1; i < num_contributions; i++) {
+ /* This will be the boundary between previous contribution and
the current one */
+ BU_GET(intersection, struct intersection_t);
+ /* Compute the vector that connects the two points */
+ VSUB2(curr_to_prev, contrib_array[i-1]->density_point->point,
+ contrib_array[i]->density_point->point);
+ /* Divide it in half and add it to the current to get the
middle point between them */
+ VSCALE(curr_to_prev, curr_to_prev, 0.5);
+ VADD2(plane_point, contrib_array[i]->density_point->point,
curr_to_prev);
+ /* Get the intersection and check if it's correct */
+ compute_intersection(intersection, inhit, outhit, plane_point,
curr_to_prev);
+ if (intersection->intersect == 0) {
+ bu_log("No intersection happened while calculating
boundaries, aborting.\n ");
+ return 0;
+ }
+ /* Copy the data to out contributions structures */
+ VMOVE(contrib_array[i]->left_boundary,
intersection->intersection_point);
+ VMOVE(contrib_array[i - 1]->right_boundary,
intersection->intersection_point);
+ }
+ /* Add the outhit point as right boundary for the last point */
+ VMOVE(contrib_array[i-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<num_contributions; i++) {
+ VSUB2(contribution_segment, contrib_array[i]->left_boundary,
contrib_array[i]->right_boundary);
+ contribution_segment_len = MAGNITUDE(contribution_segment);
+ average_density += contrib_array[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));
+ 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, -2000, -2000, -2000);
+ VSET(ap.a_ray.r_dir, 1.0, 1.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;
- /* 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);
+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;
+}
- /* 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
