Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
moab::BVHTree Class Reference

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< HandleDataHandleDataVec

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< TreeNodemyTree
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;
}

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:
entitiesEntities with which to build the tree
tree_rootRoot set for tree (see function description)
optsOptions 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_pointPoint to be located in tree
distanceDistance within which to query
result_listLeaves within distance or containing point
iter_tolTolerance for convergence of point search
inside_tolTolerance for inside element calculation
result_distsIf non-NULL, will contain distsances to leaves
result_paramsIf non-NULL, will contain parameters of the point in the ents in leaves_out
tree_rootStart 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:
boxThe box for this tree
tree_nodeIf 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() );
}

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

Parameters:
optsOptions

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:
    pointPoint to be located in tree
    leaf_outLeaf containing point
    iter_tolTolerance for convergence of point search
    inside_tolTolerance for inside element calculation
    multiple_leavesSome 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_nodeStart 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;
}

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:
rootIf 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;
}

Member Data Documentation

Definition at line 323 of file BVHTree.hpp.

Definition at line 325 of file BVHTree.hpp.

Referenced by establish_buckets(), find_split(), order_elements(), and parse_options().

const char * moab::BVHTree::treeName = "BVHTree" [static, private]

Definition at line 327 of file BVHTree.hpp.

Referenced by BVHTree().

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines