Bounding Volume Hierarchy (sorta like a tree), for sorting and searching entities spatially. More...

#include <BVHTree.hpp>

Inheritance diagram for moab::BVHTree:

Collaboration diagram for moab::BVHTree:

Classes
class	Bucket
class	HandleData
class	HandleData_comparator
class	Node
class	Split_comparator
class	SplitData
class	TreeNode
Public Member Functions
	BVHTree (Interface *impl)
	~BVHTree ()
virtual ErrorCode	reset_tree ()
	Destroy the tree maintained by this object, optionally checking we have the right root.
virtual ErrorCode	parse_options (FileOptions &opts)
	Parse options for this tree, including common options for all trees.
virtual ErrorCode	get_bounding_box (BoundBox &box, EntityHandle *tree_node=NULL) const
	Get bounding box for tree below tree_node, or entire tree If no tree has been built yet, returns +/- DBL_MAX for all dimensions. Note for some tree types, boxes are not available for non-root nodes, and this function will return failure if non-root is passed in.
virtual ErrorCode	build_tree (const Range &entities, EntityHandle tree_root_set=NULL, FileOptions options=NULL)
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)
	Get leaf containing input position.
virtual ErrorCode	distance_search (const double from_point[3], const double distance, std::vector< EntityHandle > &result_list, const double iter_tol=1.0e-10, const double inside_tol=1.0e-6, std::vector< double > result_dists=NULL, std::vector< CartVect > result_params=NULL, EntityHandle *tree_root=NULL)
	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.
virtual ErrorCode	print ()
	print various things about this tree
Private Types
typedef std::vector< HandleData >	HandleDataVec
Private Member Functions
	BVHTree (const BVHTree &s)
void	establish_buckets (HandleDataVec::const_iterator begin, HandleDataVec::const_iterator end, const BoundBox &interval, std::vector< std::vector< Bucket > > &buckets) const
unsigned int	set_interval (BoundBox &interval, std::vector< Bucket >::const_iterator begin, std::vector< Bucket >::const_iterator end) const
void	initialize_splits (std::vector< std::vector< SplitData > > &splits, const std::vector< std::vector< Bucket > > &buckets, const SplitData &data) const
void	order_elements (HandleDataVec::iterator &begin, HandleDataVec::iterator &end, SplitData &data) const
void	median_order (HandleDataVec::iterator &begin, HandleDataVec::iterator &end, SplitData &data) const
void	choose_best_split (const std::vector< std::vector< SplitData > > &splits, SplitData &data) const
void	find_split (HandleDataVec::iterator &begin, HandleDataVec::iterator &end, SplitData &data) const
ErrorCode	find_point (const std::vector< double > &point, const unsigned int &index, const double iter_tol, const double inside_tol, std::pair< EntityHandle, CartVect > &result)
EntityHandle	bruteforce_find (const double *point, const double iter_tol=1.0e-10, const double inside_tol=1.0e-6)
int	local_build_tree (std::vector< Node > &tree_nodes, HandleDataVec::iterator begin, HandleDataVec::iterator end, const int index, const BoundBox &box, const int depth=0)
ErrorCode	construct_element_vec (std::vector< HandleData > &handle_data_vec, const Range &elements, BoundBox &bounding_box)
ErrorCode	convert_tree (std::vector< Node > &tree_nodes)
ErrorCode	print_nodes (std::vector< Node > &nodes)
Private Attributes
Range	entityHandles
std::vector< TreeNode >	myTree
int	splitsPerDir
EntityHandle	startSetHandle
Static Private Attributes
static MOAB_EXPORT const char *	treeName = "BVHTree"

Detailed Description

Bounding Volume Hierarchy (sorta like a tree), for sorting and searching entities spatially.

Definition at line 31 of file BVHTree.hpp.

Member Typedef Documentation

typedef std::vector< HandleData > moab::BVHTree::HandleDataVec [private]

Definition at line 139 of file BVHTree.hpp.

Constructor & Destructor Documentation

moab::BVHTree::BVHTree ( Interface * impl ) [inline]

Definition at line 357 of file BVHTree.hpp.

References moab::Tree::boxTagName, and treeName.

                                         : Tree( impl ), splitsPerDir( 3 ), startSetHandle( 0 )
{
    boxTagName = treeName;
}

moab::BVHTree::~BVHTree ( ) [inline]

Definition at line 36 of file BVHTree.hpp.

References reset_tree().

    {
        reset_tree();
    }

moab::BVHTree::BVHTree ( const BVHTree & s ) [private]

Member Function Documentation

EntityHandle moab::BVHTree::bruteforce_find	(	const double *	point,
		const double	iter_tol = `1.0e-10`,
		const double	inside_tol = `1.0e-6`
	)		`[private]`

Definition at line 547 of file BVHTree.cpp.

References moab::CartVect::array(), dim, ErrorCode, moab::ElemEvaluator::find_containing_entity(), moab::TreeStats::leavesVisited, MB_SUCCESS, moab::Tree::myEval, myTree, moab::TreeStats::numTraversals, startSetHandle, moab::TreeStats::traversalLeafObjectTests, and moab::Tree::treeStats.

{
    treeStats.numTraversals++;
    CartVect params;
    for( unsigned int i = 0; i < myTree.size(); i++ )
    {
        if( myTree[i].dim != 3 || !myTree[i].box.contains_point( point, iter_tol ) ) continue;
        if( myEval )
        {
            EntityHandle entity = 0;
            treeStats.leavesVisited++;
            ErrorCode rval = myEval->find_containing_entity( startSetHandle + i, point, iter_tol, inside_tol, entity,
                                                             params.array(), &treeStats.traversalLeafObjectTests );
            if( entity )
                return entity;
            else if( MB_SUCCESS != rval )
                return 0;
        }
        else
            return startSetHandle + i;
    }
    return 0;
}

