Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
BVHTree.hpp
Go to the documentation of this file.
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,
00070                                   EntityHandle* tree_root_set = NULL,
00071                                   FileOptions* options        = NULL );
00072 
00073     /** \brief Get leaf containing input position.
00074      *
00075      * Does not take into account global bounding box of tree.
00076      * - Therefore there is always one leaf containing the point.
00077      * - If caller wants to account for global bounding box, then
00078      * caller can test against that box and not call this method
00079      * at all if the point is outside the box, as there is no leaf
00080      * containing the point in that case.
00081      * \param point Point to be located in tree
00082      * \param leaf_out Leaf containing point
00083      * \param iter_tol Tolerance for convergence of point search
00084      * \param inside_tol Tolerance for inside element calculation
00085      * \param multiple_leaves Some tree types can have multiple leaves containing a point;
00086      *          if non-NULL, this parameter is returned true if multiple leaves contain
00087      *          the input point
00088      * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
00089      * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
00090      */
00091     virtual ErrorCode point_search( const double* point,
00092                                     EntityHandle& leaf_out,
00093                                     const double iter_tol    = 1.0e-10,
00094                                     const double inside_tol  = 1.0e-6,
00095                                     bool* multiple_leaves    = NULL,
00096                                     EntityHandle* start_node = NULL,
00097                                     CartVect* params         = NULL );
00098 
00099     /** \brief Find all leaves within a given distance from point
00100      * If dists_out input non-NULL, also returns distances from each leaf; if
00101      * point i is inside leaf, 0 is given as dists_out[i].
00102      * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
00103      * in tree nodes and return containing entities in leaves_out.  In those cases, if params_out
00104      * is also non-NULL, will return parameters in those elements in that vector.
00105      * \param from_point Point to be located in tree
00106      * \param distance Distance within which to query
00107      * \param result_list Leaves within distance or containing point
00108      * \param iter_tol Tolerance for convergence of point search
00109      * \param inside_tol Tolerance for inside element calculation
00110      * \param result_dists If non-NULL, will contain distsances to leaves
00111      * \param result_params If non-NULL, will contain parameters of the point in the ents in
00112      * leaves_out \param tree_root Start from this tree node (non-NULL) instead of tree root (NULL)
00113      */
00114     virtual ErrorCode distance_search( const double from_point[3],
00115                                        const double distance,
00116                                        std::vector< EntityHandle >& result_list,
00117                                        const double iter_tol                  = 1.0e-10,
00118                                        const double inside_tol                = 1.0e-6,
00119                                        std::vector< double >* result_dists    = NULL,
00120                                        std::vector< CartVect >* result_params = NULL,
00121                                        EntityHandle* tree_root                = NULL );
00122 
00123     //! print various things about this tree
00124     virtual ErrorCode print();
00125 
00126   private:
00127     // don't allow copy constructor, too complicated
00128     BVHTree( const BVHTree& s );
00129 
00130     class HandleData
00131     {
00132       public:
00133         HandleData( EntityHandle h, const BoundBox& bx, const double dp ) : myHandle( h ), myBox( bx ), myDim( dp ) {}
00134         HandleData() : myHandle( 0 ), myDim( -1 ) {}
00135         EntityHandle myHandle;
00136         BoundBox myBox;
00137         double myDim;
00138     };
00139     typedef std::vector< HandleData > HandleDataVec;
00140 
00141     class SplitData
00142     {
00143       public:
00144         SplitData()
00145             : dim( UINT_MAX ), nl( UINT_MAX ), nr( UINT_MAX ), split( DBL_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX )
00146         {
00147         }
00148         SplitData( const SplitData& f )
00149             : dim( f.dim ), nl( f.nl ), nr( f.nr ), split( f.split ), Lmax( f.Lmax ), Rmin( f.Rmin ),
00150               boundingBox( f.boundingBox ), leftBox( f.leftBox ), rightBox( f.rightBox )
00151         {
00152         }
00153         unsigned int dim, nl, nr;
00154         double split;
00155         double Lmax, Rmin;
00156         BoundBox boundingBox, leftBox, rightBox;
00157         SplitData& operator=( const SplitData& f )
00158         {
00159             dim         = f.dim;
00160             nl          = f.nl;
00161             nr          = f.nr;
00162             split       = f.split;
00163             Lmax        = f.Lmax;
00164             Rmin        = f.Rmin;
00165             boundingBox = f.boundingBox;
00166             leftBox     = f.leftBox;
00167             rightBox    = f.rightBox;
00168             return *this;
00169         }
00170     };  // SplitData
00171 
00172     class Split_comparator
00173     {
00174         // deprecation of binary_function
00175         typedef SplitData first_argument_type;
00176         typedef SplitData second_argument_type;
00177         typedef bool result_type;
00178 
00179         inline double objective( const SplitData& a ) const
00180         {
00181             return a.Lmax * a.nl - a.Rmin * a.nr;
00182         }
00183 
00184       public:
00185         bool operator()( const SplitData& a, const SplitData& b ) const
00186         {
00187             return objective( a ) < objective( b );
00188         }
00189     };  // Split_comparator
00190 
00191     class HandleData_comparator
00192     {
00193         // deprecation of binary_function
00194         typedef HandleData first_argument_type;
00195         typedef HandleData second_argument_type;
00196         typedef bool result_type;
00197 
00198       public:
00199         bool operator()( const HandleData& a, const HandleData& b )
00200         {
00201             return a.myDim < b.myDim || ( a.myDim == b.myDim && a.myHandle < b.myHandle );
00202         }
00203     };  // Iterator_comparator
00204 
00205     class Bucket
00206     {
00207       public:
00208         Bucket() : mySize( 0 ) {}
00209         Bucket( const Bucket& f ) : mySize( f.mySize ), boundingBox( f.boundingBox ) {}
00210         Bucket( const unsigned int sz ) : mySize( sz ) {}
00211         static unsigned int bucket_index( int num_splits,
00212                                           const BoundBox& box,
00213                                           const BoundBox& interval,
00214                                           const unsigned int dim );
00215         unsigned int mySize;
00216         BoundBox boundingBox;
00217         Bucket& operator=( const Bucket& f )
00218         {
00219             boundingBox = f.boundingBox;
00220             mySize      = f.mySize;
00221             return *this;
00222         }
00223     };  // Bucket
00224 
00225     class Node
00226     {
00227       public:
00228         HandleDataVec entities;
00229         unsigned int dim, child;
00230         double Lmax, Rmin;
00231         BoundBox box;
00232         Node() : dim( UINT_MAX ), child( UINT_MAX ), Lmax( -DBL_MAX ), Rmin( DBL_MAX ) {}
00233         Node( const Node& f ) : entities( f.entities ), dim( f.dim ), child( f.child ), Lmax( f.Lmax ), Rmin( f.Rmin )
00234         {
00235         }
00236         Node& operator=( const Node& f )
00237         {
00238             dim      = f.dim;
00239             child    = f.child;
00240             Lmax     = f.Lmax;
00241             Rmin     = f.Rmin;
00242             entities = f.entities;
00243             return *this;
00244         }
00245     };  // Node
00246 
00247     class TreeNode
00248     {
00249       public:
00250         unsigned int dim, child;
00251         double Lmax, Rmin;
00252         BoundBox box;
00253         TreeNode( int dm, int chld, double lmx, double rmn, BoundBox& bx )
00254             : dim( dm ), child( chld ), Lmax( lmx ), Rmin( rmn ), box( bx )
00255         {
00256         }
00257         TreeNode( const TreeNode& f )
00258         {
00259             dim   = f.dim;
00260             child = f.child;
00261             Lmax  = f.Lmax;
00262             Rmin  = f.Rmin;
00263         }
00264         TreeNode& operator=( const TreeNode& f )
00265         {
00266             dim   = f.dim;
00267             child = f.child;
00268             Lmax  = f.Lmax;
00269             Rmin  = f.Rmin;
00270             return *this;
00271         }
00272     };  // TreeNode
00273 
00274     void establish_buckets( HandleDataVec::const_iterator begin,
00275                             HandleDataVec::const_iterator end,
00276                             const BoundBox& interval,
00277                             std::vector< std::vector< Bucket > >& buckets ) const;
00278 
00279     unsigned int set_interval( BoundBox& interval,
00280                                std::vector< Bucket >::const_iterator begin,
00281                                std::vector< Bucket >::const_iterator end ) const;
00282 
00283     void initialize_splits( std::vector< std::vector< SplitData > >& splits,
00284                             const std::vector< std::vector< Bucket > >& buckets,
00285                             const SplitData& data ) const;
00286 
00287     void order_elements( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
00288 
00289     void median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
00290 
00291     void choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const;
00292 
00293     void find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const;
00294 
00295     ErrorCode find_point( const std::vector< double >& point,
00296                           const unsigned int& index,
00297                           const double iter_tol,
00298                           const double inside_tol,
00299                           std::pair< EntityHandle, CartVect >& result );
00300 
00301     EntityHandle bruteforce_find( const double* point,
00302                                   const double iter_tol   = 1.0e-10,
00303                                   const double inside_tol = 1.0e-6 );
00304 
00305     int local_build_tree( std::vector< Node >& tree_nodes,
00306                           HandleDataVec::iterator begin,
00307                           HandleDataVec::iterator end,
00308                           const int index,
00309                           const BoundBox& box,
00310                           const int depth = 0 );
00311 
00312     // builds up vector of HandleData, which caches elements' bounding boxes
00313     ErrorCode construct_element_vec( std::vector< HandleData >& handle_data_vec,
00314                                      const Range& elements,
00315                                      BoundBox& bounding_box );
00316 
00317     // convert the std::vector<Node> to myTree and a bunch of entity sets
00318     ErrorCode convert_tree( std::vector< Node >& tree_nodes );
00319 
00320     // print tree nodes
00321     ErrorCode print_nodes( std::vector< Node >& nodes );
00322 
00323     Range entityHandles;
00324     std::vector< TreeNode > myTree;
00325     int splitsPerDir;
00326     EntityHandle startSetHandle;
00327     static MOAB_EXPORT const char* treeName;
00328 };  // class Bvh_tree
00329 
00330 inline unsigned int BVHTree::Bucket::bucket_index( int num_splits,
00331                                                    const BoundBox& box,
00332                                                    const BoundBox& interval,
00333                                                    const unsigned int dim )
00334 {
00335     // see FastMemoryEfficientCellLocationinUnstructuredGridsForVisualization.pdf
00336     // around page 9
00337 
00338     // Paper arithmetic is over-optimized.. this is safer.
00339     const double min    = interval.bMin[dim];
00340     const double length = ( interval.bMax[dim] - min ) / ( num_splits + 1 );
00341     const double center = ( ( box.bMax[dim] + box.bMin[dim] ) / 2.0 ) - min;
00342 #ifndef NDEBUG
00343 #ifdef BVH_SHOW_INDEX
00344     std::cout << "[" << min << " , " << interval.max[dim] << " ]" << std::endl;
00345     std::cout << "[" << box.bMin[dim] << " , " << box.bMax[dim] << " ]" << std::endl;
00346     std::cout << "Length of bucket" << length << std::endl;
00347     std::cout << "Center: " << ( box.bMax[dim] + box.bMin[dim] ) / 2.0 << std::endl;
00348     std::cout << "Distance of center from min:  " << center << std::endl;
00349     std::cout << "ratio: " << center / length << std::endl;
00350     std::cout << "index: " << std::ceil( center / length ) - 1 << std::endl;
00351 #endif
00352 #endif
00353     unsigned int cl = std::ceil( center / length );
00354     return ( cl > 0 ? cl - 1 : 0 );
00355 }
00356 
00357 inline BVHTree::BVHTree( Interface* impl ) : Tree( impl ), splitsPerDir( 3 ), startSetHandle( 0 )
00358 {
00359     boxTagName = treeName;
00360 }
00361 
00362 inline unsigned int BVHTree::set_interval( BoundBox& interval,
00363                                            std::vector< Bucket >::const_iterator begin,
00364                                            std::vector< Bucket >::const_iterator end ) const
00365 {
00366     bool first         = true;
00367     unsigned int count = 0;
00368     for( std::vector< Bucket >::const_iterator b = begin; b != end; ++b )
00369     {
00370         const BoundBox& box = b->boundingBox;
00371         count += b->mySize;
00372         if( b->mySize != 0 )
00373         {
00374             if( first )
00375             {
00376                 interval = box;
00377                 first    = false;
00378             }
00379             else
00380                 interval.update( box );
00381         }
00382     }
00383     return count;
00384 }
00385 
00386 inline void BVHTree::order_elements( HandleDataVec::iterator& begin,
00387                                      HandleDataVec::iterator& end,
00388                                      SplitData& data ) const
00389 {
00390     for( HandleDataVec::iterator i = begin; i != end; ++i )
00391     {
00392         const int index = Bucket::bucket_index( splitsPerDir, i->myBox, data.boundingBox, data.dim );
00393         i->myDim        = ( index <= data.split ) ? 0 : 1;
00394     }
00395     std::sort( begin, end, HandleData_comparator() );
00396 }
00397 
00398 inline void BVHTree::choose_best_split( const std::vector< std::vector< SplitData > >& splits, SplitData& data ) const
00399 {
00400     std::vector< SplitData > best_splits;
00401     for( std::vector< std::vector< SplitData > >::const_iterator i = splits.begin(); i != splits.end(); ++i )
00402     {
00403         std::vector< SplitData >::const_iterator j =
00404             std::min_element( ( *i ).begin(), ( *i ).end(), Split_comparator() );
00405         best_splits.push_back( *j );
00406     }
00407     data = *std::min_element( best_splits.begin(), best_splits.end(), Split_comparator() );
00408 }
00409 
00410 inline ErrorCode BVHTree::construct_element_vec( std::vector< HandleData >& handle_data_vec,
00411                                                  const Range& elements,
00412                                                  BoundBox& bounding_box )
00413 {
00414     std::vector< double > coordinate( 3 * CN::MAX_NODES_PER_ELEMENT );
00415     int num_conn;
00416     ErrorCode rval = MB_SUCCESS;
00417     std::vector< EntityHandle > storage;
00418     for( Range::iterator i = elements.begin(); i != elements.end(); ++i )
00419     {
00420         // TODO: not generic enough. Why dim != 3
00421         // Commence un-necessary deep copying.
00422         const EntityHandle* connect;
00423         rval = mbImpl->get_connectivity( *i, connect, num_conn, false, &storage );
00424         if( MB_SUCCESS != rval ) return rval;
00425         rval = mbImpl->get_coords( connect, num_conn, &coordinate[0] );
00426         if( MB_SUCCESS != rval ) return rval;
00427         BoundBox box;
00428         for( int j = 0; j < num_conn; j++ )
00429             box.update( &coordinate[3 * j] );
00430         if( i == elements.begin() )
00431             bounding_box = box;
00432         else
00433             bounding_box.update( box );
00434         handle_data_vec.push_back( HandleData( *i, box, 0.0 ) );
00435     }
00436 
00437     return rval;
00438 }
00439 
00440 inline ErrorCode BVHTree::reset_tree()
00441 {
00442     return delete_tree_sets();
00443 }
00444 
00445 }  // namespace moab
00446 
00447 #endif  // BVH_TREE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines