Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
BVHTree.hpp
Go to the documentation of this file.
1 /**\file BVHTree.hpp
2  * \class moab::BVHTree
3  * \brief Bounding Volume Hierarchy (sorta like a tree), for sorting and searching entities
4  * spatially
5  */
6 
7 #ifndef BVH_TREE_HPP
8 #define BVH_TREE_HPP
9 
10 #include "moab/Interface.hpp"
11 #include "moab/CartVect.hpp"
12 #include "moab/BoundBox.hpp"
13 #include "moab/Tree.hpp"
14 #include "moab/Range.hpp"
15 #include "moab/CN.hpp"
16 
17 #include <vector>
18 #include <cfloat>
19 #include <climits>
20 #include <map>
21 #include <set>
22 #include <iostream>
23 
24 #include "moab/win32_config.h"
25 
26 namespace moab
27 {
28 
29 class ElemEvaluator;
30 
31 class BVHTree : public Tree
32 {
33  public:
34  BVHTree( Interface* impl );
35 
37  {
38  reset_tree();
39  }
40 
41  /** \brief Destroy the tree maintained by this object, optionally checking we have the right
42  * root. \param root If non-NULL, check that this is the root, return failure if not
43  */
44  virtual ErrorCode reset_tree();
45 
46  virtual ErrorCode parse_options( FileOptions& opts );
47 
48  /** \brief Get bounding box for tree below tree_node, or entire tree
49  * If no tree has been built yet, returns +/- DBL_MAX for all dimensions. Note for some tree
50  * types, boxes are not available for non-root nodes, and this function will return failure if
51  * non-root is passed in \param box The box for this tree \param tree_node If non-NULL, node for
52  * which box is requested, tree root if NULL \return Only returns error on fatal condition
53  */
54  virtual ErrorCode get_bounding_box( BoundBox& box, EntityHandle* tree_node = NULL ) const;
55 
56  /** Build the tree
57  * Build a tree with the entities input. If a non-NULL tree_root_set pointer is input,
58  * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct
59  * a new root set and pass its handle back in *tree_root_set. Options vary by tree type;
60  * see Tree.hpp for common options; options specific to AdaptiveKDTree:
61  * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
62  * CANDIDATE_PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
63  * for possible values; default = 1 (SUBDIVISION_SNAP)
64  * \param entities Entities with which to build the tree
65  * \param tree_root Root set for tree (see function description)
66  * \param opts Options for tree (see function description)
67  * \return Error is returned only on build failure
68  */
69  virtual ErrorCode build_tree( const Range& entities,
70  EntityHandle* tree_root_set = NULL,
71  FileOptions* options = NULL );
72 
73  /** \brief Get leaf containing input position.
74  *
75  * Does not take into account global bounding box of tree.
76  * - Therefore there is always one leaf containing the point.
77  * - If caller wants to account for global bounding box, then
78  * caller can test against that box and not call this method
79  * at all if the point is outside the box, as there is no leaf
80  * containing the point in that case.
81  * \param point Point to be located in tree
82  * \param leaf_out Leaf containing point
83  * \param iter_tol Tolerance for convergence of point search
84  * \param inside_tol Tolerance for inside element calculation
85  * \param multiple_leaves Some tree types can have multiple leaves containing a point;
86  * if non-NULL, this parameter is returned true if multiple leaves contain
87  * the input point
88  * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
89  * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
90  */
91  virtual ErrorCode point_search( const double* point,
92  EntityHandle& leaf_out,
93  const double iter_tol = 1.0e-10,
94  const double inside_tol = 1.0e-6,
95  bool* multiple_leaves = NULL,
96  EntityHandle* start_node = NULL,
97  CartVect* params = NULL );
98 
99  /** \brief Find all leaves within a given distance from point
100  * If dists_out input non-NULL, also returns distances from each leaf; if
101  * point i is inside leaf, 0 is given as dists_out[i].
102  * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
103  * in tree nodes and return containing entities in leaves_out. In those cases, if params_out
104  * is also non-NULL, will return parameters in those elements in that vector.
105  * \param from_point Point to be located in tree
106  * \param distance Distance within which to query
107  * \param result_list Leaves within distance or containing point
108  * \param iter_tol Tolerance for convergence of point search
109  * \param inside_tol Tolerance for inside element calculation
110  * \param result_dists If non-NULL, will contain distsances to leaves
111  * \param result_params If non-NULL, will contain parameters of the point in the ents in
112  * leaves_out \param tree_root Start from this tree node (non-NULL) instead of tree root (NULL)
113  */
114  virtual ErrorCode distance_search( const double from_point[3],
115  const double distance,
116  std::vector< EntityHandle >& result_list,
117  const double iter_tol = 1.0e-10,
118  const double inside_tol = 1.0e-6,
119  std::vector< double >* result_dists = NULL,
120  std::vector< CartVect >* result_params = NULL,
121  EntityHandle* tree_root = NULL );
122 
123  //! print various things about this tree
124  virtual ErrorCode print();
125 
126  private:
127  // don't allow copy constructor, too complicated
128  BVHTree( const BVHTree& s );
129 
131  {
132  public:
133  HandleData( EntityHandle h, const BoundBox& bx, const double dp ) : myHandle( h ), myBox( bx ), myDim( dp ) {}
134  HandleData() : myHandle( 0 ), myDim( -1 ) {}
137  double myDim;
138  };
139  typedef std::vector< HandleData > HandleDataVec;
140 
141  class SplitData
142  {
143  public:
145  : dim( UINT_MAX ), nl( UINT_MAX ), nr( UINT_MAX ), split( DBL_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX )
146  {
147  }
148  SplitData( const SplitData& f )
149  : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ),
151  {
152  }
153  unsigned int dim, nl, nr;
154  double split;
155  double Lmax, Rmin;
158  {
159  dim = f.dim;
160  nl = f.nl;
161  nr = f.nr;
162  split = f.split;
163  Lmax = f.Lmax;
164  Rmin = f.Rmin;
166  leftBox = f.leftBox;
167  rightBox = f.rightBox;
168  return *this;
169  }
170  }; // SplitData
171 
173  {
174  // deprecation of binary_function
177  typedef bool result_type;
178 
179  inline double objective( const SplitData& a ) const
180  {
181  return a.Lmax * a.nl - a.Rmin * a.nr;
182  }
183 
184  public:
185  bool operator()( const SplitData& a, const SplitData& b ) const
186  {
187  return objective( a ) < objective( b );
188  }
189  }; // Split_comparator
190 
192  {
193  // deprecation of binary_function
196  typedef bool result_type;
197 
198  public:
199  bool operator()( const HandleData& a, const HandleData& b )
200  {
201  return a.myDim < b.myDim || ( a.myDim == b.myDim && a.myHandle < b.myHandle );
202  }
203  }; // Iterator_comparator
204 
205  class Bucket
206  {
207  public:
208  Bucket() : mySize( 0 ) {}
209  Bucket( const Bucket& f ) : mySize( f.mySize ), boundingBox( f.boundingBox ) {}
210  Bucket( const unsigned int sz ) : mySize( sz ) {}
211  static unsigned int bucket_index( int num_splits,
212  const BoundBox& box,
213  const BoundBox& interval,
214  const unsigned int dim );
215  unsigned int mySize;
217  Bucket& operator=( const Bucket& f )
218  {
220  mySize = f.mySize;
221  return *this;
222  }
223  }; // Bucket
224 
225  class Node
226  {
227  public:
229  unsigned int dim, child;
230  double Lmax, Rmin;
232  Node() : dim( UINT_MAX ), child( UINT_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX ) {}
233  Node( const Node& f ) : entities( f.entities ), dim( f.dim ), child( f.child ), Lmax( f.Lmax ), Rmin( f.Rmin )
234  {
235  }
236  Node& operator=( const Node& f )
237  {
238  dim = f.dim;
239  child = f.child;
240  Lmax = f.Lmax;
241  Rmin = f.Rmin;
242  entities = f.entities;
243  return *this;
244  }
245  }; // Node
246 
247  class TreeNode
248  {
249  public:
250  unsigned int dim, child;
251  double Lmax, Rmin;
253  TreeNode( int dm, int chld, double lmx, double rmn, BoundBox& bx )
254  : dim( dm ), child( chld ), Lmax( lmx ), Rmin( rmn ), box( bx )
255  {
256  }
257  TreeNode( const TreeNode& f )
258  {
259  dim = f.dim;
260  child = f.child;
261  Lmax = f.Lmax;
262  Rmin = f.Rmin;
263  }
265  {
266  dim = f.dim;
267  child = f.child;
268  Lmax = f.Lmax;
269  Rmin = f.Rmin;
270  return *this;
271  }
272  }; // TreeNode
273 
274  void establish_buckets( HandleDataVec::const_iterator begin,
275  HandleDataVec::const_iterator end,
276  const BoundBox& interval,
277  std::vector< std::vector< Bucket > >& buckets ) const;
278 
279  unsigned int set_interval( BoundBox& interval,
280  std::vector< Bucket >::const_iterator begin,
281  std::vector< Bucket >::const_iterator end ) const;
282 
283  void initialize_splits( std::vector< std::vector< SplitData > >& splits,
284  const std::vector< std::vector< Bucket > >& buckets,
285  const SplitData& data ) const;
286 
287  void order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
288 
289  void median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
290 
291  void choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const;
292 
293  void find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
294 
295  ErrorCode find_point( const std::vector< double >& point,
296  const unsigned int& index,
297  const double iter_tol,
298  const double inside_tol,
299  std::pair< EntityHandle, CartVect >& result );
300 
301  EntityHandle bruteforce_find( const double* point,
302  const double iter_tol = 1.0e-10,
303  const double inside_tol = 1.0e-6 );
304 
305  int local_build_tree( std::vector< Node >& tree_nodes,
306  HandleDataVec::iterator begin,
307  HandleDataVec::iterator end,
308  const int index,
309  const BoundBox& box,
310  const int depth = 0 );
311 
312  // builds up vector of HandleData, which caches elements' bounding boxes
313  ErrorCode construct_element_vec( std::vector< HandleData >& handle_data_vec,
314  const Range& elements,
316 
317  // convert the std::vector<Node> to myTree and a bunch of entity sets
318  ErrorCode convert_tree( std::vector< Node >& tree_nodes );
319 
320  // print tree nodes
321  ErrorCode print_nodes( std::vector< Node >& nodes );
322 
324  std::vector< TreeNode > myTree;
327  static MOAB_EXPORT const char* treeName;
328 }; // class Bvh_tree
329 
330 inline unsigned int BVHTree::Bucket::bucket_index( int num_splits,
331  const BoundBox& box,
332  const BoundBox& interval,
333  const unsigned int dim )
334 {
335  // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf
336  // around page 9
337 
338  // Paper arithmetic is over-optimized.. this is safer.
339  const double min = interval.bMin[dim];
340  const double length = ( interval.bMax[dim] - min ) / ( num_splits + 1 );
341  const double center = ( ( box.bMax[dim] + box.bMin[dim] ) / 2.0 ) - min;
342 #ifndef NDEBUG
343 #ifdef BVH_SHOW_INDEX
344  std::cout << "[" << min << " , " << interval.max[dim] << " ]" << std::endl;
345  std::cout << "[" << box.bMin[dim] << " , " << box.bMax[dim] << " ]" << std::endl;
346  std::cout << "Length of bucket" << length << std::endl;
347  std::cout << "Center: " << ( box.bMax[dim] + box.bMin[dim] ) / 2.0 << std::endl;
348  std::cout << "Distance of center from min: " << center << std::endl;
349  std::cout << "ratio: " << center / length << std::endl;
350  std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl;
351 #endif
352 #endif
353  unsigned int cl = std::ceil( center / length );
354  return ( cl > 0 ? cl - 1 : 0 );
355 }
356 
357 inline BVHTree::BVHTree( Interface* impl ) : Tree( impl ), splitsPerDir( 3 ), startSetHandle( 0 )
358 {
360 }
361 
362 inline unsigned int BVHTree::set_interval( BoundBox& interval,
363  std::vector< Bucket >::const_iterator begin,
364  std::vector< Bucket >::const_iterator end ) const
365 {
366  bool first = true;
367  unsigned int count = 0;
368  for( std::vector< Bucket >::const_iterator b = begin; b != end; ++b )
369  {
370  const BoundBox& box = b->boundingBox;
371  count += b->mySize;
372  if( b->mySize != 0 )
373  {
374  if( first )
375  {
376  interval = box;
377  first = false;
378  }
379  else
380  interval.update( box );
381  }
382  }
383  return count;
384 }
385 
386 inline void BVHTree::order_elements( HandleDataVec::iterator& begin,
387  HandleDataVec::iterator& end,
388  SplitData& data ) const
389 {
390  for( HandleDataVec::iterator i = begin; i != end; ++i )
391  {
392  const int index = Bucket::bucket_index( splitsPerDir, i->myBox, data.boundingBox, data.dim );
393  i->myDim = ( index <= data.split ) ? 0 : 1;
394  }
395  std::sort( begin, end, HandleData_comparator() );
396 }
397 
398 inline void BVHTree::choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const
399 {
400  std::vector< SplitData > best_splits;
401  for( std::vector< std::vector< SplitData > >::const_iterator i = splits.begin(); i != splits.end(); ++i )
402  {
403  std::vector< SplitData >::const_iterator j =
404  std::min_element( ( *i ).begin(), ( *i ).end(), Split_comparator() );
405  best_splits.push_back( *j );
406  }
407  data = *std::min_element( best_splits.begin(), best_splits.end(), Split_comparator() );
408 }
409 
410 inline ErrorCode BVHTree::construct_element_vec( std::vector< HandleData >& handle_data_vec,
411  const Range& elements,
413 {
414  std::vector< double > coordinate( 3 * CN::MAX_NODES_PER_ELEMENT );
415  int num_conn;
416  ErrorCode rval = MB_SUCCESS;
417  std::vector< EntityHandle > storage;
418  for( Range::iterator i = elements.begin(); i != elements.end(); ++i )
419  {
420  // TODO: not generic enough. Why dim != 3
421  // Commence un-necessary deep copying.
422  const EntityHandle* connect;
423  rval = mbImpl->get_connectivity( *i, connect, num_conn, false, &storage );
424  if( MB_SUCCESS != rval ) return rval;
425  rval = mbImpl->get_coords( connect, num_conn, &coordinate[0] );
426  if( MB_SUCCESS != rval ) return rval;
427  BoundBox box;
428  for( int j = 0; j < num_conn; j++ )
429  box.update( &coordinate[3 * j] );
430  if( i == elements.begin() )
431  bounding_box = box;
432  else
433  bounding_box.update( box );
434  handle_data_vec.push_back( HandleData( *i, box, 0.0 ) );
435  }
436 
437  return rval;
438 }
439 
441 {
442  return delete_tree_sets();
443 }
444 
445 } // namespace moab
446 
447 #endif // BVH_TREE_HPP