HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
UT_KDTree.h
Go to the documentation of this file.
1 /*
2  * PROPRIETARY INFORMATION. This software is proprietary to
3  * Side Effects Software Inc., and is not to be reproduced,
4  * transmitted, or disclosed in any way without written permission.
5  *
6  * NAME: UT_KDTree.h ( UT Library, C++)
7  *
8  * COMMENTS: KD Tree Implementation
9  *
10  * The template class representing the nodes must have the following methods
11  * defined:
12  * const UT_Vector3/UT_Vector4 getP(int idx) const;
13  * bool isValid(int idx) const;
14  * It's ok for the area function to return 0.
15  * The valid method is used when searching for valid nodes.
16  * The splits are the dimension used which the node is split in. Thus, the
17  * split direction can be encoded using a minimal 2 bits.
18  *
19  * TODO: We may want to convert this to a template class
20  */
21 
22 #ifndef __UT_KDTree__
23 #define __UT_KDTree__
24 
25 #include "UT_API.h"
26 #include "UT_BoundingBox.h"
27 #include "UT_BoundingRect.h"
28 #include "UT_IntArray.h"
29 #include "UT_Array.h"
30 #include "UT_Lock.h"
31 #include "UT_VectorTypes.h"
32 #include <SYS/SYS_Math.h>
33 
34 
35 class ut_KDPQueue;
36 class UT_Filter;
37 
38 // The maximum dimension must be less than 255 (unless you change the splits
39 // variable to a ushort)
40 #define UT_KD_MAXDIM 254
41 
42 /// Lookup point information to be passed to the query functions
43 struct UT_KDQueryPt {
44  UT_KDQueryPt(const float *p)
45  { P = p; myData = 0; }
46  UT_KDQueryPt(const float *p, const void *data)
47  { P = p; myData = data; }
48 
49  float boxDist(int axis, float bmin, float bmax) const
50  {
51  return SYSmax(P[axis] - bmax, bmin - P[axis], 0.0F);
52  }
53  float dist(int axis, float off) const
54  {
55  return P[axis] - off;
56  }
57  float dist(const float *P, int dim) const
58  {
59  float d2 = 0;
60  float d;
61 
62  switch (dim) {
63  default:
64  for (int i = dim; i-- > 4; ) {
65  d = dist(i, P[i]); d2 += d*d;
66  }
67  case 4: d = dist(3, P[3]); d2 += d*d;
68  case 3: d = dist(2, P[2]); d2 += d*d;
69  case 2: d = dist(1, P[1]); d2 += d*d;
70  case 1: d = dist(0, P[0]); d2 += d*d;
71  }
72  return d2;
73  }
74 
75  const float *P;
76  const void *myData;
77 };
78 
79 /// This query point considers space to wrap between 0 and 1
80 /// in all dimensions. It only supports up to 3 dimensions
81 /// due to UT_BoundingBox only having 3 dimensions.
83  UT_KDQueryPtUnitWrap(const float *p, int ndims=3)
84  : myP(p)
85  , myNDims(ndims)
86  {}
87 
88  /// This can be an underestimate, but not an overestimate
89  /// of the distance squared.
90  float boxDist(const UT_BoundingBox &box) const
91  {
92  float maxdist = 0;
93  for (int axis = 0; axis < myNDims; ++axis)
94  {
95  float p = myP[axis];
96  float c = box.centerAxis(axis);
97  float centredist;
98  if (c < p)
99  centredist = SYSmin(p - c, 1 - p + c);
100  else
101  centredist = SYSmin(c - p, 1 - c + p);
102  float axisdist = centredist - 0.5f*box.sizeAxis(axis);
103  maxdist = SYSmax(maxdist, axisdist);
104  }
105  return ((maxdist >= 0) ? 1 : -1) * maxdist*maxdist;
106  }
107 
108  /// This distance squared needs to be exact.
109  float dist(const float *P, int) const
110  {
111  float d2 = 0;
112  for (int axis = 0; axis < myNDims; ++axis)
113  {
114  float p = myP[axis];
115  float q = P[axis];
116  float axisdist;
117  if (q < p)
118  axisdist = SYSmin(p - q, 1 - p + q);
119  else
120  axisdist = SYSmin(q - p, 1 - q + p);
121 
122  d2 += axisdist*axisdist;
123  }
124  return d2;
125  }
126 
127  const float *myP;
128  const int myNDims;
129 };
130 
131 /// Lookup information for 2D-tree triangle queries
132 /// NOTE: Distance squared here is not Euclidean distance squared, but
133 /// edge-perpendicular distance squared, i.e. distance is from
134 /// the incentre out, perpendicular to one of the edges, minus the
135 /// incircle's radius. This avoids the need to have special
136 /// cases for what would be the circular sections around each vertex.
137 /// It basically indicates how far the triangle would have to be
138 /// expanded (or contracted) relative to the incentre in order to
139 /// hit the query point.
141 public:
143  const UT_Vector2 &a,
144  const UT_Vector2 &b,
145  const UT_Vector2 &c)
146  {
147  UT_Vector2 sidea = c-b;
148  UT_Vector2 sideb = a-c;
149  UT_Vector2 sidec = b-a;
150  float la = sidea.length();
151  float lb = sideb.length();
152  float lc = sidec.length();
153  float p = la + lb + lc;
154  myIncentre = (la*a + lb*b + lc*c)/p;
155  float s = 0.5f*p;
156  myDist = sqrt((s-la)*(s-lb)*(s-lc)/s);
157  myDirs[0] = UT_Vector2(sidea.y(), -sidea.x());
158  myDirs[1] = UT_Vector2(sideb.y(), -sideb.x());
159  myDirs[2] = UT_Vector2(sidec.y(), -sidec.x());
160  myDirs[0].normalize();
161  myDirs[1].normalize();
162  myDirs[2].normalize();
163  // Winding order matters for myDirs.
164  // We need to make sure that they point toward their corresponding
165  // sides, instead of away from them, else we'll get the triangle
166  // rotated a half turn. Either point on side a works for checking
167  // dir 0, and either all of the dirs are inverted, or none of them are.
168  if (dot(myDirs[0], b - myIncentre) < 0)
169  {
170  myDirs[0] = -myDirs[0];
171  myDirs[1] = -myDirs[1];
172  myDirs[2] = -myDirs[2];
173  }
174  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[0], c - myIncentre), 0));
175  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[1], c - myIncentre), 0));
176  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[1], a - myIncentre), 0));
177  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[2], a - myIncentre), 0));
178  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[2], b - myIncentre), 0));
179  }
180 
181  /// This can be an underestimate, but not an overestimate
182  /// of the distance squared.
183  float boxDist(const UT_BoundingRect &box) const
184  {
185  // Expand the box to the minimum circle fully containing it
186  // and then compute the maximum of the 3 distances.
187  float boxradius = 0.5f*SYSsqrt(box.sizeX()*box.sizeX()
188  + box.sizeY()*box.sizeY());
189  UT_Vector2 boxcentre(box.centerX(), box.centerY());
190  UT_Vector2 dir = boxcentre - myIncentre;
191  float d0 = dot(dir, myDirs[0]);
192  float d1 = dot(dir, myDirs[1]);
193  float d2 = dot(dir, myDirs[2]);
194  float d = SYSmax(d0, d1, d2) - myDist - boxradius;
195  float dsquared = d*d;
196  return (d < 0) ? -dsquared : dsquared;
197  }
198  /// This distance squared needs to be exact.
199  float dist(const float *P, int) const
200  {
201  // Compute the maximum of the 3 distances.
202  UT_Vector2 dir = UT_Vector2(P[0], P[1]) - myIncentre;
203  float d0 = dot(dir, myDirs[0]);
204  float d1 = dot(dir, myDirs[1]);
205  float d2 = dot(dir, myDirs[2]);
206  float d = SYSmax(d0, d1, d2) - myDist;
207  float dsquared = d*d;
208  return (d < 0) ? -dsquared : dsquared;
209  }
210 
211 private:
212  UT_Vector2 myDirs[3];
213  UT_Vector2 myIncentre;
214  float myDist;
215 };
216 
217 /// Lookup information for 3D-tree tetrahedron queries
218 /// NOTE: Distance squared here is not Euclidean distance squared, but
219 /// face-perpendicular distance squared, i.e. distance is from
220 /// the incentre out, perpendicular to one of the faces, minus the
221 /// insphere's radius. This avoids the need to have special
222 /// cases for what would be the spherical sections around each vertex.
223 /// It basically indicates how far the tetrahedron would have to be
224 /// expanded (or contracted) relative to the incentre in order to
225 /// hit the query point.
227 public:
229  const UT_Vector3 &a,
230  const UT_Vector3 &b,
231  const UT_Vector3 &c,
232  const UT_Vector3 &d)
233  {
234  UT_Vector3 edgeda = a-d;
235  UT_Vector3 edgedb = b-d;
236  UT_Vector3 edgedc = c-d;
237  UT_Vector3 edgeab = b-a;
238  UT_Vector3 edgebc = c-b;
239  myDirs[0] = cross(edgedb, edgedc);
240  float volume = SYSabs(dot(edgeda, myDirs[0]))/6;
241  myDirs[1] = cross(edgedc, edgeda);
242  myDirs[2] = cross(edgeda, edgedb);
243  myDirs[3] = cross(edgebc, edgeab);
244  float areaa = myDirs[0].normalize()/2;
245  float areab = myDirs[1].normalize()/2;
246  float areac = myDirs[2].normalize()/2;
247  float aread = myDirs[3].normalize()/2;
248  float area = areaa + areab + areac + aread;
249  myDist = (3*volume)/area;
250  myIncentre = (areaa*a + areab*b + areac*c + aread*d)/area;
251 
252  // Winding order matters for myDirs.
253  // We need to make sure that they point toward their corresponding
254  // sides, instead of away from them, else we'll get the tetrahedron
255  // inverted through the incentre. Any point on side a works for checking
256  // dir 0, and either all of the dirs are inverted, or none of them are.
257  if (dot(myDirs[0], b - myIncentre) < 0)
258  {
259  myDirs[0] = -myDirs[0];
260  myDirs[1] = -myDirs[1];
261  myDirs[2] = -myDirs[2];
262  myDirs[3] = -myDirs[3];
263  }
264  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[0], c - myIncentre), 0));
265  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[0], d - myIncentre), 0));
266  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[1], c - myIncentre), 0));
267  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[1], d - myIncentre), 0));
268  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[1], a - myIncentre), 0));
269  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[2], d - myIncentre), 0));
270  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[2], a - myIncentre), 0));
271  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[2], b - myIncentre), 0));
272  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[3], a - myIncentre), 0));
273  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[3], b - myIncentre), 0));
274  UT_ASSERT_P(SYSisGreaterOrEqual(dot(myDirs[3], c - myIncentre), 0));
275  }
276 
277  /// This can be an underestimate, but not an overestimate
278  /// of the distance squared.
279  float boxDist(const UT_BoundingBox &box) const
280  {
281  // Expand the box to the minimum sphere fully containing it
282  // and then compute the maximum of the 4 distances.
283  float boxradius = box.getRadius();
284  UT_Vector3 boxcentre = box.center();
285  UT_Vector3 dir = boxcentre - myIncentre;
286  float d0 = dot(dir, myDirs[0]);
287  float d1 = dot(dir, myDirs[1]);
288  float d2 = dot(dir, myDirs[2]);
289  float d3 = dot(dir, myDirs[3]);
290  float d = SYSmax(d0, d1, d2, d3) - myDist - boxradius;
291  float dsquared = d*d;
292  return (d < 0) ? -dsquared : dsquared;
293  }
294  /// This distance squared needs to be exact.
295  float dist(const float *P, int) const
296  {
297  // Compute the maximum of the 4 distances.
298  UT_Vector3 dir = UT_Vector3(P[0], P[1], P[2]) - myIncentre;
299  float d0 = dot(dir, myDirs[0]);
300  float d1 = dot(dir, myDirs[1]);
301  float d2 = dot(dir, myDirs[2]);
302  float d3 = dot(dir, myDirs[3]);
303  float d = SYSmax(d0, d1, d2, d3) - myDist;
304  float dsquared = d*d;
305  return (d < 0) ? -dsquared : dsquared;
306  }
307 
308 private:
309  UT_Vector3 myDirs[4];
310  UT_Vector3 myIncentre;
311  float myDist;
312 };
313 
314 /// Lookup information for line queries
316 public:
318  const UT_Vector3 &dir)
319  {
320  myP0 = orig;
321  myP1 = orig+dir;
322  myDir = dir;
323  }
324 
325  /// This can be an underestimate, but not an overestimate
326  /// of the distance squared.
327  float boxDist(const UT_BoundingBox &box) const
328  {
329  return box.approxLineDist2(myP0, myDir);
330  }
331  /// This distance squared needs to be exact.
332  float dist(const float *P, int) const
333  {
334  return segmentPointDist2(UT_Vector3(P), myP0, myP1);
335  }
336 
337 private:
338  UT_Vector3 myP0;
339  UT_Vector3 myP1;
340  UT_Vector3 myDir;
341 };
342 
343 /// Queries for infinite lines (infinite tubes)
345 public:
347  const UT_Vector3 &dir)
348  {
349  myP0 = orig;
350  myDir = dir;
351  }
352 
353  /// This can be an underestimate, but not an overestimate
354  /// of the distance squared.
355  float boxDist(const UT_BoundingBox &box) const
356  {
357  UT_Vector3 pos = box.center();
358  pos -= myP0;
359  pos -= myDir.project(pos);
360  float dist = SYSmax(pos.length() - box.getRadius(), 0.0);
361  return dist*dist;
362  }
363  /// This distance squared needs to be exact.
364  float dist(const float *P, int) const
365  {
366  UT_Vector3 pos(P);
367  pos -= myP0;
368  pos -= myDir.project(pos);
369  return pos.length2();
370  }
371 
372 private:
373  UT_Vector3 myP0;
374  UT_Vector3 myDir;
375 };
376 
377 #if 0
378 /// Queries for spherical cones
379 class UT_KDConeQuery {
380 public:
381  UT_KDConeQuery(
382  const UT_Vector3 &orig,
383  const UT_Vector3 &dir,
384  const float dotthreshold,
385  const float maxdistancesquared = 1e37)
386  {
387  myP0 = orig;
388  myDir = dir;
389  myDotThreshold = dotthreshold;
390  mySine = SYSsqrt(1-myDotThreshold*myDotThreshold);
391  myMaxDistaceSquared = maxdistancesquared;
392  myMaxDistace = SYSsqrt(myMaxDistaceSquared);
393  }
394 
395  /// This can be an underestimate, but not an overestimate
396  /// of the distance squared.
397  float boxDist(const UT_BoundingBox &box) const
398  {
399  if (box.isInside(myP0))
400  return 0;
401  UT_ASSERT_MSG(0, "Unimplemented!!! Note: if box is completely outside cone, should return something strictly > myMaxDistanceSquared");
402  return 2*myMaxDistaceSquared;
403  }
404  /// This distance squared needs to be exact.
405  float dist(const float *P, int) const
406  {
407  UT_Vector3 pos(P);
408  pos -= myP0;
409 
410  // Counterintuitively, this just uses regular distance squared.
411  // Points outside the cone are rejected by isValid().
412  float actualdist2 = pos.length2();
413 
414  return actualdist2;
415  }
416  float distOutsideCone(const float *P) const
417  {
418  UT_Vector3 pos(P);
419  pos -= myP0;
420  UT_Vector3 dir = pos;
421  dir.normalize();
422 
423  float actualdist2 = pos.length2();
424 
425  // If inside cone, return -1
426  float dt = dir.dot(myDir);
427  if (dt >= myDotThreshold)
428  return -1;
429 
430  // If outside, use max of actualdist2 and (maxdist + dist to cone side)^2
431 
432  UT_Vector3 planex = pos - myDir.project(pos);
433  UT_Vector3 sideraydir = myDotThreshold*myDir + mySine*planex;
434  float raydist = sideraydir.dot(pos);
435  float disttoconeside;
436  if (raydist <= 0)
437  disttoconeside = SYSsqrt(actualdist2);
438  else
439  disttoconeside = (pos - raydist*sideraydir).length();
440 
441  return disttoconeside;
442  }
443 
444 private:
445  UT_Vector3 myP0;
446  UT_Vector3 myDir;
447  float myDotThreshold;
448  float mySine;
449  float myMaxDistace;
450  float myMaxDistaceSquared;
451 };
452 #endif
453 
455 {
456 public:
457  /// KD Tree balancing algorithms. See setBalancer.
458  typedef enum {
461  UT_KD_CENTROID
462  } ut_KDBalancer;
463 
464 public:
465  UT_KDTree(int dim=3, size_t size=0)
466  {
467  myDim = dim;
468  myList = 0;
469  myMaxLeafNodes = 25;
470  myRebalanceCount = 128;
471  myBalancer = UT_KD_MEDIAN;
472  myHasRadius = false;
473  setEntries(size);
474  }
475  virtual ~UT_KDTree();
476 
477  /// This must be called before querying, to build and balance the tree.
478  virtual void buildIfNeeded(bool enable_multithreading = true)
479  {
480  balance(enable_multithreading);
481  }
482 
483  int64 getMemoryUsage(bool inclusive) const
484  {
485  int64 mem = inclusive ? sizeof(*this) : 0;
486  mem += mySplits.getMemoryUsage(false);
487  mem += myLock.getMemoryUsage(false);
488  if (myList)
489  mem += sizeof(myList[0])*myEntries;
490  return mem;
491  }
492 
493  /// setEntries instructs the KD Tree about the number of points
494  /// that it should read back from its callbacks.
495  /// Note that each call to setEntries has O(size) cost! Thus,
496  /// for (i = 0; i < n; i++) tree.setEntries(i);
497  /// is O(n^2). There isn't a good reason for this (as one could instead
498  /// just mark the tree as dirty until the next query, much like
499  /// already happens with rebalancing) but it is the way things are so be
500  // /careful.
501  void setEntries(size_t size);
502  size_t getEntries() const { return myFullEntries; }
503 
504  /// Sometimes, a tree might have new points added since it's construction.
505  /// Rather than re-balancing the tree (which can be relatively expensive),
506  /// we have allow N points to be added to the tree without re-balancing.
507  /// To add entries to the tree, simply call "growEntries" with the number
508  /// of points being added. If the size of the unsorted pool exceeds the
509  /// re-balance count, then the tree will be re-balanced.
510  ///
511  /// The rebalance count is automatically scaled as the tree is constructed,
512  /// so it shouldn't be required to set the count (unless you know what
513  /// you're doing).
514  ///
515  /// To force a re-build of the tree, simply call setEntries(), or flag the
516  /// tree as dirty.
517  void growEntries(size_t amount);
518  size_t getRebalanceCount() const
519  {
520  return myRebalanceCount;
521  }
522  void setRebalanceCount(size_t count);
523 
524  /// Sets the KD-tree balancing algorithm.
525  /// - UT_KD_MEDIAN: splits at the median point
526  /// - UT_KD_CENTROID: splits at the spatial centroid
527  /// - UT_KD_SAH: splits at SAH optimal splitting point
528  /// The UT_KD_SAH (surface area heuristic) algorithm builds more
529  /// efficient trees but they are slower to construct. The default is
530  /// UT_KD_MEDIAN.
532  {
533  myBalancer = balance;
534  }
535  /// Sets the maximum number of nodes to be stored in a leaf. Smaller
536  /// values will produce deeper and more memory-hungry trees.
537  void setMaxLeafNodes(int max_leaf_nodes)
538  {
539  UT_ASSERT_P(max_leaf_nodes >= 3);
540  myMaxLeafNodes = SYSmax(3, max_leaf_nodes);
541  }
542 
543  /// Marks the tree as dirty. Note, however, that this has
544  /// O(size) time cost.
545  void flagDirty()
546  {
547  // Re-build the lists
548  setEntries(myFullEntries);
549  }
550 
551  /// The compare function should compare the distance between two points
552  /// (given by idx0 and idx1) in the dimension specified. The dimension will
553  /// be between 0 and 3. For example,
554  /// comparePosition(...) {
555  /// fpreal delta = myP[idx0](dim) - myP[idx1](dim);
556  /// return delta < 0 ? -1 : (delta > 0 ? 1 : 0);
557  /// }
558  virtual int comparePosition(int idx0,
559  int idx1, int dim) const = 0;
560  /// Return the position associated with the given point
561  virtual const float *getP(int idx) const = 0;
562 
563  /// If each point in the KD Tree has a radius associated with it, then this
564  /// method should return true. Adding a per-point radius can affect
565  /// performance adversely.
566  virtual bool pointsHaveRadius() const { return false; }
567  /// Return the radius associated with the point in question
568  virtual fpreal getRadius(int /*idx*/) const { return 0; }
569 
570  /// Returns whether the given index should be considered in searches
571  virtual bool isValid(int /*idx*/) const
572  { return true; }
573  virtual bool isValid(int idx, const UT_KDQueryPt & /*pt*/) const
574  { return isValid(idx); }
575 
576  /// Return the maximum number of invalid indices that should be
577  /// considered before bailing. Return -1 for no limit.
578  /// maxn - the search limit
579  virtual int getInvalidLimit(int maxn) const { return -1; }
580 
581  /// Find the closest node to the position specified. This method
582  /// returns -1 if no photon is found.
583  /// NOTE THAT max_distance_squared IS SQUARED DISTANCE.
584  template <typename QueryPoint>
585  int findClosest(const QueryPoint &pt,
586  fpreal max_distance_squared);
587  template <typename QueryPoint>
588  int findClosestQueue(const QueryPoint &pt,
589  ut_KDPQueue &queue,
590  fpreal max_distance_squared);
591 
592  /// Find the closest N (max_nodes) nodes to the position given. Only points
593  /// within the sphere defined by max_distance_squared will be considered.
594  /// NOTE THAT max_distance_squared IS SQUARED DISTANCE.
595  template <typename QueryPoint>
596  int findClosest(UT_IntArray &list,
597  const QueryPoint &pt,
598  fpreal max_distance_squared,
599  int max_nodes,
600  bool sorted = true);
601  template <typename QueryPoint>
602  int findClosestQueue(UT_IntArray &list,
603  const QueryPoint &pt,
604  ut_KDPQueue &q,
605  fpreal max_distance_squared,
606  int max_nodes,
607  bool sorted = true);
608 
609  /// Find the closest N (max_nodes) nodes to the position given. Only points
610  /// within the sphere defined by max_distance_squared will be considered.
611  /// NOTE THAT max_distance_squared IS SQUARED DISTANCE.
612  template <typename QueryPoint>
613  int findClosest(UT_IntArray &list,
615  const QueryPoint &pt,
616  fpreal max_distance_squared,
617  int max_nodes,
618  bool sorted = true);
619  template <typename QueryPoint>
620  int findClosestQueue(UT_IntArray &list,
622  const QueryPoint &pt,
623  ut_KDPQueue &q,
624  fpreal max_distance_squared,
625  int max_nodes,
626  bool sorted = true);
627  template <typename QueryPoint>
629  const QueryPoint &pt,
630  int max_nodes)
631  {
632  return findClosest(list, pt, 1e37, max_nodes);
633  }
634  template <typename QueryPoint>
636  const QueryPoint &pt,
637  fpreal max_distance_squared)
638  {
639  return findClosest(list, pt, max_distance_squared,
640  myFullEntries, false);
641  }
642 
643  /// Interface for aggregate data used by volume filtering
644  class VolumeData {
645  public:
646  virtual ~VolumeData() {}
647  virtual void sum(const VolumeData &data, float weight) = 0;
648  virtual void scale(float scale) = 0;
649  };
650 
651  /// Interface for tree traversals
652  class Visitor {
653  public:
654  virtual ~Visitor() {}
655 
656  /// Before traversal, this method is called to indicate the total
657  /// number of nodes (internal + leaf) that will be traversed.
658  /// nodeid passed to the visit functions will take on values in the
659  /// range [0, num_nodes-1].
660  virtual void setNumNodes(int num_nodes)
661  {
662  }
663 
664  /// Traverses a leaf node in the tree. Leaf nodes contain a list
665  /// of points (specified by point_list and size) that is disjoint
666  /// from other leaf nodes.
667  virtual void visitLeaf(int nodeid,
668  const int *point_list, int size,
669  const UT_BoundingBox &box) = 0;
670 
671  /// Post-traverses an internal node (left and right are guaranteed
672  /// to have been traversed before this method is called)
673  virtual void postVisitInternal(int nodeid, int leftid, int rightid,
674  const UT_BoundingBox &box,
675  int axis, float split)
676  {
677  }
678  };
679 
680  /// Traverse the KD Tree in depth first order. This routine only works
681  /// on 3D trees.
682  /// NOTE: This will call balance() if the tree isn't balanced.
683  void traverse(Visitor &visitor);
684 
685  /// This class knows how to compute aggregate volume data (using a tree
686  /// traversal) and also provides an accessor to the computed data.
687  class AggregateVolume : public Visitor {
688  public:
689  virtual const VolumeData &getData(int nodeid) const = 0;
690  };
691 
692  ///
693  /// Filtered evaluation for 3D volumes. This method interprets the
694  /// KD-Tree as a volume, and filters aggregate point data using the
695  /// volumetric coverage for aggregate nodes in the tree. The max_nodes
696  /// parameter controls the approximate number of points that should be
697  /// averaged to produce the filtered result. The actual filter radius
698  /// used will depend on max_nodes and on the size of KD-Tree nodes
699  /// close to the query point. Since this is an approximate query, the
700  /// actual number of points that contribute to the result may be larger
701  /// than max_nodes.
702  ///
703  /// This method relies on implemented subclasses for AggregateVolume
704  /// and for VolumeData. You must compute aggregate volume data before
705  /// calling this method.
706  ///
707  /// UT_KDTree tree;
708  /// AggregateVolumeSubclass aggdata;
709  /// VolumeDataSubclass data;
710  ///
711  /// // Store aggregate values (do this once)
712  /// kdtree.traverse(aggdata);
713  ///
714  /// // Filter aggregate data (do this many times)
715  /// kdtree.filterVolume(data, pos, filter, aggdata, max_nodes);
716  ///
717  void filterVolume(VolumeData &data,
718  const UT_Vector3 &pos,
719  const UT_Filter &filter,
720  const AggregateVolume &aggdata,
721  int max_nodes);
722 
723  const fpreal *getBoxMin();
724  const fpreal *getBoxMax();
725 
726  /// This allows you to create a search queue so you can maintain
727  /// it over time, avoiding the cost of reallocing for each search.
728  /// NOTE THAT max_distance_squared IS SQUARED DISTANCE.
729  static ut_KDPQueue *createQueue();
730  static void destroyQueue(ut_KDPQueue *q);
731 
732 protected:
733  //
734  // KD tree splits
735  //
736 
737  class KDSplit {
738  public:
740  {
741  myLeft = myRight = -1;
742  mySplitIdx = 0;
743  myAxis = 0;
744  }
746  {
747  }
748  void init(int left, int right,
749  fpreal off, fpreal radius, int mid, int axis)
750  {
751  myLeft = left;
752  myRight = right;
753  myOffset = off;
754  myRadius = radius;
755  mySplitIdx = mid;
756  myAxis = axis;
757  }
758 
759  int left() const { return myLeft; }
760  int right() const { return myRight; }
761  fpreal offset() const { return myOffset; }
762  fpreal radius() const { return myRadius; }
763  int splitIdx() const { return mySplitIdx; }
764  int axis() const { return myAxis; }
765 
766  private:
767  int myLeft;
768  int myRight;
769  fpreal myOffset;
770 
771  /// This stores the maximum radius of any children of this node. The
772  /// idea is that this is the maximum distance you can be on the wrong
773  /// side of a splitting plane and still be legal.
774  fpreal myRadius;
775  int mySplitIdx;
776  int myAxis;
777  };
778 
779  int getHead() const { return mySplits.entries() ? 0 : -1; }
780 
781  // Dummy definition required for tube, triangle, and tetrahedron searches
782  bool isValid(int node, const UT_KDTubeQuery &) const
783  { return true; }
784  bool isValid(int node, const UT_KDLineQuery &) const
785  { return true; }
786  bool isValid(int node, const UT_KDTriQuery &) const
787  { return true; }
788  bool isValid(int node, const UT_KDTetQuery &) const
789  { return true; }
790  bool isValid(int node, const UT_KDQueryPtUnitWrap &) const
791  { return true; }
792 #if 0
793  bool isValid(int node, const UT_KDConeQuery &q) const
794  {
795  if (!isValid(node))
796  return false;
797 
798  float dist = q.distOutsideCone(getP(node));
799  if (dist <= 0)
800  return true;
801  if (!myHasRadius)
802  return false;
803  return getRadius(node) > dist;
804  }
805 #endif
806 
807  template <typename QueryPoint>
808  int findClosest(ut_KDPQueue &list, const QueryPoint &pt);
809  template <typename QueryPoint>
810  void recurseFind(ut_KDPQueue &list, const QueryPoint &pt,
811  int split, int lo, int hi) const;
812  template <typename QueryPoint>
813  void recurseFindTube(
814  ut_KDPQueue &list, const QueryPoint &pt,
815  int split, int lo, int hi) const;
816  template <typename QueryPoint>
817  void recurseFindTri(
818  ut_KDPQueue &list, const QueryPoint &pt,
819  int split, int lo, int hi) const;
820  template <typename QueryPoint>
821  void findInLeaf(ut_KDPQueue &list, const QueryPoint &pt,
822  int lo, int hi, int invalid_limit,
823  int &invalid) const;
824 
825  bool isBalanced() const { return myBalanced; }
826 
827  void traverseRecursive(Visitor &visitor,
828  int split, int nodeid,
829  UT_BoundingBox &box,
830  int lo, int hi);
831 
832  void computeBox(int start_index=0);
833  bool isBoxClose(const fpreal *P, fpreal qd, fpreal r) const;
834 
835  /// Creates a KD tree from an unsorted set of points.
836  /// Given a list of points, we:
837  /// 1) Find the dimension in which there is most change.
838  /// 2) Use the balancing algorithm to choose a splitting plane and
839  /// partition myList into left and right portions.
840  /// 3) Store the splitting information in mySplits
841  /// 4) Recurse on the two halves.
842  /// This is repeated until myMaxLeafNodes is reached where we just
843  /// leave it as a linear list.
844  void balance(bool enable_multithreading=true);
845  /// Balances a subset, returning the maximum radius of that subset.
846  /// If splitlock is NULL, the balancing will all be done serially
847  /// instead of in parallel.
848  void balanceSet(int &split, fpreal &radius, int *list, int entries,
849  fpreal* boxmin, fpreal* boxmax, UT_Lock* splitlock);
850 
851  /// Stores the bounding box of all of the points in the KD tree, including
852  /// their individual radii, in each of the possible dimensions.
853  /// Note that during balanceSet these values are temporarily altered
854  /// before recursion, so don't be too shocked.
857 
858  /// Nodes in a KD tree have two indices. One is the index that
859  /// they were added in. This is the index used in getP, for example.
860  /// The other is there location in the binary heap. myList takes
861  /// a binary heap location and returns the index useable for
862  /// getP().
863  int *myList;
864 
865  /// Split information stored as a tree structure. This stores the
866  /// midpoints, splitting offset, axis, and max radius.
868 
869  /// For protecting tree balancing and bounding box computation
871 
872  /// This marks the number of entries that have been added to our tree
873  size_t myEntries;
874 
875  /// This marks the total number of entries in this data structure.
876  /// The difference between this and myEntries is the unsorted
877  /// chaff at the end.
880  int myDim;
884 
885  /// A faster way to evaluate pointsHaveRadius() at runtime.
888 };
889 
890 #endif
#define SYSmax(a, b)
Definition: SYS_Math.h:1365
void setBalancer(ut_KDBalancer balance)
Definition: UT_KDTree.h:531
int myMaxLeafNodes
Definition: UT_KDTree.h:881
int findNClosest(UT_IntArray &list, const QueryPoint &pt, int max_nodes)
Definition: UT_KDTree.h:628
GA_API const UT_StringHolder dist
virtual void setNumNodes(int num_nodes)
Definition: UT_KDTree.h:660
bool isValid(int node, const UT_KDLineQuery &) const
Definition: UT_KDTree.h:784
bool isValid(int node, const UT_KDTriQuery &) const
Definition: UT_KDTree.h:786
bool isValid(int node, const UT_KDTubeQuery &) const
Definition: UT_KDTree.h:782
UT_Vector3T< T > project(const UT_Vector3T< T > &u) const
Calculates the orthogonal projection of a vector u on the *this vector.
Queries for infinite lines (infinite tubes)
Definition: UT_KDTree.h:344
Lookup information for line queries.
Definition: UT_KDTree.h:315
int getHead() const
Definition: UT_KDTree.h:779
UT_Vector2T< float > UT_Vector2
GLint left
Definition: glcorearb.h:2004
float boxDist(const UT_BoundingBox &box) const
Definition: UT_KDTree.h:327
virtual bool isValid(int) const
Returns whether the given index should be considered in searches.
Definition: UT_KDTree.h:571
void setMaxLeafNodes(int max_leaf_nodes)
Definition: UT_KDTree.h:537
UT_KDTriQuery(const UT_Vector2 &a, const UT_Vector2 &b, const UT_Vector2 &c)
Definition: UT_KDTree.h:142
bool myBoxComputed
Definition: UT_KDTree.h:883
virtual ~VolumeData()
Definition: UT_KDTree.h:646
bool isBalanced() const
Definition: UT_KDTree.h:825
UT_Vector3T< float > UT_Vector3
size_t getEntries() const
Definition: UT_KDTree.h:502
T approxLineDist2(const UT_Vector3T< T > &v0, const UT_Vector3T< T > &dir) const
size_t getRebalanceCount() const
Definition: UT_KDTree.h:518
UT_KDQueryPtUnitWrap(const float *p, int ndims=3)
Definition: UT_KDTree.h:83
UT_Array< KDSplit > mySplits
Definition: UT_KDTree.h:867
T & x(void)
Definition: UT_Vector2.h:284
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
#define SYSabs(a)
Definition: SYS_Math.h:1367
float dist(int axis, float off) const
Definition: UT_KDTree.h:53
void flagDirty()
Definition: UT_KDTree.h:545
#define UT_API
Definition: UT_API.h:12
int left() const
Definition: UT_KDTree.h:759
int64 getMemoryUsage(bool inclusive) const
Definition: UT_KDTree.h:483
T sizeAxis(int axis) const
png_uint_32 i
Definition: png.h:2877
float dist(const float *P, int) const
This distance squared needs to be exact.
Definition: UT_KDTree.h:109
GLsizeiptr size
Definition: glcorearb.h:663
int findAllClosest(UT_IntArray &list, const QueryPoint &pt, fpreal max_distance_squared)
Definition: UT_KDTree.h:635
UT_Lock myLock
For protecting tree balancing and bounding box computation.
Definition: UT_KDTree.h:870
Interface for tree traversals.
Definition: UT_KDTree.h:652
T segmentPointDist2(const UT_Vector3T< T > &pos, const UT_Vector3T< T > &pt1, const UT_Vector3T< T > &pt2)
Definition: UT_Vector3.h:988
#define UT_ASSERT_P(ZZ)
Definition: UT_Assert.h:101
SYS_FORCE_INLINE UT_StorageNum< typename UT_StorageBetter< T, S >::Type >::AtLeast32Bit dot(const UT_FixedVector< S, SIZE, S_INSTANTIATED > &that) const
virtual ~Visitor()
Definition: UT_KDTree.h:654
long long int64
Definition: SYS_Types.h:100
virtual void postVisitInternal(int nodeid, int leftid, int rightid, const UT_BoundingBox &box, int axis, float split)
Definition: UT_KDTree.h:673
GA_API const UT_StringHolder scale
UT_KDQueryPt(const float *p)
Definition: UT_KDTree.h:44
UT_Vector3T< T > center() const
Lookup point information to be passed to the query functions.
Definition: UT_KDTree.h:43
UT_KDQueryPt(const float *p, const void *data)
Definition: UT_KDTree.h:46
float boxDist(const UT_BoundingBox &box) const
Definition: UT_KDTree.h:279
fpreal64 dot(const CE_VectorT< T > &a, const CE_VectorT< T > &b)
Definition: CE_Vector.h:218
float boxDist(int axis, float bmin, float bmax) const
Definition: UT_KDTree.h:49
virtual void buildIfNeeded(bool enable_multithreading=true)
This must be called before querying, to build and balance the tree.
Definition: UT_KDTree.h:478
Interface for aggregate data used by volume filtering.
Definition: UT_KDTree.h:644
UT_KDTetQuery(const UT_Vector3 &a, const UT_Vector3 &b, const UT_Vector3 &c, const UT_Vector3 &d)
Definition: UT_KDTree.h:228
float boxDist(const UT_BoundingBox &box) const
Definition: UT_KDTree.h:90
int right() const
Definition: UT_KDTree.h:760
int splitIdx() const
Definition: UT_KDTree.h:763
const float * myP
Definition: UT_KDTree.h:127
GLboolean * data
Definition: glcorearb.h:130
virtual bool pointsHaveRadius() const
Definition: UT_KDTree.h:566
virtual bool isValid(int idx, const UT_KDQueryPt &) const
Definition: UT_KDTree.h:573
UT_KDTree(int dim=3, size_t size=0)
Definition: UT_KDTree.h:465
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
fpreal offset() const
Definition: UT_KDTree.h:761
float dist(const float *P, int) const
This distance squared needs to be exact.
Definition: UT_KDTree.h:295
GLint GLsizei count
Definition: glcorearb.h:404
ut_KDBalancer myBalancer
Definition: UT_KDTree.h:887
UT_KDLineQuery(const UT_Vector3 &orig, const UT_Vector3 &dir)
Definition: UT_KDTree.h:346
double fpreal
Definition: SYS_Types.h:263
int myDim
Definition: UT_KDTree.h:880
bool myBalanced
Definition: UT_KDTree.h:882
png_infop png_sPLT_tpp entries
Definition: png.h:2481
int axis() const
Definition: UT_KDTree.h:764
T getRadius() const
Returns the radius of a sphere that would fully enclose the box.
float dist(const float *P, int) const
This distance squared needs to be exact.
Definition: UT_KDTree.h:199
float dist(const float *P, int) const
This distance squared needs to be exact.
Definition: UT_KDTree.h:364
#define UT_KD_MAXDIM
Definition: UT_KDTree.h:40
size_t myRebalanceCount
Definition: UT_KDTree.h:879
const float * P
Definition: UT_KDTree.h:75
ut_KDBalancer
KD Tree balancing algorithms. See setBalancer.
Definition: UT_KDTree.h:458
virtual int getInvalidLimit(int maxn) const
Definition: UT_KDTree.h:579
bool myHasRadius
A faster way to evaluate pointsHaveRadius() at runtime.
Definition: UT_KDTree.h:886
const void * myData
Definition: UT_KDTree.h:76
fpreal radius() const
Definition: UT_KDTree.h:762
virtual fpreal getRadius(int) const
Return the radius associated with the point in question.
Definition: UT_KDTree.h:568
SYS_FORCE_INLINE Storage::MathFloat normalize()
int isInside(const UT_Vector3T< T > &pt) const
size_t myEntries
This marks the number of entries that have been added to our tree.
Definition: UT_KDTree.h:873
float boxDist(const UT_BoundingBox &box) const
Definition: UT_KDTree.h:355
float boxDist(const UT_BoundingRect &box) const
Definition: UT_KDTree.h:183
T & y(void)
Definition: UT_Vector2.h:286
SYS_FORCE_INLINE Storage::AtLeast32Bit length2() const noexcept
GLboolean r
Definition: glcorearb.h:1221
void init(int left, int right, fpreal off, fpreal radius, int mid, int axis)
Definition: UT_KDTree.h:748
float dist(const float *P, int) const
This distance squared needs to be exact.
Definition: UT_KDTree.h:332
bool isValid(int node, const UT_KDQueryPtUnitWrap &) const
Definition: UT_KDTree.h:790
SYS_FORCE_INLINE Storage::MathFloat length() const
#define SYSmin(a, b)
Definition: SYS_Math.h:1366
size_t myFullEntries
Definition: UT_KDTree.h:878
#define UT_ASSERT_MSG(ZZ, MM)
Definition: UT_Assert.h:105
int * myList
Definition: UT_KDTree.h:863
T centerAxis(int axis) const
SIM_DerVector3 cross(const SIM_DerVector3 &lhs, const SIM_DerVector3 &rhs)
bool isValid(int node, const UT_KDTetQuery &) const
Definition: UT_KDTree.h:788
UT_KDTubeQuery(const UT_Vector3 &orig, const UT_Vector3 &dir)
Definition: UT_KDTree.h:317
GA_API const UT_StringHolder area
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:794
GLint GLint GLint GLint GLint GLint GLint GLbitfield GLenum filter
Definition: glcorearb.h:1296
float dist(const float *P, int dim) const
Definition: UT_KDTree.h:57