ErrorCode moab::BVHTree::build_tree	(	const Range &	entities,
		EntityHandle *	tree_root_set = `NULL`,
		FileOptions *	options = `NULL`
	)		`[virtual]`

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 CANDIDATE_PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below) for possible values; default = 1 (SUBDIVISION_SNAP)

Parameters:

entities	Entities with which to build the tree
tree_root	Root set for tree (see function description)
opts	Options for tree (see function description)

Returns:: Error is returned only on build failure

Implements moab::Tree.

Definition at line 11 of file BVHTree.cpp.

References moab::Interface::add_entities(), moab::FileOptions::all_seen(), moab::CartVect::array(), moab::Range::begin(), moab::BoundBox::bMax, moab::BoundBox::bMin, moab::Tree::boundBox, moab::TreeStats::compute_stats(), construct_element_vec(), convert_tree(), moab::Tree::create_root(), moab::Range::end(), ErrorCode, moab::TreeStats::initTime, moab::BoundBox::intersects_box(), local_build_tree(), moab::Tree::maxPerLeaf, MB_SUCCESS, moab::Tree::mbImpl, moab::Tree::moab(), parse_options(), moab::TreeStats::reset(), moab::Range::size(), startSetHandle, moab::CpuTimer::time_elapsed(), moab::Tree::treeDepth, moab::Tree::treeStats, and moab::BoundBox::update().

{
    ErrorCode rval;
    CpuTimer cp;

    if( options )
    {
        rval = parse_options( *options );
        if( MB_SUCCESS != rval ) return rval;

        if( !options->all_seen() ) return MB_FAILURE;
    }

    // calculate bounding box of elements
    BoundBox box;
    rval = box.update( *moab(), entities );
    if( MB_SUCCESS != rval ) return rval;

    // create tree root
    EntityHandle tmp_root;
    if( !tree_root_set ) tree_root_set = &tmp_root;
    rval = create_root( box.bMin.array(), box.bMax.array(), *tree_root_set );
    if( MB_SUCCESS != rval ) return rval;
    rval = mbImpl->add_entities( *tree_root_set, entities );
    if( MB_SUCCESS != rval ) return rval;

    // a fully balanced tree will have 2*_entities.size()
    // which is one doubling away..
    std::vector< Node > tree_nodes;
    tree_nodes.reserve( entities.size() / maxPerLeaf );
    std::vector< HandleData > handle_data_vec;
    rval = construct_element_vec( handle_data_vec, entities, boundBox );
    if( MB_SUCCESS != rval ) return rval;

#ifndef NDEBUG
    for( std::vector< HandleData >::const_iterator i = handle_data_vec.begin(); i != handle_data_vec.end(); ++i )
    {
        if( !boundBox.intersects_box( i->myBox, 0 ) )
        {
            std::cerr << "BB:" << boundBox << "EB:" << i->myBox << std::endl;
            return MB_FAILURE;
        }
    }
#endif
    // We only build nonempty trees
    if( !handle_data_vec.empty() )
    {
        // initially all bits are set
        tree_nodes.push_back( Node() );
        const int depth = local_build_tree( tree_nodes, handle_data_vec.begin(), handle_data_vec.end(), 0, boundBox );
#ifndef NDEBUG
        std::set< EntityHandle > entity_handles;
        for( std::vector< Node >::iterator n = tree_nodes.begin(); n != tree_nodes.end(); ++n )
        {
            for( HandleDataVec::const_iterator j = n->entities.begin(); j != n->entities.end(); ++j )
            {
                entity_handles.insert( j->myHandle );
            }
        }
        if( entity_handles.size() != entities.size() )
        {
            std::cout << "Entity Handle Size Mismatch!" << std::endl;
        }
        for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
        {
            if( entity_handles.find( *i ) == entity_handles.end() )
                std::cout << "Tree is missing an entity! " << std::endl;
        }
#endif
        treeDepth = std::max( depth, treeDepth );
    }

    // convert vector of Node's to entity sets and vector of TreeNode's
    rval = convert_tree( tree_nodes );
    if( MB_SUCCESS != rval ) return rval;

    treeStats.reset();
    rval               = treeStats.compute_stats( mbImpl, startSetHandle );
    treeStats.initTime = cp.time_elapsed();

    return rval;
}

void moab::BVHTree::choose_best_split	(	const std::vector< std::vector< SplitData > > &	splits,
		SplitData &	data
	)		const `[inline, private]`

Definition at line 398 of file BVHTree.hpp.

Referenced by find_split().

{
    std::vector< SplitData > best_splits;
    for( std::vector< std::vector< SplitData > >::const_iterator i = splits.begin(); i != splits.end(); ++i )
    {
        std::vector< SplitData >::const_iterator j =
            std::min_element( ( *i ).begin(), ( *i ).end(), Split_comparator() );
        best_splits.push_back( *j );
    }
    data = *std::min_element( best_splits.begin(), best_splits.end(), Split_comparator() );
}

ErrorCode moab::BVHTree::construct_element_vec	(	std::vector< HandleData > &	handle_data_vec,
		const Range &	elements,
		BoundBox &	bounding_box
	)		`[inline, private]`

Definition at line 410 of file BVHTree.hpp.

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::CN::MAX_NODES_PER_ELEMENT, MB_SUCCESS, moab::Tree::mbImpl, and moab::BoundBox::update().

Referenced by build_tree().

{
    std::vector< double > coordinate( 3 * CN::MAX_NODES_PER_ELEMENT );
    int num_conn;
    ErrorCode rval = MB_SUCCESS;
    std::vector< EntityHandle > storage;
    for( Range::iterator i = elements.begin(); i != elements.end(); ++i )
    {
        // TODO: not generic enough. Why dim != 3
        // Commence un-necessary deep copying.
        const EntityHandle* connect;
        rval = mbImpl->get_connectivity( *i, connect, num_conn, false, &storage );
        if( MB_SUCCESS != rval ) return rval;
        rval = mbImpl->get_coords( connect, num_conn, &coordinate[0] );
        if( MB_SUCCESS != rval ) return rval;
        BoundBox box;
        for( int j = 0; j < num_conn; j++ )
            box.update( &coordinate[3 * j] );
        if( i == elements.begin() )
            bounding_box = box;
        else
            bounding_box.update( box );
        handle_data_vec.push_back( HandleData( *i, box, 0.0 ) );
    }

    return rval;
}

ErrorCode moab::BVHTree::convert_tree ( std::vector< Node > & tree_nodes ) [private]

Definition at line 94 of file BVHTree.cpp.

References moab::Interface::add_child_meshset(), moab::Interface::add_entities(), moab::ReadUtilIface::create_entity_sets(), ErrorCode, moab::Range::insert(), MB_SUCCESS, moab::Tree::mbImpl, moab::Tree::meshsetFlags, myTree, moab::Interface::query_interface(), moab::Interface::release_interface(), and startSetHandle.

Referenced by build_tree().

{
    // first construct the proper number of entity sets
    ReadUtilIface* read_util;
    ErrorCode rval = mbImpl->query_interface( read_util );
    if( MB_SUCCESS != rval ) return rval;

    {  // isolate potentially-large std::vector so it gets deleted earlier
        std::vector< unsigned int > tmp_flags( tree_nodes.size(), meshsetFlags );
        rval = read_util->create_entity_sets( tree_nodes.size(), &tmp_flags[0], 0, startSetHandle );
        if( MB_SUCCESS != rval ) return rval;
        rval = mbImpl->release_interface( read_util );
        if( MB_SUCCESS != rval ) return rval;
    }

    // populate the sets and the TreeNode vector
    EntityHandle set_handle = startSetHandle;
    std::vector< Node >::iterator it;
    myTree.reserve( tree_nodes.size() );
    for( it = tree_nodes.begin(); it != tree_nodes.end(); ++it, set_handle++ )
    {
        if( it != tree_nodes.begin() && !it->entities.empty() )
        {
            Range range;
            for( HandleDataVec::iterator hit = it->entities.begin(); hit != it->entities.end(); ++hit )
                range.insert( hit->myHandle );
            rval = mbImpl->add_entities( set_handle, range );
            if( MB_SUCCESS != rval ) return rval;
        }
        myTree.push_back( TreeNode( it->dim, it->child, it->Lmax, it->Rmin, it->box ) );

        if( it->dim != 3 )
        {
            rval = mbImpl->add_child_meshset( set_handle, startSetHandle + it->child );
            if( MB_SUCCESS != rval ) return rval;
            rval = mbImpl->add_child_meshset( set_handle, startSetHandle + it->child + 1 );
            if( MB_SUCCESS != rval ) return rval;
        }
    }

    return MB_SUCCESS;
}

ErrorCode moab::BVHTree::distance_search	(	const double	from_point[3],
		const double	distance,
		std::vector< EntityHandle > &	result_list,
		const double	iter_tol = `1.0e-10`,
		const double	inside_tol = `1.0e-6`,
		std::vector< double > *	result_dists = `NULL`,
		std::vector< CartVect > *	result_params = `NULL`,
		EntityHandle *	tree_root = `NULL`
	)		`[virtual]`

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.

Parameters:

from_point	Point to be located in tree
distance	Distance within which to query
result_list	Leaves within distance or containing point
iter_tol	Tolerance for convergence of point search
inside_tol	Tolerance for inside element calculation
result_dists	If non-NULL, will contain distsances to leaves
result_params	If non-NULL, will contain parameters of the point in the ents in leaves_out
tree_root	Start from this tree node (non-NULL) instead of tree root (NULL)

Definition at line 646 of file BVHTree.cpp.

References moab::CartVect::array(), child, dim, moab::BoundBox::distance_squared(), ErrorCode, moab::ElemEvaluator::find_containing_entity(), get_bounding_box(), moab::TreeStats::leavesVisited, moab::Tree::maxDepth, MB_SUCCESS, moab::Tree::myEval, moab::Tree::myRoot, myTree, moab::TreeStats::nodesVisited, moab::TreeStats::numTraversals, startSetHandle, moab::TreeStats::traversalLeafObjectTests, and moab::Tree::treeStats.

{
    // non-NULL root should be in tree
    // convoluted check because the root is different from startSetHandle
    EntityHandle this_set = ( tree_root ? *tree_root : startSetHandle );
    if( this_set != myRoot && ( this_set < startSetHandle || this_set >= startSetHandle + myTree.size() ) )
        return MB_FAILURE;
    else if( this_set == myRoot )
        this_set = startSetHandle;

    treeStats.numTraversals++;

    const double dist_sqr = distance * distance;
    const CartVect from( from_point );
    std::vector< EntityHandle > candidates;  // list of subtrees to traverse
                                             // pre-allocate space for default max tree depth
    candidates.reserve( maxDepth );

    // misc temporary values
    ErrorCode rval;
    BoundBox box;

    candidates.push_back( this_set - startSetHandle );

    while( !candidates.empty() )
    {

        EntityHandle ind = candidates.back();
        this_set         = startSetHandle + ind;
        candidates.pop_back();
        treeStats.nodesVisited++;
        if( myTree[ind].dim == 3 ) treeStats.leavesVisited++;

        // test box of this node
        rval = get_bounding_box( box, &this_set );
        if( MB_SUCCESS != rval ) return rval;
        double d_sqr = box.distance_squared( from_point );

        // if greater than cutoff, continue
        if( d_sqr > dist_sqr ) continue;

        // else if not a leaf, test children & put on list
        else if( myTree[ind].dim != 3 )
        {
            candidates.push_back( myTree[ind].child );
            candidates.push_back( myTree[ind].child + 1 );
            continue;
        }

        if( myEval && result_params )
        {
            EntityHandle ent;
            CartVect params;
            rval = myEval->find_containing_entity( startSetHandle + ind, from_point, iter_tol, inside_tol, ent,
                                                   params.array(), &treeStats.traversalLeafObjectTests );
            if( MB_SUCCESS != rval )
                return rval;
            else if( ent )
            {
                result_list.push_back( ent );
                result_params->push_back( params );
                if( result_dists ) result_dists->push_back( 0.0 );
            }
        }
        else
        {
            // leaf node within distance; return in list
            result_list.push_back( this_set );
            if( result_dists ) result_dists->push_back( sqrt( d_sqr ) );
        }
    }

    return MB_SUCCESS;
}

void moab::BVHTree::establish_buckets	(	HandleDataVec::const_iterator	begin,
		HandleDataVec::const_iterator	end,
		const BoundBox &	interval,
		std::vector< std::vector< Bucket > > &	buckets
	)		const `[private]`

Definition at line 150 of file BVHTree.cpp.

References moab::BVHTree::Bucket::boundingBox, moab::BVHTree::Bucket::bucket_index(), dim, moab::BoundBox::intersects_box(), moab::BVHTree::Bucket::mySize, size, splitsPerDir, and moab::BoundBox::update().

Referenced by find_split().

{
    // put each element into its bucket
    for( HandleDataVec::const_iterator i = begin; i != end; ++i )
    {
        const BoundBox& box = i->myBox;
        for( unsigned int dim = 0; dim < 3; ++dim )
        {
            const unsigned int index = Bucket::bucket_index( splitsPerDir, box, interval, dim );
            assert( index < buckets[dim].size() );
            Bucket& bucket = buckets[dim][index];
            if( bucket.mySize > 0 )
                bucket.boundingBox.update( box );
            else
                bucket.boundingBox = box;
            bucket.mySize++;
        }
    }

#ifndef NDEBUG
    BoundBox elt_union = begin->myBox;
    for( HandleDataVec::const_iterator i = begin; i != end; ++i )
    {
        const BoundBox& box = i->myBox;
        elt_union.update( box );
        for( unsigned int dim = 0; dim < 3; ++dim )
        {
            const unsigned int index = Bucket::bucket_index( splitsPerDir, box, interval, dim );
            Bucket& bucket           = buckets[dim][index];
            if( !bucket.boundingBox.intersects_box( box ) ) std::cerr << "Buckets not covering elements!" << std::endl;
        }
    }
    if( !elt_union.intersects_box( interval ) )
    {
        std::cout << "element union: " << std::endl << elt_union;
        std::cout << "intervals: " << std::endl << interval;
        std::cout << "union of elts does not contain original box!" << std::endl;
    }
    if( !interval.intersects_box( elt_union ) )
    {
        std::cout << "original box does not contain union of elts" << std::endl;
        std::cout << interval << std::endl << elt_union << std::endl;
    }
    for( unsigned int d = 0; d < 3; ++d )
    {
        std::vector< unsigned int > nonempty;
        const std::vector< Bucket >& buckets_ = buckets[d];
        unsigned int j                        = 0;
        for( std::vector< Bucket >::const_iterator i = buckets_.begin(); i != buckets_.end(); ++i, ++j )
        {
            if( i->mySize > 0 )
            {
                nonempty.push_back( j );
            }
        }
        BoundBox test_box = buckets_[nonempty.front()].boundingBox;
        for( unsigned int i = 0; i < nonempty.size(); ++i )
            test_box.update( buckets_[nonempty[i]].boundingBox );

        if( !test_box.intersects_box( interval ) )
            std::cout << "union of buckets in dimension: " << d << "does not contain original box!" << std::endl;
        if( !interval.intersects_box( test_box ) )
        {
            std::cout << "original box does "
                      << "not contain union of buckets"
                      << "in dimension: " << d << std::endl;
            std::cout << interval << std::endl << test_box << std::endl;
        }
    }
#endif
}

ErrorCode moab::BVHTree::find_point	(	const std::vector< double > &	point,
		const unsigned int &	index,
		const double	iter_tol,
		const double	inside_tol,
		std::pair< EntityHandle, CartVect > &	result
	)		`[private]`

Definition at line 466 of file BVHTree.cpp.

References moab::CartVect::array(), moab::Range::begin(), children, moab::BVHTree::TreeNode::dim, moab::Range::end(), entities, ErrorCode, moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_handle(), moab::TreeStats::leavesVisited, moab::BVHTree::TreeNode::Lmax, MB_SUCCESS, moab::Tree::mbImpl, moab::Tree::myEval, myTree, moab::TreeStats::nodesVisited, moab::TreeStats::numTraversals, moab::ElemEvaluator::reverse_eval(), moab::BVHTree::TreeNode::Rmin, moab::ElemEvaluator::set_ent_handle(), startSetHandle, moab::TreeStats::traversalLeafObjectTests, and moab::Tree::treeStats.

{
    if( index == 0 ) treeStats.numTraversals++;
    const TreeNode& node = myTree[index];
    treeStats.nodesVisited++;
    CartVect params;
    int is_inside;
    ErrorCode rval = MB_SUCCESS;
    if( node.dim == 3 )
    {
        treeStats.leavesVisited++;
        Range entities;
        rval = mbImpl->get_entities_by_handle( startSetHandle + index, entities );
        if( MB_SUCCESS != rval ) return rval;

        for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
        {
            treeStats.traversalLeafObjectTests++;
            myEval->set_ent_handle( *i );
            myEval->reverse_eval( &point[0], iter_tol, inside_tol, params.array(), &is_inside );
            if( is_inside )
            {
                result.first  = *i;
                result.second = params;
                return MB_SUCCESS;
            }
        }
        result.first = 0;
        return MB_SUCCESS;
    }
    // the extra tol here considers the case where
    // 0 < Rmin - Lmax < 2tol
    std::vector< EntityHandle > children;
    rval = mbImpl->get_child_meshsets( startSetHandle + index, children );
    if( MB_SUCCESS != rval || children.size() != 2 ) return rval;

    if( ( node.Lmax + iter_tol ) < ( node.Rmin - iter_tol ) )
    {
        if( point[node.dim] <= ( node.Lmax + iter_tol ) )
            return find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
        else if( point[node.dim] >= ( node.Rmin - iter_tol ) )
            return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
        result = std::make_pair( 0, CartVect( &point[0] ) );  // point lies in empty space.
        return MB_SUCCESS;
    }

    // Boxes overlap
    // left of Rmin, you must be on the left
    // we can't be sure about the boundaries since the boxes overlap
    // this was a typo in the paper which caused pain.
    if( point[node.dim] < ( node.Rmin - iter_tol ) )
        return find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
    // if you are on the right Lmax, you must be on the right
    else if( point[node.dim] > ( node.Lmax + iter_tol ) )
        return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );

    /* pg5 of paper
     * However, instead of always traversing either subtree
     * first (e.g. left always before right), we first traverse
     * the subtree whose bounding plane has the larger distance to the
     * sought point. This results in less overall traversal, and the correct
     * cell is identified more quickly.
     */
    // So far all testing confirms that this 'heuristic' is
    // significantly slower.
    // I conjecture this is because it gets improperly
    // branch predicted..
    // bool dir = (point[ node.dim] - node.Rmin) <=
    //              (node.Lmax - point[ node.dim]);
    find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
    if( result.first == 0 )
    {
        return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
    }
    return MB_SUCCESS;
}

void moab::BVHTree::find_split	(	HandleDataVec::iterator &	begin,
		HandleDataVec::iterator &	end,
		SplitData &	data
	)		const `[private]`

Definition at line 320 of file BVHTree.cpp.

References moab::BVHTree::SplitData::boundingBox, choose_best_split(), moab::BVHTree::SplitData::dim, establish_buckets(), initialize_splits(), moab::BoundBox::intersects_box(), left_box, moab::BVHTree::SplitData::leftBox, moab::BVHTree::SplitData::Lmax, median_order(), moab::BVHTree::SplitData::nl, moab::BVHTree::SplitData::nr, order_elements(), right_box, moab::BVHTree::SplitData::rightBox, moab::BVHTree::SplitData::Rmin, splitsPerDir, and moab::BoundBox::update().

Referenced by local_build_tree().

{
    std::vector< std::vector< Bucket > > buckets( 3, std::vector< Bucket >( splitsPerDir + 1 ) );
    std::vector< std::vector< SplitData > > splits( 3, std::vector< SplitData >( splitsPerDir, data ) );

    const BoundBox interval = data.boundingBox;
    establish_buckets( begin, end, interval, buckets );
    initialize_splits( splits, buckets, data );
    choose_best_split( splits, data );
    const bool use_median = ( 0 == data.nl ) || ( data.nr == 0 );
    if( !use_median )
        order_elements( begin, end, data );
    else
        median_order( begin, end, data );

#ifndef NDEBUG
    bool seen_one = false, issue = false;
    bool first_left = true, first_right = true;
    unsigned int count_left = 0, count_right = 0;
    BoundBox left_box, right_box;
    for( HandleDataVec::iterator i = begin; i != end; ++i )
    {
        int order = i->myDim;
        if( order != 0 && order != 1 )
        {
            std::cerr << "Invalid order element !";
            std::cerr << order << std::endl;
            std::exit( -1 );
        }
        if( order == 1 )
        {
            seen_one = 1;
            count_right++;
            if( first_right )
            {
                right_box   = i->myBox;
                first_right = false;
            }
            else
            {
                right_box.update( i->myBox );
            }
            if( !right_box.intersects_box( i->myBox ) )
            {
                if( !issue )
                {
                    std::cerr << "Bounding right box issue!" << std::endl;
                }
                issue = true;
            }
        }
        if( order == 0 )
        {
            count_left++;
            if( first_left )
            {
                left_box   = i->myBox;
                first_left = false;
            }
            else
            {
                left_box.update( i->myBox );
            }
            if( !data.leftBox.intersects_box( i->myBox ) )
            {
                if( !issue )
                {
                    std::cerr << "Bounding left box issue!" << std::endl;
                }
                issue = true;
            }
            if( seen_one )
            {
                std::cerr << "Invalid ordering!" << std::endl;
                std::cout << ( i - 1 )->myDim << order << std::endl;
                exit( -1 );
            }
        }
    }
    if( !left_box.intersects_box( data.leftBox ) ) std::cout << "left elts do not contain left box" << std::endl;
    if( !data.leftBox.intersects_box( left_box ) ) std::cout << "left box does not contain left elts" << std::endl;
    if( !right_box.intersects_box( data.rightBox ) ) std::cout << "right elts do not contain right box" << std::endl;
    if( !data.rightBox.intersects_box( right_box ) ) std::cout << "right box do not contain right elts" << std::endl;

    if( count_left != data.nl || count_right != data.nr )
    {
        std::cerr << "counts are off!" << std::endl;
        std::cerr << "total: " << std::distance( begin, end ) << std::endl;
        std::cerr << "Dim: " << data.dim << std::endl;
        std::cerr << data.Lmax << " , " << data.Rmin << std::endl;
        std::cerr << "Right box: " << std::endl << data.rightBox << "Left box: " << std::endl << data.leftBox;
        std::cerr << "supposed to be: " << data.nl << " " << data.nr << std::endl;
        std::cerr << "accountant says: " << count_left << " " << count_right << std::endl;
        std::exit( -1 );
    }
#endif
}

ErrorCode moab::BVHTree::get_bounding_box	(	BoundBox &	box,
		EntityHandle *	tree_node = `NULL`
	)		const `[virtual]`

Get bounding box for tree below tree_node, or entire tree If no tree has been built yet, returns +/- DBL_MAX for all dimensions. Note for some tree types, boxes are not available for non-root nodes, and this function will return failure if non-root is passed in.

Parameters:

box	The box for this tree
tree_node	If non-NULL, node for which box is requested, tree root if NULL

Returns:: Only returns error on fatal condition

Reimplemented from moab::Tree.

Definition at line 571 of file BVHTree.cpp.

References moab::Tree::boundBox, MB_SUCCESS, myTree, and startSetHandle.

Referenced by distance_search(), and point_search().

{
    if( !tree_node || *tree_node == startSetHandle )
    {
        box = boundBox;
        return MB_SUCCESS;
    }
    else if( ( tree_node && !startSetHandle ) || *tree_node < startSetHandle ||
             *tree_node - startSetHandle > myTree.size() )
        return MB_FAILURE;

    box = myTree[*tree_node - startSetHandle].box;
    return MB_SUCCESS;
}

void moab::BVHTree::initialize_splits	(	std::vector< std::vector< SplitData > > &	splits,
		const std::vector< std::vector< Bucket > > &	buckets,
		const SplitData &	data
	)		const `[private]`

Definition at line 225 of file BVHTree.cpp.

References moab::BoundBox::bMax, moab::BoundBox::bMin, moab::BVHTree::SplitData::boundingBox, moab::BoundBox::intersects_box(), set_interval(), and moab::BoundBox::update().

Referenced by find_split().

{
    for( unsigned int d = 0; d < 3; ++d )
    {
        std::vector< SplitData >::iterator splits_begin  = splits[d].begin();
        std::vector< SplitData >::iterator splits_end    = splits[d].end();
        std::vector< Bucket >::const_iterator left_begin = buckets[d].begin();
        std::vector< Bucket >::const_iterator _end       = buckets[d].end();
        std::vector< Bucket >::const_iterator left_end   = buckets[d].begin() + 1;
        for( std::vector< SplitData >::iterator s = splits_begin; s != splits_end; ++s, ++left_end )
        {
            s->nl = set_interval( s->leftBox, left_begin, left_end );
            if( s->nl == 0 )
            {
                s->leftBox         = data.boundingBox;
                s->leftBox.bMax[d] = s->leftBox.bMin[d];
            }
            s->nr = set_interval( s->rightBox, left_end, _end );
            if( s->nr == 0 )
            {
                s->rightBox         = data.boundingBox;
                s->rightBox.bMin[d] = s->rightBox.bMax[d];
            }
            s->Lmax  = s->leftBox.bMax[d];
            s->Rmin  = s->rightBox.bMin[d];
            s->split = std::distance( splits_begin, s );
            s->dim   = d;
        }
#ifndef NDEBUG
        for( std::vector< SplitData >::iterator s = splits_begin; s != splits_end; ++s )
        {
            BoundBox test_box = s->leftBox;
            test_box.update( s->rightBox );
            if( !data.boundingBox.intersects_box( test_box ) )
            {
                std::cout << "nr: " << s->nr << std::endl;
                std::cout << "Test box: " << std::endl << test_box;
                std::cout << "Left box: " << std::endl << s->leftBox;
                std::cout << "Right box: " << std::endl << s->rightBox;
                std::cout << "Interval: " << std::endl << data.boundingBox;
                std::cout << "Split boxes larger than bb" << std::endl;
            }
            if( !test_box.intersects_box( data.boundingBox ) )
            {
                std::cout << "bb larger than union of split boxes" << std::endl;
            }
        }
#endif
    }
}

int moab::BVHTree::local_build_tree	(	std::vector< Node > &	tree_nodes,
		HandleDataVec::iterator	begin,
		HandleDataVec::iterator	end,
		const int	index,
		const BoundBox &	box,
		const int	depth = `0`
	)		`[private]`

Definition at line 418 of file BVHTree.cpp.

