HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UT_SpatialTree.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_SpatialTree.h (UT Library, C++)
7  *
8  * COMMENTS:
9  * An n-dimensional tree implementation (e.g. quadtree, octree, etc.)
10  * that stores arbitrary data and supports arbitrary bounding regions
11  * and subdivision strategies.
12  *
13  * A bounding region provider must implement getDistance2(), getMin(),
14  * and getMax() as follows:
15  *
16  * /// Returns the distance (squared) from the bounding region
17  * /// corresponding to the supplied object to the specified point.
18  * fpreal getDistance2(T object, const fpreal *pt) const;
19  *
20  * /// Returns the maximum coordinates of the bounding box for the
21  * /// bounding region corresponding to the supplied object. The
22  * /// supplied buffer may be used to return the result if neccessary.
23  * /// Its size is equal to sizeof(fpreal) * ndimensions, where
24  * /// ndimensions is the dimension of the tree.
25  * const fpreal *getMax(T object, fpreal *buf) const;
26  *
27  * /// Returns the mnimum coordinates of the bounding box for the
28  * /// bounding region corresponding to the supplied object. The
29  * /// supplied buffer may be used to return the result if neccessary.
30  * /// Its size is equal to sizeof(fpreal) * ndimensions, where
31  * /// ndimensions is the dimension of the tree.
32  * const fpreal *getMin(T object, fpreal *buf) const;
33  *
34  * A subdivision strategy must implement subdivide() as follows:
35  *
36  * /// Returns true if the node at the specified depth with the
37  * /// specified bounds (min/max), containing the specified
38  * /// objects, whose bounds are given, should be subdivided.
39  * template<class T, class B>
40  * bool subdivide(
41  * const B *bounds_provider,
42  * const fpreal *min,
43  * const fpreal *max,
44  * const UT_ValArray<T> &objects,
45  * int depth,
46  * int ndimensions);
47  */
48 
49 #ifndef __UT_SpatialTree__
50 #define __UT_SpatialTree__
51 
52 #include "UT_API.h"
53 
54 #include "UT_RingBuffer.h"
55 #include "UT_SmallObject.h"
56 #include "UT_StackAlloc.h"
57 #include "UT_ValArray.h"
58 
59 #include <iostream>
60 #include <string.h>
61 #include <stack>
62 
63 #define MIN_NODE_SIZE 1e-3
64 
66 {
67 public:
68  UT_STAccepter(bool accept) : myAccept(accept) {}
69  bool accept(int objectno) const { return myAccept; }
70 
71 private:
72  bool myAccept;
73 };
74 
75 template<class T, class S, class B, int N>
77 {
78  enum { BranchFactor = 1 << N };
79 public:
80  /// Construct a tree with the specified subdivision strategy
81  /// (subd_strategy must not be null), the bounds provider (must not
82  /// be null), the specified bounds (min/max), and with the specified
83  /// dimensions.
85  S *subd_strategy,
86  B *bounds_provider,
87  const fpreal *min,
88  const fpreal *max);
89 
90  /// Construct a tree with the specified subdivision strategy
91  /// (subd_strategy must not be null), the specified bounds provider
92  /// (must not be null), and with the specified dimensions.
94  S *subd_strategy,
95  B *bounds_provider);
96 
97  /// Destructor.
99 
100  /// Add the specified object to the tree.
101  void addObject(T object);
102 
103  /// Returns the bounds provider used by this tree.
104  const B *getBoundsProvider() const { return myBoundsProvider; }
105 
106  /// Return the maximum number of dimensions allowed.
107  static int getMaxDimensions() { return 31; }
108 
109  /// Return the objects in the tree whose bounds contain the specified
110  /// point. Objects are returned in the supplied array.
111  template<class A>
112  void getObjects(
113  const fpreal *pt,
114  UT_ValArray<T> &objects,
115  const A &accepter) const;
116 
118  const fpreal *pt,
119  UT_ValArray<T> &objects) const
120  {
121  getObjects(pt, objects, UT_STAccepter(true));
122  }
123 
124  /// Return the objects in the tree whose bounds are withing the
125  /// specified radius of the specified point.
126  template<class A>
127  int getObjects(
128  const fpreal *pt,
129  UT_ValArray<T> &objects,
130  UT_FloatArray &distances2,
131  fpreal radius,
132  int max_objects,
133  const A &accepter) const;
134 
136  const fpreal *pt,
137  UT_ValArray<T> &objects,
138  UT_FloatArray &distances2,
139  fpreal radius,
140  int max_objects) const
141  {
142  return getObjects(
143  pt, objects, distances2, radius, max_objects,
144  UT_STAccepter(true));
145  }
146 
147  /// Returns the subdivision strategy used by this tree.
148  const S *getSubdivisionStrategy() const
149  { return mySubDStrategy; }
150 
151 private:
152  class Node : public UT_SmallObject<Node, UT_SMALLOBJECT_CLEANPAGES_ON, 1024>
153  {
154  public:
155  Node();
156  ~Node();
157 
158  bool isLeaf() const { return !myChildren; }
159  void steal(Node &node);
160 
161  fpreal myCenter[N];
162  UT_ValArray<T> myObjects;
163  Node *myChildren;
164  };
165 
166  struct NodeChildPointer
167  {
168  NodeChildPointer(const Node &node) :
169  myNode(node), myChildIndex(0)
170  {}
171 
172  const Node &myNode;
173  int myChildIndex;
174  };
175 
176 public:
177  void dump(std::ostream &os)
178  {
179  os << "UT_SpatialTree:\n";
180  if (!myRoot)
181  {
182  os << "<<Empty>>\n";
183  return;
184  }
185  std::stack<NodeChildPointer> nodes;
186  nodes.push(NodeChildPointer(*myRoot));
187  while(!nodes.empty())
188  {
189  const Node &node = nodes.top().myNode;
190  int &child_idx = nodes.top().myChildIndex;
191  if (child_idx == 0)
192  {
193  for (size_t i = 0; i < nodes.size(); i++) os << " ";
194  os << "Node [" << &node << "] - Object count: "
195  << node.myObjects.entries() << "\n";
196  }
197  if (child_idx < BranchFactor && !node.isLeaf())
198  nodes.push(NodeChildPointer(node.myChildren[child_idx++]));
199  else
200  nodes.pop();
201  }
202  }
203 
205  {
206  public:
207  bool atEnd() const { return myNodeStack.empty(); }
208  void advance()
209  {
210  // Advance, depth-first
211  while(!myNodeStack.empty())
212  {
213  NodeChildPointer &top = myNodeStack.top();
214  if (top.myNode.isLeaf() ||
215  top.myChildIndex >= BranchFactor)
216  {
217  myNodeStack.pop();
218  }
219  else
220  {
221  myNodeStack.push(NodeChildPointer(
222  top.myNode.myChildren[top.myChildIndex++]));
223  break;
224  }
225 
226  }
227  }
228 
229  const UT_ValArray<T> &objects() const
230  {
231  static UT_ValArray<T> empty;
232  if (myNodeStack.empty())
233  return empty;
234 
235  const NodeChildPointer &top = myNodeStack.top();
236  return top.myNode.myObjects;
237  }
238 
239  NodeIterator &operator++() { advance(); return *this; }
240 
241  bool operator==(const NodeIterator &other)
242  {
243  if (myNodeStack.empty() || other.myNodeStack.empty())
244  return myNodeStack.empty() && other.myNodeStack.empty();
245 
246  NodeChildPointer &top = myNodeStack.top();
247  NodeChildPointer &other_top = other.myNodeStack.top();
248  return top.myNode == other_top.myNode &&
249  top.myChildIndex == other_top.myChildIndex;
250  }
251 
252  bool operator!=(const NodeIterator &other)
253  {
254  return !(*this == other);
255  }
256 
257  protected:
258  friend class UT_SpatialTree;
259 
260  NodeIterator(const Node *root)
261  {
262  if (root)
263  myNodeStack.push(NodeChildPointer(*root));
264  }
265 
266  private:
267  std::stack<NodeChildPointer> myNodeStack;
268  };
269 
271  {
272  return NodeIterator(myRoot);
273  }
274 
276  {
277  return NodeIterator(NULL);
278  }
279 
280 
281 private:
282  /// Expand the tree so that it contains the specified bounds.
283  void expand(const fpreal *min, const fpreal *max);
284 
285  /// Compute the union of the two bounds, storing the result in
286  /// min1/max1.
287  void expand(
288  fpreal *min1,
289  fpreal *max1,
290  const fpreal *min2,
291  const fpreal *max2) const;
292 
293  /// Return the center coordinates, in center, of the specified bounds.
294  void getCenter(
295  const fpreal *min,
296  const fpreal *max,
297  fpreal *center) const;
298 
299  /// Return the bounding rectangle of the specified quadrant (childno)
300  /// of the specified bounds.
301  void getChildBounds(
302  int childno,
303  const fpreal *min,
304  const fpreal *max,
305  const fpreal *center,
306  fpreal *child_min,
307  fpreal *child_max) const;
308 
309  /// Return the index of the child node that contains the specified
310  /// point. min/max specify the bounds of the parent node.
311  int getChildContainer(
312  const fpreal *center,
313  const fpreal *pt) const;
314 
315  /// Given two masks, less and more, that indicate which children
316  /// of a node we would like to iterate over, return the first such
317  /// child. Use this, along with getNextChild(), to iterate over
318  /// all children indicated by the masks. Each bit of the masks
319  /// possibly filters out half of a certain dimension. For example,
320  /// if bit 0 of less is set to 0, the first half of dimension 0
321  /// is filtered out; otherwise, it is included in the iteration.
322  /// If a certain bit of less and the same bit of more are both 1,
323  /// no filtering happens for that dimension.
324  int getFirstChild(int less, int more) const;
325 
326  /// Returns the next child index in the supplied childno parameter
327  /// and returns true if the end of the iteration has not yet been
328  /// reached.
329  bool getNextChild(int &childno, int less, int more) const;
330 
331  /// Initialize the children of the specified node, computing their
332  /// center coordinates from the supplied bounds (min/max).
333  void initChildren(
334  Node *node,
335  const fpreal *min,
336  const fpreal *max) const;
337 
338  /// Returns true if the bounding region specified by min1/max1 contains
339  /// the bounding region specified by min2/max2.
340  bool isInside(
341  const fpreal *min1,
342  const fpreal *max1,
343  const fpreal *pt) const;
344 
345  /// Returns true if the bounding region specified by min1/max1 contains
346  /// the bounding region specified by min2/max2.
347  bool isInside(
348  const fpreal *min1,
349  const fpreal *max1,
350  const fpreal *min2,
351  const fpreal *max2) const;
352 
353  /// Store the ordering (ascending) of the n specified values in indices.
354  /// For example, values[indices[0]] will be the smallest value. This
355  /// function does not move elements in values.
356  static void sort(const fpreal *values, int n, int *indices);
357 
358  /// Subdivide the specified node whose depth and bounds are given.
359  void subdivide(
360  Node *node,
361  int depth,
362  const fpreal *min,
363  const fpreal *max);
364 
365  S *mySubDStrategy;
366  B *myBoundsProvider;
367  Node *myRoot;
368  fpreal myMin[N];
369  fpreal myMax[N];
370  int myDepth;
371 };
372 
373 /// A basic subdivision strategy for use with UT_SpatialTree that subdivides
374 /// based on three criteria: node size, node object count, and node depth.
376 {
377 public:
378  /// Construct an instance of UT_STBasicSubD that subdivides a node
379  /// when all of the following are true:
380  /// * the node's width and height are greater than min_size
381  /// * the node has more than min_objects children
382  /// * the node's depth is less than max_depth
384  fpreal min_size,
385  int min_objects,
386  int max_depth)
387  : myMinSize(min_size)
388  , myMinObjects(min_objects)
389  , myMaxDepth(max_depth)
390  {}
391 
392  /// Returns true if the node at the specified depth with the
393  /// specified bounds provider, containing the specified
394  /// objects, whose bounds are given, should be subdivided.
395  template<class T, class B>
396  bool subdivide(
397  B *bounds_provider,
398  const fpreal *min,
399  const fpreal *max,
400  const UT_ValArray<T> &objects,
401  int depth,
402  int ndimensions)
403  {
404  int idimension;
405 
406  if (depth >= myMaxDepth ||
407  objects.entries() <= myMinObjects)
408  return false;
409 
410  for (idimension = 0; idimension < ndimensions;++idimension)
411  {
412  if (max[idimension] - min[idimension] < myMinSize)
413  return false;
414  }
415 
416  return true;
417  }
418 
419 private:
420  fpreal myMinSize;
421  int myMinObjects;
422  int myMaxDepth;
423 };
424 
425 template<class T, class S, class B, int N>
427  : myChildren(0)
428 {
429 }
430 
431 template<class T, class S, class B, int N>
433 {
434  delete[] myChildren;
435 }
436 
437 template<class T, class S, class B, int N>
438 void
440 {
441  if (this != &node)
442  {
443  int iobject;
444 
445  delete[] myChildren;
446  myObjects.entries(0);
447 
448  ::memcpy(myCenter, node.myCenter, N * sizeof(fpreal));
449 
450  myChildren = node.myChildren;
451  node.myChildren = 0;
452 
453  for (iobject = 0; iobject < node.myObjects.entries(); ++iobject)
454  myObjects.append(node.myObjects(iobject));
455  }
456 }
457 
458 template<class T, class S, class B, int N>
460  B *bounds_provider,
461  const fpreal *min,
462  const fpreal *max)
463  : myRoot(0)
464  , mySubDStrategy(subd_strategy)
465  , myBoundsProvider(bounds_provider)
466  , myDepth(0)
467 {
468  UT_ASSERT(subd_strategy);
469  UT_ASSERT(0 < N && N <= getMaxDimensions());
470 
471  ::memcpy(myMin, min, N * sizeof(fpreal));
472  ::memcpy(myMax, max, N * sizeof(fpreal));
473 }
474 
475 template<class T, class S, class B, int N>
476 UT_SpatialTree<T, S, B, N>::UT_SpatialTree(S *subd_strategy, B *bounds_provider)
477  : myRoot(0)
478  , mySubDStrategy(subd_strategy)
479  , myBoundsProvider(bounds_provider)
480  , myDepth(0)
481 {
482  UT_ASSERT(subd_strategy);
483  UT_ASSERT(0 < N && N <= getMaxDimensions());
484 
485  ::memset(myMin, 0, N * sizeof(fpreal));
486  ::memset(myMax, 0, N * sizeof(fpreal));
487 }
488 
489 template<class T, class S, class B, int N>
491 {
492  delete myRoot;
493 }
494 
495 template<class T, class S, class B, int N>
496 void
498 {
499  Node *node;
500  const fpreal *bounds_min, *bounds_max;
501  fpreal bounds_min_buf[N], bounds_max_buf[N];
502  fpreal min[N], max[N];
503  fpreal child_min[N], child_max[N];
504  int childno, depth;
505 
506  // Get the object's bounds.
507  bounds_min = myBoundsProvider->getMin(object, bounds_min_buf);
508  bounds_max = myBoundsProvider->getMax(object, bounds_max_buf);
509 
510  // Expand the octree until it contains the bounds.
511  if (!isInside(myMin, myMax, bounds_min, bounds_max))
512  expand(bounds_min, bounds_max);
513 
514  UT_ASSERT(isInside(myMin, myMax, bounds_min, bounds_max));
515 
516  // Initialize min/max with the octree's bounds.
517  ::memcpy(min, myMin, N * sizeof(fpreal));
518  ::memcpy(max, myMax, N * sizeof(fpreal));
519 
520  depth = 1;
521 
522  // If we don't have a root node yet, create one.
523  if (!myRoot)
524  {
525  myRoot = node = new Node();
526  getCenter(myMin, myMax, myRoot->myCenter);
527  myDepth = 1;
528  }
529  else
530  {
531  // Search the tree until we find a leaf node or a non-leaf node
532  // that contains the entire bounds but has no children that
533  // contain the entire bounds.
534  node = myRoot;
535 
536  while (!node->isLeaf())
537  {
538  childno = getChildContainer(node->myCenter, bounds_min);
539  if (childno != getChildContainer(node->myCenter, bounds_max))
540  break;
541 
542  getChildBounds(
543  childno, min, max, node->myCenter, child_min, child_max);
544  node = &node->myChildren[childno];
545  UT_ASSERT(isInside(child_min,child_max,bounds_min,bounds_max));
546  ::memcpy(min, child_min, N * sizeof(fpreal));
547  ::memcpy(max, child_max, N * sizeof(fpreal));
548  depth++;
549  }
550  }
551 
552  UT_ASSERT(node);
553 
554  // Add the object and its bounds to the node.
555  node->myObjects.append(object);
556 
557  // Possibly subdivide the node.
558  if (node->isLeaf())
559  subdivide(node, depth, min, max);
560 }
561 
562 template<class T, class S, class B, int N>
563 void
565 {
566  Node *root;
567  fpreal center[2];
568  fpreal size;
569  int childno, idimension;
570 
571  // If we haven't created a tree yet, simply expand myMin/myMax until
572  // they contain min/max.
573  if (!myRoot)
574  {
575  ::memcpy(myMin, min, N * sizeof(fpreal));
576  ::memcpy(myMax, max, N * sizeof(fpreal));
577  }
578  else if (myRoot->isLeaf())
579  {
580  expand(myMin, myMax, min, max);
581  getCenter(myMin, myMax, myRoot->myCenter);
582  }
583  else
584  {
585  // Add one level at a time to our octree until it contains min/max.
586  while (!isInside(myMin, myMax, min, max))
587  {
588  root = new Node();
589  childno = 0;
590 
591  for (idimension = 0; idimension < N; ++idimension)
592  {
593  // Compute the center of myMin/myMax and min/max for the
594  // current dimension so that we can determine which direction
595  // to expand the octree in.
596  size = myMax[idimension] - myMin[idimension];
597  center[0] = (myMin[idimension] + myMax[idimension]) / 2.0;
598  center[1] = (min[idimension] + max[idimension]) / 2.0;
599 
600  if (center[0] < center[1])
601  {
602  // Expand to the right.
603  childno |= 1 << idimension;
604  root->myCenter[idimension] = myMax[idimension];
605  myMax[idimension] += size;
606  }
607  else
608  {
609  // Expand to the left.
610  root->myCenter[idimension] = myMin[idimension];
611  myMin[idimension] -= size;
612  }
613  }
614 
615  // Initialize the children of the new root node and assign
616  // the data from the old root node to the appropriate child.
617  initChildren(root, myMin, myMax);
618  root->myChildren[childno].steal(*myRoot);
619  delete myRoot;
620  myRoot = root;
621  myDepth++;
622  }
623  }
624 }
625 
626 template<class T, class S, class B, int N>
627 void
629  fpreal *max1,
630  const fpreal *min2,
631  const fpreal *max2) const
632 {
633  int idimension;
634 
635  for (idimension = 0; idimension < N; ++idimension)
636  {
637  min1[idimension] = SYSmin(min1[idimension], min2[idimension]);
638  max1[idimension] = SYSmax(max1[idimension], max2[idimension]);
639  }
640 }
641 
642 template<class T, class S, class B, int N>
643 void
645  const fpreal *max,
646  fpreal *center) const
647 {
648  int idimension;
649 
650  for (idimension = 0; idimension < N; ++idimension)
651  center[idimension] = (min[idimension] + max[idimension]) / 2.0;
652 }
653 
654 template<class T, class S, class B, int N>
655 void
657  const fpreal *min,
658  const fpreal *max,
659  const fpreal *center,
660  fpreal *child_min,
661  fpreal *child_max) const
662 {
663  int idimension;
664 
665  for (idimension = 0; idimension < N; ++idimension)
666  {
667  if ((childno & (1 << idimension)) != 0)
668  {
669  child_min[idimension] = min[idimension];
670  child_max[idimension] = center[idimension];
671  }
672  else
673  {
674  child_min[idimension] = center[idimension];
675  child_max[idimension] = max[idimension];
676  }
677  }
678 }
679 
680 template<class T, class S, class B, int N>
681 int
683  const fpreal *pt) const
684 {
685  int childno, idimension;
686 
687  childno = 0;
688 
689  for (idimension = 0; idimension < N; ++idimension)
690  {
691  if (pt[idimension] < center[idimension])
692  childno |= 1 << idimension;
693  }
694 
695  UT_ASSERT(0 <= childno && childno < BranchFactor);
696 
697  return childno;
698 }
699 
700 template<class T, class S, class B, int N>
701 int
703 {
704  int free, childno;
705 
706  // Find all dimensions that may be less or more. Every other dimension
707  // is either less or more (but must be at least one).
708  free = less & more;
709 
710  // Remove the free dimensions from both less and more.
711  less &= ~free;
712  more &= ~free;
713 
714  // Ensure that the non-free dimensions are set correctly. Set all the
715  // free dimensions to 0.
716  childno = 0;
717  childno |= less;
718  childno &= ~more;
719 
720  UT_ASSERT(childno < BranchFactor);
721 
722  return childno;
723 }
724 
725 template<class T, class S, class B, int N>
726 bool
727 UT_SpatialTree<T, S, B, N>::getNextChild(int &childno, int less, int more) const
728 {
729  int free, mask, idimension;
730 
731  // Find all dimensions that may be less or more. Every other dimension
732  // is either less or more (but must be at least one).
733  free = less & more;
734 
735  // Find the first free dimension that is set to 0 and flip it. Flip
736  // all free dimensions occuring before this back to 0.
737  for (idimension = 0; idimension < N; ++idimension)
738  {
739  mask = 1 << idimension;
740  UT_ASSERT((less & mask) != 0 || (more & mask) != 0);
741 
742  if ((free & mask) != 0)
743  {
744  childno ^= mask;
745  if ((childno & mask) != 0)
746  break;
747  }
748  }
749 
750  UT_ASSERT(childno < BranchFactor || idimension == N);
751 
752  return idimension < N;
753 }
754 
755 template<class T, class S, class B, int N>
756 template<class A>
757 void
759  UT_ValArray<T> &objects,
760  const A &accepter) const
761 {
762  Node *node;
763  const fpreal *bounds_min, *bounds_max;
764  fpreal bounds_min_buf[N], bounds_max_buf[N];
765  int iobject, childno;
766 
767  // Search the tree for objects containing the specified point.
768  if (!myRoot)
769  return;
770 
771  node = myRoot;
772 
773  while (node)
774  {
775  // Check each object in the current node to see if it contains
776  // the query position.
777  for (iobject = 0; iobject < node->myObjects.entries(); ++iobject)
778  {
779  // Get the object's bounds.
780  bounds_min = myBoundsProvider->getMin(
781  node->myObjects(iobject), bounds_min_buf);
782  bounds_max = myBoundsProvider->getMax(
783  node->myObjects(iobject), bounds_max_buf);
784 
785  if (isInside(bounds_min, bounds_max, pt) &&
786  accepter.accept(node->myObjects(iobject)))
787  objects.append(node->myObjects(iobject));
788  }
789 
790  if (node->isLeaf())
791  break;
792 
793  // Get the next node.
794  childno = getChildContainer(node->myCenter, pt);
795  node = &node->myChildren[childno];
796  }
797 }
798 
799 template<class T, class S, class B, int N>
800 template<class A>
801 int
803  UT_ValArray<T> &objects,
804  UT_FloatArray &distances2,
805  fpreal radius,
806  int max_objects,
807  const A &accepter) const
808 {
809  Node **node_stack;
810  fpreal *dist_stack;
811  int *childnos;
812  fpreal *child_dists;
813  int *indices;
814  Node *node;
815  const fpreal *center;
816  fpreal dim_dists2[N];
817  fpreal radius2, dist2;
818  int iobject, idimension, ichild;
819  int childno, ptchildno, index, nchildren;
820  int max_stack_size, stack_size, less, more, diff;
821 
822  objects.entries(0);
823  distances2.entries(0);
824 
825  // If the caller requested 0 (or fewer) objects, we are done.
826  if (!myRoot || max_objects <= 0)
827  return 0;
828 
829  // Allocate stack buffers.
830  max_stack_size = (myDepth - 1) * BranchFactor + 1;
831  node_stack = (Node **) UTstackAlloc(max_stack_size * sizeof(Node *));
832  dist_stack = (fpreal *) UTstackAlloc(max_stack_size * sizeof(fpreal));
833  child_dists = (fpreal *) UTstackAlloc(BranchFactor * sizeof(fpreal));
834  childnos = (int *) UTstackAlloc(2 * BranchFactor * sizeof(int));
835  indices = childnos + BranchFactor;
836 
837  // Validate the radius and compute the squared radius.
838  radius = SYSmin(SYSmax(0.0, radius), 1e18);
839  radius2 = radius * radius;
840 
841  // Push the root node onto the node stack.
842  node_stack[0] = myRoot;
843  dist_stack[0] = 0.0;
844  stack_size = 1;
845 
846  while (stack_size > 0)
847  {
848  // Retrieve the next node. We perform a distance check again
849  // (we did this previously when we inserted the node into the
850  // stack) because our radius may have shrunk. This can drastically
851  // improve performance because we can avoid touching many nodes.
852  while (stack_size > 0 && dist_stack[--stack_size] >= radius2)
853  ;
854 
855  if (stack_size < 0)
856  break;
857 
858  node = node_stack[stack_size];
859  UT_ASSERT(node);
860  center = node->myCenter;
861 
862  // Check to see if the objects in the current node are within
863  // the specified radius.
864  for (iobject = 0; iobject < node->myObjects.entries(); ++iobject)
865  {
866  dist2 = myBoundsProvider->getDistance2(
867  node->myObjects(iobject), pt);
868 
869  if (dist2 < radius2 && accepter.accept(node->myObjects(iobject)))
870  {
871  // Determine the position of the object in the object list.
872  index = distances2.sortedInsert(dist2);
873  UT_ASSERT(index >= 0);
874  objects.insert(node->myObjects(iobject), index);
875  objects.truncate(max_objects);
876  distances2.truncate(max_objects);
877 
878  // If we already have the number of requested objects
879  // (max_objects), then we can come up with a new, lower bound
880  // for the maximum distance that an object can be from
881  // the query position.
882  if (objects.entries() == max_objects)
883  radius2 = distances2.last();
884  }
885  }
886 
887  if (!node->isLeaf())
888  {
889  // Determine which child nodes may possibly be within the
890  // specified radius of the specified point.
891  less = more = 0;
892 
893  for (idimension = 0; idimension < N; ++idimension)
894  {
895  dim_dists2[idimension] = pt[idimension] - center[idimension];
896  dim_dists2[idimension] *= dim_dists2[idimension];
897 
898  if (dim_dists2[idimension] <= radius2)
899  {
900  less |= 1 << idimension;
901  more |= 1 << idimension;
902  }
903  else if (pt[idimension] < center[idimension])
904  less |= 1 << idimension;
905  else
906  more |= 1 << idimension;
907  }
908 
909  // Sort the children by their distance to the supplied point
910  // (measured from the center of the child).
911  ptchildno = getChildContainer(center, pt);
912  childno = getFirstChild(less, more);
913  nchildren = 0;
914 
915  do
916  {
917  diff = ptchildno ^ childno;
918  dist2 = 0.0;
919 
920  for (idimension = 0; idimension < N; ++idimension)
921  {
922  if ((diff & (1 << idimension)) != 0)
923  dist2 += dim_dists2[idimension];
924  }
925 
926  childnos[nchildren] = childno;
927  child_dists[nchildren] = dist2;
928  nchildren++;
929  }
930  while (getNextChild(childno, less, more));
931 
932  sort(child_dists, nchildren, indices);
933 
934  // Add the children to the stack starting with the ones furthest
935  // from the supplied point.
936  for (ichild = nchildren - 1; ichild >= 0; --ichild)
937  {
938  childno = childnos[indices[ichild]];
939  node_stack[stack_size] = &node->myChildren[childno];
940  dist_stack[stack_size] = child_dists[indices[ichild]];
941  stack_size++;
942  UT_ASSERT(stack_size <= max_stack_size);
943  }
944  }
945  }
946 
947  if (node_stack)
948  {
949  UTstackFree(node_stack);
950  UTstackFree(dist_stack);
951  UTstackFree(child_dists);
952  UTstackFree(childnos);
953  }
954 
955  return objects.entries();
956 }
957 
958 template<class T, class S, class B, int N>
959 void
961  const fpreal *min,
962  const fpreal *max) const
963 {
964  fpreal child_min[N], child_max[N];
965  int ichild;
966 
967  UT_ASSERT(node);
968  UT_ASSERT(!node->myChildren);
969  UT_ASSERT(min);
970  UT_ASSERT(max);
971 
972  node->myChildren = new Node[BranchFactor];
973 
974  // Compute the center coordinates for each child.
975  for (ichild = 0; ichild < BranchFactor; ++ichild)
976  {
977  getChildBounds(ichild, min, max, node->myCenter, child_min, child_max);
978  getCenter(child_min, child_max, node->myChildren[ichild].myCenter);
979  }
980 }
981 
982 template<class T, class S, class B, int N>
983 bool
985  const fpreal *max,
986  const fpreal *pt) const
987 {
988  int idimension;
989 
990  for (idimension = 0; idimension < N; ++idimension)
991  {
992  if (pt[idimension] < min[idimension] ||
993  pt[idimension] > max[idimension])
994  return false;
995  }
996 
997  return true;
998 }
999 
1000 template<class T, class S, class B, int N>
1001 bool
1003  const fpreal *max1,
1004  const fpreal *min2,
1005  const fpreal *max2) const
1006 {
1007  int idimension;
1008 
1009  for (idimension = 0; idimension < N; ++idimension)
1010  {
1011  if (min2[idimension] < min1[idimension] ||
1012  max2[idimension] > max1[idimension])
1013  return false;
1014  }
1015 
1016  return true;
1017 }
1018 
1019 template<class T, class S, class B, int N>
1020 void
1022 {
1023  int i, j, tmp;
1024 
1025  for (i = 0; i < n; ++i)
1026  indices[i] = i;
1027 
1028  for (i = 0; i < n; ++i)
1029  {
1030  for (j = i + 1; j < n; ++j)
1031  {
1032  if (values[indices[j]] < values[indices[i]])
1033  {
1034  tmp = indices[i];
1035  indices[i] = indices[j];
1036  indices[j] = tmp;
1037  }
1038  }
1039  }
1040 }
1041 
1042 template<class T, class S, class B, int N>
1043 void
1045  const fpreal *min, const fpreal *max)
1046 {
1047  UT_ASSERT(node);
1048  UT_ASSERT(node->isLeaf());
1049 
1050  Node *child;
1051  const fpreal *bounds_min, *bounds_max;
1052  fpreal bounds_min_buf[N], bounds_max_buf[N];
1053  fpreal child_min[N], child_max[N];
1054  int childno, ichild, iobject;
1055  bool subd;
1056 
1057  // Query the subdivision strategy to determine if the node should
1058  // be subdivided.
1059  UT_ASSERT(mySubDStrategy);
1060  subd = mySubDStrategy->subdivide(
1061  myBoundsProvider,
1062  min,
1063  max,
1064  (const UT_ValArray<T> &) node->myObjects,
1065  depth,
1066  N);
1067 
1068  if (!subd)
1069  return;
1070 
1071  // Increase the tree's depth.
1072  myDepth = SYSmax(myDepth, depth + 1);
1073 
1074  // Initialize the node's children.
1075  initChildren(node, min, max);
1076 
1077  // Test each of the node's objects to see if they can be placed
1078  // in one of its children.
1079  for (iobject = 0; iobject < node->myObjects.entries(); )
1080  {
1081  // Get the object's bounds.
1082  bounds_min = myBoundsProvider->getMin(
1083  node->myObjects(iobject), bounds_min_buf);
1084  bounds_max = myBoundsProvider->getMax(
1085  node->myObjects(iobject), bounds_max_buf);
1086 
1087  // Test for containment.
1088  childno = getChildContainer(node->myCenter, bounds_min);
1089  if (childno == getChildContainer(node->myCenter, bounds_max))
1090  {
1091  child = &node->myChildren[childno];
1092  child->myObjects.append(node->myObjects(iobject));
1093  node->myObjects.removeIndex(iobject);
1094  }
1095  else
1096  iobject++;
1097  }
1098 
1099  // Possibly subdivide each of the children.
1100  for (ichild = 0; ichild < BranchFactor; ++ichild)
1101  {
1102  getChildBounds(ichild, min, max, node->myCenter, child_min, child_max);
1103  subdivide(&node->myChildren[ichild], depth + 1, child_min, child_max);
1104  }
1105 }
1106 
1107 #undef MIN_NODE_SIZE
1108 
1109 #endif
1110 
#define SYSmax(a, b)
Definition: SYS_Math.h:1570
T & last()
Definition: UT_Array.h:816
GLsizei GLenum const void * indices
Definition: glcorearb.h:406
const UT_ValArray< T > & objects() const
const B * getBoundsProvider() const
Returns the bounds provider used by this tree.
NodeIterator leaf_end() const
Definition: Node.h:52
UT_STAccepter(bool accept)
NodeIterator(const Node *root)
void getObjects(const fpreal *pt, UT_ValArray< T > &objects) const
#define UT_API
Definition: UT_API.h:14
bool operator!=(const NodeIterator &other)
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
~UT_SpatialTree()
Destructor.
UT_SpatialTree(S *subd_strategy, B *bounds_provider, const fpreal *min, const fpreal *max)
NodeIterator begin() const
#define UTstackAlloc(size)
Definition: UT_StackAlloc.h:34
bool accept(int objectno) const
GLdouble n
Definition: glcorearb.h:2008
void dump(std::ostream &os)
GLint GLuint mask
Definition: glcorearb.h:124
exint sortedInsert(const T &t, Comparator compare)
Definition: UT_ArrayImpl.h:819
bool operator==(const NodeIterator &other)
exint append()
Definition: UT_Array.h:142
bool subdivide(B *bounds_provider, const fpreal *min, const fpreal *max, const UT_ValArray< T > &objects, int depth, int ndimensions)
GLint GLint GLsizei GLsizei GLsizei depth
Definition: glcorearb.h:476
UT_STBasicSubD(fpreal min_size, int min_objects, int max_depth)
exint entries() const
Alias of size(). size() is preferred.
Definition: UT_Array.h:648
GLint j
Definition: glad.h:2733
GLsizeiptr size
Definition: glcorearb.h:664
static int getMaxDimensions()
Return the maximum number of dimensions allowed.
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1602
fpreal64 fpreal
Definition: SYS_Types.h:277
GLuint index
Definition: glcorearb.h:786
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
void getObjects(const fpreal *pt, UT_ValArray< T > &objects, const A &accepter) const
GA_API const UT_StringHolder N
int getObjects(const fpreal *pt, UT_ValArray< T > &objects, UT_FloatArray &distances2, fpreal radius, int max_objects) const
void truncate(exint maxsize)
Decreases, but never expands, to the given maxsize.
Definition: UT_Array.h:723
GLdouble GLdouble GLdouble top
Definition: glad.h:2817
#define UT_ASSERT(ZZ)
Definition: UT_Assert.h:156
void addObject(T object)
Add the specified object to the tree.
#define SYSmin(a, b)
Definition: SYS_Math.h:1571
void sort(I begin, I end, const Pred &pred)
Definition: pugixml.cpp:7334
const S * getSubdivisionStrategy() const
Returns the subdivision strategy used by this tree.
exint insert(exint index)
Definition: UT_ArrayImpl.h:721
myRoot
Definition: UT_RTreeImpl.h:716
#define UTstackFree(ptr)
Definition: UT_StackAlloc.h:38