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