AdaptiveKDTree.hpp

Go to the documentation of this file.

/**\file AdaptiveKDTree.hpp
 * \class moab::AdaptiveKDTree
 * \brief Adaptive KD tree, for sorting and searching entities spatially
 */
 
#ifndef MOAB_ADAPTIVE_KD_TREE_HPP
#define MOAB_ADAPTIVE_KD_TREE_HPP
 
#include "moab/Types.hpp"
#include "moab/Tree.hpp"
 
#include <string>
#include <vector>
#include <cmath>
 
namespace moab
{
 
class AdaptiveKDTreeIter;
class Interface;
class Range;
 
class AdaptiveKDTree : public Tree
{
 public:
 AdaptiveKDTree( Interface* iface );
 
 /** \brief Constructor (build the tree on construction)
 * Construct a tree object, and build the tree with entities input. See comments
 * for build_tree() for detailed description of arguments.
 * \param iface MOAB instance
 * \param entities Entities to build tree around
 * \param tree_root Root set for tree (see function description)
 * \param opts Options for tree (see function description)
 */
 AdaptiveKDTree( Interface* iface,
 const Range& entities,
 EntityHandle* tree_root_set = NULL,
 FileOptions* opts = NULL );
 
 ~AdaptiveKDTree();
 
 /** \brief Parse options for tree creation
 * \param options Options passed in by application
 * \return Failure is returned if any options were passed in and not interpreted; could mean
 * inappropriate options for a particular tree type
 */
 ErrorCode parse_options( FileOptions& options );
 
 /** Build the tree
 * Build a tree with the entities input. If a non-NULL tree_root_set pointer is input,
 * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct
 * a new root set and pass its handle back in *tree_root_set. Options vary by tree type;
 * see Tree.hpp for common options; options specific to AdaptiveKDTree:
 * SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
 * PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
 * for possible values; default = 1 (SUBDIVISION_SNAP)
 * \param entities Entities with which to build the tree
 * \param tree_root Root set for tree (see function description)
 * \param opts Options for tree (see function description)
 * \return Error is returned only on build failure
 */
 virtual ErrorCode build_tree( const Range& entities,
 EntityHandle* tree_root_set = NULL,
 FileOptions* options = NULL );
 
 //! Reset the tree, optionally checking we have the right root
 virtual ErrorCode reset_tree();
 
 /** \brief Get leaf containing input position.
 *
 * Does not take into account global bounding box of tree.
 * - Therefore there is always one leaf containing the point.
 * - If caller wants to account for global bounding box, then
 * caller can test against that box and not call this method
 * at all if the point is outside the box, as there is no leaf
 * containing the point in that case.
 * \param point Point to be located in tree
 * \param leaf_out Leaf containing point
 * \param iter_tol Tolerance for convergence of point search
 * \param inside_tol Tolerance for inside element calculation
 * \param multiple_leaves Some tree types can have multiple leaves containing a point;
 * if non-NULL, this parameter is returned true if multiple leaves contain
 * the input point
 * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
 * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
 */
 virtual ErrorCode point_search( const double* point,
 EntityHandle& leaf_out,
 const double iter_tol = 1.0e-10,
 const double inside_tol = 1.0e-6,
 bool* multiple_leaves = NULL,
 EntityHandle* start_node = NULL,
 CartVect* params = NULL );
 
 /** \brief Get leaf containing input position.
 *
 * Does not take into account global bounding box of tree.
 * - Therefore there is always one leaf containing the point.
 * - If caller wants to account for global bounding box, then
 * caller can test against that box and not call this method
 * at all if the point is outside the box, as there is no leaf
 * containing the point in that case.
 * \param point Point to be located in tree
 * \param leaf_it Iterator to leaf containing point
 * \param iter_tol Tolerance for convergence of point search
 * \param inside_tol Tolerance for inside element calculation
 * \param multiple_leaves Some tree types can have multiple leaves containing a point;
 * if non-NULL, this parameter is returned true if multiple leaves contain
 * the input point
 * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
 * \return Non-success returned only in case of failure; not-found indicated by leaf_out=0
 */
 ErrorCode point_search( const double* point,
 AdaptiveKDTreeIter& leaf_it,
 const double iter_tol = 1.0e-10,
 const double inside_tol = 1.0e-6,
 bool* multiple_leaves = NULL,
 EntityHandle* start_node = NULL );
 
