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