MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /**\file BVHTree.hpp 00002 * \class moab::BVHTree 00003 * \brief Bounding Volume Hierarchy (sorta like a tree), for sorting and searching entities 00004 * spatially 00005 */ 00006 00007 #ifndef BVH_TREE_HPP 00008 #define BVH_TREE_HPP 00009 00010 #include "moab/Interface.hpp" 00011 #include "moab/CartVect.hpp" 00012 #include "moab/BoundBox.hpp" 00013 #include "moab/Tree.hpp" 00014 #include "moab/Range.hpp" 00015 #include "moab/CN.hpp" 00016 00017 #include <vector> 00018 #include <cfloat> 00019 #include <climits> 00020 #include <map> 00021 #include <set> 00022 #include <iostream> 00023 00024 #include "moab/win32_config.h" 00025 00026 namespace moab 00027 { 00028 00029 class ElemEvaluator; 00030 00031 class BVHTree : public Tree 00032 { 00033 public: 00034 BVHTree( Interface* impl ); 00035 00036 ~BVHTree() 00037 { 00038 reset_tree(); 00039 } 00040 00041 /** \brief Destroy the tree maintained by this object, optionally checking we have the right 00042 * root. \param root If non-NULL, check that this is the root, return failure if not 00043 */ 00044 virtual ErrorCode reset_tree(); 00045 00046 virtual ErrorCode parse_options( FileOptions& opts ); 00047 00048 /** \brief Get bounding box for tree below tree_node, or entire tree 00049 * If no tree has been built yet, returns +/- DBL_MAX for all dimensions. Note for some tree 00050 * types, boxes are not available for non-root nodes, and this function will return failure if 00051 * non-root is passed in \param box The box for this tree \param tree_node If non-NULL, node for 00052 * which box is requested, tree root if NULL \return Only returns error on fatal condition 00053 */ 00054 virtual ErrorCode get_bounding_box( BoundBox& box, EntityHandle* tree_node = NULL ) const; 00055 00056 /** Build the tree 00057 * Build a tree with the entities input. If a non-NULL tree_root_set pointer is input, 00058 * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct 00059 * a new root set and pass its handle back in *tree_root_set. Options vary by tree type; 00060 * see Tree.hpp for common options; options specific to AdaptiveKDTree: 00061 * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3 00062 * CANDIDATE_PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below) 00063 * for possible values; default = 1 (SUBDIVISION_SNAP) 00064 * \param entities Entities with which to build the tree 00065 * \param tree_root Root set for tree (see function description) 00066 * \param opts Options for tree (see function description) 00067 * \return Error is returned only on build failure 00068 */ 00069 virtual ErrorCode build_tree( const Range& entities, EntityHandle* tree_root_set = NULL, 00070 FileOptions* options = NULL ); 00071 00072 /** \brief Get leaf containing input position. 00073 * 00074 * Does not take into account global bounding box of tree. 00075 * - Therefore there is always one leaf containing the point. 00076 * - If caller wants to account for global bounding box, then 00077 * caller can test against that box and not call this method 00078 * at all if the point is outside the box, as there is no leaf 00079 * containing the point in that case. 00080 * \param point Point to be located in tree 00081 * \param leaf_out Leaf containing point 00082 * \param iter_tol Tolerance for convergence of point search 00083 * \param inside_tol Tolerance for inside element calculation 00084 * \param multiple_leaves Some tree types can have multiple leaves containing a point; 00085 * if non-NULL, this parameter is returned true if multiple leaves contain 00086 * the input point 00087 * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL) 00088 * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0 00089 */ 00090 virtual ErrorCode point_search( const double* point, EntityHandle& leaf_out, const double iter_tol = 1.0e-10, 00091 const double inside_tol = 1.0e-6, bool* multiple_leaves = NULL, 00092 EntityHandle* start_node = NULL, CartVect* params = NULL ); 00093 00094 /** \brief Find all leaves within a given distance from point 00095 * If dists_out input non-NULL, also returns distances from each leaf; if 00096 * point i is inside leaf, 0 is given as dists_out[i]. 00097 * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities 00098 * in tree nodes and return containing entities in leaves_out. In those cases, if params_out 00099 * is also non-NULL, will return parameters in those elements in that vector. 00100 * \param from_point Point to be located in tree 00101 * \param distance Distance within which to query 00102 * \param result_list Leaves within distance or containing point 00103 * \param iter_tol Tolerance for convergence of point search 00104 * \param inside_tol Tolerance for inside element calculation 00105 * \param result_dists If non-NULL, will contain distsances to leaves 00106 * \param result_params If non-NULL, will contain parameters of the point in the ents in 00107 * leaves_out \param tree_root Start from this tree node (non-NULL) instead of tree root (NULL) 00108 */ 00109 virtual ErrorCode distance_search( const double from_point[3], const double distance, 00110 std::vector< EntityHandle >& result_list, const double iter_tol = 1.0e-10, 00111 const double inside_tol = 1.0e-6, std::vector< double >* result_dists = NULL, 00112 std::vector< CartVect >* result_params = NULL, EntityHandle* tree_root = NULL ); 00113 00114 //! print various things about this tree 00115 virtual ErrorCode print(); 00116 00117 private: 00118 // don't allow copy constructor, too complicated 00119 BVHTree( const BVHTree& s ); 00120 00121 class HandleData 00122 { 00123 public: 00124 HandleData( EntityHandle h, const BoundBox& bx, const double dp ) : myHandle( h ), myBox( bx ), myDim( dp ) {} 00125 HandleData() : myHandle( 0 ), myDim( -1 ) {} 00126 EntityHandle myHandle; 00127 BoundBox myBox; 00128 double myDim; 00129 }; 00130 typedef std::vector< HandleData > HandleDataVec; 00131 00132 class SplitData 00133 { 00134 public: 00135 SplitData() 00136 : dim( UINT_MAX ), nl( UINT_MAX ), nr( UINT_MAX ), split( DBL_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX ) 00137 { 00138 } 00139 SplitData( const SplitData& f ) 00140 : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ), 00141 boundingBox( f.boundingBox ), leftBox( f.leftBox ), rightBox( f.rightBox ) 00142 { 00143 } 00144 unsigned int dim, nl, nr; 00145 double split; 00146 double Lmax, Rmin; 00147 BoundBox boundingBox, leftBox, rightBox; 00148 SplitData& operator=( const SplitData& f ) 00149 { 00150 dim = f.dim; 00151 nl = f.nl; 00152 nr = f.nr; 00153 split = f.split; 00154 Lmax = f.Lmax; 00155 Rmin = f.Rmin; 00156 boundingBox = f.boundingBox; 00157 leftBox = f.leftBox; 00158 rightBox = f.rightBox; 00159 return *this; 00160 } 00161 }; // SplitData 00162 00163 class Split_comparator : public std::binary_function< SplitData, SplitData, bool > 00164 { 00165 inline double objective( const SplitData& a ) const 00166 { 00167 return a.Lmax * a.nl - a.Rmin * a.nr; 00168 } 00169 00170 public: 00171 bool operator()( const SplitData& a, const SplitData& b ) const 00172 { 00173 return objective( a ) < objective( b ); 00174 } 00175 }; // Split_comparator 00176 00177 class HandleData_comparator : public std::binary_function< HandleData, HandleData, bool > 00178 { 00179 public: 00180 bool operator()( const HandleData& a, const HandleData& b ) 00181 { 00182 return a.myDim < b.myDim || ( a.myDim == b.myDim && a.myHandle < b.myHandle ); 00183 } 00184 }; // Iterator_comparator 00185 00186 class Bucket 00187 { 00188 public: 00189 Bucket() : mySize( 0 ) {} 00190 Bucket( const Bucket& f ) : mySize( f.mySize ), boundingBox( f.boundingBox ) {} 00191 Bucket( const unsigned int sz ) : mySize( sz ) {} 00192 static unsigned int bucket_index( int num_splits, const BoundBox& box, const BoundBox& interval, 00193 const unsigned int dim ); 00194 unsigned int mySize; 00195 BoundBox boundingBox; 00196 Bucket& operator=( const Bucket& f ) 00197 { 00198 boundingBox = f.boundingBox; 00199 mySize = f.mySize; 00200 return *this; 00201 } 00202 }; // Bucket 00203 00204 class Node 00205 { 00206 public: 00207 HandleDataVec entities; 00208 unsigned int dim, child; 00209 double Lmax, Rmin; 00210 BoundBox box; 00211 Node() : dim( UINT_MAX ), child( UINT_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX ) {} 00212 Node& operator=( const Node& f ) 00213 { 00214 dim = f.dim; 00215 child = f.child; 00216 Lmax = f.Lmax; 00217 Rmin = f.Rmin; 00218 entities = f.entities; 00219 return *this; 00220 } 00221 }; // Node 00222 00223 class TreeNode 00224 { 00225 public: 00226 unsigned int dim, child; 00227 double Lmax, Rmin; 00228 BoundBox box; 00229 TreeNode( int dm, int chld, double lmx, double rmn, BoundBox& bx ) 00230 : dim( dm ), child( chld ), Lmax( lmx ), Rmin( rmn ), box( bx ) 00231 { 00232 } 00233 TreeNode& operator=( const TreeNode& f ) 00234 { 00235 dim = f.dim; 00236 child = f.child; 00237 Lmax = f.Lmax; 00238 Rmin = f.Rmin; 00239 return *this; 00240 } 00241 }; // TreeNode 00242 00243 void establish_buckets( HandleDataVec::const_iterator begin, HandleDataVec::const_iterator end, 00244 const BoundBox& interval, std::vector< std::vector< Bucket > >& buckets ) const; 00245 00246 unsigned int set_interval( BoundBox& interval, std::vector< Bucket >::const_iterator begin, 00247 std::vector< Bucket >::const_iterator end ) const; 00248 00249 void initialize_splits( std::vector< std::vector< SplitData > >& splits, 00250 const std::vector< std::vector< Bucket > >& buckets, const SplitData& data ) const; 00251 00252 void order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const; 00253 00254 void median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const; 00255 00256 void choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const; 00257 00258 void find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const; 00259 00260 ErrorCode find_point( const std::vector< double >& point, const unsigned int& index, const double iter_tol, 00261 const double inside_tol, std::pair< EntityHandle, CartVect >& result ); 00262 00263 EntityHandle bruteforce_find( const double* point, const double iter_tol = 1.0e-10, 00264 const double inside_tol = 1.0e-6 ); 00265 00266 int local_build_tree( std::vector< Node >& tree_nodes, HandleDataVec::iterator begin, HandleDataVec::iterator end, 00267 const int index, const BoundBox& box, const int depth = 0 ); 00268 00269 // builds up vector of HandleData, which caches elements' bounding boxes 00270 ErrorCode construct_element_vec( std::vector< HandleData >& handle_data_vec, const Range& elements, 00271 BoundBox& bounding_box ); 00272 00273 // convert the std::vector<Node> to myTree and a bunch of entity sets 00274 ErrorCode convert_tree( std::vector< Node >& tree_nodes ); 00275 00276 // print tree nodes 00277 ErrorCode print_nodes( std::vector< Node >& nodes ); 00278 00279 Range entityHandles; 00280 std::vector< TreeNode > myTree; 00281 int splitsPerDir; 00282 EntityHandle startSetHandle; 00283 static MOAB_EXPORT const char* treeName; 00284 }; // class Bvh_tree 00285 00286 inline unsigned int BVHTree::Bucket::bucket_index( int num_splits, const BoundBox& box, const BoundBox& interval, 00287 const unsigned int dim ) 00288 { 00289 // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf 00290 // around page 9 00291 00292 // Paper arithmetic is over-optimized.. this is safer. 00293 const double min = interval.bMin[dim]; 00294 const double length = ( interval.bMax[dim] - min ) / ( num_splits + 1 ); 00295 const double center = ( ( box.bMax[dim] + box.bMin[dim] ) / 2.0 ) - min; 00296 #ifndef NDEBUG 00297 #ifdef BVH_SHOW_INDEX 00298 std::cout << "[" << min << " , " << interval.max[dim] << " ]" << std::endl; 00299 std::cout << "[" << box.bMin[dim] << " , " << box.bMax[dim] << " ]" << std::endl; 00300 std::cout << "Length of bucket" << length << std::endl; 00301 std::cout << "Center: " << ( box.bMax[dim] + box.bMin[dim] ) / 2.0 << std::endl; 00302 std::cout << "Distance of center from min: " << center << std::endl; 00303 std::cout << "ratio: " << center / length << std::endl; 00304 std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl; 00305 #endif 00306 #endif 00307 unsigned int cl = std::ceil( center / length ); 00308 return ( cl > 0 ? cl - 1 : 0 ); 00309 } 00310 00311 inline BVHTree::BVHTree( Interface* impl ) : Tree( impl ), splitsPerDir( 3 ), startSetHandle( 0 ) 00312 { 00313 boxTagName = treeName; 00314 } 00315 00316 inline unsigned int BVHTree::set_interval( BoundBox& interval, std::vector< Bucket >::const_iterator begin, 00317 std::vector< Bucket >::const_iterator end ) const 00318 { 00319 bool first = true; 00320 unsigned int count = 0; 00321 for( std::vector< Bucket >::const_iterator b = begin; b != end; ++b ) 00322 { 00323 const BoundBox& box = b->boundingBox; 00324 count += b->mySize; 00325 if( b->mySize != 0 ) 00326 { 00327 if( first ) 00328 { 00329 interval = box; 00330 first = false; 00331 } 00332 else 00333 interval.update( box ); 00334 } 00335 } 00336 return count; 00337 } 00338 00339 inline void BVHTree::order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, 00340 SplitData& data ) const 00341 { 00342 for( HandleDataVec::iterator i = begin; i != end; ++i ) 00343 { 00344 const int index = Bucket::bucket_index( splitsPerDir, i->myBox, data.boundingBox, data.dim ); 00345 i->myDim = ( index <= data.split ) ? 0 : 1; 00346 } 00347 std::sort( begin, end, HandleData_comparator() ); 00348 } 00349 00350 inline void BVHTree::choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const 00351 { 00352 std::vector< SplitData > best_splits; 00353 for( std::vector< std::vector< SplitData > >::const_iterator i = splits.begin(); i != splits.end(); ++i ) 00354 { 00355 std::vector< SplitData >::const_iterator j = 00356 std::min_element( ( *i ).begin(), ( *i ).end(), Split_comparator() ); 00357 best_splits.push_back( *j ); 00358 } 00359 data = *std::min_element( best_splits.begin(), best_splits.end(), Split_comparator() ); 00360 } 00361 00362 inline ErrorCode BVHTree::construct_element_vec( std::vector< HandleData >& handle_data_vec, const Range& elements, 00363 BoundBox& bounding_box ) 00364 { 00365 std::vector< double > coordinate( 3 * CN::MAX_NODES_PER_ELEMENT ); 00366 int num_conn; 00367 ErrorCode rval = MB_SUCCESS; 00368 std::vector< EntityHandle > storage; 00369 for( Range::iterator i = elements.begin(); i != elements.end(); ++i ) 00370 { 00371 // TODO: not generic enough. Why dim != 3 00372 // Commence un-necessary deep copying. 00373 const EntityHandle* connect; 00374 rval = mbImpl->get_connectivity( *i, connect, num_conn, false, &storage ); 00375 if( MB_SUCCESS != rval ) return rval; 00376 rval = mbImpl->get_coords( connect, num_conn, &coordinate[0] ); 00377 if( MB_SUCCESS != rval ) return rval; 00378 BoundBox box; 00379 for( int j = 0; j < num_conn; j++ ) 00380 box.update( &coordinate[3 * j] ); 00381 if( i == elements.begin() ) 00382 bounding_box = box; 00383 else 00384 bounding_box.update( box ); 00385 handle_data_vec.push_back( HandleData( *i, box, 0.0 ) ); 00386 } 00387 00388 return rval; 00389 } 00390 00391 inline ErrorCode BVHTree::reset_tree() 00392 { 00393 return delete_tree_sets(); 00394 } 00395 00396 } // namespace moab 00397 00398 #endif // BVH_TREE_HPP