BRL-CAD
tie_kdtree.c
Go to the documentation of this file.
1 /* T I E _ K D T R E 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_kdtree.c
21  * tie_kdtree.c
22  *
23  */
24 
25 #include "tie.h"
26 
27 #include <math.h>
28 #include <stdlib.h>
29 #include <string.h>
30 #ifdef HAVE_STDINT_H
31 # include <stdint.h>
32 #endif
33 
34 #include "bio.h"
35 
36 #include "vmath.h"
37 #include "rtgeom.h"
38 #include "raytrace.h"
39 
40 #include "tieprivate.h"
41 
42 #define MAX_SLICES 100
43 #define MIN_SLICES 35
44 #define MIN_DENSITY 0.01
45 #define MIN_SPAN 0.15
46 #define SCALE_COEF 1.80
47 
48 #define TIE_KDTREE_NODE_MAX 4 /* Maximum number of triangles that can reside in a given node until it should be split */
49 #define TIE_KDTREE_DEPTH_K1 1.4 /* K1 Depth Constant Coefficient */
50 #define TIE_KDTREE_DEPTH_K2 1 /* K2 Constant */
51 
52 #define _MIN(a, b) (a)<(b)?(a):(b)
53 #define _MAX(a, b) (a)>(b)?(a):(b)
54 #define MATH_MIN3(_a, _b, _c, _d) _a = _MIN((_b), _MIN((_c), (_d)))
55 #define MATH_MAX3(_a, _b, _c, _d) _a = _MAX((_b), _MAX((_c), (_d)))
56 
57 #define MATH_BBOX(_a, _b, _c, _d, _e) { \
58  MATH_MIN3(_a.v[0], _c.v[0], _d.v[0], _e.v[0]); \
59  MATH_MIN3(_a.v[1], _c.v[1], _d.v[1], _e.v[1]); \
60  MATH_MIN3(_a.v[2], _c.v[2], _d.v[2], _e.v[2]); \
61  MATH_MAX3(_b.v[0], _c.v[0], _d.v[0], _e.v[0]); \
62  MATH_MAX3(_b.v[1], _c.v[1], _d.v[1], _e.v[1]); \
63  MATH_MAX3(_b.v[2], _c.v[2], _d.v[2], _e.v[2]); }
64 
65 /* ======================== X-tests ======================== */
66 #define AXISTEST_X01(a, b, fa, fb) \
67  p.v[0] = a*v0.v[1] - b*v0.v[2]; \
68  p.v[2] = a*v2.v[1] - b*v2.v[2]; \
69  if (p.v[0] < p.v[2]) { min = p.v[0]; max = p.v[2]; } else { min = p.v[2]; max = p.v[0]; } \
70  rad = fa * half_size->v[1] + fb * half_size->v[2]; \
71  if (min > rad || max < -rad) return 0; \
72 
73 #define AXISTEST_X2(a, b, fa, fb) \
74  p.v[0] = a*v0.v[1] - b*v0.v[2]; \
75  p.v[1] = a*v1.v[1] - b*v1.v[2]; \
76  if (p.v[0] < p.v[1]) { min = p.v[0]; max = p.v[1]; } else { min = p.v[1]; max = p.v[0]; } \
77  rad = fa * half_size->v[1] + fb * half_size->v[2]; \
78  if (min > rad || max < -rad) return 0;
79 
80 /* ======================== Y-tests ======================== */
81 #define AXISTEST_Y02(a, b, fa, fb) \
82  p.v[0] = -a*v0.v[0] + b*v0.v[2]; \
83  p.v[2] = -a*v2.v[0] + b*v2.v[2]; \
84  if (p.v[0] < p.v[2]) { min = p.v[0]; max = p.v[2]; } else { min = p.v[2]; max = p.v[0]; } \
85  rad = fa * half_size->v[0] + fb * half_size->v[2]; \
86  if (min > rad || max < -rad) return 0;
87 
88 #define AXISTEST_Y1(a, b, fa, fb) \
89  p.v[0] = -a*v0.v[0] + b*v0.v[2]; \
90  p.v[1] = -a*v1.v[0] + b*v1.v[2]; \
91  if (p.v[0] < p.v[1]) { min = p.v[0]; max = p.v[1]; } else { min = p.v[1]; max = p.v[0]; } \
92  rad = fa * half_size->v[0] + fb * half_size->v[2]; \
93  if (min > rad || max < -rad) return 0;
94 
95 /* ======================== Z-tests ======================== */
96 #define AXISTEST_Z12(a, b, fa, fb) \
97  p.v[1] = a*v1.v[0] - b*v1.v[1]; \
98  p.v[2] = a*v2.v[0] - b*v2.v[1]; \
99  if (p.v[2] < p.v[1]) { min = p.v[2]; max = p.v[1]; } else { min = p.v[1]; max = p.v[2]; } \
100  rad = fa * half_size->v[0] + fb * half_size->v[1]; \
101  if (min > rad || max < -rad) return 0;
102 
103 #define AXISTEST_Z0(a, b, fa, fb) \
104  p.v[0] = a*v0.v[0] - b*v0.v[1]; \
105  p.v[1] = a*v1.v[0] - b*v1.v[1]; \
106  if (p.v[0] < p.v[1]) { min = p.v[0]; max = p.v[1]; } else { min = p.v[1]; max = p.v[0]; } \
107  rad = fa * half_size->v[0] + fb * half_size->v[1]; \
108  if (min > rad || max < -rad) return 0;
109 
110 
111 /*************************************************************
112  **************** PRIVATE FUNCTIONS **************************
113  *************************************************************/
114 
115 
116 static void
117 tie_kdtree_free_node(struct tie_kdtree_s *node)
118 {
119  if (node && TIE_HAS_CHILDREN(node->b)) {
120  /* Node Data is KDTREE Children, Recurse */
121  tie_kdtree_free_node(&((struct tie_kdtree_s *)(node->data))[0]);
122  tie_kdtree_free_node(&((struct tie_kdtree_s *)(node->data))[1]);
123  } else {
124  /* This node points to a geometry node, free it */
125  struct tie_geom_s *tmp;
126  tmp = (struct tie_geom_s *)node->data;
127  if (tmp) {
128  if (tmp->tri_num > 0) {
129  bu_free(tmp->tri_list, "tri_list");
130  }
131  bu_free(tmp, "data");
132  }
133  }
134 }
135 
136 
137 static void
138 tie_kdtree_prep_head(struct tie_s *tie, struct tie_tri_s *tri_list, unsigned int tri_num)
139 {
140  struct tie_geom_s *g;
141  TIE_3 min, max;
142  unsigned int i;
143 
144  if (tri_num == 0 || tie->kdtree != NULL)
145  return;
146 
147  /* Insert all triangles into the Head Node */
148  BU_ALLOC(tie->kdtree, struct tie_kdtree_s);
149  BU_ALLOC(tie->kdtree->data, struct tie_geom_s);
150 
151  g = ((struct tie_geom_s *)(tie->kdtree->data));
152  g->tri_num = 0;
153  VSETALL(tie->min, +INFINITY);
154  VSETALL(tie->max, -INFINITY);
155 
156  g->tri_list = (struct tie_tri_s **)bu_calloc(tri_num, sizeof(struct tie_tri_s *), __FUNCTION__);
157 
158  /* form bounding box of scene */
159  for (i = 0; i < tri_num; i++) {
160  g->tri_list[i] = &tri_list[i];
161 
162  /* Get Bounding Box of Triangle */
163  MATH_BBOX(min, max, tri_list[i].data[0], tri_list[i].data[1], tri_list[i].data[2]);
164  /* Check to see if defines a new Max or Min point */
165  VMIN(tie->min, min.v);
166  VMAX(tie->max, max.v);
167  }
168  VADD2SCALE(tie->mid, tie->min, tie->max, 0.5);
169  tie->radius = DIST_PT_PT(tie->max, tie->mid);
170  g->tri_num = tri_num;
171 }
172 
173 static int
174 tie_kdtree_tri_box_overlap(TIE_3 *center, TIE_3 *half_size, TIE_3 *triverts)
175 {
176  /*
177  * use separating axis theorem to test overlap between triangle and box
178  * need to test for overlap in these directions:
179  * 1) the {x, y, z}-directions (actually, since we use the AABB of the triangle
180  * we do not even need to test these)
181  * 2) normal of the triangle
182  * 3) crossproduct(edge from tri, {x, y, z}-direction)
183  * this gives 3x3=9 more tests
184  */
185  TIE_3 v0, v1, v2, normal, e0, e1, e2, fe, p;
186  tfloat min, max, d, t, rad;
187 
188  /* move everything so that the boxcenter is in (0, 0, 0) */
189  VSUB2(v0.v, triverts[0].v, (*center).v);
190  VSUB2(v1.v, triverts[1].v, (*center).v);
191  VSUB2(v2.v, triverts[2].v, (*center).v);
192 
193  /*
194  * First test overlap in the {x, y, z}-directions
195  * find min, max of the triangle each direction, and test for overlap in
196  * that direction -- this is equivalent to testing a minimal AABB around
197  * the triangle against the AABB
198  */
199 
200  /* test in X-direction */
201  MATH_MIN3(min, v0.v[0], v1.v[0], v2.v[0]);
202  MATH_MAX3(max, v0.v[0], v1.v[0], v2.v[0]);
203  if (min > half_size->v[0] || max < -half_size->v[0])
204  return 0;
205 
206  /* test in Y-direction */
207  MATH_MIN3(min, v0.v[1], v1.v[1], v2.v[1]);
208  MATH_MAX3(max, v0.v[1], v1.v[1], v2.v[1]);
209  if (min > half_size->v[1] || max < -half_size->v[1])
210  return 0;
211 
212  /* test in Z-direction */
213  MATH_MIN3(min, v0.v[2], v1.v[2], v2.v[2]);
214  MATH_MAX3(max, v0.v[2], v1.v[2], v2.v[2]);
215  if (min > half_size->v[2] || max < -half_size->v[2])
216  return 0;
217 
218  /* compute triangle edges */
219  VSUB2(e0.v, v1.v, v0.v); /* tri edge 0 */
220  VSUB2(e1.v, v2.v, v1.v); /* tri edge 1 */
221  VSUB2(e2.v, v0.v, v2.v); /* tri edge 2 */
222 
223  /* Perform the 9 tests */
224  fe.v[0] = fabs(e0.v[0]);
225  fe.v[1] = fabs(e0.v[1]);
226  fe.v[2] = fabs(e0.v[2]);
227 
228  AXISTEST_X01(e0.v[2], e0.v[1], fe.v[2], fe.v[1]);
229  AXISTEST_Y02(e0.v[2], e0.v[0], fe.v[2], fe.v[0]);
230  AXISTEST_Z12(e0.v[1], e0.v[0], fe.v[1], fe.v[0]);
231 
232  fe.v[0] = fabs(e1.v[0]);
233  fe.v[1] = fabs(e1.v[1]);
234  fe.v[2] = fabs(e1.v[2]);
235 
236  AXISTEST_X01(e1.v[2], e1.v[1], fe.v[2], fe.v[1]);
237  AXISTEST_Y02(e1.v[2], e1.v[0], fe.v[2], fe.v[0]);
238  AXISTEST_Z0(e1.v[1], e1.v[0], fe.v[1], fe.v[0]);
239 
240  fe.v[0] = fabs(e2.v[0]);
241  fe.v[1] = fabs(e2.v[1]);
242  fe.v[2] = fabs(e2.v[2]);
243 
244  AXISTEST_X2(e2.v[2], e2.v[1], fe.v[2], fe.v[1]);
245  AXISTEST_Y1(e2.v[2], e2.v[0], fe.v[2], fe.v[0]);
246  AXISTEST_Z12(e2.v[1], e2.v[0], fe.v[1], fe.v[0]);
247 
248  /*
249  * Test if the box intersects the plane of the triangle
250  * compute plane equation of triangle: normal*x+d=0
251  */
252  VCROSS(normal.v, e0.v, e1.v);
253  d = VDOT( normal.v, v0.v); /* plane eq: normal . x + d = 0 */
254 
255  p.v[0] = normal.v[0] > 0 ? half_size->v[0] : -half_size->v[0];
256  p.v[1] = normal.v[1] > 0 ? half_size->v[1] : -half_size->v[1];
257  p.v[2] = normal.v[2] > 0 ? half_size->v[2] : -half_size->v[2];
258 
259  t = VDOT( normal.v, p.v);
260  return t >= d ? 1 : 0;
261 }
262 
263 static unsigned int
264 find_split_fast(struct tie_kdtree_s *node, TIE_3 *cmin, TIE_3 *cmax)
265 {
266  /**********************
267  * MID-SPLIT ALGORITHM *
268  ***********************/
269  TIE_3 vec, center[2];
270 
271  VADD2(center[0].v, cmax[0].v, cmin[0].v);
272  VSCALE(center[0].v, center[0].v, 0.5);
273 
274  /* Split along largest Axis to keep node sizes relatively cube-like (Naive) */
275  VSUB2(vec.v, cmax[0].v, cmin[0].v);
276 
277  /* Determine the largest Axis */
278  if (vec.v[0] >= vec.v[1] && vec.v[0] >= vec.v[2]) {
279  cmax[0].v[0] = center[0].v[0];
280  cmin[1].v[0] = center[0].v[0];
281  node->axis = center[0].v[0];
282  return 0;
283  }
284 
285  if (vec.v[1] >= vec.v[0] && vec.v[1] >= vec.v[2]) {
286  cmax[0].v[1] = center[0].v[1];
287  cmin[1].v[1] = center[0].v[1];
288  node->axis = center[0].v[1];
289  return 1;
290  }
291 
292  cmax[0].v[2] = center[0].v[2];
293  cmin[1].v[2] = center[0].v[2];
294  node->axis = center[0].v[2];
295  return 2;
296 }
297 
298 static unsigned int
299 find_split_optimal(struct tie_s *tie, struct tie_kdtree_s *node, TIE_3 *cmin, TIE_3 *cmax)
300 {
301  /****************************************
302  * Justin's Aggressive KD-Tree Algorithm *
303  *****************************************/
304  unsigned int slice[3][MAX_SLICES+MIN_SLICES], gap[3][2], active, split_slice = 0, split;
305  unsigned int side[3][MAX_SLICES+MIN_SLICES][2], i, j, d, s, n, k, smax[3], smin, slice_num;
306  tfloat coef[3][MAX_SLICES+MIN_SLICES], split_coef, beg, end, d_min = 0.0, d_max = 0.0;
307  struct tie_tri_s *tri;
308  struct tie_geom_s *node_gd = (struct tie_geom_s *)(node->data);
309  TIE_3 min, max;
310  TIE_3 center[2], half_size[2];
311  VSETALL(min.v, 0.0);
312  VSETALL(max.v, 0.0);
313 
314  /*
315  * Calculate number of slices to use as a function of triangle density.
316  * Setting slices as a function of relative node size does not work so well.
317  */
318  /* slice_num = MIN_SLICES + MAX_SLICES * ((tfloat)node_gd->tri_num / (tfloat)tie->tri_num); */
319  slice_num = 1*node_gd->tri_num > MAX_SLICES ? MAX_SLICES : 1*node_gd->tri_num;
320 
321  for (d = 0; d < 3; d++) {
322  /*
323  * Optimization: Walk each triangle and find the min and max for the given dimension
324  * of the complete triangle list. This will tell us what slices we needn't bother
325  * doing any computations for.
326  */
327  for (i = 0; i < node_gd->tri_num; i++) {
328  tri = node_gd->tri_list[i];
329  /* Set min anx max */
330  MATH_MIN3(tri->v[0], tri->data[0].v[d], tri->data[1].v[d], tri->data[2].v[d]);
331  MATH_MAX3(tri->v[1], tri->data[0].v[d], tri->data[1].v[d], tri->data[2].v[d]);
332 
333  /* Clamp to node AABB */
334  if (tri->v[0] < min.v[d])
335  tri->v[0] = min.v[d];
336  if (tri->v[1] > max.v[d])
337  tri->v[1] = max.v[d];
338 
339  if (i == 0 || tri->v[0] < d_min)
340  d_min = tri->v[0];
341 
342  if (i == 0 || tri->v[1] > d_max)
343  d_max = tri->v[1];
344  }
345 
346  for (k = 0; k < slice_num; k++) {
347  slice[d][k] = 0;
348  side[d][k][0] = 0;
349  side[d][k][1] = 0;
350 
351  /* Left Child */
352  cmin[0] = min;
353  cmax[0] = max;
354 
355  /* Right Child */
356  cmin[1] = min;
357  cmax[1] = max;
358 
359  /* construct slices so as not to use the boundaries as slices */
360  coef[d][k] = ((tfloat)k / (tfloat)(slice_num-1)) * (tfloat)(slice_num-2) / (tfloat)slice_num + (tfloat)1 / (tfloat)slice_num;
361  cmax[0].v[d] = min.v[d]*(1.0-coef[d][k]) + max.v[d]*coef[d][k];
362  cmin[1].v[d] = cmax[0].v[d];
363 
364  /* determine whether to bother with this slice */
365  if (cmax[0].v[d] < d_min || cmax[0].v[d] > d_max)
366  continue;
367 
368  for (n = 0; n < 2; n++) {
369  VADD2(center[n].v, cmax[n].v, cmin[n].v);
370  VSCALE(center[n].v, center[n].v, 0.5);
371  VSUB2(half_size[n].v, cmax[n].v, cmin[n].v);
372  VSCALE(half_size[n].v, half_size[n].v, 0.5);
373  }
374 
375  for (i = 0; i < node_gd->tri_num; i++) {
376  /*
377  * Optimization: If the points for the triangle of the dimension being tested
378  * do not span the cutting plane, then do not bother with the next test.
379  */
380  if ((node_gd->tri_list[i]->data[0].v[d] > cmax[0].v[d] &&
381  node_gd->tri_list[i]->data[1].v[d] > cmax[0].v[d] &&
382  node_gd->tri_list[i]->data[2].v[d] > cmax[0].v[d])||
383  (node_gd->tri_list[i]->data[0].v[d] < cmax[0].v[d] &&
384  node_gd->tri_list[i]->data[1].v[d] < cmax[0].v[d] &&
385  node_gd->tri_list[i]->data[2].v[d] < cmax[0].v[d]))
386  continue;
387 
388  /* Check that the triangle is in both node A and B for it to span. */
389  s = 0;
390  for (n = 0; n < 2; n++) {
391  /*
392  * Check to see if any triangle points are inside of the node before
393  * spending a lot of cycles on the full blown triangle box overlap
394  */
395  for (j = 0; j < 3; j++)
396  if (node_gd->tri_list[i]->data[j].v[0] > cmin[n].v[0] &&
397  node_gd->tri_list[i]->data[j].v[0] < cmax[n].v[0] &&
398  node_gd->tri_list[i]->data[j].v[1] > cmin[n].v[1] &&
399  node_gd->tri_list[i]->data[j].v[1] < cmax[n].v[1] &&
400  node_gd->tri_list[i]->data[j].v[2] > cmin[n].v[2] &&
401  node_gd->tri_list[i]->data[j].v[2] < cmax[n].v[2]) {
402  j = 4;
403  }
404 
405  if (j == 5) {
406  s++;
407  side[d][k][n]++;
408  } else {
409  if (tie_kdtree_tri_box_overlap(&center[n], &half_size[n], node_gd->tri_list[i]->data)) {
410  s++;
411  side[d][k][n]++;
412  }
413  }
414  }
415 
416  if (s == 2)
417  slice[d][k]++;
418  }
419  }
420  }
421 
422  /* Store the max value from each of the 3 Slice arrays */
423  for (d = 0; d < 3; d++) {
424  smax[d] = 0;
425  for (k = 0; k < slice_num; k++)
426  if (slice[d][k] > smax[d])
427  smax[d] = slice[d][k];
428  }
429 
430  /*
431  * To eliminate "empty" areas, build a list of spans where geometric complexity is
432  * less than MIN_SPAN of the overall nodes size and then selecting the splitting plane
433  * that corresponds to the span slice domain nearest the center to bias towards a balanced tree
434  */
435 
436  for (d = 0; d < 3; d++) {
437  gap[d][0] = 0;
438  gap[d][1] = 0;
439  beg = 0;
440  end = 0;
441  active = 0;
442 
443  for (k = 0; k < slice_num; k++) {
444  if (slice[d][k] < (unsigned int)(MIN_DENSITY * (tfloat)smax[d])) {
445  if (!active) {
446  active = 1;
447  beg = k;
448  end = k;
449  } else
450  end = k;
451  } else {
452  if (active)
453  if (end - beg > gap[d][1] - gap[d][0]) {
454  gap[d][0] = beg;
455  gap[d][1] = end;
456  }
457  active = 0;
458  }
459  }
460 
461  if (active)
462  if (end - beg > gap[d][1] - gap[d][0]) {
463  gap[d][0] = beg;
464  gap[d][1] = end;
465  }
466  }
467 
468  /*
469  * If there is a gap at least MIN_SPAN in side w/r/t the nodes dimension size
470  * then use the nearest edge of the gap to 0.5 as the splitting plane.
471  * Use the gap with the largest span.
472  * If no gaps are found meeting the criteria then weight the span values to
473  * bias towards a balanced kd-tree and choose the minima of that weighted curve.
474  */
475 
476  /* Largest gap */
477  d = 0;
478  if (gap[1][1] - gap[1][0] > gap[d][1] - gap[d][0])
479  d = 1;
480  if (gap[2][1] - gap[2][0] > gap[d][1] - gap[d][0])
481  d = 2;
482 
483  /*
484  * Largest gap found must meet MIN_SPAN requirements
485  * There must be at least 500 triangles or we don't bother.
486  * Lower triangle numbers means there is a higher probability that
487  * triangles lack any sort of coherent structure.
488  */
489  if ((tfloat)(gap[d][1] - gap[d][0]) / (tfloat)slice_num > MIN_SPAN && node_gd->tri_num > 500) {
490  split = d;
491  if (abs(gap[d][0] - slice_num/2) < abs(gap[d][1] - slice_num/2)) {
492  /* choose gap[d][0] as splitting plane */
493  split_coef = ((tfloat)gap[d][0] / (tfloat)(slice_num-1)) * (tfloat)(slice_num-2) / (tfloat)slice_num + (tfloat)1 / (tfloat)slice_num;
494  split_slice = gap[d][0];
495  } else {
496  /* choose gap[d][1] as splitting plane */
497  split_coef = ((tfloat)gap[d][1] / (tfloat)(slice_num-1)) * (tfloat)(slice_num-2) / (tfloat)slice_num + (tfloat)1 / (tfloat)slice_num;
498  split_slice = gap[d][1];
499  }
500  } else {
501  /*
502  * Weight the slices based on a heuristic driven linear scaling function to bias values
503  * towards the center as more desirable. This solves the case of a partially linear graph
504  * to prevent marching in order to determine a desirable splitting point. If this section of code
505  * is being executed it's typically because most 'empty space' has now been eliminated
506  * and/or the resulting geometry is now losing structure as the smaller cells are being
507  * created, i.e. dividing a fraction of a wing-nut instead of an engine-block.
508  */
509  for (d = 0; d < 3; d++)
510  for (k = 0; k < slice_num; k++)
511  slice[d][k] += fabs(coef[d][k]-0.5) * SCALE_COEF * smax[d];
512 
513  /* Choose the slice with the graphs minima as the splitting plane. */
514  split = 0;
515  smin = tie->tri_num;
516  split_coef = 0.5;
517  for (d = 0; d < 3; d++) {
518  for (k = 0; k < slice_num; k++) {
519  if (slice[d][k] < smin) {
520  split_coef = coef[d][k];
521  split = d;
522  split_slice = k;
523  smin = slice[d][k];
524  }
525  }
526  }
527 
528  /*
529  * If the dimension chosen to split along has a value of 0 for the maximum value
530  * then the geometry was aligned such that it fell undetectable between the slices
531  * and therefore was not picked up by the marching slices. In the event that this
532  * happens, choose to naively split along the middle as this last ditch decision
533  * will give better results than the algorithm naively picking the first of the
534  * the slices forming these irregular, short followed by a long box, splits.
535  */
536  if (smax[split] == 0) {
537  split_slice = slice_num / 2;
538  split_coef = coef[split][split_slice];
539  }
540  }
541 
542  /*
543  * Lastly, after we have supposedly chosen the ideal splitting point,
544  * check to see that the subdivision that is about to take place is worth
545  * doing. In other words, if one of the children have the same number of triangles
546  * as the parent does then stop.
547  */
548  if (side[split][split_slice][0] == node_gd->tri_num || side[split][split_slice][1] == node_gd->tri_num) {
549  tie->stat += node_gd->tri_num;
550  return split;
551  }
552 
553  /* Based on the winner, construct the two child nodes */
554  VMOVE(cmin[0].v, min.v);
555  VMOVE(cmax[0].v, max.v);
556 
557  VMOVE(cmin[1].v, min.v);
558  VMOVE(cmax[1].v, max.v);
559 
560  cmax[0].v[split] = min.v[split]*(1.0-split_coef) + max.v[split]*split_coef;
561  cmin[1].v[split] = cmax[0].v[split];
562  node->axis = cmax[0].v[split];
563  return split;
564 }
565 
566 static void
567 tie_kdtree_build(struct tie_s *tie, struct tie_kdtree_s *node, unsigned int depth, point_t min, point_t max)
568 {
569  struct tie_geom_s *child[2], *node_gd = (struct tie_geom_s *)(node->data);
570  TIE_3 cmin[2], cmax[2], center[2], half_size[2];
571  unsigned int i, j, n, split = 0, cnt[2];
572 
573  if (node_gd == NULL) {
574  bu_log("null geom, aborting\n");
575  return;
576  }
577 
578  /* initialize cmax to make the compiler happy */
579  VMOVE(cmax[0].v, max);
580  VMOVE(cmin[0].v, min);
581  VMOVE(cmax[1].v, max);
582  VMOVE(cmin[1].v, min);
583 
584  /* Terminating criteria for KDTREE subdivision */
585  if (node_gd->tri_num <= TIE_KDTREE_NODE_MAX || depth > tie->max_depth) {
586  tie->stat += node_gd->tri_num;
587  return;
588  }
589 
590  if (tie->kdmethod == TIE_KDTREE_FAST)
591  split = find_split_fast(node, &cmin[0], &cmax[0]);
592  else if (tie->kdmethod == TIE_KDTREE_OPTIMAL)
593  split = find_split_optimal(tie, node, &cmin[0], &cmax[0]);
594  else
595  bu_bomb("Illegal tie kdtree method\n");
596 
597  /* Allocate 2 children nodes for the parent node */
598  node->data = bu_calloc(2, sizeof(struct tie_kdtree_s), __FUNCTION__);
599  node->b = 0;
600 
601  BU_ALLOC(((struct tie_kdtree_s *)(node->data))[0].data, struct tie_geom_s);
602  BU_ALLOC(((struct tie_kdtree_s *)(node->data))[1].data, struct tie_geom_s);
603 
604  /* Initialize Triangle List */
605  child[0] = ((struct tie_geom_s *)(((struct tie_kdtree_s *)(node->data))[0].data));
606  child[1] = ((struct tie_geom_s *)(((struct tie_kdtree_s *)(node->data))[1].data));
607 
608  child[0]->tri_list = (struct tie_tri_s **)bu_calloc(node_gd->tri_num, sizeof(struct tie_tri_s *), "child[0]->tri_list");
609  child[0]->tri_num = 0;
610 
611  child[1]->tri_list = (struct tie_tri_s **)bu_calloc(node_gd->tri_num, sizeof(struct tie_tri_s *), "child[1]->tri_list");
612  child[1]->tri_num = 0;
613 
614 
615  /*
616  * Determine if the triangles touch either of the two children nodes,
617  * if it does insert it into them respectively.
618  */
619  for (n = 0; n < 2; n++) {
620  cnt[n] = 0;
621 
622  VADD2(center[n].v, cmax[n].v, cmin[n].v);
623  VSCALE(center[n].v, center[n].v, 0.5);
624  VSUB2(half_size[n].v, cmax[n].v, cmin[n].v);
625  VSCALE(half_size[n].v, half_size[n].v, 0.5);
626 
627  for (i = 0; i < node_gd->tri_num; i++) {
628  /*
629  * Check to see if any triangle points are inside of the node before
630  * spending a lot of cycles on the full blown triangle box overlap
631  */
632  for (j = 0; j < 3; j++)
633  if (node_gd->tri_list[i]->data[j].v[0] > cmin[n].v[0] &&
634  node_gd->tri_list[i]->data[j].v[0] < cmax[n].v[0] &&
635  node_gd->tri_list[i]->data[j].v[1] > cmin[n].v[1] &&
636  node_gd->tri_list[i]->data[j].v[1] < cmax[n].v[1] &&
637  node_gd->tri_list[i]->data[j].v[2] > cmin[n].v[2] &&
638  node_gd->tri_list[i]->data[j].v[2] < cmax[n].v[2]) {
639  j = 4;
640  }
641 
642  if (j == 5) {
643  child[n]->tri_list[child[n]->tri_num++] = node_gd->tri_list[i];
644  cnt[n]++;
645  } else
646  if (tie_kdtree_tri_box_overlap(&center[n], &half_size[n], node_gd->tri_list[i]->data)) {
647  child[n]->tri_list[child[n]->tri_num++] = node_gd->tri_list[i];
648  cnt[n]++;
649  }
650  }
651 
652  /* Resize Tri List to actual amount of memory used */
653  /* TODO: examine if this is correct. A 0 re-alloc is probably a very bad
654  * thing. */
655  if ( child[n]->tri_num == 0 ) {
656  if ( child[n]->tri_list ) {
657  bu_free( child[n]->tri_list, "tri_list");
658  child[n]->tri_list = NULL;
659  }
660  } else {
661  child[n]->tri_list = (struct tie_tri_s **)bu_realloc(child[n]->tri_list, sizeof(struct tie_tri_s *)*child[n]->tri_num, __FUNCTION__);
662  }
663  }
664 
665  /*
666  * Now that the triangles have been propagated to the appropriate child nodes,
667  * free the triangle list on this node.
668  */
669  node_gd->tri_num = 0;
670  bu_free(node_gd->tri_list, "tri_list");
671  bu_free(node_gd, "node");
672 
673  /* Push each child through the same process. */
674  {
675  point_t lmin[2], lmax[2];
676  VMOVE(lmin[0], cmin[0].v);
677  VMOVE(lmin[1], cmin[1].v);
678  VMOVE(lmax[0], cmax[0].v);
679  VMOVE(lmax[1], cmax[1].v);
680  tie_kdtree_build(tie, &((struct tie_kdtree_s *)(node->data))[0], depth+1, lmin[0], lmax[0]);
681  tie_kdtree_build(tie, &((struct tie_kdtree_s *)(node->data))[1], depth+1, lmin[1], lmax[1]);
682  }
683 
684  /* Assign the splitting dimension to the node */
685  /* If we've come this far then YES, this node DOES have child nodes, MARK it as so. */
686  node->b = TIE_SET_HAS_CHILDREN(node->b) + split;
687 }
688 
689 /*************************************************************
690  **************** EXPORTED FUNCTIONS *************************
691  *************************************************************/
692 
693 void
694 TIE_VAL(tie_kdtree_free)(struct tie_s *tie)
695 {
696  /* Free KDTREE Nodes */
697  /* prevent tie from crashing when a tie_free() is called right after a tie_init() */
698  if (tie->kdtree) {
699  tie_kdtree_free_node(tie->kdtree);
700  bu_free(tie->kdtree, "kdtree");
701  }
702 }
703 
704 void
705 TIE_VAL(tie_kdtree_prep)(struct tie_s *tie)
706 {
707  TIE_3 delta;
708  int already_built;
709  struct tie_geom_s *g;
710 
711  already_built = tie->kdtree ? 1 : 0;
712 
713  /* Set bounding volume and make head node a geometry node */
714  if (!already_built)
715  tie_kdtree_prep_head(tie, tie->tri_list, tie->tri_num);
716 
717  if (!tie->kdtree)
718  return;
719 
720  g = (struct tie_geom_s *)tie->kdtree->data;
721 
722  /* Trim KDTREE to number of actual triangles if it's not that size already. */
723  /* TODO: examine if this is correct. A 0 re-alloc is probably a very bad
724  * thing. */
725 
726  if (!already_built) {
727  if (g->tri_num)
728  g->tri_list = (struct tie_tri_s **)bu_realloc(g->tri_list, sizeof(struct tie_tri_s *) * g->tri_num, "prep tri_list");
729  } else {
730  bu_free(g->tri_list, "tri_list");
731  }
732 
733  /*
734  * Compute Floating Fuzz Precision Value
735  * For now, take largest dimension as basis for TIE_PREC
736  */
737  VMOVE(tie->amin, tie->min);
738  VMOVE(tie->amax, tie->max);
739  VSUB2(delta.v, tie->max, tie->min);
740  MATH_MAX3(TIE_PREC, delta.v[0], delta.v[1], delta.v[2]);
741 #if defined(TIE_PRECISION) && defined(TIE_PRECISION_SINGLE) && TIE_PRECISION == TIE_PRECISION_SINGLE
742  TIE_PREC *= 0.000000001;
743 #else
744  TIE_PREC *= 0.000000000001;
745 #endif
746 
747  /* Grow the head node to avoid floating point fuzz in the building process with edges */
748  VSCALE(delta.v, delta.v, 1.0);
749  VSUB2(tie->min, tie->min, delta.v);
750  VADD2(tie->max, tie->max, delta.v);
751 
752  /* Compute Max Depth to allow the KD-Tree to grow to */
753  tie->max_depth = (int)(TIE_KDTREE_DEPTH_K1 * (log(tie->tri_num) / log(2)) + TIE_KDTREE_DEPTH_K2);
754 
755  /* Build the KDTREE */
756  if (!already_built)
757  tie_kdtree_build(tie, tie->kdtree, 0, tie->min, tie->max);
758 
759  tie->stat = 0;
760 }
761 
762 /*
763  * Local Variables:
764  * tab-width: 8
765  * mode: C
766  * indent-tabs-mode: t
767  * c-file-style: "stroustrup"
768  * End:
769  * ex: shiftwidth=4 tabstop=8
770  */
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define AXISTEST_X01(a, b, fa, fb)
Definition: tie_kdtree.c:66
#define AXISTEST_Z12(a, b, fa, fb)
Definition: tie_kdtree.c:96
#define TIE_KDTREE_NODE_MAX
Definition: tie_kdtree.c:48
if lu s
Definition: nmg_mod.c:3860
struct tie_tri_s ** tri_list
Definition: tieprivate.h:42
#define VSETALL(a, s)
Definition: color.c:54
void TIE_VAL() tie_kdtree_free(struct tie_s *tie)
Definition: tie_kdtree.c:694
#define AXISTEST_X2(a, b, fa, fb)
Definition: tie_kdtree.c:73
#define AXISTEST_Y02(a, b, fa, fb)
Definition: tie_kdtree.c:81
if(share_geom)
Definition: nmg_mod.c:3829
COMPLEX data[64]
Definition: fftest.c:34
#define TIE_HAS_CHILDREN(bits)
Definition: tieprivate.h:38
#define TIE_KDTREE_DEPTH_K1
Definition: tie_kdtree.c:49
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
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 MIN_SPAN
Definition: tie_kdtree.c:45
void * bu_realloc(void *ptr, size_t siz, const char *str)
#define MAX_SLICES
Definition: tie_kdtree.c:42
#define MATH_MIN3(_a, _b, _c, _d)
Definition: tie_kdtree.c:54
#define MATH_MAX3(_a, _b, _c, _d)
Definition: tie_kdtree.c:55
#define MATH_BBOX(_a, _b, _c, _d, _e)
Definition: tie_kdtree.c:57
#define TIE_SET_HAS_CHILDREN(bits)
Definition: tieprivate.h:39
#define MIN_SLICES
Definition: tie_kdtree.c:43
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
#define MIN_DENSITY
Definition: tie_kdtree.c:44
#define AXISTEST_Z0(a, b, fa, fb)
Definition: tie_kdtree.c:103
#define SCALE_COEF
Definition: tie_kdtree.c:46
#define TIE_KDTREE_DEPTH_K2
Definition: tie_kdtree.c:50
#define AXISTEST_Y1(a, b, fa, fb)
Definition: tie_kdtree.c:88
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
HIDDEN const point_t delta
Definition: sh_prj.c:618
fastf_t TIE_PREC
Definition: btg.c:41
void TIE_VAL() tie_kdtree_prep(struct tie_s *tie)
Definition: tie_kdtree.c:705