MOAB: Mesh Oriented datABase  (version 5.2.1)
OrientedBoxTreeTool.cpp
Go to the documentation of this file.
00001 /*
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 /**\file OrientedBox.hpp
00017  *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
00018  *\date 2006-07-18
00019  */
00020 
00021 #include "moab/Interface.hpp"
00022 #include "Internals.hpp"
00023 #include "moab/OrientedBoxTreeTool.hpp"
00024 #include "moab/Range.hpp"
00025 #include "moab/CN.hpp"
00026 #include "moab/GeomUtil.hpp"
00027 #include "MBTagConventions.hpp"
00028 #include <iostream>
00029 #include <iomanip>
00030 #include <algorithm>
00031 #include <limits>
00032 #include <assert.h>
00033 #include <math.h>
00034 
00035 //#define MB_OBB_USE_VECTOR_QUERIES
00036 //#define MB_OBB_USE_TYPE_QUERIES
00037 
00038 namespace moab
00039 {
00040 
00041 #if defined( MB_OBB_USE_VECTOR_QUERIES ) && defined( MB_OBB_USE_TYPE_QUERIES )
00042 #undef MB_OBB_USE_TYPE_QUERIES
00043 #endif
00044 
00045 const char DEFAULT_TAG_NAME[] = "OBB";
00046 
00047 OrientedBoxTreeTool::Op::~Op() {}
00048 
00049 OrientedBoxTreeTool::OrientedBoxTreeTool( Interface* i, const char* tag_name, bool destroy_created_trees )
00050     : instance( i ), cleanUpTrees( destroy_created_trees )
00051 {
00052     if( !tag_name ) tag_name = DEFAULT_TAG_NAME;
00053     ErrorCode rval = OrientedBox::tag_handle( tagHandle, instance, tag_name );
00054     if( MB_SUCCESS != rval ) tagHandle = 0;
00055 }
00056 
00057 OrientedBoxTreeTool::~OrientedBoxTreeTool()
00058 {
00059     if( !cleanUpTrees ) return;
00060 
00061     while( !createdTrees.empty() )
00062     {
00063         EntityHandle tree = createdTrees.back();
00064         // make sure this is a tree (rather than some other, stale handle)
00065         const void* data_ptr = 0;
00066         ErrorCode rval       = instance->tag_get_by_ptr( tagHandle, &tree, 1, &data_ptr );
00067         if( MB_SUCCESS == rval ) rval = delete_tree( tree );
00068         if( MB_SUCCESS != rval ) createdTrees.pop_back();
00069     }
00070 }
00071 
00072 OrientedBoxTreeTool::Settings::Settings()
00073     : max_leaf_entities( 8 ), max_depth( 0 ), worst_split_ratio( 0.95 ), best_split_ratio( 0.4 ),
00074       set_options( MESHSET_SET )
00075 {
00076 }
00077 
00078 bool OrientedBoxTreeTool::Settings::valid() const
00079 {
00080     return max_leaf_entities > 0 && max_depth >= 0 && worst_split_ratio <= 1.0 && best_split_ratio >= 0.0 &&
00081            worst_split_ratio >= best_split_ratio;
00082 }
00083 
00084 /********************** Simple Tree Access Methods ****************************/
00085 
00086 ErrorCode OrientedBoxTreeTool::box( EntityHandle set, OrientedBox& obb )
00087 {
00088     return instance->tag_get_data( tagHandle, &set, 1, &obb );
00089 }
00090 
00091 ErrorCode OrientedBoxTreeTool::box( EntityHandle set, double center[3], double axis1[3], double axis2[3],
00092                                     double axis3[3] )
00093 {
00094     OrientedBox obb;
00095     ErrorCode rval = this->box( set, obb );
00096     obb.center.get( center );
00097     obb.scaled_axis( 0 ).get( axis1 );
00098     obb.scaled_axis( 1 ).get( axis2 );
00099     obb.scaled_axis( 2 ).get( axis3 );
00100     return rval;
00101 }
00102 
00103 /********************** Tree Construction Code ****************************/
00104 
00105 struct OrientedBoxTreeTool::SetData
00106 {
00107     EntityHandle handle;
00108     OrientedBox::CovarienceData box_data;
00109     // Range vertices;
00110 };
00111 
00112 ErrorCode OrientedBoxTreeTool::build( const Range& entities, EntityHandle& set_handle_out, const Settings* settings )
00113 {
00114     if( !entities.all_of_dimension( 2 ) ) return MB_TYPE_OUT_OF_RANGE;
00115     if( settings && !settings->valid() ) return MB_FAILURE;
00116 
00117     return build_tree( entities, set_handle_out, 0, settings ? *settings : Settings() );
00118 }
00119 
00120 ErrorCode OrientedBoxTreeTool::join_trees( const Range& sets, EntityHandle& set_handle_out, const Settings* settings )
00121 {
00122     if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
00123     if( settings && !settings->valid() ) return MB_FAILURE;
00124 
00125     // Build initial set data list.
00126     std::list< SetData > data;
00127     for( Range::const_iterator i = sets.begin(); i != sets.end(); ++i )
00128     {
00129         Range elements;
00130         ErrorCode rval = instance->get_entities_by_dimension( *i, 2, elements, true );
00131         if( MB_SUCCESS != rval ) return rval;
00132         if( elements.empty() ) continue;
00133 
00134         data.push_back( SetData() );
00135         SetData& set_data = data.back();
00136         set_data.handle   = *i;
00137         rval              = OrientedBox::covariance_data_from_tris( set_data.box_data, instance, elements );
00138         if( MB_SUCCESS != rval ) return rval;
00139     }
00140 
00141     ErrorCode result = build_sets( data, set_handle_out, 0, settings ? *settings : Settings() );
00142     if( MB_SUCCESS != result ) return result;
00143 
00144     for( Range::reverse_iterator i = sets.rbegin(); i != sets.rend(); ++i )
00145     {
00146         createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), *i ), createdTrees.end() );
00147     }
00148     createdTrees.push_back( set_handle_out );
00149     return MB_SUCCESS;
00150 }
00151 
00152 /**\brief Split triangles by which side of a plane they are on
00153  *
00154  * Given a plane specified as a bisecting plane normal to one
00155  * of the axes of a box, split triangles based on which side
00156  * of the plane they are on.
00157  *\param instance   MOAB instance
00158  *\param box        The oriented box containing all the entities
00159  *\param axis       The axis for which the split plane is orthogonal
00160  *\param left_list  Output, entities to the left of the plane
00161  *\param right_list Output, entities to the right of the plane
00162  *\param num_intersecting Output, number entities intersecting plane
00163  */
00164 static ErrorCode split_box( Interface* instance, const OrientedBox& box, int axis, const Range& entities,
00165                             Range& left_list, Range& right_list )
00166 {
00167     ErrorCode rval;
00168     left_list.clear();
00169     right_list.clear();
00170 
00171     std::vector< CartVect > coords;
00172     for( Range::reverse_iterator i = entities.rbegin(); i != entities.rend(); ++i )
00173     {
00174         const EntityHandle* conn = NULL;
00175         int conn_len             = 0;
00176         rval                     = instance->get_connectivity( *i, conn, conn_len );
00177         if( MB_SUCCESS != rval ) return rval;
00178 
00179         coords.resize( conn_len );
00180         rval = instance->get_coords( conn, conn_len, coords[0].array() );
00181         if( MB_SUCCESS != rval ) return rval;
00182 
00183         CartVect centroid( 0.0 );
00184         for( int j = 0; j < conn_len; ++j )
00185             centroid += coords[j];
00186         centroid /= conn_len;
00187 
00188         if( ( box.axis( axis ) % ( centroid - box.center ) ) < 0.0 )
00189             left_list.insert( *i );
00190         else
00191             right_list.insert( *i );
00192     }
00193 
00194     return MB_SUCCESS;
00195 }
00196 
00197 ErrorCode OrientedBoxTreeTool::build_tree( const Range& entities, EntityHandle& set, int depth,
00198                                            const Settings& settings )
00199 {
00200     OrientedBox tmp_box;
00201     ErrorCode rval;
00202 
00203     if( entities.empty() )
00204     {
00205         Matrix3 axis;
00206         tmp_box = OrientedBox( axis, CartVect( 0. ) );
00207     }
00208     else
00209     {
00210         rval = OrientedBox::compute_from_2d_cells( tmp_box, instance, entities );
00211         if( MB_SUCCESS != rval ) return rval;
00212     }
00213 
00214     // create an entity set for the tree node
00215     rval = instance->create_meshset( settings.set_options, set );
00216     if( MB_SUCCESS != rval ) return rval;
00217 
00218     rval = instance->tag_set_data( tagHandle, &set, 1, &tmp_box );
00219     if( MB_SUCCESS != rval )
00220     {
00221         delete_tree( set );
00222         return rval;
00223     }
00224 
00225     // check if should create children
00226     bool leaf = true;
00227     ++depth;
00228     if( ( !settings.max_depth || depth < settings.max_depth ) &&
00229         entities.size() > (unsigned)settings.max_leaf_entities )
00230     {
00231         // try splitting with planes normal to each axis of the box
00232         // until we find an acceptable split
00233         double best_ratio = settings.worst_split_ratio;  // worst case ratio
00234         Range best_left_list, best_right_list;
00235         // Axes are sorted from shortest to longest, so search backwards
00236         for( int axis = 2; best_ratio > settings.best_split_ratio && axis >= 0; --axis )
00237         {
00238             Range left_list, right_list;
00239 
00240             rval = split_box( instance, tmp_box, axis, entities, left_list, right_list );
00241             if( MB_SUCCESS != rval )
00242             {
00243                 delete_tree( set );
00244                 return rval;
00245             }
00246 
00247             double ratio = fabs( (double)right_list.size() - left_list.size() ) / entities.size();
00248 
00249             if( ratio < best_ratio )
00250             {
00251                 best_ratio = ratio;
00252                 best_left_list.swap( left_list );
00253                 best_right_list.swap( right_list );
00254             }
00255         }
00256 
00257         // create children
00258         if( !best_left_list.empty() )
00259         {
00260             EntityHandle child = 0;
00261 
00262             rval = build_tree( best_left_list, child, depth, settings );
00263             if( MB_SUCCESS != rval )
00264             {
00265                 delete_tree( set );
00266                 return rval;
00267             }
00268             rval = instance->add_child_meshset( set, child );
00269             if( MB_SUCCESS != rval )
00270             {
00271                 delete_tree( set );
00272                 delete_tree( child );
00273                 return rval;
00274             }
00275 
00276             rval = build_tree( best_right_list, child, depth, settings );
00277             if( MB_SUCCESS != rval )
00278             {
00279                 delete_tree( set );
00280                 return rval;
00281             }
00282             rval = instance->add_child_meshset( set, child );
00283             if( MB_SUCCESS != rval )
00284             {
00285                 delete_tree( set );
00286                 delete_tree( child );
00287                 return rval;
00288             }
00289 
00290             leaf = false;
00291         }
00292     }
00293 
00294     if( leaf )
00295     {
00296         rval = instance->add_entities( set, entities );
00297         if( MB_SUCCESS != rval )
00298         {
00299             delete_tree( set );
00300             return rval;
00301         }
00302     }
00303 
00304     createdTrees.push_back( set );
00305     return MB_SUCCESS;
00306 }
00307 
00308 static ErrorCode split_sets( Interface*, const OrientedBox& box, int axis,
00309                              const std::list< OrientedBoxTreeTool::SetData >& sets,
00310                              std::list< OrientedBoxTreeTool::SetData >& left,
00311                              std::list< OrientedBoxTreeTool::SetData >& right )
00312 {
00313     left.clear();
00314     right.clear();
00315 
00316     std::list< OrientedBoxTreeTool::SetData >::const_iterator i;
00317     for( i = sets.begin(); i != sets.end(); ++i )
00318     {
00319         CartVect centroid( i->box_data.center / i->box_data.area );
00320         if( ( box.axis( axis ) % ( centroid - box.center ) ) < 0.0 )
00321             left.push_back( *i );
00322         else
00323             right.push_back( *i );
00324     }
00325 
00326     return MB_SUCCESS;
00327 }
00328 
00329 ErrorCode OrientedBoxTreeTool::build_sets( std::list< SetData >& sets, EntityHandle& node_set, int depth,
00330                                            const Settings& settings )
00331 {
00332     ErrorCode rval;
00333     int count = sets.size();
00334     if( 0 == count ) return MB_FAILURE;
00335 
00336     // calculate box
00337     OrientedBox obox;
00338 
00339     // make vector go out of scope when done, so memory is released
00340     {
00341         Range elems;
00342         std::vector< OrientedBox::CovarienceData > data( sets.size() );
00343         data.clear();
00344         for( std::list< SetData >::iterator i = sets.begin(); i != sets.end(); ++i )
00345         {
00346             data.push_back( i->box_data );
00347             rval = instance->get_entities_by_dimension( i->handle, 2, elems, true );
00348             if( MB_SUCCESS != rval ) return rval;
00349         }
00350 
00351         Range points;
00352         rval = instance->get_adjacencies( elems, 0, false, points, Interface::UNION );
00353         if( MB_SUCCESS != rval ) return rval;
00354 
00355         rval = OrientedBox::compute_from_covariance_data( obox, instance, &data[0], data.size(), points );
00356         if( MB_SUCCESS != rval ) return rval;
00357     }
00358 
00359     // If only one set in list...
00360     if( count == 1 )
00361     {
00362         node_set = sets.front().handle;
00363         return instance->tag_set_data( tagHandle, &node_set, 1, &obox );
00364     }
00365 
00366     // create an entity set for the tree node
00367     rval = instance->create_meshset( settings.set_options, node_set );
00368     if( MB_SUCCESS != rval ) return rval;
00369 
00370     rval = instance->tag_set_data( tagHandle, &node_set, 1, &obox );
00371     if( MB_SUCCESS != rval )
00372     {
00373         delete_tree( node_set );
00374         return rval;
00375     }
00376 
00377     double best_ratio = 2.0;
00378     std::list< SetData > best_left_list, best_right_list;
00379     for( int axis = 0; axis < 2; ++axis )
00380     {
00381         std::list< SetData > left_list, right_list;
00382         rval = split_sets( instance, obox, axis, sets, left_list, right_list );
00383         if( MB_SUCCESS != rval )
00384         {
00385             delete_tree( node_set );
00386             return rval;
00387         }
00388 
00389         double ratio = fabs( (double)right_list.size() - left_list.size() ) / sets.size();
00390 
00391         if( ratio < best_ratio )
00392         {
00393             best_ratio = ratio;
00394             best_left_list.swap( left_list );
00395             best_right_list.swap( right_list );
00396         }
00397     }
00398 
00399     // We must subdivide the list of sets, because we want to guarantee that
00400     // there is a node in the tree corresponding to each of the sets.  So if
00401     // we couldn't find a usable split plane, just split them in an arbitrary
00402     // manner.
00403     if( best_left_list.empty() || best_right_list.empty() )
00404     {
00405         best_left_list.clear();
00406         best_right_list.clear();
00407         std::list< SetData >* lists[2] = { &best_left_list, &best_right_list };
00408         int i                          = 0;
00409         while( !sets.empty() )
00410         {
00411             lists[i]->push_back( sets.front() );
00412             sets.pop_front();
00413             i = 1 - i;
00414         }
00415     }
00416     else
00417     {
00418         sets.clear();  // release memory before recursion
00419     }
00420 
00421     // Create child sets
00422 
00423     EntityHandle child = 0;
00424 
00425     rval = build_sets( best_left_list, child, depth + 1, settings );
00426     if( MB_SUCCESS != rval )
00427     {
00428         delete_tree( node_set );
00429         return rval;
00430     }
00431     rval = instance->add_child_meshset( node_set, child );
00432     if( MB_SUCCESS != rval )
00433     {
00434         delete_tree( node_set );
00435         delete_tree( child );
00436         return rval;
00437     }
00438 
00439     rval = build_sets( best_right_list, child, depth + 1, settings );
00440     if( MB_SUCCESS != rval )
00441     {
00442         delete_tree( node_set );
00443         return rval;
00444     }
00445     rval = instance->add_child_meshset( node_set, child );
00446     if( MB_SUCCESS != rval )
00447     {
00448         delete_tree( node_set );
00449         delete_tree( child );
00450         return rval;
00451     }
00452 
00453     return MB_SUCCESS;
00454 }
00455 
00456 ErrorCode OrientedBoxTreeTool::delete_tree( EntityHandle set )
00457 {
00458     std::vector< EntityHandle > children;
00459     ErrorCode rval = instance->get_child_meshsets( set, children, 0 );
00460     if( MB_SUCCESS != rval ) return rval;
00461 
00462     createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), set ), createdTrees.end() );
00463     children.insert( children.begin(), set );
00464     return instance->delete_entities( &children[0], children.size() );
00465 }
00466 
00467 ErrorCode OrientedBoxTreeTool::remove_root( EntityHandle root )
00468 {
00469     std::vector< EntityHandle >::iterator i = find( createdTrees.begin(), createdTrees.end(), root );
00470     if( i != createdTrees.end() )
00471     {
00472         createdTrees.erase( i );
00473         return MB_SUCCESS;
00474     }
00475     else
00476         return MB_ENTITY_NOT_FOUND;
00477 }
00478 
00479 /********************** Generic Tree Traversal ****************************/
00480 struct Data
00481 {
00482     EntityHandle set;
00483     int depth;
00484 };
00485 ErrorCode OrientedBoxTreeTool::preorder_traverse( EntityHandle set, Op& operation, TrvStats* accum )
00486 {
00487     ErrorCode rval;
00488     std::vector< EntityHandle > children;
00489     std::vector< Data > the_stack;
00490     Data data = { set, 0 };
00491     the_stack.push_back( data );
00492     int max_depth = -1;
00493 
00494     while( !the_stack.empty() )
00495     {
00496         data = the_stack.back();
00497         the_stack.pop_back();
00498 
00499         // increment traversal statistics
00500         if( accum )
00501         {
00502             accum->increment( data.depth );
00503             max_depth = std::max( max_depth, data.depth );
00504         }
00505 
00506         bool descend = true;
00507         rval         = operation.visit( data.set, data.depth, descend );
00508         assert( MB_SUCCESS == rval );
00509         if( MB_SUCCESS != rval ) return rval;
00510 
00511         if( !descend ) continue;
00512 
00513         children.clear();
00514         rval = instance->get_child_meshsets( data.set, children );
00515         assert( MB_SUCCESS == rval );
00516         if( MB_SUCCESS != rval ) return rval;
00517         if( children.empty() )
00518         {
00519             if( accum ) { accum->increment_leaf( data.depth ); }
00520             rval = operation.leaf( data.set );
00521             assert( MB_SUCCESS == rval );
00522             if( MB_SUCCESS != rval ) return rval;
00523         }
00524         else if( children.size() == 2 )
00525         {
00526             data.depth++;
00527             data.set = children[0];
00528             the_stack.push_back( data );
00529             data.set = children[1];
00530             the_stack.push_back( data );
00531         }
00532         else
00533             return MB_MULTIPLE_ENTITIES_FOUND;
00534     }
00535 
00536     if( accum ) { accum->end_traversal( max_depth ); }
00537 
00538     return MB_SUCCESS;
00539 }
00540 
00541 /********************** General Sphere/Triangle Intersection ***************/
00542 
00543 struct OBBTreeSITFrame
00544 {
00545     OBBTreeSITFrame( EntityHandle n, EntityHandle s, int dp ) : node( n ), surf( s ), depth( dp ) {}
00546     EntityHandle node;
00547     EntityHandle surf;
00548     int depth;
00549 };
00550 
00551 ErrorCode OrientedBoxTreeTool::sphere_intersect_triangles( const double* center_v, double radius,
00552                                                            EntityHandle tree_root,
00553                                                            std::vector< EntityHandle >& facets_out,
00554                                                            std::vector< EntityHandle >* sets_out, TrvStats* accum )
00555 {
00556     const double radsqr = radius * radius;
00557     OrientedBox b;
00558     ErrorCode rval;
00559     Range sets;
00560     const CartVect center( center_v );
00561     CartVect closest, coords[3];
00562     const EntityHandle* conn;
00563     int num_conn;
00564 #ifndef MB_OBB_USE_VECTOR_QUERIES
00565     Range tris;
00566     Range::const_iterator t;
00567 #else
00568     std::vector< EntityHandle > tris;
00569     std::vector< EntityHandle >::const_iterator t;
00570 #endif
00571 
00572     std::vector< OBBTreeSITFrame > stack;
00573     std::vector< EntityHandle > children;
00574     stack.reserve( 30 );
00575     stack.push_back( OBBTreeSITFrame( tree_root, 0, 0 ) );
00576 
00577     int max_depth = -1;
00578 
00579     while( !stack.empty() )
00580     {
00581         EntityHandle surf = stack.back().surf;
00582         EntityHandle node = stack.back().node;
00583         int current_depth = stack.back().depth;
00584         stack.pop_back();
00585 
00586         // increment traversal statistics.
00587         if( accum )
00588         {
00589             accum->increment( current_depth );
00590             max_depth = std::max( max_depth, current_depth );
00591         }
00592 
00593         if( !surf && sets_out )
00594         {
00595             rval = get_moab_instance()->get_entities_by_type( node, MBENTITYSET, sets );
00596             if( !sets.empty() ) surf = sets.front();
00597             sets.clear();
00598         }
00599 
00600         // check if sphere intersects box
00601         rval = box( node, b );
00602         if( MB_SUCCESS != rval ) return rval;
00603         b.closest_location_in_box( center, closest );
00604         closest -= center;
00605         if( ( closest % closest ) > radsqr ) continue;
00606 
00607         // push child boxes on stack
00608         children.clear();
00609         rval = instance->get_child_meshsets( node, children );
00610         if( MB_SUCCESS != rval ) return rval;
00611         if( !children.empty() )
00612         {
00613             assert( children.size() == 2 );
00614             stack.push_back( OBBTreeSITFrame( children[0], surf, current_depth + 1 ) );
00615             stack.push_back( OBBTreeSITFrame( children[1], surf, current_depth + 1 ) );
00616             continue;
00617         }
00618 
00619         if( accum ) { accum->increment_leaf( current_depth ); }
00620         // if leaf, intersect sphere with triangles
00621 #ifndef MB_OBB_USE_VECTOR_QUERIES
00622 #ifdef MB_OBB_USE_TYPE_QUERIES
00623         rval = get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
00624 #else
00625         rval = get_moab_instance()->get_entities_by_handle( node, tris );
00626 #endif
00627         t = tris.begin();
00628 #else
00629         rval = get_moab_instance()->get_entities_by_handle( node, tris );
00630         t    = tris.lower_bound( MBTRI );
00631 #endif
00632         if( MB_SUCCESS != rval ) return rval;
00633 
00634         for( t = tris.begin(); t != tris.end(); ++t )
00635         {
00636 #ifndef MB_OBB_USE_VECTOR_QUERIES
00637             if( TYPE_FROM_HANDLE( *t ) != MBTRI ) break;
00638 #elif !defined( MB_OBB_USE_TYPE_QUERIES )
00639             if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
00640 #endif
00641             rval = get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
00642             if( MB_SUCCESS != rval ) return rval;
00643             if( num_conn != 3 ) continue;
00644 
00645             rval = get_moab_instance()->get_coords( conn, num_conn, coords[0].array() );
00646             if( MB_SUCCESS != rval ) return rval;
00647 
00648             GeomUtil::closest_location_on_tri( center, coords, closest );
00649             closest -= center;
00650             if( ( closest % closest ) <= radsqr &&
00651                 std::find( facets_out.begin(), facets_out.end(), *t ) == facets_out.end() )
00652             {
00653                 facets_out.push_back( *t );
00654                 if( sets_out ) sets_out->push_back( surf );
00655             }
00656         }
00657     }
00658 
00659     if( accum ) { accum->end_traversal( max_depth ); }
00660 
00661     return MB_SUCCESS;
00662 }
00663 
00664 /********************** General Ray/Tree and Ray/Triangle Intersection ***************/
00665 
00666 class RayIntersector : public OrientedBoxTreeTool::Op
00667 {
00668   private:
00669     OrientedBoxTreeTool* tool;
00670     const CartVect b, m;
00671     const double* len;
00672     const double tol;
00673     Range& boxes;
00674 
00675   public:
00676     RayIntersector( OrientedBoxTreeTool* tool_ptr, const double* ray_point, const double* unit_ray_dir,
00677                     const double* ray_length, double tolerance, Range& leaf_boxes )
00678         : tool( tool_ptr ), b( ray_point ), m( unit_ray_dir ), len( ray_length ), tol( tolerance ), boxes( leaf_boxes )
00679     {
00680     }
00681 
00682     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
00683     virtual ErrorCode leaf( EntityHandle node );
00684 };
00685 
00686 //#include <stdio.h>
00687 // inline void dump_fragmentation( const Range& range ) {
00688 //  static FILE* file = fopen( "fragmentation", "w" );
00689 //  unsigned ranges = 0, entities = 0;
00690 //  for (Range::const_pair_iterator i = range.const_pair_begin(); i != range.const_pair_end(); ++i)
00691 //  {
00692 //    ++ranges;
00693 //    entities += i->second - i->first + 1;
00694 //  }
00695 //  fprintf( file, "%u %u\n", ranges, entities );
00696 //}
00697 
00698 ErrorCode OrientedBoxTreeTool::ray_intersect_triangles( std::vector< double >& intersection_distances_out,
00699                                                         std::vector< EntityHandle >& intersection_facets_out,
00700                                                         const Range& boxes, double /*tolerance*/,
00701                                                         const double ray_point[3], const double unit_ray_dir[3],
00702                                                         const double* ray_length, unsigned int* raytri_test_count )
00703 {
00704     ErrorCode rval;
00705     intersection_distances_out.clear();
00706 #ifdef MB_OBB_USE_VECTOR_QUERIES
00707     std::vector< EntityHandle > tris;
00708 #endif
00709 
00710     const CartVect point( ray_point );
00711     const CartVect dir( unit_ray_dir );
00712 
00713     for( Range::iterator b = boxes.begin(); b != boxes.end(); ++b )
00714     {
00715 #ifndef MB_OBB_USE_VECTOR_QUERIES
00716         Range tris;
00717 #ifdef MB_OBB_USE_TYPE_QUERIES
00718         rval = instance->get_entities_by_type( *b, MBTRI, tris );
00719 #else
00720         rval = instance->get_entities_by_handle( *b, tris );
00721 #endif
00722 #else
00723         tris.clear();
00724         rval = instance->get_entities_by_handle( *b, tris );
00725 #endif
00726         if( MB_SUCCESS != rval ) return rval;
00727             // dump_fragmentation( tris );
00728 
00729 #ifndef MB_OBB_USE_VECTOR_QUERIES
00730         for( Range::iterator t = tris.begin(); t != tris.end(); ++t )
00731 #else
00732         for( std::vector< EntityHandle >::iterator t = tris.begin(); t != tris.end(); ++t )
00733 #endif
00734         {
00735 #ifndef MB_OBB_USE_TYPE_QUERIES
00736             if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
00737 #endif
00738 
00739             const EntityHandle* conn = NULL;
00740             int len                  = 0;
00741             rval                     = instance->get_connectivity( *t, conn, len, true );
00742             if( MB_SUCCESS != rval ) return rval;
00743 
00744             CartVect coords[3];
00745             rval = instance->get_coords( conn, 3, coords[0].array() );
00746             if( MB_SUCCESS != rval ) return rval;
00747 
00748             if( raytri_test_count ) *raytri_test_count += 1;
00749 
00750             double td;
00751             if( GeomUtil::plucker_ray_tri_intersect( coords, point, dir, td, ray_length ) )
00752             {
00753                 intersection_distances_out.push_back( td );
00754                 intersection_facets_out.push_back( *t );
00755             }
00756         }
00757     }
00758 
00759     return MB_SUCCESS;
00760 }
00761 
00762 ErrorCode OrientedBoxTreeTool::ray_intersect_triangles( std::vector< double >& intersection_distances_out,
00763                                                         std::vector< EntityHandle >& intersection_facets_out,
00764                                                         EntityHandle root_set, double tolerance,
00765                                                         const double ray_point[3], const double unit_ray_dir[3],
00766                                                         const double* ray_length, TrvStats* accum )
00767 {
00768     Range boxes;
00769     ErrorCode rval;
00770 
00771     rval = ray_intersect_boxes( boxes, root_set, tolerance, ray_point, unit_ray_dir, ray_length, accum );
00772     if( MB_SUCCESS != rval ) return rval;
00773 
00774     return ray_intersect_triangles( intersection_distances_out, intersection_facets_out, boxes, tolerance, ray_point,
00775                                     unit_ray_dir, ray_length, accum ? &( accum->ray_tri_tests_count ) : NULL );
00776 }
00777 
00778 ErrorCode OrientedBoxTreeTool::ray_intersect_boxes( Range& boxes_out, EntityHandle root_set, double tolerance,
00779                                                     const double ray_point[3], const double unit_ray_dir[3],
00780                                                     const double* ray_length, TrvStats* accum )
00781 {
00782     RayIntersector op( this, ray_point, unit_ray_dir, ray_length, tolerance, boxes_out );
00783     return preorder_traverse( root_set, op, accum );
00784 }
00785 
00786 ErrorCode RayIntersector::visit( EntityHandle node, int, bool& descend )
00787 {
00788     OrientedBox box;
00789     ErrorCode rval = tool->box( node, box );
00790     if( MB_SUCCESS != rval ) return rval;
00791 
00792     descend = box.intersect_ray( b, m, tol, len );
00793     return MB_SUCCESS;
00794 }
00795 
00796 ErrorCode RayIntersector::leaf( EntityHandle node )
00797 {
00798     boxes.insert( node );
00799     return MB_SUCCESS;
00800 }
00801 
00802 /********************** Ray/Set Intersection ****************************/
00803 
00804 ErrorCode OrientedBoxTreeTool::get_close_tris( CartVect int_pt, double tol, const EntityHandle* rootSet,
00805                                                const EntityHandle* geomVol, const Tag* senseTag,
00806                                                std::vector< EntityHandle >& close_tris,
00807                                                std::vector< int >& close_senses )
00808 {
00809     std::vector< EntityHandle > close_surfs;
00810     ErrorCode rval = sphere_intersect_triangles( int_pt.array(), tol, *rootSet, close_tris, &close_surfs );
00811     assert( MB_SUCCESS == rval );
00812     if( MB_SUCCESS != rval ) return rval;
00813 
00814     // for each surface, get the surf sense wrt parent volume
00815     close_senses.resize( close_surfs.size() );
00816     for( unsigned i = 0; i < close_surfs.size(); ++i )
00817     {
00818         EntityHandle vols[2];
00819         rval = get_moab_instance()->tag_get_data( *senseTag, &( close_surfs[i] ), 1, vols );
00820         assert( MB_SUCCESS == rval );
00821         if( MB_SUCCESS != rval ) return rval;
00822         if( vols[0] == vols[1] )
00823         {
00824             std::cerr << "error: surf has positive and negative sense wrt same volume" << std::endl;
00825             return MB_FAILURE;
00826         }
00827         if( *geomVol == vols[0] ) { close_senses[i] = 1; }
00828         else if( *geomVol == vols[1] )
00829         {
00830             close_senses[i] = -1;
00831         }
00832         else
00833         {
00834             return MB_FAILURE;
00835         }
00836     }
00837 
00838     return MB_SUCCESS;
00839 }
00840 
00841 class RayIntersectSets : public OrientedBoxTreeTool::Op
00842 {
00843   private:
00844     // Input
00845     OrientedBoxTreeTool* tool;
00846     const CartVect ray_origin;
00847     const CartVect ray_direction;
00848     OrientedBoxTreeTool::IntersectSearchWindow& search_win; /* length to search ahead/behind of ray origin */
00849     const double tol;                                       /* used for box.intersect_ray, radius of
00850                                                                neighborhood for adjacent triangles,
00851                                                                and old mode of add_intersection */
00852     OrientedBoxTreeTool::IntRegCtxt& int_reg_callback;
00853 
00854     // Optional Input - to screen RTIs by orientation and edge/node intersection
00855     int* surfTriOrient; /* holds desired orientation of tri wrt surface */
00856     int surfTriOrient_val;
00857 
00858     // Other Variables
00859     unsigned int* raytri_test_count;
00860     EntityHandle lastSet;
00861     int lastSetDepth;
00862 
00863   public:
00864     RayIntersectSets( OrientedBoxTreeTool* tool_ptr, const double* ray_point, const double* unit_ray_dir,
00865                       const double tolerance, OrientedBoxTreeTool::IntersectSearchWindow& win,
00866                       unsigned int* ray_tri_test_count, OrientedBoxTreeTool::IntRegCtxt& intRegCallback )
00867         : tool( tool_ptr ), ray_origin( ray_point ), ray_direction( unit_ray_dir ), search_win( win ), tol( tolerance ),
00868           int_reg_callback( intRegCallback ), surfTriOrient_val( 0 ), raytri_test_count( ray_tri_test_count ),
00869           lastSet( 0 ), lastSetDepth( 0 )
00870     {
00871 
00872         // specified orientation should be 1 or -1, indicating ray and surface
00873         // normal in the same or opposite directions, respectively.
00874         if( int_reg_callback.getDesiredOrient() ) { surfTriOrient = &surfTriOrient_val; }
00875         else
00876         {
00877             surfTriOrient = NULL;
00878         }
00879 
00880         // check the limits
00881         if( search_win.first ) { assert( 0 <= *( search_win.first ) ); }
00882         if( search_win.second ) { assert( 0 >= *( search_win.second ) ); }
00883     };
00884 
00885     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
00886     virtual ErrorCode leaf( EntityHandle node );
00887 };
00888 
00889 ErrorCode RayIntersectSets::visit( EntityHandle node, int depth, bool& descend )
00890 {
00891     OrientedBox box;
00892     ErrorCode rval = tool->box( node, box );
00893     assert( MB_SUCCESS == rval );
00894     if( MB_SUCCESS != rval ) return rval;
00895 
00896     descend = box.intersect_ray( ray_origin, ray_direction, tol, search_win.first, search_win.second );
00897 
00898     if( lastSet && depth <= lastSetDepth ) lastSet = 0;
00899 
00900     if( descend && !lastSet )
00901     {
00902         Range tmp_sets;
00903         rval = tool->get_moab_instance()->get_entities_by_type( node, MBENTITYSET, tmp_sets );
00904         assert( MB_SUCCESS == rval );
00905         if( MB_SUCCESS != rval ) return rval;
00906 
00907         if( !tmp_sets.empty() )
00908         {
00909             if( tmp_sets.size() > 1 ) return MB_FAILURE;
00910             lastSet      = *tmp_sets.begin();
00911             lastSetDepth = depth;
00912 
00913             rval = int_reg_callback.update_orient( lastSet, surfTriOrient );
00914             if( MB_SUCCESS != rval ) return rval;
00915         }
00916     }
00917 
00918     return MB_SUCCESS;
00919 }
00920 
00921 ErrorCode RayIntersectSets::leaf( EntityHandle node )
00922 {
00923     assert( lastSet );
00924     if( !lastSet )  // if no surface has been visited yet, something's messed up.
00925         return MB_FAILURE;
00926 
00927 #ifndef MB_OBB_USE_VECTOR_QUERIES
00928     Range tris;
00929 #ifdef MB_OBB_USE_TYPE_QUERIES
00930     ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
00931 #else
00932     ErrorCode rval = tool->get_moab_instance()->get_entities_by_handle( node, tris );
00933 #endif
00934 #else
00935     std::vector< EntityHandle > tris;
00936     ErrorCode rval = tool->get_moab_instance()->get_entities_by_handle( node, tris );
00937 #endif
00938     assert( MB_SUCCESS == rval );
00939     if( MB_SUCCESS != rval ) return rval;
00940 
00941 #ifndef MB_OBB_USE_VECTOR_QUERIES
00942     for( Range::iterator t = tris.begin(); t != tris.end(); ++t )
00943 #else
00944     for( std::vector< EntityHandle >::iterator t = tris.begin(); t != tris.end(); ++t )
00945 #endif
00946     {
00947 #ifndef MB_OBB_USE_TYPE_QUERIES
00948         if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
00949 #endif
00950 
00951         const EntityHandle* conn;
00952         int num_conn;
00953         rval = tool->get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
00954         assert( MB_SUCCESS == rval );
00955         if( MB_SUCCESS != rval ) return rval;
00956 
00957         CartVect coords[3];
00958         rval = tool->get_moab_instance()->get_coords( conn, 3, coords[0].array() );
00959         assert( MB_SUCCESS == rval );
00960         if( MB_SUCCESS != rval ) return rval;
00961 
00962         if( raytri_test_count ) *raytri_test_count += 1;
00963 
00964         double int_dist;
00965         GeomUtil::intersection_type int_type = GeomUtil::NONE;
00966 
00967         if( GeomUtil::plucker_ray_tri_intersect( coords, ray_origin, ray_direction, int_dist, search_win.first,
00968                                                  search_win.second, surfTriOrient, &int_type ) )
00969         { int_reg_callback.register_intersection( lastSet, *t, int_dist, search_win, int_type ); }
00970     }
00971     return MB_SUCCESS;
00972 }
00973 
00974 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
00975                                                    std::vector< EntityHandle >& sets_out,
00976                                                    std::vector< EntityHandle >& facets_out, EntityHandle root_set,
00977                                                    const double tolerance, const double ray_point[3],
00978                                                    const double unit_ray_dir[3], IntersectSearchWindow& search_win,
00979                                                    IntRegCtxt& int_reg_callback, TrvStats* accum )
00980 
00981 {
00982     RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
00983                          accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
00984     ErrorCode rval = preorder_traverse( root_set, op, accum );
00985 
00986     distances_out = int_reg_callback.get_intersections();
00987     sets_out      = int_reg_callback.get_sets();
00988     facets_out    = int_reg_callback.get_facets();
00989 
00990     return rval;
00991 }
00992 
00993 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
00994                                                    std::vector< EntityHandle >& sets_out,
00995                                                    std::vector< EntityHandle >& facets_out, EntityHandle root_set,
00996                                                    double tolerance, const double ray_point[3],
00997                                                    const double unit_ray_dir[3], const double* ray_length,
00998                                                    TrvStats* accum )
00999 {
01000     IntRegCtxt int_reg_ctxt;
01001 
01002     OrientedBoxTreeTool::IntersectSearchWindow search_win( ray_length, (double*)0 );
01003 
01004     RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
01005                          accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_ctxt );
01006 
01007     ErrorCode rval = preorder_traverse( root_set, op, accum );
01008 
01009     if( MB_SUCCESS != rval ) { return rval; }
01010 
01011     distances_out = int_reg_ctxt.get_intersections();
01012     sets_out      = int_reg_ctxt.get_sets();
01013     facets_out    = int_reg_ctxt.get_facets();
01014 
01015     return MB_SUCCESS;
01016 }
01017 
01018 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( EntityHandle root_set, const double tolerance,
01019                                                    const double ray_point[3], const double unit_ray_dir[3],
01020                                                    IntersectSearchWindow& search_win, IntRegCtxt& int_reg_callback,
01021                                                    TrvStats* accum )
01022 
01023 {
01024     RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
01025                          accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
01026     return preorder_traverse( root_set, op, accum );
01027 }
01028 
01029 /********************** Closest Point code ***************/
01030 
01031 struct OBBTreeCPFrame
01032 {
01033     OBBTreeCPFrame( double d, EntityHandle n, EntityHandle s, int dp )
01034         : dist_sqr( d ), node( n ), mset( s ), depth( dp )
01035     {
01036     }
01037     double dist_sqr;
01038     EntityHandle node;
01039     EntityHandle mset;
01040     int depth;
01041 };
01042 
01043 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point, EntityHandle root, double* point_out,
01044                                                     EntityHandle& facet_out, EntityHandle* set_out, TrvStats* accum )
01045 {
01046     ErrorCode rval;
01047     const CartVect loc( point );
01048     double smallest_dist_sqr = std::numeric_limits< double >::max();
01049 
01050     EntityHandle current_set = 0;
01051     Range sets;
01052     std::vector< EntityHandle > children( 2 );
01053     std::vector< double > coords;
01054     std::vector< OBBTreeCPFrame > stack;
01055     int max_depth = -1;
01056 
01057     stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
01058 
01059     while( !stack.empty() )
01060     {
01061 
01062         // pop from top of stack
01063         EntityHandle node = stack.back().node;
01064         double dist_sqr   = stack.back().dist_sqr;
01065         current_set       = stack.back().mset;
01066         int current_depth = stack.back().depth;
01067         stack.pop_back();
01068 
01069         // If current best result is closer than the box, skip this tree node.
01070         if( dist_sqr > smallest_dist_sqr ) continue;
01071 
01072         // increment traversal statistics.
01073         if( accum )
01074         {
01075             accum->increment( current_depth );
01076             max_depth = std::max( max_depth, current_depth );
01077         }
01078 
01079         // Check if this node has a set associated with it
01080         if( set_out && !current_set )
01081         {
01082             sets.clear();
01083             rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
01084             if( MB_SUCCESS != rval ) return rval;
01085             if( !sets.empty() )
01086             {
01087                 if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
01088                 current_set = sets.front();
01089             }
01090         }
01091 
01092         // Get child boxes
01093         children.clear();
01094         rval = instance->get_child_meshsets( node, children );
01095         if( MB_SUCCESS != rval ) return rval;
01096 
01097         // if not a leaf node
01098         if( !children.empty() )
01099         {
01100             if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
01101 
01102             // get boxes
01103             OrientedBox box1, box2;
01104             rval = box( children[0], box1 );
01105             if( MB_SUCCESS != rval ) return rval;
01106             rval = box( children[1], box2 );
01107             if( MB_SUCCESS != rval ) return rval;
01108 
01109             // get distance from each box
01110             CartVect pt1, pt2;
01111             box1.closest_location_in_box( loc, pt1 );
01112             box2.closest_location_in_box( loc, pt2 );
01113             pt1 -= loc;
01114             pt2 -= loc;
01115             const double dsqr1 = pt1 % pt1;
01116             const double dsqr2 = pt2 % pt2;
01117 
01118             // push children on tree such that closer one is on top
01119             if( dsqr1 < dsqr2 )
01120             {
01121                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01122                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01123             }
01124             else
01125             {
01126                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01127                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01128             }
01129         }
01130         else
01131         {  // LEAF NODE
01132             if( accum ) { accum->increment_leaf( current_depth ); }
01133 
01134             Range facets;
01135             rval = instance->get_entities_by_dimension( node, 2, facets );
01136             if( MB_SUCCESS != rval ) return rval;
01137 
01138             const EntityHandle* conn = NULL;
01139             int len                  = 0;
01140             CartVect tmp, diff;
01141             for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
01142             {
01143                 rval = instance->get_connectivity( *i, conn, len, true );
01144                 if( MB_SUCCESS != rval ) return rval;
01145 
01146                 coords.resize( 3 * len );
01147                 rval = instance->get_coords( conn, len, &coords[0] );
01148                 if( MB_SUCCESS != rval ) return rval;
01149 
01150                 if( len == 3 )
01151                     GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
01152                 else
01153                     GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
01154 
01155                 diff     = tmp - loc;
01156                 dist_sqr = diff % diff;
01157                 if( dist_sqr < smallest_dist_sqr )
01158                 {
01159                     smallest_dist_sqr = dist_sqr;
01160                     facet_out         = *i;
01161                     tmp.get( point_out );
01162                     if( set_out ) *set_out = current_set;
01163                 }
01164             }
01165         }  // LEAF NODE
01166     }
01167 
01168     if( accum ) { accum->end_traversal( max_depth ); }
01169 
01170     return MB_SUCCESS;
01171 }
01172 
01173 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point, EntityHandle root, double tolerance,
01174                                                     std::vector< EntityHandle >& facets_out,
01175                                                     std::vector< EntityHandle >* sets_out, TrvStats* accum )
01176 {
01177     ErrorCode rval;
01178     const CartVect loc( point );
01179     double smallest_dist_sqr = std::numeric_limits< double >::max();
01180     double smallest_dist     = smallest_dist_sqr;
01181 
01182     EntityHandle current_set = 0;
01183     Range sets;
01184     std::vector< EntityHandle > children( 2 );
01185     std::vector< double > coords;
01186     std::vector< OBBTreeCPFrame > stack;
01187     int max_depth = -1;
01188 
01189     stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
01190 
01191     while( !stack.empty() )
01192     {
01193 
01194         // pop from top of stack
01195         EntityHandle node = stack.back().node;
01196         double dist_sqr   = stack.back().dist_sqr;
01197         current_set       = stack.back().mset;
01198         int current_depth = stack.back().depth;
01199         stack.pop_back();
01200 
01201         // If current best result is closer than the box, skip this tree node.
01202         if( dist_sqr > smallest_dist_sqr + tolerance ) continue;
01203 
01204         // increment traversal statistics.
01205         if( accum )
01206         {
01207             accum->increment( current_depth );
01208             max_depth = std::max( max_depth, current_depth );
01209         }
01210 
01211         // Check if this node has a set associated with it
01212         if( sets_out && !current_set )
01213         {
01214             sets.clear();
01215             rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
01216             if( MB_SUCCESS != rval ) return rval;
01217             if( !sets.empty() )
01218             {
01219                 if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
01220                 current_set = *sets.begin();
01221             }
01222         }
01223 
01224         // Get child boxes
01225         children.clear();
01226         rval = instance->get_child_meshsets( node, children );
01227         if( MB_SUCCESS != rval ) return rval;
01228 
01229         // if not a leaf node
01230         if( !children.empty() )
01231         {
01232             if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
01233 
01234             // get boxes
01235             OrientedBox box1, box2;
01236             rval = box( children[0], box1 );
01237             if( MB_SUCCESS != rval ) return rval;
01238             rval = box( children[1], box2 );
01239             if( MB_SUCCESS != rval ) return rval;
01240 
01241             // get distance from each box
01242             CartVect pt1, pt2;
01243             box1.closest_location_in_box( loc, pt1 );
01244             box2.closest_location_in_box( loc, pt2 );
01245             pt1 -= loc;
01246             pt2 -= loc;
01247             const double dsqr1 = pt1 % pt1;
01248             const double dsqr2 = pt2 % pt2;
01249 
01250             // push children on tree such that closer one is on top
01251             if( dsqr1 < dsqr2 )
01252             {
01253                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01254                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01255             }
01256             else
01257             {
01258                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01259                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01260             }
01261         }
01262         else
01263         {  // LEAF NODE
01264             if( accum ) { accum->increment_leaf( current_depth ); }
01265 
01266             Range facets;
01267             rval = instance->get_entities_by_dimension( node, 2, facets );
01268             if( MB_SUCCESS != rval ) return rval;
01269 
01270             const EntityHandle* conn = NULL;
01271             int len                  = 0;
01272             CartVect tmp, diff;
01273             for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
01274             {
01275                 rval = instance->get_connectivity( *i, conn, len, true );
01276                 if( MB_SUCCESS != rval ) return rval;
01277 
01278                 coords.resize( 3 * len );
01279                 rval = instance->get_coords( conn, len, &coords[0] );
01280                 if( MB_SUCCESS != rval ) return rval;
01281 
01282                 if( len == 3 )
01283                     GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
01284                 else
01285                     GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
01286 
01287                 diff     = tmp - loc;
01288                 dist_sqr = diff % diff;
01289                 if( dist_sqr < smallest_dist_sqr )
01290                 {
01291                     if( 0.5 * dist_sqr < 0.5 * smallest_dist_sqr + tolerance * ( 0.5 * tolerance - smallest_dist ) )
01292                     {
01293                         facets_out.clear();
01294                         if( sets_out ) sets_out->clear();
01295                     }
01296                     smallest_dist_sqr = dist_sqr;
01297                     smallest_dist     = sqrt( smallest_dist_sqr );
01298                     // put closest result at start of lists
01299                     facets_out.push_back( *i );
01300                     std::swap( facets_out.front(), facets_out.back() );
01301                     if( sets_out )
01302                     {
01303                         sets_out->push_back( current_set );
01304                         std::swap( sets_out->front(), sets_out->back() );
01305                     }
01306                 }
01307                 else if( dist_sqr <= smallest_dist_sqr + tolerance * ( tolerance + 2 * smallest_dist ) )
01308                 {
01309                     facets_out.push_back( *i );
01310                     if( sets_out ) sets_out->push_back( current_set );
01311                 }
01312             }
01313         }  // LEAF NODE
01314     }
01315 
01316     if( accum ) { accum->end_traversal( max_depth ); }
01317 
01318     return MB_SUCCESS;
01319 }
01320 
01321 /********************** Tree Printing Code ****************************/
01322 
01323 class TreeLayoutPrinter : public OrientedBoxTreeTool::Op
01324 {
01325   public:
01326     TreeLayoutPrinter( std::ostream& stream, Interface* instance );
01327 
01328     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
01329     virtual ErrorCode leaf( EntityHandle node );
01330 
01331   private:
01332     Interface* instance;
01333     std::ostream& outputStream;
01334     std::vector< bool > path;
01335 };
01336 
01337 TreeLayoutPrinter::TreeLayoutPrinter( std::ostream& stream, Interface* interface )
01338     : instance( interface ), outputStream( stream )
01339 {
01340 }
01341 
01342 ErrorCode TreeLayoutPrinter::visit( EntityHandle node, int depth, bool& descend )
01343 {
01344     descend = true;
01345 
01346     if( (unsigned)depth > path.size() )
01347     {
01348         // assert(depth+1 == path.size); // preorder traversal
01349         path.push_back( true );
01350     }
01351     else
01352     {
01353         path.resize( depth );
01354         if( depth ) path.back() = false;
01355     }
01356 
01357     for( unsigned i = 0; i + 1 < path.size(); ++i )
01358     {
01359         if( path[i] )
01360             outputStream << "|   ";
01361         else
01362             outputStream << "    ";
01363     }
01364     if( depth )
01365     {
01366         if( path.back() )
01367             outputStream << "+---";
01368         else
01369             outputStream << "\\---";
01370     }
01371     outputStream << instance->id_from_handle( node ) << std::endl;
01372     return MB_SUCCESS;
01373 }
01374 
01375 ErrorCode TreeLayoutPrinter::leaf( EntityHandle )
01376 {
01377     return MB_SUCCESS;
01378 }
01379 
01380 class TreeNodePrinter : public OrientedBoxTreeTool::Op
01381 {
01382   public:
01383     TreeNodePrinter( std::ostream& stream, bool list_contents, bool list_box, const char* id_tag_name,
01384                      OrientedBoxTreeTool* tool_ptr );
01385 
01386     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
01387 
01388     virtual ErrorCode leaf( EntityHandle )
01389     {
01390         return MB_SUCCESS;
01391     }
01392 
01393   private:
01394     ErrorCode print_geometry( EntityHandle node );
01395     ErrorCode print_contents( EntityHandle node );
01396     ErrorCode print_counts( EntityHandle node );
01397 
01398     bool printContents;
01399     bool printGeometry;
01400     bool haveTag;
01401     Tag tag, gidTag, geomTag;
01402     Interface* instance;
01403     OrientedBoxTreeTool* tool;
01404     std::ostream& outputStream;
01405 };
01406 
01407 TreeNodePrinter::TreeNodePrinter( std::ostream& stream, bool list_contents, bool list_box, const char* id_tag_name,
01408                                   OrientedBoxTreeTool* tool_ptr )
01409     : printContents( list_contents ), printGeometry( list_box ), haveTag( false ), tag( 0 ), gidTag( 0 ), geomTag( 0 ),
01410       instance( tool_ptr->get_moab_instance() ), tool( tool_ptr ), outputStream( stream )
01411 {
01412     ErrorCode rval;
01413     if( id_tag_name )
01414     {
01415         rval = instance->tag_get_handle( id_tag_name, 1, MB_TYPE_INTEGER, tag );
01416         if( !rval )
01417         {
01418             std::cerr << "Could not get tag \"" << id_tag_name << "\"\n";
01419             stream << "Could not get tag \"" << id_tag_name << "\"\n";
01420         }
01421         else
01422         {
01423             haveTag = true;
01424         }
01425     }
01426 
01427     gidTag = instance->globalId_tag();
01428 
01429     rval = instance->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );
01430     if( MB_SUCCESS != rval ) geomTag = 0;
01431 }
01432 
01433 ErrorCode TreeNodePrinter::visit( EntityHandle node, int, bool& descend )
01434 {
01435     descend            = true;
01436     EntityHandle setid = instance->id_from_handle( node );
01437     outputStream << setid << ":" << std::endl;
01438 
01439     Range surfs;
01440     ErrorCode r3 = MB_SUCCESS;
01441     if( geomTag )
01442     {
01443         const int two         = 2;
01444         const void* tagdata[] = { &two };
01445         r3 = instance->get_entities_by_type_and_tag( node, MBENTITYSET, &geomTag, tagdata, 1, surfs );
01446 
01447         if( MB_SUCCESS == r3 && surfs.size() == 1 )
01448         {
01449             EntityHandle surf = *surfs.begin();
01450             int id;
01451             if( gidTag && MB_SUCCESS == instance->tag_get_data( gidTag, &surf, 1, &id ) )
01452                 outputStream << "  Surface " << id << std::endl;
01453             else
01454                 outputStream << "  Surface w/ unknown ID (" << surf << ")" << std::endl;
01455         }
01456     }
01457 
01458     ErrorCode r1 = printGeometry ? print_geometry( node ) : MB_SUCCESS;
01459     ErrorCode r2 = printContents ? print_contents( node ) : print_counts( node );
01460     outputStream << std::endl;
01461 
01462     if( MB_SUCCESS != r1 )
01463         return r1;
01464     else if( MB_SUCCESS != r2 )
01465         return r2;
01466     else
01467         return r3;
01468 }
01469 
01470 ErrorCode TreeNodePrinter::print_geometry( EntityHandle node )
01471 {
01472     OrientedBox box;
01473     ErrorCode rval = tool->box( node, box );
01474     if( MB_SUCCESS != rval ) return rval;
01475 
01476     CartVect length = box.dimensions();
01477 
01478     outputStream << box.center << "  Radius: " << box.inner_radius() << " - " << box.outer_radius() << std::endl
01479                  << '+' << box.axis( 0 ) << " : " << length[0] << std::endl
01480                  << 'x' << box.axis( 1 ) << " : " << length[1] << std::endl
01481                  << 'x' << box.axis( 2 ) << " : " << length[2] << std::endl;
01482     return MB_SUCCESS;
01483 }
01484 
01485 ErrorCode TreeNodePrinter::print_counts( EntityHandle node )
01486 {
01487     for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
01488     {
01489         int count      = 0;
01490         ErrorCode rval = instance->get_number_entities_by_type( node, type, count );
01491         if( MB_SUCCESS != rval ) return rval;
01492         if( count > 0 ) outputStream << " " << count << " " << CN::EntityTypeName( type ) << std::endl;
01493     }
01494     return MB_SUCCESS;
01495 }
01496 
01497 ErrorCode TreeNodePrinter::print_contents( EntityHandle node )
01498 {
01499     // list contents
01500     for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
01501     {
01502         Range range;
01503         ErrorCode rval = instance->get_entities_by_type( node, type, range );
01504         if( MB_SUCCESS != rval ) return rval;
01505         if( range.empty() ) continue;
01506         outputStream << " " << CN::EntityTypeName( type ) << " ";
01507         std::vector< int > ids( range.size() );
01508         if( haveTag )
01509         {
01510             rval = instance->tag_get_data( tag, range, &ids[0] );
01511             std::sort( ids.begin(), ids.end() );
01512         }
01513         else
01514         {
01515             Range::iterator ri              = range.begin();
01516             std::vector< int >::iterator vi = ids.begin();
01517             while( ri != range.end() )
01518             {
01519                 *vi = instance->id_from_handle( *ri );
01520                 ++ri;
01521                 ++vi;
01522             }
01523         }
01524 
01525         unsigned i = 0;
01526         for( ;; )
01527         {
01528             unsigned beg = i, end;
01529             do
01530             {
01531                 end = i++;
01532             } while( i < ids.size() && ids[end] + 1 == ids[i] );
01533             if( end == beg )
01534                 outputStream << ids[end];
01535             else if( end == beg + 1 )
01536                 outputStream << ids[beg] << ", " << ids[end];
01537             else
01538                 outputStream << ids[beg] << "-" << ids[end];
01539 
01540             if( i == ids.size() )
01541             {
01542                 outputStream << std::endl;
01543                 break;
01544             }
01545             else
01546                 outputStream << ", ";
01547         }
01548     }
01549 
01550     return MB_SUCCESS;
01551 }
01552 
01553 void OrientedBoxTreeTool::print( EntityHandle set, std::ostream& str, bool list, const char* tag )
01554 {
01555     TreeLayoutPrinter op1( str, instance );
01556     TreeNodePrinter op2( str, list, true, tag, this );
01557     ErrorCode r1 = preorder_traverse( set, op1 );
01558     str << std::endl;
01559     ErrorCode r2 = preorder_traverse( set, op2 );
01560     if( r1 != MB_SUCCESS || r2 != MB_SUCCESS )
01561     {
01562         std::cerr << "Errors encountered while printing tree\n";
01563         str << "Errors encountered while printing tree\n";
01564     }
01565 }
01566 
01567 /********************* Traversal Metrics Code  **************************/
01568 
01569 void OrientedBoxTreeTool::TrvStats::reset()
01570 {
01571     nodes_visited_count.clear();
01572     leaves_visited_count.clear();
01573     traversals_ended_count.clear();
01574     ray_tri_tests_count = 0;
01575 }
01576 
01577 void OrientedBoxTreeTool::TrvStats::increment( unsigned depth )
01578 {
01579 
01580     while( nodes_visited_count.size() <= depth )
01581     {
01582         nodes_visited_count.push_back( 0 );
01583         leaves_visited_count.push_back( 0 );
01584         traversals_ended_count.push_back( 0 );
01585     }
01586     nodes_visited_count[depth] += 1;
01587 }
01588 
01589 void OrientedBoxTreeTool::TrvStats::increment_leaf( unsigned depth )
01590 {
01591     // assume safe depth, because increment is called first
01592     leaves_visited_count[depth] += 1;
01593 }
01594 
01595 void OrientedBoxTreeTool::TrvStats::end_traversal( unsigned depth )
01596 {
01597     // assume safe depth, because increment is always called on a given
01598     // tree level first
01599     traversals_ended_count[depth] += 1;
01600 }
01601 
01602 void OrientedBoxTreeTool::TrvStats::print( std::ostream& str ) const
01603 {
01604 
01605     const std::string h1 = "OBBTree Depth";
01606     const std::string h2 = " - NodesVisited";
01607     const std::string h3 = " - LeavesVisited";
01608     const std::string h4 = " - TraversalsEnded";
01609 
01610     str << h1 << h2 << h3 << h4 << std::endl;
01611 
01612     unsigned num_visited = 0, num_leaves = 0, num_traversals = 0;
01613     for( unsigned i = 0; i < traversals_ended_count.size(); ++i )
01614     {
01615 
01616         num_visited += nodes_visited_count[i];
01617         num_leaves += leaves_visited_count[i];
01618         num_traversals += traversals_ended_count[i];
01619 
01620         str << std::setw( h1.length() ) << i << std::setw( h2.length() ) << nodes_visited_count[i]
01621             << std::setw( h3.length() ) << leaves_visited_count[i] << std::setw( h4.length() )
01622             << traversals_ended_count[i] << std::endl;
01623     }
01624 
01625     str << std::setw( h1.length() ) << "---- Totals:" << std::setw( h2.length() ) << num_visited
01626         << std::setw( h3.length() ) << num_leaves << std::setw( h4.length() ) << num_traversals << std::endl;
01627 
01628     if( ray_tri_tests_count )
01629     { str << std::setw( h1.length() ) << "---- Total ray-tri tests: " << ray_tri_tests_count << std::endl; }
01630 }
01631 
01632 /********************** Tree Statistics Code ****************************/
01633 
01634 class StatData
01635 {
01636   public:
01637     struct Ratio
01638     {
01639         double min, max, sum, sqr;
01640         int hist[10];
01641         Ratio()
01642             : min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() ), sum( 0.0 ),
01643               sqr( 0.0 )
01644         {
01645             hist[0] = hist[1] = hist[2] = hist[3] = hist[4] = hist[5] = hist[6] = hist[7] = hist[8] = hist[9] = 0;
01646         }
01647         void accum( double v )
01648         {
01649             if( v < min ) min = v;
01650             if( v > max ) max = v;
01651             sum += v;
01652             sqr += v * v;
01653             int i = (int)( 10 * v );
01654             if( i < 0 )
01655                 i = 0;
01656             else if( i > 9 )
01657                 i = 9;
01658             ++hist[i];
01659         }
01660     };
01661 
01662     template < typename T >
01663     struct Stat
01664     {
01665         T min, max;
01666         double sum, sqr;
01667         Stat() : sum( 0.0 ), sqr( 0.0 )
01668         {
01669             std::numeric_limits< T > lim;
01670             min = lim.max();
01671             if( lim.is_integer )
01672                 max = lim.min();
01673             else
01674                 max = -lim.max();
01675         }
01676         void accum( T v )
01677         {
01678             if( v < min ) min = v;
01679             if( v > max ) max = v;
01680             sum += v;
01681             sqr += (double)v * v;
01682         }
01683     };
01684 
01685     StatData() : count( 0 ) {}
01686 
01687     Ratio volume;
01688     Ratio entities;
01689     Ratio radius;
01690     Stat< unsigned > leaf_ent;
01691     Stat< double > vol;
01692     Stat< double > area;
01693     std::vector< unsigned > leaf_depth;
01694     unsigned count;
01695 };
01696 
01697 static int measure( const CartVect& v, double& result )
01698 {
01699     const double tol = 1e-6;
01700     int dims         = 0;
01701     result           = 1;
01702     for( int i = 0; i < 3; ++i )
01703         if( v[i] > tol )
01704         {
01705             ++dims;
01706             result *= v[i];
01707         }
01708     return dims;
01709 }
01710 
01711 ErrorCode OrientedBoxTreeTool::recursive_stats( OrientedBoxTreeTool* tool, Interface* inst, EntityHandle set, int depth,
01712                                                 StatData& data, unsigned& count_out, CartVect& dimensions_out )
01713 {
01714     ErrorCode rval;
01715     OrientedBox tmp_box;
01716     std::vector< EntityHandle > children( 2 );
01717     unsigned counts[2];
01718     bool isleaf;
01719 
01720     ++data.count;
01721 
01722     rval = tool->box( set, tmp_box );
01723     if( MB_SUCCESS != rval ) return rval;
01724     children.clear();
01725     rval = inst->get_child_meshsets( set, children );
01726     if( MB_SUCCESS != rval ) return rval;
01727     isleaf = children.empty();
01728     if( !isleaf && children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
01729 
01730     dimensions_out = tmp_box.dimensions();
01731     data.radius.accum( tmp_box.inner_radius() / tmp_box.outer_radius() );
01732     data.vol.accum( tmp_box.volume() );
01733     data.area.accum( tmp_box.area() );
01734 
01735     if( isleaf )
01736     {
01737         if( data.leaf_depth.size() <= (unsigned)depth ) data.leaf_depth.resize( depth + 1, 0 );
01738         ++data.leaf_depth[depth];
01739 
01740         int count;
01741         rval = inst->get_number_entities_by_handle( set, count );
01742         if( MB_SUCCESS != rval ) return rval;
01743         count_out = count;
01744         data.leaf_ent.accum( count_out );
01745     }
01746     else
01747     {
01748         for( int i = 0; i < 2; ++i )
01749         {
01750             CartVect dims;
01751             rval = recursive_stats( tool, inst, children[i], depth + 1, data, counts[i], dims );
01752             if( MB_SUCCESS != rval ) return rval;
01753             double this_measure, chld_measure;
01754             int this_dim = measure( dimensions_out, this_measure );
01755             int chld_dim = measure( dims, chld_measure );
01756             double ratio;
01757             if( chld_dim < this_dim )
01758                 ratio = 0;
01759             else
01760                 ratio = chld_measure / this_measure;
01761 
01762             data.volume.accum( ratio );
01763         }
01764         count_out = counts[0] + counts[1];
01765         data.entities.accum( (double)counts[0] / count_out );
01766         data.entities.accum( (double)counts[1] / count_out );
01767     }
01768     return MB_SUCCESS;
01769 }
01770 
01771 static inline double std_dev( double sqr, double sum, double count )
01772 {
01773     sum /= count;
01774     sqr /= count;
01775     return sqrt( sqr - sum * sum );
01776 }
01777 
01778 //#define WW <<std::setw(10)<<std::fixed<<
01779 #define WE << std::setw( 10 ) <<
01780 #define WW WE
01781 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set, unsigned& total_entities, double& rv, double& tot_node_volume,
01782                                       double& tot_to_root_volume, unsigned& tree_height, unsigned& node_count,
01783                                       unsigned& num_leaves )
01784 {
01785     StatData d;
01786     ErrorCode rval;
01787     unsigned i;
01788     CartVect total_dim;
01789 
01790     rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
01791     if( MB_SUCCESS != rval ) return rval;
01792 
01793     tree_height                 = d.leaf_depth.size();
01794     unsigned min_leaf_depth     = tree_height;
01795     num_leaves                  = 0;
01796     unsigned max_leaf_per_depth = 0;
01797     // double sum_leaf_depth = 0, sqr_leaf_depth = 0;
01798     for( i = 0; i < d.leaf_depth.size(); ++i )
01799     {
01800         unsigned val = d.leaf_depth[i];
01801         num_leaves += val;
01802         // sum_leaf_depth += (double)val*i;
01803         // sqr_leaf_depth += (double)val*i*i;
01804         if( val && i < min_leaf_depth ) min_leaf_depth = i;
01805         if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
01806     }
01807     rv                 = total_dim[0] * total_dim[1] * total_dim[2];
01808     tot_node_volume    = d.vol.sum;
01809     tot_to_root_volume = d.vol.sum / rv;
01810     node_count         = d.count;
01811 
01812     return MB_SUCCESS;
01813 }
01814 
01815 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set, std::ostream& s )
01816 {
01817     StatData d;
01818     ErrorCode rval;
01819     unsigned total_entities, i;
01820     CartVect total_dim;
01821 
01822     rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
01823     if( MB_SUCCESS != rval ) return rval;
01824 
01825     unsigned tree_height    = d.leaf_depth.size();
01826     unsigned min_leaf_depth = tree_height, num_leaves = 0;
01827     unsigned max_leaf_per_depth = 0;
01828     double sum_leaf_depth = 0, sqr_leaf_depth = 0;
01829     for( i = 0; i < d.leaf_depth.size(); ++i )
01830     {
01831         unsigned val = d.leaf_depth[i];
01832         num_leaves += val;
01833         sum_leaf_depth += (double)val * i;
01834         sqr_leaf_depth += (double)val * i * i;
01835         if( val && i < min_leaf_depth ) min_leaf_depth = i;
01836         if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
01837     }
01838     unsigned num_non_leaf = d.count - num_leaves;
01839 
01840     double rv = total_dim[0] * total_dim[1] * total_dim[2];
01841     s << "entities in tree:  " << total_entities << std::endl
01842       << "root volume:       " << rv << std::endl
01843       << "total node volume: " << d.vol.sum << std::endl
01844       << "total/root volume: " << d.vol.sum / rv << std::endl
01845       << "tree height:       " << tree_height << std::endl
01846       << "node count:        " << d.count << std::endl
01847       << "leaf count:        " << num_leaves << std::endl
01848       << std::endl;
01849 
01850     double avg_leaf_depth = sum_leaf_depth / num_leaves;
01851     double rms_leaf_depth = sqrt( sqr_leaf_depth / num_leaves );
01852     double std_leaf_depth = std_dev( sqr_leaf_depth, sum_leaf_depth, num_leaves );
01853 
01854     double avg_leaf_ent = d.leaf_ent.sum / num_leaves;
01855     double rms_leaf_ent = sqrt( d.leaf_ent.sqr / num_leaves );
01856     double std_leaf_ent = std_dev( d.leaf_ent.sqr, d.leaf_ent.sum, num_leaves );
01857 
01858     unsigned num_child = 2 * num_non_leaf;
01859 
01860     double avg_vol_ratio = d.volume.sum / num_child;
01861     double rms_vol_ratio = sqrt( d.volume.sqr / num_child );
01862     double std_vol_ratio = std_dev( d.volume.sqr, d.volume.sum, num_child );
01863 
01864     double avg_ent_ratio = d.entities.sum / num_child;
01865     double rms_ent_ratio = sqrt( d.entities.sqr / num_child );
01866     double std_ent_ratio = std_dev( d.entities.sqr, d.entities.sum, num_child );
01867 
01868     double avg_rad_ratio = d.radius.sum / d.count;
01869     double rms_rad_ratio = sqrt( d.radius.sqr / d.count );
01870     double std_rad_ratio = std_dev( d.radius.sqr, d.radius.sum, d.count );
01871 
01872     double avg_vol = d.vol.sum / d.count;
01873     double rms_vol = sqrt( d.vol.sqr / d.count );
01874     double std_vol = std_dev( d.vol.sqr, d.vol.sum, d.count );
01875 
01876     double avg_area = d.area.sum / d.count;
01877     double rms_area = sqrt( d.area.sqr / d.count );
01878     double std_area = std_dev( d.area.sqr, d.area.sum, d.count );
01879 
01880     int prec = s.precision();
01881     s << "                   " WW "Minimum" WW "Average" WW "RMS" WW "Maximum" WW "Std.Dev." << std::endl;
01882     s << std::setprecision( 1 )
01883       << "Leaf Depth         " WW min_leaf_depth WW avg_leaf_depth WW rms_leaf_depth WW d.leaf_depth.size() -
01884              1 WW std_leaf_depth
01885       << std::endl;
01886     s << std::setprecision( 0 )
01887       << "Entities/Leaf      " WW d.leaf_ent.min WW avg_leaf_ent WW rms_leaf_ent WW d.leaf_ent.max WW std_leaf_ent
01888       << std::endl;
01889     s << std::setprecision( 3 )
01890       << "Child Volume Ratio " WW d.volume.min WW avg_vol_ratio WW rms_vol_ratio WW d.volume.max WW std_vol_ratio
01891       << std::endl;
01892     s << std::setprecision( 3 )
01893       << "Child Entity Ratio " WW d.entities.min WW avg_ent_ratio WW rms_ent_ratio WW d.entities.max WW std_ent_ratio
01894       << std::endl;
01895     s << std::setprecision( 3 )
01896       << "Box Radius Ratio   " WW d.radius.min WW avg_rad_ratio WW rms_rad_ratio WW d.radius.max WW std_rad_ratio
01897       << std::endl;
01898     s << std::setprecision( 0 ) << "Box volume         " WE d.vol.min WE avg_vol WE rms_vol WE d.vol.max WE std_vol
01899       << std::endl;
01900     s << std::setprecision( 0 ) << "Largest side area  " WE d.area.min WE avg_area WE rms_area WE d.area.max WE std_area
01901       << std::endl;
01902     s << std::setprecision( prec ) << std::endl;
01903 
01904     s << "Leaf Depth Histogram (Root depth is 0)" << std::endl;
01905     double f = 60.0 / max_leaf_per_depth;
01906     for( i = min_leaf_depth; i < d.leaf_depth.size(); ++i )
01907         s << std::setw( 2 ) << i << " " << std::setw( 5 ) << d.leaf_depth[i] << " |" << std::setfill( '*' )
01908           << std::setw( (int)floor( f * d.leaf_depth[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01909     s << std::endl;
01910 
01911     s << "Child/Parent Volume Ratio Histogram" << std::endl;
01912     f = 60.0 / *( std::max_element( d.volume.hist, d.volume.hist + 10 ) );
01913     for( i = 0; i < 10u; ++i )
01914         s << "0." << i << " " << std::setw( 5 ) << d.volume.hist[i] << " |" << std::setfill( '*' )
01915           << std::setw( (int)floor( f * d.volume.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01916     s << std::endl;
01917 
01918     s << "Child/Parent Entity Count Ratio Histogram" << std::endl;
01919     f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 10 ) );
01920     for( i = 0; i < 10u; ++i )
01921         s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
01922           << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01923     s << std::endl;
01924 
01925     s << "Inner/Outer Radius Ratio Histogram (~0.70 for cube)" << std::endl;
01926     // max radius ratio for a box is about 0.7071.  Move any boxes
01927     // in the .7 bucket into .6 and print .0 to .6.
01928     d.radius.hist[6] += d.radius.hist[7];
01929     f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 7 ) );
01930     for( i = 0; i < 7u; ++i )
01931         s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
01932           << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01933     s << std::endl;
01934 
01935     return MB_SUCCESS;
01936 }
01937 
01938 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines