MOAB: Mesh Oriented datABase  (version 5.4.0)
BVHTree.cpp
Go to the documentation of this file.
00001 #include "moab/BVHTree.hpp"
00002 #include "moab/Interface.hpp"
00003 #include "moab/ElemEvaluator.hpp"
00004 #include "moab/ReadUtilIface.hpp"
00005 #include "moab/CpuTimer.hpp"
00006 
00007 namespace moab
00008 {
00009 const char* BVHTree::treeName = "BVHTree";
00010 
00011 ErrorCode BVHTree::build_tree( const Range& entities, EntityHandle* tree_root_set, FileOptions* options )
00012 {
00013     ErrorCode rval;
00014     CpuTimer cp;
00015 
00016     if( options )
00017     {
00018         rval = parse_options( *options );
00019         if( MB_SUCCESS != rval ) return rval;
00020 
00021         if( !options->all_seen() ) return MB_FAILURE;
00022     }
00023 
00024     // calculate bounding box of elements
00025     BoundBox box;
00026     rval = box.update( *moab(), entities );
00027     if( MB_SUCCESS != rval ) return rval;
00028 
00029     // create tree root
00030     EntityHandle tmp_root;
00031     if( !tree_root_set ) tree_root_set = &tmp_root;
00032     rval = create_root( box.bMin.array(), box.bMax.array(), *tree_root_set );
00033     if( MB_SUCCESS != rval ) return rval;
00034     rval = mbImpl->add_entities( *tree_root_set, entities );
00035     if( MB_SUCCESS != rval ) return rval;
00036 
00037     // a fully balanced tree will have 2*_entities.size()
00038     // which is one doubling away..
00039     std::vector< Node > tree_nodes;
00040     tree_nodes.reserve( entities.size() / maxPerLeaf );
00041     std::vector< HandleData > handle_data_vec;
00042     rval = construct_element_vec( handle_data_vec, entities, boundBox );
00043     if( MB_SUCCESS != rval ) return rval;
00044 
00045 #ifndef NDEBUG
00046     for( std::vector< HandleData >::const_iterator i = handle_data_vec.begin(); i != handle_data_vec.end(); ++i )
00047     {
00048         if( !boundBox.intersects_box( i->myBox, 0 ) )
00049         {
00050             std::cerr << "BB:" << boundBox << "EB:" << i->myBox << std::endl;
00051             return MB_FAILURE;
00052         }
00053     }
00054 #endif
00055     // We only build nonempty trees
00056     if( !handle_data_vec.empty() )
00057     {
00058         // initially all bits are set
00059         tree_nodes.push_back( Node() );
00060         const int depth = local_build_tree( tree_nodes, handle_data_vec.begin(), handle_data_vec.end(), 0, boundBox );
00061 #ifndef NDEBUG
00062         std::set< EntityHandle > entity_handles;
00063         for( std::vector< Node >::iterator n = tree_nodes.begin(); n != tree_nodes.end(); ++n )
00064         {
00065             for( HandleDataVec::const_iterator j = n->entities.begin(); j != n->entities.end(); ++j )
00066             {
00067                 entity_handles.insert( j->myHandle );
00068             }
00069         }
00070         if( entity_handles.size() != entities.size() )
00071         {
00072             std::cout << "Entity Handle Size Mismatch!" << std::endl;
00073         }
00074         for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
00075         {
00076             if( entity_handles.find( *i ) == entity_handles.end() )
00077                 std::cout << "Tree is missing an entity! " << std::endl;
00078         }
00079 #endif
00080         treeDepth = std::max( depth, treeDepth );
00081     }
00082 
00083     // convert vector of Node's to entity sets and vector of TreeNode's
00084     rval = convert_tree( tree_nodes );
00085     if( MB_SUCCESS != rval ) return rval;
00086 
00087     treeStats.reset();
00088     rval               = treeStats.compute_stats( mbImpl, startSetHandle );
00089     treeStats.initTime = cp.time_elapsed();
00090 
00091     return rval;
00092 }
00093 
00094 ErrorCode BVHTree::convert_tree( std::vector< Node >& tree_nodes )
00095 {
00096     // first construct the proper number of entity sets
00097     ReadUtilIface* read_util;
00098     ErrorCode rval = mbImpl->query_interface( read_util );
00099     if( MB_SUCCESS != rval ) return rval;
00100 
00101     {  // isolate potentially-large std::vector so it gets deleted earlier
00102         std::vector< unsigned int > tmp_flags( tree_nodes.size(), meshsetFlags );
00103         rval = read_util->create_entity_sets( tree_nodes.size(), &tmp_flags[0], 0, startSetHandle );
00104         if( MB_SUCCESS != rval ) return rval;
00105         rval = mbImpl->release_interface( read_util );
00106         if( MB_SUCCESS != rval ) return rval;
00107     }
00108 
00109     // populate the sets and the TreeNode vector
00110     EntityHandle set_handle = startSetHandle;
00111     std::vector< Node >::iterator it;
00112     myTree.reserve( tree_nodes.size() );
00113     for( it = tree_nodes.begin(); it != tree_nodes.end(); ++it, set_handle++ )
00114     {
00115         if( it != tree_nodes.begin() && !it->entities.empty() )
00116         {
00117             Range range;
00118             for( HandleDataVec::iterator hit = it->entities.begin(); hit != it->entities.end(); ++hit )
00119                 range.insert( hit->myHandle );
00120             rval = mbImpl->add_entities( set_handle, range );
00121             if( MB_SUCCESS != rval ) return rval;
00122         }
00123         myTree.push_back( TreeNode( it->dim, it->child, it->Lmax, it->Rmin, it->box ) );
00124 
00125         if( it->dim != 3 )
00126         {
00127             rval = mbImpl->add_child_meshset( set_handle, startSetHandle + it->child );
00128             if( MB_SUCCESS != rval ) return rval;
00129             rval = mbImpl->add_child_meshset( set_handle, startSetHandle + it->child + 1 );
00130             if( MB_SUCCESS != rval ) return rval;
00131         }
00132     }
00133 
00134     return MB_SUCCESS;
00135 }
00136 
00137 ErrorCode BVHTree::parse_options( FileOptions& opts )
00138 {
00139     ErrorCode rval = parse_common_options( opts );
00140     if( MB_SUCCESS != rval ) return rval;
00141 
00142     //  SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
00143     int tmp_int;
00144     rval = opts.get_int_option( "SPLITS_PER_DIR", tmp_int );
00145     if( MB_SUCCESS == rval ) splitsPerDir = tmp_int;
00146 
00147     return MB_SUCCESS;
00148 }
00149 
00150 void BVHTree::establish_buckets( HandleDataVec::const_iterator begin,
00151                                  HandleDataVec::const_iterator end,
00152                                  const BoundBox& interval,
00153                                  std::vector< std::vector< Bucket > >& buckets ) const
00154 {
00155     // put each element into its bucket
00156     for( HandleDataVec::const_iterator i = begin; i != end; ++i )
00157     {
00158         const BoundBox& box = i->myBox;
00159         for( unsigned int dim = 0; dim < 3; ++dim )
00160         {
00161             const unsigned int index = Bucket::bucket_index( splitsPerDir, box, interval, dim );
00162             assert( index < buckets[dim].size() );
00163             Bucket& bucket = buckets[dim][index];
00164             if( bucket.mySize > 0 )
00165                 bucket.boundingBox.update( box );
00166             else
00167                 bucket.boundingBox = box;
00168             bucket.mySize++;
00169         }
00170     }
00171 
00172 #ifndef NDEBUG
00173     BoundBox elt_union = begin->myBox;
00174     for( HandleDataVec::const_iterator i = begin; i != end; ++i )
00175     {
00176         const BoundBox& box = i->myBox;
00177         elt_union.update( box );
00178         for( unsigned int dim = 0; dim < 3; ++dim )
00179         {
00180             const unsigned int index = Bucket::bucket_index( splitsPerDir, box, interval, dim );
00181             Bucket& bucket           = buckets[dim][index];
00182             if( !bucket.boundingBox.intersects_box( box ) ) std::cerr << "Buckets not covering elements!" << std::endl;
00183         }
00184     }
00185     if( !elt_union.intersects_box( interval ) )
00186     {
00187         std::cout << "element union: " << std::endl << elt_union;
00188         std::cout << "intervals: " << std::endl << interval;
00189         std::cout << "union of elts does not contain original box!" << std::endl;
00190     }
00191     if( !interval.intersects_box( elt_union ) )
00192     {
00193         std::cout << "original box does not contain union of elts" << std::endl;
00194         std::cout << interval << std::endl << elt_union << std::endl;
00195     }
00196     for( unsigned int d = 0; d < 3; ++d )
00197     {
00198         std::vector< unsigned int > nonempty;
00199         const std::vector< Bucket >& buckets_ = buckets[d];
00200         unsigned int j                        = 0;
00201         for( std::vector< Bucket >::const_iterator i = buckets_.begin(); i != buckets_.end(); ++i, ++j )
00202         {
00203             if( i->mySize > 0 )
00204             {
00205                 nonempty.push_back( j );
00206             }
00207         }
00208         BoundBox test_box = buckets_[nonempty.front()].boundingBox;
00209         for( unsigned int i = 0; i < nonempty.size(); ++i )
00210             test_box.update( buckets_[nonempty[i]].boundingBox );
00211 
00212         if( !test_box.intersects_box( interval ) )
00213             std::cout << "union of buckets in dimension: " << d << "does not contain original box!" << std::endl;
00214         if( !interval.intersects_box( test_box ) )
00215         {
00216             std::cout << "original box does "
00217                       << "not contain union of buckets"
00218                       << "in dimension: " << d << std::endl;
00219             std::cout << interval << std::endl << test_box << std::endl;
00220         }
00221     }
00222 #endif
00223 }
00224 
00225 void BVHTree::initialize_splits( std::vector< std::vector< SplitData > >& splits,
00226                                  const std::vector< std::vector< Bucket > >& buckets,
00227                                  const SplitData& data ) const
00228 {
00229     for( unsigned int d = 0; d < 3; ++d )
00230     {
00231         std::vector< SplitData >::iterator splits_begin  = splits[d].begin();
00232         std::vector< SplitData >::iterator splits_end    = splits[d].end();
00233         std::vector< Bucket >::const_iterator left_begin = buckets[d].begin();
00234         std::vector< Bucket >::const_iterator _end       = buckets[d].end();
00235         std::vector< Bucket >::const_iterator left_end   = buckets[d].begin() + 1;
00236         for( std::vector< SplitData >::iterator s = splits_begin; s != splits_end; ++s, ++left_end )
00237         {
00238             s->nl = set_interval( s->leftBox, left_begin, left_end );
00239             if( s->nl == 0 )
00240             {
00241                 s->leftBox         = data.boundingBox;
00242                 s->leftBox.bMax[d] = s->leftBox.bMin[d];
00243             }
00244             s->nr = set_interval( s->rightBox, left_end, _end );
00245             if( s->nr == 0 )
00246             {
00247                 s->rightBox         = data.boundingBox;
00248                 s->rightBox.bMin[d] = s->rightBox.bMax[d];
00249             }
00250             s->Lmax  = s->leftBox.bMax[d];
00251             s->Rmin  = s->rightBox.bMin[d];
00252             s->split = std::distance( splits_begin, s );
00253             s->dim   = d;
00254         }
00255 #ifndef NDEBUG
00256         for( std::vector< SplitData >::iterator s = splits_begin; s != splits_end; ++s )
00257         {
00258             BoundBox test_box = s->leftBox;
00259             test_box.update( s->rightBox );
00260             if( !data.boundingBox.intersects_box( test_box ) )
00261             {
00262                 std::cout << "nr: " << s->nr << std::endl;
00263                 std::cout << "Test box: " << std::endl << test_box;
00264                 std::cout << "Left box: " << std::endl << s->leftBox;
00265                 std::cout << "Right box: " << std::endl << s->rightBox;
00266                 std::cout << "Interval: " << std::endl << data.boundingBox;
00267                 std::cout << "Split boxes larger than bb" << std::endl;
00268             }
00269             if( !test_box.intersects_box( data.boundingBox ) )
00270             {
00271                 std::cout << "bb larger than union of split boxes" << std::endl;
00272             }
00273         }
00274 #endif
00275     }
00276 }
00277 
00278 void BVHTree::median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const
00279 {
00280     int dim = data.dim;
00281     for( HandleDataVec::iterator i = begin; i != end; ++i )
00282     {
00283         i->myDim = 0.5 * ( i->myBox.bMin[dim], i->myBox.bMax[dim] );
00284     }
00285     std::sort( begin, end, BVHTree::HandleData_comparator() );
00286     const unsigned int total       = std::distance( begin, end );
00287     HandleDataVec::iterator middle = begin + ( total / 2 );
00288     double middle_center           = middle->myDim;
00289     middle_center += ( ++middle )->myDim;
00290     middle_center *= 0.5;
00291     data.split = middle_center;
00292     data.nl    = std::distance( begin, middle ) + 1;
00293     data.nr    = total - data.nl;
00294     ++middle;
00295     data.leftBox  = begin->myBox;
00296     data.rightBox = middle->myBox;
00297     for( HandleDataVec::iterator i = begin; i != middle; ++i )
00298     {
00299         i->myDim = 0;
00300         data.leftBox.update( i->myBox );
00301     }
00302     for( HandleDataVec::iterator i = middle; i != end; ++i )
00303     {
00304         i->myDim = 1;
00305         data.rightBox.update( i->myBox );
00306     }
00307     data.Rmin = data.rightBox.bMin[data.dim];
00308     data.Lmax = data.leftBox.bMax[data.dim];
00309 #ifndef NDEBUG
00310     BoundBox test_box( data.rightBox );
00311     if( !data.boundingBox.intersects_box( test_box ) )
00312     {
00313         std::cerr << "MEDIAN: BB Does not contain splits" << std::endl;
00314         std::cerr << "test_box:         " << test_box << std::endl;
00315         std::cerr << "data.boundingBox: " << data.boundingBox << std::endl;
00316     }
00317 #endif
00318 }
00319 
00320 void BVHTree::find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const
00321 {
00322     std::vector< std::vector< Bucket > > buckets( 3, std::vector< Bucket >( splitsPerDir + 1 ) );
00323     std::vector< std::vector< SplitData > > splits( 3, std::vector< SplitData >( splitsPerDir, data ) );
00324 
00325     const BoundBox interval = data.boundingBox;
00326     establish_buckets( begin, end, interval, buckets );
00327     initialize_splits( splits, buckets, data );
00328     choose_best_split( splits, data );
00329     const bool use_median = ( 0 == data.nl ) || ( data.nr == 0 );
00330     if( !use_median )
00331         order_elements( begin, end, data );
00332     else
00333         median_order( begin, end, data );
00334 
00335 #ifndef NDEBUG
00336     bool seen_one = false, issue = false;
00337     bool first_left = true, first_right = true;
00338     unsigned int count_left = 0, count_right = 0;
00339     BoundBox left_box, right_box;
00340     for( HandleDataVec::iterator i = begin; i != end; ++i )
00341     {
00342         int order = i->myDim;
00343         if( order != 0 && order != 1 )
00344         {
00345             std::cerr << "Invalid order element !";
00346             std::cerr << order << std::endl;
00347             std::exit( -1 );
00348         }
00349         if( order == 1 )
00350         {
00351             seen_one = 1;
00352             count_right++;
00353             if( first_right )
00354             {
00355                 right_box   = i->myBox;
00356                 first_right = false;
00357             }
00358             else
00359             {
00360                 right_box.update( i->myBox );
00361             }
00362             if( !right_box.intersects_box( i->myBox ) )
00363             {
00364                 if( !issue )
00365                 {
00366                     std::cerr << "Bounding right box issue!" << std::endl;
00367                 }
00368                 issue = true;
00369             }
00370         }
00371         if( order == 0 )
00372         {
00373             count_left++;
00374             if( first_left )
00375             {
00376                 left_box   = i->myBox;
00377                 first_left = false;
00378             }
00379             else
00380             {
00381                 left_box.update( i->myBox );
00382             }
00383             if( !data.leftBox.intersects_box( i->myBox ) )
00384             {
00385                 if( !issue )
00386                 {
00387                     std::cerr << "Bounding left box issue!" << std::endl;
00388                 }
00389                 issue = true;
00390             }
00391             if( seen_one )
00392             {
00393                 std::cerr << "Invalid ordering!" << std::endl;
00394                 std::cout << ( i - 1 )->myDim << order << std::endl;
00395                 exit( -1 );
00396             }
00397         }
00398     }
00399     if( !left_box.intersects_box( data.leftBox ) ) std::cout << "left elts do not contain left box" << std::endl;
00400     if( !data.leftBox.intersects_box( left_box ) ) std::cout << "left box does not contain left elts" << std::endl;
00401     if( !right_box.intersects_box( data.rightBox ) ) std::cout << "right elts do not contain right box" << std::endl;
00402     if( !data.rightBox.intersects_box( right_box ) ) std::cout << "right box do not contain right elts" << std::endl;
00403 
00404     if( count_left != data.nl || count_right != data.nr )
00405     {
00406         std::cerr << "counts are off!" << std::endl;
00407         std::cerr << "total: " << std::distance( begin, end ) << std::endl;
00408         std::cerr << "Dim: " << data.dim << std::endl;
00409         std::cerr << data.Lmax << " , " << data.Rmin << std::endl;
00410         std::cerr << "Right box: " << std::endl << data.rightBox << "Left box: " << std::endl << data.leftBox;
00411         std::cerr << "supposed to be: " << data.nl << " " << data.nr << std::endl;
00412         std::cerr << "accountant says: " << count_left << " " << count_right << std::endl;
00413         std::exit( -1 );
00414     }
00415 #endif
00416 }
00417 
00418 int BVHTree::local_build_tree( std::vector< Node >& tree_nodes,
00419                                HandleDataVec::iterator begin,
00420                                HandleDataVec::iterator end,
00421                                const int index,
00422                                const BoundBox& box,
00423                                const int depth )
00424 {
00425 #ifndef NDEBUG
00426     for( HandleDataVec::const_iterator i = begin; i != end; ++i )
00427     {
00428         if( !box.intersects_box( i->myBox, 0 ) )
00429         {
00430             std::cerr << "depth: " << depth << std::endl;
00431             std::cerr << "BB:" << box << "EB:" << i->myBox << std::endl;
00432             std::exit( -1 );
00433         }
00434     }
00435 #endif
00436 
00437     const unsigned int total_num_elements = std::distance( begin, end );
00438     tree_nodes[index].box                 = box;
00439 
00440     // logic for splitting conditions
00441     if( (int)total_num_elements > maxPerLeaf && depth < maxDepth )
00442     {
00443         SplitData data;
00444         data.boundingBox = box;
00445         find_split( begin, end, data );
00446         // assign data to node
00447         tree_nodes[index].Lmax  = data.Lmax;
00448         tree_nodes[index].Rmin  = data.Rmin;
00449         tree_nodes[index].dim   = data.dim;
00450         tree_nodes[index].child = tree_nodes.size();
00451         // insert left, right children;
00452         tree_nodes.push_back( Node() );
00453         tree_nodes.push_back( Node() );
00454         const int left_depth =
00455             local_build_tree( tree_nodes, begin, begin + data.nl, tree_nodes[index].child, data.leftBox, depth + 1 );
00456         const int right_depth =
00457             local_build_tree( tree_nodes, begin + data.nl, end, tree_nodes[index].child + 1, data.rightBox, depth + 1 );
00458         return std::max( left_depth, right_depth );
00459     }
00460 
00461     tree_nodes[index].dim = 3;
00462     std::copy( begin, end, std::back_inserter( tree_nodes[index].entities ) );
00463     return depth;
00464 }
00465 
00466 ErrorCode BVHTree::find_point( const std::vector< double >& point,
00467                                const unsigned int& index,
00468                                const double iter_tol,
00469                                const double inside_tol,
00470                                std::pair< EntityHandle, CartVect >& result )
00471 {
00472     if( index == 0 ) treeStats.numTraversals++;
00473     const TreeNode& node = myTree[index];
00474     treeStats.nodesVisited++;
00475     CartVect params;
00476     int is_inside;
00477     ErrorCode rval = MB_SUCCESS;
00478     if( node.dim == 3 )
00479     {
00480         treeStats.leavesVisited++;
00481         Range entities;
00482         rval = mbImpl->get_entities_by_handle( startSetHandle + index, entities );
00483         if( MB_SUCCESS != rval ) return rval;
00484 
00485         for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
00486         {
00487             treeStats.traversalLeafObjectTests++;
00488             myEval->set_ent_handle( *i );
00489             myEval->reverse_eval( &point[0], iter_tol, inside_tol, params.array(), &is_inside );
00490             if( is_inside )
00491             {
00492                 result.first  = *i;
00493                 result.second = params;
00494                 return MB_SUCCESS;
00495             }
00496         }
00497         result.first = 0;
00498         return MB_SUCCESS;
00499     }
00500     // the extra tol here considers the case where
00501     // 0 < Rmin - Lmax < 2tol
00502     std::vector< EntityHandle > children;
00503     rval = mbImpl->get_child_meshsets( startSetHandle + index, children );
00504     if( MB_SUCCESS != rval || children.size() != 2 ) return rval;
00505 
00506     if( ( node.Lmax + iter_tol ) < ( node.Rmin - iter_tol ) )
00507     {
00508         if( point[node.dim] <= ( node.Lmax + iter_tol ) )
00509             return find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
00510         else if( point[node.dim] >= ( node.Rmin - iter_tol ) )
00511             return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
00512         result = std::make_pair( 0, CartVect( &point[0] ) );  // point lies in empty space.
00513         return MB_SUCCESS;
00514     }
00515 
00516     // Boxes overlap
00517     // left of Rmin, you must be on the left
00518     // we can't be sure about the boundaries since the boxes overlap
00519     // this was a typo in the paper which caused pain.
00520     if( point[node.dim] < ( node.Rmin - iter_tol ) )
00521         return find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
00522     // if you are on the right Lmax, you must be on the right
00523     else if( point[node.dim] > ( node.Lmax + iter_tol ) )
00524         return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
00525 
00526     /* pg5 of paper
00527      * However, instead of always traversing either subtree
00528      * first (e.g. left always before right), we first traverse
00529      * the subtree whose bounding plane has the larger distance to the
00530      * sought point. This results in less overall traversal, and the correct
00531      * cell is identified more quickly.
00532      */
00533     // So far all testing confirms that this 'heuristic' is
00534     // significantly slower.
00535     // I conjecture this is because it gets improperly
00536     // branch predicted..
00537     // bool dir = (point[ node.dim] - node.Rmin) <=
00538     //              (node.Lmax - point[ node.dim]);
00539     find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
00540     if( result.first == 0 )
00541     {
00542         return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
00543     }
00544     return MB_SUCCESS;
00545 }
00546 
00547 EntityHandle BVHTree::bruteforce_find( const double* point, const double iter_tol, const double inside_tol )
00548 {
00549     treeStats.numTraversals++;
00550     CartVect params;
00551     for( unsigned int i = 0; i < myTree.size(); i++ )
00552     {
00553         if( myTree[i].dim != 3 || !myTree[i].box.contains_point( point, iter_tol ) ) continue;
00554         if( myEval )
00555         {
00556             EntityHandle entity = 0;
00557             treeStats.leavesVisited++;
00558             ErrorCode rval = myEval->find_containing_entity( startSetHandle + i, point, iter_tol, inside_tol, entity,
00559                                                              params.array(), &treeStats.traversalLeafObjectTests );
00560             if( entity )
00561                 return entity;
00562             else if( MB_SUCCESS != rval )
00563                 return 0;
00564         }
00565         else
00566             return startSetHandle + i;
00567     }
00568     return 0;
00569 }
00570 
00571 ErrorCode BVHTree::get_bounding_box( BoundBox& box, EntityHandle* tree_node ) const
00572 {
00573     if( !tree_node || *tree_node == startSetHandle )
00574     {
00575         box = boundBox;
00576         return MB_SUCCESS;
00577     }
00578     else if( ( tree_node && !startSetHandle ) || *tree_node < startSetHandle ||
00579              *tree_node - startSetHandle > myTree.size() )
00580         return MB_FAILURE;
00581 
00582     box = myTree[*tree_node - startSetHandle].box;
00583     return MB_SUCCESS;
00584 }
00585 
00586 ErrorCode BVHTree::point_search( const double* point,
00587                                  EntityHandle& leaf_out,
00588                                  const double iter_tol,
00589                                  const double inside_tol,
00590                                  bool* multiple_leaves,
00591                                  EntityHandle* start_node,
00592                                  CartVect* params )
00593 {
00594     treeStats.numTraversals++;
00595 
00596     EntityHandle this_set = ( start_node ? *start_node : startSetHandle );
00597     // convoluted check because the root is different from startSetHandle
00598     if( this_set != myRoot && ( this_set < startSetHandle || this_set >= startSetHandle + myTree.size() ) )
00599         return MB_FAILURE;
00600     else if( this_set == myRoot )
00601         this_set = startSetHandle;
00602 
00603     std::vector< EntityHandle > candidates,
00604         result_list;  // list of subtrees to traverse, and results
00605     candidates.push_back( this_set - startSetHandle );
00606 
00607     BoundBox box;
00608     while( !candidates.empty() )
00609     {
00610         EntityHandle ind = candidates.back();
00611         treeStats.nodesVisited++;
00612         if( myTree[ind].dim == 3 ) treeStats.leavesVisited++;
00613         this_set = startSetHandle + ind;
00614         candidates.pop_back();
00615 
00616         // test box of this node
00617         ErrorCode rval = get_bounding_box( box, &this_set );
00618         if( MB_SUCCESS != rval ) return rval;
00619         if( !box.contains_point( point, iter_tol ) ) continue;
00620 
00621         // else if not a leaf, test children & put on list
00622         else if( myTree[ind].dim != 3 )
00623         {
00624             candidates.push_back( myTree[ind].child );
00625             candidates.push_back( myTree[ind].child + 1 );
00626             continue;
00627         }
00628         else if( myTree[ind].dim == 3 && myEval && params )
00629         {
00630             rval = myEval->find_containing_entity( startSetHandle + ind, point, iter_tol, inside_tol, leaf_out,
00631                                                    params->array(), &treeStats.traversalLeafObjectTests );
00632             if( leaf_out || MB_SUCCESS != rval ) return rval;
00633         }
00634         else
00635         {
00636             // leaf node within distance; return in list
00637             result_list.push_back( this_set );
00638         }
00639     }
00640 
00641     if( !result_list.empty() ) leaf_out = result_list[0];
00642     if( multiple_leaves && result_list.size() > 1 ) *multiple_leaves = true;
00643     return MB_SUCCESS;
00644 }
00645 
00646 ErrorCode BVHTree::distance_search( const double from_point[3],
00647                                     const double distance,
00648                                     std::vector< EntityHandle >& result_list,
00649                                     const double iter_tol,
00650                                     const double inside_tol,
00651                                     std::vector< double >* result_dists,
00652                                     std::vector< CartVect >* result_params,
00653                                     EntityHandle* tree_root )
00654 {
00655     // non-NULL root should be in tree
00656     // convoluted check because the root is different from startSetHandle
00657     EntityHandle this_set = ( tree_root ? *tree_root : startSetHandle );
00658     if( this_set != myRoot && ( this_set < startSetHandle || this_set >= startSetHandle + myTree.size() ) )
00659         return MB_FAILURE;
00660     else if( this_set == myRoot )
00661         this_set = startSetHandle;
00662 
00663     treeStats.numTraversals++;
00664 
00665     const double dist_sqr = distance * distance;
00666     const CartVect from( from_point );
00667     std::vector< EntityHandle > candidates;  // list of subtrees to traverse
00668                                              // pre-allocate space for default max tree depth
00669     candidates.reserve( maxDepth );
00670 
00671     // misc temporary values
00672     ErrorCode rval;
00673     BoundBox box;
00674 
00675     candidates.push_back( this_set - startSetHandle );
00676 
00677     while( !candidates.empty() )
00678     {
00679 
00680         EntityHandle ind = candidates.back();
00681         this_set         = startSetHandle + ind;
00682         candidates.pop_back();
00683         treeStats.nodesVisited++;
00684         if( myTree[ind].dim == 3 ) treeStats.leavesVisited++;
00685 
00686         // test box of this node
00687         rval = get_bounding_box( box, &this_set );
00688         if( MB_SUCCESS != rval ) return rval;
00689         double d_sqr = box.distance_squared( from_point );
00690 
00691         // if greater than cutoff, continue
00692         if( d_sqr > dist_sqr ) continue;
00693 
00694         // else if not a leaf, test children & put on list
00695         else if( myTree[ind].dim != 3 )
00696         {
00697             candidates.push_back( myTree[ind].child );
00698             candidates.push_back( myTree[ind].child + 1 );
00699             continue;
00700         }
00701 
00702         if( myEval && result_params )
00703         {
00704             EntityHandle ent;
00705             CartVect params;
00706             rval = myEval->find_containing_entity( startSetHandle + ind, from_point, iter_tol, inside_tol, ent,
00707                                                    params.array(), &treeStats.traversalLeafObjectTests );
00708             if( MB_SUCCESS != rval )
00709                 return rval;
00710             else if( ent )
00711             {
00712                 result_list.push_back( ent );
00713                 result_params->push_back( params );
00714                 if( result_dists ) result_dists->push_back( 0.0 );
00715             }
00716         }
00717         else
00718         {
00719             // leaf node within distance; return in list
00720             result_list.push_back( this_set );
00721             if( result_dists ) result_dists->push_back( sqrt( d_sqr ) );
00722         }
00723     }
00724 
00725     return MB_SUCCESS;
00726 }
00727 
00728 ErrorCode BVHTree::print_nodes( std::vector< Node >& nodes )
00729 {
00730     int i;
00731     std::vector< Node >::iterator it;
00732     for( it = nodes.begin(), i = 0; it != nodes.end(); ++it, i++ )
00733     {
00734         std::cout << "Node " << i << ": dim = " << it->dim << ", child = " << it->child << ", Lmax/Rmin = " << it->Lmax
00735                   << "/" << it->Rmin << ", box = " << it->box << std::endl;
00736     }
00737     return MB_SUCCESS;
00738 }
00739 
00740 ErrorCode BVHTree::print()
00741 {
00742     int i;
00743     std::vector< TreeNode >::iterator it;
00744     for( it = myTree.begin(), i = 0; it != myTree.end(); ++it, i++ )
00745     {
00746         std::cout << "Node " << i << ": dim = " << it->dim << ", child = " << it->child << ", Lmax/Rmin = " << it->Lmax
00747                   << "/" << it->Rmin << ", box = " << it->box << std::endl;
00748     }
00749     return MB_SUCCESS;
00750 }
00751 
00752 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines