opennurbs_rtree.cpp
Go to the documentation of this file.
00001 /* $NoKeywords: $ */
00002 /*
00003 //
00004 // Copyright (c) 1993-2012 Robert McNeel & Associates. All rights reserved.
00005 // OpenNURBS, Rhinoceros, and Rhino3D are registered trademarks of Robert
00006 // McNeel & Associates.
00007 //
00008 // THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY.
00009 // ALL IMPLIED WARRANTIES OF FITNESS FOR ANY PARTICULAR PURPOSE AND OF
00010 // MERCHANTABILITY ARE HEREBY DISCLAIMED.
00011 //                              
00012 // For complete openNURBS copyright information see <http://www.opennurbs.org>.
00013 //
00015 */
00016 
00017 #include "pcl/surface/3rdparty/opennurbs/opennurbs.h"
00018 
00019 // Dimension of tree bounding boxes
00020 #define ON_RTree_NODE_DIM 3
00021 
00023 struct ON_RTreeListNode
00024 {
00025   struct ON_RTreeListNode* m_next;      
00026   struct ON_RTreeNode* m_node;  
00027 };
00028 
00030 struct ON_RTreePartitionVars
00031 {
00032   int m_partition[ON_RTree_MAX_NODE_COUNT+1];
00033   int m_total;
00034   int m_minFill;
00035   int m_taken[ON_RTree_MAX_NODE_COUNT+1];
00036   int m_count[2];
00037   ON_RTreeBBox m_cover[2];
00038   double m_area[2];
00039 
00040   ON_RTreeBranch m_branchBuf[ON_RTree_MAX_NODE_COUNT+1];
00041   int m_branchCount;
00042   ON_RTreeBBox m_coverSplit;
00043   double m_coverSplitArea;
00044 }; 
00045 
00046 typedef bool (*ON_RTreeSearchCallback)(void* a_context, ON__INT_PTR a_id );
00047 
00048 struct ON_RTreeSearchResultCallback
00049 {
00050   void* m_context;
00051   ON_RTreeSearchCallback m_resultCallback;
00052 };
00053 
00054 static void ChoosePartition(ON_RTreePartitionVars* a_parVars, int a_minFill);
00055 static void CountRec(ON_RTreeNode* a_node, int& a_count);
00056 static void PickSeeds(ON_RTreePartitionVars* a_parVars);
00057 static void InitParVars(ON_RTreePartitionVars* a_parVars, int a_maxRects, int a_minFill);
00058 static void GetBranches(ON_RTreeNode* a_node, ON_RTreeBranch* a_branch, ON_RTreePartitionVars* a_parVars);
00059 static int PickBranch(ON_RTreeBBox* a_rect, ON_RTreeNode* a_node);
00060 static void DisconnectBranch(ON_RTreeNode* a_node, int a_index);
00061 static ON_RTreeBBox NodeCover(ON_RTreeNode* a_node);
00062 static void InitRect(ON_RTreeBBox* a_rect);
00063 static ON_RTreeBBox CombineRectHelper(const ON_RTreeBBox* a_rectA, const ON_RTreeBBox* a_rectB);
00064 static double CalcRectVolumeHelper(const ON_RTreeBBox* a_rect);
00065 static bool OverlapHelper(const ON_RTreeBBox* a_rectA, const ON_RTreeBBox* a_rectB);
00066 static double DistanceToCapsuleAxisHelper(const struct ON_RTreeCapsule* a_capsule, const ON_RTreeBBox* a_rect);
00067 static void ClassifyHelper(int a_index, int a_group, struct ON_RTreePartitionVars* a_parVars);
00068 static bool SearchHelper(const ON_RTreeNode* a_node, ON_RTreeBBox* a_rect, ON_RTreeSearchResultCallback& a_result );
00069 static bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_RTreeSearchResult& a_result );
00070 static bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_SimpleArray<ON_RTreeLeaf> &a_result );
00071 static bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_SimpleArray<int> &a_result );
00072 static bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_SimpleArray<void*> &a_result );
00073 static bool SearchHelper(const ON_RTreeNode* a_node, struct ON_RTreeSphere* a_sphere, ON_RTreeSearchResultCallback& a_result );
00074 static bool SearchHelper(const ON_RTreeNode* a_node, struct ON_RTreeCapsule* a_capsule, ON_RTreeSearchResultCallback& a_result );
00075 
00077 //
00078 // ON_RTreeMemPool
00079 //
00080 
00081 size_t ON_MemoryPageSize()
00082 {
00083 #if defined(ON_COMPILER_MSC)
00084   static size_t pagesize = 0;
00085   if ( 0 == pagesize )
00086   {
00087     SYSTEM_INFO system_info;
00088     memset(&system_info,0,sizeof(system_info));
00089     ::GetSystemInfo(&system_info);
00090     pagesize = system_info.dwPageSize;
00091     if ( pagesize <= 0 )
00092       pagesize = 4096;
00093   }
00094   return pagesize;
00095 #else
00096   return 4096;
00097 #endif
00098 }
00099 
00100 
00101 static size_t SizeofBlkLink()
00102 {
00103   // Room to reserve for struct ON_RTreeMemPool::Blk.
00104   // This number should be >= sizeof(struct ON_RTreeMemPool::Blk)
00105   // and a multiple of 16 to insure proper alignment on all CPUs.
00106   return 16;
00107 }
00108 
00109 static size_t MemPoolBlkSize( size_t leaf_count )
00110 {
00111   // sizeof_blklink = number of bytes to reserve for the
00112   // struct ON_RTreeMemPool::Blk header used to keep track
00113   // of our allocations.
00114   const size_t sizeof_blklink = SizeofBlkLink();
00115 
00116   // pagesize = OS memory manager page size.  We want the
00117   // allocated blocks to be some smallish mulitples of pagesize
00118   // to the active sections of the tree will end up in CPU cache 
00119   // when the tree is being repeatedly searched.
00120   size_t pagesize = ON_MemoryPageSize();
00121   if ( pagesize <= sizeof_blklink )
00122     pagesize = 4096;
00123 
00124   size_t node_count = 32;
00125   if ( leaf_count > 0 )
00126   {
00127     // When the user can estimate the number of leaves,
00128     // node_count is an estimate of the number of nodes needed
00129     // to build the tree.  When node_count is small, we avoid
00130     // allocating too much memory for a tiny tree.  The goal
00131     // is to avoid wasting lots of memory when there are thousands
00132     // of individual trees, most with a less than six leaves.
00133     // If there are only a few trees, or most of the trees have
00134     // lots of leaves, then hard coding node_count = 32 works
00135     // just fine.
00136     if ( 5*leaf_count < 4*ON_RTree_MAX_NODE_COUNT )
00137       node_count = 3;
00138     else if ( 5*leaf_count < 4*ON_RTree_MAX_NODE_COUNT*ON_RTree_MAX_NODE_COUNT )
00139       node_count = ON_RTree_MAX_NODE_COUNT+1;
00140   }
00141   
00142   // Set "sz" to an multiple of pagesize that is big enough
00143   // to hold node_count nodes.
00144   const size_t sizeof_node = sizeof(ON_RTreeNode);
00145   size_t sz = pagesize;
00146   size_t nodes_per_blk = ( node_count < 32 )
00147                        ? node_count
00148                        : (sz-sizeof_blklink)/sizeof_node;
00149   while ( nodes_per_blk < node_count )
00150   {
00151     sz += pagesize;
00152     nodes_per_blk = (sz-sizeof_blklink)/sizeof_node;
00153   }
00154 
00155   // Some lame allocators pad each allocation request and use the extra space
00156   // to store allocation bookkeeping information.  The "+ 2*sizeof(void*) allows 
00157   // for up to two pointers of "lame" system overhead per allocation.  The goal
00158   // is to prevent the "real" allocation from being just a hair bigger that a 
00159   // multiple of pagesize.
00160   if ( sz < sizeof_blklink + nodes_per_blk*sizeof_node + 2*sizeof(void*)  ) 
00161     nodes_per_blk--; // prevent memory manager overhead from allocating another page
00162 
00163   // Return the minimum number of bytes we need for each block. An decent
00164   // OS should assign this allocation an efficient set of pages.
00165   return (sizeof_blklink + nodes_per_blk*sizeof_node);
00166 }
00167 
00168 ON_RTreeMemPool::ON_RTreeMemPool( ON_MEMORY_POOL* heap, size_t leaf_count  )
00169 : m_nodes(0)
00170 , m_list_nodes(0)
00171 , m_buffer(0)
00172 , m_buffer_capacity(0)
00173 , m_blk_list(0)
00174 , m_sizeof_blk(0)
00175 , m_heap(heap)
00176 , m_sizeof_heap(0)
00177 {
00178   m_sizeof_blk = MemPoolBlkSize(leaf_count); 
00179 }
00180 
00181 ON_RTreeMemPool::~ON_RTreeMemPool()
00182 {
00183   DeallocateAll();
00184 }
00185 
00186 void ON_RTreeMemPool::GrowBuffer()
00187 {
00188   if ( (0 == m_sizeof_blk) || (0 != m_blk_list && 0 == m_blk_list->m_next) )
00189   {
00190     // Setting m_sizeof_blk happens twice per ON_RTreeMemPool.
00191     // The first time is typically at construction and not here.
00192     // The (0 == m_sizeof_blk) check is here to support cases 
00193     // where the ON_RTreeMemPool class is initialized by a "memset(...,0)
00194     // instead of calling its constructor.  The first block can be small
00195     // if the caller passed in a leaf_count estimate.  For the second
00196     // (0 != m_blk_list && 0 == m_blk_list->m_next) and subsequent calls
00197     // to GrowBuffer(), we use the default block size.
00198     m_sizeof_blk = MemPoolBlkSize(0); 
00199   }
00200 
00201   struct Blk* blk = (struct Blk*)onmalloc_from_pool(m_heap,m_sizeof_blk);
00202   if ( blk )
00203   {
00204     m_sizeof_heap += m_sizeof_blk;
00205     blk->m_next = m_blk_list;
00206     m_blk_list = blk;
00207     const size_t sizeof_blklink = SizeofBlkLink();
00208     m_buffer = ((unsigned char*)m_blk_list) + sizeof_blklink;
00209     m_buffer_capacity = m_sizeof_blk - sizeof_blklink;
00210   }
00211   else
00212   {
00213     m_buffer = 0;
00214     m_buffer_capacity = 0;
00215     ON_ERROR("ON_RTreeMemPool::GrowBuffer - out of memory");
00216   }
00217 }
00218 
00219 ON_RTreeNode* ON_RTreeMemPool::AllocNode()
00220 {
00221   ON_RTreeNode* node = (ON_RTreeNode*)m_nodes;
00222   if ( node )
00223   {
00224     m_nodes = m_nodes->m_next;
00225   }
00226   else
00227   {
00228     const size_t node_sz = sizeof(*node);
00229     if ( m_buffer_capacity < node_sz )
00230       GrowBuffer();
00231 
00232     if ( 0 == (node = (ON_RTreeNode*)m_buffer) )
00233     {
00234       ON_ERROR("ON_RTreeMemPool::AllocNode() - out of memory");
00235       return 0;
00236     }
00237 
00238     m_buffer += node_sz;
00239     m_buffer_capacity -= node_sz;
00240   }
00241 
00242   // initialize node
00243   node->m_count = 0;
00244   node->m_level = -1;
00245 
00246   return node;
00247 }
00248 
00249 void ON_RTreeMemPool::FreeNode(ON_RTreeNode* node)
00250 {
00251   if ( node )
00252   {
00253     struct Blk* blk = (struct Blk*)node;
00254     blk->m_next = m_nodes;
00255     m_nodes = blk;
00256   }
00257 }
00258 
00259 struct ON_RTreeListNode* ON_RTreeMemPool::AllocListNode()
00260 {
00261   struct ON_RTreeListNode* list_node = (struct ON_RTreeListNode*)m_list_nodes;
00262   if ( list_node )
00263   {
00264     m_list_nodes = m_list_nodes->m_next;
00265   }
00266   else
00267   {
00268     size_t list_node_sz = sizeof(*list_node);
00269     if ( m_buffer_capacity < list_node_sz )
00270     {
00271       GrowBuffer();
00272     }
00273     list_node = (struct ON_RTreeListNode*)m_buffer;
00274     if ( list_node )
00275     {
00276       m_buffer += list_node_sz;
00277       m_buffer_capacity -= list_node_sz;
00278     }
00279   }
00280   return list_node;
00281 }
00282 
00283 void ON_RTreeMemPool::FreeListNode(struct ON_RTreeListNode* list_node)
00284 {
00285   if ( list_node )
00286   {
00287     struct Blk* blk = (struct Blk*)list_node;
00288     blk->m_next = m_list_nodes;
00289     m_list_nodes = blk;
00290   }
00291 }
00292 
00293 size_t ON_RTreeMemPool::SizeOf() const
00294 {
00295   return m_sizeof_heap;
00296 }
00297 
00298 size_t ON_RTreeMemPool::SizeOfUnusedBuffer() const
00299 {
00300   const struct Blk* blk;
00301   size_t sz = m_buffer_capacity;
00302   for ( blk = m_nodes; blk; blk = blk->m_next )
00303   {
00304     sz += sizeof(struct ON_RTreeNode);
00305   }
00306   for ( blk = m_list_nodes; blk; blk = blk->m_next )
00307   {
00308     sz += sizeof(struct ON_RTreeListNode);
00309   }
00310   return sz;
00311 }
00312 
00313 void ON_RTreeMemPool::DeallocateAll()
00314 {
00315   struct Blk* p = m_blk_list;
00316 
00317   m_nodes = 0;
00318   m_list_nodes = 0;
00319   m_buffer = 0;
00320   m_buffer_capacity = 0;
00321   m_blk_list = 0;
00322   m_sizeof_blk = 0;
00323   m_sizeof_heap = 0;
00324 
00325   while(p)
00326   {
00327     struct Blk* next = p->m_next; 
00328     onfree(p);
00329     p = next; 
00330   }
00331 }
00332 
00334 //
00335 // ON_RTreeIterator
00336 //
00337 
00338 ON_RTreeIterator::ON_RTreeIterator()
00339 { 
00340   Initialize(0);
00341 }
00342 
00343 ON_RTreeIterator::ON_RTreeIterator(const class ON_RTree& rtree)
00344 { 
00345   Initialize(rtree);
00346 }
00347 
00348 ON_RTreeIterator::~ON_RTreeIterator()
00349 { 
00350 }
00351 
00352 const ON_RTreeBranch* ON_RTreeIterator::Value() const
00353 {
00354   return ( 0 != m_sp )
00355          ? &m_sp->m_node->m_branch[m_sp->m_branchIndex]
00356          : 0;
00357 }
00358 
00359 bool ON_RTreeIterator::Initialize(const ON_RTree& a_rtree)
00360 {
00361   return Initialize(a_rtree.Root());
00362 }
00363 
00364 bool ON_RTreeIterator::Initialize(const ON_RTreeNode* a_node)
00365 { 
00366   m_sp = 0;
00367   m_root = ( 0 != a_node && a_node->m_count > 0 ) ? a_node : 0;
00368   return First();
00369 }
00370 
00371 bool ON_RTreeIterator::PushChildren(StackElement* sp, bool bFirstChild )
00372 { 
00373   StackElement* spmax = &m_stack[0] + MAX_STACK;
00374   const ON_RTreeNode* node = sp->m_node;
00375   m_sp = 0;
00376   // push first leaf coverted by this node onto the stack
00377   while( 0 != node && node->m_level >= 0 && node->m_count > 0 )
00378   {
00379     if ( 0 == node->m_level )
00380     {
00381       m_sp = sp;
00382       return true;
00383     }
00384     node = node->m_branch[sp->m_branchIndex].m_child;
00385     if( ++sp == spmax )
00386     {
00387       // Either this is a GIGANTIC R-tree, or, more likely, there is
00388       // a bug in the code that creates the R-tree and this R-tree
00389       // is horribly unbalanced. If the case is valid, then we must 
00390       // increase MAX_STACK and ship a service release.
00391       ON_ERROR("ON_RTreeIterator::PushFirstChild - stack overflow");
00392       return false;
00393     }
00394     sp->m_node = node;
00395     sp->m_branchIndex = bFirstChild ? 0 : node->m_count-1; // 0 for first child
00396   }
00397   return false;
00398 }
00399 
00400 bool ON_RTreeIterator::First()
00401 { 
00402   m_sp = 0;
00403   if ( 0 == m_root || m_root->m_level < 0 || m_root->m_count <= 0 )
00404     return false;
00405   m_stack[0].m_node = m_root;
00406   m_stack[0].m_branchIndex = 0;
00407   return PushChildren(&m_stack[0],true);
00408 }
00409 
00410 bool ON_RTreeIterator::Last()
00411 { 
00412   m_sp = 0;
00413   if ( 0 == m_root || m_root->m_level < 0 || m_root->m_count <= 0 )
00414     return false;
00415   m_stack[0].m_node = m_root;
00416   m_stack[0].m_branchIndex = m_root->m_count - 1;
00417   return PushChildren(&m_stack[0],false);
00418 }
00419 
00420 bool ON_RTreeIterator::Next()
00421 { 
00422   if ( 0 == m_sp )
00423     return false; // invalid iterator
00424 
00425   if ( ++(m_sp->m_branchIndex) < m_sp->m_node->m_count )
00426     return true; // m_sp->m_node is always at leaf level
00427 
00428   // pop the stack until we find an element with room to move over.
00429   StackElement* sp0 = &m_stack[0];
00430   StackElement* sp = m_sp;
00431   m_sp = 0;
00432   while ( sp > sp0 )
00433   {        
00434     sp--; // pop the stack
00435 
00436     // ++(sp->m_branchIndex) moves to the next element in sp->m_node.
00437     if ( ++(sp->m_branchIndex) >= sp->m_node->m_count )
00438       continue; // this element is used up
00439 
00440     // Since we've popped the stack, we cannot be at the leaf level.
00441     // PushFirst() pushes the first child onto the stack until 
00442     // it reaches the leaf level.
00443     return PushChildren(sp,true);
00444   }
00445   return false; // we were at the last element and now there are no more.
00446 }
00447 
00448 bool ON_RTreeIterator::Prev()
00449 { 
00450   if ( 0 == m_sp )
00451     return false; // invalid iterator
00452 
00453   if ( --(m_sp->m_branchIndex) >= 0 )
00454     return true; // m_sp->m_node is always at leaf level
00455 
00456   // pop the stack until we find an element with room to move over.
00457   StackElement* sp0 = &m_stack[0];
00458   StackElement* sp = m_sp;
00459   m_sp = 0;
00460   while ( sp > sp0 )
00461   {        
00462     sp--; // pop the stack
00463 
00464     // --(sp->m_branchIndex) moves to the previous element in sp->m_node.
00465     if ( --(sp->m_branchIndex) < 0 )
00466       continue; // this element is used up
00467 
00468     // Since we've popped the stack, we cannot be at the leaf level.
00469     // PushFirst() pushes the first child onto the stack until 
00470     // it reaches the leaf level.
00471     return PushChildren(sp,false);
00472   }
00473   return false; // we were at the last element and now there are no more.
00474 }
00475 
00476 
00478 //
00479 // ON_RTree
00480 //
00481 
00482 
00483 ON_RTree::ON_RTree( ON_MEMORY_POOL* heap, size_t leaf_count )
00484 : m_root(0)
00485 , m_reserved(0)
00486 , m_mem_pool(heap,leaf_count)
00487 {
00488 }
00489 
00490 
00491 
00492 ON_RTree::~ON_RTree()
00493 {
00494   RemoveAll();
00495 }
00496 
00497 bool ON_RTree::CreateMeshFaceTree( const ON_Mesh* mesh )
00498 {
00499   double fmin[3], fmax[3];
00500   ON_3dPoint V;
00501   unsigned int fi, fcount;
00502   const int* fvi;
00503   const ON_MeshFace* meshF;
00504   const ON_3fPoint* meshfV;
00505   const ON_3dPoint* meshdV;
00506 
00507   RemoveAll();
00508 
00509   if ( 0 == mesh )
00510     return false;
00511 
00512   fcount = mesh->m_F.UnsignedCount();
00513   if ( fcount <= 0 )
00514     return false;
00515 
00516   meshF = mesh->m_F.Array();
00517   if ( 0 == meshF )
00518     return false;
00519 
00520   meshfV = mesh->m_V.Array();
00521 
00522   meshdV = mesh->HasDoublePrecisionVertices() 
00523          ? mesh->DoublePrecisionVertices().Array() 
00524          : 0;
00525 
00526   if ( 0 != meshfV )
00527   {
00528     if ( 0 != meshdV )
00529     {
00530       for ( fi = 0; fi < fcount; fi++ )
00531       {
00532         fvi = meshF[fi].vi;
00533 
00534         V = meshfV[fvi[0]];
00535         fmin[0] = fmax[0] = V.x;
00536         fmin[1] = fmax[1] = V.y;
00537         fmin[2] = fmax[2] = V.z;
00538         V = meshdV[fvi[0]];
00539         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00540         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00541         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00542 
00543         V = meshfV[fvi[1]];
00544         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00545         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00546         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00547         V = meshdV[fvi[1]];
00548         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00549         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00550         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00551 
00552         V = meshfV[fvi[2]];
00553         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00554         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00555         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00556         V = meshdV[fvi[2]];
00557         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00558         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00559         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00560 
00561         if ( fvi[2] != fvi[3] )
00562         {
00563           V = meshfV[fvi[3]];
00564           if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00565           if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00566           if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;      
00567           V = meshdV[fvi[3]];
00568           if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00569           if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00570           if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00571         }
00572 
00573         if ( !Insert(fmin,fmax,fi) )
00574         {
00575           RemoveAll();
00576           return false;
00577         }
00578       }
00579     }
00580     else
00581     {
00582       for ( fi = 0; fi < fcount; fi++ )
00583       {
00584         fvi = meshF[fi].vi;
00585 
00586         V = meshfV[fvi[0]];
00587         fmin[0] = fmax[0] = V.x;
00588         fmin[1] = fmax[1] = V.y;
00589         fmin[2] = fmax[2] = V.z;
00590 
00591         V = meshfV[fvi[1]];
00592         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00593         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00594         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00595 
00596         V = meshfV[fvi[2]];
00597         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00598         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00599         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00600 
00601         if ( fvi[2] != fvi[3] )
00602         {
00603           V = meshfV[fvi[3]];
00604           if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00605           if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00606           if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;      
00607         }
00608 
00609         if ( !Insert(fmin,fmax,fi) )
00610         {
00611           RemoveAll();
00612           return false;
00613         }
00614       }
00615     }
00616   }
00617   else if ( 0 != meshdV )
00618   {
00619     for ( fi = 0; fi < fcount; fi++ )
00620     {
00621       fvi = meshF[fi].vi;
00622 
00623       V = meshdV[fvi[0]];
00624       fmin[0] = fmax[0] = V.x;
00625       fmin[1] = fmax[1] = V.y;
00626       fmin[2] = fmax[2] = V.z;
00627 
00628       V = meshdV[fvi[1]];
00629       if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00630       if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00631       if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00632 
00633       V = meshdV[fvi[2]];
00634       if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00635       if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00636       if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;  
00637 
00638       if ( fvi[2] != fvi[3] )
00639       {
00640         V = meshdV[fvi[3]];
00641         if ( V.x < fmin[0] ) fmin[0] = V.x; else if ( V.x > fmax[0] ) fmax[0] = V.x;
00642         if ( V.y < fmin[1] ) fmin[1] = V.y; else if ( V.y > fmax[1] ) fmax[1] = V.y;
00643         if ( V.z < fmin[2] ) fmin[2] = V.z; else if ( V.z > fmax[2] ) fmax[2] = V.z;      
00644       }
00645 
00646       if ( !Insert(fmin,fmax,fi) )
00647       {
00648         RemoveAll();
00649         return false;
00650       }
00651     }
00652   }
00653   else
00654   {
00655     // no vertices
00656     return false;
00657   }
00658 
00659   return (0 != m_root);
00660 }
00661 
00662 bool ON_RTree::Insert2d(const double a_min[2], const double a_max[2], int a_element_id)
00663 {
00664   const double min3d[3] = {a_min[0],a_min[1],0.0};
00665   const double max3d[3] = {a_max[0],a_max[1],0.0};
00666   return Insert(min3d,max3d,a_element_id);
00667 }
00668 
00669 bool ON_RTree::Insert(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM], int a_element_id)
00670 {
00671   bool rc;
00672   ON_RTreeBBox rect;
00673   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00674   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00675   if ( rect.m_min[0] <= rect.m_max[0] && rect.m_min[1] <= rect.m_max[1] && rect.m_min[2] <= rect.m_max[2] )
00676   {  
00677     if ( 0 == m_root )
00678     {
00679       m_root = m_mem_pool.AllocNode();
00680       m_root->m_level = 0;
00681     }
00682     InsertRect(&rect, a_element_id, &m_root, 0);
00683     rc = true;
00684   }
00685   else
00686   {
00687     // invalid bounding box - don't let this corrupt the tree
00688     rc = false;
00689     ON_ERROR("ON_RTree::Insert - invalid a_min[] or a_max[] input.");
00690   }
00691   return rc;
00692 }
00693 
00694 bool ON_RTree::Insert2d(const double a_min[2], const double a_max[2], void* a_element_id)
00695 {
00696   const double min3d[3] = {a_min[0],a_min[1],0.0};
00697   const double max3d[3] = {a_max[0],a_max[1],0.0};
00698   return Insert(min3d,max3d,a_element_id);
00699 }
00700 
00701 bool ON_RTree::Insert(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM], void* a_element_id)
00702 {
00703   bool rc;
00704   ON_RTreeBBox rect;
00705   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00706   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00707   if ( rect.m_min[0] <= rect.m_max[0] && rect.m_min[1] <= rect.m_max[1] && rect.m_min[2] <= rect.m_max[2] )
00708   {  
00709     if ( 0 == m_root )
00710     {
00711       m_root = m_mem_pool.AllocNode();
00712       m_root->m_level = 0;
00713     }
00714 
00715     // The ON__INT_PTR cast is safe because ON__INT_PTR == sizeof(void*)
00716 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
00717 #pragma warning( push )
00718 // Disable warning C4311: 'type cast' : pointer truncation from 'void *' to 'ON__INT_PTR'
00719 #pragma warning( disable : 4311 )
00720 #endif
00721     InsertRect(&rect, (ON__INT_PTR)a_element_id, &m_root, 0);
00722 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
00723 #pragma warning( pop )
00724 #endif
00725 
00726     rc = true;
00727   }
00728   else
00729   {
00730     // invalid bounding box - don't let this corrupt the tree
00731     rc = false;
00732     ON_ERROR("ON_RTree::Insert - invalid a_min[] or a_max[] input.");
00733   }
00734   return rc;
00735 }
00736 
00737 bool ON_RTree::Remove2d(const double a_min[2], const double a_max[2], int a_dataId)
00738 {
00739   const double min3d[3] = {a_min[0],a_min[1],0.0};
00740   const double max3d[3] = {a_max[0],a_max[1],0.0};
00741   return Remove(min3d,max3d,a_dataId);
00742 }
00743 
00744 bool ON_RTree::Remove(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM], int a_dataId)
00745 {
00746   bool rc = false;
00747   if ( 0 != m_root )
00748   {
00749     ON_RTreeBBox rect;
00750     memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00751     memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00752     if ( rect.m_min[0] <= rect.m_max[0] && rect.m_min[1] <= rect.m_max[1] && rect.m_min[2] <= rect.m_max[2] )
00753     {  
00754       // RemoveRect() returns 0 on success
00755       rc = (0 ==  RemoveRect(&rect, a_dataId, &m_root));
00756     }
00757     else
00758     {
00759       // invalid bounding box - don't let this corrupt the tree
00760       ON_ERROR("ON_RTree::Remove - invalid a_min[] or a_max[] input.");
00761     }
00762   }
00763   return rc;
00764 }
00765 
00766 bool ON_RTree::Remove2d(const double a_min[2], const double a_max[2],  void* a_dataId)
00767 {
00768   const double min3d[3] = {a_min[0],a_min[1],0.0};
00769   const double max3d[3] = {a_max[0],a_max[1],0.0};
00770   return Remove(min3d,max3d,a_dataId);
00771 }
00772 
00773 bool ON_RTree::Remove(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM], void* a_dataId)
00774 {
00775   bool rc = false;
00776   if ( 0 != m_root )
00777   {
00778     ON_RTreeBBox rect;
00779     memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00780     memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00781     if ( rect.m_min[0] <= rect.m_max[0] && rect.m_min[1] <= rect.m_max[1] && rect.m_min[2] <= rect.m_max[2] )
00782     {  
00783       // RemoveRect() returns 0 on success
00784       // The ON__INT_PTR cast is save because ON__INT_PTR == sizeof(void*)
00785 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
00786 #pragma warning( push )
00787 // Disable warning C4311: 'type cast' : pointer truncation from 'void *' to 'ON__INT_PTR'
00788 #pragma warning( disable : 4311 )
00789 #endif
00790       rc = (0 ==  RemoveRect(&rect, (ON__INT_PTR)a_dataId, &m_root));
00791 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
00792 #pragma warning( pop )
00793 #endif
00794 
00795     }
00796     else
00797     {
00798       // invalid bounding box - don't let this corrupt the tree
00799       ON_ERROR("ON_RTree::Remove - invalid a_min[] or a_max[] input.");
00800     }
00801   }
00802   return rc;
00803 }
00804 
00805 
00806 bool ON_RTree::Search2d(const double a_min[2], const double a_max[2], 
00807                           bool ON_MSC_CDECL a_resultCallback(void* a_context, ON__INT_PTR a_data), 
00808                           void* a_context
00809                           ) const
00810 {
00811   if ( 0 == m_root )
00812     return false;
00813 
00814   ON_RTreeBBox rect;
00815   memcpy(rect.m_min,a_min,2*sizeof(a_min[0]));
00816   rect.m_min[2] = 0.0;
00817   memcpy(rect.m_max,a_max,2*sizeof(a_max[0]));
00818   rect.m_max[2] = 0.0;
00819 
00820   ON_RTreeSearchResultCallback result;
00821   result.m_context = a_context;
00822   result.m_resultCallback = a_resultCallback;
00823   return SearchHelper(m_root, &rect, result);
00824 }
00825 
00826 bool ON_RTree::Search(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM], 
00827                           bool ON_MSC_CDECL a_resultCallback(void* a_context, ON__INT_PTR a_data), 
00828                           void* a_context
00829                           ) const
00830 {
00831   if ( 0 == m_root )
00832     return false;
00833 
00834   ON_RTreeBBox rect;
00835   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00836   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00837   return Search( &rect, a_resultCallback, a_context );
00838 }
00839 
00840 bool ON_RTree::Search( ON_RTreeBBox* a_rect,
00841                        bool ON_MSC_CDECL a_resultCallback(void* a_context, ON__INT_PTR a_data), 
00842                        void* a_context
00843                       ) const
00844 {
00845   if ( 0 == m_root || 0 == a_rect )
00846     return false;
00847 
00848   ON_RTreeSearchResultCallback result;
00849   result.m_context = a_context;
00850   result.m_resultCallback = a_resultCallback;
00851   return SearchHelper(m_root, a_rect, result);
00852 }
00853 
00854 bool ON_RTree::Search( 
00855     struct ON_RTreeSphere* a_sphere,
00856     bool ON_MSC_CDECL a_resultCallback(void* a_context, ON__INT_PTR a_id), 
00857     void* a_context
00858     ) const
00859 {
00860   if ( 0 == m_root || 0 == a_sphere )
00861     return false;
00862 
00863   ON_RTreeSearchResultCallback result;
00864   result.m_context = a_context;
00865   result.m_resultCallback = a_resultCallback;
00866 
00867   return SearchHelper(m_root, a_sphere, result);
00868 }
00869 
00870 bool ON_RTree::Search( 
00871   struct ON_RTreeCapsule* a_capsule,
00872   bool ON_MSC_CDECL a_resultCallback(void* a_context, ON__INT_PTR a_id), 
00873   void* a_context
00874   ) const
00875 {
00876   
00877   if ( 0 == m_root || 0 == a_capsule )
00878     return false;
00879 
00880   ON_RTreeSearchResultCallback result;
00881   result.m_context = a_context;
00882   result.m_resultCallback = a_resultCallback;
00883 
00884   return SearchHelper(m_root, a_capsule, result);
00885 }
00886 
00887 bool ON_RTree::Search2d(const double a_min[2], const double a_max[2],
00888                           ON_SimpleArray<ON_RTreeLeaf>& a_result 
00889                           ) const
00890 {
00891   if ( 0 == m_root )
00892     return false;
00893 
00894   ON_RTreeBBox rect;
00895   memcpy(rect.m_min,a_min,2*sizeof(a_min[0]));
00896   rect.m_min[2] = 0.0;
00897   memcpy(rect.m_max,a_max,2*sizeof(a_max[0]));
00898   rect.m_max[2] = 0.0;
00899 
00900   return SearchHelper(m_root, &rect, a_result);
00901 }
00902 
00903 
00904 bool ON_RTree::Search(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM],
00905                           ON_SimpleArray<ON_RTreeLeaf>& a_result 
00906                           ) const
00907 {
00908   if ( 0 == m_root )
00909     return false;
00910 
00911   ON_RTreeBBox rect;
00912   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00913   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00914 
00915   return SearchHelper(m_root, &rect, a_result);
00916 }
00917 
00918 
00919 bool ON_RTree::Search2d(const double a_min[2], const double a_max[2],
00920                           ON_SimpleArray<void*>& a_result 
00921                           ) const
00922 {
00923   if ( 0 == m_root )
00924     return false;
00925 
00926   ON_RTreeBBox rect;
00927   memcpy(rect.m_min,a_min,2*sizeof(a_min[0]));
00928   rect.m_min[2] = 0.0;
00929   memcpy(rect.m_max,a_max,2*sizeof(a_max[0]));
00930   rect.m_max[2] = 0.0;
00931 
00932   return SearchHelper(m_root, &rect,  a_result);
00933 }
00934 
00935 
00936 bool ON_RTree::Search(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM],
00937                           ON_SimpleArray<void*>& a_result 
00938                           ) const
00939 {
00940   if ( 0 == m_root )
00941     return false;
00942 
00943   ON_RTreeBBox rect;
00944   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00945   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00946 
00947   return SearchHelper(m_root, &rect,  a_result);
00948 }
00949 
00950 
00951 bool ON_RTree::Search2d(const double a_min[2], const double a_max[2],
00952                           ON_SimpleArray<int>& a_result 
00953                           ) const
00954 {
00955   if ( 0 == m_root )
00956     return false;
00957 
00958   ON_RTreeBBox rect;
00959   memcpy(rect.m_min,a_min,2*sizeof(a_min[0]));
00960   rect.m_min[2] = 0.0;
00961   memcpy(rect.m_max,a_max,2*sizeof(a_max[0]));
00962   rect.m_max[2] = 0.0;
00963 
00964   return SearchHelper(m_root, &rect,  a_result);
00965 }
00966 
00967 bool ON_RTree::Search(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM],
00968                           ON_SimpleArray<int>& a_result 
00969                           ) const
00970 {
00971   if ( 0 == m_root )
00972     return false;
00973 
00974   ON_RTreeBBox rect;
00975   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
00976   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
00977 
00978   return SearchHelper(m_root, &rect,  a_result);
00979 }
00980 
00981 bool ON_RTree::Search2d(const double a_min[2], const double a_max[2],
00982                           ON_RTreeSearchResult& a_result ) const
00983 {
00984   if ( 0 == m_root )
00985     return false;
00986 
00987   ON_RTreeBBox rect;
00988   memcpy(rect.m_min,a_min,2*sizeof(a_min[0]));
00989   rect.m_min[2] = 0.0;
00990   memcpy(rect.m_max,a_max,2*sizeof(a_max[0]));
00991   rect.m_max[2] = 0.0;
00992 
00993   return SearchHelper(m_root, &rect, a_result);
00994 }
00995 
00996 bool ON_RTree::Search(const double a_min[ON_RTree_NODE_DIM], const double a_max[ON_RTree_NODE_DIM],
00997                           ON_RTreeSearchResult& a_result ) const
00998 {
00999   if ( 0 == m_root )
01000     return false;
01001 
01002   ON_RTreeBBox rect;
01003   memcpy(rect.m_min,a_min,sizeof(rect.m_min));
01004   memcpy(rect.m_max,a_max,sizeof(rect.m_max));
01005 
01006   return SearchHelper(m_root, &rect, a_result);
01007 }
01008 
01009 struct ON_RTreePairSearchResult
01010 {
01011   double m_tolerance;
01012   ON_SimpleArray<ON_2dex>* m_result;
01013 };
01014 
01015 
01016 static bool PairSearchOverlapHelper( const ON_RTreeBBox* a_rectA, const ON_RTreeBBox* a_rectB, double tolerance )
01017 {
01018   double dx,dy,dz,d;
01019   const double* mn;
01020   const double* mx;
01021 
01022 
01023   mx = a_rectA->m_max;
01024   mn = a_rectB->m_min;
01025   dx = *mn++ - *mx++;
01026   if ( dx > tolerance ) return false;
01027   dy = *mn++ - *mx++;
01028   if ( dy > tolerance ) return false;
01029   dz = *mn - *mx;
01030   if ( dz > tolerance ) return false;
01031 
01032   mx = a_rectB->m_max;
01033   mn = a_rectA->m_min;
01034   d = *mn++ - *mx++;
01035   if ( d > tolerance ) return false;
01036   if ( d > dx ) dx = d;
01037   d = *mn++ - *mx++;
01038   if ( d > tolerance ) return false;
01039   if ( d > dy ) dy = d;
01040   d = *mn++ - *mx++;
01041   if ( d > tolerance ) return false;
01042   if ( d > dz ) dz = d;
01043 
01044   d  = (dx > 0.0) ? dx*dx : 0.0;
01045   d += (dy > 0.0) ? dy*dy : 0.0;
01046   d += (dz > 0.0) ? dz*dz : 0.0;
01047 
01048   return (d <= tolerance*tolerance);
01049 }
01050 
01051 
01052 static void PairSearchHelper( const ON_RTreeBranch* a_branchA, const ON_RTreeNode* a_nodeB, ON_RTreePairSearchResult* a_result )
01053 {
01054   // DO NOT ADD ANYTHING TO THIS FUNCTION
01055   const ON_RTreeBranch *branchB, *branchBmax;
01056 
01057   branchB = a_nodeB->m_branch;
01058   branchBmax = branchB + a_nodeB->m_count;
01059   while(branchB < branchBmax)
01060   {
01061     if ( PairSearchOverlapHelper( &a_branchA->m_rect, &branchB->m_rect, a_result->m_tolerance ) )
01062     {
01063       if ( a_nodeB->m_level > 0 )
01064       {
01065         PairSearchHelper(a_branchA,branchB->m_child,a_result);
01066       }
01067       else
01068       {
01069         ON_2dex& r = a_result->m_result->AppendNew();
01070         r.i = (int)a_branchA->m_id;
01071         r.j = (int)branchB->m_id;
01072       }
01073     }
01074     branchB++;
01075   }
01076 }
01077 
01078 static void PairSearchHelper( const ON_RTreeNode* a_nodeA, const ON_RTreeBranch* a_branchB, ON_RTreePairSearchResult* a_result )
01079 {
01080   // DO NOT ADD ANYTHING TO THIS FUNCTION
01081   const ON_RTreeBranch *branchA, *branchAmax;
01082 
01083   branchA = a_nodeA->m_branch;
01084   branchAmax = branchA + a_nodeA->m_count;
01085   while(branchA < branchAmax)
01086   {
01087     if ( PairSearchOverlapHelper( &branchA->m_rect, &a_branchB->m_rect, a_result->m_tolerance ) )
01088     {
01089       if ( a_nodeA->m_level > 0 )
01090       {
01091         PairSearchHelper(branchA->m_child,a_branchB,a_result);
01092       }
01093       else
01094       {
01095         ON_2dex& r = a_result->m_result->AppendNew();
01096         r.i = (int)branchA->m_id;
01097         r.j = (int)a_branchB->m_id;
01098       }
01099     }
01100     branchA++;
01101   }
01102 }
01103 
01104 
01105 static void PairSearchHelper( const ON_RTreeNode* a_nodeA, const ON_RTreeNode* a_nodeB, ON_RTreePairSearchResult* a_result )
01106 {
01107   // DO NOT ADD ANYTHING TO THIS FUNCTION
01108   const ON_RTreeBranch *branchA, *branchAmax, *branchB, *branchBmax;
01109 
01110   branchA = a_nodeA->m_branch;
01111   branchAmax = branchA + a_nodeA->m_count;
01112   branchBmax = a_nodeB->m_branch + a_nodeB->m_count;
01113   while(branchA < branchAmax)
01114   {
01115     for ( branchB = a_nodeB->m_branch; branchB < branchBmax; branchB++ )
01116     {
01117       if ( PairSearchOverlapHelper( &branchA->m_rect, &branchB->m_rect, a_result->m_tolerance ) )
01118       {
01119         if ( a_nodeA->m_level > 0 )
01120         {
01121           if ( a_nodeB->m_level > 0 )
01122             PairSearchHelper(branchA->m_child,branchB->m_child,a_result);
01123           else
01124             PairSearchHelper(branchA->m_child,branchB,a_result);
01125         }
01126         else if ( a_nodeB->m_level > 0 )
01127         {
01128           PairSearchHelper(branchA,branchB->m_child,a_result);
01129         }
01130         else
01131         {
01132           ON_2dex& r = a_result->m_result->AppendNew();
01133           r.i = (int)branchA->m_id;
01134           r.j = (int)branchB->m_id;
01135         }
01136       }
01137     }
01138     branchA++;
01139   }
01140 }
01141 
01142 
01143 bool ON_RTree::Search( 
01144           const ON_RTree& a_rtreeA,
01145           const ON_RTree& a_rtreeB, 
01146           double tolerance,
01147           ON_SimpleArray<ON_2dex>& a_result
01148           )
01149 {
01150   if ( 0 == a_rtreeA.m_root )
01151     return false;
01152   if ( 0 == a_rtreeB.m_root )
01153     return false;
01154   ON_RTreePairSearchResult r;
01155   r.m_tolerance = (ON_IsValid(tolerance) && tolerance > 0.0) ? tolerance : 0.0;
01156   r.m_result = &a_result;
01157   PairSearchHelper(a_rtreeA.m_root,a_rtreeB.m_root,&r);
01158   return true;
01159 }
01160 
01161 typedef void (*ON_RTreePairSearchCallback)(void*, ON__INT_PTR, ON__INT_PTR);
01162 
01163 struct ON_RTreePairSearchCallbackResult
01164 {
01165   double m_tolerance;
01166   void* m_context;
01167   ON_RTreePairSearchCallback m_resultCallback;
01168 };
01169 
01170 typedef bool (*ON_RTreePairSearchCallbackBool)(void*, ON__INT_PTR, ON__INT_PTR);
01171 
01172 struct ON_RTreePairSearchCallbackResultBool
01173 {
01174   double m_tolerance;
01175   void* m_context;
01176   ON_RTreePairSearchCallbackBool m_resultCallbackBool;
01177 };
01178 
01179 static void PairSearchHelper( const ON_RTreeBranch* a_branchA, const ON_RTreeNode* a_nodeB, ON_RTreePairSearchCallbackResult* a_result )
01180 {
01181   // DO NOT ADD ANYTHING TO THIS FUNCTION
01182   const ON_RTreeBranch *branchB, *branchBmax;
01183 
01184   branchB = a_nodeB->m_branch;
01185   branchBmax = branchB + a_nodeB->m_count;
01186   while(branchB < branchBmax)
01187   {
01188     if ( PairSearchOverlapHelper( &a_branchA->m_rect, &branchB->m_rect, a_result->m_tolerance ) )
01189     {
01190       if ( a_nodeB->m_level > 0 )
01191       {
01192         PairSearchHelper(a_branchA,branchB->m_child,a_result);
01193       }
01194       else
01195       {
01196         a_result->m_resultCallback(a_result->m_context,a_branchA->m_id,branchB->m_id);
01197       }
01198     }
01199     branchB++;
01200   }
01201 }
01202 
01203 static bool PairSearchHelperBool( const ON_RTreeBranch* a_branchA, const ON_RTreeNode* a_nodeB, ON_RTreePairSearchCallbackResultBool* a_result )
01204 {
01205   // DO NOT ADD ANYTHING TO THIS FUNCTION
01206   const ON_RTreeBranch *branchB, *branchBmax;
01207 
01208   branchB = a_nodeB->m_branch;
01209   branchBmax = branchB + a_nodeB->m_count;
01210   while(branchB < branchBmax)
01211   {
01212     if ( PairSearchOverlapHelper( &a_branchA->m_rect, &branchB->m_rect, a_result->m_tolerance ) )
01213     {
01214       if ( a_nodeB->m_level > 0 )
01215       {
01216         if ( !PairSearchHelperBool(a_branchA,branchB->m_child,a_result) )
01217           return false;
01218       }
01219       else
01220       {
01221         if ( !a_result->m_resultCallbackBool(a_result->m_context,a_branchA->m_id,branchB->m_id) )
01222           return false;
01223       }
01224     }
01225     branchB++;
01226   }
01227   return true;
01228 }
01229 
01230 static void PairSearchHelper( const ON_RTreeNode* a_nodeA, const ON_RTreeBranch* a_branchB, ON_RTreePairSearchCallbackResult* a_result )
01231 {
01232   // DO NOT ADD ANYTHING TO THIS FUNCTION
01233   const ON_RTreeBranch *branchA, *branchAmax;
01234 
01235   branchA = a_nodeA->m_branch;
01236   branchAmax = branchA + a_nodeA->m_count;
01237   while(branchA < branchAmax)
01238   {
01239     if ( PairSearchOverlapHelper( &branchA->m_rect, &a_branchB->m_rect, a_result->m_tolerance ) )
01240     {
01241       if ( a_nodeA->m_level > 0 )
01242       {
01243         PairSearchHelper(branchA->m_child,a_branchB,a_result);
01244       }
01245       else
01246       {
01247         a_result->m_resultCallback(a_result->m_context,branchA->m_id,a_branchB->m_id);
01248       }
01249     }
01250     branchA++;
01251   }
01252 }
01253 
01254 static bool PairSearchHelperBool( const ON_RTreeNode* a_nodeA, const ON_RTreeBranch* a_branchB, ON_RTreePairSearchCallbackResultBool* a_result )
01255 {
01256   // DO NOT ADD ANYTHING TO THIS FUNCTION
01257   const ON_RTreeBranch *branchA, *branchAmax;
01258 
01259   branchA = a_nodeA->m_branch;
01260   branchAmax = branchA + a_nodeA->m_count;
01261   while(branchA < branchAmax)
01262   {
01263     if ( PairSearchOverlapHelper( &branchA->m_rect, &a_branchB->m_rect, a_result->m_tolerance ) )
01264     {
01265       if ( a_nodeA->m_level > 0 )
01266       {
01267         if ( !PairSearchHelperBool(branchA->m_child,a_branchB,a_result) )
01268           return false;
01269       }
01270       else
01271       {
01272         if ( !a_result->m_resultCallbackBool(a_result->m_context,branchA->m_id,a_branchB->m_id) )
01273           return false;
01274       }
01275     }
01276     branchA++;
01277   }
01278   return true;
01279 }
01280 
01281 
01282 static void PairSearchHelper( const ON_RTreeNode* a_nodeA, const ON_RTreeNode* a_nodeB, ON_RTreePairSearchCallbackResult* a_result )
01283 {
01284   // DO NOT ADD ANYTHING TO THIS FUNCTION
01285   const ON_RTreeBranch *branchA, *branchAmax, *branchB, *branchBmax;
01286 
01287   branchA = a_nodeA->m_branch;
01288   branchAmax = branchA + a_nodeA->m_count;
01289   branchBmax = a_nodeB->m_branch + a_nodeB->m_count;
01290   while(branchA < branchAmax)
01291   {
01292     for ( branchB = a_nodeB->m_branch; branchB < branchBmax; branchB++ )
01293     {
01294       if ( PairSearchOverlapHelper( &branchA->m_rect, &branchB->m_rect, a_result->m_tolerance ) )
01295       {
01296         if ( a_nodeA->m_level > 0 )
01297         {
01298           if ( a_nodeB->m_level > 0 )
01299             PairSearchHelper(branchA->m_child,branchB->m_child,a_result);
01300           else
01301             PairSearchHelper(branchA->m_child,branchB,a_result);
01302         }
01303         else if ( a_nodeB->m_level > 0 )
01304         {
01305           PairSearchHelper(branchA,branchB->m_child,a_result);
01306         }
01307         else
01308         {
01309           a_result->m_resultCallback(a_result->m_context,branchA->m_id,branchB->m_id);
01310         }
01311       }
01312     }
01313     branchA++;
01314   }
01315 }
01316 
01317 
01318 static bool PairSearchHelperBool( const ON_RTreeNode* a_nodeA, const ON_RTreeNode* a_nodeB, ON_RTreePairSearchCallbackResultBool* a_result )
01319 {
01320   // DO NOT ADD ANYTHING TO THIS FUNCTION
01321   const ON_RTreeBranch *branchA, *branchAmax, *branchB, *branchBmax;
01322 
01323   branchA = a_nodeA->m_branch;
01324   branchAmax = branchA + a_nodeA->m_count;
01325   branchBmax = a_nodeB->m_branch + a_nodeB->m_count;
01326   while(branchA < branchAmax)
01327   {
01328     for ( branchB = a_nodeB->m_branch; branchB < branchBmax; branchB++ )
01329     {
01330       if ( PairSearchOverlapHelper( &branchA->m_rect, &branchB->m_rect, a_result->m_tolerance ) )
01331       {
01332         if ( a_nodeA->m_level > 0 )
01333         {
01334           if ( a_nodeB->m_level > 0 )
01335           {
01336             if ( !PairSearchHelperBool(branchA->m_child,branchB->m_child,a_result) )
01337               return false;
01338           }
01339           else
01340           {
01341             if ( !PairSearchHelperBool(branchA->m_child,branchB,a_result) )
01342               return false;
01343           }
01344         }
01345         else if ( a_nodeB->m_level > 0 )
01346         {
01347           if ( !PairSearchHelperBool(branchA,branchB->m_child,a_result) )
01348             return false;
01349         }
01350         else
01351         {
01352           if ( !a_result->m_resultCallbackBool(a_result->m_context,branchA->m_id,branchB->m_id) )
01353             return false;
01354         }
01355       }
01356     }
01357     branchA++;
01358   }
01359   return true;
01360 }
01361 
01362 bool ON_RTree::Search( 
01363           const ON_RTree& a_rtreeA,
01364           const ON_RTree& a_rtreeB, 
01365           double tolerance,
01366           void ON_MSC_CDECL resultCallback(void* a_context,ON__INT_PTR a_idA, ON__INT_PTR a_idB),
01367           void* a_context
01368           )
01369 {
01370   if ( 0 == a_rtreeA.m_root )
01371     return false;
01372   if ( 0 == a_rtreeB.m_root )
01373     return false;
01374   ON_RTreePairSearchCallbackResult r;
01375   r.m_tolerance = (ON_IsValid(tolerance) && tolerance > 0.0) ? tolerance : 0.0;
01376   r.m_context = a_context;
01377   r.m_resultCallback = resultCallback;
01378   PairSearchHelper(a_rtreeA.m_root,a_rtreeB.m_root,&r);
01379   return true;
01380 }
01381 
01382 bool ON_RTree::Search( 
01383           const ON_RTree& a_rtreeA,
01384           const ON_RTree& a_rtreeB, 
01385           double tolerance,
01386           bool ON_MSC_CDECL resultCallback(void* a_context,ON__INT_PTR a_idA, ON__INT_PTR a_idB),
01387           void* a_context
01388           )
01389 {
01390   if ( 0 == a_rtreeA.m_root )
01391     return false;
01392   if ( 0 == a_rtreeB.m_root )
01393     return false;
01394   ON_RTreePairSearchCallbackResultBool r;
01395   r.m_tolerance = (ON_IsValid(tolerance) && tolerance > 0.0) ? tolerance : 0.0;
01396   r.m_context = a_context;
01397   r.m_resultCallbackBool = resultCallback;
01398 
01399   // Do not return false if PairSearchHelperBool() returns false.  The only reason
01400   // PairSearchHelperBool() returns false is that the user specified resultCallback()
01401   // terminated the search. This way a programmer with the ability to reason can
01402   // distinguish between a terminiation and a failure to start because input is
01403   // missing.
01404   PairSearchHelperBool(a_rtreeA.m_root,a_rtreeB.m_root,&r);
01405 
01406   return true;
01407 }
01408 
01409 
01410 
01411 
01412 int ON_RTree::ElementCount()
01413 {
01414   int count = 0;
01415 
01416   if ( 0 != m_root )
01417     CountRec(m_root, count);
01418   
01419   return count;
01420 }
01421 
01422 const ON_RTreeNode* ON_RTree::Root() const
01423 {
01424   return m_root;
01425 }
01426 
01427 ON_BoundingBox ON_RTree::BoundingBox() const
01428 {
01429   ON_BoundingBox bbox;
01430   if ( 0 != m_root && m_root->m_count > 0 )
01431   {
01432     bbox.m_min = m_root->m_branch[0].m_rect.m_min;
01433     bbox.m_max = m_root->m_branch[0].m_rect.m_max;
01434     for ( int i = 1; i < m_root->m_count; i++ )
01435     {
01436       if ( m_root->m_branch[i].m_rect.m_min[0] < bbox.m_min.x )
01437         bbox.m_min.x = m_root->m_branch[i].m_rect.m_min[0];
01438       if ( m_root->m_branch[i].m_rect.m_min[1] < bbox.m_min.y )
01439         bbox.m_min.y = m_root->m_branch[i].m_rect.m_min[1];
01440       if ( m_root->m_branch[i].m_rect.m_min[2] < bbox.m_min.z )
01441         bbox.m_min.z = m_root->m_branch[i].m_rect.m_min[2];
01442 
01443       if ( m_root->m_branch[i].m_rect.m_max[0] > bbox.m_max.x )
01444         bbox.m_max.x = m_root->m_branch[i].m_rect.m_max[0];
01445       if ( m_root->m_branch[i].m_rect.m_max[1] > bbox.m_max.y )
01446         bbox.m_max.y = m_root->m_branch[i].m_rect.m_max[1];
01447       if ( m_root->m_branch[i].m_rect.m_max[2] > bbox.m_max.z )
01448         bbox.m_max.z = m_root->m_branch[i].m_rect.m_max[2];
01449     }
01450   }
01451   return bbox;
01452 }
01453 
01454 static void CountRec(ON_RTreeNode* a_node, int& a_count)
01455 {
01456   if(a_node->IsInternalNode())  // not a leaf node
01457   {
01458     for(int index = 0; index < a_node->m_count; ++index)
01459     {
01460       CountRec(a_node->m_branch[index].m_child, a_count);
01461     }
01462   }
01463   else // A leaf node
01464   {
01465     a_count += a_node->m_count;
01466   }
01467 }
01468 
01469 size_t ON_RTree::SizeOf() const
01470 {
01471   return m_mem_pool.SizeOf();
01472 }
01473 
01474 
01475 static void NodeCountHelper( const ON_RTreeNode* node, size_t& node_count, size_t& wasted_branch_count, size_t& leaf_count )
01476 {
01477   if ( 0 == node )
01478     return;
01479   node_count++;
01480   wasted_branch_count += (ON_RTree_MAX_NODE_COUNT - node->m_count);
01481   if ( node->m_level > 0 )
01482   {
01483     for ( int i = 0; i < node->m_count; i++ )
01484     {
01485       NodeCountHelper(node->m_branch[i].m_child,node_count,wasted_branch_count,leaf_count);
01486     }
01487   }
01488   else
01489     leaf_count += node->m_count;
01490 }
01491 
01492 void ON_RTree::RemoveAll()
01493 {
01494   m_root = 0;
01495   m_mem_pool.DeallocateAll();
01496 }
01497 
01498 void ON_RTree::RemoveAllRec(ON_RTreeNode* a_node)
01499 {
01500   if(a_node->IsInternalNode()) // This is an internal node in the tree
01501   {
01502     for(int index=0; index < a_node->m_count; ++index)
01503     {
01504       RemoveAllRec(a_node->m_branch[index].m_child);
01505     }
01506   }
01507   m_mem_pool.FreeNode(a_node); 
01508 }
01509 
01510 
01511 
01512 static void InitRect(ON_RTreeBBox* a_rect)
01513 {
01514   for(int index = 0; index < ON_RTree_NODE_DIM; ++index)
01515   {
01516     a_rect->m_min[index] = (double)0;
01517     a_rect->m_max[index] = (double)0;
01518   }
01519 }
01520 
01521 
01522 // Inserts a new data rectangle into the index structure.
01523 // Recursively descends tree, propagates splits back up.
01524 // Returns 0 if node was not split.  Old node updated.
01525 // If node was split, returns 1 and sets the pointer pointed to by
01526 // new_node to point to the new node.  Old node updated to become one of two.
01527 // The level argument specifies the number of steps up from the leaf
01528 // level to insert; e.g. a data rectangle goes in at level = 0.
01529 
01530 bool ON_RTree::InsertRectRec(ON_RTreeBBox* a_rect, ON__INT_PTR a_id, ON_RTreeNode* a_node, ON_RTreeNode** a_newNode, int a_level)
01531 {
01532   int index;
01533   ON_RTreeBranch branch;
01534   ON_RTreeNode* otherNode;
01535 
01536   // Still above level for insertion, go down tree recursively
01537   if(a_node->m_level > a_level)
01538   {
01539     index = PickBranch(a_rect, a_node);
01540     if ( index < 0 )
01541     {
01542       return false;
01543     }
01544     if (!InsertRectRec(a_rect, a_id, a_node->m_branch[index].m_child, &otherNode, a_level))
01545     {
01546       // Child was not split
01547       a_node->m_branch[index].m_rect = CombineRectHelper(a_rect, &(a_node->m_branch[index].m_rect));
01548       return false;
01549     }
01550     else // Child was split
01551     {
01552       a_node->m_branch[index].m_rect = NodeCover(a_node->m_branch[index].m_child);
01553       branch.m_child = otherNode;
01554       branch.m_rect = NodeCover(otherNode);
01555       return AddBranch(&branch, a_node, a_newNode);
01556     }
01557   }
01558   else if(a_node->m_level == a_level) // Have reached level for insertion. Add rect, split if necessary
01559   {
01560     branch.m_rect = *a_rect;
01561 
01562     // The (ON_RTreeNode*) cast is safe because ON__INT_PTR == sizeof(void*)
01563 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
01564 #pragma warning( push )
01565 // Disable warning C4312: 'type cast' : conversion from 'ON__INT_PTR' to 'ON_RTreeNode *' of greater size
01566 #pragma warning( disable : 4312 )
01567 #endif
01568     branch.m_child = (ON_RTreeNode*)a_id;
01569 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
01570 #pragma warning( pop )
01571 #endif
01572 
01573     // Child field of leaves contains id of data record
01574     return AddBranch(&branch, a_node, a_newNode);
01575   }
01576 
01577   // We should never get here
01578   ON_ERROR("ON_RTree::InsertRectRec - bug in algorithm");
01579   return false;
01580 }
01581 
01582 
01583 // Insert a data rectangle into an index structure.
01584 // InsertRect provides for splitting the root;
01585 // returns 1 if root was split, 0 if it was not.
01586 // The level argument specifies the number of steps up from the leaf
01587 // level to insert; e.g. a data rectangle goes in at level = 0.
01588 // InsertRect2 does the recursion.
01589 //
01590 
01591 bool ON_RTree::InsertRect(ON_RTreeBBox* a_rect, ON__INT_PTR a_id, ON_RTreeNode** a_root, int a_level)
01592 {
01593   ON_RTreeNode* newRoot;
01594   ON_RTreeNode* newNode;
01595   ON_RTreeBranch branch;
01596 
01597   if(InsertRectRec(a_rect, a_id, *a_root, &newNode, a_level))  // Root split
01598   {
01599     newRoot = m_mem_pool.AllocNode();  // Grow tree taller and new root
01600     newRoot->m_level = (*a_root)->m_level + 1;
01601     branch.m_rect = NodeCover(*a_root);
01602     branch.m_child = *a_root;
01603     AddBranch(&branch, newRoot, NULL);
01604     branch.m_rect = NodeCover(newNode);
01605     branch.m_child = newNode;
01606     AddBranch(&branch, newRoot, NULL);
01607     *a_root = newRoot;
01608     return true;
01609   }
01610 
01611   return false;
01612 }
01613 
01614 
01615 // Find the smallest rectangle that includes all rectangles in branches of a node.
01616 
01617 static ON_RTreeBBox NodeCover(ON_RTreeNode* a_node)
01618 {
01619   int i;
01620   const ON_RTreeBranch* branch;
01621   ON_RTreeBBox rect;
01622 
01623   if ( (i = a_node->m_count) > 0 )
01624   {
01625     rect = a_node->m_branch[--i].m_rect;
01626     for ( branch = a_node->m_branch; i; i--, branch++ )
01627     {
01628 #if (3 == ON_RTree_NODE_DIM)
01629       if ( rect.m_min[0] > branch->m_rect.m_min[0] )
01630         rect.m_min[0] = branch->m_rect.m_min[0];
01631       if ( rect.m_min[1] > branch->m_rect.m_min[1] )
01632         rect.m_min[1] = branch->m_rect.m_min[1];
01633       if ( rect.m_min[2] > branch->m_rect.m_min[2] )
01634         rect.m_min[2] = branch->m_rect.m_min[2];
01635       if ( rect.m_max[0] < branch->m_rect.m_max[0] )
01636         rect.m_max[0] = branch->m_rect.m_max[0];
01637       if ( rect.m_max[1] < branch->m_rect.m_max[1] )
01638         rect.m_max[1] = branch->m_rect.m_max[1];
01639       if ( rect.m_max[2] < branch->m_rect.m_max[2] )
01640         rect.m_max[2] = branch->m_rect.m_max[2];
01641 #else
01642       for( int j = 0; j < ON_RTree_NODE_DIM; j++ )
01643       {
01644         if ( rect.m_min[j] > branch->m_rect.m_min[j] )
01645           rect.m_min[j] = branch->m_rect.m_min[j];
01646         if ( rect.m_max[j] < branch->m_rect.m_max[j] )
01647           rect.m_max[j] = branch->m_rect.m_max[j];
01648       }
01649 #endif
01650     }
01651   }
01652   else
01653   {
01654     InitRect(&rect);
01655   }
01656 
01657   return rect;
01658 }
01659 
01660 
01661 // Add a branch to a node.  Split the node if necessary.
01662 // Returns 0 if node not split.  Old node updated.
01663 // Returns 1 if node split, sets *new_node to address of new node.
01664 // Old node updated, becomes one of two.
01665 
01666 bool ON_RTree::AddBranch(ON_RTreeBranch* a_branch, ON_RTreeNode* a_node, ON_RTreeNode** a_newNode)
01667 {
01668   if(a_node->m_count < ON_RTree_MAX_NODE_COUNT)  // Split won't be necessary
01669   {
01670     a_node->m_branch[a_node->m_count] = *a_branch;
01671     ++a_node->m_count;
01672 
01673     return false;
01674   }
01675   else
01676   {
01677     SplitNode(a_node, a_branch, a_newNode);
01678     return true;
01679   }
01680 }
01681 
01682 
01683 // Disconnect a dependent node.
01684 // Caller must return (or stop using iteration index) after this as count has changed
01685 
01686 static void DisconnectBranch(ON_RTreeNode* a_node, int a_index)
01687 {
01688   // Remove element by swapping with the last element to prevent gaps in array
01689   a_node->m_branch[a_index] = a_node->m_branch[a_node->m_count - 1];
01690   
01691   --a_node->m_count;
01692 }
01693 
01694 
01695 // Pick a branch.  Pick the one that will need the smallest increase
01696 // in area to accomodate the new rectangle.  This will result in the
01697 // least total area for the covering rectangles in the current node.
01698 // In case of a tie, pick the one which was smaller before, to get
01699 // the best resolution when searching.
01700 
01701 static int PickBranch(ON_RTreeBBox* a_rect, ON_RTreeNode* a_node)
01702 {
01703   bool firstTime = true;
01704   double increase;
01705   double bestIncr = -1.0;
01706   double area;
01707   double bestArea;
01708   int best;
01709   ON_RTreeBBox tempRect;
01710 
01711   best = -1;
01712   bestArea = -1.0;
01713 
01714   for(int index=0; index < a_node->m_count; ++index)
01715   {
01716     ON_RTreeBBox* curRect = &a_node->m_branch[index].m_rect;
01717     area = CalcRectVolumeHelper(curRect);
01718     tempRect = CombineRectHelper(a_rect, curRect);
01719     increase = CalcRectVolumeHelper(&tempRect) - area;
01720     if((increase < bestIncr) || firstTime)
01721     {
01722       best = index;
01723       bestArea = area;
01724       bestIncr = increase;
01725       firstTime = false;
01726     }
01727     else if((increase == bestIncr) && (area <=bestArea))
01728     {
01729       best = index;
01730       bestArea = area;
01731       bestIncr = increase;
01732     }
01733   }
01734   return best;
01735 }
01736 
01737 
01738 // Combine two rectangles into larger one containing both
01739 
01740 ON_RTreeBBox CombineRectHelper(const ON_RTreeBBox* a_rectA, const ON_RTreeBBox* a_rectB)
01741 {
01742   ON_RTreeBBox rect = *a_rectA;
01743 
01744 #if (3 == ON_RTree_NODE_DIM)
01745   if ( rect.m_min[0] > a_rectB->m_min[0] )
01746     rect.m_min[0] = a_rectB->m_min[0];
01747   if ( rect.m_min[1] > a_rectB->m_min[1] )
01748     rect.m_min[1] = a_rectB->m_min[1];
01749   if ( rect.m_min[2] > a_rectB->m_min[2] )
01750     rect.m_min[2] = a_rectB->m_min[2];
01751   if ( rect.m_max[0] < a_rectB->m_max[0] )
01752     rect.m_max[0] = a_rectB->m_max[0];
01753   if ( rect.m_max[1] < a_rectB->m_max[1] )
01754     rect.m_max[1] = a_rectB->m_max[1];
01755   if ( rect.m_max[2] < a_rectB->m_max[2] )
01756     rect.m_max[2] = a_rectB->m_max[2];
01757 #else
01758   for( int j = 0; j < ON_RTree_NODE_DIM; j++ )
01759   {
01760     if ( rect.m_min[j] > a_rectB->m_min[j] )
01761       rect.m_min[j] = a_rectB->m_min[j];
01762     if ( rect.m_max[j] < a_rectB->m_max[j] )
01763       rect.m_max[j] = a_rectB->m_max[j];
01764   }
01765 #endif
01766 
01767   return rect;
01768 }
01769 
01770 
01771 
01772 // Split a node.
01773 // Divides the nodes branches and the extra one between two nodes.
01774 // Old node is one of the new ones, and one really new one is created.
01775 // Tries more than one method for choosing a partition, uses best result.
01776 
01777 void ON_RTree::SplitNode(ON_RTreeNode* a_node, ON_RTreeBranch* a_branch, ON_RTreeNode** a_newNode)
01778 {
01779   ON_RTreePartitionVars localVars;
01780   int level;
01781 
01782   // Load all the branches into a buffer, initialize a_node to be empty
01783   level = a_node->m_level; // save m_level (The InitNode() call in GetBranches will set it to -1)
01784   GetBranches(a_node, a_branch, &localVars);
01785 
01786   // Find partition
01787   ChoosePartition(&localVars, ON_RTree_MIN_NODE_COUNT);
01788 
01789   // Put branches from buffer into 2 nodes according to chosen partition
01790   *a_newNode = m_mem_pool.AllocNode();
01791   (*a_newNode)->m_level = a_node->m_level = level; // restore m_level
01792   LoadNodes(a_node, *a_newNode, &localVars);
01793 }
01794 
01795 double CalcRectVolumeHelper(const ON_RTreeBBox* a_rect)
01796 {
01797   double d, r;
01798 
01799   // Bounding sphere volume calculation is slower, but helps certain merge cases
01800 #if ( 3 == ON_RTree_NODE_DIM)
01801   // 3d bounding sphere volume
01802   d = (a_rect->m_max[0] - a_rect->m_min[0]);
01803   r = d * d;
01804   d = (a_rect->m_max[1] - a_rect->m_min[1]);
01805   r += d * d;
01806   d = (a_rect->m_max[2] - a_rect->m_min[2]);
01807   r += d * d;
01808   r = sqrt(r*0.5); // r = sqrt((dx^2 + dy^2 + dz^2)/2);
01809   return (r * r * r * 4.1887902047863909846168578443727); // 4/3 pi r^3
01810 #elif ( 2 == ON_RTree_NODE_DIM )
01811   // 2d bounding circle volume
01812   d = (a_rect->m_max[0] - a_rect->m_min[0]);
01813   r = d * d;
01814   d = (a_rect->m_max[1] - a_rect->m_min[1]);
01815   r += d * d;
01816   r = sqrt(r*0.5); // r = sqrt((dx^2 + dy^2)/2);
01817   return (r * r * ON_PI);
01818 #else
01819 
01820   // n-dim unit sphere volumes
01821   //  0.000000f, 2.000000f, 3.141593f, // Dimension  0,1,2
01822   //  4.188790f, 4.934802f, 5.263789f, // Dimension  3,4,5
01823   //  5.167713f, 4.724766f, 4.058712f, // Dimension  6,7,8
01824   //  3.298509f, 2.550164f, 1.884104f, // Dimension  9,10,11
01825   //  1.335263f, 0.910629f, 0.599265f, // Dimension  12,13,14
01826   //  0.381443f, 0.235331f, 0.140981f, // Dimension  15,16,17
01827   //  0.082146f, 0.046622f, 0.025807f, // Dimension  18,19,20 
01828   //return (unit_sphere_volume * radius^ON_RTree_NODE_DIM);
01829 
01830   // Faster rectangle volume calculation, but can cause poor merges
01831   d = a_rect->m_max[0] - a_rect->m_min[0];
01832   for(int i = 1; i < ON_RTree_NODE_DIM; ++i)
01833   {
01834     d *= a_rect->m_max[i] - a_rect->m_min[i];
01835   }
01836   return d;
01837 #endif
01838 }
01839 
01840 
01841 // Load branch buffer with branches from full node plus the extra branch.
01842 
01843 static void GetBranches(ON_RTreeNode* a_node, ON_RTreeBranch* a_branch, ON_RTreePartitionVars* a_parVars)
01844 {
01845   // Load the branch buffer
01846   for(int index=0; index < ON_RTree_MAX_NODE_COUNT; ++index)
01847   {
01848     a_parVars->m_branchBuf[index] = a_node->m_branch[index];
01849   }
01850   a_parVars->m_branchBuf[ON_RTree_MAX_NODE_COUNT] = *a_branch;
01851   a_parVars->m_branchCount = ON_RTree_MAX_NODE_COUNT + 1;
01852 
01853   // Calculate rect containing all in the set
01854   a_parVars->m_coverSplit = a_parVars->m_branchBuf[0].m_rect;
01855   for(int index=1; index < ON_RTree_MAX_NODE_COUNT+1; ++index)
01856   {
01857     a_parVars->m_coverSplit = CombineRectHelper(&a_parVars->m_coverSplit, &a_parVars->m_branchBuf[index].m_rect);
01858   }
01859   a_parVars->m_coverSplitArea = CalcRectVolumeHelper(&a_parVars->m_coverSplit);
01860 
01861   a_node->m_count = 0;
01862   a_node->m_level = -1;
01863 }
01864 
01865 
01866 // Method #0 for choosing a partition:
01867 // As the seeds for the two groups, pick the two rects that would waste the
01868 // most area if covered by a single rectangle, i.e. evidently the worst pair
01869 // to have in the same group.
01870 // Of the remaining, one at a time is chosen to be put in one of the two groups.
01871 // The one chosen is the one with the greatest difference in area expansion
01872 // depending on which group - the rect most strongly attracted to one group
01873 // and repelled from the other.
01874 // If one group gets too full (more would force other group to violate min
01875 // fill requirement) then other group gets the rest.
01876 // These last are the ones that can go in either group most easily.
01877 
01878 static void ChoosePartition(ON_RTreePartitionVars* a_parVars, int a_minFill)
01879 {
01880   double biggestDiff;
01881   int group, chosen, betterGroup;
01882   
01883   InitParVars(a_parVars, a_parVars->m_branchCount, a_minFill);
01884   PickSeeds(a_parVars);
01885 
01886   while (((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total)
01887        && (a_parVars->m_count[0] < (a_parVars->m_total - a_parVars->m_minFill))
01888        && (a_parVars->m_count[1] < (a_parVars->m_total - a_parVars->m_minFill)))
01889   {
01890     biggestDiff = -1.0;
01891     chosen = 0;
01892     betterGroup = 0;
01893     for(int index=0; index<a_parVars->m_total; ++index)
01894     {
01895       if(!a_parVars->m_taken[index])
01896       {
01897         ON_RTreeBBox* curRect = &a_parVars->m_branchBuf[index].m_rect;
01898         ON_RTreeBBox rect0 = CombineRectHelper(curRect, &a_parVars->m_cover[0]);
01899         ON_RTreeBBox rect1 = CombineRectHelper(curRect, &a_parVars->m_cover[1]);
01900         double growth0 = CalcRectVolumeHelper(&rect0) - a_parVars->m_area[0];
01901         double growth1 = CalcRectVolumeHelper(&rect1) - a_parVars->m_area[1];
01902         double diff = growth1 - growth0;
01903         if(diff >= 0)
01904         {
01905           group = 0;
01906         }
01907         else
01908         {
01909           group = 1;
01910           diff = -diff;
01911         }
01912 
01913         if(diff > biggestDiff)
01914         {
01915           biggestDiff = diff;
01916           chosen = index;
01917           betterGroup = group;
01918         }
01919         else if((diff == biggestDiff) && (a_parVars->m_count[group] < a_parVars->m_count[betterGroup]))
01920         {
01921           chosen = index;
01922           betterGroup = group;
01923         }
01924       }
01925     }
01926     ClassifyHelper(chosen, betterGroup, a_parVars);
01927   }
01928 
01929   // If one group too full, put remaining rects in the other
01930   if((a_parVars->m_count[0] + a_parVars->m_count[1]) < a_parVars->m_total)
01931   {
01932     if(a_parVars->m_count[0] >= a_parVars->m_total - a_parVars->m_minFill)
01933     {
01934       group = 1;
01935     }
01936     else
01937     {
01938       group = 0;
01939     }
01940     for(int index=0; index<a_parVars->m_total; ++index)
01941     {
01942       if(!a_parVars->m_taken[index])
01943       {
01944         ClassifyHelper(index, group, a_parVars);
01945       }
01946     }
01947   }
01948 }
01949 
01950 
01951 // Copy branches from the buffer into two nodes according to the partition.
01952 
01953 void ON_RTree::LoadNodes(ON_RTreeNode* a_nodeA, ON_RTreeNode* a_nodeB, ON_RTreePartitionVars* a_parVars)
01954 {
01955   for(int index=0; index < a_parVars->m_total; ++index)
01956   {
01957     if(a_parVars->m_partition[index] == 0)
01958     {
01959       AddBranch(&a_parVars->m_branchBuf[index], a_nodeA, NULL);
01960     }
01961     else if(a_parVars->m_partition[index] == 1)
01962     {
01963       AddBranch(&a_parVars->m_branchBuf[index], a_nodeB, NULL);
01964     }
01965   }
01966 }
01967 
01968 
01969 // Initialize a ON_RTreePartitionVars structure.
01970 
01971 static void InitParVars(ON_RTreePartitionVars* a_parVars, int a_maxRects, int a_minFill)
01972 {
01973   a_parVars->m_count[0] = a_parVars->m_count[1] = 0;
01974   a_parVars->m_area[0] = a_parVars->m_area[1] = (double)0;
01975   a_parVars->m_total = a_maxRects;
01976   a_parVars->m_minFill = a_minFill;
01977   for(int index=0; index < a_maxRects; ++index)
01978   {
01979     a_parVars->m_taken[index] = false;
01980     a_parVars->m_partition[index] = -1;
01981   }
01982 }
01983 
01984 
01985 
01986 static void PickSeeds(ON_RTreePartitionVars* a_parVars)
01987 {
01988   int seed0 = 0, seed1 = 1;
01989   double worst, waste;
01990   double area[ON_RTree_MAX_NODE_COUNT+1];
01991 
01992   for(int index=0; index<a_parVars->m_total; ++index)
01993   {
01994     area[index] = CalcRectVolumeHelper(&a_parVars->m_branchBuf[index].m_rect);
01995   }
01996 
01997   worst = -a_parVars->m_coverSplitArea - 1;
01998   for(int indexA=0; indexA < a_parVars->m_total-1; ++indexA)
01999   {
02000     for(int indexB = indexA+1; indexB < a_parVars->m_total; ++indexB)
02001     {
02002       ON_RTreeBBox oneRect = CombineRectHelper(&a_parVars->m_branchBuf[indexA].m_rect, &a_parVars->m_branchBuf[indexB].m_rect);
02003       waste = CalcRectVolumeHelper(&oneRect) - area[indexA] - area[indexB];
02004       if(waste > worst)
02005       {
02006         worst = waste;
02007         seed0 = indexA;
02008         seed1 = indexB;
02009       }
02010     }
02011   }
02012   ClassifyHelper(seed0, 0, a_parVars);
02013   ClassifyHelper(seed1, 1, a_parVars);
02014 }
02015 
02016 // Put a branch in one of the groups.
02017 
02018 void ClassifyHelper(int a_index, int a_group, ON_RTreePartitionVars* a_parVars)
02019 {
02020   a_parVars->m_partition[a_index] = a_group;
02021   a_parVars->m_taken[a_index] = true;
02022 
02023   if (a_parVars->m_count[a_group] == 0)
02024   {
02025     a_parVars->m_cover[a_group] = a_parVars->m_branchBuf[a_index].m_rect;
02026   }
02027   else
02028   {
02029     a_parVars->m_cover[a_group] = CombineRectHelper(&a_parVars->m_branchBuf[a_index].m_rect, &a_parVars->m_cover[a_group]);
02030   }
02031   a_parVars->m_area[a_group] = CalcRectVolumeHelper(&a_parVars->m_cover[a_group]);
02032   ++a_parVars->m_count[a_group];
02033 }
02034 
02035 
02036 // Delete a data rectangle from an index structure.
02037 // Pass in a pointer to a ON_RTreeBBox, the tid of the record, ptr to ptr to root node.
02038 // Returns 1 if record not found, 0 if success.
02039 // RemoveRect provides for eliminating the root.
02040 
02041 bool ON_RTree::RemoveRect(ON_RTreeBBox* a_rect, ON__INT_PTR a_id, ON_RTreeNode** a_root)
02042 {
02043   ON_RTreeNode* tempNode;
02044   ON_RTreeListNode* reInsertList = NULL;
02045 
02046   if(!RemoveRectRec(a_rect, a_id, *a_root, &reInsertList))
02047   {
02048     // Found and deleted a data item
02049     // Reinsert any branches from eliminated nodes
02050     while(reInsertList)
02051     {
02052       tempNode = reInsertList->m_node;
02053 
02054       for(int index = 0; index < tempNode->m_count; ++index)
02055       {
02056         InsertRect(&(tempNode->m_branch[index].m_rect),
02057                    tempNode->m_branch[index].m_id,
02058                    a_root,
02059                    tempNode->m_level);
02060       }
02061       
02062       ON_RTreeListNode* remLNode = reInsertList;
02063       reInsertList = reInsertList->m_next;
02064       
02065       m_mem_pool.FreeNode(remLNode->m_node);
02066       m_mem_pool.FreeListNode(remLNode);
02067     }
02068     
02069     // Check for redundant root (not leaf, 1 child) and eliminate
02070     if((*a_root)->m_count == 1 && (*a_root)->IsInternalNode())
02071     {
02072       tempNode = (*a_root)->m_branch[0].m_child;
02073       m_mem_pool.FreeNode(*a_root);
02074       *a_root = tempNode;
02075     }
02076     return false;
02077   }
02078   else
02079   {
02080     return true;
02081   }
02082 }
02083 
02084 
02085 // Delete a rectangle from non-root part of an index structure.
02086 // Called by RemoveRect.  Descends tree recursively,
02087 // merges branches on the way back up.
02088 // Returns 1 if record not found, 0 if success.
02089 
02090 bool ON_RTree::RemoveRectRec(ON_RTreeBBox* a_rect, ON__INT_PTR a_id, ON_RTreeNode* a_node, ON_RTreeListNode** a_listNode)
02091 {
02092   if(a_node->IsInternalNode())  // not a leaf node
02093   {
02094     for(int index = 0; index < a_node->m_count; ++index)
02095     {
02096       if(OverlapHelper(a_rect, &(a_node->m_branch[index].m_rect)))
02097       {
02098         if(!RemoveRectRec(a_rect, a_id, a_node->m_branch[index].m_child, a_listNode))
02099         {
02100           if(a_node->m_branch[index].m_child->m_count >= ON_RTree_MIN_NODE_COUNT)
02101           {
02102             // child removed, just resize parent rect
02103             a_node->m_branch[index].m_rect = NodeCover(a_node->m_branch[index].m_child);
02104           }
02105           else
02106           {
02107             // child removed, not enough entries in node, eliminate node
02108             ReInsert(a_node->m_branch[index].m_child, a_listNode);
02109             DisconnectBranch(a_node, index); // Must return after this call as count has changed
02110           }
02111           return false;
02112         }
02113       }
02114     }
02115     return true;
02116   }
02117   else // A leaf node
02118   {
02119     for(int index = 0; index < a_node->m_count; ++index)
02120     {
02121       if(a_node->m_branch[index].m_id == a_id)
02122       {
02123         DisconnectBranch(a_node, index); // Must return after this call as count has changed
02124         return false;
02125       }
02126     }
02127     return true;
02128   }
02129 }
02130 
02131 
02132 // Decide whether two rectangles overlap.
02133 bool OverlapHelper(const ON_RTreeBBox* a_rectA, const ON_RTreeBBox* a_rectB)
02134 {
02135   const double* mn;
02136   const double* mx;
02137 
02138   mx = a_rectA->m_max;
02139   mn = a_rectB->m_min;
02140   if ( *mx++ < *mn++ ) return false;
02141   if ( *mx++ < *mn++ ) return false;
02142   if ( *mx   < *mn   ) return false;
02143 
02144   mx = a_rectB->m_max;
02145   mn = a_rectA->m_min;
02146   if ( *mx++ < *mn++ ) return false;
02147   if ( *mx++ < *mn++ ) return false;
02148   if ( *mx   < *mn   ) return false;
02149 
02150   return true;
02151 }
02152 
02153 //static bool OverlapHelper(const struct ON_RTreeSphere* a_sphere, const ON_RTreeBBox* a_rect)
02154 //{
02155 //  double d[3], t, r;
02156 //  const double* mn;
02157 //  const double* mx;
02158 //  const double* pt;
02159 //
02160 //  pt = a_sphere->m_point;
02161 //  r  = a_sphere->m_radius;
02162 //  mn = a_rect->m_min;
02163 //  mx = a_rect->m_max;
02164 //
02165 //  if ( *pt < *mn )
02166 //  {
02167 //    d[0] = *mn - *pt;
02168 //    if ( d[0] > r )
02169 //      return false;
02170 //  }
02171 //  else if ( *pt > *mx )
02172 //  {
02173 //    d[0] = *pt - *mx;
02174 //    if ( d[0] > r )
02175 //      return false;
02176 //  }
02177 //  else
02178 //  {
02179 //    d[0] = 0.0;
02180 //  }
02181 //
02182 //  mn++;
02183 //  mx++;
02184 //  pt++;
02185 //  if ( *pt < *mn )
02186 //  {
02187 //    d[1] = *mn - *pt;
02188 //    if ( d[1] > r )
02189 //      return false;
02190 //    if ( d[1] > d[0] )
02191 //    {
02192 //      t = d[0]; d[0] = d[1]; d[1] = t;
02193 //    }
02194 //  }
02195 //  else if ( *pt > *mx )
02196 //  {
02197 //    d[1] = *pt - *mx;
02198 //    if ( d[1] > r )
02199 //      return false;
02200 //    if ( d[1] > d[0] )
02201 //    {
02202 //      t = d[0]; d[0] = d[1]; d[1] = t;
02203 //    }
02204 //  }
02205 //  else
02206 //  {
02207 //    d[1] = 0.0;
02208 //  }
02209 //
02210 //  mn++;
02211 //  mx++;
02212 //  pt++;
02213 //  if ( *pt < *mn )
02214 //  {
02215 //    d[2] = *mn - *pt;
02216 //    if ( d[2] > r )
02217 //      return false;
02218 //    if ( d[2] > d[0] )
02219 //    {
02220 //      t = d[0]; d[0] = d[2]; d[2] = t;
02221 //    }
02222 //  }
02223 //  else if ( *pt > *mx )
02224 //  {
02225 //    d[2] = *pt - *mx;
02226 //    if ( d[2] > r )
02227 //      return false;
02228 //    if ( d[2] > d[0] )
02229 //    {
02230 //      t = d[0]; d[0] = d[2]; d[2] = t;
02231 //    }
02232 //  }
02233 //  else
02234 //  {
02235 //    d[2] = 0.0;
02236 //  }
02237 //
02238 //  if ( d[0] > 0.0 )
02239 //  {
02240 //    d[1] /= d[0];
02241 //    d[2] /= d[0];
02242 //    d[0] *= sqrt(1.0 + d[1]*d[1] + d[2]*d[2]);
02243 //    return (d[0] <= r);
02244 //  }
02245 //
02246 //  return true;
02247 //}
02248 
02249 static double DistanceToBoxHelper( 
02250   const double* pt, 
02251   double r, 
02252   const ON_RTreeBBox* a_rect
02253   )
02254 {
02255   // If the sphere with center at pt and radius r intersects a_rect, then
02256   // the distance from pt to a_rect is returned. A value of 0.0 indicates
02257   // that pt is inside a_rect.  If the distance from pt to a_rect is
02258   // greater than r, then some number > r and <= actual distance from
02259   // pt to a_rect is returned as quickly as possible.
02260 
02261   double d[3], t;
02262   const double* mn;
02263   const double* mx;
02264 
02265   mn = a_rect->m_min;
02266   mx = a_rect->m_max;
02267 
02268   if ( *pt < *mn )
02269   {
02270     d[0] = *mn - *pt;
02271     if ( d[0] > r )
02272       return d[0];
02273   }
02274   else if ( *pt > *mx )
02275   {
02276     d[0] = *pt - *mx;
02277     if ( d[0] > r )
02278       return d[0];
02279   }
02280   else
02281   {
02282     d[0] = 0.0;
02283   }
02284 
02285   mn++;
02286   mx++;
02287   pt++;
02288   if ( *pt < *mn )
02289   {
02290     d[1] = *mn - *pt;
02291     if ( d[1] > r )
02292       return d[1];
02293     if ( d[1] > d[0] )
02294     {
02295       t = d[0]; d[0] = d[1]; d[1] = t;
02296     }
02297   }
02298   else if ( *pt > *mx )
02299   {
02300     d[1] = *pt - *mx;
02301     if ( d[1] > r )
02302       return d[1];
02303     if ( d[1] > d[0] )
02304     {
02305       t = d[0]; d[0] = d[1]; d[1] = t;
02306     }
02307   }
02308   else
02309   {
02310     d[1] = 0.0;
02311   }
02312 
02313   mn++;
02314   mx++;
02315   pt++;
02316   if ( *pt < *mn )
02317   {
02318     d[2] = *mn - *pt;
02319     if ( d[2] > r )
02320       return d[2];
02321     if ( d[2] > d[0] )
02322     {
02323       t = d[0]; d[0] = d[2]; d[2] = t;
02324     }
02325   }
02326   else if ( *pt > *mx )
02327   {
02328     d[2] = *pt - *mx;
02329     if ( d[2] > r )
02330       return d[2];
02331     if ( d[2] > d[0] )
02332     {
02333       t = d[0]; d[0] = d[2]; d[2] = t;
02334     }
02335   }
02336   else
02337   {
02338     d[2] = 0.0;
02339   }
02340 
02341   if ( d[0] > 0.0 )
02342   {
02343     d[1] /= d[0];
02344     d[2] /= d[0];
02345     d[0] *= sqrt(1.0 + d[1]*d[1] + d[2]*d[2]);
02346   }
02347 
02348   return d[0];
02349 }
02350 
02351 static double DistanceToCapsuleAxisHelper(const struct ON_RTreeCapsule* a_capsule, const ON_RTreeBBox* a_rect)
02352 {
02353   double L[2][3], s[2];
02354   if ( 0.0 == a_capsule->m_domain[0] && 1.0 == a_capsule->m_domain[1] )
02355     return ((const ON_BoundingBox*)a_rect->m_min)->MinimumDistanceTo( *((const ON_Line*)a_capsule->m_point[0]) );
02356 
02357   if ( 0.0 == a_capsule->m_domain[0] )
02358   {
02359     L[0][0] = a_capsule->m_point[0][0];
02360     L[0][1] = a_capsule->m_point[0][1];
02361     L[0][2] = a_capsule->m_point[0][2];
02362   }
02363   else
02364   {
02365     s[0] = 1.0 - a_capsule->m_domain[0];
02366     s[1] = a_capsule->m_domain[0];
02367     L[0][0] = s[0]*a_capsule->m_point[0][0] + s[1]*a_capsule->m_point[1][0];
02368     L[0][1] = s[0]*a_capsule->m_point[0][1] + s[1]*a_capsule->m_point[1][1];
02369     L[0][2] = s[0]*a_capsule->m_point[0][2] + s[1]*a_capsule->m_point[1][2];
02370   }
02371 
02372   if ( 0.0 == a_capsule->m_domain[1] )
02373   {
02374     L[1][0] = a_capsule->m_point[1][0];
02375     L[1][1] = a_capsule->m_point[1][1];
02376     L[1][2] = a_capsule->m_point[1][2];
02377   }
02378   else
02379   {
02380     s[0] = 1.0 - a_capsule->m_domain[1];
02381     s[1] = a_capsule->m_domain[1];
02382     L[1][0] = s[0]*a_capsule->m_point[0][0] + s[1]*a_capsule->m_point[1][0];
02383     L[1][1] = s[0]*a_capsule->m_point[0][1] + s[1]*a_capsule->m_point[1][1];
02384     L[1][2] = s[0]*a_capsule->m_point[0][2] + s[1]*a_capsule->m_point[1][2];
02385   }
02386 
02387   return ((const ON_BoundingBox*)a_rect->m_min)->MinimumDistanceTo( *((const ON_Line*)L[0]) );
02388 }
02389 
02390 // Add a node to the reinsertion list.  All its branches will later
02391 // be reinserted into the index structure.
02392 
02393 void ON_RTree::ReInsert(ON_RTreeNode* a_node, ON_RTreeListNode** a_listNode)
02394 {
02395   ON_RTreeListNode* newListNode;
02396 
02397   newListNode = m_mem_pool.AllocListNode();
02398   newListNode->m_node = a_node;
02399   newListNode->m_next = *a_listNode;
02400   *a_listNode = newListNode;
02401 }
02402 
02403 
02404 static
02405 bool OverlapBoundedPlaneXYZHelper( const double* a_bounded_plane, const ON_RTreeBBox* a_rect )
02406 {
02407   unsigned char flag = 0;
02408   double x, y, z, d, v;
02409   
02410   // check the 8 corners of the box minimizing the number of evaluations
02411   // and unrolling the loop for speed
02412 
02413   // corner = (min, min, min)
02414   x = a_bounded_plane[0]*a_rect->m_min[0]; 
02415   y = a_bounded_plane[1]*a_rect->m_min[1];
02416   z = a_bounded_plane[2]*a_rect->m_min[2];
02417   d = a_bounded_plane[3];
02418   v = x + y + z + d;
02419   if ( v < a_bounded_plane[4] )
02420     flag = 1;
02421   else if ( v > a_bounded_plane[5] )
02422     flag = 2;
02423   else
02424     return true;
02425   
02426   // corner = (max, min, min)
02427   x = a_bounded_plane[0]*a_rect->m_max[0]; 
02428   v = x + y + z + d;
02429   if ( v < a_bounded_plane[4] )
02430   {
02431     flag |= 1;
02432     if ( 3 == flag )
02433       return true;
02434   }
02435   else if ( v > a_bounded_plane[5] )
02436   {
02437     flag |= 2;
02438     if ( 3 == flag )
02439       return true;
02440   }
02441   else
02442     return true;
02443   
02444   // corner = (max, max, min)
02445   y = a_bounded_plane[1]*a_rect->m_max[1];
02446   v = x + y + z + d;
02447   if ( v < a_bounded_plane[4] )
02448   {
02449     flag |= 1;
02450     if ( 3 == flag )
02451       return true;
02452   }
02453   else if ( v > a_bounded_plane[5] )
02454   {
02455     flag |= 2;
02456     if ( 3 == flag )
02457       return true;
02458   }
02459   else
02460     return true;
02461     
02462   // corner = (max, max, max)
02463   z = a_bounded_plane[2]*a_rect->m_max[2];
02464   v = x + y + z + d;
02465   if ( v < a_bounded_plane[4] )
02466   {
02467     flag |= 1;
02468     if ( 3 == flag )
02469       return true;
02470   }
02471   else if ( v > a_bounded_plane[5] )
02472   {
02473     flag |= 2;
02474     if ( 3 == flag )
02475       return true;
02476   }
02477   else
02478     return true;
02479 
02480   // corner = (min, max, max)
02481   x = a_bounded_plane[0]*a_rect->m_min[0]; 
02482   v = x + y + z + d;
02483   if ( v < a_bounded_plane[4] )
02484   {
02485     flag |= 1;
02486     if ( 3 == flag )
02487       return true;
02488   }
02489   else if ( v > a_bounded_plane[5] )
02490   {
02491     flag |= 2;
02492     if ( 3 == flag )
02493       return true;
02494   }
02495   else
02496     return true;
02497   
02498   // corner = (min, min, max)
02499   y = a_bounded_plane[1]*a_rect->m_min[1];
02500   v = x + y + z + d;
02501   if ( v < a_bounded_plane[4] )
02502   {
02503     flag |= 1;
02504     if ( 3 == flag )
02505       return true;
02506   }
02507   else if ( v > a_bounded_plane[5] )
02508   {
02509     flag |= 2;
02510     if ( 3 == flag )
02511       return true;
02512   }
02513   else
02514     return true;
02515 
02516   // corner = (max, min, max)
02517   x = a_bounded_plane[0]*a_rect->m_max[0]; 
02518   v = x + y + z + d;
02519   if ( v < a_bounded_plane[4] )
02520   {
02521     flag |= 1;
02522     if ( 3 == flag )
02523       return true;
02524   }
02525   else if ( v > a_bounded_plane[5] )
02526   {
02527     flag |= 2;
02528     if ( 3 == flag )
02529       return true;
02530   }
02531   else
02532     return true;
02533 
02534   // corner = (min, max, min)
02535   x = a_bounded_plane[0]*a_rect->m_min[0]; 
02536   y = a_bounded_plane[1]*a_rect->m_max[1];
02537   z = a_bounded_plane[2]*a_rect->m_min[2];
02538   v = x + y + z + d;
02539   if ( v < a_bounded_plane[4] )
02540   {
02541     flag |= 1;
02542     if ( 3 == flag )
02543       return true;
02544   }
02545   else if ( v > a_bounded_plane[5] )
02546   {
02547     flag |= 2;
02548     if ( 3 == flag )
02549       return true;
02550   }
02551   else
02552     return true;
02553   
02554   // Either all 8 box corners 
02555   // are below the min plane (flag=1)
02556   // or above the max plane (flag=2).
02557   return false;
02558 }
02559 
02560 static
02561 bool SearchBoundedPlaneXYZHelper(const ON_RTreeNode* a_node, const double* a_bounded_plane, ON_RTreeSearchResultCallback& a_result ) 
02562 {
02563   int i, count;
02564 
02565   if ( (count = a_node->m_count) > 0 )
02566   {
02567     const ON_RTreeBranch* branch = a_node->m_branch;
02568     if(a_node->IsInternalNode()) 
02569     {
02570       // a_node is an internal node - search m_branch[].m_child as needed
02571       for( i=0; i < count; ++i )
02572       {
02573         if(OverlapBoundedPlaneXYZHelper(a_bounded_plane, &branch[i].m_rect))
02574         {
02575           if(!SearchBoundedPlaneXYZHelper(branch[i].m_child, a_bounded_plane, a_result) )
02576           {
02577             return false; // Don't continue searching
02578           }
02579         }
02580       }
02581     }
02582     else
02583     {
02584       // a_node is a leaf node - return m_branch[].m_id values
02585       for(i=0; i < count; ++i)
02586       {
02587         if(OverlapBoundedPlaneXYZHelper(a_bounded_plane, &branch[i].m_rect))
02588         {
02589           if ( !a_result.m_resultCallback( a_result.m_context, branch[i].m_id ) )
02590           {
02591             // callback canceled search
02592             return false;
02593           }
02594         }
02595       }
02596     }
02597   }
02598 
02599   return true; // Continue searching
02600 }
02601 
02602 bool ON_RTree::Search(
02603   const double a_plane_eqn[4],
02604   double a_min,
02605   double a_max,
02606   bool ON_MSC_CDECL a_resultCallback(void* a_context, ON__INT_PTR a_id), 
02607   void* a_context
02608   ) const
02609 {
02610   if (    0 == m_root 
02611        || 0 == a_plane_eqn 
02612        || !(a_min <= a_max) 
02613        || (0.0 == a_plane_eqn[0] && 0.0 == a_plane_eqn[1] && 0.0 == a_plane_eqn[2])
02614      )
02615     return false;
02616 
02617   double bounded_plane[6];
02618   bounded_plane[0] = a_plane_eqn[0];
02619   bounded_plane[1] = a_plane_eqn[1];
02620   bounded_plane[2] = a_plane_eqn[2];
02621   bounded_plane[3] = a_plane_eqn[3];
02622   bounded_plane[4] = a_min;
02623   bounded_plane[5] = a_max;
02624 
02625   ON_RTreeSearchResultCallback result;
02626   result.m_context = a_context;
02627   result.m_resultCallback = a_resultCallback;
02628 
02629   return SearchBoundedPlaneXYZHelper(m_root, bounded_plane, result);
02630 }
02631 
02632 // Search in an index tree or subtree for all data retangles that overlap the argument rectangle.
02633 
02634 static
02635 bool SearchHelper(const ON_RTreeNode* a_node, ON_RTreeBBox* a_rect, ON_RTreeSearchResultCallback& a_result ) 
02636 {
02637   // NOTE: 
02638   //  Some versions of ON_RTree::Search shrink a_rect as the search progresses.
02639   int i, count;
02640 
02641   if ( (count = a_node->m_count) > 0 )
02642   {
02643     const ON_RTreeBranch* branch = a_node->m_branch;
02644     if(a_node->IsInternalNode()) 
02645     {
02646       // a_node is an internal node - search m_branch[].m_child as needed
02647       for( i=0; i < count; ++i )
02648       {
02649         if(OverlapHelper(a_rect, &branch[i].m_rect))
02650         {
02651           if(!SearchHelper(branch[i].m_child, a_rect, a_result) )
02652           {
02653             return false; // Don't continue searching
02654           }
02655         }
02656       }
02657     }
02658     else
02659     {
02660       // a_node is a leaf node - return m_branch[].m_id values
02661       for(i=0; i < count; ++i)
02662       {
02663         if(OverlapHelper(a_rect, &branch[i].m_rect))
02664         {
02665           if ( !a_result.m_resultCallback( a_result.m_context, branch[i].m_id ) )
02666           {
02667             // callback canceled search
02668             return false;
02669           }
02670         }
02671       }
02672     }
02673   }
02674 
02675   return true; // Continue searching
02676 }
02677 
02678 static
02679 bool SearchHelper(const ON_RTreeNode* a_node, struct ON_RTreeSphere* a_sphere, ON_RTreeSearchResultCallback& a_result ) 
02680 {
02681   // NOTE: 
02682   //  Some versions of ON_RTree::Search shrink a_sphere as the search progresses.
02683   int i, closest_i, count;
02684   const double* sphere_center;
02685   const ON_RTreeBranch* branch;
02686   double r[ON_RTree_MAX_NODE_COUNT], sphere_radius, closest_d;
02687   
02688   if ( (count = a_node->m_count) > 0 )
02689   {
02690     branch = a_node->m_branch;
02691     sphere_center = a_sphere->m_point;
02692     sphere_radius = a_sphere->m_radius;
02693     closest_i = -1;
02694     closest_d = sphere_radius;
02695 
02696     for( i = 0; i < count; ++i )
02697     {
02698       // The radius parameter passed to DistanceToBoxHelper()
02699       // needs to be sphere_radius and not closest_d in order 
02700       // for the for() loops below to work correctly.
02701       r[i] = DistanceToBoxHelper( sphere_center, sphere_radius, &branch[i].m_rect );
02702       if ( r[i] <= closest_d )
02703       {
02704         closest_d = r[i];
02705         closest_i = i;
02706       }
02707     }
02708 
02709     // If all of the branches rectangles do not intersect the sphere,
02710     // then closest_i = -1.
02711     if ( closest_i >= 0 )
02712     {
02713       if(a_node->IsInternalNode()) 
02714       {
02715         // a_node is an internal node - search m_branch[].m_child as needed.
02716         // Search a closer node first to avoid worst case search times
02717         // in calculations where the calls to a_result.m_resultCallback()
02718         // reduce a_sphere->m_radius as results are found.  Closest point
02719         // calculations are an example.
02720         if ( !SearchHelper(branch[closest_i].m_child, a_sphere, a_result) )
02721         {
02722           // callback canceled search
02723           return false;
02724         }
02725 
02726         for( i = 0; i < count; ++i )
02727         {
02728           // Note that the calls to SearchHelper() can reduce the
02729           // value of a_sphere->m_radius.
02730           if( i != closest_i && r[i] <= a_sphere->m_radius )
02731           {
02732             if(!SearchHelper(branch[i].m_child, a_sphere, a_result) )
02733             {
02734               return false; // Don't continue searching
02735             }
02736           }
02737         }
02738       }
02739       else
02740       {
02741         // a_node is a leaf node - return m_branch[].m_id values
02742         // Search a closer node first to avoid worst case search times
02743         // in calculations where the calls to a_result.m_resultCallback()
02744         // reduce a_sphere->m_radius as results are found.  Closest point
02745         // calculations are an example.
02746         if ( !a_result.m_resultCallback( a_result.m_context, branch[closest_i].m_id ) ) 
02747         {
02748           // callback canceled search
02749           return false;
02750         }
02751 
02752         for( i = 0; i < count; ++i)
02753         {
02754           // Note that the calls to a_result.m_resultCallback() can reduce
02755           // the value of a_sphere->m_radius.
02756           if( i != closest_i && r[i] <= a_sphere->m_radius )
02757           {
02758             if ( !a_result.m_resultCallback( a_result.m_context, branch[i].m_id ) )
02759             {
02760               // callback canceled search
02761               return false;
02762             }
02763           }
02764         }
02765       }
02766     }
02767   }
02768 
02769   return true; // Continue searching
02770 }
02771 
02772 
02773 static
02774 bool SearchHelper(const ON_RTreeNode* a_node, struct ON_RTreeCapsule* a_capsule, ON_RTreeSearchResultCallback& a_result ) 
02775 {
02776   // NOTE: 
02777   //  Some versions of ON_RTree::Search shrink a_sphere as the search progresses.
02778   int i, count;
02779   double r[2];
02780   
02781   if ( (count = a_node->m_count) > 0 )
02782   {
02783     const ON_RTreeBranch* branch = a_node->m_branch;
02784     if(a_node->IsInternalNode()) 
02785     {
02786       // a_node is an internal node - search m_branch[].m_child as needed
02787       if ( count > 1 )
02788       {
02789         // search a closer node first to avoid worst case search times
02790         // in closest point style calculations
02791         r[0] = DistanceToCapsuleAxisHelper( a_capsule, &branch[0].m_rect );
02792         r[1] = DistanceToCapsuleAxisHelper( a_capsule, &branch[count-1].m_rect );
02793         i = ( r[0] <= r[1] ) ? 0 : count-1;
02794         if (    (r[i?1:0] <= a_capsule->m_radius && !SearchHelper(branch[i].m_child, a_capsule, a_result)) 
02795              || (r[i?0:1] <= a_capsule->m_radius && !SearchHelper(branch[count-1 - i].m_child, a_capsule, a_result))
02796           )
02797         {
02798           // callback canceled search
02799           return false;
02800         }
02801         count -= 2;
02802         branch++;
02803       }
02804 
02805       r[1] = a_capsule->m_radius;
02806       for( i = 0; i < count; ++i )
02807       {
02808         r[0] = DistanceToCapsuleAxisHelper(a_capsule, &branch[i].m_rect);
02809         if(r[0] <= r[1])
02810         {
02811           if(!SearchHelper(branch[i].m_child, a_capsule, a_result) )
02812           {
02813             return false; // Don't continue searching
02814           }
02815           // a_result.m_resultCallback can shrink the capsule
02816           r[1] = a_capsule->m_radius;
02817         }
02818       }
02819     }
02820     else
02821     {
02822       // a_node is a leaf node - return m_branch[].m_id values
02823       if ( count > 1 )
02824       {
02825         // search a closer node first to avoid worst case search times
02826         // in closest point style calculations
02827         r[0] = DistanceToCapsuleAxisHelper( a_capsule, &branch[0].m_rect );
02828         r[1] = DistanceToCapsuleAxisHelper( a_capsule, &branch[count-1].m_rect );
02829         i = ( r[0] <= r[1] ) ? 0 : count-1;
02830         if (    (r[i?1:0] <= a_capsule->m_radius && !a_result.m_resultCallback( a_result.m_context, branch[i].m_id )) 
02831              || (r[i?0:1] <= a_capsule->m_radius && !a_result.m_resultCallback( a_result.m_context, branch[count-1 - i].m_id ))
02832           )
02833         {
02834           // callback canceled search
02835           return false;
02836         }
02837         count -= 2;
02838         branch++;
02839       }
02840 
02841       r[1] = a_capsule->m_radius;
02842       for( i = 0; i < count; ++i)
02843       {
02844         r[0] = DistanceToCapsuleAxisHelper(a_capsule, &branch[i].m_rect);
02845         if(r[0] <= r[1])
02846         {
02847           if ( !a_result.m_resultCallback( a_result.m_context, branch[i].m_id ) )
02848           {
02849             // callback canceled search
02850             return false;
02851           }
02852           // a_result.m_resultCallback can shrink the capsule
02853           r[1] = a_capsule->m_radius;
02854         }
02855       }
02856     }
02857   }
02858 
02859   return true; // Continue searching
02860 }
02861 
02862 
02863 // Search in an index tree or subtree for all data retangles that overlap the argument rectangle.
02864 
02865 static bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_SimpleArray<ON_RTreeLeaf> &a_result)
02866 {
02867   int i, count;
02868 
02869   if ( (count = a_node->m_count) > 0 )
02870   {
02871     const ON_RTreeBranch* branch = a_node->m_branch;
02872     if(a_node->IsInternalNode()) 
02873     {
02874       // a_node is an internal node - search m_branch[].m_child as needed
02875       for( i=0; i < count; ++i )
02876       {
02877         if(OverlapHelper(a_rect, &branch[i].m_rect))
02878         {
02879           if(!SearchHelper(branch[i].m_child, a_rect, a_result) )
02880           {
02881             return false; // Don't continue searching
02882           }
02883         }
02884       }
02885     }
02886     else
02887     {
02888       // a_node is a leaf node - return point to the branch
02889       for(i=0; i < count; ++i)
02890       {
02891         if(OverlapHelper(a_rect, &branch[i].m_rect))
02892         {
02893           ON_RTreeLeaf& leaf = a_result.AppendNew();
02894           leaf.m_rect = branch[i].m_rect;
02895           leaf.m_id = branch[i].m_id;
02896         }
02897       }
02898     }
02899   }
02900 
02901   return true; // Continue searching
02902 }
02903 
02904 
02905 static
02906 bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_SimpleArray<void*> &a_result)
02907 {
02908   int i, count;
02909 
02910   if ( (count = a_node->m_count) > 0 )
02911   {
02912     const ON_RTreeBranch* branch = a_node->m_branch;
02913     if(a_node->IsInternalNode()) 
02914     {
02915       // a_node is an internal node - search m_branch[].m_child as needed
02916       for( i=0; i < count; ++i )
02917       {
02918         if(OverlapHelper(a_rect, &branch[i].m_rect))
02919         {
02920           if(!SearchHelper(branch[i].m_child, a_rect, a_result) )
02921           {
02922             return false; // Don't continue searching
02923           }
02924         }
02925       }
02926     }
02927     else
02928     {
02929       // a_node is a leaf node - return m_branch[].m_id values
02930       for(i=0; i < count; ++i)
02931       {
02932         if(OverlapHelper(a_rect, &branch[i].m_rect))
02933         {
02934           // The (void*) cast is safe because branch[i].m_id is an ON__INT_PTR
02935 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
02936 #pragma warning( push )
02937 // Disable warning C4312: 'type cast' : conversion from 'const ON__INT_PTR' to 'void *' of greater size
02938 #pragma warning( disable : 4312 )
02939 #endif
02940           a_result.Append( (void*)branch[i].m_id );
02941 #if defined(ON_COMPILER_MSC) && 4 == ON_SIZEOF_POINTER
02942 #pragma warning( pop )
02943 #endif
02944         }
02945       }
02946     }
02947   }
02948 
02949   return true; // Continue searching
02950 }
02951 
02952 static
02953 bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_SimpleArray<int> &a_result)
02954 {
02955   int i, count;
02956 
02957   if ( (count = a_node->m_count) > 0 )
02958   {
02959     const ON_RTreeBranch* branch = a_node->m_branch;
02960     if(a_node->IsInternalNode()) 
02961     {
02962       // a_node is an internal node - search m_branch[].m_child as needed
02963       for( i=0; i < count; ++i )
02964       {
02965         if(OverlapHelper(a_rect, &branch[i].m_rect))
02966         {
02967           if(!SearchHelper(branch[i].m_child, a_rect, a_result) )
02968           {
02969             return false; // Don't continue searching
02970           }
02971         }
02972       }
02973     }
02974     else
02975     {
02976       // a_node is a leaf node - return m_branch[].m_id values
02977       for(i=0; i < count; ++i)
02978       {
02979         if(OverlapHelper(a_rect, &branch[i].m_rect))
02980         {
02981           a_result.Append( (int)branch[i].m_id );
02982         }
02983       }
02984     }
02985   }
02986 
02987   return true; // Continue searching
02988 }
02989 
02990 static
02991 bool SearchHelper(const ON_RTreeNode* a_node, const ON_RTreeBBox* a_rect, ON_RTreeSearchResult& a_result)
02992 {
02993   int i, count;
02994 
02995   if ( (count = a_node->m_count) > 0 )
02996   {
02997     const ON_RTreeBranch* branch = a_node->m_branch;
02998     if(a_node->IsInternalNode()) 
02999     {
03000       // a_node is an internal node - search m_branch[].m_child as needed
03001       for( i=0; i < count; ++i )
03002       {
03003         if(OverlapHelper(a_rect, &branch[i].m_rect))
03004         {
03005           if(!SearchHelper(branch[i].m_child, a_rect, a_result) )
03006           {
03007             return false; // Don't continue searching
03008           }
03009         }
03010       }
03011     }
03012     else
03013     {
03014       // a_node is a leaf node - return m_branch[].m_id values
03015       for(i=0; i < count; ++i)
03016       {
03017         if(OverlapHelper(a_rect, &branch[i].m_rect))
03018         {
03019           if ( a_result.m_count >= a_result.m_capacity )
03020             return false; // No more space for results
03021           a_result.m_id[a_result.m_count++] = branch[i].m_id;
03022         }
03023       }
03024     }
03025   }
03026 
03027   return true; // Continue searching
03028 }
03029 


pcl
Author(s): Open Perception
autogenerated on Wed Aug 26 2015 15:27:03