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