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