AdaptiveKDTree.cpp

Go to the documentation of this file.

/**\file AdaptiveKDTree.cpp
 */
 
#include "moab/AdaptiveKDTree.hpp"
#include "moab/Interface.hpp"
#include "moab/GeomUtil.hpp"
#include "moab/Range.hpp"
#include "moab/ElemEvaluator.hpp"
#include "moab/CpuTimer.hpp"
#include "Internals.hpp"
#include "moab/Util.hpp"
#include <cmath>
 
#include <cassert>
#include <algorithm>
#include <limits>
#include <iostream>
#include <cstdio>
 
namespace moab
{
 
const char* AdaptiveKDTree::treeName = "AKDTree";
 
#define MB_AD_KD_TREE_DEFAULT_TAG_NAME
 
// If defined, use single tag for both axis and location of split plane
#define MB_AD_KD_TREE_USE_SINGLE_TAG
 
// No effect if MB_AD_KD_TREE_USE_SINGLE_TAG is not defined.
// If defined, store plane axis as double so tag has consistent
// type (doubles for both location and axis). If not defined,
// store struct Plane as opaque.
#define MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG
 
AdaptiveKDTree::AdaptiveKDTree( Interface* iface )
 : Tree( iface ), planeTag( 0 ), axisTag( 0 ), splitsPerDir( 3 ), planeSet( SUBDIVISION_SNAP ), spherical( false ),
 radius( 1.0 )
{
 boxTagName = treeName;
 
 ErrorCode rval = init();
 if( MB_SUCCESS != rval ) throw rval;
}
 
AdaptiveKDTree::AdaptiveKDTree( Interface* iface,
 const Range& entities,
 EntityHandle* tree_root_set,
 FileOptions* opts )
 : Tree( iface ), planeTag( 0 ), axisTag( 0 ), splitsPerDir( 3 ), planeSet( SUBDIVISION_SNAP ), spherical( false ),
 radius( 1.0 )
{
 boxTagName = treeName;
 
 ErrorCode rval;
 if( opts )
 {
 rval = parse_options( *opts );
 if( MB_SUCCESS != rval ) throw rval;
 }
 
 rval = init();
 if( MB_SUCCESS != rval ) throw rval;
 
 rval = build_tree( entities, tree_root_set, opts );
 if( MB_SUCCESS != rval ) throw rval;
}
 
AdaptiveKDTree::~AdaptiveKDTree()
{
 if( !cleanUp ) return;
 
 if( myRoot )
 {
 reset_tree();
 myRoot = 0;
 }
}
 
ErrorCode AdaptiveKDTree::build_tree( const Range& entities, EntityHandle* tree_root_set, FileOptions* options )
{
 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, spherical, radius );
 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 = moab()->add_entities( *tree_root_set, entities );
 if( MB_SUCCESS != rval ) return rval;
 
 AdaptiveKDTreeIter iter;
 iter.initialize( this, *tree_root_set, box.bMin.array(), box.bMax.array(), AdaptiveKDTreeIter::LEFT );
 
 std::vector< double > tmp_data;
 std::vector< EntityHandle > tmp_data2;
 for( ;; )
 {
 
 int pcount;
 rval = moab()->get_number_entities_by_handle( iter.handle(), pcount );
 if( MB_SUCCESS != rval ) break;
 
 const size_t p_count = pcount;
 Range best_left, best_right, best_both;
 Plane best_plane = { HUGE_VAL, -1 };
 if( (int)p_count > maxPerLeaf && (int)iter.depth() < maxDepth )
 {
 switch( planeSet )
 {
 case AdaptiveKDTree::SUBDIVISION:
 rval = best_subdivision_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
 minWidth );
 break;
 case AdaptiveKDTree::SUBDIVISION_SNAP:
 rval = best_subdivision_snap_plane( splitsPerDir, iter, best_left, best_right, best_both,
 best_plane, tmp_data, minWidth );
 break;
 case AdaptiveKDTree::VERTEX_MEDIAN:
 rval = best_vertex_median_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
 tmp_data, minWidth );
 break;
 case AdaptiveKDTree::VERTEX_SAMPLE:
 rval = best_vertex_sample_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
 tmp_data, tmp_data2, minWidth );
 break;
 default:
 rval = MB_FAILURE;
 }
 
 if( MB_SUCCESS != rval ) return rval;
 }
 
 if( best_plane.norm >= 0 )
 {
 best_left.merge( best_both );
 best_right.merge( best_both );
 rval = split_leaf( iter, best_plane, best_left, best_right );
 if( MB_SUCCESS != rval ) return rval;
 }
 else
 {
 rval = iter.step();
 if( MB_ENTITY_NOT_FOUND == rval )
 {
 rval = treeStats.compute_stats( mbImpl, myRoot );
 treeStats.initTime = cp.time_elapsed();
 return rval; // at end
 }
 else if( MB_SUCCESS != rval )
 break;
 }
 }
 
 reset_tree();
 
 treeStats.reset();
 
 return rval;
}
 
ErrorCode AdaptiveKDTree::parse_options( FileOptions& opts )
{
 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;
 
 // PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
 // for possible values; default = 1 (SUBDIVISION_SNAP)
 rval = opts.get_int_option( "PLANE_SET", tmp_int );
 if( MB_SUCCESS == rval && ( tmp_int < SUBDIVISION || tmp_int > VERTEX_SAMPLE ) )
 return MB_FAILURE;
 else if( MB_ENTITY_NOT_FOUND == rval )
 planeSet = SUBDIVISION;
 else
 planeSet = (CandidatePlaneSet)( tmp_int );
 
 rval = opts.get_toggle_option( "SPHERICAL", false, spherical );
 if( MB_SUCCESS != rval ) spherical = false;
 
 double tmp = 1.0;
 rval = opts.get_real_option( "RADIUS", tmp );
 if( MB_SUCCESS != rval )
 radius = 1.0;
 else
 radius = tmp;
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::make_tag( Interface* iface,
 std::string name,
 TagType storage,
 DataType type,
 int count,
 void* default_val,
 Tag& tag_handle,
 std::vector< Tag >& created_tags )
{
 ErrorCode rval =
 iface->tag_get_handle( name.c_str(), count, type, tag_handle, MB_TAG_CREAT | storage, default_val );
 
 if( MB_SUCCESS == rval )
 {
 if( std::find( created_tags.begin(), created_tags.end(), tag_handle ) == created_tags.end() )
 created_tags.push_back( tag_handle );
 }
 else
 {
 while( !created_tags.empty() )
 {
 iface->tag_delete( created_tags.back() );
 created_tags.pop_back();
 }
 
 planeTag = axisTag = (Tag)-1;
 }
 
 return rval;
}
 
ErrorCode AdaptiveKDTree::init()
{
 std::vector< Tag > ctl;
 
#ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
 // create two tags, one for axis direction and one for axis coordinate
 std::string n1( treeName ), n2( treeName );
 n1 += "_coord";
 n2 += "_norm";
 ErrorCode rval = make_tag( moab(), n1, MB_TAG_DENSE, MB_TYPE_DOUBLE, 1, 0, planeTag, ctl );
 if( MB_SUCCESS != rval ) return rval;
 rval = make_tag( moab(), n2, MB_TAG_DENSE, MB_TYPE_INT, 1, 0, axisTag, ctl );
 if( MB_SUCCESS != rval ) return rval;
 
#elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
 // create tag to hold two doubles, one for location and one for axis
 std::string double_tag_name = std::string( treeName ) + std::string( "_coord_norm" );
 ErrorCode rval = make_tag( moab(), double_tag_name, MB_TAG_DENSE, MB_TYPE_DOUBLE, 2, 0, planeTag, ctl );
 if( MB_SUCCESS != rval ) return rval;
#else
 // create opaque tag to hold struct Plane
 ErrorCode rval = make_tag( moab(), tagname, MB_TAG_DENSE, MB_TYPE_OPAQUE, sizeof( Plane ), 0, planeTag, ctl );
 if( MB_SUCCESS != rval ) return rval;
 
#ifdef MOAB_HAVE_HDF5
 // create a mesh tag holding the HDF5 type for a struct Plane
 Tag type_tag;
 std::string type_tag_name = "__hdf5_tag_type_";
 type_tag_name += boxTagName;
 rval = make_tag( moab(), type_tag_name, MB_TAG_MESH, MB_TYPE_OPAQUE, sizeof( hid_t ), 0, type_tag, ctl );
 if( MB_SUCCESS != rval ) return rval;
 // create HDF5 type object describing struct Plane
 Plane p;
 hid_t handle = H5Tcreate( H5T_COMPOUND, sizeof( Plane ) );
 H5Tinsert( handle, "coord", &( p.coord ) - &p, H5T_NATIVE_DOUBLE );
 H5Tinsert( handle, "norm", &( p.axis ) - &p, H5T_NATIVE_INT );
 EntityHandle root = 0;
 rval = mbImpl->tag_set_data( type_tag, &root, 1, &handle );
 if( MB_SUCCESS != rval ) return rval;
#endif
#endif
 
 return rval;
}
 
ErrorCode AdaptiveKDTree::get_split_plane( EntityHandle entity, Plane& plane )
{
#ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
 ErrorCode r1, r2;
 r1 = moab()->tag_get_data( planeTag, &entity, 1, &plane.coord );
 r2 = moab()->tag_get_data( axisTag, &entity, 1, &plane.norm );
 return MB_SUCCESS == r1 ? r2 : r1;
#elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
 double values[2];
 ErrorCode rval = moab()->tag_get_data( planeTag, &entity, 1, values );
 plane.coord = values[0];
 plane.norm = (int)values[1];
 return rval;
#else
 return moab()->tag_get_data( planeTag, &entity, 1, &plane );
#endif
}
 
ErrorCode AdaptiveKDTree::set_split_plane( EntityHandle entity, const Plane& plane )
{
#ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
 ErrorCode r1, r2;
 r1 = moab()->tag_set_data( planeTag, &entity, 1, &plane.coord );
 r2 = moab()->tag_set_data( axisTag, &entity, 1, &plane.norm );
 return MB_SUCCESS == r1 ? r2 : r1;
#elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
 double values[2] = { plane.coord, static_cast< double >( plane.norm ) };
 return moab()->tag_set_data( planeTag, &entity, 1, values );
#else
 return moab()->tag_set_data( planeTag, &entity, 1, &plane );
#endif
}
 
ErrorCode AdaptiveKDTree::get_tree_iterator( EntityHandle root, AdaptiveKDTreeIter& iter )
{
 double box[6];
 ErrorCode rval = moab()->tag_get_data( boxTag, &root, 1, box );
 if( MB_SUCCESS != rval ) return rval;
 
 return get_sub_tree_iterator( root, box, box + 3, iter );
}
 
ErrorCode AdaptiveKDTree::get_last_iterator( EntityHandle root, AdaptiveKDTreeIter& iter )
{
 double box[6];
 ErrorCode rval = moab()->tag_get_data( boxTag, &root, 1, box );
 if( MB_SUCCESS != rval ) return rval;
 
 return iter.initialize( this, root, box, box + 3, AdaptiveKDTreeIter::RIGHT );
}
 
ErrorCode AdaptiveKDTree::get_sub_tree_iterator( EntityHandle root,
 const double min[3],
 const double max[3],
 AdaptiveKDTreeIter& result )
{
 return result.initialize( this, root, min, max, AdaptiveKDTreeIter::LEFT );
}
 
ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, EntityHandle& left, EntityHandle& right )
{
 ErrorCode rval;
 
 rval = moab()->create_meshset( meshsetFlags, left );
 if( MB_SUCCESS != rval ) return rval;
 
 rval = moab()->create_meshset( meshsetFlags, right );
 if( MB_SUCCESS != rval )
 {
 moab()->delete_entities( &left, 1 );
 return rval;
 }
 
 if( MB_SUCCESS != set_split_plane( leaf.handle(), plane ) ||
 MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), left ) ||
 MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), right ) ||
 MB_SUCCESS != leaf.step_to_first_leaf( AdaptiveKDTreeIter::LEFT ) )
 {
 EntityHandle children[] = { left, right };
 moab()->delete_entities( children, 2 );
 return MB_FAILURE;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf, Plane plane )
{
 EntityHandle left, right;
 return split_leaf( leaf, plane, left, right );
}
 
ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf,
 Plane plane,
 const Range& left_entities,
 const Range& right_entities )
{
 EntityHandle left, right, parent = leaf.handle();
 ErrorCode rval = split_leaf( leaf, plane, left, right );
 if( MB_SUCCESS != rval ) return rval;
 
 if( MB_SUCCESS == moab()->add_entities( left, left_entities ) &&
 MB_SUCCESS == moab()->add_entities( right, right_entities ) &&
 MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
 return MB_SUCCESS;
 
 moab()->remove_child_meshset( parent, left );
 moab()->remove_child_meshset( parent, right );
 EntityHandle children[] = { left, right };
 moab()->delete_entities( children, 2 );
 return MB_FAILURE;
}
 
ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf,
 Plane plane,
 const std::vector< EntityHandle >& left_entities,
 const std::vector< EntityHandle >& right_entities )
{
 EntityHandle left, right, parent = leaf.handle();
 ErrorCode rval = split_leaf( leaf, plane, left, right );
 if( MB_SUCCESS != rval ) return rval;
 
 if( MB_SUCCESS == moab()->add_entities( left, &left_entities[0], left_entities.size() ) &&
 MB_SUCCESS == moab()->add_entities( right, &right_entities[0], right_entities.size() ) &&
 MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
 return MB_SUCCESS;
 
 moab()->remove_child_meshset( parent, left );
 moab()->remove_child_meshset( parent, right );
 EntityHandle children[] = { left, right };
 moab()->delete_entities( children, 2 );
 return MB_FAILURE;
}
 
ErrorCode AdaptiveKDTree::merge_leaf( AdaptiveKDTreeIter& iter )
{
 ErrorCode rval;
 if( iter.depth() == 1 ) // at root
 return MB_FAILURE;
 
 // Move iter to parent
 
 AdaptiveKDTreeIter::StackObj node = iter.mStack.back();
 iter.mStack.pop_back();
 
 iter.childVect.clear();
 rval = moab()->get_child_meshsets( iter.mStack.back().entity, iter.childVect );
 if( MB_SUCCESS != rval ) return rval;
 Plane plane;
 rval = get_split_plane( iter.mStack.back().entity, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 int child_idx = iter.childVect[0] == node.entity ? 0 : 1;
 assert( iter.childVect[child_idx] == node.entity );
 iter.mBox[1 - child_idx][plane.norm] = node.coord;
 
 // Get all entities from children and put them in parent
 EntityHandle parent = iter.handle();
 moab()->remove_child_meshset( parent, iter.childVect[0] );
 moab()->remove_child_meshset( parent, iter.childVect[1] );
 std::vector< EntityHandle > stack( iter.childVect );
 
 Range range;
 while( !stack.empty() )
 {
 EntityHandle h = stack.back();
 stack.pop_back();
 range.clear();
 rval = moab()->get_entities_by_handle( h, range );
 if( MB_SUCCESS != rval ) return rval;
 rval = moab()->add_entities( parent, range );
 if( MB_SUCCESS != rval ) return rval;
 
 iter.childVect.clear();
 rval = moab()->get_child_meshsets( h, iter.childVect );MB_CHK_ERR( rval );
 if( !iter.childVect.empty() )
 {
 moab()->remove_child_meshset( h, iter.childVect[0] );
 moab()->remove_child_meshset( h, iter.childVect[1] );
 stack.push_back( iter.childVect[0] );
 stack.push_back( iter.childVect[1] );
 }
 
 rval = moab()->delete_entities( &h, 1 );
 if( MB_SUCCESS != rval ) return rval;
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTreeIter::initialize( AdaptiveKDTree* ttool,
 EntityHandle root,
 const double bmin[3],
 const double bmax[3],
 Direction direction )
{
 mStack.clear();
 treeTool = ttool;
 mBox[BMIN][0] = bmin[0];
 mBox[BMIN][1] = bmin[1];
 mBox[BMIN][2] = bmin[2];
 mBox[BMAX][0] = bmax[0];
 mBox[BMAX][1] = bmax[1];
 mBox[BMAX][2] = bmax[2];
 mStack.push_back( StackObj( root, 0 ) );
 return step_to_first_leaf( direction );
}
 
ErrorCode AdaptiveKDTreeIter::step_to_first_leaf( Direction direction )
{
 ErrorCode rval;
 AdaptiveKDTree::Plane plane;
 const Direction opposite = static_cast< Direction >( 1 - direction );
 
 for( ;; )
 {
 childVect.clear();
 treeTool->treeStats.nodesVisited++; // not sure whether this is the visit or the push_back below
 rval = treeTool->moab()->get_child_meshsets( mStack.back().entity, childVect );
 if( MB_SUCCESS != rval ) return rval;
 if( childVect.empty() )
 { // leaf
 treeTool->treeStats.leavesVisited++;
 break;
 }
 
 rval = treeTool->get_split_plane( mStack.back().entity, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 mStack.push_back( StackObj( childVect[direction], mBox[opposite][plane.norm] ) );
 mBox[opposite][plane.norm] = plane.coord;
 }
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTreeIter::step( Direction direction )
{
 StackObj node, parent;
 ErrorCode rval;
 AdaptiveKDTree::Plane plane;
 const Direction opposite = static_cast< Direction >( 1 - direction );
 
 // If stack is empty, then either this iterator is uninitialized
 // or we reached the end of the iteration (and return
 // MB_ENTITY_NOT_FOUND) already.
 if( mStack.empty() ) return MB_FAILURE;
 
 // Pop the current node from the stack.
 // The stack should then contain the parent of the current node.
 // If the stack is empty after this pop, then we've reached the end.
 node = mStack.back();
 mStack.pop_back();
 treeTool->treeStats.nodesVisited++;
 if( mStack.empty() ) treeTool->treeStats.leavesVisited++;
 
 while( !mStack.empty() )
 {
 // Get data for parent entity
 parent = mStack.back();
 childVect.clear();
 rval = treeTool->moab()->get_child_meshsets( parent.entity, childVect );
 if( MB_SUCCESS != rval ) return rval;
 rval = treeTool->get_split_plane( parent.entity, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 // If we're at the left child
 if( childVect[opposite] == node.entity )
 {
 // change from box of left child to box of parent
 mBox[direction][plane.norm] = node.coord;
 // push right child on stack
 node.entity = childVect[direction];
 treeTool->treeStats.nodesVisited++; // changed node
 node.coord = mBox[opposite][plane.norm];
 mStack.push_back( node );
 // change from box of parent to box of right child
 mBox[opposite][plane.norm] = plane.coord;
 // descend to left-most leaf of the right child
 return step_to_first_leaf( opposite );
 }
 
 // The current node is the right child of the parent,
 // continue up the tree.
 assert( childVect[direction] == node.entity );
 mBox[opposite][plane.norm] = node.coord;
 node = parent;
 treeTool->treeStats.nodesVisited++;
 mStack.pop_back();
 }
 
 return MB_ENTITY_NOT_FOUND;
}
 
ErrorCode AdaptiveKDTreeIter::get_neighbors( AdaptiveKDTree::Axis norm,
 bool neg,
 std::vector< AdaptiveKDTreeIter >& results,
 double epsilon ) const
{
 StackObj node, parent;
 ErrorCode rval;
 AdaptiveKDTree::Plane plane;
 int child_idx;
 
 // Find tree node at which the specified side of the box
 // for this node was created.
 AdaptiveKDTreeIter iter( *this ); // temporary iterator (don't modify *this)
 node = iter.mStack.back();
 iter.mStack.pop_back();
 for( ;; )
 {
 // reached the root - original node was on boundary (no neighbors)
 if( iter.mStack.empty() ) return MB_SUCCESS;
 
 // get parent node data
 parent = iter.mStack.back();
 iter.childVect.clear();
 rval = treeTool->moab()->get_child_meshsets( parent.entity, iter.childVect );
 if( MB_SUCCESS != rval ) return rval;
 rval = treeTool->get_split_plane( parent.entity, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 child_idx = iter.childVect[0] == node.entity ? 0 : 1;
 assert( iter.childVect[child_idx] == node.entity );
 
 // if we found the split plane for the desired side
 // push neighbor on stack and stop
 if( plane.norm == norm && (int)neg == child_idx )
 {
 // change from box of previous child to box of parent
 iter.mBox[1 - child_idx][plane.norm] = node.coord;
 // push other child of parent onto stack
 node.entity = iter.childVect[1 - child_idx];
 node.coord = iter.mBox[child_idx][plane.norm];
 iter.mStack.push_back( node );
 // change from parent box to box of new child
 iter.mBox[child_idx][plane.norm] = plane.coord;
 break;
 }
 
 // continue up the tree
 iter.mBox[1 - child_idx][plane.norm] = node.coord;
 node = parent;
 iter.mStack.pop_back();
 }
 
 // now move down tree, searching for adjacent boxes
 std::vector< AdaptiveKDTreeIter > list;
 // loop over all potential paths to neighbors (until list is empty)
 for( ;; )
 {
 // follow a single path to a leaf, append any other potential
 // paths to neighbors to 'list'
 node = iter.mStack.back();
 for( ;; )
 {
 iter.childVect.clear();
 rval = treeTool->moab()->get_child_meshsets( node.entity, iter.childVect );
 if( MB_SUCCESS != rval ) return rval;
 
 // if leaf
 if( iter.childVect.empty() )
 {
 results.push_back( iter );
 break;
 }
 
 rval = treeTool->get_split_plane( node.entity, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 // if split parallel to side
 if( plane.norm == norm )
 {
 // continue with whichever child is on the correct side of the split
 node.entity = iter.childVect[neg];
 node.coord = iter.mBox[1 - neg][plane.norm];
 iter.mStack.push_back( node );
 iter.mBox[1 - neg][plane.norm] = plane.coord;
 }
 // if left child is adjacent
 else if( this->mBox[BMIN][plane.norm] - plane.coord <= epsilon )
 {
 // if right child is also adjacent, add to list
 if( plane.coord - this->mBox[BMAX][plane.norm] <= epsilon )
 {
 list.push_back( iter );
 list.back().mStack.push_back( StackObj( iter.childVect[1], iter.mBox[BMIN][plane.norm] ) );
 list.back().mBox[BMIN][plane.norm] = plane.coord;
 }
 // continue with left child
 node.entity = iter.childVect[0];
 node.coord = iter.mBox[BMAX][plane.norm];
 iter.mStack.push_back( node );
 iter.mBox[BMAX][plane.norm] = plane.coord;
 }
 // right child is adjacent
 else
 {
 // if left child is not adjacent, right must be or something
 // is really messed up.
 assert( plane.coord - this->mBox[BMAX][plane.norm] <= epsilon );
 // continue with left child
 node.entity = iter.childVect[1];
 node.coord = iter.mBox[BMIN][plane.norm];
 iter.mStack.push_back( node );
 iter.mBox[BMIN][plane.norm] = plane.coord;
 }
 }
 
 if( list.empty() ) break;
 
 iter = list.back();
 list.pop_back();
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTreeIter::sibling_side( AdaptiveKDTree::Axis& axis_out, bool& neg_out ) const
{
 if( mStack.size() < 2 ) // at tree root
 return MB_ENTITY_NOT_FOUND;
 
 EntityHandle parent = mStack[mStack.size() - 2].entity;
 AdaptiveKDTree::Plane plane;
 ErrorCode rval = tool()->get_split_plane( parent, plane );
 if( MB_SUCCESS != rval ) return MB_FAILURE;
 
 childVect.clear();
 rval = tool()->moab()->get_child_meshsets( parent, childVect );
 if( MB_SUCCESS != rval || childVect.size() != 2 ) return MB_FAILURE;
 
 axis_out = static_cast< AdaptiveKDTree::Axis >( plane.norm );
 neg_out = ( childVect[1] == handle() );
 assert( childVect[neg_out] == handle() );
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTreeIter::get_parent_split_plane( AdaptiveKDTree::Plane& plane ) const
{
 if( mStack.size() < 2 ) // at tree root
 return MB_ENTITY_NOT_FOUND;
 
 EntityHandle parent = mStack[mStack.size() - 2].entity;
 return tool()->get_split_plane( parent, plane );
}
 
bool AdaptiveKDTreeIter::is_sibling( const AdaptiveKDTreeIter& other_leaf ) const
{
 const size_t s = mStack.size();
 return ( s > 1 ) && ( s == other_leaf.mStack.size() ) &&
 ( other_leaf.mStack[s - 2].entity == mStack[s - 2].entity ) && other_leaf.handle() != handle();
}
 
bool AdaptiveKDTreeIter::is_sibling( EntityHandle other_leaf ) const
{
 if( mStack.size() < 2 || other_leaf == handle() ) return false;
 EntityHandle parent = mStack[mStack.size() - 2].entity;
 childVect.clear();
 ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
 if( MB_SUCCESS != rval || childVect.size() != 2 )
 {
 assert( false );
 return false;
 }
 return childVect[0] == other_leaf || childVect[1] == other_leaf;
}
 
bool AdaptiveKDTreeIter::sibling_is_forward() const
{
 if( mStack.size() < 2 ) // if root
 return false;
 EntityHandle parent = mStack[mStack.size() - 2].entity;
 childVect.clear();
 ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
 if( MB_SUCCESS != rval || childVect.size() != 2 )
 {
 assert( false );
 return false;
 }
 return childVect[0] == handle();
}
 
bool AdaptiveKDTreeIter::intersect_ray( const double ray_point[3],
 const double ray_vect[3],
 double& t_enter,
 double& t_exit ) const
{
 treeTool->treeStats.traversalLeafObjectTests++;
 return GeomUtil::ray_box_intersect( CartVect( box_min() ), CartVect( box_max() ), CartVect( ray_point ),
 CartVect( ray_vect ), t_enter, t_exit );
}
 
ErrorCode AdaptiveKDTree::intersect_children_with_elems( const Range& elems,
 AdaptiveKDTree::Plane plane,
 double eps,
 CartVect box_min,
 CartVect box_max,
 Range& left_tris,
 Range& right_tris,
 Range& both_tris,
 double& metric_value )
{
 left_tris.clear();
 right_tris.clear();
 both_tris.clear();
 CartVect coords[16];
 
 // get extents of boxes for left and right sides
 BoundBox left_box( box_min, box_max ), right_box( box_min, box_max );
 right_box.bMin = box_min;
 left_box.bMax = box_max;
 right_box.bMin[plane.norm] = left_box.bMax[plane.norm] = plane.coord;
 const CartVect left_cen = 0.5 * ( left_box.bMax + box_min );
 const CartVect left_dim = 0.5 * ( left_box.bMax - box_min );
 const CartVect right_cen = 0.5 * ( box_max + right_box.bMin );
 const CartVect right_dim = 0.5 * ( box_max - right_box.bMin );
 const CartVect dim = box_max - box_min;
 const double max_tol = std::max( dim[0], std::max( dim[1], dim[2] ) ) / 10;
 
 // test each entity
 ErrorCode rval;
 int count, count2;
 const EntityHandle *conn, *conn2;
 
 const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
 const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
 const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
 Range::iterator left_ins = left_tris.begin();
 Range::iterator right_ins = right_tris.begin();
 Range::iterator both_ins = both_tris.begin();
 Range::const_iterator i;
 
 // vertices
 for( i = elems.begin(); i != elem_begin; ++i )
 {
 tree_stats().constructLeafObjectTests++;
 rval = moab()->get_coords( &*i, 1, coords[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 
 bool lo = false, ro = false;
 if( coords[0][plane.norm] <= plane.coord ) lo = true;
 if( coords[0][plane.norm] >= plane.coord ) ro = true;
 
 if( lo && ro )
 both_ins = both_tris.insert( both_ins, *i, *i );
 else if( lo )
 left_ins = left_tris.insert( left_ins, *i, *i );
 else // if (ro)
 right_ins = right_tris.insert( right_ins, *i, *i );
 }
 
 // non-polyhedron elements
 std::vector< EntityHandle > dum_vector;
 for( i = elem_begin; i != poly_begin; ++i )
 {
 tree_stats().constructLeafObjectTests++;
 rval = moab()->get_connectivity( *i, conn, count, true, &dum_vector );
 if( MB_SUCCESS != rval ) return rval;
 if( count > (int)( sizeof( coords ) / sizeof( coords[0] ) ) ) return MB_FAILURE;
 rval = moab()->get_coords( &conn[0], count, coords[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 
 bool lo = false, ro = false;
 for( int j = 0; j < count; ++j )
 {
 if( coords[j][plane.norm] <= plane.coord ) lo = true;
 if( coords[j][plane.norm] >= plane.coord ) ro = true;
 }
 
 // Triangle must be in at least one leaf. If test against plane
 // identified that leaf, then we're done. If triangle is on both
 // sides of plane, do more precise test to ensure that it is really
 // in both.
 // BoundBox box;
 // box.update(*moab(), *i);
 if( lo && ro )
 {
 double tol = eps;
 lo = ro = false;
 while( !lo && !ro && tol <= max_tol )
 {
 tree_stats().boxElemTests += 2;
 lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE( *i ), left_cen, left_dim + CartVect( tol ),
 count );
 ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE( *i ), right_cen, right_dim + CartVect( tol ),
 count );
 
 tol *= 10.0;
 }
 }
 if( lo && ro )
 both_ins = both_tris.insert( both_ins, *i, *i );
 else if( lo )
 left_ins = left_tris.insert( left_ins, *i, *i );
 else if( ro )
 right_ins = right_tris.insert( right_ins, *i, *i );
 }
 
 // polyhedra
 for( i = poly_begin; i != set_begin; ++i )
 {
 tree_stats().constructLeafObjectTests++;
 rval = moab()->get_connectivity( *i, conn, count, true );
 if( MB_SUCCESS != rval ) return rval;
 
 // just check the bounding box of the polyhedron
 bool lo = false, ro = false;
 for( int j = 0; j < count; ++j )
 {
 rval = moab()->get_connectivity( conn[j], conn2, count2, true );
 if( MB_SUCCESS != rval ) return rval;
 
 for( int k = 0; k < count2; ++k )
 {
 rval = moab()->get_coords( conn2 + k, 1, coords[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 if( coords[0][plane.norm] <= plane.coord ) lo = true;
 if( coords[0][plane.norm] >= plane.coord ) ro = true;
 }
 }
 
 if( lo && ro )
 both_ins = both_tris.insert( both_ins, *i, *i );
 else if( lo )
 left_ins = left_tris.insert( left_ins, *i, *i );
 else if( ro )
 right_ins = right_tris.insert( right_ins, *i, *i );
 }
 
 // sets
 BoundBox tbox;
 for( i = set_begin; i != elems.end(); ++i )
 {
 tree_stats().constructLeafObjectTests++;
 rval = tbox.update( *moab(), *i, spherical, radius );
 if( MB_SUCCESS != rval ) return rval;
 
 bool lo = false, ro = false;
 if( tbox.bMin[plane.norm] <= plane.coord ) lo = true;
 if( tbox.bMax[plane.norm] >= plane.coord ) ro = true;
 
 if( lo && ro )
 both_ins = both_tris.insert( both_ins, *i, *i );
 else if( lo )
 left_ins = left_tris.insert( left_ins, *i, *i );
 else // if (ro)
 right_ins = right_tris.insert( right_ins, *i, *i );
 }
 
 CartVect box_dim = box_max - box_min;
 double area_left = left_dim[0] * left_dim[1] + left_dim[1] * left_dim[2] + left_dim[2] * left_dim[0];
 double area_right = right_dim[0] * right_dim[1] + right_dim[1] * right_dim[2] + right_dim[2] * right_dim[0];
 double area_both = box_dim[0] * box_dim[1] + box_dim[1] * box_dim[2] + box_dim[2] * box_dim[0];
 metric_value = ( area_left * left_tris.size() + area_right * right_tris.size() ) / area_both + both_tris.size();
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::best_subdivision_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 double eps )
{
 double metric_val = std::numeric_limits< unsigned >::max();
 
 ErrorCode r;
 const CartVect box_min( iter.box_min() );
 const CartVect box_max( iter.box_max() );
 const CartVect diff( box_max - box_min );
 
 Range entities;
 r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
 if( MB_SUCCESS != r ) return r;
 const size_t p_count = entities.size();
 
 for( int axis = 0; axis < 3; ++axis )
 {
 int plane_count = num_planes;
 if( ( num_planes + 1 ) * eps >= diff[axis] ) plane_count = (int)( diff[axis] / eps ) - 1;
 
 for( int p = 1; p <= plane_count; ++p )
 {
 AdaptiveKDTree::Plane plane = { box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis], axis };
 Range left, right, both;
 double val;
 r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
 if( MB_SUCCESS != r ) return r;
 const size_t sdiff = p_count - both.size();
 if( left.size() == sdiff || right.size() == sdiff ) continue;
 
 if( val >= metric_val ) continue;
 
 metric_val = val;
 best_plane = plane;
 best_left.swap( left );
 best_right.swap( right );
 best_both.swap( both );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::best_subdivision_snap_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 std::vector< double >& tmp_data,
 double eps )
{
 double metric_val = std::numeric_limits< unsigned >::max();
 
 ErrorCode r;
 // const CartVect tol(eps*diff);
 
 Range entities, vertices;
 r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
 if( MB_SUCCESS != r ) return r;
 const size_t p_count = entities.size();
 r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
 if( MB_SUCCESS != r ) return r;
 
 unsigned int nverts = vertices.size();
 tmp_data.resize( 3 * nverts );
 r = iter.tool()->moab()->get_coords( vertices, &tmp_data[0], &tmp_data[nverts], &tmp_data[2 * nverts] );
 if( MB_SUCCESS != r ) return r;
 
 // calculate bounding box of vertices
 // decide based on the actual box the splitting plane
 // do not decide based on iterator box.
 // it could be too big
 // BoundBox box;
 // r = box.update(*moab(), vertices);
 CartVect box_min;
 CartVect box_max;
 for( int dir = 0; dir < 3; dir++ )
 {
 double amin = tmp_data[dir * nverts];
 double amax = amin;
 double* p = &tmp_data[dir * nverts + 1];
 for( unsigned int i = 1; i < nverts; i++ )
 {
 if( *p < amin ) amin = *p;
 if( *p > amax ) amax = *p;
 p++;
 }
 box_min[dir] = amin;
 box_max[dir] = amax;
 }
 CartVect diff( box_max - box_min );
 
 for( int axis = 0; axis < 3; ++axis )
 {
 int plane_count = num_planes;
 
 // if num_planes results in width < eps, reset the plane count
 if( ( num_planes + 1 ) * eps >= diff[axis] ) plane_count = (int)( diff[axis] / eps ) - 1;
 
 for( int p = 1; p <= plane_count; ++p )
 {
 
 // coord of this plane on axis
 double coord = box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis];
 
 // find closest vertex coordinate to this plane position
 unsigned int istrt = axis * nverts;
 double closest_coord = tmp_data[istrt];
 for( unsigned i = 1; i < nverts; ++i )
 if( fabs( coord - tmp_data[istrt + i] ) < fabs( coord - closest_coord ) )
 closest_coord = tmp_data[istrt + i];
 if( closest_coord - box_min[axis] <= eps || box_max[axis] - closest_coord <= eps ) continue;
 
 // seprate elems into left/right/both, and compute separating metric
 AdaptiveKDTree::Plane plane = { closest_coord, axis };
 Range left, right, both;
 double val;
 r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
 if( MB_SUCCESS != r ) return r;
 const size_t d = p_count - both.size();
 if( left.size() == d || right.size() == d ) continue;
 
 if( val >= metric_val ) continue;
 
 metric_val = val;
 best_plane = plane;
 best_left.swap( left );
 best_right.swap( right );
 best_both.swap( both );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::best_vertex_median_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 std::vector< double >& coords,
 double eps )
{
 double metric_val = std::numeric_limits< unsigned >::max();
 
 ErrorCode r;
 const CartVect box_min( iter.box_min() );
 const CartVect box_max( iter.box_max() );
 
 Range entities, vertices;
 r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
 if( MB_SUCCESS != r ) return r;
 const size_t p_count = entities.size();
 r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
 if( MB_SUCCESS != r ) return r;
 
 coords.resize( vertices.size() );
 for( int axis = 0; axis < 3; ++axis )
 {
 if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
 
 double* ptrs[] = { 0, 0, 0 };
 ptrs[axis] = &coords[0];
 r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
 if( MB_SUCCESS != r ) return r;
 
 std::sort( coords.begin(), coords.end() );
 std::vector< double >::iterator citer;
 citer = std::upper_bound( coords.begin(), coords.end(), box_min[axis] + eps );
 const size_t count = std::upper_bound( citer, coords.end(), box_max[axis] - eps ) - citer;
 size_t step;
 int np = num_planes;
 if( count < 2 * (size_t)num_planes )
 {
 step = 1;
 np = count - 1;
 }
 else
 {
 step = count / ( num_planes + 1 );
 }
 
 for( int p = 1; p <= np; ++p )
 {
 
 citer += step;
 AdaptiveKDTree::Plane plane = { *citer, axis };
 Range left, right, both;
 double val;
 r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
 if( MB_SUCCESS != r ) return r;
 const size_t diff = p_count - both.size();
 if( left.size() == diff || right.size() == diff ) continue;
 
 if( val >= metric_val ) continue;
 
 metric_val = val;
 best_plane = plane;
 best_left.swap( left );
 best_right.swap( right );
 best_both.swap( both );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::best_vertex_sample_plane( int num_planes,
 const AdaptiveKDTreeIter& iter,
 Range& best_left,
 Range& best_right,
 Range& best_both,
 AdaptiveKDTree::Plane& best_plane,
 std::vector< double >& coords,
 std::vector< EntityHandle >& indices,
 double eps )
{
 const size_t random_elem_threshold = 20 * num_planes;
 double metric_val = std::numeric_limits< unsigned >::max();
 
 ErrorCode r;
 const CartVect box_min( iter.box_min() );
 const CartVect box_max( iter.box_max() );
 
 Range entities, vertices;
 r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
 if( MB_SUCCESS != r ) return r;
 
 // We are selecting random vertex coordinates to use for candidate split
 // planes. So if element list is large, begin by selecting random elements.
 const size_t p_count = entities.size();
 coords.resize( 3 * num_planes );
 if( p_count < random_elem_threshold )
 {
 r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
 if( MB_SUCCESS != r ) return r;
 }
 else
 {
 indices.resize( random_elem_threshold );
 const int num_rand = p_count / RAND_MAX + 1;
 for( size_t j = 0; j < random_elem_threshold; ++j )
 {
 size_t rnd = rand();
 for( int i = num_rand; i > 1; --i )
 rnd *= rand();
 rnd %= p_count;
 indices[j] = entities[rnd];
 }
 r = iter.tool()->moab()->get_adjacencies( &indices[0], random_elem_threshold, 0, false, vertices,
 Interface::UNION );
 if( MB_SUCCESS != r ) return r;
 }
 
 coords.resize( vertices.size() );
 for( int axis = 0; axis < 3; ++axis )
 {
 if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
 
 double* ptrs[] = { 0, 0, 0 };
 ptrs[axis] = &coords[0];
 r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
 if( MB_SUCCESS != r ) return r;
 
 size_t num_valid_coords = 0;
 for( size_t i = 0; i < coords.size(); ++i )
 if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) ++num_valid_coords;
 
 if( 2 * (size_t)num_planes > num_valid_coords )
 {
 indices.clear();
 for( size_t i = 0; i < coords.size(); ++i )
 if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) indices.push_back( i );
 }
 else
 {
 indices.resize( num_planes );
 // make sure random indices are sufficient to cover entire range
 const int num_rand = coords.size() / RAND_MAX + 1;
 for( int j = 0; j < num_planes; ++j )
 {
 size_t rnd;
 do
 {
 rnd = rand();
 for( int i = num_rand; i > 1; --i )
 rnd *= rand();
 rnd %= coords.size();
 } while( coords[rnd] <= box_min[axis] + eps || coords[rnd] >= box_max[axis] - eps );
 indices[j] = rnd;
 }
 }
 
 for( unsigned p = 0; p < indices.size(); ++p )
 {
 
 AdaptiveKDTree::Plane plane = { coords[indices[p]], axis };
 Range left, right, both;
 double val;
 r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
 if( MB_SUCCESS != r ) return r;
 const size_t diff = p_count - both.size();
 if( left.size() == diff || right.size() == diff ) continue;
 
 if( val >= metric_val ) continue;
 
 metric_val = val;
 best_plane = plane;
 best_left.swap( left );
 best_right.swap( right );
 best_both.swap( both );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::point_search( const double* point,
 EntityHandle& leaf_out,
 const double iter_tol,
 const double inside_tol,
 bool* multiple_leaves,
 EntityHandle* start_node,
 CartVect* params )
{
 std::vector< EntityHandle > children;
 Plane plane;
 
 treeStats.numTraversals++;
 leaf_out = 0;
 BoundBox box;
 // kdtrees never have multiple leaves containing a pt
 if( multiple_leaves ) *multiple_leaves = false;
 
 EntityHandle node = ( start_node ? *start_node : myRoot );
 
 treeStats.nodesVisited++;
 ErrorCode rval = get_bounding_box( box, &node );
 if( MB_SUCCESS != rval ) return rval;
 if( !box.contains_point( point, iter_tol ) ) return MB_SUCCESS;
 
 rval = moab()->get_child_meshsets( node, children );
 if( MB_SUCCESS != rval ) return rval;
 
 while( !children.empty() )
 {
 treeStats.nodesVisited++;
 
 rval = get_split_plane( node, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 const double d = point[plane.norm] - plane.coord;
 node = children[( d > 0.0 )];
 
 children.clear();
 rval = moab()->get_child_meshsets( node, children );
 if( MB_SUCCESS != rval ) return rval;
 }
 
 treeStats.leavesVisited++;
 if( myEval && params )
 {
 rval = myEval->find_containing_entity( node, point, iter_tol, inside_tol, leaf_out, params->array(),
 &treeStats.traversalLeafObjectTests );
 if( MB_SUCCESS != rval ) return rval;
 }
 else
 leaf_out = node;
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::point_search( const double* point,
 AdaptiveKDTreeIter& leaf_it,
 const double iter_tol,
 const double /*inside_tol*/,
 bool* multiple_leaves,
 EntityHandle* start_node )
{
 ErrorCode rval;
 treeStats.numTraversals++;
 
 // kdtrees never have multiple leaves containing a pt
 if( multiple_leaves ) *multiple_leaves = false;
 
 leaf_it.mBox[0] = boundBox.bMin;
 leaf_it.mBox[1] = boundBox.bMax;
 
 // test that point is inside tree
 if( !boundBox.contains_point( point, iter_tol ) )
 {
 treeStats.nodesVisited++;
 return MB_ENTITY_NOT_FOUND;
 }
 
 // initialize iterator at tree root
 leaf_it.treeTool = this;
 leaf_it.mStack.clear();
 leaf_it.mStack.push_back( AdaptiveKDTreeIter::StackObj( ( start_node ? *start_node : myRoot ), 0 ) );
 
 // loop until we reach a leaf
 AdaptiveKDTree::Plane plane;
 for( ;; )
 {
 treeStats.nodesVisited++;
 
 // get children
 leaf_it.childVect.clear();
 rval = moab()->get_child_meshsets( leaf_it.handle(), leaf_it.childVect );
 if( MB_SUCCESS != rval ) return rval;
 
 // if no children, then at leaf (done)
 if( leaf_it.childVect.empty() )
 {
 treeStats.leavesVisited++;
 break;
 }
 
 // get split plane
 rval = get_split_plane( leaf_it.handle(), plane );
 if( MB_SUCCESS != rval ) return rval;
 
 // step iterator to appropriate child
 // idx: 0->left, 1->right
 const int idx = ( point[plane.norm] > plane.coord );
 leaf_it.mStack.push_back(
 AdaptiveKDTreeIter::StackObj( leaf_it.childVect[idx], leaf_it.mBox[1 - idx][plane.norm] ) );
 leaf_it.mBox[1 - idx][plane.norm] = plane.coord;
 }
 
 return MB_SUCCESS;
}
 
struct NodeDistance
{
 EntityHandle handle;
 CartVect dist; // from_point - closest_point_on_box
};
 
ErrorCode AdaptiveKDTree::distance_search( const double from_point[3],
 const double distance,
 std::vector< EntityHandle >& result_list,
 const double iter_tol,
 const double inside_tol,
 std::vector< double >* result_dists,
 std::vector< CartVect >* result_params,
 EntityHandle* tree_root )
{
 treeStats.numTraversals++;
 const double dist_sqr = distance * distance;
 const CartVect from( from_point );
 std::vector< NodeDistance > list,
 result_list_nodes; // list of subtrees to traverse, and results
 // pre-allocate space for default max tree depth
 list.reserve( maxDepth );
 
 // misc temporary values
 Plane plane;
 NodeDistance node;
 ErrorCode rval;
 std::vector< EntityHandle > children;
 
 // Get distance from input position to bounding box of tree
 // (zero if inside box)
 BoundBox box;
 rval = get_bounding_box( box );
 if( MB_SUCCESS == rval && !box.contains_point( from_point, iter_tol ) )
 {
 treeStats.nodesVisited++;
 return MB_SUCCESS;
 }
 
 // if bounding box is not available (e.g. not starting from true root)
 // just start with zero. Less efficient, but will work.
 node.dist = CartVect( 0.0 );
 if( MB_SUCCESS == rval )
 {
 for( int i = 0; i < 3; ++i )
 {
 if( from_point[i] < box.bMin[i] )
 node.dist[i] = box.bMin[i] - from_point[i];
 else if( from_point[i] > box.bMax[i] )
 node.dist[i] = from_point[i] - box.bMax[i];
 }
 if( node.dist % node.dist > dist_sqr )
 {
 treeStats.nodesVisited++;
 return MB_SUCCESS;
 }
 }
 
 // begin with root in list
 node.handle = ( tree_root ? *tree_root : myRoot );
 list.push_back( node );
 
 while( !list.empty() )
 {
 
 node = list.back();
 list.pop_back();
 treeStats.nodesVisited++;
 
 // If leaf node, test contained triangles
 children.clear();
 rval = moab()->get_child_meshsets( node.handle, children );
 if( children.empty() )
 {
 treeStats.leavesVisited++;
 if( myEval && result_params )
 {
 EntityHandle ent;
 CartVect params;
 rval = myEval->find_containing_entity( node.handle, 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
 {
 result_list_nodes.push_back( node );
 continue;
 }
 }
 
 // If not leaf node, add children to working list
 rval = get_split_plane( node.handle, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 const double d = from[plane.norm] - plane.coord;
 
 // right of plane?
 if( d > 0 )
 {
 node.handle = children[1];
 list.push_back( node );
 // if the split plane is close to the input point, add
 // the left child also (we'll check the exact distance
 /// when we pop it from the list.)
 if( d <= distance )
 {
 node.dist[plane.norm] = d;
 if( node.dist % node.dist <= dist_sqr )
 {
 node.handle = children[0];
 list.push_back( node );
 }
 }
 }
 // left of plane
 else
 {
 node.handle = children[0];
 list.push_back( node );
 // if the split plane is close to the input point, add
 // the right child also (we'll check the exact distance
 /// when we pop it from the list.)
 if( -d <= distance )
 {
 node.dist[plane.norm] = -d;
 if( node.dist % node.dist <= dist_sqr )
 {
 node.handle = children[1];
 list.push_back( node );
 }
 }
 }
 }
 
 if( myEval && result_params ) return MB_SUCCESS;
 
 // separate loops to avoid if test inside loop
 
 result_list.reserve( result_list_nodes.size() );
 for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end(); ++vit )
 result_list.push_back( ( *vit ).handle );
 
 if( result_dists && distance > 0.0 )
 {
 result_dists->reserve( result_list_nodes.size() );
 for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end();
 ++vit )
 result_dists->push_back( ( *vit ).dist.length() );
 }
 
 return MB_SUCCESS;
}
 
static ErrorCode closest_to_triangles( Interface* moab,
 const Range& tris,
 const CartVect& from,
 double& shortest_dist_sqr,
 CartVect& closest_pt,
 EntityHandle& closest_tri )
{
 ErrorCode rval;
 CartVect pos, diff, verts[3];
 const EntityHandle* conn = NULL;
 int len = 0;
 
 for( Range::iterator i = tris.begin(); i != tris.end(); ++i )
 {
 rval = moab->get_connectivity( *i, conn, len );
 if( MB_SUCCESS != rval ) return rval;
 
 rval = moab->get_coords( conn, 3, verts[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 
 GeomUtil::closest_location_on_tri( from, verts, pos );
 diff = pos - from;
 double dist_sqr = diff % diff;
 if( dist_sqr < shortest_dist_sqr )
 {
 // new closest location
 shortest_dist_sqr = dist_sqr;
 closest_pt = pos;
 closest_tri = *i;
 }
 }
 
 return MB_SUCCESS;
}
 
static ErrorCode closest_to_triangles( Interface* moab,
 EntityHandle set_handle,
 const CartVect& from,
 double& shortest_dist_sqr,
 CartVect& closest_pt,
 EntityHandle& closest_tri )
{
 ErrorCode rval;
 Range tris;
 
 rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
 if( MB_SUCCESS != rval ) return rval;
 
 return closest_to_triangles( moab, tris, from, shortest_dist_sqr, closest_pt, closest_tri );
}
 
ErrorCode AdaptiveKDTree::find_close_triangle( EntityHandle root,
 const double from[3],
 double pt[3],
 EntityHandle& triangle )
{
 ErrorCode rval;
 Range tris;
 Plane split;
 std::vector< EntityHandle > stack;
 std::vector< EntityHandle > children( 2 );
 stack.reserve( 30 );
 assert( root );
 stack.push_back( root );
 
 while( !stack.empty() )
 {
 EntityHandle node = stack.back();
 stack.pop_back();
 
 for( ;; )
 { // loop until we find a leaf
 
 children.clear();
 rval = moab()->get_child_meshsets( node, children );
 if( MB_SUCCESS != rval ) return rval;
 
 // loop termination criterion
 if( children.empty() ) break;
 
 // if not a leaf, get split plane
 rval = get_split_plane( node, split );
 if( MB_SUCCESS != rval ) return rval;
 
 // continue down the side that contains the point,
 // and push the other side onto the stack in case
 // we need to check it later.
 int rs = split.right_side( from );
 node = children[rs];
 stack.push_back( children[1 - rs] );
 }
 
 // We should now be at a leaf.
 // If it has some triangles, we're done.
 // If not, continue searching for another leaf.
 tris.clear();
 rval = moab()->get_entities_by_type( node, MBTRI, tris );
 if( !tris.empty() )
 {
 double dist_sqr = HUGE_VAL;
 CartVect point( pt );
 rval = closest_to_triangles( moab(), tris, CartVect( from ), dist_sqr, point, triangle );
 point.get( pt );
 return rval;
 }
 }
 
 // If we got here, then we traversed the entire tree
 // and all the leaves were empty.
 return MB_ENTITY_NOT_FOUND;
}
 
/** Find the triangles in a set that are closer to the input
 * position than any triangles in the 'closest_tris' list.
 *
 * closest_tris is assumed to contain a list of triangles for
 * which the first is the closest known triangle to the input
 * position and the first entry in 'closest_pts' is the closest
 * location on that triangle. Any other values in the lists must
 * be other triangles for which the closest point is within the
 * input tolerance of the closest closest point. This function
 * will update the lists as appropriate if any closer triangles
 * or triangles within the tolerance of the current closest location
 * are found. The first entry is maintained as the closest of the
 * list of triangles.
 */
/*
 static ErrorCode closest_to_triangles( Interface* moab,
 EntityHandle set_handle,
 double tolerance,
 const CartVect& from,
 std::vector<EntityHandle>& closest_tris,
 std::vector<CartVect>& closest_pts )
 {
 ErrorCode rval;
 Range tris;
 CartVect pos, diff, verts[3];
 const EntityHandle* conn;
 int len;
 double shortest_dist_sqr = HUGE_VAL;
 if (!closest_pts.empty()) {
 diff = from - closest_pts.front();
 shortest_dist_sqr = diff % diff;
 }
 
 rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
 if (MB_SUCCESS != rval)
 return rval;
 
 for (Range::iterator i = tris.begin(); i != tris.end(); ++i) {
 rval = moab->get_connectivity( *i, conn, len );
 if (MB_SUCCESS != rval)
 return rval;
 
 rval = moab->get_coords( conn, 3, verts[0].array() );
 if (MB_SUCCESS != rval)
 return rval;
 
 GeomUtil::closest_location_on_tri( from, verts, pos );
 diff = pos - from;
 double dist_sqr = diff % diff;
 if (dist_sqr < shortest_dist_sqr) {
 // new closest location
 shortest_dist_sqr = dist_sqr;
 
 if (closest_pts.empty()) {
 closest_tris.push_back( *i );
 closest_pts.push_back( pos );
 }
 // if have a previous closest location
 else {
 // if previous closest is more than 2*tolerance away
 // from new closest, then nothing in the list can
 // be within tolerance of new closest point.
 diff = pos - closest_pts.front();
 dist_sqr = diff % diff;
 if (dist_sqr > 4.0 * tolerance * tolerance) {
 closest_tris.clear();
 closest_pts.clear();
 closest_tris.push_back( *i );
 closest_pts.push_back( pos );
 }
 // otherwise need to remove any triangles that are
 // not within tolerance of the new closest point.
 else {
 unsigned r = 0, w = 0;
 for (r = 0; r < closest_pts.size(); ++r) {
 diff = pos - closest_pts[r];
 if (diff % diff <= tolerance*tolerance) {
 closest_pts[w] = closest_pts[r];
 closest_tris[w] = closest_tris[r];
 ++w;
 }
 }
 closest_pts.resize( w + 1 );
 closest_tris.resize( w + 1 );
 // always put the closest one in the front
 if (w > 0) {
 closest_pts.back() = closest_pts.front();
 closest_tris.back() = closest_tris.front();
 }
 closest_pts.front() = pos;
 closest_tris.front() = *i;
 }
 }
 }
 else {
 // If within tolerance of old closest triangle,
 // add this one to the list.
 diff = closest_pts.front() - pos;
 if (diff % diff <= tolerance*tolerance) {
 closest_pts.push_back( pos );
 closest_tris.push_back( *i );
 }
 }
 }
 
 return MB_SUCCESS;
 }
*/
 
ErrorCode AdaptiveKDTree::closest_triangle( EntityHandle tree_root,
 const double from_coords[3],
 double closest_point_out[3],
 EntityHandle& triangle_out )
{
 ErrorCode rval;
 double shortest_dist_sqr = HUGE_VAL;
 std::vector< EntityHandle > leaves;
 const CartVect from( from_coords );
 CartVect closest_pt;
 
 // Find the leaf containing the input point
 // This search does not take into account any bounding box for the
 // tree, so it always returns one leaf.
 assert( tree_root );
 rval = find_close_triangle( tree_root, from_coords, closest_pt.array(), triangle_out );
 if( MB_SUCCESS != rval ) return rval;
 
 // Find any other leaves for which the bounding box is within
 // the same distance from the input point as the current closest
 // point is.
 CartVect diff = closest_pt - from;
 rval = distance_search( from_coords, sqrt( diff % diff ), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
 if( MB_SUCCESS != rval ) return rval;
 
 // Check any close leaves to see if they contain triangles that
 // are as close to or closer than the current closest triangle(s).
 for( unsigned i = 0; i < leaves.size(); ++i )
 {
 rval = closest_to_triangles( moab(), leaves[i], from, shortest_dist_sqr, closest_pt, triangle_out );
 if( MB_SUCCESS != rval ) return rval;
 }
 
 // pass back resulting position
 closest_pt.get( closest_point_out );
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::sphere_intersect_triangles( EntityHandle tree_root,
 const double center[3],
 double rad,
 std::vector< EntityHandle >& triangles )
{
 ErrorCode rval;
 std::vector< EntityHandle > leaves;
 const CartVect from( center );
 CartVect closest_pt;
 const EntityHandle* conn;
 CartVect coords[3];
 int conn_len;
 
 // get leaves of tree that intersect sphere
 assert( tree_root );
 rval = distance_search( center, rad, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
 if( MB_SUCCESS != rval ) return rval;
 
 // search each leaf for triangles intersecting sphere
 for( unsigned i = 0; i < leaves.size(); ++i )
 {
 Range tris;
 rval = moab()->get_entities_by_type( leaves[i], MBTRI, tris );
 if( MB_SUCCESS != rval ) return rval;
 
 for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
 {
 rval = moab()->get_connectivity( *j, conn, conn_len );
 if( MB_SUCCESS != rval ) return rval;
 rval = moab()->get_coords( conn, 3, coords[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 GeomUtil::closest_location_on_tri( from, coords, closest_pt );
 closest_pt -= from;
 if( ( closest_pt % closest_pt ) <= ( rad * rad ) ) triangles.push_back( *j );
 }
 }
 
 // remove duplicates from triangle list
 std::sort( triangles.begin(), triangles.end() );
 triangles.erase( std::unique( triangles.begin(), triangles.end() ), triangles.end() );
 return MB_SUCCESS;
}
 
struct NodeSeg
{
 NodeSeg( EntityHandle h, double b, double e ) : handle( h ), beg( b ), end( e ) {}
 EntityHandle handle;
 double beg, end;
};
 
ErrorCode AdaptiveKDTree::ray_intersect_triangles( EntityHandle root,
 const double tol,
 const double ray_dir_in[3],
 const double ray_pt_in[3],
 std::vector< EntityHandle >& tris_out,
 std::vector< double >& dists_out,
 int max_ints,
 double ray_end )
{
 ErrorCode rval;
 double ray_beg = 0.0;
 if( ray_end < 0.0 ) ray_end = HUGE_VAL;
 
 // if root has bounding box, trim ray to that box
 CartVect tvec( tol );
 BoundBox box;
 const CartVect ray_pt( ray_pt_in ), ray_dir( ray_dir_in );
 rval = get_bounding_box( box );
 if( MB_SUCCESS == rval )
 {
 if( !GeomUtil::segment_box_intersect( box.bMin - tvec, box.bMax + tvec, ray_pt, ray_dir, ray_beg, ray_end ) )
 return MB_SUCCESS; // ray misses entire tree.
 }
 
 Range tris;
 Range::iterator iter;
 CartVect tri_coords[3];
 const EntityHandle* tri_conn;
 int conn_len;
 double tri_t;
 
 Plane plane;
 std::vector< EntityHandle > children;
 std::vector< NodeSeg > list;
 NodeSeg seg( root, ray_beg, ray_end );
 list.push_back( seg );
 
 while( !list.empty() )
 {
 seg = list.back();
 list.pop_back();
 
 // If we are limited to a certain number of intersections
 // (max_ints != 0), then ray_end will contain the distance
 // to the furthest intersection we have so far. If the
 // tree node is further than that, skip it.
 if( seg.beg > ray_end ) continue;
 
 // Check if at a leaf
 children.clear();
 rval = moab()->get_child_meshsets( seg.handle, children );
 if( MB_SUCCESS != rval ) return rval;
 if( children.empty() )
 { // leaf
 
 tris.clear();
 rval = moab()->get_entities_by_type( seg.handle, MBTRI, tris );
 if( MB_SUCCESS != rval ) return rval;
 
 for( iter = tris.begin(); iter != tris.end(); ++iter )
 {
 rval = moab()->get_connectivity( *iter, tri_conn, conn_len );
 if( MB_SUCCESS != rval ) return rval;
 rval = moab()->get_coords( tri_conn, 3, tri_coords[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 
 if( GeomUtil::ray_tri_intersect( tri_coords, ray_pt, ray_dir, tri_t, &ray_end ) )
 {
 if( !max_ints )
 {
 if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
 {
 tris_out.push_back( *iter );
 dists_out.push_back( tri_t );
 }
 }
 else if( tri_t < ray_end )
 {
 if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
 {
 if( tris_out.size() < (unsigned)max_ints )
 {
 tris_out.resize( tris_out.size() + 1 );
 dists_out.resize( dists_out.size() + 1 );
 }
 int w = tris_out.size() - 1;
 for( ; w > 0 && tri_t < dists_out[w - 1]; --w )
 {
 tris_out[w] = tris_out[w - 1];
 dists_out[w] = dists_out[w - 1];
 }
 tris_out[w] = *iter;
 dists_out[w] = tri_t;
 if( tris_out.size() >= (unsigned)max_ints )
 // when we have already reached the max intx points, we cans safely
 // reset ray_end, because we will accept new points only "closer"
 // than the last one
 ray_end = dists_out.back();
 }
 }
 }
 }
 
 continue;
 }
 
 rval = get_split_plane( seg.handle, plane );
 if( MB_SUCCESS != rval ) return rval;
 
 // Consider two planes that are the split plane +/- the tolerance.
 // Calculate the segment parameter at which the line segment intersects
 // the true plane, and also the difference between that value and the
 // intersection with either of the +/- tol planes.
 const double inv_dir = 1.0 / ray_dir[plane.norm]; // only do division once
 const double t = ( plane.coord - ray_pt[plane.norm] ) * inv_dir; // intersection with plane
 const double diff = tol * inv_dir; // t adjustment for +tol plane
 // const double t0 = t - diff; // intersection with -tol plane
 // const double t1 = t + diff; // intersection with +tol plane
 
 // The index of the child tree node (0 or 1) that is on the
 // side of the plane to which the ray direction points. That is,
 // if the ray direction is opposite the plane normal, the index
 // of the child corresponding to the side beneath the plane. If
 // the ray direction is the same as the plane normal, the index
 // of the child corresponding to the side above the plane.
 const int fwd_child = ( ray_dir[plane.norm] > 0.0 );
 
 // Note: we maintain seg.beg <= seg.end at all times, so assume that here.
 
 // If segment is parallel to plane
 if( !Util::is_finite( t ) )
 {
 if( ray_pt[plane.norm] - tol <= plane.coord ) list.push_back( NodeSeg( children[0], seg.beg, seg.end ) );
 if( ray_pt[plane.norm] + tol >= plane.coord ) list.push_back( NodeSeg( children[1], seg.beg, seg.end ) );
 }
 // If segment is entirely to one side of plane such that the
 // intersection with the split plane is past the end of the segment
 else if( seg.end + diff < t )
 {
 // If segment direction is opposite that of plane normal, then
 // being past the end of the segment means that we are to the
 // right (or above) the plane and what the right child (index == 1).
 // Otherwise we want the left child (index == 0);
 list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
 }
 // If the segment is entirely to one side of the plane such that
 // the intersection with the split plane is before the start of the
 // segment
 else if( seg.beg - diff > t )
 {
 // If segment direction is opposite that of plane normal, then
 // being before the start of the segment means that we are to the
 // left (or below) the plane and what the left child (index == 0).
 // Otherwise we want the right child (index == 1);
 list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
 }
 // Otherwise we must intersect the plane.
 // Note: be careful not to grow the segment if t is slightly
 // outside the current segment, as doing so would effectively
 // increase the tolerance as we descend the tree.
 else if( t <= seg.beg )
 {
 list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.beg ) );
 list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
 }
 else if( t >= seg.end )
 {
 list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
 list.push_back( NodeSeg( children[fwd_child], seg.end, seg.end ) );
 }
 else
 {
 list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, t ) );
 list.push_back( NodeSeg( children[fwd_child], t, seg.end ) );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::compute_depth( EntityHandle root, unsigned int& min_depth, unsigned int& max_depth )
{
 AdaptiveKDTreeIter iter;
 get_tree_iterator( root, iter );
 iter.step_to_first_leaf( AdaptiveKDTreeIter::LEFT );
 min_depth = max_depth = iter.depth();
 
 int num_of_elements = 0, max, min;
 moab()->get_number_entities_by_handle( iter.handle(), num_of_elements );
 max = min = num_of_elements;
 int k = 0;
 while( MB_SUCCESS == iter.step() )
 {
 int temp = 0;
 moab()->get_number_entities_by_handle( iter.handle(), temp );
 max = std::max( max, temp );
 min = std::min( min, temp );
 if( iter.depth() > max_depth )
 max_depth = iter.depth();
 else if( iter.depth() < min_depth )
 min_depth = iter.depth();
 ++k;
 }
 return MB_SUCCESS;
}
 
ErrorCode AdaptiveKDTree::get_info( EntityHandle root, double bmin[3], double bmax[3], unsigned int& dep )
{
 BoundBox box;
 ErrorCode result = get_bounding_box( box, &root );
 if( MB_SUCCESS != result ) return result;
 box.bMin.get( bmin );
 box.bMax.get( bmax );
 
 unsigned min_depth;
 return compute_depth( root, min_depth, dep );
}
 
static std::string mem_to_string( unsigned long mem )
{
 char unit[3] = "B";
 if( mem > 9 * 1024 )
 {
 mem = ( mem + 512 ) / 1024;
 strcpy( unit, "kB" );
 }
 if( mem > 9 * 1024 )
 {
 mem = ( mem + 512 ) / 1024;
 strcpy( unit, "MB" );
 }
 if( mem > 9 * 1024 )
 {
 mem = ( mem + 512 ) / 1024;
 strcpy( unit, "GB" );
 }
 char buffer[256];
 sprintf( buffer, "%lu %s", mem, unit );
 return buffer;
}
 
template < typename T >
struct SimpleStat
{
 T min, max, sum, sqr;
 size_t count;
 SimpleStat();
 void add( T value );
 double avg() const
 {
 return (double)sum / count;
 }
 double rms() const
 {
 return sqrt( (double)sqr / count );
 }
 double dev() const
 {
 return ( count > 1
 ? sqrt( ( count * (double)sqr - (double)sum * (double)sum ) / ( (double)count * ( count - 1 ) ) )
 : 0.0 );
 }
};
 
template < typename T >
SimpleStat< T >::SimpleStat()
 : min( std::numeric_limits< T >::max() ), max( std::numeric_limits< T >::min() ), sum( 0 ), sqr( 0 ), count( 0 )
{
}
 
template < typename T >
void SimpleStat< T >::add( T value )
{
 if( value < min ) min = value;
 if( value > max ) max = value;
 sum += value;
 sqr += value * value;
 ++count;
}
 
ErrorCode AdaptiveKDTree::print()
{
 Range range;
 
 Range tree_sets, elem2d, elem3d, verts, all;
 moab()->get_child_meshsets( myRoot, tree_sets, 0 );
 for( Range::iterator rit = tree_sets.begin(); rit != tree_sets.end(); ++rit )
 {
 moab()->get_entities_by_dimension( *rit, 2, elem2d );
 moab()->get_entities_by_dimension( *rit, 3, elem3d );
 moab()->get_entities_by_type( *rit, MBVERTEX, verts );
 }
 all.merge( verts );
 all.merge( elem2d );
 all.merge( elem3d );
 tree_sets.insert( myRoot );
 unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized,
 elem_used, elem_amortized;
 moab()->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0, 0,
 0, &set_tag_used, &set_tag_amortized );
 moab()->estimated_memory_use( all, &elem_used, &elem_amortized );
 
 int num_2d = 0, num_3d = 0;
 ;
 moab()->get_number_entities_by_dimension( 0, 2, num_2d );
 moab()->get_number_entities_by_dimension( 0, 3, num_3d );
 
 BoundBox box;
 ErrorCode rval = get_bounding_box( box, &myRoot );
 if( MB_SUCCESS != rval || box == BoundBox() ) throw rval;
 double diff[3] = { box.bMax[0] - box.bMin[0], box.bMax[1] - box.bMin[1], box.bMax[2] - box.bMin[2] };
 double tree_vol = diff[0] * diff[1] * diff[2];
 double tree_surf_area = 2 * ( diff[0] * diff[1] + diff[1] * diff[2] + diff[2] * diff[0] );
 
 SimpleStat< unsigned > depth, size;
 SimpleStat< double > vol, surf;
 
 AdaptiveKDTreeIter iter;
 get_tree_iterator( myRoot, iter );
 do
 {
 depth.add( iter.depth() );
 
 int num_leaf_elem;
 moab()->get_number_entities_by_handle( iter.handle(), num_leaf_elem );
 size.add( num_leaf_elem );
 
 const double* n = iter.box_min();
 const double* x = iter.box_max();
 double dims[3] = { x[0] - n[0], x[1] - n[1], x[2] - n[2] };
 
 double leaf_vol = dims[0] * dims[1] * dims[2];
 vol.add( leaf_vol );
 
 double area = 2.0 * ( dims[0] * dims[1] + dims[1] * dims[2] + dims[2] * dims[0] );
 surf.add( area );
 
 } while( MB_SUCCESS == iter.step() );
 
 printf( "------------------------------------------------------------------\n" );
 printf( "tree volume: %f\n", tree_vol );
 printf( "total elements: %d\n", num_2d + num_3d );
 printf( "number of leaves: %lu\n", (unsigned long)depth.count );
 printf( "number of nodes: %lu\n", (unsigned long)tree_sets.size() );
 printf( "volume ratio: %0.2f%%\n", 100 * ( vol.sum / tree_vol ) );
 printf( "surface ratio: %0.2f%%\n", 100 * ( surf.sum / tree_surf_area ) );
 printf( "\nmemory: used amortized\n" );
 printf( " ---------- ----------\n" );
 printf( "elements %10s %10s\n", mem_to_string( elem_used ).c_str(), mem_to_string( elem_amortized ).c_str() );
 printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() );
 printf( "sets %10s %10s\n", mem_to_string( set_store_used ).c_str(),
 mem_to_string( set_store_amortized ).c_str() );
 printf( "set tags %10s %10s\n", mem_to_string( set_tag_used ).c_str(),
 mem_to_string( set_tag_amortized ).c_str() );
 printf( "\nleaf stats: min avg rms max std.dev\n" );
 printf( " ---------- ---------- ---------- ---------- ----------\n" );
 printf( "depth %10u %10.1f %10.1f %10u %10.2f\n", depth.min, depth.avg(), depth.rms(), depth.max,
 depth.dev() );
 printf( "triangles %10u %10.1f %10.1f %10u %10.2f\n", size.min, size.avg(), size.rms(), size.max, size.dev() );
 printf( "volume %10.2g %10.2g %10.2g %10.2g %10.2g\n", vol.min, vol.avg(), vol.rms(), vol.max, vol.dev() );
 printf( "surf. area %10.2g %10.2g %10.2g %10.2g %10.2g\n", surf.min, surf.avg(), surf.rms(), surf.max,
 surf.dev() );
 printf( "------------------------------------------------------------------\n" );
 
 return MB_SUCCESS;
}
 
} // namespace moab