1 #ifndef RTREE_H
2 #define RTREE_H
3 
4 // NOTE This file compiles under MSVC 6 SP5 and MSVC .Net 2003 it may not work on other compilers without modification.
5 
6 // NOTE These next few lines may be win32 specific, you may need to modify them to compile on other platform
7 #include <stdio.h>
8 #include <math.h>
9 #include <assert.h>
10 #include <stdlib.h>
11 
12 #define ASSERT assert // RTree uses ASSERT( condition )
13 #ifndef Min
14   #define Min __min
15 #endif //Min
16 #ifndef Max
17   #define Max __max
18 #endif //Max
19 
20 #define rtree_min(a,b) (a<b?a:b)
21 #define rtree_max(a,b) (a>b?a:b)
22 
23 
24 
25 //
26 // RTree.h
27 //
28 
29 #define RTREE_TEMPLATE template<class DATATYPE, class DATATYPENP, class ELEMTYPE, int NUMDIMS, class CONTEXT, class ELEMTYPEREAL, int TMAXNODES, int TMINNODES>
30 #define RTREE_QUAL RTree<DATATYPE, DATATYPENP, ELEMTYPE, NUMDIMS, CONTEXT, ELEMTYPEREAL, TMAXNODES, TMINNODES>
31 
32 #define RTREE_DONT_USE_MEMPOOLS // This version does not contain a fixed memory allocator, fill in lines with EXAMPLE to implement one.
33 #define RTREE_USE_SPHERICAL_VOLUME // Better split classification, may be slower on some systems
34 
35 #undef RTREE_WANTS_IO
36 
37 // Fwd decl
38 #ifdef RTREE_WANTS_IO
39 class RTFileStream;  // File I/O helper class, look below for implementation and notes.
40 #endif
41 
42 
43 /// \class RTree
44 /// Implementation of RTree, a multidimensional bounding rectangle tree.
45 /// Example usage: For a 3-dimensional tree use RTree<Object*, float, 3> myTree;
46 ///
47 /// This modified, templated C++ version by Greg Douglas at Auran (http://www.auran.com)
48 ///
49 /// DATATYPE Referenced data, should be int, void*, obj* etc. no larger than sizeof<void*> and simple type
50 /// ELEMTYPE Type of element such as int or float
51 /// NUMDIMS Number of dimensions such as 2 or 3
52 /// ELEMTYPEREAL Type of element that allows fractional and large values such as float or double, for use in volume calcs
53 ///
54 /// NOTES: Inserting and removing data requires the knowledge of its constant Minimal Bounding Rectangle.
55 ///        This version uses new/delete for nodes, I recommend using a fixed size allocator for efficiency.
56 ///        Instead of using a callback function for returned results, I recommend and efficient pre-sized, grow-only memory
57 ///        array similar to MFC CArray or STL Vector for returning search query result.
58 ///
59 template<class DATATYPE, class DATATYPENP, class ELEMTYPE, int NUMDIMS, class CONTEXT,
60          class ELEMTYPEREAL = ELEMTYPE, int TMAXNODES = 8, int TMINNODES = TMAXNODES / 2>
61 class RTree
62 {
63 protected:
64 
65   struct Node;  // Fwd decl.  Used by other internal structs and iterator
66 
67 public:
68 
69   // These constant must be declared after Branch and before Node struct
70   // Stuck up here for MSVC 6 compiler.  NSVC .NET 2003 is much happier.
71   enum
72   {
73     MAXNODES = TMAXNODES,                         ///< Max elements in node
74     MINNODES = TMINNODES                         ///< Min elements in node
75   };
76 
77 
78 public:
79   typedef void(DATATYPENP::* Operation)(const CONTEXT &) const;
80 
81   RTree(Operation operation);
82   virtual ~RTree();
83 
84   /// Insert entry
85   /// \param a_min Min of bounding rect
86   /// \param a_max Max of bounding rect
87   /// \param a_dataId Positive Id of data.  Maybe zero, but negative numbers not allowed.
88   virtual void Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE& a_dataId);
89 
90   /// Remove entry
91   /// \param a_min Min of bounding rect
92   /// \param a_max Max of bounding rect
93   /// \param a_dataId Positive Id of data.  Maybe zero, but negative numbers not allowed.
94   virtual void Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE& a_dataId);
95 
96 
97   /// DK 15.10.2008 - begin
98   //typedef void(DATATYPENP::* Operation)(const CONTEXT &) const;
99 
100   /// Find all within search rectangle
101   /// \param a_min Min of search bounding rect
102   /// \param a_max Max of search bounding rect
103   /// \param a_searchResult Search result array.  Caller should set grow size. Function will reset, not append to array.
104   /// \param a_resultCallback Callback function to return result.  Callback should return 'true' to continue searching
105   /// \param a_context User context to pass as parameter to a_resultCallback
106   /// \return Returns the number of entries found
107   virtual int Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const CONTEXT &c) const;
108 
109   /// DK 15.10.2008 - end
110 
111   /// Remove all entries from tree
112   void RemoveAll();
113 
114   /// Count the data elements in this container.  This is slow as no internal counter is maintained.
115   int Count();
116 
117 #ifdef RTREE_WANTS_IO
118   /// Load tree contents from file
119   bool Load(const char* a_fileName);
120   /// Load tree contents from stream
121   bool Load(RTFileStream& a_stream);
122 
123 
124   /// Save tree contents to file
125   bool Save(const char* a_fileName);
126   /// Save tree contents to stream
127   bool Save(RTFileStream& a_stream);
128 #endif
129 
130   /// Iterator is not remove safe.
131   class Iterator
132   {
133   private:
134 
135     enum { MAX_STACK = 32 }; //  Max stack size. Allows almost n^32 where n is number of branches in node
136 
137     struct StackElement
138     {
139       Node* m_node;
140       int m_branchIndex;
141     };
142 
143   public:
144 
Iterator()145     Iterator()                                    { Init(); }
146 
~Iterator()147     ~Iterator()                                   { }
148 
149     /// Is iterator invalid
IsNull()150     bool IsNull()                                 { return (m_tos <= 0); }
151 
152     /// Is iterator pointing to valid data
IsNotNull()153     bool IsNotNull()                              { return (m_tos > 0); }
154 
155     /// Access the current data element. Caller must be sure iterator is not NULL first.
156     DATATYPE& operator*()
157     {
158       ASSERT(IsNotNull());
159       StackElement& curTos = m_stack[m_tos - 1];
160       return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;
161     }
162 
163     /// Access the current data element. Caller must be sure iterator is not NULL first.
164     const DATATYPE& operator*() const
165     {
166       ASSERT(IsNotNull());
167       StackElement& curTos = m_stack[m_tos - 1];
168       return curTos.m_node->m_branch[curTos.m_branchIndex].m_data;
169     }
170 
171     /// Find the next data element
172     bool operator++()                             { return FindNextData(); }
173 
174   private:
175 
176     /// Reset iterator
Init()177     void Init()                                   { m_tos = 0; }
178 
179     /// Find the next data element in the tree (For internal use only)
FindNextData()180     bool FindNextData()
181     {
182       for(;;)
183       {
184         if(m_tos <= 0)
185         {
186           return false;
187         }
188         StackElement curTos = Pop(); // Copy stack top cause it may change as we use it
189 
190         if(curTos.m_node->IsLeaf())
191         {
192           // Keep walking through data while we can
193           if(curTos.m_branchIndex+1 < curTos.m_node->m_count)
194           {
195             // There is more data, just point to the next one
196             Push(curTos.m_node, curTos.m_branchIndex + 1);
197             return true;
198           }
199           // No more data, so it will fall back to previous level
200         }
201         else
202         {
203           if(curTos.m_branchIndex+1 < curTos.m_node->m_count)
204           {
205             // Push sibling on for future tree walk
206             // This is the 'fall back' node when we finish with the current level
207             Push(curTos.m_node, curTos.m_branchIndex + 1);
208           }
209           // Since cur node is not a leaf, push first of next level to get deeper into the tree
210           Node* nextLevelnode = curTos.m_node->m_branch[curTos.m_branchIndex].m_child;
211           Push(nextLevelnode, 0);
212 
213           // If we pushed on a new leaf, exit as the data is ready at TOS
214           if(nextLevelnode->IsLeaf())
215           {
216             return true;
217           }
218         }
219       }
220     }
221 
222     /// Push node and branch onto iteration stack (For internal use only)
Push(Node * a_node,int a_branchIndex)223     void Push(Node* a_node, int a_branchIndex)
224     {
225       m_stack[m_tos].m_node = a_node;
226       m_stack[m_tos].m_branchIndex = a_branchIndex;
227       ++m_tos;
228       ASSERT(m_tos <= MAX_STACK);
229     }
230 
231     /// Pop element off iteration stack (For internal use only)
Pop()232     StackElement& Pop()
233     {
234       ASSERT(m_tos > 0);
235       --m_tos;
236       return m_stack[m_tos];
237     }
238 
239     StackElement m_stack[MAX_STACK];              ///< Stack as we are doing iteration instead of recursion
240     int m_tos;                                    ///< Top Of Stack index
241 
242     friend class RTree<DATATYPE, DATATYPENP, ELEMTYPE, NUMDIMS, CONTEXT, ELEMTYPEREAL, TMAXNODES, TMINNODES>; // Allow hiding of non-public functions while allowing manipulation by logical owner
243   };
244 
245   /// Get 'first' for iteration
GetFirst(Iterator & a_it)246   void GetFirst(Iterator& a_it)
247   {
248     a_it.Init();
249     if(m_root && (m_root->m_count > 0))
250     {
251       a_it.Push(m_root, 0);
252       a_it.FindNextData();
253     }
254   }
255 
256   /// Get Next for iteration
GetNext(Iterator & a_it)257   void GetNext(Iterator& a_it)                    { ++a_it; }
258 
259   /// Is iterator NULL, or at end?
IsNull(Iterator & a_it)260   bool IsNull(Iterator& a_it)                     { return a_it.IsNull(); }
261 
262   /// Get object at iterator position
GetAt(Iterator & a_it)263   DATATYPE& GetAt(Iterator& a_it)                 { return *a_it; }
264 
265 protected:
266 
267   /// Minimal bounding rectangle (n-dimensional)
268   struct Rect
269   {
270     ELEMTYPE m_min[NUMDIMS];                      ///< Min dimensions of bounding box
271     ELEMTYPE m_max[NUMDIMS];                      ///< Max dimensions of bounding box
272   };
273 
274   /// May be data or may be another subtree
275   /// The parents level determines this.
276   /// If the parents level is 0, then this is data
277   struct Branch
278   {
279     Rect m_rect;                                  ///< Bounds
280     union
281     {
282       Node* m_child;                              ///< Child node
283       DATATYPE m_data;                            ///< Data Id or Ptr
284     };
285   };
286 
287   /// Node for each branch level
288   struct Node
289   {
IsInternalNodeNode290     bool IsInternalNode() const                   { return (m_level > 0); } // Not a leaf, but a internal node
IsLeafNode291     bool IsLeaf() const                           { return (m_level == 0); } // A leaf, contains data
292 
293     int m_count;                                  ///< Count
294     int m_level;                                  ///< Leaf is zero, others positive
295     Branch m_branch[MAXNODES];                    ///< Branch
296   };
297 
298   /// A link list of nodes for reinsertion after a delete operation
299   struct ListNode
300   {
301     ListNode* m_next;                             ///< Next in list
302     Node* m_node;                                 ///< Node
303   };
304 
305   /// Variables for finding a split partition
306   struct PartitionVars
307   {
308     int m_partition[MAXNODES+1];
309     int m_total;
310     int m_minFill;
311     int m_taken[MAXNODES+1];
312     int m_count[2];
313     Rect m_cover[2];
314     ELEMTYPEREAL m_area[2];
315 
316     Branch m_branchBuf[MAXNODES+1];
317     int m_branchCount;
318     Rect m_coverSplit;
319     ELEMTYPEREAL m_coverSplitArea;
320   };
321 
322   Node* AllocNode();
323   void FreeNode(Node* a_node);
324   void InitNode(Node* a_node);
325   void InitRect(Rect* a_rect);
326   bool InsertRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node, Node** a_newNode, int a_level);
327   bool InsertRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root, int a_level);
328   Rect NodeCover(Node* a_node);
329   bool AddBranch(Branch* a_branch, Node* a_node, Node** a_newNode);
330   void DisconnectBranch(Node* a_node, int a_index);
331   int PickBranch(Rect* a_rect, Node* a_node);
332   Rect CombineRect(Rect* a_rectA, Rect* a_rectB);
333   void SplitNode(Node* a_node, Branch* a_branch, Node** a_newNode);
334   ELEMTYPEREAL RectSphericalVolume(Rect* a_rect);
335   ELEMTYPEREAL RectVolume(Rect* a_rect);
336   ELEMTYPEREAL CalcRectVolume(Rect* a_rect);
337   void GetBranches(Node* a_node, Branch* a_branch, PartitionVars* a_parVars);
338   void ChoosePartition(PartitionVars* a_parVars, int a_minFill);
339   void LoadNodes(Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVars);
340   void InitParVars(PartitionVars* a_parVars, int a_maxRects, int a_minFill);
341   void PickSeeds(PartitionVars* a_parVars);
342   void Classify(int a_index, int a_group, PartitionVars* a_parVars);
343   bool RemoveRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root);
344   bool RemoveRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node, ListNode** a_listNode);
345   ListNode* AllocListNode();
346   void FreeListNode(ListNode* a_listNode);
347   bool Overlap(Rect* a_rectA, Rect* a_rectB) const;
348   void ReInsert(Node* a_node, ListNode** a_listNode);
349   bool Search(Node* a_node, Rect* a_rect, int& a_foundCount, const CONTEXT &c) const;
350   void RemoveAllRec(Node* a_node);
351   void Reset();
352   void CountRec(Node* a_node, int& a_count);
353 
354 #ifdef RTREE_WANTS_IO
355   bool SaveRec(Node* a_node, RTFileStream& a_stream);
356   bool LoadRec(Node* a_node, RTFileStream& a_stream);
357 #endif
358 
359   Node* m_root;                                    ///< Root of tree
360   ELEMTYPEREAL m_unitSphereVolume;                 ///< Unit sphere constant for required number of dimensions
361   Operation myOperation;
362 };
363 
364 
365 #ifdef RTREE_WANTS_IO
366 // Because there is not stream support, this is a quick and dirty file I/O helper.
367 // Users will likely replace its usage with a Stream implementation from their favorite API.
368 class RTFileStream
369 {
370   FILE* m_file;
371 
372 public:
373 
374 
RTFileStream()375   RTFileStream()
376   {
377     m_file = NULL;
378   }
379 
~RTFileStream()380   ~RTFileStream()
381   {
382     Close();
383   }
384 
OpenRead(const char * a_fileName)385   bool OpenRead(const char* a_fileName)
386   {
387     m_file = fopen(a_fileName, "rb");
388     if(!m_file)
389     {
390       return false;
391     }
392     return true;
393   }
394 
OpenWrite(const char * a_fileName)395   bool OpenWrite(const char* a_fileName)
396   {
397     m_file = fopen(a_fileName, "wb");
398     if(!m_file)
399     {
400       return false;
401     }
402     return true;
403   }
404 
Close()405   void Close()
406   {
407     if(m_file)
408     {
409       fclose(m_file);
410       m_file = NULL;
411     }
412   }
413 
414   template< typename TYPE >
Write(const TYPE & a_value)415   size_t Write(const TYPE& a_value)
416   {
417     ASSERT(m_file);
418     return fwrite((void*)&a_value, sizeof(a_value), 1, m_file);
419   }
420 
421   template< typename TYPE >
WriteArray(const TYPE * a_array,int a_count)422   size_t WriteArray(const TYPE* a_array, int a_count)
423   {
424     ASSERT(m_file);
425     return fwrite((void*)a_array, sizeof(TYPE) * a_count, 1, m_file);
426   }
427 
428   template< typename TYPE >
Read(TYPE & a_value)429   size_t Read(TYPE& a_value)
430   {
431     ASSERT(m_file);
432     return fread((void*)&a_value, sizeof(a_value), 1, m_file);
433   }
434 
435   template< typename TYPE >
ReadArray(TYPE * a_array,int a_count)436   size_t ReadArray(TYPE* a_array, int a_count)
437   {
438     ASSERT(m_file);
439     return fread((void*)a_array, sizeof(TYPE) * a_count, 1, m_file);
440   }
441 };
442 #endif
443 
444 
445 RTREE_TEMPLATE
RTree(Operation operation)446 RTREE_QUAL::RTree(Operation operation)
447 : myOperation(operation)
448 {
449   ASSERT(MAXNODES > MINNODES);
450   ASSERT(MINNODES > 0);
451 
452 
453   // We only support machine word size simple data type eg. integer index or object pointer.
454   // Since we are storing as union with non data branch
455   ASSERT(sizeof(DATATYPE) == sizeof(void*) || sizeof(DATATYPE) == sizeof(int));
456 
457   // Precomputed volumes of the unit spheres for the first few dimensions
458   const float UNIT_SPHERE_VOLUMES[] = {
459     0.000000f, 2.000000f, 3.141593f, // Dimension  0,1,2
460     4.188790f, 4.934802f, 5.263789f, // Dimension  3,4,5
461     5.167713f, 4.724766f, 4.058712f, // Dimension  6,7,8
462     3.298509f, 2.550164f, 1.884104f, // Dimension  9,10,11
463     1.335263f, 0.910629f, 0.599265f, // Dimension  12,13,14
464     0.381443f, 0.235331f, 0.140981f, // Dimension  15,16,17
465     0.082146f, 0.046622f, 0.025807f, // Dimension  18,19,20
466   };
467 
468   m_root = AllocNode();
469   m_root->m_level = 0;
470   m_unitSphereVolume = (ELEMTYPEREAL)UNIT_SPHERE_VOLUMES[NUMDIMS];
471 }
472 
473 
474 RTREE_TEMPLATE
~RTree()475 RTREE_QUAL::~RTree()
476 {
477   Reset(); // Free, or reset node memory
478 }
479 
480 
481 RTREE_TEMPLATE
Insert(const ELEMTYPE a_min[NUMDIMS],const ELEMTYPE a_max[NUMDIMS],const DATATYPE & a_dataId)482 void RTREE_QUAL::Insert(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE& a_dataId)
483 {
484 #ifdef _DEBUG
485   for(int index=0; index<NUMDIMS; ++index)
486   {
487     ASSERT(a_min[index] <= a_max[index]);
488   }
489 #endif //_DEBUG
490 
491   Rect rect;
492 
493   for(int axis=0; axis<NUMDIMS; ++axis)
494   {
495     rect.m_min[axis] = a_min[axis];
496     rect.m_max[axis] = a_max[axis];
497   }
498 
499   InsertRect(&rect, a_dataId, &m_root, 0);
500 }
501 
502 
503 RTREE_TEMPLATE
Remove(const ELEMTYPE a_min[NUMDIMS],const ELEMTYPE a_max[NUMDIMS],const DATATYPE & a_dataId)504 void RTREE_QUAL::Remove(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const DATATYPE& a_dataId)
505 {
506 #ifdef _DEBUG
507   for(int index=0; index<NUMDIMS; ++index)
508   {
509     ASSERT(a_min[index] <= a_max[index]);
510   }
511 #endif //_DEBUG
512 
513   Rect rect;
514 
515   for(int axis=0; axis<NUMDIMS; ++axis)
516   {
517     rect.m_min[axis] = a_min[axis];
518     rect.m_max[axis] = a_max[axis];
519   }
520 
521   RemoveRect(&rect, a_dataId, &m_root);
522 }
523 
524 
525 
526 RTREE_TEMPLATE
Search(const ELEMTYPE a_min[NUMDIMS],const ELEMTYPE a_max[NUMDIMS],const CONTEXT & c)527 int RTREE_QUAL::Search(const ELEMTYPE a_min[NUMDIMS], const ELEMTYPE a_max[NUMDIMS], const CONTEXT &c) const
528 {
529 #ifdef _DEBUG
530   for(int index=0; index<NUMDIMS; ++index)
531   {
532     ASSERT(a_min[index] <= a_max[index]);
533   }
534 #endif //_DEBUG
535 
536   Rect rect;
537 
538   for(int axis=0; axis<NUMDIMS; ++axis)
539   {
540     rect.m_min[axis] = a_min[axis];
541     rect.m_max[axis] = a_max[axis];
542   }
543 
544   // NOTE: May want to return search result another way, perhaps returning the number of found elements here.
545 
546   int foundCount = 0;
547   Search(m_root, &rect, foundCount, c);
548 
549   return foundCount;
550 }
551 
552 
553 
554 
555 
556 
557 RTREE_TEMPLATE
Count()558 int RTREE_QUAL::Count()
559 {
560   int count = 0;
561   CountRec(m_root, count);
562 
563   return count;
564 }
565 
566 
567 
568 RTREE_TEMPLATE
CountRec(Node * a_node,int & a_count)569 void RTREE_QUAL::CountRec(Node* a_node, int& a_count)
570 {
571   if(a_node->IsInternalNode())  // not a leaf node
572   {
573     for(int index = 0; index < a_node->m_count; ++index)
574     {
575       CountRec(a_node->m_branch[index].m_child, a_count);
576     }
577   }
578   else // A leaf node
579   {
580     a_count += a_node->m_count;
581   }
582 }
583 
584 
585 #ifdef RTREE_WANTS_IO
586 RTREE_TEMPLATE
Load(const char * a_fileName)587 bool RTREE_QUAL::Load(const char* a_fileName)
588 {
589   RemoveAll(); // Clear existing tree
590 
591   RTFileStream stream;
592   if(!stream.OpenRead(a_fileName))
593   {
594     return false;
595   }
596 
597   bool result = Load(stream);
598 
599   stream.Close();
600 
601   return result;
602 }
603 
604 
605 
606 RTREE_TEMPLATE
Load(RTFileStream & a_stream)607 bool RTREE_QUAL::Load(RTFileStream& a_stream)
608 {
609   // Write some kind of header
610   int _dataFileId = ('R'<<0)|('T'<<8)|('R'<<16)|('E'<<24);
611   int _dataSize = sizeof(DATATYPE);
612   int _dataNumDims = NUMDIMS;
613   int _dataElemSize = sizeof(ELEMTYPE);
614   int _dataElemRealSize = sizeof(ELEMTYPEREAL);
615   int _dataMaxNodes = TMAXNODES;
616   int _dataMinNodes = TMINNODES;
617 
618   int dataFileId = 0;
619   int dataSize = 0;
620   int dataNumDims = 0;
621   int dataElemSize = 0;
622   int dataElemRealSize = 0;
623   int dataMaxNodes = 0;
624   int dataMinNodes = 0;
625 
626   a_stream.Read(dataFileId);
627   a_stream.Read(dataSize);
628   a_stream.Read(dataNumDims);
629   a_stream.Read(dataElemSize);
630   a_stream.Read(dataElemRealSize);
631   a_stream.Read(dataMaxNodes);
632   a_stream.Read(dataMinNodes);
633 
634   bool result = false;
635 
636   // Test if header was valid and compatible
637   if(    (dataFileId == _dataFileId)
638       && (dataSize == _dataSize)
639       && (dataNumDims == _dataNumDims)
640       && (dataElemSize == _dataElemSize)
641       && (dataElemRealSize == _dataElemRealSize)
642       && (dataMaxNodes == _dataMaxNodes)
643       && (dataMinNodes == _dataMinNodes)
644     )
645   {
646     // Recursively load tree
647     result = LoadRec(m_root, a_stream);
648   }
649 
650   return result;
651 }
652 
653 
654 RTREE_TEMPLATE
LoadRec(Node * a_node,RTFileStream & a_stream)655 bool RTREE_QUAL::LoadRec(Node* a_node, RTFileStream& a_stream)
656 {
657   a_stream.Read(a_node->m_level);
658   a_stream.Read(a_node->m_count);
659 
660   if(a_node->IsInternalNode())  // not a leaf node
661   {
662     for(int index = 0; index < a_node->m_count; ++index)
663     {
664       Branch* curBranch = &a_node->m_branch[index];
665 
666       a_stream.ReadArray(curBranch->m_rect.m_min, NUMDIMS);
667       a_stream.ReadArray(curBranch->m_rect.m_max, NUMDIMS);
668 
669       curBranch->m_child = AllocNode();
670       LoadRec(curBranch->m_child, a_stream);
671     }
672   }
673   else // A leaf node
674   {
675     for(int index = 0; index < a_node->m_count; ++index)
676     {
677       Branch* curBranch = &a_node->m_branch[index];
678 
679       a_stream.ReadArray(curBranch->m_rect.m_min, NUMDIMS);
680       a_stream.ReadArray(curBranch->m_rect.m_max, NUMDIMS);
681 
682       a_stream.Read(curBranch->m_data);
683     }
684   }
685 
686   return true; // Should do more error checking on I/O operations
687 }
688 
689 
690 RTREE_TEMPLATE
Save(const char * a_fileName)691 bool RTREE_QUAL::Save(const char* a_fileName)
692 {
693   RTFileStream stream;
694   if(!stream.OpenWrite(a_fileName))
695   {
696     return false;
697   }
698 
699   bool result = Save(stream);
700 
701   stream.Close();
702 
703   return result;
704 }
705 
706 
707 RTREE_TEMPLATE
Save(RTFileStream & a_stream)708 bool RTREE_QUAL::Save(RTFileStream& a_stream)
709 {
710   // Write some kind of header
711   int dataFileId = ('R'<<0)|('T'<<8)|('R'<<16)|('E'<<24);
712   int dataSize = sizeof(DATATYPE);
713   int dataNumDims = NUMDIMS;
714   int dataElemSize = sizeof(ELEMTYPE);
715   int dataElemRealSize = sizeof(ELEMTYPEREAL);
716   int dataMaxNodes = TMAXNODES;
717   int dataMinNodes = TMINNODES;
718 
719   a_stream.Write(dataFileId);
720   a_stream.Write(dataSize);
721   a_stream.Write(dataNumDims);
722   a_stream.Write(dataElemSize);
723   a_stream.Write(dataElemRealSize);
724   a_stream.Write(dataMaxNodes);
725   a_stream.Write(dataMinNodes);
726 
727   // Recursively save tree
728   bool result = SaveRec(m_root, a_stream);
729 
730   return result;
731 }
732 
733 
734 RTREE_TEMPLATE
SaveRec(Node * a_node,RTFileStream & a_stream)735 bool RTREE_QUAL::SaveRec(Node* a_node, RTFileStream& a_stream)
736 {
737   a_stream.Write(a_node->m_level);
738   a_stream.Write(a_node->m_count);
739 
740   if(a_node->IsInternalNode())  // not a leaf node
741   {
742     for(int index = 0; index < a_node->m_count; ++index)
743     {
744       Branch* curBranch = &a_node->m_branch[index];
745 
746       a_stream.WriteArray(curBranch->m_rect.m_min, NUMDIMS);
747       a_stream.WriteArray(curBranch->m_rect.m_max, NUMDIMS);
748 
749       SaveRec(curBranch->m_child, a_stream);
750     }
751   }
752   else // A leaf node
753   {
754     for(int index = 0; index < a_node->m_count; ++index)
755     {
756       Branch* curBranch = &a_node->m_branch[index];
757 
758       a_stream.WriteArray(curBranch->m_rect.m_min, NUMDIMS);
759       a_stream.WriteArray(curBranch->m_rect.m_max, NUMDIMS);
760 
761       a_stream.Write(curBranch->m_data);
762     }
763   }
764 
765   return true; // Should do more error checking on I/O operations
766 }
767 #endif
768 
769 
770 RTREE_TEMPLATE
RemoveAll()771 void RTREE_QUAL::RemoveAll()
772 {
773   // Delete all existing nodes
774   Reset();
775 
776   m_root = AllocNode();
777   m_root->m_level = 0;
778 }
779 
780 
781 RTREE_TEMPLATE
Reset()782 void RTREE_QUAL::Reset()
783 {
784 #ifdef RTREE_DONT_USE_MEMPOOLS
785   // Delete all existing nodes
786   RemoveAllRec(m_root);
787 #else // RTREE_DONT_USE_MEMPOOLS
788   // Just reset memory pools.  We are not using complex types
789   // EXAMPLE
790 #endif // RTREE_DONT_USE_MEMPOOLS
791 }
792 
793 
794 RTREE_TEMPLATE
RemoveAllRec(Node * a_node)795 void RTREE_QUAL::RemoveAllRec(Node* a_node)
796 {
797   ASSERT(a_node);
798   ASSERT(a_node->m_level >= 0);
799 
800   if(a_node->IsInternalNode()) // This is an internal node in the tree
801   {
802     for(int index=0; index < a_node->m_count; ++index)
803     {
804       RemoveAllRec(a_node->m_branch[index].m_child);
805     }
806   }
807   FreeNode(a_node);
808 }
809 
810 
811 RTREE_TEMPLATE
AllocNode()812 typename RTREE_QUAL::Node* RTREE_QUAL::AllocNode()
813 {
814   Node* newNode;
815 #ifdef RTREE_DONT_USE_MEMPOOLS
816   newNode = new Node;
817 #else // RTREE_DONT_USE_MEMPOOLS
818   // EXAMPLE
819 #endif // RTREE_DONT_USE_MEMPOOLS
820   InitNode(newNode);
821   return newNode;
822 }
823 
824 
825 RTREE_TEMPLATE
FreeNode(Node * a_node)826 void RTREE_QUAL::FreeNode(Node* a_node)
827 {
828   ASSERT(a_node);
829 
830 #ifdef RTREE_DONT_USE_MEMPOOLS
831   delete a_node;
832 #else // RTREE_DONT_USE_MEMPOOLS
833   // EXAMPLE
834 #endif // RTREE_DONT_USE_MEMPOOLS
835 }
836 
837 
838 // Allocate space for a node in the list used in DeletRect to
839 // store Nodes that are too empty.
840 RTREE_TEMPLATE
AllocListNode()841 typename RTREE_QUAL::ListNode* RTREE_QUAL::AllocListNode()
842 {
843 #ifdef RTREE_DONT_USE_MEMPOOLS
844   return new ListNode;
845 #else // RTREE_DONT_USE_MEMPOOLS
846   // EXAMPLE
847 #endif // RTREE_DONT_USE_MEMPOOLS
848 }
849 
850 
851 RTREE_TEMPLATE
FreeListNode(ListNode * a_listNode)852 void RTREE_QUAL::FreeListNode(ListNode* a_listNode)
853 {
854 #ifdef RTREE_DONT_USE_MEMPOOLS
855   delete a_listNode;
856 #else // RTREE_DONT_USE_MEMPOOLS
857   // EXAMPLE
858 #endif // RTREE_DONT_USE_MEMPOOLS
859 }
860 
861 
862 RTREE_TEMPLATE
InitNode(Node * a_node)863 void RTREE_QUAL::InitNode(Node* a_node)
864 {
865   a_node->m_count = 0;
866   a_node->m_level = -1;
867 }
868 
869 
870 RTREE_TEMPLATE
InitRect(Rect * a_rect)871 void RTREE_QUAL::InitRect(Rect* a_rect)
872 {
873   for(int index = 0; index < NUMDIMS; ++index)
874   {
875     a_rect->m_min[index] = (ELEMTYPE)0;
876     a_rect->m_max[index] = (ELEMTYPE)0;
877   }
878 }
879 
880 
881 // Inserts a new data rectangle into the index structure.
882 // Recursively descends tree, propagates splits back up.
883 // Returns 0 if node was not split.  Old node updated.
884 // If node was split, returns 1 and sets the pointer pointed to by
885 // new_node to point to the new node.  Old node updated to become one of two.
886 // The level argument specifies the number of steps up from the leaf
887 // level to insert; e.g. a data rectangle goes in at level = 0.
888 RTREE_TEMPLATE
InsertRectRec(Rect * a_rect,const DATATYPE & a_id,Node * a_node,Node ** a_newNode,int a_level)889 bool RTREE_QUAL::InsertRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node, Node** a_newNode, int a_level)
890 {
891   ASSERT(a_rect && a_node && a_newNode);
892   ASSERT(a_level >= 0 && a_level <= a_node->m_level);
893 
894   int index;
895   Branch branch;
896   Node* otherNode;
897 
898   // Still above level for insertion, go down tree recursively
899   if(a_node->m_level > a_level)
900   {
901     index = PickBranch(a_rect, a_node);
902     if (!InsertRectRec(a_rect, a_id, a_node->m_branch[index].m_child, &otherNode, a_level))
903     {
904       // Child was not split
905       a_node->m_branch[index].m_rect = CombineRect(a_rect, &(a_node->m_branch[index].m_rect));
906       return false;
907     }
908     else // Child was split
909     {
910       a_node->m_branch[index].m_rect = NodeCover(a_node->m_branch[index].m_child);
911       branch.m_child = otherNode;
912       branch.m_rect = NodeCover(otherNode);
913       return AddBranch(&branch, a_node, a_newNode);
914     }
915   }
916   else if(a_node->m_level == a_level) // Have reached level for insertion. Add rect, split if necessary
917   {
918     branch.m_rect = *a_rect;
919     branch.m_child = (Node*) a_id;
920     // Child field of leaves contains id of data record
921     return AddBranch(&branch, a_node, a_newNode);
922   }
923   else
924   {
925     // Should never occur
926     ASSERT(0);
927     return false;
928   }
929 }
930 
931 
932 // Insert a data rectangle into an index structure.
933 // InsertRect provides for splitting the root;
934 // returns 1 if root was split, 0 if it was not.
935 // The level argument specifies the number of steps up from the leaf
936 // level to insert; e.g. a data rectangle goes in at level = 0.
937 // InsertRect2 does the recursion.
938 //
939 RTREE_TEMPLATE
InsertRect(Rect * a_rect,const DATATYPE & a_id,Node ** a_root,int a_level)940 bool RTREE_QUAL::InsertRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root, int a_level)
941 {
942   ASSERT(a_rect && a_root);
943   ASSERT(a_level >= 0 && a_level <= (*a_root)->m_level);
944 #ifdef _DEBUG
945   for(int index=0; index < NUMDIMS; ++index)
946   {
947     ASSERT(a_rect->m_min[index] <= a_rect->m_max[index]);
948   }
949 #endif //_DEBUG
950 
951   Node* newRoot;
952   Node* newNode;
953   Branch branch;
954 
955   if(InsertRectRec(a_rect, a_id, *a_root, &newNode, a_level))  // Root split
956   {
957     newRoot = AllocNode();  // Grow tree taller and new root
958     newRoot->m_level = (*a_root)->m_level + 1;
959     branch.m_rect = NodeCover(*a_root);
960     branch.m_child = *a_root;
961     AddBranch(&branch, newRoot, NULL);
962     branch.m_rect = NodeCover(newNode);
963     branch.m_child = newNode;
964     AddBranch(&branch, newRoot, NULL);
965     *a_root = newRoot;
966     return true;
967   }
968 
969   return false;
970 }
971 
972 
973 // Find the smallest rectangle that includes all rectangles in branches of a node.
974 RTREE_TEMPLATE
NodeCover(Node * a_node)975 typename RTREE_QUAL::Rect RTREE_QUAL::NodeCover(Node* a_node)
976 {
977   ASSERT(a_node);
978 
979   int firstTime = true;
980   Rect rect;
981   InitRect(&rect);
982 
983   for(int index = 0; index < a_node->m_count; ++index)
984   {
985     if(firstTime)
986     {
987       rect = a_node->m_branch[index].m_rect;
988       firstTime = false;
989     }
990     else
991     {
992       rect = CombineRect(&rect, &(a_node->m_branch[index].m_rect));
993     }
994   }
995 
996   return rect;
997 }
998 
999 
1000 // Add a branch to a node.  Split the node if necessary.
1001 // Returns 0 if node not split.  Old node updated.
1002 // Returns 1 if node split, sets *new_node to address of new node.
1003 // Old node updated, becomes one of two.
1004 RTREE_TEMPLATE
AddBranch(Branch * a_branch,Node * a_node,Node ** a_newNode)1005 bool RTREE_QUAL::AddBranch(Branch* a_branch, Node* a_node, Node** a_newNode)
1006 {
1007   ASSERT(a_branch);
1008   ASSERT(a_node);
1009 
1010   if(a_node->m_count < MAXNODES)  // Split won't be necessary
1011   {
1012     a_node->m_branch[a_node->m_count] = *a_branch;
1013     ++a_node->m_count;
1014 
1015     return false;
1016   }
1017   else
1018   {
1019     ASSERT(a_newNode);
1020 
1021     SplitNode(a_node, a_branch, a_newNode);
1022     return true;
1023   }
1024 }
1025 
1026 
1027 // Disconnect a dependent node.
1028 // Caller must return (or stop using iteration index) after this as count has changed
1029 RTREE_TEMPLATE
DisconnectBranch(Node * a_node,int a_index)1030 void RTREE_QUAL::DisconnectBranch(Node* a_node, int a_index)
1031 {
1032   ASSERT(a_node && (a_index >= 0) && (a_index < MAXNODES));
1033   ASSERT(a_node->m_count > 0);
1034 
1035   // Remove element by swapping with the last element to prevent gaps in array
1036   a_node->m_branch[a_index] = a_node->m_branch[a_node->m_count - 1];
1037 
1038   --a_node->m_count;
1039 }
1040 
1041 
1042 // Pick a branch.  Pick the one that will need the smallest increase
1043 // in area to accomodate the new rectangle.  This will result in the
1044 // least total area for the covering rectangles in the current node.
1045 // In case of a tie, pick the one which was smaller before, to get
1046 // the best resolution when searching.
1047 RTREE_TEMPLATE
PickBranch(Rect * a_rect,Node * a_node)1048 int RTREE_QUAL::PickBranch(Rect* a_rect, Node* a_node)
1049 {
1050   ASSERT(a_rect && a_node);
1051 
1052   bool firstTime = true;
1053   ELEMTYPEREAL increase;
1054   ELEMTYPEREAL bestIncr = (ELEMTYPEREAL)-1;
1055   ELEMTYPEREAL area;
1056   ELEMTYPEREAL bestArea = 0;
1057   int best = 0;
1058   Rect tempRect;
1059 
1060   for(int index=0; index < a_node->m_count; ++index)
1061   {
1062     Rect* curRect = &a_node->m_branch[index].m_rect;
1063     area = CalcRectVolume(curRect);
1064     tempRect = CombineRect(a_rect, curRect);
1065     increase = CalcRectVolume(&tempRect) - area;
1066     if((increase < bestIncr) || firstTime)
1067     {
1068       best = index;
1069       bestArea = area;
1070       bestIncr = increase;
1071       firstTime = false;
1072     }
1073     else if((increase == bestIncr) && (area < bestArea))
1074     {
1075       best = index;
1076       bestArea = area;
1077       bestIncr = increase;
1078     }
1079   }
1080   return best;
1081 }
1082 
1083 
1084 // Combine two rectangles into larger one containing both
1085 RTREE_TEMPLATE
CombineRect(Rect * a_rectA,Rect * a_rectB)1086 typename RTREE_QUAL::Rect RTREE_QUAL::CombineRect(Rect* a_rectA, Rect* a_rectB)
1087 {
1088   ASSERT(a_rectA && a_rectB);
1089 
1090   Rect newRect;
1091 
1092   for(int index = 0; index < NUMDIMS; ++index)
1093   {
1094     newRect.m_min[index] = rtree_min(a_rectA->m_min[index], a_rectB->m_min[index]);
1095     newRect.m_max[index] = rtree_max(a_rectA->m_max[index], a_rectB->m_max[index]);
1096   }
1097 
1098   return newRect;
1099 }
1100 
1101 
1102 
1103 // Split a node.
1104 // Divides the nodes branches and the extra one between two nodes.
1105 // Old node is one of the new ones, and one really new one is created.
1106 // Tries more than one method for choosing a partition, uses best result.
1107 RTREE_TEMPLATE
SplitNode(Node * a_node,Branch * a_branch,Node ** a_newNode)1108 void RTREE_QUAL::SplitNode(Node* a_node, Branch* a_branch, Node** a_newNode)
1109 {
1110   ASSERT(a_node);
1111   ASSERT(a_branch);
1112 
1113   // Could just use local here, but member or external is faster since it is reused
1114   PartitionVars localVars;
1115   PartitionVars* parVars = &localVars;
1116   int level;
1117 
1118   // Load all the branches into a buffer, initialize old node
1119   level = a_node->m_level;
1120   GetBranches(a_node, a_branch, parVars);
1121 
1122   // Find partition
1123   ChoosePartition(parVars, MINNODES);
1124 
1125   // Put branches from buffer into 2 nodes according to chosen partition
1126   *a_newNode = AllocNode();
1127   (*a_newNode)->m_level = a_node->m_level = level;
1128   LoadNodes(a_node, *a_newNode, parVars);
1129 
1130   ASSERT((a_node->m_count + (*a_newNode)->m_count) == parVars->m_total);
1131 }
1132 
1133 
1134 // Calculate the n-dimensional volume of a rectangle
1135 RTREE_TEMPLATE
RectVolume(Rect * a_rect)1136 ELEMTYPEREAL RTREE_QUAL::RectVolume(Rect* a_rect)
1137 {
1138   ASSERT(a_rect);
1139 
1140   ELEMTYPEREAL volume = (ELEMTYPEREAL)1;
1141 
1142   for(int index=0; index<NUMDIMS; ++index)
1143   {
1144     volume *= a_rect->m_max[index] - a_rect->m_min[index];
1145   }
1146 
1147   ASSERT(volume >= (ELEMTYPEREAL)0);
1148 
1149   return volume;
1150 }
1151 
1152 
1153 // The exact volume of the bounding sphere for the given Rect
1154 RTREE_TEMPLATE
RectSphericalVolume(Rect * a_rect)1155 ELEMTYPEREAL RTREE_QUAL::RectSphericalVolume(Rect* a_rect)
1156 {
1157   ASSERT(a_rect);
1158 
1159   ELEMTYPEREAL sumOfSquares = (ELEMTYPEREAL)0;
1160   ELEMTYPEREAL radius;
1161 
1162   for(int index=0; index < NUMDIMS; ++index)
1163   {
1164     ELEMTYPEREAL halfExtent = ((ELEMTYPEREAL)a_rect->m_max[index] - (ELEMTYPEREAL)a_rect->m_min[index]) * 0.5f;
1165     sumOfSquares += halfExtent * halfExtent;
1166   }
1167 
1168   radius = (ELEMTYPEREAL)sqrt(sumOfSquares);
1169 
1170   // Pow maybe slow, so test for common dims like 2,3 and just use x*x, x*x*x.
1171   if(NUMDIMS == 3)
1172   {
1173     return (radius * radius * radius * m_unitSphereVolume);
1174   }
1175   else if(NUMDIMS == 2)
1176   {
1177     return (radius * radius * m_unitSphereVolume);
1178   }
1179   else
1180   {
1181     return (ELEMTYPEREAL)(pow(radius, NUMDIMS) * m_unitSphereVolume);
1182   }
1183 }
1184 
1185 
1186 // Use one of the methods to calculate retangle volume
1187 RTREE_TEMPLATE
CalcRectVolume(Rect * a_rect)1188 ELEMTYPEREAL RTREE_QUAL::CalcRectVolume(Rect* a_rect)
1189 {
1190 #ifdef RTREE_USE_SPHERICAL_VOLUME
1191   return RectSphericalVolume(a_rect); // Slower but helps certain merge cases
1192 #else // RTREE_USE_SPHERICAL_VOLUME
1193   return RectVolume(a_rect); // Faster but can cause poor merges
1194 #endif // RTREE_USE_SPHERICAL_VOLUME
1195 }
1196 
1197 
1198 // Load branch buffer with branches from full node plus the extra branch.
1199 RTREE_TEMPLATE
GetBranches(Node * a_node,Branch * a_branch,PartitionVars * a_parVars)1200 void RTREE_QUAL::GetBranches(Node* a_node, Branch* a_branch, PartitionVars* a_parVars)
1201 {
1202   ASSERT(a_node);
1203   ASSERT(a_branch);
1204 
1205   ASSERT(a_node->m_count == MAXNODES);
1206 
1207   // Load the branch buffer
1208   int index;
1209   for(index=0; index < MAXNODES; ++index)
1210   {
1211     a_parVars->m_branchBuf[index] = a_node->m_branch[index];
1212   }
1213   a_parVars->m_branchBuf[MAXNODES] = *a_branch;
1214   a_parVars->m_branchCount = MAXNODES + 1;
1215 
1216   // Calculate rect containing all in the set
1217   a_parVars->m_coverSplit = a_parVars->m_branchBuf[0].m_rect;
1218   for(index=1; index < MAXNODES+1; ++index)
1219   {
1220     a_parVars->m_coverSplit = CombineRect(&a_parVars->m_coverSplit, &a_parVars->m_branchBuf[index].m_rect);
1221   }
1222   a_parVars->m_coverSplitArea = CalcRectVolume(&a_parVars->m_coverSplit);
1223 
1224   InitNode(a_node);
1225 }
1226 
1227 
1228 // Method #0 for choosing a partition:
1229 // As the seeds for the two groups, pick the two rects that would waste the
1230 // most area if covered by a single rectangle, i.e. evidently the worst pair
1231 // to have in the same group.
1232 // Of the remaining, one at a time is chosen to be put in one of the two groups.
1233 // The one chosen is the one with the greatest difference in area expansion
1234 // depending on which group - the rect most strongly attracted to one group
1235 // and repelled from the other.
1236 // If one group gets too full (more would force other group to violate min
1237 // fill requirement) then other group gets the rest.
1238 // These last are the ones that can go in either group most easily.
1239 RTREE_TEMPLATE
ChoosePartition(PartitionVars * a_parVars,int a_minFill)1240 void RTREE_QUAL::ChoosePartition(PartitionVars* a_parVars, int a_minFill)
1241 {
1242   ASSERT(a_parVars);
1243 
1244   ELEMTYPEREAL biggestDiff;
1245   int group, chosen = 0, betterGroup = 0;
1246 
1247   InitParVars(a_parVars, a_parVars->m_branchCount, a_minFill);
1248   PickSeeds(a_parVars);
1249 
1250   while (((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total)
1251        && (a_parVars->m_count[0] < (a_parVars->m_total - a_parVars->m_minFill))
1252        && (a_parVars->m_count[1] < (a_parVars->m_total - a_parVars->m_minFill)))
1253   {
1254     biggestDiff = (ELEMTYPEREAL) -1;
1255     for(int index=0; index<a_parVars->m_total; ++index)
1256     {
1257       if(!a_parVars->m_taken[index])
1258       {
1259         Rect* curRect = &a_parVars->m_branchBuf[index].m_rect;
1260         Rect rect0 = CombineRect(curRect, &a_parVars->m_cover[0]);
1261         Rect rect1 = CombineRect(curRect, &a_parVars->m_cover[1]);
1262         ELEMTYPEREAL growth0 = CalcRectVolume(&rect0) - a_parVars->m_area[0];
1263         ELEMTYPEREAL growth1 = CalcRectVolume(&rect1) - a_parVars->m_area[1];
1264         ELEMTYPEREAL diff = growth1 - growth0;
1265         if(diff >= 0)
1266         {
1267           group = 0;
1268         }
1269         else
1270         {
1271           group = 1;
1272           diff = -diff;
1273         }
1274 
1275         if(diff > biggestDiff)
1276         {
1277           biggestDiff = diff;
1278           chosen = index;
1279           betterGroup = group;
1280         }
1281         else if((diff == biggestDiff) && (a_parVars->m_count[group] < a_parVars->m_count[betterGroup]))
1282         {
1283           chosen = index;
1284           betterGroup = group;
1285         }
1286       }
1287     }
1288     Classify(chosen, betterGroup, a_parVars);
1289   }
1290 
1291   // If one group too full, put remaining rects in the other
1292   if((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total)
1293   {
1294     if(a_parVars->m_count[0] >= a_parVars->m_total - a_parVars->m_minFill)
1295     {
1296       group = 1;
1297     }
1298     else
1299     {
1300       group = 0;
1301     }
1302     for(int index=0; index<a_parVars->m_total; ++index)
1303     {
1304       if(!a_parVars->m_taken[index])
1305       {
1306         Classify(index, group, a_parVars);
1307       }
1308     }
1309   }
1310 
1311   ASSERT((a_parVars->m_count[0] + a_parVars->m_count[1]) == a_parVars->m_total);
1312   ASSERT((a_parVars->m_count[0] >= a_parVars->m_minFill) &&
1313         (a_parVars->m_count[1] >= a_parVars->m_minFill));
1314 }
1315 
1316 
1317 // Copy branches from the buffer into two nodes according to the partition.
1318 RTREE_TEMPLATE
LoadNodes(Node * a_nodeA,Node * a_nodeB,PartitionVars * a_parVars)1319 void RTREE_QUAL::LoadNodes(Node* a_nodeA, Node* a_nodeB, PartitionVars* a_parVars)
1320 {
1321   ASSERT(a_nodeA);
1322   ASSERT(a_nodeB);
1323   ASSERT(a_parVars);
1324 
1325   for(int index=0; index < a_parVars->m_total; ++index)
1326   {
1327     ASSERT(a_parVars->m_partition[index] == 0 || a_parVars->m_partition[index] == 1);
1328 
1329     if(a_parVars->m_partition[index] == 0)
1330     {
1331       AddBranch(&a_parVars->m_branchBuf[index], a_nodeA, NULL);
1332     }
1333     else if(a_parVars->m_partition[index] == 1)
1334     {
1335       AddBranch(&a_parVars->m_branchBuf[index], a_nodeB, NULL);
1336     }
1337   }
1338 }
1339 
1340 
1341 // Initialize a PartitionVars structure.
1342 RTREE_TEMPLATE
InitParVars(PartitionVars * a_parVars,int a_maxRects,int a_minFill)1343 void RTREE_QUAL::InitParVars(PartitionVars* a_parVars, int a_maxRects, int a_minFill)
1344 {
1345   ASSERT(a_parVars);
1346 
1347   a_parVars->m_count[0] = a_parVars->m_count[1] = 0;
1348   a_parVars->m_area[0] = a_parVars->m_area[1] = (ELEMTYPEREAL)0;
1349   a_parVars->m_total = a_maxRects;
1350   a_parVars->m_minFill = a_minFill;
1351   for(int index=0; index < a_maxRects; ++index)
1352   {
1353     a_parVars->m_taken[index] = false;
1354     a_parVars->m_partition[index] = -1;
1355   }
1356 }
1357 
1358 
1359 RTREE_TEMPLATE
PickSeeds(PartitionVars * a_parVars)1360 void RTREE_QUAL::PickSeeds(PartitionVars* a_parVars)
1361 {
1362   int seed0=0, seed1=1;
1363   ELEMTYPEREAL worst, waste;
1364   ELEMTYPEREAL area[MAXNODES+1];
1365 
1366   for(int index=0; index<a_parVars->m_total; ++index)
1367   {
1368     area[index] = CalcRectVolume(&a_parVars->m_branchBuf[index].m_rect);
1369   }
1370 
1371   worst = -a_parVars->m_coverSplitArea - 1;
1372   for(int indexA=0; indexA < a_parVars->m_total-1; ++indexA)
1373   {
1374     for(int indexB = indexA+1; indexB < a_parVars->m_total; ++indexB)
1375     {
1376       Rect oneRect = CombineRect(&a_parVars->m_branchBuf[indexA].m_rect, &a_parVars->m_branchBuf[indexB].m_rect);
1377       waste = CalcRectVolume(&oneRect) - area[indexA] - area[indexB];
1378       if(waste > worst)
1379       {
1380         worst = waste;
1381         seed0 = indexA;
1382         seed1 = indexB;
1383       }
1384     }
1385   }
1386   Classify(seed0, 0, a_parVars);
1387   Classify(seed1, 1, a_parVars);
1388 }
1389 
1390 
1391 // Put a branch in one of the groups.
1392 RTREE_TEMPLATE
Classify(int a_index,int a_group,PartitionVars * a_parVars)1393 void RTREE_QUAL::Classify(int a_index, int a_group, PartitionVars* a_parVars)
1394 {
1395   ASSERT(a_parVars);
1396   ASSERT(!a_parVars->m_taken[a_index]);
1397 
1398   a_parVars->m_partition[a_index] = a_group;
1399   a_parVars->m_taken[a_index] = true;
1400 
1401   if (a_parVars->m_count[a_group] == 0)
1402   {
1403     a_parVars->m_cover[a_group] = a_parVars->m_branchBuf[a_index].m_rect;
1404   }
1405   else
1406   {
1407     a_parVars->m_cover[a_group] = CombineRect(&a_parVars->m_branchBuf[a_index].m_rect, &a_parVars->m_cover[a_group]);
1408   }
1409   a_parVars->m_area[a_group] = CalcRectVolume(&a_parVars->m_cover[a_group]);
1410   ++a_parVars->m_count[a_group];
1411 }
1412 
1413 
1414 // Delete a data rectangle from an index structure.
1415 // Pass in a pointer to a Rect, the tid of the record, ptr to ptr to root node.
1416 // Returns 1 if record not found, 0 if success.
1417 // RemoveRect provides for eliminating the root.
1418 RTREE_TEMPLATE
RemoveRect(Rect * a_rect,const DATATYPE & a_id,Node ** a_root)1419 bool RTREE_QUAL::RemoveRect(Rect* a_rect, const DATATYPE& a_id, Node** a_root)
1420 {
1421   ASSERT(a_rect && a_root);
1422   ASSERT(*a_root);
1423 
1424   Node* tempNode;
1425   ListNode* reInsertList = NULL;
1426 
1427   if(!RemoveRectRec(a_rect, a_id, *a_root, &reInsertList))
1428   {
1429     // Found and deleted a data item
1430     // Reinsert any branches from eliminated nodes
1431     while(reInsertList)
1432     {
1433       tempNode = reInsertList->m_node;
1434 
1435       for(int index = 0; index < tempNode->m_count; ++index)
1436       {
1437         InsertRect(&(tempNode->m_branch[index].m_rect),
1438                    tempNode->m_branch[index].m_data,
1439                    a_root,
1440                    tempNode->m_level);
1441       }
1442 
1443       ListNode* remLNode = reInsertList;
1444       reInsertList = reInsertList->m_next;
1445 
1446       FreeNode(remLNode->m_node);
1447       FreeListNode(remLNode);
1448     }
1449 
1450     // Check for redundant root (not leaf, 1 child) and eliminate
1451     if((*a_root)->m_count == 1 && (*a_root)->IsInternalNode())
1452     {
1453       tempNode = (*a_root)->m_branch[0].m_child;
1454 
1455       ASSERT(tempNode);
1456       FreeNode(*a_root);
1457       *a_root = tempNode;
1458     }
1459     return false;
1460   }
1461   else
1462   {
1463     return true;
1464   }
1465 }
1466 
1467 
1468 // Delete a rectangle from non-root part of an index structure.
1469 // Called by RemoveRect.  Descends tree recursively,
1470 // merges branches on the way back up.
1471 // Returns 1 if record not found, 0 if success.
1472 RTREE_TEMPLATE
RemoveRectRec(Rect * a_rect,const DATATYPE & a_id,Node * a_node,ListNode ** a_listNode)1473 bool RTREE_QUAL::RemoveRectRec(Rect* a_rect, const DATATYPE& a_id, Node* a_node, ListNode** a_listNode)
1474 {
1475   ASSERT(a_rect && a_node && a_listNode);
1476   ASSERT(a_node->m_level >= 0);
1477 
1478   if(a_node->IsInternalNode())  // not a leaf node
1479   {
1480     for(int index = 0; index < a_node->m_count; ++index)
1481     {
1482       if(Overlap(a_rect, &(a_node->m_branch[index].m_rect)))
1483       {
1484         if(!RemoveRectRec(a_rect, a_id, a_node->m_branch[index].m_child, a_listNode))
1485         {
1486           if(a_node->m_branch[index].m_child->m_count >= MINNODES)
1487           {
1488             // child removed, just resize parent rect
1489             a_node->m_branch[index].m_rect = NodeCover(a_node->m_branch[index].m_child);
1490           }
1491           else
1492           {
1493             // child removed, not enough entries in node, eliminate node
1494             ReInsert(a_node->m_branch[index].m_child, a_listNode);
1495             DisconnectBranch(a_node, index); // Must return after this call as count has changed
1496           }
1497           return false;
1498         }
1499       }
1500     }
1501     return true;
1502   }
1503   else // A leaf node
1504   {
1505     for(int index = 0; index < a_node->m_count; ++index)
1506     {
1507       if(a_node->m_branch[index].m_child == (Node*)a_id)
1508       {
1509         DisconnectBranch(a_node, index); // Must return after this call as count has changed
1510         return false;
1511       }
1512     }
1513     return true;
1514   }
1515 }
1516 
1517 
1518 // Decide whether two rectangles overlap.
1519 RTREE_TEMPLATE
Overlap(Rect * a_rectA,Rect * a_rectB)1520 bool RTREE_QUAL::Overlap(Rect* a_rectA, Rect* a_rectB) const
1521 {
1522   ASSERT(a_rectA && a_rectB);
1523 
1524   for(int index=0; index < NUMDIMS; ++index)
1525   {
1526     if (a_rectA->m_min[index] > a_rectB->m_max[index] ||
1527         a_rectB->m_min[index] > a_rectA->m_max[index])
1528     {
1529       return false;
1530     }
1531   }
1532   return true;
1533 }
1534 
1535 
1536 // Add a node to the reinsertion list.  All its branches will later
1537 // be reinserted into the index structure.
1538 RTREE_TEMPLATE
ReInsert(Node * a_node,ListNode ** a_listNode)1539 void RTREE_QUAL::ReInsert(Node* a_node, ListNode** a_listNode)
1540 {
1541   ListNode* newListNode;
1542 
1543   newListNode = AllocListNode();
1544   newListNode->m_node = a_node;
1545   newListNode->m_next = *a_listNode;
1546   *a_listNode = newListNode;
1547 }
1548 
1549 
1550 
1551 // Search in an index tree or subtree for all data retangles that overlap the argument rectangle.
1552 RTREE_TEMPLATE
Search(Node * a_node,Rect * a_rect,int & a_foundCount,const CONTEXT & c)1553 bool RTREE_QUAL::Search(Node* a_node, Rect* a_rect, int& a_foundCount, const CONTEXT &c) const
1554 {
1555   ASSERT(a_node);
1556   ASSERT(a_node->m_level >= 0);
1557   ASSERT(a_rect);
1558 
1559   if(a_node->IsInternalNode()) // This is an internal node in the tree
1560   {
1561     for(int index=0; index < a_node->m_count; ++index)
1562     {
1563       if(Overlap(a_rect, &a_node->m_branch[index].m_rect))
1564       {
1565         if(!Search(a_node->m_branch[index].m_child, a_rect, a_foundCount, c))
1566         {
1567           return false; // Don't continue searching
1568         }
1569       }
1570     }
1571   }
1572   else // This is a leaf node
1573   {
1574     for(int index=0; index < a_node->m_count; ++index)
1575     {
1576       if(Overlap(a_rect, &a_node->m_branch[index].m_rect))
1577       {
1578         DATATYPE& id = a_node->m_branch[index].m_data;
1579 
1580         // NOTE: There are different ways to return results.  Here's where to modify
1581 /*        if(&a_resultCallback)
1582         {*/
1583           ++a_foundCount;
1584           (id->*myOperation)(c);
1585           /*
1586           if(!a_resultCallback(id, a_context))
1587           {
1588             return false; // Don't continue searching
1589           }
1590           */
1591 //        }
1592       }
1593     }
1594   }
1595 
1596   return true; // Continue searching
1597 }
1598 
1599 
1600 
1601 
1602 #undef RTREE_TEMPLATE
1603 #undef RTREE_QUAL
1604 
1605 #endif //RTREE_H
1606 
1607 
1608 
1609