 /** \brief Find all leaves within a given distance from point
 * If dists_out input non-NULL, also returns distances from each leaf; if
 * point i is inside leaf, 0 is given as dists_out[i].
 * If params_out is non-NULL and myEval is non-NULL, will evaluate individual entities
 * in tree nodes and return containing entities in leaves_out. In those cases, if params_out
 * is also non-NULL, will return parameters in those elements in that vector.
 * \param point Point to be located in tree
 * \param distance Distance within which to query
 * \param leaves_out Leaves within distance or containing point
 * \param iter_tol Tolerance for convergence of point search
 * \param inside_tol Tolerance for inside element calculation
 * \param dists_out If non-NULL, will contain distsances to leaves
 * \param params_out If non-NULL, will contain parameters of the point in the ents in leaves_out
 * \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
 */
 virtual ErrorCode distance_search( const double* point,
 const double distance,
 std::vector< EntityHandle >& leaves_out,
 const double iter_tol = 1.0e-10,
 const double inside_tol = 1.0e-6,
 std::vector< double >* dists_out = NULL,
 std::vector< CartVect >* params_out = NULL,
 EntityHandle* start_node = NULL );
 
 ErrorCode get_info( EntityHandle root, double min[3], double max[3], unsigned int& dep );
 
 //! Enumeriate split plane directions
 enum Axis
 {
 X = 0,
 Y = 1,
 Z = 2
 };
 
 //! Split plane
 struct Plane
 {
 double coord; //!< Location of plane as coordinate on normal axis
 int norm; //!< The principal axis that is the normal of the plane;
 
 /** return true if point is below/to the left of the split plane */
 bool left_side( const double point[3] )
 {
 return point[norm] < coord;
 }
 /** return true if point is above/to the right of the split plane */
 bool right_side( const double point[3] )
 {
 return point[norm] > coord;
 }
 /** return distance from point to plane */
 double distance( const double point[3] ) const
 {
 return fabs( point[norm] - coord );
 }
 };
 
 //! Get split plane for tree node
 ErrorCode get_split_plane( EntityHandle node, Plane& plane );
 
 //! Set split plane for tree node
 ErrorCode set_split_plane( EntityHandle node, const Plane& plane );
 
 //! Get iterator for tree
 ErrorCode get_tree_iterator( EntityHandle tree_root, AdaptiveKDTreeIter& result );
 
 //! Get iterator at right-most ('last') leaf.
 ErrorCode get_last_iterator( EntityHandle tree_root, AdaptiveKDTreeIter& result );
 
 //! Get iterator for tree or subtree
 ErrorCode get_sub_tree_iterator( EntityHandle tree_root,
 const double box_min[3],
 const double box_max[3],
 AdaptiveKDTreeIter& result );
 
 //! Split leaf of tree
 //! Updates iterator location to point to first new leaf node.
 ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane );
 
 //! Split leaf of tree
 //! Updates iterator location to point to first new leaf node.
 ErrorCode split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, EntityHandle& left_child, EntityHandle& right_child );
 //! Split leaf of tree
 //! Updates iterator location to point to first new leaf node.
 ErrorCode split_leaf( AdaptiveKDTreeIter& leaf,
 Plane plane,
 const Range& left_entities,
 const Range& right_entities );
 
 //! Split leaf of tree
 //! Updates iterator location to point to first new leaf node.
 ErrorCode split_leaf( AdaptiveKDTreeIter& leaf,
 Plane plane,
 const std::vector< EntityHandle >& left_entities,
 const std::vector< EntityHandle >& right_entities );
 
 //! Merge the leaf pointed to by the current iterator with it's
 //! sibling. If the sibling is not a leaf, multiple merges may
 //! be done.
 ErrorCode merge_leaf( AdaptiveKDTreeIter& iter );
 
 //! Find triangle closest to input position.
 //!\param from_coords The input position to test against
 //!\param closest_point_out The closest point on the set of triangles in the tree
 //!\param triangle_out The triangle closest to the input position
 ErrorCode closest_triangle( EntityHandle tree_root,
 const double from_coords[3],
 double closest_point_out[3],
 EntityHandle& triangle_out );
 
 ErrorCode sphere_intersect_triangles( EntityHandle tree_root,
 const double center[3],
 double radius,
 std::vector< EntityHandle >& triangles );
 
 ErrorCode ray_intersect_triangles( EntityHandle tree_root,
 const double tolerance,
 const double ray_unit_dir[3],
 const double ray_base_pt[3],
 std::vector< EntityHandle >& triangles_out,
 std::vector< double >& distance_out,
 int result_count_limit = 0,
 double distance_limit = -1.0 );
 
 ErrorCode compute_depth( EntityHandle root, unsigned int& min_depth, unsigned int& max_depth );
 
 //! methods for selecting candidate split planes
 enum CandidatePlaneSet
 {
 //! Candidiate planes at evenly spaced intervals
 SUBDIVISION = 0,
 //! Like SUBDIVISION, except snap to closest vertex coordinate
 SUBDIVISION_SNAP, // = 1
 //! Median vertex coodinate values
 VERTEX_MEDIAN, // = 2
 //! Random sampling of vertex coordinate values
 VERTEX_SAMPLE // = 3
 };
 
 //! print various things about this tree
 virtual ErrorCode print();
 
 private:
 friend class AdaptiveKDTreeIter;
 
 ErrorCode init();
 
 /**\brief find a triangle near the input point */
 ErrorCode find_close_triangle( EntityHandle root,
 const double from_point[3],
 double pt[3],
 EntityHandle& triangle );
 
 ErrorCode make_tag( Interface* iface,
 std::string name,
 TagType storage,
 DataType type,
 int count,
 void* default_val,
 Tag& tag_handle,
 std::vector< Tag >& created_tags );
 
 ErrorCode intersect_children_with_elems( const Range& elems,
 AdaptiveKDTree::Plane plane,
 double eps,
 CartVect box_min,
 CartVect box_max,
 Range& left_tris,
 Range& right_tris,
 Range& both_tris,
 double& metric_value );
 
 ErrorCode best_subdivision_snap_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 std::vector< double >& tmp_data,
 double eps );
 
 ErrorCode best_subdivision_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 double eps );
 
 ErrorCode best_vertex_median_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 std::vector< double >& coords,
 double eps );
 
 ErrorCode best_vertex_sample_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 std::vector< double >& coords,
 std::vector< EntityHandle >& indices,
 double eps );
 
 static const char* treeName;
 
 Tag planeTag, axisTag;
 
 unsigned splitsPerDir;
 
 CandidatePlaneSet planeSet;
 
 bool spherical;
 double radius;
};
 
//! Iterate over leaves of an adaptive kD-tree
class AdaptiveKDTreeIter
{
 public:
 enum Direction
 {
 LEFT = 0,
 RIGHT = 1
 };
 
 private:
 struct StackObj
 {
 StackObj( EntityHandle e, double c ) : entity( e ), coord( c ) {}
 StackObj() : entity( 0 ), coord( 0.0 ) {}
 EntityHandle entity; //!< handle for tree node
 double coord; //!< box coordinate of parent
 };
 
 enum
 {
 BMIN = 0,
 BMAX = 1
 }; //!< indices into mBox and child list
 
 CartVect mBox[2]; //!< min and max corners of bounding box
 AdaptiveKDTree* treeTool; //!< tool for tree
 std::vector< StackObj > mStack; //!< stack storing path through tree
 mutable std::vector< EntityHandle > childVect; //!< temporary storage of child handles
 
 //! Descend tree to left most leaf from current position
 //! No-op if at leaf.
 ErrorCode step_to_first_leaf( Direction direction );
 
 friend class AdaptiveKDTree;
 
 public:
 AdaptiveKDTreeIter() : treeTool( 0 ), childVect( 2 ) {}
 
 ErrorCode initialize( AdaptiveKDTree* tool,
 EntityHandle root,
 const double box_min[3],
 const double box_max[3],
 Direction direction );
 
 AdaptiveKDTree* tool() const
 {
 return treeTool;
 }
 
 //! Get handle for current leaf
 EntityHandle handle() const
 {
 return mStack.back().entity;
 }
 
 //! Get min corner of axis-aligned box for current leaf
 const double* box_min() const
 {
 return mBox[BMIN].array();
 }
 
 //! Get max corner of axis-aligned box for current leaf
 const double* box_max() const
 {
 return mBox[BMAX].array();
 }
 
 double volume() const
 {
 return ( mBox[BMAX][0] - mBox[BMIN][0] ) * ( mBox[BMAX][1] - mBox[BMIN][1] ) *
 ( mBox[BMAX][2] - mBox[BMIN][2] );
 }
 
 //! test if a plane intersects the leaf box
 bool intersects( const AdaptiveKDTree::Plane& plane ) const
 {
 return mBox[BMIN][plane.norm] <= plane.coord && mBox[BMAX][plane.norm] >= plane.coord;
 }
 
 //! Get depth in tree. root is at depth of 1.
 unsigned depth() const
 {
 return mStack.size();
 }
 
 //! Advance the iterator either left or right in the tree
 //! Note: stepping past the end of the tree will invalidate
 //! the iterator. It will *not* be work step the
 //! other direction.
 ErrorCode step( Direction direction );
 
 //! Advance to next leaf
 //! Returns MB_ENTITY_NOT_FOUND if at end.
 //! Note: steping past the end of the tree will invalidate
 //! the iterator. Calling back() will not work.
 ErrorCode step()
 {
 return step( RIGHT );
 }
 
 //! Move back to previous leaf
 //! Returns MB_ENTITY_NOT_FOUND if at beginning.
 //! Note: steping past the start of the tree will invalidate
 //! the iterator. Calling step() will not work.
 ErrorCode back()
 {
 return step( LEFT );
 }
 
 //! Return the side of the box bounding this tree node
 //! that is shared with the immediately adjacent sibling
 //! (the tree node that shares a common parent node with
 //! this node in the binary tree.)
 //!
 //!\param axis_out The principal axis orthogonal to the side of the box
 //!\param neg_out true if the side of the box is toward the decreasing
 //! direction of the principal axis indicated by axis_out,
 //! false if it is toward the increasing direction.
 //!\return MB_ENTITY_NOT FOUND if root node.
 //! MB_FAILURE if internal error.
 //! MB_SUCCESS otherwise.
 ErrorCode sibling_side( AdaptiveKDTree::Axis& axis_out, bool& neg_out ) const;
 
 //! Get adjacent leaf nodes on side indicated by norm and neg.
 //!
 //! E.g. if norm == X and neg == true, then get neighbor(s)
 //! adjacent to the side of the box contained in the plane
 //! with normal to the X axis and with the x coordinate equal
 //! to the minimum x of the bounding box.
 //!
 //! E.g. if norm == Y and neg == false, then get neighbor(s)
 //! adjacent to the side of the box with y = maximum y of bounding box.
 //!
 //!\param norm Normal vector for box side (X, Y, or Z)
 //!\param neg Which of two planes with norm (true->smaller coord,
 //! false->larger coord)
 //!\param results List to which to append results. This function does
 //! *not* clear existing values in list.
 //!\param epsilon Tolerance on overlap. A positive value E will
 //! result in nodes that are separated by as much as E
 //! to be considered touching. A negative value -E will
 //! cause leaves that do not overlap by at least E to be
 //! considered non-overlapping. Amongst other things,
 //! this value can be used to control whether or not
 //! leaves adjacent at only their edges or corners are
 //! returned.
 ErrorCode get_neighbors( AdaptiveKDTree::Axis norm,
 bool neg,
 std::vector< AdaptiveKDTreeIter >& results,
 double epsilon = 0.0 ) const;
 
 //! Get split plane that separates this node from its immediate sibling.
 ErrorCode get_parent_split_plane( AdaptiveKDTree::Plane& plane ) const;
 
 //! Return true if thos node and the passed node share the
 //! same immediate parent.
 bool is_sibling( const AdaptiveKDTreeIter& other_leaf ) const;
 
 //! Return true if thos node and the passed node share the
 //! same immediate parent.
 bool is_sibling( EntityHandle other_leaf ) const;
 
 //! Returns true if calling step() will advance to the
 //! immediate sibling of the current node. Returns false
 //! if current node is root or back() will move to the
 //! immediate sibling.
 bool sibling_is_forward() const;
 
 //! Find range of overlap between ray and leaf.
 //!
 //!\param ray_point Coordinates of start point of ray
 //!\param ray_vect Directionion vector for ray such that
 //! the ray is defined by r(t) = ray_point + t * ray_vect
 //! for t > 0.
 //!\param t_enter Output: if return value is true, this value
 //! is the parameter location along the ray at which
 //! the ray entered the leaf. If return value is false,
 //! then this value is undefined.
 //!\param t_exit Output: if return value is true, this value
 //! is the parameter location along the ray at which
 //! the ray exited the leaf. If return value is false,
 //! then this value is undefined.
 //!\return true if ray intersects leaf, false otherwise.
 bool intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter, double& t_exit ) const;
};
 
inline ErrorCode AdaptiveKDTree::reset_tree()
{
 return delete_tree_sets();
}
 
} // namespace moab
 
#endif