References moab::BVHTree::SplitData::boundingBox, moab::BVHTree::SplitData::dim, entities, find_split(), moab::BoundBox::intersects_box(), moab::BVHTree::SplitData::leftBox, moab::BVHTree::SplitData::Lmax, moab::Tree::maxDepth, moab::Tree::maxPerLeaf, moab::BVHTree::SplitData::nl, moab::BVHTree::SplitData::rightBox, and moab::BVHTree::SplitData::Rmin.

Referenced by build_tree().

{
#ifndef NDEBUG
    for( HandleDataVec::const_iterator i = begin; i != end; ++i )
    {
        if( !box.intersects_box( i->myBox, 0 ) )
        {
            std::cerr << "depth: " << depth << std::endl;
            std::cerr << "BB:" << box << "EB:" << i->myBox << std::endl;
            std::exit( -1 );
        }
    }
#endif

    const unsigned int total_num_elements = std::distance( begin, end );
    tree_nodes[index].box                 = box;

    // logic for splitting conditions
    if( (int)total_num_elements > maxPerLeaf && depth < maxDepth )
    {
        SplitData data;
        data.boundingBox = box;
        find_split( begin, end, data );
        // assign data to node
        tree_nodes[index].Lmax  = data.Lmax;
        tree_nodes[index].Rmin  = data.Rmin;
        tree_nodes[index].dim   = data.dim;
        tree_nodes[index].child = tree_nodes.size();
        // insert left, right children;
        tree_nodes.push_back( Node() );
        tree_nodes.push_back( Node() );
        const int left_depth =
            local_build_tree( tree_nodes, begin, begin + data.nl, tree_nodes[index].child, data.leftBox, depth + 1 );
        const int right_depth =
            local_build_tree( tree_nodes, begin + data.nl, end, tree_nodes[index].child + 1, data.rightBox, depth + 1 );
        return std::max( left_depth, right_depth );
    }

    tree_nodes[index].dim = 3;
    std::copy( begin, end, std::back_inserter( tree_nodes[index].entities ) );
    return depth;
}

void moab::BVHTree::median_order	(	HandleDataVec::iterator &	begin,
		HandleDataVec::iterator &	end,
		SplitData &	data
	)		const `[private]`

Definition at line 278 of file BVHTree.cpp.

References moab::BoundBox::bMax, moab::BoundBox::bMin, moab::BVHTree::SplitData::boundingBox, dim, moab::BVHTree::SplitData::dim, moab::BoundBox::intersects_box(), moab::BVHTree::SplitData::leftBox, moab::BVHTree::SplitData::Lmax, moab::BVHTree::SplitData::nl, moab::BVHTree::SplitData::nr, moab::BVHTree::SplitData::rightBox, moab::BVHTree::SplitData::Rmin, moab::BVHTree::SplitData::split, and moab::BoundBox::update().

Referenced by find_split().

{
    int dim = data.dim;
    for( HandleDataVec::iterator i = begin; i != end; ++i )
    {
        i->myDim = 0.5 * ( i->myBox.bMin[dim], i->myBox.bMax[dim] );
    }
    std::sort( begin, end, BVHTree::HandleData_comparator() );
    const unsigned int total       = std::distance( begin, end );
    HandleDataVec::iterator middle = begin + ( total / 2 );
    double middle_center           = middle->myDim;
    middle_center += ( ++middle )->myDim;
    middle_center *= 0.5;
    data.split = middle_center;
    data.nl    = std::distance( begin, middle ) + 1;
    data.nr    = total - data.nl;
    ++middle;
    data.leftBox  = begin->myBox;
    data.rightBox = middle->myBox;
    for( HandleDataVec::iterator i = begin; i != middle; ++i )
    {
        i->myDim = 0;
        data.leftBox.update( i->myBox );
    }
    for( HandleDataVec::iterator i = middle; i != end; ++i )
    {
        i->myDim = 1;
        data.rightBox.update( i->myBox );
    }
    data.Rmin = data.rightBox.bMin[data.dim];
    data.Lmax = data.leftBox.bMax[data.dim];
#ifndef NDEBUG
    BoundBox test_box( data.rightBox );
    if( !data.boundingBox.intersects_box( test_box ) )
    {
        std::cerr << "MEDIAN: BB Does not contain splits" << std::endl;
        std::cerr << "test_box:         " << test_box << std::endl;
        std::cerr << "data.boundingBox: " << data.boundingBox << std::endl;
    }
#endif
}

void moab::BVHTree::order_elements	(	HandleDataVec::iterator &	begin,
		HandleDataVec::iterator &	end,
		SplitData &	data
	)		const `[inline, private]`

Definition at line 386 of file BVHTree.hpp.

References moab::BVHTree::SplitData::boundingBox, moab::BVHTree::Bucket::bucket_index(), moab::BVHTree::SplitData::dim, moab::BVHTree::SplitData::split, and splitsPerDir.

Referenced by find_split().

{
    for( HandleDataVec::iterator i = begin; i != end; ++i )
    {
        const int index = Bucket::bucket_index( splitsPerDir, i->myBox, data.boundingBox, data.dim );
        i->myDim        = ( index <= data.split ) ? 0 : 1;
    }
    std::sort( begin, end, HandleData_comparator() );
}

ErrorCode moab::BVHTree::parse_options ( FileOptions & opts ) [virtual]

Parse options for this tree, including common options for all trees.

Parameters:

opts Options

Implements moab::Tree.

Definition at line 137 of file BVHTree.cpp.

References ErrorCode, moab::FileOptions::get_int_option(), MB_SUCCESS, moab::Tree::parse_common_options(), and splitsPerDir.

Referenced by build_tree().

