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