Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
bvh_tree.hpp
Go to the documentation of this file.
1 /**
2  * bvh_tree.hpp
3  * Ryan H. Lewis
4  * (C) 2012
5  *
6  * An element tree partitions a mesh composed of elements.
7  * We subdivide the bounding box of a me, by putting boxes
8  * on the left if there center is on the left of a split line
9  * and vice versa.
10  */
11 #include <vector>
12 #include <set>
13 #include <iostream>
14 #include <map>
15 #include <algorithm>
16 #include <bitset>
17 #include <numeric>
18 #include <cmath>
19 #include <tr1/unordered_map>
20 #include <limits>
21 #include "common_tree.hpp"
22 
23 //#define BVH_TREE_DEBUG
24 #ifndef BVH_TREE_HPP
25 #define BVH_TREE_HPP
26 
27 namespace ct = moab::common_tree;
28 
29 namespace moab
30 {
31 
32 // forward declarations
33 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
34 class Bvh_tree;
35 
36 // non-exported functionality
37 namespace
38 {
39  namespace _bvh
40  {
41  template < typename Box, typename Entity_handle >
42  struct _Node
43  {
44  typedef typename std::vector< std::pair< Box, Entity_handle > > Entities;
45  std::size_t dim;
46  std::size_t child;
47  double Lmax, Rmin;
48  Entities entities;
49  _Node& operator=( const _Node& f )
50  {
51  dim = f.dim;
52  child = f.child;
53  Lmax = f.Lmax;
54  Rmin = f.Rmin;
55  entities = f.entities;
56  return *this;
57  }
58  }; // _Node
59 
60  template < typename Split >
61  class Split_comparator
62  {
63  // deprecation of binary_function
64  typedef Split first_argument_type;
65  typedef Split second_argument_type;
66  typedef bool result_type;
67 
68  inline double objective( const Split& a ) const
69  {
70  return a.Lmax * a.nl - a.Rmin * a.nr;
71  }
72 
73  public:
74  bool operator()( const Split& a, const Split& b ) const
75  {
76  return objective( a ) < objective( b );
77  }
78  }; // Split_comparator
79 
80  template < typename Iterator >
81  class Iterator_comparator
82  {
83  // deprecation of binary_function
84  typedef Iterator first_argument_type;
85  typedef Iterator second_argument_type;
86  typedef bool result_type;
87 
88  public:
89  bool operator()( const Iterator& a, const Iterator& b ) const
90  {
91  return a->second.second < b->second.second ||
92  ( !( b->second.second < a->second.second ) && a->first < b->first );
93  }
94  }; // Split_comparator
95 
96  class _Split_data
97  {
98  public:
99  typedef ct::Box< double > Box;
100  _Split_data()
101  : dim( 0 ), nl( 0 ), nr( 0 ), split( 0.0 ), Lmax( 0.0 ), Rmin( 0.0 ), bounding_box(), left_box(),
102  right_box()
103  {
104  }
105  _Split_data( const _Split_data& f )
106  : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ),
108  {
109  }
110  std::size_t dim;
111  std::size_t nl;
112  std::size_t nr;
113  double split;
114  double Lmax, Rmin;
116  Box left_box;
118  _Split_data& operator=( const _Split_data& f )
119  {
120  dim = f.dim;
121  nl = f.nl;
122  nr = f.nr;
123  split = f.split;
124  Lmax = f.Lmax;
125  Rmin = f.Rmin;
126  bounding_box = f.bounding_box;
127  left_box = f.left_box;
128  right_box = f.right_box;
129  return *this;
130  }
131  }; //_Split_data
132 
133  class _Bucket
134  {
135  public:
136  _Bucket() : size( 0 ), bounding_box() {}
137  _Bucket( const _Bucket& f ) : size( f.size ), bounding_box( f.bounding_box ) {}
138  _Bucket( const std::size_t size_ ) : size( size_ ), bounding_box() {}
139  std::size_t size;
141  _Bucket& operator=( const _Bucket& f )
142  {
143  bounding_box = f.bounding_box;
144  size = f.size;
145  return *this;
146  }
147  }; //_Split_data
148  } // namespace _bvh
149 } // namespace
150 
151 template < typename _Entity_handles, typename _Box, typename _Moab, typename _Parametrizer >
152 class Bvh_tree
153 {
154  // public types
155  public:
156  typedef _Entity_handles Entity_handles;
157  typedef _Box Box;
158  typedef _Moab Moab;
159  typedef _Parametrizer Parametrizer;
160  typedef typename Entity_handles::value_type Entity_handle;
161 
162  // private types
163  private:
165  typedef typename std::pair< Box, Entity_handle > Leaf_element;
166  typedef _bvh::_Node< Box, Entity_handle > Node;
167  typedef typename std::vector< Node > Nodes;
168  // public methods
169  public:
170  // Constructor
171  Bvh_tree( Entity_handles& _entities, Moab& _moab, Box& _bounding_box, Parametrizer& _entity_contains )
172  : entity_handles_( _entities ), tree_(), moab( _moab ), bounding_box( _bounding_box ),
173  entity_contains( _entity_contains )
174  {
175  typedef typename Entity_handles::iterator Entity_handle_iterator;
176  typedef ct::_Element_data< const _Box, double > Element_data;
177  typedef typename std::tr1::unordered_map< Entity_handle, Element_data > Entity_map;
178  typedef typename Entity_map::iterator Entity_map_iterator;
179  typedef std::vector< Entity_map_iterator > Vector;
180  // a fully balanced tree will have 2*_entities.size()
181  // which is one doubling away..
182  tree_.reserve( entity_handles_.size() );
183  Entity_map entity_map( entity_handles_.size() );
185 #ifdef BVH_TREE_DEBUG
186  for( Entity_map_iterator i = entity_map.begin(); i != entity_map.end(); ++i )
187  {
188  if( !box_contains_box( bounding_box, i->second.first, 0 ) )
189  {
190  std::cerr << "BB:" << bounding_box << "EB:" << i->second.first << std::endl;
191  std::exit( -1 );
192  }
193  }
194 #endif
195  //_bounding_box = bounding_box;
196  Vector entity_ordering;
197  construct_ordering( entity_map, entity_ordering );
198  // We only build nonempty trees
199  if( entity_ordering.size() )
200  {
201  // initially all bits are set
202  tree_.push_back( Node() );
203  const int depth = build_tree( entity_ordering.begin(), entity_ordering.end(), 0, bounding_box );
204 #ifdef BVH_TREE_DEBUG
205  typedef typename Nodes::iterator Node_iterator;
206  typedef typename Node::Entities::iterator Entity_iterator;
207  std::set< Entity_handle > entity_handles;
208  for( Node_iterator i = tree_.begin(); i != tree_.end(); ++i )
209  {
210  for( Entity_iterator j = i->entities.begin(); j != i->entities.end(); ++j )
211  {
212  entity_handles.insert( j->second );
213  }
214  }
215  if( entity_handles.size() != entity_handles_.size() )
216  {
217  std::cout << "Entity Handle Size Mismatch!" << std::endl;
218  }
219  typedef typename Entity_handles::iterator Entity_iterator_;
220  for( Entity_iterator_ i = entity_handles_.begin(); i != entity_handles_.end(); ++i )
221  {
222  if( entity_handles.find( *i ) == entity_handles.end() )
223  {
224  std::cout << "Tree is missing an entity! " << std::endl;
225  }
226  }
227 
228 #endif
229  std::cout << "max tree depth: " << depth << std::endl;
230  }
231  }
232 
233  // Copy constructor
237  {
238  }
239 
240 // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf
241 // around page 9
242 #define NUM_SPLITS 4
243 #define NUM_BUCKETS ( NUM_SPLITS + 1 ) // NUM_SPLITS+1
244 #define SMAX 5
245  // Paper arithmetic is over-optimized.. this is safer.
246  template < typename Box >
247  std::size_t bucket_index( const Box& box, const Box& interval, const std::size_t dim ) const
248  {
249  const double min = interval.min[dim];
250  const double length = ( interval.max[dim] - min ) / NUM_BUCKETS;
251  const double center = ( ( box.max[dim] + box.min[dim] ) / 2.0 ) - min;
252 #ifdef BVH_TREE_DEBUG
253 #ifdef BVH_SHOW_INDEX
254  std::cout << "[ " << min << " , " << interval.max[dim] << " ]" << std::endl;
255  std::cout << "[ " << box.min[dim] << " , " << box.max[dim] << " ]" << std::endl;
256  std::cout << "Length of bucket" << length << std::endl;
257  std::cout << "Center: " << ( box.max[dim] + box.min[dim] ) / 2.0 << std::endl;
258  std::cout << "Distance of center from min: " << center << std::endl;
259  std::cout << "ratio: " << center / length << std::endl;
260  std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl;
261 #endif
262 #endif
263  return std::ceil( center / length ) - 1;
264  }
265 
266  template < typename Iterator, typename Bounding_box, typename Buckets >
267  void establish_buckets( const Iterator begin,
268  const Iterator end,
269  const Bounding_box& interval,
270  Buckets& buckets ) const
271  {
272  // put each element into its bucket
273  for( Iterator i = begin; i != end; ++i )
274  {
275  const Bounding_box& box = ( *i )->second.first;
276  for( std::size_t dim = 0; dim < NUM_DIM; ++dim )
277  {
278  const std::size_t index = bucket_index( box, interval, dim );
279  _bvh::_Bucket& bucket = buckets[dim][index];
280  if( bucket.size > 0 )
281  {
282  ct::update_bounding_box( bucket.bounding_box, box );
283  }
284  else
285  {
286  bucket.bounding_box = box;
287  }
288  bucket.size++;
289  }
290  }
291 #ifdef BVH_TREE_DEBUG
292  Bounding_box elt_union = ( *begin )->second.first;
293  for( Iterator i = begin; i != end; ++i )
294  {
295  const Bounding_box& box = ( *i )->second.first;
296  ct::update_bounding_box( elt_union, box );
297  for( std::size_t dim = 0; dim < NUM_DIM; ++dim )
298  {
299  const std::size_t index = bucket_index( box, interval, dim );
300  _bvh::_Bucket& bucket = buckets[dim][index];
301  if( !box_contains_box( bucket.bounding_box, box ) )
302  {
303  std::cerr << "Buckets not covering elements!" << std::endl;
304  }
305  }
306  }
307  if( !box_contains_box( elt_union, interval ) )
308  {
309  std::cout << "element union: " << std::endl << elt_union;
310  std::cout << "intervals: " << std::endl << interval;
311  std::cout << "union of elts does not contain original box!" << std::endl;
312  }
313  if( !box_contains_box( interval, elt_union ) )
314  {
315  std::cout << "original box does not contain union of elts" << std::endl;
316  std::cout << interval << std::endl;
317  std::cout << elt_union << std::endl;
318  }
319  typedef typename Buckets::value_type Bucket_list;
320  typedef typename Bucket_list::const_iterator Bucket_iterator;
321  for( std::size_t d = 0; d < NUM_DIM; ++d )
322  {
323  std::vector< std::size_t > nonempty;
324  const Bucket_list& buckets_ = buckets[d];
325  std::size_t j = 0;
326  for( Bucket_iterator i = buckets_.begin(); i != buckets_.end(); ++i, ++j )
327  {
328  if( i->size > 0 )
329  {
330  nonempty.push_back( j );
331  }
332  }
333  Bounding_box test_box = buckets_[nonempty.front()].bounding_box;
334  for( std::size_t i = 0; i < nonempty.size(); ++i )
335  {
336  ct::update_bounding_box( test_box, buckets_[nonempty[i]].bounding_box );
337  }
338  if( !box_contains_box( test_box, interval ) )
339  {
340  std::cout << "union of buckets in dimension: " << d << "does not contain original box!" << std::endl;
341  }
342  if( !box_contains_box( interval, test_box ) )
343  {
344  std::cout << "original box does "
345  << "not contain union of buckets"
346  << "in dimension: " << d << std::endl;
347  std::cout << interval << std::endl;
348  std::cout << test_box << std::endl;
349  }
350  }
351 #endif
352  }
353 
354  template < typename Box, typename Iterator >
355  std::size_t set_interval( Box& interval, const Iterator begin, const Iterator end ) const
356  {
357  bool first = true;
358  std::size_t count = 0;
359  for( Iterator b = begin; b != end; ++b )
360  {
361  const Box& box = b->bounding_box;
362  count += b->size;
363  if( b->size != 0 )
364  {
365  if( first )
366  {
367  interval = box;
368  first = false;
369  }
370  else
371  {
372  ct::update_bounding_box( interval, box );
373  }
374  }
375  }
376  return count;
377  }
378 
379  template < typename Splits, typename Buckets, typename Split_data >
380  void initialize_splits( Splits& splits, const Buckets& buckets, const Split_data& data ) const
381  {
382  typedef typename Buckets::value_type Bucket_list;
383  typedef typename Bucket_list::value_type Bucket;
384  typedef typename Bucket_list::const_iterator Bucket_iterator;
385  typedef typename Splits::value_type Split_list;
386  typedef typename Split_list::value_type Split;
387  typedef typename Split_list::iterator Split_iterator;
388  for( std::size_t d = 0; d < NUM_DIM; ++d )
389  {
390  const Split_iterator splits_begin = splits[d].begin();
391  const Split_iterator splits_end = splits[d].end();
392  const Bucket_iterator left_begin = buckets[d].begin();
393  const Bucket_iterator _end = buckets[d].end();
394  Bucket_iterator left_end = buckets[d].begin() + 1;
395  for( Split_iterator s = splits_begin; s != splits_end; ++s, ++left_end )
396  {
397  s->nl = set_interval( s->left_box, left_begin, left_end );
398  if( s->nl == 0 )
399  {
400  s->left_box = data.bounding_box;
401  s->left_box.max[d] = s->left_box.min[d];
402  }
403  s->nr = set_interval( s->right_box, left_end, _end );
404  if( s->nr == 0 )
405  {
406  s->right_box = data.bounding_box;
407  s->right_box.min[d] = s->right_box.max[d];
408  }
409  s->Lmax = s->left_box.max[d];
410  s->Rmin = s->right_box.min[d];
411  s->split = std::distance( splits_begin, s );
412  s->dim = d;
413  }
414 #ifdef BVH_TREE_DEBUG
415  for( Split_iterator s = splits_begin; s != splits_end; ++s, ++left_end )
416  {
417  typename Split::Box test_box = s->left_box;
418  ct::update_bounding_box( test_box, s->right_box );
419  if( !box_contains_box( data.bounding_box, test_box ) )
420  {
421  std::cout << "nr: " << s->nr << std::endl;
422  std::cout << "Test box: " << std::endl << test_box;
423  std::cout << "Left box: " << std::endl << s->left_box;
424  std::cout << "Right box: " << std::endl << s->right_box;
425  std::cout << "Interval: " << std::endl << data.bounding_box;
426  std::cout << "Split boxes larger than bb" << std::endl;
427  }
428  if( !box_contains_box( test_box, data.bounding_box ) )
429  {
430  std::cout << "bb larger than union "
431  << "of split boxes" << std::endl;
432  }
433  }
434 #endif
435  }
436  }
437 
438  template < typename Iterator, typename Split_data >
439  void order_elements( const Iterator& begin, const Iterator& end, const Split_data& data ) const
440  {
441  typedef typename Iterator::value_type Map_iterator;
442  for( Iterator i = begin; i != end; ++i )
443  {
444  const int index = bucket_index( ( *i )->second.first, data.bounding_box, data.dim );
445  ( *i )->second.second = ( index <= data.split ) ? 0 : 1;
446  }
447  std::sort( begin, end, _bvh::Iterator_comparator< Map_iterator >() );
448  }
449 
450  template < typename Iterator, typename Split_data >
451  void median_order( const Iterator& begin, const Iterator& end, Split_data& data ) const
452  {
453  typedef typename Iterator::value_type Map_iterator;
454  for( Iterator i = begin; i != end; ++i )
455  {
456  const double center = compute_box_center( ( *i )->second.first, data.dim );
457  ( *i )->second.second = center;
458  }
459  std::sort( begin, end, _bvh::Iterator_comparator< Map_iterator >() );
460  const std::size_t total = std::distance( begin, end );
461  Iterator middle = begin + ( total / 2 );
462  double middle_center = ( *middle )->second.second;
463  middle_center += ( *( ++middle ) )->second.second;
464  middle_center /= 2.0;
465  data.split = middle_center;
466  data.nl = std::distance( begin, middle ) + 1;
467  data.nr = total - data.nl;
468  middle++;
469  data.left_box = ( *begin )->second.first;
470  data.right_box = ( *middle )->second.first;
471  for( Iterator i = begin; i != middle; ++i )
472  {
473  ( *i )->second.second = 0;
474  update_bounding_box( data.left_box, ( *i )->second.first );
475  }
476  for( Iterator i = middle; i != end; ++i )
477  {
478  ( *i )->second.second = 1;
479  update_bounding_box( data.right_box, ( *i )->second.first );
480  }
481  data.Rmin = data.right_box.min[data.dim];
482  data.Lmax = data.left_box.max[data.dim];
483 #ifdef BVH_TREE_DEBUG
484  typename Split_data::Box test_box = data.left_box;
485  ct::update_bounding_box( test_box, data.right_box );
486  if( !box_contains_box( data.bounding_box, test_box ) )
487  {
488  std::cerr << "MEDIAN: BB Does not contain splits" << std::endl;
489  }
490  if( !box_contains_box( test_box, data.bounding_box ) )
491  {
492  std::cerr << "MEDIAN: splits do not contain BB" << std::endl;
493  }
494 #endif
495  }
496 
497  template < typename Splits, typename Split_data >
498  void choose_best_split( const Splits& splits, Split_data& data ) const
499  {
500  typedef typename Splits::const_iterator List_iterator;
501  typedef typename List_iterator::value_type::const_iterator Split_iterator;
502  typedef typename Split_iterator::value_type Split;
503  std::vector< Split > best_splits;
504  typedef typename _bvh::Split_comparator< Split > Comparator;
505  Comparator compare;
506  for( List_iterator i = splits.begin(); i != splits.end(); ++i )
507  {
508  Split_iterator j = std::min_element( i->begin(), i->end(), compare );
509  best_splits.push_back( *j );
510  }
511  data = *std::min_element( best_splits.begin(), best_splits.end(), compare );
512  }
513 
514  template < typename Iterator, typename Split_data >
515  void find_split( const Iterator& begin, const Iterator& end, Split_data& data ) const
516  {
517  typedef typename Iterator::value_type Map_iterator;
518  typedef typename Map_iterator::value_type::second_type Box_data;
519  typedef typename Box_data::first_type _Bounding_box; // Note, not global typedef moab::common_tree::Box<
520  // double> Bounding_box;
521  typedef typename std::vector< Split_data > Split_list;
522  typedef typename std::vector< Split_list > Splits;
523  typedef typename Splits::iterator Split_iterator;
524  typedef typename std::vector< _bvh::_Bucket > Bucket_list;
525  typedef typename std::vector< Bucket_list > Buckets;
526  Buckets buckets( NUM_DIM, Bucket_list( NUM_BUCKETS ) );
527  Splits splits( NUM_DIM, Split_list( NUM_SPLITS, data ) );
528 
529  const _Bounding_box interval = data.bounding_box;
530  establish_buckets( begin, end, interval, buckets );
531  initialize_splits( splits, buckets, data );
532  choose_best_split( splits, data );
533  const bool use_median = ( 0 == data.nl ) || ( data.nr == 0 );
534  if( !use_median )
535  {
536  order_elements( begin, end, data );
537  }
538  else
539  {
540  median_order( begin, end, data );
541  }
542 #ifdef BVH_TREE_DEBUG
543  bool seen_one = false, issue = false;
544  bool first_left = true, first_right = true;
545  std::size_t count_left = 0, count_right = 0;
546  typename Split_data::Box left_box, right_box;
547  for( Iterator i = begin; i != end; ++i )
548  {
549  double order = ( *i )->second.second;
550  if( order != 0 && order != 1 )
551  {
552  std::cerr << "Invalid order element !";
553  std::cerr << order << std::endl;
554  std::exit( -1 );
555  }
556  if( order == 1 )
557  {
558  seen_one = 1;
559  count_right++;
560  if( first_right )
561  {
562  right_box = ( *i )->second.first;
563  first_right = false;
564  }
565  else
566  {
567  ct::update_bounding_box( right_box, ( *i )->second.first );
568  }
569  if( !box_contains_box( data.right_box, ( *i )->second.first ) )
570  {
571  if( !issue )
572  {
573  std::cerr << "Bounding right box issue!" << std::endl;
574  }
575  issue = true;
576  }
577  }
578  if( order == 0 )
579  {
580  count_left++;
581  if( first_left )
582  {
583  left_box = ( *i )->second.first;
584  first_left = false;
585  }
586  else
587  {
588  ct::update_bounding_box( left_box, ( *i )->second.first );
589  }
590  if( !box_contains_box( data.left_box, ( *i )->second.first ) )
591  {
592  if( !issue )
593  {
594  std::cerr << "Bounding left box issue!" << std::endl;
595  }
596  issue = true;
597  }
598  if( seen_one )
599  {
600  std::cerr << "Invalid ordering!" << std::endl;
601  std::cout << ( *( i - 1 ) )->second.second << order << std::endl;
602  exit( -1 );
603  }
604  }
605  }
606  if( !box_contains_box( left_box, data.left_box ) )
607  {
608  std::cout << "left elts do not contain left box" << std::endl;
609  }
610  if( !box_contains_box( data.left_box, left_box ) )
611  {
612  std::cout << "left box does not contain left elts" << std::endl;
613  }
614  if( !box_contains_box( right_box, data.right_box ) )
615  {
616  std::cout << "right elts do not contain right box" << std::endl;
617  }
618  if( !box_contains_box( data.right_box, right_box ) )
619  {
620  std::cout << "right box do not contain right elts" << std::endl;
621  }
622  if( count_left != data.nl || count_right != data.nr )
623  {
624  std::cerr << "counts are off!" << std::endl;
625  std::cerr << "total: " << std::distance( begin, end ) << std::endl;
626  std::cerr << "Dim: " << data.dim << std::endl;
627  std::cerr << data.Lmax << " , " << data.Rmin << std::endl;
628  std::cerr << "Right box: " << std::endl << data.right_box << "Left box: " << std::endl << data.left_box;
629  std::cerr << "supposed to be: " << data.nl << " " << data.nr << std::endl;
630  std::cerr << "accountant says: " << count_left << " " << count_right << std::endl;
631  std::exit( -1 );
632  }
633 #endif
634  }
635 
636  // private functionality
637  private:
638  template < typename Iterator >
639  int build_tree( const Iterator begin, const Iterator end, const int index, const Box& box, const int depth = 0 )
640  {
641 #ifdef BVH_TREE_DEBUG
642  for( Iterator i = begin; i != end; ++i )
643  {
644  if( !box_contains_box( box, ( *i )->second.first, 0 ) )
645  {
646  std::cerr << "depth: " << depth << std::endl;
647  std::cerr << "BB:" << box << "EB:" << ( *i )->second.first << std::endl;
648  std::exit( -1 );
649  }
650  }
651 #endif
652 
653  const std::size_t total_num_elements = std::distance( begin, end );
654  Node& node = tree_[index];
655  // logic for splitting conditions
656  if( total_num_elements > SMAX )
657  {
658  _bvh::_Split_data data;
659  data.bounding_box = box;
660  find_split( begin, end, data );
661  // assign data to node
662  node.Lmax = data.Lmax;
663  node.Rmin = data.Rmin;
664  node.dim = data.dim;
665  node.child = tree_.size();
666  // insert left, right children;
667  tree_.push_back( Node() );
668  tree_.push_back( Node() );
669  const int left_depth = build_tree( begin, begin + data.nl, node.child, data.left_box, depth + 1 );
670  const int right_depth = build_tree( begin + data.nl, end, node.child + 1, data.right_box, depth + 1 );
671  return std::max( left_depth, right_depth );
672  }
673  node.dim = 3;
674  ct::assign_entities( node.entities, begin, end );
675  return depth;
676  }
677 
678  template < typename Vector, typename Node_index, typename Result >
679  Result& _find_point( const Vector& point, const Node_index& index, const double tol, Result& result ) const
680  {
681  typedef typename Node::Entities::const_iterator Entity_iterator;
682  const Node& node = tree_[index];
683  if( node.dim == 3 )
684  {
685  for( Entity_iterator i = node.entities.begin(); i != node.entities.end(); ++i )
686  {
687  if( ct::box_contains_point( i->first, point, tol ) )
688  {
689  const std::pair< bool, Vector > r = entity_contains( moab, i->second, point );
690  if( r.first )
691  {
692  return result = std::make_pair( i->second, r.second );
693  }
694  }
695  }
696  result = Result( 0, point );
697  return result;
698  }
699  // the extra tol here considers the case where
700  // 0 < Rmin - Lmax < 2tol
701  if( ( node.Lmax + tol ) < ( node.Rmin - tol ) )
702  {
703  if( point[node.dim] <= ( node.Lmax + tol ) )
704  {
705  return _find_point( point, node.child, tol, result );
706  }
707  else if( point[node.dim] >= ( node.Rmin - tol ) )
708  {
709  return _find_point( point, node.child + 1, tol, result );
710  }
711  result = Result( 0, point );
712  return result; // point lies in empty space.
713  }
714  // Boxes overlap
715  // left of Rmin, you must be on the left
716  // we can't be sure about the boundaries since the boxes overlap
717  // this was a typo in the paper which caused pain.
718  if( point[node.dim] < ( node.Rmin - tol ) )
719  {
720  return _find_point( point, node.child, tol, result );
721  // if you are on the right Lmax, you must be on the right
722  }
723  else if( point[node.dim] > ( node.Lmax + tol ) )
724  {
725  return _find_point( point, node.child + 1, tol, result );
726  }
727  /* pg5 of paper
728  * However, instead of always traversing either subtree
729  * first (e.g. left always before right), we first traverse
730  * the subtree whose bounding plane has the larger distance to the
731  * sought point. This results in less overall traversal, and the correct
732  * cell is identified more quickly.
733  */
734  // Sofar all testing confirms that this 'heuristic' is
735  // significantly slower.
736  // I conjecture this is because it gets improperly
737  // branch predicted..
738  // bool dir = (point[ node.dim] - node.Rmin) <=
739  // (node.Lmax - point[ node.dim]);
740  bool dir = false;
741  _find_point( point, node.child + dir, tol, result );
742  if( result.first == 0 )
743  {
744  return _find_point( point, node.child + ( !dir ), tol, result );
745  }
746  return result;
747  }
748 
749  // public functionality
750  public:
751  template < typename Vector, typename Result >
752  Result& find( const Vector& point, const double tol, Result& result ) const
753  {
754  typedef typename Vector::const_iterator Point_iterator;
755  return _find_point( point, 0, tol, result );
756  }
757 
758  // public functionality
759  public:
760  template < typename Vector >
761  Entity_handle bruteforce_find( const Vector& point, const double tol ) const
762  {
763  typedef typename Vector::const_iterator Point_iterator;
764  typedef typename Nodes::value_type Node;
765  typedef typename Nodes::const_iterator Node_iterator;
766  typedef typename Node::Entities::const_iterator Entity_iterator;
767  for( Node_iterator i = tree_.begin(); i != tree_.end(); ++i )
768  {
769  if( i->dim == 3 )
770  {
771  for( Entity_iterator j = i->entities.begin(); j != i->entities.end(); ++j )
772  {
773  if( ct::box_contains_point( j->first, point, tol ) )
774  {
775  const std::pair< bool, Vector > result = entity_contains( moab, j->second, point );
776  if( result.first )
777  {
778  return j->second;
779  }
780  }
781  }
782  }
783  }
784  return 0;
785  }
786 
787  // public accessor methods
788  public:
789  // private data members
790  private:
796 
797 }; // class Bvh_tree
798 
799 } // namespace moab
800 
801 #endif // BVH_TREE_HPP