{
    ErrorCode rval = parse_common_options( opts );
    if( MB_SUCCESS != rval ) return rval;

    //  SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
    int tmp_int;
    rval = opts.get_int_option( "SPLITS_PER_DIR", tmp_int );
    if( MB_SUCCESS == rval ) splitsPerDir = tmp_int;

    return MB_SUCCESS;
}

ErrorCode moab::BVHTree::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`
	)		`[virtual]`

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.

Parameters:

point	Point to be located in tree
leaf_out	Leaf containing point
iter_tol	Tolerance for convergence of point search
inside_tol	Tolerance for inside element calculation
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
start_node	Start from this tree node (non-NULL) instead of tree root (NULL)

Returns:: Non-success returned only in case of failure; not-found indicated by leaf_out=0

Implements moab::Tree.

Definition at line 586 of file BVHTree.cpp.

References moab::CartVect::array(), child, moab::BoundBox::contains_point(), dim, ErrorCode, moab::ElemEvaluator::find_containing_entity(), get_bounding_box(), moab::TreeStats::leavesVisited, MB_SUCCESS, moab::Tree::myEval, moab::Tree::myRoot, myTree, moab::TreeStats::nodesVisited, moab::TreeStats::numTraversals, startSetHandle, moab::TreeStats::traversalLeafObjectTests, and moab::Tree::treeStats.

{
    treeStats.numTraversals++;

    EntityHandle this_set = ( start_node ? *start_node : startSetHandle );
    // convoluted check because the root is different from startSetHandle
    if( this_set != myRoot && ( this_set < startSetHandle || this_set >= startSetHandle + myTree.size() ) )
        return MB_FAILURE;
    else if( this_set == myRoot )
        this_set = startSetHandle;

    std::vector< EntityHandle > candidates,
        result_list;  // list of subtrees to traverse, and results
    candidates.push_back( this_set - startSetHandle );

    BoundBox box;
    while( !candidates.empty() )
    {
        EntityHandle ind = candidates.back();
        treeStats.nodesVisited++;
        if( myTree[ind].dim == 3 ) treeStats.leavesVisited++;
        this_set = startSetHandle + ind;
        candidates.pop_back();

        // test box of this node
        ErrorCode rval = get_bounding_box( box, &this_set );
        if( MB_SUCCESS != rval ) return rval;
        if( !box.contains_point( point, iter_tol ) ) continue;

        // else if not a leaf, test children & put on list
        else if( myTree[ind].dim != 3 )
        {
            candidates.push_back( myTree[ind].child );
            candidates.push_back( myTree[ind].child + 1 );
            continue;
        }
        else if( myTree[ind].dim == 3 && myEval && params )
        {
            rval = myEval->find_containing_entity( startSetHandle + ind, point, iter_tol, inside_tol, leaf_out,
                                                   params->array(), &treeStats.traversalLeafObjectTests );
            if( leaf_out || MB_SUCCESS != rval ) return rval;
        }
        else
        {
            // leaf node within distance; return in list
            result_list.push_back( this_set );
        }
    }

    if( !result_list.empty() ) leaf_out = result_list[0];
    if( multiple_leaves && result_list.size() > 1 ) *multiple_leaves = true;
    return MB_SUCCESS;
}

ErrorCode moab::BVHTree::print ( ) [virtual]

print various things about this tree

Implements moab::Tree.

Definition at line 740 of file BVHTree.cpp.

References MB_SUCCESS, and myTree.

{
    int i;
    std::vector< TreeNode >::iterator it;
    for( it = myTree.begin(), i = 0; it != myTree.end(); ++it, i++ )
    {
        std::cout << "Node " << i << ": dim = " << it->dim << ", child = " << it->child << ", Lmax/Rmin = " << it->Lmax
                  << "/" << it->Rmin << ", box = " << it->box << std::endl;
    }
    return MB_SUCCESS;
}

ErrorCode moab::BVHTree::print_nodes ( std::vector< Node > & nodes ) [private]

Definition at line 728 of file BVHTree.cpp.

References MB_SUCCESS.

{
    int i;
    std::vector< Node >::iterator it;
    for( it = nodes.begin(), i = 0; it != nodes.end(); ++it, i++ )
    {
        std::cout << "Node " << i << ": dim = " << it->dim << ", child = " << it->child << ", Lmax/Rmin = " << it->Lmax
                  << "/" << it->Rmin << ", box = " << it->box << std::endl;
    }
    return MB_SUCCESS;
}

ErrorCode moab::BVHTree::reset_tree ( ) [inline, virtual]

Destroy the tree maintained by this object, optionally checking we have the right root.

Parameters:

root	If non-NULL, check that this is the root, return failure if not

Implements moab::Tree.

Definition at line 440 of file BVHTree.hpp.

References moab::Tree::delete_tree_sets().

Referenced by ~BVHTree().

{
    return delete_tree_sets();
}

unsigned int moab::BVHTree::set_interval	(	BoundBox &	interval,
		std::vector< Bucket >::const_iterator	begin,
		std::vector< Bucket >::const_iterator	end
	)		const `[inline, private]`

Definition at line 362 of file BVHTree.hpp.

References moab::GeomUtil::first(), and moab::BoundBox::update().

Referenced by initialize_splits().

{
    bool first         = true;
    unsigned int count = 0;
    for( std::vector< Bucket >::const_iterator b = begin; b != end; ++b )
    {
        const BoundBox& box = b->boundingBox;
        count += b->mySize;
        if( b->mySize != 0 )
        {
            if( first )
            {
                interval = box;
                first    = false;
            }
            else
                interval.update( box );
        }
    }
    return count;
}