BSPTree.hpp

Go to the documentation of this file.

/*
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
/**\file BSPTree.hpp
 *\author Jason Kraftcheck ([email protected])
 *\date 2008-05-12
 */
 
#ifndef MOAB_BSP_TREE_HPP
#define MOAB_BSP_TREE_HPP
 
#include "moab/Types.hpp"
#include "moab/Interface.hpp"
 
#include <cmath>
#include <string>
#include <vector>
 
namespace moab
{
 
class BSPTreeBoxIter;
class BSPTreeIter;
class Range;
class BSPTreePoly;
 
/** \class BSPTree
 * \brief BSP tree, for sorting and searching entities spatially
 */
class BSPTree
{
 private:
 Interface* mbInstance;
 Tag planeTag, rootTag;
 unsigned meshSetFlags;
 bool cleanUpTrees;
 std::vector< EntityHandle > createdTrees;
 
 ErrorCode init_tags( const char* tagname = 0 );
 
 public:
 static double epsilon()
 {
 return 1e-6;
 }
 
 BSPTree( Interface* iface, const char* tagname = 0, unsigned meshset_creation_flags = MESHSET_SET );
 
 BSPTree( Interface* iface,
 bool destroy_created_trees,
 const char* tagname = 0,
 unsigned meshset_creation_flags = MESHSET_SET );
 
 ~BSPTree();
 
 //! Enumerate split plane directions
 enum Axis
 {
 X = 0,
 Y = 1,
 Z = 2
 };
 
 /**\brief struct to store a plane
 *
 * If plane is defined as Ax+By+Cz+D=0, then
 * norm={A,B,C} and coeff=-D.
 */
 struct Plane
 {
 Plane() : coeff( 0.0 ) {}
 Plane( const double n[3], double d ) : coeff( d )
 {
 norm[0] = n[0];
 norm[1] = n[1];
 norm[2] = n[2];
 }
 /** a x + b y + c z + d = 0 */
 Plane( double a, double b, double c, double d ) : coeff( d )
 {
 norm[0] = a;
 norm[1] = b;
 norm[2] = c;
 }
 /** Create Y = 1 plane by doing Plane( Y, 1.0 ); */
 Plane( Axis normal, double point_on_axis ) : coeff( -point_on_axis )
 {
 norm[0] = norm[1] = norm[2] = 0;
 norm[normal] = 1.0;
 }
 
 double norm[3]; //!< Unit normal of plane
 double coeff; //!< norm[0]*x + norm[1]*y + norm[2]*z + coeff = 0
 
 /** Signed distance from point to plane:
 * absolute value is distance from point to plane
 * positive if 'above' plane, negative if 'below'
 * Note: assumes unit-length normal.
 */
 double signed_distance( const double point[3] ) const
 {
 return point[0] * norm[0] + point[1] * norm[1] + point[2] * norm[2] + coeff;
 }
 
 /** return true if point is below the plane */
 bool below( const double point[3] ) const
 {
 return signed_distance( point ) <= 0.0;
 }
 
 /** return true if point is above the plane */
 bool above( const double point[3] ) const
 {
 return signed_distance( point ) >= 0.0;
 }
 
 double distance( const double point[3] ) const
 {
 return fabs( signed_distance( point ) );
 }
 
 /** reverse plane normal */
 void flip()
 {
 norm[0] = -norm[0];
 norm[1] = -norm[1];
 norm[2] = -norm[2];
 coeff = -coeff;
 }
 
 void set( const double normal[3], const double point[3] )
 {
 const double dot = normal[0] * point[0] + normal[1] * point[1] + normal[2] * point[2];
 *this = Plane( normal, -dot );
 }
 
 void set( const double pt1[3], const double pt2[3], const double pt3[3] );
 
 void set( double i, double j, double k, double cff )
 {
 *this = Plane( i, j, k, cff );
 }
 
 /** Create Y = 1 plane by doing set( Y, 1.0 ); */
 void set( Axis normal, double point_on_axis )
 {
 coeff = -point_on_axis;
 norm[0] = norm[1] = norm[2] = 0;
 norm[normal] = 1.0;
 }
 };
 
 //! Get split plane for tree node
 ErrorCode get_split_plane( EntityHandle node, Plane& plane )
 {
 return moab()->tag_get_data( planeTag, &node, 1, &plane );
 }
 
 //! Set split plane for tree node
 ErrorCode set_split_plane( EntityHandle node, const Plane& plane );
 
 //! Get bounding box for entire tree
 ErrorCode get_tree_box( EntityHandle root_node, double corner_coords[8][3] );
 
 //! Get bounding box for entire tree
 ErrorCode get_tree_box( EntityHandle root_node, double corner_coords[24] );
 
 //! Set bounding box for entire tree
 ErrorCode set_tree_box( EntityHandle root_node, const double box_min[3], const double box_max[3] );
 ErrorCode set_tree_box( EntityHandle root_node, const double corner_coords[8][3] );
 
 //! Create tree root node
 ErrorCode create_tree( const double box_min[3], const double box_max[3], EntityHandle& root_handle );
 ErrorCode create_tree( const double corner_coords[8][3], EntityHandle& root_handle );
 
 //! Create tree root node
 ErrorCode create_tree( EntityHandle& root_handle );
 
 //! Find all tree roots
 ErrorCode find_all_trees( Range& results );
 
 //! Destroy a tree
 ErrorCode delete_tree( EntityHandle root_handle );
 
 Interface* moab()
 {
 return mbInstance;
 }
 
 //! Get iterator for tree
 ErrorCode get_tree_iterator( EntityHandle tree_root, BSPTreeIter& result );
 
 //! Get iterator at right-most ('last') leaf.
 ErrorCode get_tree_end_iterator( EntityHandle tree_root, BSPTreeIter& result );
 
 //! Split leaf of tree
 //! Updates iterator location to point to first new leaf node.
 ErrorCode split_leaf( BSPTreeIter& leaf, Plane plane );
 
 //! Split leaf of tree
 //! Updates iterator location to point to first new leaf node.
 ErrorCode split_leaf( BSPTreeIter& 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( BSPTreeIter& 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( BSPTreeIter& 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( BSPTreeIter& iter );
 
 //! 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.
 ErrorCode leaf_containing_point( EntityHandle tree_root, const double point[3], EntityHandle& leaf_out );
 
 //! Get iterator at leaf containing input position.
 //!
 //! Returns MB_ENTITY_NOT_FOUND if point is not within
 //! bounding box of tree.
 ErrorCode leaf_containing_point( EntityHandle tree_root, const double xyz[3], BSPTreeIter& result );
};
 
/** \class BSPTreeIter
 * \brief Iterate over leaves of a BSPTree
 */
class BSPTreeIter
{
 public:
 enum Direction
 {
 LEFT = 0,
 RIGHT = 1
 };
 
 private:
 friend class BSPTree;
 
 BSPTree* treeTool;
 
 protected:
 std::vector< EntityHandle > mStack;
 mutable std::vector< EntityHandle > childVect;
 
 virtual ErrorCode step_to_first_leaf( Direction direction );
 
 virtual ErrorCode up();
 virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction );
 
 virtual ErrorCode initialize( BSPTree* tool, EntityHandle root, const double* point = 0 );
 
 public:
 BSPTreeIter() : treeTool( 0 ), childVect( 2 ) {}
 virtual ~BSPTreeIter() {}
 
 BSPTree* tool() const
 {
 return treeTool;
 }
 
 //! Get handle for current leaf
 EntityHandle handle() const
 {
 return mStack.back();
 }
 
 //! 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.
 virtual ErrorCode step( Direction direction );
 
 //! Advance to next leaf
 //! Returns MB_ENTITY_NOT_FOUND if at end.
 //! Note: stepping 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: stepping past the start of the tree will invalidate
 //! the iterator. Calling step() will not work.
 ErrorCode back()
 {
 return step( LEFT );
 }
 
 //! Get split plane that separates this node from its immediate sibling.
 ErrorCode get_parent_split_plane( BSPTree::Plane& plane ) const;
 
 //! Get volume of leaf polyhedron
 virtual double volume() 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.
 virtual bool intersect_ray( const double ray_point[3],
 const double ray_vect[3],
 double& t_enter,
 double& t_exit ) const;
 
 //! Return true if thos node and the passed node share the
 //! same immediate parent.
 bool is_sibling( const BSPTreeIter& 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;
 
 //! Calculate the convex polyhedron bounding this leaf.
 virtual ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
};
 
/** \class BSPTreeBoxIter
 * \brief Iterate over leaves of a BSPTree
 */
class BSPTreeBoxIter : public BSPTreeIter
{
 private:
 double leafCoords[8][3];
 struct Corners
 {
 double coords[4][3];
 };
 std::vector< Corners > stackData;
 
 protected:
 virtual ErrorCode step_to_first_leaf( Direction direction );
 
 virtual ErrorCode up();
 virtual ErrorCode down( const BSPTree::Plane& plane, Direction direction );
 
 virtual ErrorCode initialize( BSPTree* tool, EntityHandle root, const double* point = 0 );
 
 public:
 BSPTreeBoxIter() {}
 virtual ~BSPTreeBoxIter() {}
 
 //! Faces of a hex : corner bitmap
 enum SideBits
 {
 B0154 = 0x33, //!< Face defined by corners {0,1,5,4}: 1st Exodus side
 B1265 = 0x66, //!< Face defined by corners {1,2,6,5}: 2nd Exodus side
 B2376 = 0xCC, //!< Face defined by corners {2,3,7,6}: 3rd Exodus side
 B3047 = 0x99, //!< Face defined by corners {3,0,4,7}: 4th Exodus side
 B3210 = 0x0F, //!< Face defined by corners {3,2,1,0}: 5th Exodus side
 B4567 = 0xF0 //!< Face defined by corners {4,5,6,7}: 6th Exodus side
 };
 
 static SideBits side_above_plane( const double hex_coords[8][3], const BSPTree::Plane& plane );
 
 static SideBits side_on_plane( const double hex_coords[8][3], const BSPTree::Plane& plane );
 
 static SideBits opposite_face( const SideBits& bits )
 {
 return (SideBits)( ( ~bits ) & 0xFF );
 }
 
 static ErrorCode face_corners( const SideBits face, const double hex_corners[8][3], double face_corners_out[4][3] );
 
 //! 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* work to subsequently
 //! step the other direction.
 virtual ErrorCode step( Direction direction );
 
 //! Advance to next leaf
 //! Returns MB_ENTITY_NOT_FOUND if at end.
 //! Note: stepping past the end of the tree will invalidate
 //! the iterator. Calling back() will not work.
 ErrorCode step()
 {
 return BSPTreeIter::step();
 }
 
 //! Move back to previous leaf
 //! Returns MB_ENTITY_NOT_FOUND if at beginning.
 //! Note: stepping past the start of the tree will invalidate
 //! the iterator. Calling step() will not work.
 ErrorCode back()
 {
 return BSPTreeIter::back();
 }
 
 //! Get coordinates of box corners, in Exodus II hexahedral ordering.
 ErrorCode get_box_corners( double coords[8][3] ) const;
 
 //! Get volume of leaf box
 double volume() const;
 
 //! test if a plane intersects the leaf box
 enum XSect
 {
 MISS = 0,
 SPLIT = 1,
 NONHEX = -1
 };
 XSect splits( const BSPTree::Plane& plane ) const;
 
 //! test if a plane intersects the leaf box
 bool intersects( const BSPTree::Plane& plane ) 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;
 
 //! 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.)
 //!
 //!\return MB_ENTITY_NOT FOUND if root node.
 //! MB_FAILURE if internal error.
 //! MB_SUCCESS otherwise.
 ErrorCode sibling_side( SideBits& side_out ) const;
 
 //! Get adjacent leaf nodes on indicated side
 //!
 //!\param side Face of box for which to retrieve neighbors
 //!\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( SideBits side, std::vector< BSPTreeBoxIter >& results, double epsilon = 0.0 ) const;
 
 //! Calculate the convex polyhedron bounding this leaf.
 ErrorCode calculate_polyhedron( BSPTreePoly& polyhedron_out ) const;
};
 
} // namespace moab
 
#endif