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