MOAB: Mesh Oriented datABase  (version 5.2.1)
AdaptiveKDTree.hpp
Go to the documentation of this file.
00001 /**\file AdaptiveKDTree.hpp
00002  * \class moab::AdaptiveKDTree
00003  * \brief Adaptive KD tree, for sorting and searching entities spatially
00004  */
00005 
00006 #ifndef MOAB_ADAPTIVE_KD_TREE_HPP
00007 #define MOAB_ADAPTIVE_KD_TREE_HPP
00008 
00009 #include "moab/Types.hpp"
00010 #include "moab/Tree.hpp"
00011 
00012 #include <string>
00013 #include <vector>
00014 #include <cmath>
00015 
00016 namespace moab
00017 {
00018 
00019 class AdaptiveKDTreeIter;
00020 class Interface;
00021 class Range;
00022 
00023 class AdaptiveKDTree : public Tree
00024 {
00025   public:
00026     AdaptiveKDTree( Interface* iface );
00027 
00028     /** \brief Constructor (build the tree on construction)
00029      * Construct a tree object, and build the tree with entities input.  See comments
00030      * for build_tree() for detailed description of arguments.
00031      * \param iface MOAB instance
00032      * \param entities Entities to build tree around
00033      * \param tree_root Root set for tree (see function description)
00034      * \param opts Options for tree (see function description)
00035      */
00036     AdaptiveKDTree( Interface* iface, const Range& entities, EntityHandle* tree_root_set = NULL,
00037                     FileOptions* opts = NULL );
00038 
00039     ~AdaptiveKDTree();
00040 
00041     /** \brief Parse options for tree creation
00042      * \param options Options passed in by application
00043      * \return Failure is returned if any options were passed in and not interpreted; could mean
00044      * inappropriate options for a particular tree type
00045      */
00046     ErrorCode parse_options( FileOptions& options );
00047 
00048     /** Build the tree
00049      * Build a tree with the entities input.  If a non-NULL tree_root_set pointer is input,
00050      * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct
00051      * a new root set and pass its handle back in *tree_root_set.  Options vary by tree type;
00052      * see Tree.hpp for common options; options specific to AdaptiveKDTree:
00053      * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
00054      * PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
00055      *          for possible values; default = 1 (SUBDIVISION_SNAP)
00056      * \param entities Entities with which to build the tree
00057      * \param tree_root Root set for tree (see function description)
00058      * \param opts Options for tree (see function description)
00059      * \return Error is returned only on build failure
00060      */
00061     virtual ErrorCode build_tree( const Range& entities, EntityHandle* tree_root_set = NULL,
00062                                   FileOptions* options = NULL );
00063 
00064     //! Reset the tree, optionally checking we have the right root
00065     virtual ErrorCode reset_tree();
00066 
00067     /** \brief Get leaf containing input position.
00068      *
00069      * Does not take into account global bounding box of tree.
00070      * - Therefore there is always one leaf containing the point.
00071      * - If caller wants to account for global bounding box, then
00072      * caller can test against that box and not call this method
00073      * at all if the point is outside the box, as there is no leaf
00074      * containing the point in that case.
00075      * \param point Point to be located in tree
00076      * \param leaf_out Leaf containing point
00077      * \param iter_tol Tolerance for convergence of point search
00078      * \param inside_tol Tolerance for inside element calculation
00079      * \param multiple_leaves Some tree types can have multiple leaves containing a point;
00080      *          if non-NULL, this parameter is returned true if multiple leaves contain
00081      *          the input point
00082      * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
00083      * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
00084      */
00085     virtual ErrorCode point_search( const double* point, EntityHandle& leaf_out, const double iter_tol = 1.0e-10,
00086                                     const double inside_tol = 1.0e-6, bool* multiple_leaves = NULL,
00087                                     EntityHandle* start_node = NULL, CartVect* params = NULL );
00088 
00089     /** \brief Get leaf containing input position.
00090      *
00091      * Does not take into account global bounding box of tree.
00092      * - Therefore there is always one leaf containing the point.
00093      * - If caller wants to account for global bounding box, then
00094      * caller can test against that box and not call this method
00095      * at all if the point is outside the box, as there is no leaf
00096      * containing the point in that case.
00097      * \param point Point to be located in tree
00098      * \param leaf_it Iterator to leaf containing point
00099      * \param iter_tol Tolerance for convergence of point search
00100      * \param inside_tol Tolerance for inside element calculation
00101      * \param multiple_leaves Some tree types can have multiple leaves containing a point;
00102      *          if non-NULL, this parameter is returned true if multiple leaves contain
00103      *          the input point
00104      * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
00105      * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
00106      */
00107     ErrorCode point_search( const double* point, AdaptiveKDTreeIter& leaf_it, const double iter_tol = 1.0e-10,
00108                             const double inside_tol = 1.0e-6, bool* multiple_leaves = NULL,
00109                             EntityHandle* start_node = NULL );
00110 
00111     /** \brief Find all leaves within a given distance from point
00112      * If dists_out input non-NULL, also returns distances from each leaf; if
00113      * point i is inside leaf, 0 is given as dists_out[i].
00114      * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
00115      * in tree nodes and return containing entities in leaves_out.  In those cases, if params_out
00116      * is also non-NULL, will return parameters in those elements in that vector.
00117      * \param point Point to be located in tree
00118      * \param distance Distance within which to query
00119      * \param leaves_out Leaves within distance or containing point
00120      * \param iter_tol Tolerance for convergence of point search
00121      * \param inside_tol Tolerance for inside element calculation
00122      * \param dists_out If non-NULL, will contain distsances to leaves
00123      * \param params_out If non-NULL, will contain parameters of the point in the ents in leaves_out
00124      * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
00125      */
00126     virtual ErrorCode distance_search( const double* point, const double distance,
00127                                        std::vector< EntityHandle >& leaves_out, const double iter_tol = 1.0e-10,
00128                                        const double inside_tol = 1.0e-6, std::vector< double >* dists_out = NULL,
00129                                        std::vector< CartVect >* params_out = NULL, EntityHandle* start_node = NULL );
00130 
00131     ErrorCode get_info( EntityHandle root, double min[3], double max[3], unsigned int& dep );
00132 
00133     //! Enumeriate split plane directions
00134     enum Axis
00135     {
00136         X = 0,
00137         Y = 1,
00138         Z = 2
00139     };
00140 
00141     //! Split plane
00142     struct Plane
00143     {
00144         double coord;  //!< Location of plane as coordinate on normal axis
00145         int norm;      //!< The principal axis that is the normal of the plane;
00146 
00147         /** return true if point is below/to the left of the split plane */
00148         bool left_side( const double point[3] )
00149         {
00150             return point[norm] < coord;
00151         }
00152         /** return true if point is above/to the right of the split plane */
00153         bool right_side( const double point[3] )
00154         {
00155             return point[norm] > coord;
00156         }
00157         /** return distance from point to plane */
00158         double distance( const double point[3] ) const
00159         {
00160             return fabs( point[norm] - coord );
00161         }
00162     };
00163 
00164     //! Get split plane for tree node
00165     ErrorCode get_split_plane( EntityHandle node, Plane& plane );
00166 
00167     //! Set split plane for tree node
00168     ErrorCode set_split_plane( EntityHandle node, const Plane& plane );
00169 
00170     //! Get iterator for tree
00171     ErrorCode get_tree_iterator( EntityHandle tree_root, AdaptiveKDTreeIter& result );
00172 
00173     //! Get iterator at right-most ('last') leaf.
00174     ErrorCode get_last_iterator( EntityHandle tree_root, AdaptiveKDTreeIter& result );
00175 
00176     //! Get iterator for tree or subtree
00177     ErrorCode get_sub_tree_iterator( EntityHandle tree_root, const double box_min[3], const double box_max[3],
00178                                      AdaptiveKDTreeIter& result );
00179 
00180     //! Split leaf of tree
00181     //! Updates iterator location to point to first new leaf node.
00182     ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane );
00183 
00184     //! Split leaf of tree
00185     //! Updates iterator location to point to first new leaf node.
00186     ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, EntityHandle& left_child, EntityHandle& right_child );
00187     //! Split leaf of tree
00188     //! Updates iterator location to point to first new leaf node.
00189     ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, const Range& left_entities,
00190                           const Range& right_entities );
00191 
00192     //! Split leaf of tree
00193     //! Updates iterator location to point to first new leaf node.
00194     ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, const std::vector< EntityHandle >& left_entities,
00195                           const std::vector< EntityHandle >& right_entities );
00196 
00197     //! Merge the leaf pointed to by the current iterator with it's
00198     //! sibling.  If the sibling is not a leaf, multiple merges may
00199     //! be done.
00200     ErrorCode merge_leaf( AdaptiveKDTreeIter& iter );
00201 
00202     //! Find triangle closest to input position.
00203     //!\param from_coords  The input position to test against
00204     //!\param closest_point_out  The closest point on the set of triangles in the tree
00205     //!\param triangle_out The triangle closest to the input position
00206     ErrorCode closest_triangle( EntityHandle tree_root, const double from_coords[3], double closest_point_out[3],
00207                                 EntityHandle& triangle_out );
00208 
00209     ErrorCode sphere_intersect_triangles( EntityHandle tree_root, const double center[3], double radius,
00210                                           std::vector< EntityHandle >& triangles );
00211 
00212     ErrorCode ray_intersect_triangles( EntityHandle tree_root, const double tolerance, const double ray_unit_dir[3],
00213                                        const double ray_base_pt[3], std::vector< EntityHandle >& triangles_out,
00214                                        std::vector< double >& distance_out, int result_count_limit = 0,
00215                                        double distance_limit = -1.0 );
00216 
00217     ErrorCode compute_depth( EntityHandle root, unsigned int& min_depth, unsigned int& max_depth );
00218 
00219     //! methods for selecting candidate split planes
00220     enum CandidatePlaneSet
00221     {
00222         //! Candidiate planes at evenly spaced intervals
00223         SUBDIVISION = 0,
00224         //! Like SUBDIVISION, except snap to closest vertex coordinate
00225         SUBDIVISION_SNAP,  // = 1
00226                            //! Median vertex coodinate values
00227         VERTEX_MEDIAN,     // = 2
00228                            //! Random sampling of vertex coordinate values
00229         VERTEX_SAMPLE      // = 3
00230     };
00231 
00232     //! print various things about this tree
00233     virtual ErrorCode print();
00234 
00235   private:
00236     friend class AdaptiveKDTreeIter;
00237 
00238     ErrorCode init();
00239 
00240     /**\brief find a triangle near the input point */
00241     ErrorCode find_close_triangle( EntityHandle root, const double from_point[3], double pt[3],
00242                                    EntityHandle& triangle );
00243 
00244     ErrorCode make_tag( Interface* iface, std::string name, TagType storage, DataType type, int count,
00245                         void* default_val, Tag& tag_handle, std::vector< Tag >& created_tags );
00246 
00247     ErrorCode intersect_children_with_elems( const Range& elems, AdaptiveKDTree::Plane plane, double eps,
00248                                              CartVect box_min, CartVect box_max, Range& left_tris, Range& right_tris,
00249                                              Range& both_tris, double& metric_value );
00250 
00251     ErrorCode best_subdivision_snap_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
00252                                            Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
00253                                            std::vector< double >& tmp_data, double eps );
00254 
00255     ErrorCode best_subdivision_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
00256                                       Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
00257                                       double eps );
00258 
00259     ErrorCode best_vertex_median_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
00260                                         Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
00261                                         std::vector< double >& coords, double eps );
00262 
00263     ErrorCode best_vertex_sample_plane( int num_planes, const AdaptiveKDTreeIter& iter, Range& best_left,
00264                                         Range& best_right, Range& best_both, AdaptiveKDTree::Plane& best_plane,
00265                                         std::vector< double >& coords, std::vector< EntityHandle >& indices,
00266                                         double eps );
00267 
00268     static const char* treeName;
00269 
00270     Tag planeTag, axisTag;
00271 
00272     unsigned splitsPerDir;
00273 
00274     CandidatePlaneSet planeSet;
00275 
00276     bool spherical;
00277     double radius;
00278 };
00279 
00280 //! Iterate over leaves of an adaptive kD-tree
00281 class AdaptiveKDTreeIter
00282 {
00283   public:
00284     enum Direction
00285     {
00286         LEFT  = 0,
00287         RIGHT = 1
00288     };
00289 
00290   private:
00291     struct StackObj
00292     {
00293         StackObj( EntityHandle e, double c ) : entity( e ), coord( c ) {}
00294         StackObj() : entity( 0 ), coord( 0.0 ) {}
00295         EntityHandle entity;  //!< handle for tree node
00296         double coord;         //!< box coordinate of parent
00297     };
00298 
00299     enum
00300     {
00301         BMIN = 0,
00302         BMAX = 1
00303     };  //!< indices into mBox and child list
00304 
00305     CartVect mBox[2];                               //!< min and max corners of bounding box
00306     AdaptiveKDTree* treeTool;                       //!< tool for tree
00307     std::vector< StackObj > mStack;                 //!< stack storing path through tree
00308     mutable std::vector< EntityHandle > childVect;  //!< temporary storage of child handles
00309 
00310     //! Descend tree to left most leaf from current position
00311     //! No-op if at leaf.
00312     ErrorCode step_to_first_leaf( Direction direction );
00313 
00314     friend class AdaptiveKDTree;
00315 
00316   public:
00317     AdaptiveKDTreeIter() : treeTool( 0 ), childVect( 2 ) {}
00318 
00319     ErrorCode initialize( AdaptiveKDTree* tool, EntityHandle root, const double box_min[3], const double box_max[3],
00320                           Direction direction );
00321 
00322     AdaptiveKDTree* tool() const
00323     {
00324         return treeTool;
00325     }
00326 
00327     //! Get handle for current leaf
00328     EntityHandle handle() const
00329     {
00330         return mStack.back().entity;
00331     }
00332 
00333     //! Get min corner of axis-aligned box for current leaf
00334     const double* box_min() const
00335     {
00336         return mBox[BMIN].array();
00337     }
00338 
00339     //! Get max corner of axis-aligned box for current leaf
00340     const double* box_max() const
00341     {
00342         return mBox[BMAX].array();
00343     }
00344 
00345     double volume() const
00346     {
00347         return ( mBox[BMAX][0] - mBox[BMIN][0] ) * ( mBox[BMAX][1] - mBox[BMIN][1] ) *
00348                ( mBox[BMAX][2] - mBox[BMIN][2] );
00349     }
00350 
00351     //! test if a plane intersects the leaf box
00352     bool intersects( const AdaptiveKDTree::Plane& plane ) const
00353     {
00354         return mBox[BMIN][plane.norm] <= plane.coord && mBox[BMAX][plane.norm] >= plane.coord;
00355     }
00356 
00357     //! Get depth in tree. root is at depth of 1.
00358     unsigned depth() const
00359     {
00360         return mStack.size();
00361     }
00362 
00363     //! Advance the iterator either left or right in the tree
00364     //! Note:  stepping past the end of the tree will invalidate
00365     //!        the iterator.  It will *not* be work step the
00366     //!        other direction.
00367     ErrorCode step( Direction direction );
00368 
00369     //! Advance to next leaf
00370     //! Returns MB_ENTITY_NOT_FOUND if at end.
00371     //! Note: steping past the end of the tree will invalidate
00372     //!       the iterator. Calling back() will not work.
00373     ErrorCode step()
00374     {
00375         return step( RIGHT );
00376     }
00377 
00378     //! Move back to previous leaf
00379     //! Returns MB_ENTITY_NOT_FOUND if at beginning.
00380     //! Note: steping past the start of the tree will invalidate
00381     //!       the iterator. Calling step() will not work.
00382     ErrorCode back()
00383     {
00384         return step( LEFT );
00385     }
00386 
00387     //! Return the side of the box bounding this tree node
00388     //! that is shared with the immediately adjacent sibling
00389     //! (the tree node that shares a common parent node with
00390     //! this node in the binary tree.)
00391     //!
00392     //!\param axis_out The principal axis orthogonal to the side of the box
00393     //!\param neg_out  true if the side of the box is toward the decreasing
00394     //!                direction of the principal axis indicated by axis_out,
00395     //!                false if it is toward the increasing direction.
00396     //!\return MB_ENTITY_NOT FOUND if root node.
00397     //!        MB_FAILURE if internal error.
00398     //!        MB_SUCCESS otherwise.
00399     ErrorCode sibling_side( AdaptiveKDTree::Axis& axis_out, bool& neg_out ) const;
00400 
00401     //! Get adjacent leaf nodes on side indicated by norm and neg.
00402     //!
00403     //! E.g. if norm == X and neg == true, then get neighbor(s)
00404     //! adjacent to the side of the box contained in the plane
00405     //! with normal to the X axis and with the x coordinate equal
00406     //! to the minimum x of the bounding box.
00407     //!
00408     //! E.g. if norm == Y and neg == false, then get neighbor(s)
00409     //! adjacent to the side of the box with y = maximum y of bounding box.
00410     //!
00411     //!\param norm  Normal vector for box side (X, Y, or Z)
00412     //!\param neg   Which of two planes with norm (true->smaller coord,
00413     //!             false->larger coord)
00414     //!\param results List to which to append results.  This function does
00415     //!             *not* clear existing values in list.
00416     //!\param epsilon Tolerance on overlap.  A positive value E will
00417     //!              result in nodes that are separated by as much as E
00418     //!              to be considered touching.  A negative value -E will
00419     //!              cause leaves that do not overlap by at least E to be
00420     //!              considered non-overlapping.  Amongst other things,
00421     //!              this value can be used to control whether or not
00422     //!              leaves adjacent at only their edges or corners are
00423     //!              returned.
00424     ErrorCode get_neighbors( AdaptiveKDTree::Axis norm, bool neg, std::vector< AdaptiveKDTreeIter >& results,
00425                              double epsilon = 0.0 ) const;
00426 
00427     //! Get split plane that separates this node from its immediate sibling.
00428     ErrorCode get_parent_split_plane( AdaptiveKDTree::Plane& plane ) const;
00429 
00430     //! Return true if thos node and the passed node share the
00431     //! same immediate parent.
00432     bool is_sibling( const AdaptiveKDTreeIter& other_leaf ) const;
00433 
00434     //! Return true if thos node and the passed node share the
00435     //! same immediate parent.
00436     bool is_sibling( EntityHandle other_leaf ) const;
00437 
00438     //! Returns true if calling step() will advance to the
00439     //! immediate sibling of the current node.  Returns false
00440     //! if current node is root or back() will move to the
00441     //! immediate sibling.
00442     bool sibling_is_forward() const;
00443 
00444     //! Find range of overlap between ray and leaf.
00445     //!
00446     //!\param ray_point Coordinates of start point of ray
00447     //!\param ray_vect  Directionion vector for ray such that
00448     //!                 the ray is defined by r(t) = ray_point + t * ray_vect
00449     //!                 for t > 0.
00450     //!\param t_enter   Output: if return value is true, this value
00451     //!                 is the parameter location along the ray at which
00452     //!                 the ray entered the leaf.  If return value is false,
00453     //!                 then this value is undefined.
00454     //!\param t_exit    Output: if return value is true, this value
00455     //!                 is the parameter location along the ray at which
00456     //!                 the ray exited the leaf.  If return value is false,
00457     //!                 then this value is undefined.
00458     //!\return true if ray intersects leaf, false otherwise.
00459     bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter, double& t_exit ) const;
00460 };
00461 
00462 inline ErrorCode AdaptiveKDTree::reset_tree()
00463 {
00464     return delete_tree_sets();
00465 }
00466 
00467 }  // namespace moab
00468 
00469 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines