MOAB: Mesh Oriented datABase  (version 5.2.1)
bvh_tree.hpp
Go to the documentation of this file.
00001 /**
00002  * bvh_tree.hpp
00003  * Ryan H. Lewis
00004  * (C) 2012
00005  *
00006  * An element tree partitions a mesh composed of elements.
00007  * We subdivide the bounding box of a me, by putting boxes
00008  * on the left if there center is on the left of a split line
00009  * and vice versa.
00010  */
00011 #include <vector>
00012 #include <set>
00013 #include <iostream>
00014 #include <map>
00015 #include <algorithm>
00016 #include <bitset>
00017 #include <numeric>
00018 #include <cmath>
00019 #include <tr1/unordered_map>
00020 #include <limits>
00021 #include "common_tree.hpp"
00022 
00023 //#define BVH_TREE_DEBUG
00024 #ifndef BVH_TREE_HPP
00025 #define BVH_TREE_HPP
00026 
00027 namespace ct = moab::common_tree;
00028 
00029 namespace moab
00030 {
00031 
00032 // forward declarations
00033 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
00034 class Bvh_tree;
00035 
00036 // non-exported functionality
00037 namespace
00038 {
00039     namespace _bvh
00040     {
00041         template < typename Box, typename Entity_handle >
00042         struct _Node
00043         {
00044             typedef typename std::vector< std::pair< Box, Entity_handle > > Entities;
00045             std::size_t dim;
00046             std::size_t child;
00047             double Lmax, Rmin;
00048             Entities entities;
00049             _Node& operator=( const _Node& f )
00050             {
00051                 dim      = f.dim;
00052                 child    = f.child;
00053                 Lmax     = f.Lmax;
00054                 Rmin     = f.Rmin;
00055                 entities = f.entities;
00056                 return *this;
00057             }
00058         };  // _Node
00059 
00060         template < typename Split >
00061         class Split_comparator : public std::binary_function< Split, Split, bool >
00062         {
00063             inline double objective( const Split& a ) const
00064             {
00065                 return a.Lmax * a.nl - a.Rmin * a.nr;
00066             }
00067 
00068           public:
00069             bool operator()( const Split& a, const Split& b ) const
00070             {
00071                 return objective( a ) < objective( b );
00072             }
00073         };  // Split_comparator
00074 
00075         template < typename Iterator >
00076         class Iterator_comparator : public std::binary_function< Iterator, Iterator, bool >
00077         {
00078           public:
00079             bool operator()( const Iterator& a, const Iterator& b ) const
00080             {
00081                 return a->second.second < b->second.second ||
00082                        ( !( b->second.second < a->second.second ) && a->first < b->first );
00083             }
00084         };  // Split_comparator
00085 
00086         class _Split_data
00087         {
00088           public:
00089             typedef ct::Box< double > Box;
00090             _Split_data()
00091                 : dim( 0 ), nl( 0 ), nr( 0 ), split( 0.0 ), Lmax( 0.0 ), Rmin( 0.0 ), bounding_box(), left_box(),
00092                   right_box()
00093             {
00094             }
00095             _Split_data( const _Split_data& f )
00096                 : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ),
00097                   bounding_box( f.bounding_box ), left_box( f.left_box ), right_box( f.right_box )
00098             {
00099             }
00100             std::size_t dim;
00101             std::size_t nl;
00102             std::size_t nr;
00103             double split;
00104             double Lmax, Rmin;
00105             Box bounding_box;
00106             Box left_box;
00107             Box right_box;
00108             _Split_data& operator=( const _Split_data& f )
00109             {
00110                 dim          = f.dim;
00111                 nl           = f.nl;
00112                 nr           = f.nr;
00113                 split        = f.split;
00114                 Lmax         = f.Lmax;
00115                 Rmin         = f.Rmin;
00116                 bounding_box = f.bounding_box;
00117                 left_box     = f.left_box;
00118                 right_box    = f.right_box;
00119                 return *this;
00120             }
00121         };  //_Split_data
00122 
00123         class _Bucket
00124         {
00125           public:
00126             _Bucket() : size( 0 ), bounding_box() {}
00127             _Bucket( const _Bucket& f ) : size( f.size ), bounding_box( f.bounding_box ) {}
00128             _Bucket( const std::size_t size_ ) : size( size_ ), bounding_box() {}
00129             std::size_t size;
00130             ct::Box< double > bounding_box;
00131             _Bucket& operator=( const _Bucket& f )
00132             {
00133                 bounding_box = f.bounding_box;
00134                 size         = f.size;
00135                 return *this;
00136             }
00137         };  //_Split_data
00138     }       // namespace _bvh
00139 }  // namespace
00140 
00141 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
00142 class Bvh_tree
00143 {
00144     // public types
00145   public:
00146     typedef _Entity_handles Entity_handles;
00147     typedef _Box Box;
00148     typedef _Moab Moab;
00149     typedef _Parametrizer Parametrizer;
00150     typedef typename Entity_handles::value_type Entity_handle;
00151 
00152     // private types
00153   private:
00154     typedef Bvh_tree< _Entity_handles, _Box, _Moab, _Parametrizer > Self;
00155     typedef typename std::pair< Box, Entity_handle > Leaf_element;
00156     typedef _bvh::_Node< Box, Entity_handle > Node;
00157     typedef typename std::vector< Node > Nodes;
00158     // public methods
00159   public:
00160     // Constructor
00161     Bvh_tree( Entity_handles& _entities, Moab& _moab, Box& _bounding_box, Parametrizer& _entity_contains )
00162         : entity_handles_( _entities ), tree_(), moab( _moab ), bounding_box( _bounding_box ),
00163           entity_contains( _entity_contains )
00164     {
00165         typedef typename Entity_handles::iterator Entity_handle_iterator;
00166         typedef ct::_Element_data< const _Box, double > Element_data;
00167         typedef typename std::tr1::unordered_map< Entity_handle, Element_data > Entity_map;
00168         typedef typename Entity_map::iterator Entity_map_iterator;
00169         typedef std::vector< Entity_map_iterator > Vector;
00170         // a fully balanced tree will have 2*_entities.size()
00171         // which is one doubling away..
00172         tree_.reserve( entity_handles_.size() );
00173         Entity_map entity_map( entity_handles_.size() );
00174         ct::construct_element_map( entity_handles_, entity_map, bounding_box, moab );
00175 #ifdef BVH_TREE_DEBUG
00176         for( Entity_map_iterator i = entity_map.begin(); i != entity_map.end(); ++i )
00177         {
00178             if( !box_contains_box( bounding_box, i->second.first, 0 ) )
00179             {
00180                 std::cerr << "BB:" << bounding_box << "EB:" << i->second.first << std::endl;
00181                 std::exit( -1 );
00182             }
00183         }
00184 #endif
00185         //_bounding_box = bounding_box;
00186         Vector entity_ordering;
00187         construct_ordering( entity_map, entity_ordering );
00188         // We only build nonempty trees
00189         if( entity_ordering.size() )
00190         {
00191             // initially all bits are set
00192             tree_.push_back( Node() );
00193             const int depth = build_tree( entity_ordering.begin(), entity_ordering.end(), 0, bounding_box );
00194 #ifdef BVH_TREE_DEBUG
00195             typedef typename Nodes::iterator Node_iterator;
00196             typedef typename Node::Entities::iterator Entity_iterator;
00197             std::set< Entity_handle > entity_handles;
00198             for( Node_iterator i = tree_.begin(); i != tree_.end(); ++i )
00199             {
00200                 for( Entity_iterator j = i->entities.begin(); j != i->entities.end(); ++j )
00201                 {
00202                     entity_handles.insert( j->second );
00203                 }
00204             }
00205             if( entity_handles.size() != entity_handles_.size() )
00206             { std::cout << "Entity Handle Size Mismatch!" << std::endl; }
00207             typedef typename Entity_handles::iterator Entity_iterator_;
00208             for( Entity_iterator_ i = entity_handles_.begin(); i != entity_handles_.end(); ++i )
00209             {
00210                 if( entity_handles.find( *i ) == entity_handles.end() )
00211                 { std::cout << "Tree is missing an entity! " << std::endl; }
00212             }
00213 
00214 #endif
00215             std::cout << "max tree depth: " << depth << std::endl;
00216         }
00217     }
00218 
00219     // Copy constructor
00220     Bvh_tree( Self& s )
00221         : entity_handles_( s.entity_handles_ ), tree_( s.tree_ ), moab( s.moab ), bounding_box( s.bounding_box ),
00222           entity_contains( s.entity_contains )
00223     {
00224     }
00225 
00226 // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf
00227 // around page 9
00228 #define NUM_SPLITS  4
00229 #define NUM_BUCKETS ( NUM_SPLITS + 1 )  // NUM_SPLITS+1
00230 #define SMAX        5
00231     // Paper arithmetic is over-optimized.. this is safer.
00232     template < typename Box >
00233     std::size_t bucket_index( const Box& box, const Box& interval, const std::size_t dim ) const
00234     {
00235         const double min    = interval.min[dim];
00236         const double length = ( interval.max[dim] - min ) / NUM_BUCKETS;
00237         const double center = ( ( box.max[dim] + box.min[dim] ) / 2.0 ) - min;
00238 #ifdef BVH_TREE_DEBUG
00239 #ifdef BVH_SHOW_INDEX
00240         std::cout << "[ " << min << " , " << interval.max[dim] << " ]" << std::endl;
00241         std::cout << "[ " << box.min[dim] << " , " << box.max[dim] << " ]" << std::endl;
00242         std::cout << "Length of bucket" << length << std::endl;
00243         std::cout << "Center: " << ( box.max[dim] + box.min[dim] ) / 2.0 << std::endl;
00244         std::cout << "Distance of center from min:  " << center << std::endl;
00245         std::cout << "ratio: " << center / length << std::endl;
00246         std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl;
00247 #endif
00248 #endif
00249         return std::ceil( center / length ) - 1;
00250     }
00251 
00252     template < typename Iterator, typename Bounding_box, typename Buckets >
00253     void establish_buckets( const Iterator begin, const Iterator end, const Bounding_box& interval,
00254                             Buckets& buckets ) const
00255     {
00256         // put each element into its bucket
00257         for( Iterator i = begin; i != end; ++i )
00258         {
00259             const Bounding_box& box = ( *i )->second.first;
00260             for( std::size_t dim = 0; dim < NUM_DIM; ++dim )
00261             {
00262                 const std::size_t index = bucket_index( box, interval, dim );
00263                 _bvh::_Bucket& bucket   = buckets[dim][index];
00264                 if( bucket.size > 0 ) { ct::update_bounding_box( bucket.bounding_box, box ); }
00265                 else
00266                 {
00267                     bucket.bounding_box = box;
00268                 }
00269                 bucket.size++;
00270             }
00271         }
00272 #ifdef BVH_TREE_DEBUG
00273         Bounding_box elt_union = ( *begin )->second.first;
00274         for( Iterator i = begin; i != end; ++i )
00275         {
00276             const Bounding_box& box = ( *i )->second.first;
00277             ct::update_bounding_box( elt_union, box );
00278             for( std::size_t dim = 0; dim < NUM_DIM; ++dim )
00279             {
00280                 const std::size_t index = bucket_index( box, interval, dim );
00281                 _bvh::_Bucket& bucket   = buckets[dim][index];
00282                 if( !box_contains_box( bucket.bounding_box, box ) )
00283                 { std::cerr << "Buckets not covering elements!" << std::endl; }
00284             }
00285         }
00286         if( !box_contains_box( elt_union, interval ) )
00287         {
00288             std::cout << "element union: " << std::endl << elt_union;
00289             std::cout << "intervals: " << std::endl << interval;
00290             std::cout << "union of elts does not contain original box!" << std::endl;
00291         }
00292         if( !box_contains_box( interval, elt_union ) )
00293         {
00294             std::cout << "original box does not contain union of elts" << std::endl;
00295             std::cout << interval << std::endl;
00296             std::cout << elt_union << std::endl;
00297         }
00298         typedef typename Buckets::value_type Bucket_list;
00299         typedef typename Bucket_list::const_iterator Bucket_iterator;
00300         for( std::size_t d = 0; d < NUM_DIM; ++d )
00301         {
00302             std::vector< std::size_t > nonempty;
00303             const Bucket_list& buckets_ = buckets[d];
00304             std::size_t j               = 0;
00305             for( Bucket_iterator i = buckets_.begin(); i != buckets_.end(); ++i, ++j )
00306             {
00307                 if( i->size > 0 ) { nonempty.push_back( j ); }
00308             }
00309             Bounding_box test_box = buckets_[nonempty.front()].bounding_box;
00310             for( std::size_t i = 0; i < nonempty.size(); ++i )
00311             {
00312                 ct::update_bounding_box( test_box, buckets_[nonempty[i]].bounding_box );
00313             }
00314             if( !box_contains_box( test_box, interval ) )
00315             { std::cout << "union of buckets in dimension: " << d << "does not contain original box!" << std::endl; }
00316             if( !box_contains_box( interval, test_box ) )
00317             {
00318                 std::cout << "original box does "
00319                           << "not contain union of buckets"
00320                           << "in dimension: " << d << std::endl;
00321                 std::cout << interval << std::endl;
00322                 std::cout << test_box << std::endl;
00323             }
00324         }
00325 #endif
00326     }
00327 
00328     template < typename Box, typename Iterator >
00329     std::size_t set_interval( Box& interval, const Iterator begin, const Iterator end ) const
00330     {
00331         bool first        = true;
00332         std::size_t count = 0;
00333         for( Iterator b = begin; b != end; ++b )
00334         {
00335             const Box& box = b->bounding_box;
00336             count += b->size;
00337             if( b->size != 0 )
00338             {
00339                 if( first )
00340                 {
00341                     interval = box;
00342                     first    = false;
00343                 }
00344                 else
00345                 {
00346                     ct::update_bounding_box( interval, box );
00347                 }
00348             }
00349         }
00350         return count;
00351     }
00352 
00353     template < typename Splits, typename Buckets, typename Split_data >
00354     void initialize_splits( Splits& splits, const Buckets& buckets, const Split_data& data ) const
00355     {
00356         typedef typename Buckets::value_type Bucket_list;
00357         typedef typename Bucket_list::value_type Bucket;
00358         typedef typename Bucket_list::const_iterator Bucket_iterator;
00359         typedef typename Splits::value_type Split_list;
00360         typedef typename Split_list::value_type Split;
00361         typedef typename Split_list::iterator Split_iterator;
00362         for( std::size_t d = 0; d < NUM_DIM; ++d )
00363         {
00364             const Split_iterator splits_begin = splits[d].begin();
00365             const Split_iterator splits_end   = splits[d].end();
00366             const Bucket_iterator left_begin  = buckets[d].begin();
00367             const Bucket_iterator _end        = buckets[d].end();
00368             Bucket_iterator left_end          = buckets[d].begin() + 1;
00369             for( Split_iterator s = splits_begin; s != splits_end; ++s, ++left_end )
00370             {
00371                 s->nl = set_interval( s->left_box, left_begin, left_end );
00372                 if( s->nl == 0 )
00373                 {
00374                     s->left_box        = data.bounding_box;
00375                     s->left_box.max[d] = s->left_box.min[d];
00376                 }
00377                 s->nr = set_interval( s->right_box, left_end, _end );
00378                 if( s->nr == 0 )
00379                 {
00380                     s->right_box        = data.bounding_box;
00381                     s->right_box.min[d] = s->right_box.max[d];
00382                 }
00383                 s->Lmax  = s->left_box.max[d];
00384                 s->Rmin  = s->right_box.min[d];
00385                 s->split = std::distance( splits_begin, s );
00386                 s->dim   = d;
00387             }
00388 #ifdef BVH_TREE_DEBUG
00389             for( Split_iterator s = splits_begin; s != splits_end; ++s, ++left_end )
00390             {
00391                 typename Split::Box test_box = s->left_box;
00392                 ct::update_bounding_box( test_box, s->right_box );
00393                 if( !box_contains_box( data.bounding_box, test_box ) )
00394                 {
00395                     std::cout << "nr: " << s->nr << std::endl;
00396                     std::cout << "Test box: " << std::endl << test_box;
00397                     std::cout << "Left box: " << std::endl << s->left_box;
00398                     std::cout << "Right box: " << std::endl << s->right_box;
00399                     std::cout << "Interval: " << std::endl << data.bounding_box;
00400                     std::cout << "Split boxes larger than bb" << std::endl;
00401                 }
00402                 if( !box_contains_box( test_box, data.bounding_box ) )
00403                 {
00404                     std::cout << "bb larger than union "
00405                               << "of split boxes" << std::endl;
00406                 }
00407             }
00408 #endif
00409         }
00410     }
00411 
00412     template < typename Iterator, typename Split_data >
00413     void order_elements( const Iterator& begin, const Iterator& end, const Split_data& data ) const
00414     {
00415         typedef typename Iterator::value_type Map_iterator;
00416         for( Iterator i = begin; i != end; ++i )
00417         {
00418             const int index       = bucket_index( ( *i )->second.first, data.bounding_box, data.dim );
00419             ( *i )->second.second = ( index <= data.split ) ? 0 : 1;
00420         }
00421         std::sort( begin, end, _bvh::Iterator_comparator< Map_iterator >() );
00422     }
00423 
00424     template < typename Iterator, typename Split_data >
00425     void median_order( const Iterator& begin, const Iterator& end, Split_data& data ) const
00426     {
00427         typedef typename Iterator::value_type Map_iterator;
00428         for( Iterator i = begin; i != end; ++i )
00429         {
00430             const double center   = compute_box_center( ( *i )->second.first, data.dim );
00431             ( *i )->second.second = center;
00432         }
00433         std::sort( begin, end, _bvh::Iterator_comparator< Map_iterator >() );
00434         const std::size_t total = std::distance( begin, end );
00435         Iterator middle         = begin + ( total / 2 );
00436         double middle_center    = ( *middle )->second.second;
00437         middle_center += ( *( ++middle ) )->second.second;
00438         middle_center /= 2.0;
00439         data.split = middle_center;
00440         data.nl    = std::distance( begin, middle ) + 1;
00441         data.nr    = total - data.nl;
00442         middle++;
00443         data.left_box  = ( *begin )->second.first;
00444         data.right_box = ( *middle )->second.first;
00445         for( Iterator i = begin; i != middle; ++i )
00446         {
00447             ( *i )->second.second = 0;
00448             update_bounding_box( data.left_box, ( *i )->second.first );
00449         }
00450         for( Iterator i = middle; i != end; ++i )
00451         {
00452             ( *i )->second.second = 1;
00453             update_bounding_box( data.right_box, ( *i )->second.first );
00454         }
00455         data.Rmin = data.right_box.min[data.dim];
00456         data.Lmax = data.left_box.max[data.dim];
00457 #ifdef BVH_TREE_DEBUG
00458         typename Split_data::Box test_box = data.left_box;
00459         ct::update_bounding_box( test_box, data.right_box );
00460         if( !box_contains_box( data.bounding_box, test_box ) )
00461         { std::cerr << "MEDIAN: BB Does not contain splits" << std::endl; }
00462         if( !box_contains_box( test_box, data.bounding_box ) )
00463         { std::cerr << "MEDIAN: splits do not contain BB" << std::endl; }
00464 #endif
00465     }
00466 
00467     template < typename Splits, typename Split_data >
00468     void choose_best_split( const Splits& splits, Split_data& data ) const
00469     {
00470         typedef typename Splits::const_iterator List_iterator;
00471         typedef typename List_iterator::value_type::const_iterator Split_iterator;
00472         typedef typename Split_iterator::value_type Split;
00473         std::vector< Split > best_splits;
00474         typedef typename _bvh::Split_comparator< Split > Comparator;
00475         Comparator compare;
00476         for( List_iterator i = splits.begin(); i != splits.end(); ++i )
00477         {
00478             Split_iterator j = std::min_element( i->begin(), i->end(), compare );
00479             best_splits.push_back( *j );
00480         }
00481         data = *std::min_element( best_splits.begin(), best_splits.end(), compare );
00482     }
00483 
00484     template < typename Iterator, typename Split_data >
00485     void find_split( const Iterator& begin, const Iterator& end, Split_data& data ) const
00486     {
00487         typedef typename Iterator::value_type Map_iterator;
00488         typedef typename Map_iterator::value_type::second_type Box_data;
00489         typedef typename Box_data::first_type _Bounding_box;  // Note, not global typedef moab::common_tree::Box<
00490                                                               // double> Bounding_box;
00491         typedef typename std::vector< Split_data > Split_list;
00492         typedef typename std::vector< Split_list > Splits;
00493         typedef typename Splits::iterator Split_iterator;
00494         typedef typename std::vector< _bvh::_Bucket > Bucket_list;
00495         typedef typename std::vector< Bucket_list > Buckets;
00496         Buckets buckets( NUM_DIM, Bucket_list( NUM_BUCKETS ) );
00497         Splits splits( NUM_DIM, Split_list( NUM_SPLITS, data ) );
00498 
00499         const _Bounding_box interval = data.bounding_box;
00500         establish_buckets( begin, end, interval, buckets );
00501         initialize_splits( splits, buckets, data );
00502         choose_best_split( splits, data );
00503         const bool use_median = ( 0 == data.nl ) || ( data.nr == 0 );
00504         if( !use_median ) { order_elements( begin, end, data ); }
00505         else
00506         {
00507             median_order( begin, end, data );
00508         }
00509 #ifdef BVH_TREE_DEBUG
00510         bool seen_one = false, issue = false;
00511         bool first_left = true, first_right = true;
00512         std::size_t count_left = 0, count_right = 0;
00513         typename Split_data::Box left_box, right_box;
00514         for( Iterator i = begin; i != end; ++i )
00515         {
00516             double order = ( *i )->second.second;
00517             if( order != 0 && order != 1 )
00518             {
00519                 std::cerr << "Invalid order element !";
00520                 std::cerr << order << std::endl;
00521                 std::exit( -1 );
00522             }
00523             if( order == 1 )
00524             {
00525                 seen_one = 1;
00526                 count_right++;
00527                 if( first_right )
00528                 {
00529                     right_box   = ( *i )->second.first;
00530                     first_right = false;
00531                 }
00532                 else
00533                 {
00534                     ct::update_bounding_box( right_box, ( *i )->second.first );
00535                 }
00536                 if( !box_contains_box( data.right_box, ( *i )->second.first ) )
00537                 {
00538                     if( !issue ) { std::cerr << "Bounding right box issue!" << std::endl; }
00539                     issue = true;
00540                 }
00541             }
00542             if( order == 0 )
00543             {
00544                 count_left++;
00545                 if( first_left )
00546                 {
00547                     left_box   = ( *i )->second.first;
00548                     first_left = false;
00549                 }
00550                 else
00551                 {
00552                     ct::update_bounding_box( left_box, ( *i )->second.first );
00553                 }
00554                 if( !box_contains_box( data.left_box, ( *i )->second.first ) )
00555                 {
00556                     if( !issue ) { std::cerr << "Bounding left box issue!" << std::endl; }
00557                     issue = true;
00558                 }
00559                 if( seen_one )
00560                 {
00561                     std::cerr << "Invalid ordering!" << std::endl;
00562                     std::cout << ( *( i - 1 ) )->second.second << order << std::endl;
00563                     exit( -1 );
00564                 }
00565             }
00566         }
00567         if( !box_contains_box( left_box, data.left_box ) )
00568         { std::cout << "left elts do not contain left box" << std::endl; }
00569         if( !box_contains_box( data.left_box, left_box ) )
00570         { std::cout << "left box does not contain left elts" << std::endl; }
00571         if( !box_contains_box( right_box, data.right_box ) )
00572         { std::cout << "right elts do not contain right box" << std::endl; }
00573         if( !box_contains_box( data.right_box, right_box ) )
00574         { std::cout << "right box do not contain right elts" << std::endl; }
00575         if( count_left != data.nl || count_right != data.nr )
00576         {
00577             std::cerr << "counts are off!" << std::endl;
00578             std::cerr << "total: " << std::distance( begin, end ) << std::endl;
00579             std::cerr << "Dim: " << data.dim << std::endl;
00580             std::cerr << data.Lmax << " , " << data.Rmin << std::endl;
00581             std::cerr << "Right box: " << std::endl << data.right_box << "Left box: " << std::endl << data.left_box;
00582             std::cerr << "supposed to be: " << data.nl << " " << data.nr << std::endl;
00583             std::cerr << "accountant says: " << count_left << " " << count_right << std::endl;
00584             std::exit( -1 );
00585         }
00586 #endif
00587     }
00588 
00589     // private functionality
00590   private:
00591     template < typename Iterator >
00592     int build_tree( const Iterator begin, const Iterator end, const int index, const Box& box, const int depth = 0 )
00593     {
00594 #ifdef BVH_TREE_DEBUG
00595         for( Iterator i = begin; i != end; ++i )
00596         {
00597             if( !box_contains_box( box, ( *i )->second.first, 0 ) )
00598             {
00599                 std::cerr << "depth: " << depth << std::endl;
00600                 std::cerr << "BB:" << box << "EB:" << ( *i )->second.first << std::endl;
00601                 std::exit( -1 );
00602             }
00603         }
00604 #endif
00605 
00606         const std::size_t total_num_elements = std::distance( begin, end );
00607         Node& node                           = tree_[index];
00608         // logic for splitting conditions
00609         if( total_num_elements > SMAX )
00610         {
00611             _bvh::_Split_data data;
00612             data.bounding_box = box;
00613             find_split( begin, end, data );
00614             // assign data to node
00615             node.Lmax  = data.Lmax;
00616             node.Rmin  = data.Rmin;
00617             node.dim   = data.dim;
00618             node.child = tree_.size();
00619             // insert left, right children;
00620             tree_.push_back( Node() );
00621             tree_.push_back( Node() );
00622             const int left_depth  = build_tree( begin, begin + data.nl, node.child, data.left_box, depth + 1 );
00623             const int right_depth = build_tree( begin + data.nl, end, node.child + 1, data.right_box, depth + 1 );
00624             return std::max( left_depth, right_depth );
00625         }
00626         node.dim = 3;
00627         ct::assign_entities( node.entities, begin, end );
00628         return depth;
00629     }
00630 
00631     template < typename Vector, typename Node_index, typename Result >
00632     Result& _find_point( const Vector& point, const Node_index& index, const double tol, Result& result ) const
00633     {
00634         typedef typename Node::Entities::const_iterator Entity_iterator;
00635         const Node& node = tree_[index];
00636         if( node.dim == 3 )
00637         {
00638             for( Entity_iterator i = node.entities.begin(); i != node.entities.end(); ++i )
00639             {
00640                 if( ct::box_contains_point( i->first, point, tol ) )
00641                 {
00642                     const std::pair< bool, Vector > r = entity_contains( moab, i->second, point );
00643                     if( r.first ) { return result = std::make_pair( i->second, r.second ); }
00644                 }
00645             }
00646             result = Result( 0, point );
00647             return result;
00648         }
00649         // the extra tol here considers the case where
00650         // 0 < Rmin - Lmax < 2tol
00651         if( ( node.Lmax + tol ) < ( node.Rmin - tol ) )
00652         {
00653             if( point[node.dim] <= ( node.Lmax + tol ) ) { return _find_point( point, node.child, tol, result ); }
00654             else if( point[node.dim] >= ( node.Rmin - tol ) )
00655             {
00656                 return _find_point( point, node.child + 1, tol, result );
00657             }
00658             result = Result( 0, point );
00659             return result;  // point lies in empty space.
00660         }
00661         // Boxes overlap
00662         // left of Rmin, you must be on the left
00663         // we can't be sure about the boundaries since the boxes overlap
00664         // this was a typo in the paper which caused pain.
00665         if( point[node.dim] < ( node.Rmin - tol ) )
00666         {
00667             return _find_point( point, node.child, tol, result );
00668             // if you are on the right Lmax, you must be on the right
00669         }
00670         else if( point[node.dim] > ( node.Lmax + tol ) )
00671         {
00672             return _find_point( point, node.child + 1, tol, result );
00673         }
00674         /* pg5 of paper
00675          * However, instead of always traversing either subtree
00676          * first (e.g. left always before right), we first traverse
00677          * the subtree whose bounding plane has the larger distance to the
00678          * sought point. This results in less overall traversal, and the correct
00679          * cell is identified more quickly.
00680          */
00681         // Sofar all testing confirms that this 'heuristic' is
00682         // significantly slower.
00683         // I conjecture this is because it gets improperly
00684         // branch predicted..
00685         // bool dir = (point[ node.dim] - node.Rmin) <=
00686         //              (node.Lmax - point[ node.dim]);
00687         bool dir = false;
00688         _find_point( point, node.child + dir, tol, result );
00689         if( result.first == 0 ) { return _find_point( point, node.child + ( !dir ), tol, result ); }
00690         return result;
00691     }
00692 
00693     // public functionality
00694   public:
00695     template < typename Vector, typename Result >
00696     Result& find( const Vector& point, const double tol, Result& result ) const
00697     {
00698         typedef typename Vector::const_iterator Point_iterator;
00699         return _find_point( point, 0, tol, result );
00700     }
00701 
00702     // public functionality
00703   public:
00704     template < typename Vector >
00705     Entity_handle bruteforce_find( const Vector& point, const double tol ) const
00706     {
00707         typedef typename Vector::const_iterator Point_iterator;
00708         typedef typename Nodes::value_type Node;
00709         typedef typename Nodes::const_iterator Node_iterator;
00710         typedef typename Node::Entities::const_iterator Entity_iterator;
00711         for( Node_iterator i = tree_.begin(); i != tree_.end(); ++i )
00712         {
00713             if( i->dim == 3 )
00714             {
00715                 for( Entity_iterator j = i->entities.begin(); j != i->entities.end(); ++j )
00716                 {
00717                     if( ct::box_contains_point( j->first, point, tol ) )
00718                     {
00719                         const std::pair< bool, Vector > result = entity_contains( moab, j->second, point );
00720                         if( result.first ) { return j->second; }
00721                     }
00722                 }
00723             }
00724         }
00725         return 0;
00726     }
00727 
00728     // public accessor methods
00729   public:
00730     // private data members
00731   private:
00732     const Entity_handles& entity_handles_;
00733     Nodes tree_;
00734     Moab& moab;
00735     Box bounding_box;
00736     Parametrizer& entity_contains;
00737 
00738 };  // class Bvh_tree
00739 
00740 }  // namespace moab
00741 
00742 #endif  // BVH_TREE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines