MOAB: Mesh Oriented datABase  (version 5.3.0)
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 <cassert>
00033 #include <cmath>
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         {
00970             int_reg_callback.register_intersection( lastSet, *t, int_dist, search_win, int_type );
00971         }
00972     }
00973     return MB_SUCCESS;
00974 }
00975 
00976 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
00977                                                    std::vector< EntityHandle >& sets_out,
00978                                                    std::vector< EntityHandle >& facets_out, EntityHandle root_set,
00979                                                    const double tolerance, const double ray_point[3],
00980                                                    const double unit_ray_dir[3], IntersectSearchWindow& search_win,
00981                                                    IntRegCtxt& int_reg_callback, TrvStats* accum )
00982 
00983 {
00984     RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
00985                          accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
00986     ErrorCode rval = preorder_traverse( root_set, op, accum );
00987 
00988     distances_out = int_reg_callback.get_intersections();
00989     sets_out      = int_reg_callback.get_sets();
00990     facets_out    = int_reg_callback.get_facets();
00991 
00992     return rval;
00993 }
00994 
00995 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
00996                                                    std::vector< EntityHandle >& sets_out,
00997                                                    std::vector< EntityHandle >& facets_out, EntityHandle root_set,
00998                                                    double tolerance, const double ray_point[3],
00999                                                    const double unit_ray_dir[3], const double* ray_length,
01000                                                    TrvStats* accum )
01001 {
01002     IntRegCtxt int_reg_ctxt;
01003 
01004     OrientedBoxTreeTool::IntersectSearchWindow search_win( ray_length, (double*)0 );
01005 
01006     RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
01007                          accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_ctxt );
01008 
01009     ErrorCode rval = preorder_traverse( root_set, op, accum );
01010 
01011     if( MB_SUCCESS != rval ) { return rval; }
01012 
01013     distances_out = int_reg_ctxt.get_intersections();
01014     sets_out      = int_reg_ctxt.get_sets();
01015     facets_out    = int_reg_ctxt.get_facets();
01016 
01017     return MB_SUCCESS;
01018 }
01019 
01020 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( EntityHandle root_set, const double tolerance,
01021                                                    const double ray_point[3], const double unit_ray_dir[3],
01022                                                    IntersectSearchWindow& search_win, IntRegCtxt& int_reg_callback,
01023                                                    TrvStats* accum )
01024 
01025 {
01026     RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
01027                          accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
01028     return preorder_traverse( root_set, op, accum );
01029 }
01030 
01031 /********************** Closest Point code ***************/
01032 
01033 struct OBBTreeCPFrame
01034 {
01035     OBBTreeCPFrame( double d, EntityHandle n, EntityHandle s, int dp )
01036         : dist_sqr( d ), node( n ), mset( s ), depth( dp )
01037     {
01038     }
01039     double dist_sqr;
01040     EntityHandle node;
01041     EntityHandle mset;
01042     int depth;
01043 };
01044 
01045 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point, EntityHandle root, double* point_out,
01046                                                     EntityHandle& facet_out, EntityHandle* set_out, TrvStats* accum )
01047 {
01048     ErrorCode rval;
01049     const CartVect loc( point );
01050     double smallest_dist_sqr = std::numeric_limits< double >::max();
01051 
01052     EntityHandle current_set = 0;
01053     Range sets;
01054     std::vector< EntityHandle > children( 2 );
01055     std::vector< double > coords;
01056     std::vector< OBBTreeCPFrame > stack;
01057     int max_depth = -1;
01058 
01059     stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
01060 
01061     while( !stack.empty() )
01062     {
01063 
01064         // pop from top of stack
01065         EntityHandle node = stack.back().node;
01066         double dist_sqr   = stack.back().dist_sqr;
01067         current_set       = stack.back().mset;
01068         int current_depth = stack.back().depth;
01069         stack.pop_back();
01070 
01071         // If current best result is closer than the box, skip this tree node.
01072         if( dist_sqr > smallest_dist_sqr ) continue;
01073 
01074         // increment traversal statistics.
01075         if( accum )
01076         {
01077             accum->increment( current_depth );
01078             max_depth = std::max( max_depth, current_depth );
01079         }
01080 
01081         // Check if this node has a set associated with it
01082         if( set_out && !current_set )
01083         {
01084             sets.clear();
01085             rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
01086             if( MB_SUCCESS != rval ) return rval;
01087             if( !sets.empty() )
01088             {
01089                 if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
01090                 current_set = sets.front();
01091             }
01092         }
01093 
01094         // Get child boxes
01095         children.clear();
01096         rval = instance->get_child_meshsets( node, children );
01097         if( MB_SUCCESS != rval ) return rval;
01098 
01099         // if not a leaf node
01100         if( !children.empty() )
01101         {
01102             if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
01103 
01104             // get boxes
01105             OrientedBox box1, box2;
01106             rval = box( children[0], box1 );
01107             if( MB_SUCCESS != rval ) return rval;
01108             rval = box( children[1], box2 );
01109             if( MB_SUCCESS != rval ) return rval;
01110 
01111             // get distance from each box
01112             CartVect pt1, pt2;
01113             box1.closest_location_in_box( loc, pt1 );
01114             box2.closest_location_in_box( loc, pt2 );
01115             pt1 -= loc;
01116             pt2 -= loc;
01117             const double dsqr1 = pt1 % pt1;
01118             const double dsqr2 = pt2 % pt2;
01119 
01120             // push children on tree such that closer one is on top
01121             if( dsqr1 < dsqr2 )
01122             {
01123                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01124                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01125             }
01126             else
01127             {
01128                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01129                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01130             }
01131         }
01132         else
01133         {  // LEAF NODE
01134             if( accum ) { accum->increment_leaf( current_depth ); }
01135 
01136             Range facets;
01137             rval = instance->get_entities_by_dimension( node, 2, facets );
01138             if( MB_SUCCESS != rval ) return rval;
01139 
01140             const EntityHandle* conn = NULL;
01141             int len                  = 0;
01142             CartVect tmp, diff;
01143             for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
01144             {
01145                 rval = instance->get_connectivity( *i, conn, len, true );
01146                 if( MB_SUCCESS != rval ) return rval;
01147 
01148                 coords.resize( 3 * len );
01149                 rval = instance->get_coords( conn, len, &coords[0] );
01150                 if( MB_SUCCESS != rval ) return rval;
01151 
01152                 if( len == 3 )
01153                     GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
01154                 else
01155                     GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
01156 
01157                 diff     = tmp - loc;
01158                 dist_sqr = diff % diff;
01159                 if( dist_sqr < smallest_dist_sqr )
01160                 {
01161                     smallest_dist_sqr = dist_sqr;
01162                     facet_out         = *i;
01163                     tmp.get( point_out );
01164                     if( set_out ) *set_out = current_set;
01165                 }
01166             }
01167         }  // LEAF NODE
01168     }
01169 
01170     if( accum ) { accum->end_traversal( max_depth ); }
01171 
01172     return MB_SUCCESS;
01173 }
01174 
01175 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point, EntityHandle root, double tolerance,
01176                                                     std::vector< EntityHandle >& facets_out,
01177                                                     std::vector< EntityHandle >* sets_out, TrvStats* accum )
01178 {
01179     ErrorCode rval;
01180     const CartVect loc( point );
01181     double smallest_dist_sqr = std::numeric_limits< double >::max();
01182     double smallest_dist     = smallest_dist_sqr;
01183 
01184     EntityHandle current_set = 0;
01185     Range sets;
01186     std::vector< EntityHandle > children( 2 );
01187     std::vector< double > coords;
01188     std::vector< OBBTreeCPFrame > stack;
01189     int max_depth = -1;
01190 
01191     stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
01192 
01193     while( !stack.empty() )
01194     {
01195 
01196         // pop from top of stack
01197         EntityHandle node = stack.back().node;
01198         double dist_sqr   = stack.back().dist_sqr;
01199         current_set       = stack.back().mset;
01200         int current_depth = stack.back().depth;
01201         stack.pop_back();
01202 
01203         // If current best result is closer than the box, skip this tree node.
01204         if( dist_sqr > smallest_dist_sqr + tolerance ) continue;
01205 
01206         // increment traversal statistics.
01207         if( accum )
01208         {
01209             accum->increment( current_depth );
01210             max_depth = std::max( max_depth, current_depth );
01211         }
01212 
01213         // Check if this node has a set associated with it
01214         if( sets_out && !current_set )
01215         {
01216             sets.clear();
01217             rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
01218             if( MB_SUCCESS != rval ) return rval;
01219             if( !sets.empty() )
01220             {
01221                 if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
01222                 current_set = *sets.begin();
01223             }
01224         }
01225 
01226         // Get child boxes
01227         children.clear();
01228         rval = instance->get_child_meshsets( node, children );
01229         if( MB_SUCCESS != rval ) return rval;
01230 
01231         // if not a leaf node
01232         if( !children.empty() )
01233         {
01234             if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
01235 
01236             // get boxes
01237             OrientedBox box1, box2;
01238             rval = box( children[0], box1 );
01239             if( MB_SUCCESS != rval ) return rval;
01240             rval = box( children[1], box2 );
01241             if( MB_SUCCESS != rval ) return rval;
01242 
01243             // get distance from each box
01244             CartVect pt1, pt2;
01245             box1.closest_location_in_box( loc, pt1 );
01246             box2.closest_location_in_box( loc, pt2 );
01247             pt1 -= loc;
01248             pt2 -= loc;
01249             const double dsqr1 = pt1 % pt1;
01250             const double dsqr2 = pt2 % pt2;
01251 
01252             // push children on tree such that closer one is on top
01253             if( dsqr1 < dsqr2 )
01254             {
01255                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01256                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01257             }
01258             else
01259             {
01260                 stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
01261                 stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
01262             }
01263         }
01264         else
01265         {  // LEAF NODE
01266             if( accum ) { accum->increment_leaf( current_depth ); }
01267 
01268             Range facets;
01269             rval = instance->get_entities_by_dimension( node, 2, facets );
01270             if( MB_SUCCESS != rval ) return rval;
01271 
01272             const EntityHandle* conn = NULL;
01273             int len                  = 0;
01274             CartVect tmp, diff;
01275             for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
01276             {
01277                 rval = instance->get_connectivity( *i, conn, len, true );
01278                 if( MB_SUCCESS != rval ) return rval;
01279 
01280                 coords.resize( 3 * len );
01281                 rval = instance->get_coords( conn, len, &coords[0] );
01282                 if( MB_SUCCESS != rval ) return rval;
01283 
01284                 if( len == 3 )
01285                     GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
01286                 else
01287                     GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
01288 
01289                 diff     = tmp - loc;
01290                 dist_sqr = diff % diff;
01291                 if( dist_sqr < smallest_dist_sqr )
01292                 {
01293                     if( 0.5 * dist_sqr < 0.5 * smallest_dist_sqr + tolerance * ( 0.5 * tolerance - smallest_dist ) )
01294                     {
01295                         facets_out.clear();
01296                         if( sets_out ) sets_out->clear();
01297                     }
01298                     smallest_dist_sqr = dist_sqr;
01299                     smallest_dist     = sqrt( smallest_dist_sqr );
01300                     // put closest result at start of lists
01301                     facets_out.push_back( *i );
01302                     std::swap( facets_out.front(), facets_out.back() );
01303                     if( sets_out )
01304                     {
01305                         sets_out->push_back( current_set );
01306                         std::swap( sets_out->front(), sets_out->back() );
01307                     }
01308                 }
01309                 else if( dist_sqr <= smallest_dist_sqr + tolerance * ( tolerance + 2 * smallest_dist ) )
01310                 {
01311                     facets_out.push_back( *i );
01312                     if( sets_out ) sets_out->push_back( current_set );
01313                 }
01314             }
01315         }  // LEAF NODE
01316     }
01317 
01318     if( accum ) { accum->end_traversal( max_depth ); }
01319 
01320     return MB_SUCCESS;
01321 }
01322 
01323 /********************** Tree Printing Code ****************************/
01324 
01325 class TreeLayoutPrinter : public OrientedBoxTreeTool::Op
01326 {
01327   public:
01328     TreeLayoutPrinter( std::ostream& stream, Interface* instance );
01329 
01330     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
01331     virtual ErrorCode leaf( EntityHandle node );
01332 
01333   private:
01334     Interface* instance;
01335     std::ostream& outputStream;
01336     std::vector< bool > path;
01337 };
01338 
01339 TreeLayoutPrinter::TreeLayoutPrinter( std::ostream& stream, Interface* interface )
01340     : instance( interface ), outputStream( stream )
01341 {
01342 }
01343 
01344 ErrorCode TreeLayoutPrinter::visit( EntityHandle node, int depth, bool& descend )
01345 {
01346     descend = true;
01347 
01348     if( (unsigned)depth > path.size() )
01349     {
01350         // assert(depth+1 == path.size); // preorder traversal
01351         path.push_back( true );
01352     }
01353     else
01354     {
01355         path.resize( depth );
01356         if( depth ) path.back() = false;
01357     }
01358 
01359     for( unsigned i = 0; i + 1 < path.size(); ++i )
01360     {
01361         if( path[i] )
01362             outputStream << "|   ";
01363         else
01364             outputStream << "    ";
01365     }
01366     if( depth )
01367     {
01368         if( path.back() )
01369             outputStream << "+---";
01370         else
01371             outputStream << "\\---";
01372     }
01373     outputStream << instance->id_from_handle( node ) << std::endl;
01374     return MB_SUCCESS;
01375 }
01376 
01377 ErrorCode TreeLayoutPrinter::leaf( EntityHandle )
01378 {
01379     return MB_SUCCESS;
01380 }
01381 
01382 class TreeNodePrinter : public OrientedBoxTreeTool::Op
01383 {
01384   public:
01385     TreeNodePrinter( std::ostream& stream, bool list_contents, bool list_box, const char* id_tag_name,
01386                      OrientedBoxTreeTool* tool_ptr );
01387 
01388     virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
01389 
01390     virtual ErrorCode leaf( EntityHandle )
01391     {
01392         return MB_SUCCESS;
01393     }
01394 
01395   private:
01396     ErrorCode print_geometry( EntityHandle node );
01397     ErrorCode print_contents( EntityHandle node );
01398     ErrorCode print_counts( EntityHandle node );
01399 
01400     bool printContents;
01401     bool printGeometry;
01402     bool haveTag;
01403     Tag tag, gidTag, geomTag;
01404     Interface* instance;
01405     OrientedBoxTreeTool* tool;
01406     std::ostream& outputStream;
01407 };
01408 
01409 TreeNodePrinter::TreeNodePrinter( std::ostream& stream, bool list_contents, bool list_box, const char* id_tag_name,
01410                                   OrientedBoxTreeTool* tool_ptr )
01411     : printContents( list_contents ), printGeometry( list_box ), haveTag( false ), tag( 0 ), gidTag( 0 ), geomTag( 0 ),
01412       instance( tool_ptr->get_moab_instance() ), tool( tool_ptr ), outputStream( stream )
01413 {
01414     ErrorCode rval;
01415     if( id_tag_name )
01416     {
01417         rval = instance->tag_get_handle( id_tag_name, 1, MB_TYPE_INTEGER, tag );
01418         if( !rval )
01419         {
01420             std::cerr << "Could not get tag \"" << id_tag_name << "\"\n";
01421             stream << "Could not get tag \"" << id_tag_name << "\"\n";
01422         }
01423         else
01424         {
01425             haveTag = true;
01426         }
01427     }
01428 
01429     gidTag = instance->globalId_tag();
01430 
01431     rval = instance->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );
01432     if( MB_SUCCESS != rval ) geomTag = 0;
01433 }
01434 
01435 ErrorCode TreeNodePrinter::visit( EntityHandle node, int, bool& descend )
01436 {
01437     descend            = true;
01438     EntityHandle setid = instance->id_from_handle( node );
01439     outputStream << setid << ":" << std::endl;
01440 
01441     Range surfs;
01442     ErrorCode r3 = MB_SUCCESS;
01443     if( geomTag )
01444     {
01445         const int two         = 2;
01446         const void* tagdata[] = { &two };
01447         r3 = instance->get_entities_by_type_and_tag( node, MBENTITYSET, &geomTag, tagdata, 1, surfs );
01448 
01449         if( MB_SUCCESS == r3 && surfs.size() == 1 )
01450         {
01451             EntityHandle surf = *surfs.begin();
01452             int id;
01453             if( gidTag && MB_SUCCESS == instance->tag_get_data( gidTag, &surf, 1, &id ) )
01454                 outputStream << "  Surface " << id << std::endl;
01455             else
01456                 outputStream << "  Surface w/ unknown ID (" << surf << ")" << std::endl;
01457         }
01458     }
01459 
01460     ErrorCode r1 = printGeometry ? print_geometry( node ) : MB_SUCCESS;
01461     ErrorCode r2 = printContents ? print_contents( node ) : print_counts( node );
01462     outputStream << std::endl;
01463 
01464     if( MB_SUCCESS != r1 )
01465         return r1;
01466     else if( MB_SUCCESS != r2 )
01467         return r2;
01468     else
01469         return r3;
01470 }
01471 
01472 ErrorCode TreeNodePrinter::print_geometry( EntityHandle node )
01473 {
01474     OrientedBox box;
01475     ErrorCode rval = tool->box( node, box );
01476     if( MB_SUCCESS != rval ) return rval;
01477 
01478     CartVect length = box.dimensions();
01479 
01480     outputStream << box.center << "  Radius: " << box.inner_radius() << " - " << box.outer_radius() << std::endl
01481                  << '+' << box.axis( 0 ) << " : " << length[0] << std::endl
01482                  << 'x' << box.axis( 1 ) << " : " << length[1] << std::endl
01483                  << 'x' << box.axis( 2 ) << " : " << length[2] << std::endl;
01484     return MB_SUCCESS;
01485 }
01486 
01487 ErrorCode TreeNodePrinter::print_counts( EntityHandle node )
01488 {
01489     for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
01490     {
01491         int count      = 0;
01492         ErrorCode rval = instance->get_number_entities_by_type( node, type, count );
01493         if( MB_SUCCESS != rval ) return rval;
01494         if( count > 0 ) outputStream << " " << count << " " << CN::EntityTypeName( type ) << std::endl;
01495     }
01496     return MB_SUCCESS;
01497 }
01498 
01499 ErrorCode TreeNodePrinter::print_contents( EntityHandle node )
01500 {
01501     // list contents
01502     for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
01503     {
01504         Range range;
01505         ErrorCode rval = instance->get_entities_by_type( node, type, range );
01506         if( MB_SUCCESS != rval ) return rval;
01507         if( range.empty() ) continue;
01508         outputStream << " " << CN::EntityTypeName( type ) << " ";
01509         std::vector< int > ids( range.size() );
01510         if( haveTag )
01511         {
01512             rval = instance->tag_get_data( tag, range, &ids[0] );
01513             std::sort( ids.begin(), ids.end() );
01514         }
01515         else
01516         {
01517             Range::iterator ri              = range.begin();
01518             std::vector< int >::iterator vi = ids.begin();
01519             while( ri != range.end() )
01520             {
01521                 *vi = instance->id_from_handle( *ri );
01522                 ++ri;
01523                 ++vi;
01524             }
01525         }
01526 
01527         unsigned i = 0;
01528         for( ;; )
01529         {
01530             unsigned beg = i, end;
01531             do
01532             {
01533                 end = i++;
01534             } while( i < ids.size() && ids[end] + 1 == ids[i] );
01535             if( end == beg )
01536                 outputStream << ids[end];
01537             else if( end == beg + 1 )
01538                 outputStream << ids[beg] << ", " << ids[end];
01539             else
01540                 outputStream << ids[beg] << "-" << ids[end];
01541 
01542             if( i == ids.size() )
01543             {
01544                 outputStream << std::endl;
01545                 break;
01546             }
01547             else
01548                 outputStream << ", ";
01549         }
01550     }
01551 
01552     return MB_SUCCESS;
01553 }
01554 
01555 void OrientedBoxTreeTool::print( EntityHandle set, std::ostream& str, bool list, const char* tag )
01556 {
01557     TreeLayoutPrinter op1( str, instance );
01558     TreeNodePrinter op2( str, list, true, tag, this );
01559     ErrorCode r1 = preorder_traverse( set, op1 );
01560     str << std::endl;
01561     ErrorCode r2 = preorder_traverse( set, op2 );
01562     if( r1 != MB_SUCCESS || r2 != MB_SUCCESS )
01563     {
01564         std::cerr << "Errors encountered while printing tree\n";
01565         str << "Errors encountered while printing tree\n";
01566     }
01567 }
01568 
01569 /********************* Traversal Metrics Code  **************************/
01570 
01571 void OrientedBoxTreeTool::TrvStats::reset()
01572 {
01573     nodes_visited_count.clear();
01574     leaves_visited_count.clear();
01575     traversals_ended_count.clear();
01576     ray_tri_tests_count = 0;
01577 }
01578 
01579 void OrientedBoxTreeTool::TrvStats::increment( unsigned depth )
01580 {
01581 
01582     while( nodes_visited_count.size() <= depth )
01583     {
01584         nodes_visited_count.push_back( 0 );
01585         leaves_visited_count.push_back( 0 );
01586         traversals_ended_count.push_back( 0 );
01587     }
01588     nodes_visited_count[depth] += 1;
01589 }
01590 
01591 void OrientedBoxTreeTool::TrvStats::increment_leaf( unsigned depth )
01592 {
01593     // assume safe depth, because increment is called first
01594     leaves_visited_count[depth] += 1;
01595 }
01596 
01597 void OrientedBoxTreeTool::TrvStats::end_traversal( unsigned depth )
01598 {
01599     // assume safe depth, because increment is always called on a given
01600     // tree level first
01601     traversals_ended_count[depth] += 1;
01602 }
01603 
01604 void OrientedBoxTreeTool::TrvStats::print( std::ostream& str ) const
01605 {
01606 
01607     const std::string h1 = "OBBTree Depth";
01608     const std::string h2 = " - NodesVisited";
01609     const std::string h3 = " - LeavesVisited";
01610     const std::string h4 = " - TraversalsEnded";
01611 
01612     str << h1 << h2 << h3 << h4 << std::endl;
01613 
01614     unsigned num_visited = 0, num_leaves = 0, num_traversals = 0;
01615     for( unsigned i = 0; i < traversals_ended_count.size(); ++i )
01616     {
01617 
01618         num_visited += nodes_visited_count[i];
01619         num_leaves += leaves_visited_count[i];
01620         num_traversals += traversals_ended_count[i];
01621 
01622         str << std::setw( h1.length() ) << i << std::setw( h2.length() ) << nodes_visited_count[i]
01623             << std::setw( h3.length() ) << leaves_visited_count[i] << std::setw( h4.length() )
01624             << traversals_ended_count[i] << std::endl;
01625     }
01626 
01627     str << std::setw( h1.length() ) << "---- Totals:" << std::setw( h2.length() ) << num_visited
01628         << std::setw( h3.length() ) << num_leaves << std::setw( h4.length() ) << num_traversals << std::endl;
01629 
01630     if( ray_tri_tests_count )
01631     {
01632         str << std::setw( h1.length() ) << "---- Total ray-tri tests: " << ray_tri_tests_count << std::endl;
01633     }
01634 }
01635 
01636 /********************** Tree Statistics Code ****************************/
01637 
01638 class StatData
01639 {
01640   public:
01641     struct Ratio
01642     {
01643         double min, max, sum, sqr;
01644         int hist[10];
01645         Ratio()
01646             : min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() ), sum( 0.0 ),
01647               sqr( 0.0 )
01648         {
01649             hist[0] = hist[1] = hist[2] = hist[3] = hist[4] = hist[5] = hist[6] = hist[7] = hist[8] = hist[9] = 0;
01650         }
01651         void accum( double v )
01652         {
01653             if( v < min ) min = v;
01654             if( v > max ) max = v;
01655             sum += v;
01656             sqr += v * v;
01657             int i = (int)( 10 * v );
01658             if( i < 0 )
01659                 i = 0;
01660             else if( i > 9 )
01661                 i = 9;
01662             ++hist[i];
01663         }
01664     };
01665 
01666     template < typename T >
01667     struct Stat
01668     {
01669         T min, max;
01670         double sum, sqr;
01671         Stat() : sum( 0.0 ), sqr( 0.0 )
01672         {
01673             std::numeric_limits< T > lim;
01674             min = lim.max();
01675             if( lim.is_integer )
01676                 max = lim.min();
01677             else
01678                 max = -lim.max();
01679         }
01680         void accum( T v )
01681         {
01682             if( v < min ) min = v;
01683             if( v > max ) max = v;
01684             sum += v;
01685             sqr += (double)v * v;
01686         }
01687     };
01688 
01689     StatData() : count( 0 ) {}
01690 
01691     Ratio volume;
01692     Ratio entities;
01693     Ratio radius;
01694     Stat< unsigned > leaf_ent;
01695     Stat< double > vol;
01696     Stat< double > area;
01697     std::vector< unsigned > leaf_depth;
01698     unsigned count;
01699 };
01700 
01701 static int measure( const CartVect& v, double& result )
01702 {
01703     const double tol = 1e-6;
01704     int dims         = 0;
01705     result           = 1;
01706     for( int i = 0; i < 3; ++i )
01707         if( v[i] > tol )
01708         {
01709             ++dims;
01710             result *= v[i];
01711         }
01712     return dims;
01713 }
01714 
01715 ErrorCode OrientedBoxTreeTool::recursive_stats( OrientedBoxTreeTool* tool, Interface* inst, EntityHandle set, int depth,
01716                                                 StatData& data, unsigned& count_out, CartVect& dimensions_out )
01717 {
01718     ErrorCode rval;
01719     OrientedBox tmp_box;
01720     std::vector< EntityHandle > children( 2 );
01721     unsigned counts[2];
01722     bool isleaf;
01723 
01724     ++data.count;
01725 
01726     rval = tool->box( set, tmp_box );
01727     if( MB_SUCCESS != rval ) return rval;
01728     children.clear();
01729     rval = inst->get_child_meshsets( set, children );
01730     if( MB_SUCCESS != rval ) return rval;
01731     isleaf = children.empty();
01732     if( !isleaf && children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
01733 
01734     dimensions_out = tmp_box.dimensions();
01735     data.radius.accum( tmp_box.inner_radius() / tmp_box.outer_radius() );
01736     data.vol.accum( tmp_box.volume() );
01737     data.area.accum( tmp_box.area() );
01738 
01739     if( isleaf )
01740     {
01741         if( data.leaf_depth.size() <= (unsigned)depth ) data.leaf_depth.resize( depth + 1, 0 );
01742         ++data.leaf_depth[depth];
01743 
01744         int count;
01745         rval = inst->get_number_entities_by_handle( set, count );
01746         if( MB_SUCCESS != rval ) return rval;
01747         count_out = count;
01748         data.leaf_ent.accum( count_out );
01749     }
01750     else
01751     {
01752         for( int i = 0; i < 2; ++i )
01753         {
01754             CartVect dims;
01755             rval = recursive_stats( tool, inst, children[i], depth + 1, data, counts[i], dims );
01756             if( MB_SUCCESS != rval ) return rval;
01757             double this_measure, chld_measure;
01758             int this_dim = measure( dimensions_out, this_measure );
01759             int chld_dim = measure( dims, chld_measure );
01760             double ratio;
01761             if( chld_dim < this_dim )
01762                 ratio = 0;
01763             else
01764                 ratio = chld_measure / this_measure;
01765 
01766             data.volume.accum( ratio );
01767         }
01768         count_out = counts[0] + counts[1];
01769         data.entities.accum( (double)counts[0] / count_out );
01770         data.entities.accum( (double)counts[1] / count_out );
01771     }
01772     return MB_SUCCESS;
01773 }
01774 
01775 static inline double std_dev( double sqr, double sum, double count )
01776 {
01777     sum /= count;
01778     sqr /= count;
01779     return sqrt( sqr - sum * sum );
01780 }
01781 
01782 //#define WW <<std::setw(10)<<std::fixed<<
01783 #define WE << std::setw( 10 ) <<
01784 #define WW WE
01785 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set, unsigned& total_entities, double& rv, double& tot_node_volume,
01786                                       double& tot_to_root_volume, unsigned& tree_height, unsigned& node_count,
01787                                       unsigned& num_leaves )
01788 {
01789     StatData d;
01790     ErrorCode rval;
01791     unsigned i;
01792     CartVect total_dim;
01793 
01794     rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
01795     if( MB_SUCCESS != rval ) return rval;
01796 
01797     tree_height                 = d.leaf_depth.size();
01798     unsigned min_leaf_depth     = tree_height;
01799     num_leaves                  = 0;
01800     unsigned max_leaf_per_depth = 0;
01801     // double sum_leaf_depth = 0, sqr_leaf_depth = 0;
01802     for( i = 0; i < d.leaf_depth.size(); ++i )
01803     {
01804         unsigned val = d.leaf_depth[i];
01805         num_leaves += val;
01806         // sum_leaf_depth += (double)val*i;
01807         // sqr_leaf_depth += (double)val*i*i;
01808         if( val && i < min_leaf_depth ) min_leaf_depth = i;
01809         if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
01810     }
01811     rv                 = total_dim[0] * total_dim[1] * total_dim[2];
01812     tot_node_volume    = d.vol.sum;
01813     tot_to_root_volume = d.vol.sum / rv;
01814     node_count         = d.count;
01815 
01816     return MB_SUCCESS;
01817 }
01818 
01819 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set, std::ostream& s )
01820 {
01821     StatData d;
01822     ErrorCode rval;
01823     unsigned total_entities, i;
01824     CartVect total_dim;
01825 
01826     rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
01827     if( MB_SUCCESS != rval ) return rval;
01828 
01829     unsigned tree_height    = d.leaf_depth.size();
01830     unsigned min_leaf_depth = tree_height, num_leaves = 0;
01831     unsigned max_leaf_per_depth = 0;
01832     double sum_leaf_depth = 0, sqr_leaf_depth = 0;
01833     for( i = 0; i < d.leaf_depth.size(); ++i )
01834     {
01835         unsigned val = d.leaf_depth[i];
01836         num_leaves += val;
01837         sum_leaf_depth += (double)val * i;
01838         sqr_leaf_depth += (double)val * i * i;
01839         if( val && i < min_leaf_depth ) min_leaf_depth = i;
01840         if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
01841     }
01842     unsigned num_non_leaf = d.count - num_leaves;
01843 
01844     double rv = total_dim[0] * total_dim[1] * total_dim[2];
01845     s << "entities in tree:  " << total_entities << std::endl
01846       << "root volume:       " << rv << std::endl
01847       << "total node volume: " << d.vol.sum << std::endl
01848       << "total/root volume: " << d.vol.sum / rv << std::endl
01849       << "tree height:       " << tree_height << std::endl
01850       << "node count:        " << d.count << std::endl
01851       << "leaf count:        " << num_leaves << std::endl
01852       << std::endl;
01853 
01854     double avg_leaf_depth = sum_leaf_depth / num_leaves;
01855     double rms_leaf_depth = sqrt( sqr_leaf_depth / num_leaves );
01856     double std_leaf_depth = std_dev( sqr_leaf_depth, sum_leaf_depth, num_leaves );
01857 
01858     double avg_leaf_ent = d.leaf_ent.sum / num_leaves;
01859     double rms_leaf_ent = sqrt( d.leaf_ent.sqr / num_leaves );
01860     double std_leaf_ent = std_dev( d.leaf_ent.sqr, d.leaf_ent.sum, num_leaves );
01861 
01862     unsigned num_child = 2 * num_non_leaf;
01863 
01864     double avg_vol_ratio = d.volume.sum / num_child;
01865     double rms_vol_ratio = sqrt( d.volume.sqr / num_child );
01866     double std_vol_ratio = std_dev( d.volume.sqr, d.volume.sum, num_child );
01867 
01868     double avg_ent_ratio = d.entities.sum / num_child;
01869     double rms_ent_ratio = sqrt( d.entities.sqr / num_child );
01870     double std_ent_ratio = std_dev( d.entities.sqr, d.entities.sum, num_child );
01871 
01872     double avg_rad_ratio = d.radius.sum / d.count;
01873     double rms_rad_ratio = sqrt( d.radius.sqr / d.count );
01874     double std_rad_ratio = std_dev( d.radius.sqr, d.radius.sum, d.count );
01875 
01876     double avg_vol = d.vol.sum / d.count;
01877     double rms_vol = sqrt( d.vol.sqr / d.count );
01878     double std_vol = std_dev( d.vol.sqr, d.vol.sum, d.count );
01879 
01880     double avg_area = d.area.sum / d.count;
01881     double rms_area = sqrt( d.area.sqr / d.count );
01882     double std_area = std_dev( d.area.sqr, d.area.sum, d.count );
01883 
01884     int prec = s.precision();
01885     s << "                   " WW "Minimum" WW "Average" WW "RMS" WW "Maximum" WW "Std.Dev." << std::endl;
01886     s << std::setprecision( 1 )
01887       << "Leaf Depth         " WW min_leaf_depth WW avg_leaf_depth WW rms_leaf_depth WW d.leaf_depth.size() -
01888              1 WW std_leaf_depth
01889       << std::endl;
01890     s << std::setprecision( 0 )
01891       << "Entities/Leaf      " WW d.leaf_ent.min WW avg_leaf_ent WW rms_leaf_ent WW d.leaf_ent.max WW std_leaf_ent
01892       << std::endl;
01893     s << std::setprecision( 3 )
01894       << "Child Volume Ratio " WW d.volume.min WW avg_vol_ratio WW rms_vol_ratio WW d.volume.max WW std_vol_ratio
01895       << std::endl;
01896     s << std::setprecision( 3 )
01897       << "Child Entity Ratio " WW d.entities.min WW avg_ent_ratio WW rms_ent_ratio WW d.entities.max WW std_ent_ratio
01898       << std::endl;
01899     s << std::setprecision( 3 )
01900       << "Box Radius Ratio   " WW d.radius.min WW avg_rad_ratio WW rms_rad_ratio WW d.radius.max WW std_rad_ratio
01901       << std::endl;
01902     s << std::setprecision( 0 ) << "Box volume         " WE d.vol.min WE avg_vol WE rms_vol WE d.vol.max WE std_vol
01903       << std::endl;
01904     s << std::setprecision( 0 ) << "Largest side area  " WE d.area.min WE avg_area WE rms_area WE d.area.max WE std_area
01905       << std::endl;
01906     s << std::setprecision( prec ) << std::endl;
01907 
01908     s << "Leaf Depth Histogram (Root depth is 0)" << std::endl;
01909     double f = 60.0 / max_leaf_per_depth;
01910     for( i = min_leaf_depth; i < d.leaf_depth.size(); ++i )
01911         s << std::setw( 2 ) << i << " " << std::setw( 5 ) << d.leaf_depth[i] << " |" << std::setfill( '*' )
01912           << std::setw( (int)floor( f * d.leaf_depth[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01913     s << std::endl;
01914 
01915     s << "Child/Parent Volume Ratio Histogram" << std::endl;
01916     f = 60.0 / *( std::max_element( d.volume.hist, d.volume.hist + 10 ) );
01917     for( i = 0; i < 10u; ++i )
01918         s << "0." << i << " " << std::setw( 5 ) << d.volume.hist[i] << " |" << std::setfill( '*' )
01919           << std::setw( (int)floor( f * d.volume.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01920     s << std::endl;
01921 
01922     s << "Child/Parent Entity Count Ratio Histogram" << std::endl;
01923     f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 10 ) );
01924     for( i = 0; i < 10u; ++i )
01925         s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
01926           << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01927     s << std::endl;
01928 
01929     s << "Inner/Outer Radius Ratio Histogram (~0.70 for cube)" << std::endl;
01930     // max radius ratio for a box is about 0.7071.  Move any boxes
01931     // in the .7 bucket into .6 and print .0 to .6.
01932     d.radius.hist[6] += d.radius.hist[7];
01933     f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 7 ) );
01934     for( i = 0; i < 7u; ++i )
01935         s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
01936           << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
01937     s << std::endl;
01938 
01939     return MB_SUCCESS;
01940 }
01941 
01942 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines