BRL-CAD
tie.c
Go to the documentation of this file.
1 /* T I E . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2008-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file primitives/bot/tie.c
21  * tie.c
22  *
23  */
24 
25 #include "tie.h"
26 
27 #include <math.h>
28 #include <stdlib.h>
29 #include <string.h>
30 
31 #ifdef TIE_SSE
32 # include <xmmintrin.h>
33 #endif
34 
35 #ifdef HAVE_STDINT_H
36 # include <stdint.h>
37 #endif
38 
39 #include "bio.h"
40 
41 
42 #include "tieprivate.h"
43 
44 #define TIE_TAB1 "\1\0\0\2\2\1" /* Triangle Index Table */
45 
46 #define TIE_DEGENERATE_THRESHOLD 0.0001
47 
48 /*************************************************************
49  **************** PRIVATE FUNCTIONS **************************
50  *************************************************************/
51 
52 static void
53 TIE_VAL(tie_tri_prep)(struct tie_s *tie)
54 {
55  TIE_3 v1, v2, u, v;
56  unsigned int i, i1, i2;
57  struct tie_tri_s *tri;
58  fastf_t mag_sq;
59 
60  for (i = 0; i < tie->tri_num; i++) {
61  tri = &tie->tri_list[i];
62 
63  v1 = tri->data[1];
64  v2 = tri->data[2];
65 
66  /* Compute Normal */
67  VSUB2(u.v, tri->data[1].v, tri->data[0].v);
68  VSUB2(v.v, tri->data[2].v, tri->data[0].v);
69  VCROSS(tri->data[1].v, u.v, v.v);
70 
71  /* Unitize Normal */
72  mag_sq = MAGSQ(tri->data[1].v);
73  if (UNLIKELY(mag_sq < SMALL_FASTF)) {
74  /* Can not unitize normal, most likely we have a zero area
75  * triangle (i.e. degenerate) so skip it.
76  */
77  continue;
78  }
79  VSCALE(tri->data[1].v, tri->data[1].v, 1.0/sqrt(mag_sq));
80 
81  /* Compute i1 and i2 */
82  u.v[0] = fabs(tri->data[1].v[0]);
83  u.v[1] = fabs(tri->data[1].v[1]);
84  u.v[2] = fabs(tri->data[1].v[2]);
85 
86  if (u.v[2] > u.v[1] && u.v[2] > u.v[0]) {
87  i1 = 0;
88  i2 = 1;
89  } else if (u.v[1] > u.v[2] && u.v[1] > u.v[0]) {
90  i1 = 0;
91  i2 = 2;
92  } else {
93  i1 = 1;
94  i2 = 2;
95  }
96 
97  /* compute u1, v2, u2, v2 */
98  tri->data[2].v[1] = v1.v[i1] - tri->data[0].v[i1];
99  tri->data[2].v[2] = v2.v[i1] - tri->data[0].v[i1];
100  tri->v[0] = v1.v[i2] - tri->data[0].v[i2];
101  tri->v[1] = v2.v[i2] - tri->data[0].v[i2];
102 
103  /* looks like this is intentionally injected to give the
104  * compiler something to do if/while v1 and tri->data are
105  * preloaded, instead of just getting set in the preceding
106  * conditional above.
107  */
108  if (i1 == 0 && i2 == 1)
109  tri->b = tri->b + 2;
110  else if (i1 == 0)
111  tri->b = tri->b + 1;
112 
113  /* Compute DotVN */
114  VSCALE(v1.v, tri->data[0].v, -1.0);
115  tri->data[2].v[0] = VDOT(v1.v, tri->data[1].v);
116  }
117 }
118 
119 
120 /*************************************************************
121  **************** EXPORTED FUNCTIONS *************************
122  *************************************************************/
123 
124 /**
125  * Initialize a struct tie_s data structure
126  *
127  * This needs to be called before any other libtie data structures are
128  * called.
129  *
130  * @param tie pointer to a struct tie_t
131  * @param tri_num initial number of triangles to allocate for.
132  * tie_push may expand the buffer, if needed.
133  * @param kdmethod Either TIE_KDTREE_FAST or TIE_KDTREE_OPTIMAL
134  * @return void
135  */
136 void
137 TIE_VAL(tie_init)(struct tie_s *tie, unsigned int tri_num, unsigned int kdmethod)
138 {
139  tie->kdtree = NULL;
140  tie->kdmethod = kdmethod;
141  tie->tri_num = 0;
142  tie->tri_num_alloc = tri_num;
143  /* set to NULL if tri_num == 0. bu_malloc(0) causes issues. */
144  tie->tri_list = tri_num?(struct tie_tri_s *)bu_calloc(tri_num, sizeof(struct tie_tri_s), "tie_init"):NULL;
145  tie->stat = 0;
146  tie->rays_fired = 0;
147 }
148 
149 
150 /**
151  * Free up all the stuff associated with libtie
152  *
153  * All of the KDTREE nodes and triangles that we have allocated need
154  * to be freed in a controlled manner. This routine does that.
155  *
156  * @param tie pointer to a struct tie_t
157  * @return void
158  */
159 void
160 TIE_VAL(tie_free)(struct tie_s *tie)
161 {
162  bu_free(tie->tri_list, "tie_free");
163 
164  /* Free KDTREE Nodes */
165  tie_kdtree_free(tie);
166 }
167 
168 
169 /**
170  * Get ready to shoot rays at triangles
171  *
172  * Build the KDTREE tree for the triangles we have
173  *
174  * @param tie pointer to a struct struct tie_s which now has all the
175  * triangles in it.
176  *
177  * @return void
178  */
179 void
180 TIE_VAL(tie_prep)(struct tie_s *tie)
181 {
182  /* Build the kd-tree */
183  tie_kdtree_prep (tie);
184 
185  /* Prep all the triangles */
186  TIE_VAL(tie_tri_prep) (tie);
187 }
188 
189 
190 /**
191  * Shoot a ray at some triangles
192  *
193  * The user-provided hitfunc is called at each ray/triangle
194  * intersection. Calls are guaranteed to be made in the
195  * ray-intersection order. The last argument (void *ptr) is passed to
196  * the hitfunc as-is, to allow application specific data to be passed
197  * to the hitfunc.
198  *
199  * @param tie a struct struct tie_s universe
200  * @param ray the ray to be intersected with the geometry
201  * @param id the intersection data for each intersection
202  * @param hitfunc the application routine to be called upon ray/triangle intersection.
203  *
204  * This function should return 0 if the ray is to continue propagating
205  * through the geometry, or non-zero if ray intersection should cease.
206  * @param ptr a pointer to be passed to the hitfunc when it is called.
207  *
208  * @return the return value from the user hitfunc() is used.
209  *
210  * In the event that the ray does not hit anything, or the ray exits
211  * the geometry space, a null value will be returned.
212  *
213  * @retval 0 ray did not hit anything, or ray was propagated through the geometry completely.
214  * @retval !0 the value returned from the last invocation of hitfunc()
215  */
216 void *
217 TIE_VAL(tie_work)(struct tie_s *tie, struct tie_ray_s *ray, struct tie_id_s *id, void *(*hitfunc)(struct tie_ray_s*, struct tie_id_s*, struct tie_tri_s*, void *ptr), void *ptr)
218 {
219  struct tie_stack_s stack[40];
220  struct tie_id_s t = {{0, 0, 0}, {0, 0, 0}, 0, 0, 0}, id_list[256];
221  struct tie_tri_s *hit_list[256], *tri;
222  struct tie_geom_s *data;
223  struct tie_kdtree_s *node, *temp[2];
224  tfloat near, far, dirinv[3], dist;
225  unsigned int i, n;
226  uint8_t hit_count;
227  int ab[3], split, stack_ind;
228  void *result;
229 
230  if (!tie->kdtree)
231  return NULL;
232 
233  ray->kdtree_depth = 0;
234 
235  /*
236  * Precompute direction inverse since it's used in a bunch of
237  * divides, this allows those divides to become fast multiplies.
238  */
239  for (i = 0; i < 3; i++) {
240  if (ZERO(ray->dir[i]))
241  ray->dir[i] = TIE_PREC;
242  dirinv[i] = 1.0 / ray->dir[i];
243  ab[i] = dirinv[i] < 0.0 ? 1.0 : 0.0;
244  }
245 
246  /* Extracting value of splitting plane from the kdtree */
247  split = tie->kdtree->b & (uint32_t)0x3L;
248 
249  /* Initialize ray segment */
250  if (ray->dir[split] < 0.0)
251  far = (tie->min[split] - ray->pos[split]) * dirinv[split];
252  else
253  far = (tie->max[split] - ray->pos[split]) * dirinv[split];
254 
255  stack_ind = 0;
256  stack[0].node = tie->kdtree;
257  stack[0].near = 0;
258  stack[0].far = far;
259 
260  /* Process items on the stack */
261  do {
262  near = stack[stack_ind].near;
263  far = stack[stack_ind].far;
264  node = stack[stack_ind].node;
265  stack_ind--;
266 
267  /*
268  * KDTREE TRAVERSAL
269  *
270  * 3 conditions can happen here:
271  * - Ray only intersects the nearest node
272  * - Ray only intersects the furthest node
273  * - Ray intersects both nodes, pushing the furthest onto the stack
274  *
275  * Gordon Stoll's Mantra - Rays are Measured in Millions :-)
276  */
277  while (node && node->data && TIE_HAS_CHILDREN(node->b)) {
278  ray->kdtree_depth++;
279 
280  /* Retrieve the splitting plane */
281  split = node->b & (uint32_t)0x3L;
282 
283  /* Calculate the projected 1d distance to splitting axis */
284  dist = (node->axis - ray->pos[split]) * dirinv[split];
285 
286  temp[0] = &((struct tie_kdtree_s *)(node->data))[ab[split]];
287  temp[1] = &((struct tie_kdtree_s *)(node->data))[1-ab[split]];
288 
289  i = near >= dist; /* Node B Only? */
290  node = temp[i];
291 
292  if (far < dist || i)
293  continue;
294 
295  /* Nearest Node and Push Furthest */
296  stack_ind++;
297  stack[stack_ind].node = temp[1];
298  stack[stack_ind].near = dist;
299  stack[stack_ind].far = far;
300  far = dist;
301  }
302 
303  if (!node || !node->data) {
304  /* bu_log("INTERNAL ERROR: Unexpected TIE tree traversal encountered\n"); */
305  continue;
306  }
307 
308  /*
309  * RAY/TRIANGLE INTERSECTION - Only gets executed on geometry
310  * nodes. This part of the function is being executed because
311  * the KDTREE Traversal is Complete.
312  */
313 
314  hit_count = 0;
315  data = (struct tie_geom_s *)(node->data);
316 
317  for (i = 0; i < data->tri_num; i++) {
318  /*
319  * Triangle Intersection Code
320  */
321  tfloat u0, v0, *v;
322  int i1, i2;
323 
324  tri = data->tri_list[i];
325  u0 = VDOT(tri->data[1].v, ray->pos);
326  v0 = VDOT(tri->data[1].v, ray->dir);
327 
328  /* skip rays that are practically perpendicular so we
329  * don't try to divide by zero and propagate NaN (or
330  * crash).
331  */
332  if (ZERO(v0))
333  continue;
334 
335  t.dist = -(tri->data[2].v[0] + u0) / v0;
336 
337  if (isnan(t.dist))
338  continue;
339 
340  /*
341  * Intersection point on triangle must lie within the
342  * kdtree node or it is rejected Apply TIE_PREC to near
343  * and far such that triangles that lie on orthogonal
344  * planes aren't in a precision fuzz boundary, thus
345  * missing something they should actually have hit.
346  */
347  if (t.dist < near-TIE_PREC || t.dist > far+TIE_PREC)
348  continue;
349 
350  /* Compute Intersection Point (P = O + Dt) */
351  VSCALE(t.pos, ray->dir, t.dist);
352  VADD2(t.pos, ray->pos, t.pos);
353 
354  /* Extract i1 and i2 indices from the 'b' field */
355  v = tri->v;
356 
357  i1 = TIE_TAB1[tri->b & (uint32_t)0x7L];
358  i2 = TIE_TAB1[3 + (tri->b & (uint32_t)0x7L)];
359 
360  /* Compute U and V */
361  u0 = t.pos[i1] - tri->data[0].v[i1];
362  v0 = t.pos[i2] - tri->data[0].v[i2];
363 
364  /*
365  * Compute the barycentric coordinates, and make sure the
366  * coordinates fall within the boundaries of the triangle
367  * plane.
368  */
369  if (fabs(tri->data[2].v[1]) <= TIE_PREC) {
370  t.beta = u0 / tri->data[2].v[2];
371  if (t.beta < 0 || t.beta > 1)
372  continue;
373  t.alpha = (v0 - t.beta*v[1]) / v[0];
374  } else {
375  t.beta = (v0*tri->data[2].v[1] - u0*v[0]) / (v[1]*tri->data[2].v[1] - tri->data[2].v[2]*v[0]);
376  if (t.beta < 0 || t.beta > 1)
377  continue;
378  t.alpha = (u0 - t.beta*tri->data[2].v[2]) / tri->data[2].v[1];
379  }
380 
381  if (t.alpha < 0 || (t.alpha + t.beta) > 1.0)
382  continue;
383 
384  /* Triangle Intersected, append it in the list */
385  if (hit_count < 0xFF) {
386  hit_list[hit_count] = tri;
387  id_list[hit_count] = t;
388  hit_count++;
389  }
390  }
391 
392 
393  /* If we hit something, then sort the hit triangles on demand */
394  for (i = 0; i < hit_count; i++) {
395  /* Sort the list so that HitList and IDList [n] is in order wrt [i] */
396  for (n = i; n < hit_count; n++) {
397  if (id_list[n].dist < id_list[i].dist) {
398  /* Swap */
399  tri = hit_list[i];
400  t = id_list[i];
401  hit_list[i] = hit_list[n];
402  id_list[i] = id_list[n];
403  hit_list[n] = tri;
404  id_list[n] = t;
405  }
406  }
407 
408  VMOVE(id_list[i].norm, hit_list[i]->data[1].v);
409  result = hitfunc(ray, &id_list[i], hit_list[i], ptr);
410 
411  if (result) {
412  *id = id_list[i];
413  return result;
414  }
415  }
416  } while (stack_ind >= 0);
417 
418  return NULL;
419 }
420 
421 
422 /**
423  * Add a triangle
424  *
425  * Add a new triangle to the universe to be raytraced.
426  *
427  * @param tie the universe
428  * @param tlist is an array of TIE_3 vertex triplets (v0, v1, v2) that
429  * form each triangle.
430  * @param tnum is the number of triangles (tlist = 3 * tnum of TIE_3's).
431  * @param plist is a list of pointer data that gets assigned to the
432  * ptr of each triangle.
433  *
434  * This will typically be 4-byte (32-bit architecture) spaced array of
435  * pointers that associate the triangle pointer with your arbitrary
436  * structure, i.e. a mesh.
437  *
438  * @param pstride is the number of bytes to increment the pointer list
439  * as it assigns a pointer to each mesh, typically a value of 4 (for
440  * 32-bit machines). If you have a single pointer that groups all
441  * triangles to a common structure then you can use a value of 0 for
442  * pstride. This will give the pointer of all triangles the pointer
443  * address of plist. @return void
444  */
445 void
446 TIE_VAL(tie_push)(struct tie_s *tie, TIE_3 **tlist, unsigned int tnum, void *plist, unsigned int pstride)
447 {
448  unsigned int i;
449 
450  /* expand the tri buffer if needed */
451  if (tnum + tie->tri_num > tie->tri_num_alloc) {
452  tie->tri_list = (struct tie_tri_s *)bu_realloc(tie->tri_list, sizeof(struct tie_tri_s) * (tie->tri_num + tnum), "tri_list during tie_push");
453  tie->tri_num_alloc += tnum;
454  }
455 
456  for (i = 0; i < tnum; i++) {
457 
458  if (tie_check_degenerate) {
459  TIE_3 u, v, w;
460 
461  VSUB2(u.v, (*tlist[i*3+1]).v, (*tlist[i*3+0]).v);
462  VSUB2(v.v, (*tlist[i*3+2]).v, (*tlist[i*3+0]).v);
463  VCROSS(w.v, u.v, v.v);
464 
465  if (MAGNITUDE(w.v) < 0.0001 * 0.0001) {
466  bu_log("WARNING: degenerate triangle found: %f %f %f | %f %f %f | %f %f %f\n",
467  V3ARGS((*tlist[i*3+0]).v), V3ARGS((*tlist[i*3+1]).v), V3ARGS((*tlist[i*3+2]).v));
468  continue;
469  }
470  }
471 
472  /* pack pack pack */
473  tie->tri_list[tie->tri_num].data[0] = *tlist[i*3+0];
474  tie->tri_list[tie->tri_num].data[1] = *tlist[i*3+1];
475  tie->tri_list[tie->tri_num].data[2] = *tlist[i*3+2];
476 
477  /* set the association pointer */
478  tie->tri_list[tie->tri_num].ptr = plist;
479  if (plist)
480  plist = (void *)((intptr_t)plist + pstride);
481 
482  V2SETALL(tie->tri_list[tie->tri_num].v, 0.0);
483  tie->tri_list[tie->tri_num].b = 0;
484  tie->tri_num++;
485  }
486  return;
487 }
488 
489 
490 /*
491  * Local Variables:
492  * tab-width: 8
493  * mode: C
494  * indent-tabs-mode: t
495  * c-file-style: "stroustrup"
496  * End:
497  * ex: shiftwidth=4 tabstop=8
498  */
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
struct tie_tri_s ** tri_list
Definition: tieprivate.h:42
int tie_check_degenerate
Definition: btg.c:40
#define TIE_TAB1
Definition: tie.c:44
void TIE_VAL() tie_kdtree_free(struct tie_s *tie)
Definition: tie_kdtree.c:694
#define SMALL_FASTF
Definition: defines.h:342
tfloat far
Definition: tieprivate.h:54
struct tie_kdtree_s * node
Definition: tieprivate.h:52
COMPLEX data[64]
Definition: fftest.c:34
#define TIE_HAS_CHILDREN(bits)
Definition: tieprivate.h:38
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
uint32_t tri_num
Definition: tieprivate.h:43
#define V3ARGS(a)
Definition: color.c:56
void * bu_realloc(void *ptr, size_t siz, const char *str)
#define ZERO(val)
Definition: units.c:38
void TIE_VAL() tie_push(struct tie_s *tie, TIE_3 **tlist, unsigned int tnum, void *plist, unsigned int pstride)
Definition: tie.c:446
void TIE_VAL() tie_prep(struct tie_s *tie)
Definition: tie.c:180
void TIE_VAL() tie_free(struct tie_s *tie)
Definition: tie.c:160
tfloat near
Definition: tieprivate.h:53
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
double fastf_t
Definition: defines.h:300
void TIE_VAL() tie_init(struct tie_s *tie, unsigned int tri_num, unsigned int kdmethod)
Definition: tie.c:137
fastf_t TIE_PREC
Definition: btg.c:41
void *TIE_VAL() tie_work(struct tie_s *tie, struct tie_ray_s *ray, struct tie_id_s *id, void *(*hitfunc)(struct tie_ray_s *, struct tie_id_s *, struct tie_tri_s *, void *ptr), void *ptr)
Definition: tie.c:217
#define UNLIKELY(expression)
Definition: common.h:282
void TIE_VAL() tie_kdtree_prep(struct tie_s *tie)
Definition: tie_kdtree.c:705