Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
OrientedBoxTreeTool.cpp
Go to the documentation of this file.
1 /*
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 /**\file OrientedBox.hpp
17  *\author Jason Kraftcheck ([email protected])
18  *\date 2006-07-18
19  */
20 
21 #include "moab/Interface.hpp"
22 #include "Internals.hpp"
23 #include "moab/OrientedBoxTreeTool.hpp"
24 #include "moab/Range.hpp"
25 #include "moab/CN.hpp"
26 #include "moab/GeomUtil.hpp"
27 #include "MBTagConventions.hpp"
28 #include <iostream>
29 #include <iomanip>
30 #include <algorithm>
31 #include <limits>
32 #include <cassert>
33 #include <cmath>
34 
35 //#define MB_OBB_USE_VECTOR_QUERIES
36 //#define MB_OBB_USE_TYPE_QUERIES
37 
38 namespace moab
39 {
40 
41 #if defined( MB_OBB_USE_VECTOR_QUERIES ) && defined( MB_OBB_USE_TYPE_QUERIES )
42 #undef MB_OBB_USE_TYPE_QUERIES
43 #endif
44 
45 const char DEFAULT_TAG_NAME[] = "OBB";
46 
47 OrientedBoxTreeTool::Op::~Op() {}
48 
49 OrientedBoxTreeTool::OrientedBoxTreeTool( Interface* i, const char* tag_name, bool destroy_created_trees )
50  : instance( i ), cleanUpTrees( destroy_created_trees )
51 {
52  if( !tag_name ) tag_name = DEFAULT_TAG_NAME;
53  ErrorCode rval = OrientedBox::tag_handle( tagHandle, instance, tag_name );
54  if( MB_SUCCESS != rval ) tagHandle = 0;
55 }
56 
57 OrientedBoxTreeTool::~OrientedBoxTreeTool()
58 {
59  if( !cleanUpTrees ) return;
60 
61  while( !createdTrees.empty() )
62  {
63  EntityHandle tree = createdTrees.back();
64  // make sure this is a tree (rather than some other, stale handle)
65  const void* data_ptr = 0;
66  ErrorCode rval = instance->tag_get_by_ptr( tagHandle, &tree, 1, &data_ptr );
67  if( MB_SUCCESS == rval ) rval = delete_tree( tree );
68  if( MB_SUCCESS != rval ) createdTrees.pop_back();
69  }
70 }
71 
72 OrientedBoxTreeTool::Settings::Settings()
73  : max_leaf_entities( 8 ), max_depth( 0 ), worst_split_ratio( 0.95 ), best_split_ratio( 0.4 ),
74  set_options( MESHSET_SET )
75 {
76 }
77 
78 bool OrientedBoxTreeTool::Settings::valid() const
79 {
80  return max_leaf_entities > 0 && max_depth >= 0 && worst_split_ratio <= 1.0 && best_split_ratio >= 0.0 &&
81  worst_split_ratio >= best_split_ratio;
82 }
83 
84 /********************** Simple Tree Access Methods ****************************/
85 
86 ErrorCode OrientedBoxTreeTool::box( EntityHandle set, OrientedBox& obb )
87 {
88  return instance->tag_get_data( tagHandle, &set, 1, &obb );
89 }
90 
91 ErrorCode OrientedBoxTreeTool::box( EntityHandle set,
92  double center[3],
93  double axis1[3],
94  double axis2[3],
95  double axis3[3] )
96 {
97  OrientedBox obb;
98  ErrorCode rval = this->box( set, obb );
99  obb.center.get( center );
100  obb.scaled_axis( 0 ).get( axis1 );
101  obb.scaled_axis( 1 ).get( axis2 );
102  obb.scaled_axis( 2 ).get( axis3 );
103  return rval;
104 }
105 
106 /********************** Tree Construction Code ****************************/
107 
108 struct OrientedBoxTreeTool::SetData
109 {
110  EntityHandle handle;
111  OrientedBox::CovarienceData box_data;
112  // Range vertices;
113 };
114 
115 ErrorCode OrientedBoxTreeTool::build( const Range& entities, EntityHandle& set_handle_out, const Settings* settings )
116 {
117  if( !entities.all_of_dimension( 2 ) ) return MB_TYPE_OUT_OF_RANGE;
118  if( settings && !settings->valid() ) return MB_FAILURE;
119 
120  return build_tree( entities, set_handle_out, 0, settings ? *settings : Settings() );
121 }
122 
123 ErrorCode OrientedBoxTreeTool::join_trees( const Range& sets, EntityHandle& set_handle_out, const Settings* settings )
124 {
125  if( !sets.all_of_type( MBENTITYSET ) ) return MB_TYPE_OUT_OF_RANGE;
126  if( settings && !settings->valid() ) return MB_FAILURE;
127 
128  // Build initial set data list.
129  std::list< SetData > data;
130  for( Range::const_iterator i = sets.begin(); i != sets.end(); ++i )
131  {
132  Range elements;
133  ErrorCode rval = instance->get_entities_by_dimension( *i, 2, elements, true );
134  if( MB_SUCCESS != rval ) return rval;
135  if( elements.empty() ) continue;
136 
137  data.push_back( SetData() );
138  SetData& set_data = data.back();
139  set_data.handle = *i;
140  rval = OrientedBox::covariance_data_from_tris( set_data.box_data, instance, elements );
141  if( MB_SUCCESS != rval ) return rval;
142  }
143 
144  ErrorCode result = build_sets( data, set_handle_out, 0, settings ? *settings : Settings() );
145  if( MB_SUCCESS != result ) return result;
146 
147  for( Range::reverse_iterator i = sets.rbegin(); i != sets.rend(); ++i )
148  {
149  createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), *i ), createdTrees.end() );
150  }
151  createdTrees.push_back( set_handle_out );
152  return MB_SUCCESS;
153 }
154 
155 /**\brief Split triangles by which side of a plane they are on
156  *
157  * Given a plane specified as a bisecting plane normal to one
158  * of the axes of a box, split triangles based on which side
159  * of the plane they are on.
160  *\param instance MOAB instance
161  *\param box The oriented box containing all the entities
162  *\param axis The axis for which the split plane is orthogonal
163  *\param entities The entities intersecting plane
164  *\param left_list Output, entities to the left of the plane
165  *\param right_list Output, entities to the right of the plane
166  */
167 static ErrorCode split_box( Interface* instance,
168  const OrientedBox& box,
169  int axis,
170  const Range& entities,
171  Range& left_list,
172  Range& right_list )
173 {
174  ErrorCode rval;
175  left_list.clear();
176  right_list.clear();
177 
178  std::vector< CartVect > coords;
179  for( Range::reverse_iterator i = entities.rbegin(); i != entities.rend(); ++i )
180  {
181  const EntityHandle* conn = NULL;
182  int conn_len = 0;
183  rval = instance->get_connectivity( *i, conn, conn_len );
184  if( MB_SUCCESS != rval ) return rval;
185 
186  coords.resize( conn_len );
187  rval = instance->get_coords( conn, conn_len, coords[0].array() );
188  if( MB_SUCCESS != rval ) return rval;
189 
190  CartVect centroid( 0.0 );
191  for( int j = 0; j < conn_len; ++j )
192  centroid += coords[j];
193  centroid /= conn_len;
194 
195  if( ( box.axis( axis ) % ( centroid - box.center ) ) < 0.0 )
196  left_list.insert( *i );
197  else
198  right_list.insert( *i );
199  }
200 
201  return MB_SUCCESS;
202 }
203 
204 ErrorCode OrientedBoxTreeTool::build_tree( const Range& entities,
205  EntityHandle& set,
206  int depth,
207  const Settings& settings )
208 {
209  OrientedBox tmp_box;
210  ErrorCode rval;
211 
212  if( entities.empty() )
213  {
214  Matrix3 axis;
215  tmp_box = OrientedBox( axis, CartVect( 0. ) );
216  }
217  else
218  {
219  rval = OrientedBox::compute_from_2d_cells( tmp_box, instance, entities );
220  if( MB_SUCCESS != rval ) return rval;
221  }
222 
223  // create an entity set for the tree node
224  rval = instance->create_meshset( settings.set_options, set );
225  if( MB_SUCCESS != rval ) return rval;
226 
227  rval = instance->tag_set_data( tagHandle, &set, 1, &tmp_box );
228  if( MB_SUCCESS != rval )
229  {
230  delete_tree( set );
231  return rval;
232  }
233 
234  // check if should create children
235  bool leaf = true;
236  ++depth;
237  if( ( !settings.max_depth || depth < settings.max_depth ) &&
238  entities.size() > (unsigned)settings.max_leaf_entities )
239  {
240  // try splitting with planes normal to each axis of the box
241  // until we find an acceptable split
242  double best_ratio = settings.worst_split_ratio; // worst case ratio
243  Range best_left_list, best_right_list;
244  // Axes are sorted from shortest to longest, so search backwards
245  for( int axis = 2; best_ratio > settings.best_split_ratio && axis >= 0; --axis )
246  {
247  Range left_list, right_list;
248 
249  rval = split_box( instance, tmp_box, axis, entities, left_list, right_list );
250  if( MB_SUCCESS != rval )
251  {
252  delete_tree( set );
253  return rval;
254  }
255 
256  double ratio = fabs( (double)right_list.size() - left_list.size() ) / entities.size();
257 
258  if( ratio < best_ratio )
259  {
260  best_ratio = ratio;
261  best_left_list.swap( left_list );
262  best_right_list.swap( right_list );
263  }
264  }
265 
266  // create children
267  if( !best_left_list.empty() )
268  {
269  EntityHandle child = 0;
270 
271  rval = build_tree( best_left_list, child, depth, settings );
272  if( MB_SUCCESS != rval )
273  {
274  delete_tree( set );
275  return rval;
276  }
277  rval = instance->add_child_meshset( set, child );
278  if( MB_SUCCESS != rval )
279  {
280  delete_tree( set );
281  delete_tree( child );
282  return rval;
283  }
284 
285  rval = build_tree( best_right_list, child, depth, settings );
286  if( MB_SUCCESS != rval )
287  {
288  delete_tree( set );
289  return rval;
290  }
291  rval = instance->add_child_meshset( set, child );
292  if( MB_SUCCESS != rval )
293  {
294  delete_tree( set );
295  delete_tree( child );
296  return rval;
297  }
298 
299  leaf = false;
300  }
301  }
302 
303  if( leaf )
304  {
305  rval = instance->add_entities( set, entities );
306  if( MB_SUCCESS != rval )
307  {
308  delete_tree( set );
309  return rval;
310  }
311  }
312 
313  createdTrees.push_back( set );
314  return MB_SUCCESS;
315 }
316 
317 static ErrorCode split_sets( Interface*,
318  const OrientedBox& box,
319  int axis,
320  const std::list< OrientedBoxTreeTool::SetData >& sets,
321  std::list< OrientedBoxTreeTool::SetData >& left,
322  std::list< OrientedBoxTreeTool::SetData >& right )
323 {
324  left.clear();
325  right.clear();
326 
327  std::list< OrientedBoxTreeTool::SetData >::const_iterator i;
328  for( i = sets.begin(); i != sets.end(); ++i )
329  {
330  CartVect centroid( i->box_data.center / i->box_data.area );
331  if( ( box.axis( axis ) % ( centroid - box.center ) ) < 0.0 )
332  left.push_back( *i );
333  else
334  right.push_back( *i );
335  }
336 
337  return MB_SUCCESS;
338 }
339 
340 ErrorCode OrientedBoxTreeTool::build_sets( std::list< SetData >& sets,
341  EntityHandle& node_set,
342  int depth,
343  const Settings& settings )
344 {
345  ErrorCode rval;
346  int count = sets.size();
347  if( 0 == count ) return MB_FAILURE;
348 
349  // calculate box
350  OrientedBox obox;
351 
352  // make vector go out of scope when done, so memory is released
353  {
354  Range elems;
355  std::vector< OrientedBox::CovarienceData > data( sets.size() );
356  data.clear();
357  for( std::list< SetData >::iterator i = sets.begin(); i != sets.end(); ++i )
358  {
359  data.push_back( i->box_data );
360  rval = instance->get_entities_by_dimension( i->handle, 2, elems, true );
361  if( MB_SUCCESS != rval ) return rval;
362  }
363 
364  Range points;
365  rval = instance->get_adjacencies( elems, 0, false, points, Interface::UNION );
366  if( MB_SUCCESS != rval ) return rval;
367 
368  rval = OrientedBox::compute_from_covariance_data( obox, instance, &data[0], data.size(), points );
369  if( MB_SUCCESS != rval ) return rval;
370  }
371 
372  // If only one set in list...
373  if( count == 1 )
374  {
375  node_set = sets.front().handle;
376  return instance->tag_set_data( tagHandle, &node_set, 1, &obox );
377  }
378 
379  // create an entity set for the tree node
380  rval = instance->create_meshset( settings.set_options, node_set );
381  if( MB_SUCCESS != rval ) return rval;
382 
383  rval = instance->tag_set_data( tagHandle, &node_set, 1, &obox );
384  if( MB_SUCCESS != rval )
385  {
386  delete_tree( node_set );
387  return rval;
388  }
389 
390  double best_ratio = 2.0;
391  std::list< SetData > best_left_list, best_right_list;
392  for( int axis = 0; axis < 2; ++axis )
393  {
394  std::list< SetData > left_list, right_list;
395  rval = split_sets( instance, obox, axis, sets, left_list, right_list );
396  if( MB_SUCCESS != rval )
397  {
398  delete_tree( node_set );
399  return rval;
400  }
401 
402  double ratio = fabs( (double)right_list.size() - left_list.size() ) / sets.size();
403 
404  if( ratio < best_ratio )
405  {
406  best_ratio = ratio;
407  best_left_list.swap( left_list );
408  best_right_list.swap( right_list );
409  }
410  }
411 
412  // We must subdivide the list of sets, because we want to guarantee that
413  // there is a node in the tree corresponding to each of the sets. So if
414  // we couldn't find a usable split plane, just split them in an arbitrary
415  // manner.
416  if( best_left_list.empty() || best_right_list.empty() )
417  {
418  best_left_list.clear();
419  best_right_list.clear();
420  std::list< SetData >* lists[2] = { &best_left_list, &best_right_list };
421  int i = 0;
422  while( !sets.empty() )
423  {
424  lists[i]->push_back( sets.front() );
425  sets.pop_front();
426  i = 1 - i;
427  }
428  }
429  else
430  {
431  sets.clear(); // release memory before recursion
432  }
433 
434  // Create child sets
435 
436  EntityHandle child = 0;
437 
438  rval = build_sets( best_left_list, child, depth + 1, settings );
439  if( MB_SUCCESS != rval )
440  {
441  delete_tree( node_set );
442  return rval;
443  }
444  rval = instance->add_child_meshset( node_set, child );
445  if( MB_SUCCESS != rval )
446  {
447  delete_tree( node_set );
448  delete_tree( child );
449  return rval;
450  }
451 
452  rval = build_sets( best_right_list, child, depth + 1, settings );
453  if( MB_SUCCESS != rval )
454  {
455  delete_tree( node_set );
456  return rval;
457  }
458  rval = instance->add_child_meshset( node_set, child );
459  if( MB_SUCCESS != rval )
460  {
461  delete_tree( node_set );
462  delete_tree( child );
463  return rval;
464  }
465 
466  return MB_SUCCESS;
467 }
468 
469 ErrorCode OrientedBoxTreeTool::delete_tree( EntityHandle set )
470 {
471  std::vector< EntityHandle > children;
472  ErrorCode rval = instance->get_child_meshsets( set, children, 0 );
473  if( MB_SUCCESS != rval ) return rval;
474 
475  createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), set ), createdTrees.end() );
476  children.insert( children.begin(), set );
477  return instance->delete_entities( &children[0], children.size() );
478 }
479 
480 ErrorCode OrientedBoxTreeTool::remove_root( EntityHandle root )
481 {
482  std::vector< EntityHandle >::iterator i = find( createdTrees.begin(), createdTrees.end(), root );
483  if( i != createdTrees.end() )
484  {
485  createdTrees.erase( i );
486  return MB_SUCCESS;
487  }
488  else
489  return MB_ENTITY_NOT_FOUND;
490 }
491 
492 /********************** Generic Tree Traversal ****************************/
493 struct Data
494 {
495  EntityHandle set;
496  int depth;
497 };
498 ErrorCode OrientedBoxTreeTool::preorder_traverse( EntityHandle set, Op& operation, TrvStats* accum )
499 {
500  ErrorCode rval;
501  std::vector< EntityHandle > children;
502  std::vector< Data > the_stack;
503  Data data = { set, 0 };
504  the_stack.push_back( data );
505  int max_depth = -1;
506 
507  while( !the_stack.empty() )
508  {
509  data = the_stack.back();
510  the_stack.pop_back();
511 
512  // increment traversal statistics
513  if( accum )
514  {
515  accum->increment( data.depth );
516  max_depth = std::max( max_depth, data.depth );
517  }
518 
519  bool descend = true;
520  rval = operation.visit( data.set, data.depth, descend );
521  assert( MB_SUCCESS == rval );
522  if( MB_SUCCESS != rval ) return rval;
523 
524  if( !descend ) continue;
525 
526  children.clear();
527  rval = instance->get_child_meshsets( data.set, children );
528  assert( MB_SUCCESS == rval );
529  if( MB_SUCCESS != rval ) return rval;
530  if( children.empty() )
531  {
532  if( accum )
533  {
534  accum->increment_leaf( data.depth );
535  }
536  rval = operation.leaf( data.set );
537  assert( MB_SUCCESS == rval );
538  if( MB_SUCCESS != rval ) return rval;
539  }
540  else if( children.size() == 2 )
541  {
542  data.depth++;
543  data.set = children[0];
544  the_stack.push_back( data );
545  data.set = children[1];
546  the_stack.push_back( data );
547  }
548  else
549  return MB_MULTIPLE_ENTITIES_FOUND;
550  }
551 
552  if( accum )
553  {
554  accum->end_traversal( max_depth );
555  }
556 
557  return MB_SUCCESS;
558 }
559 
560 /********************** General Sphere/Triangle Intersection ***************/
561 
562 struct OBBTreeSITFrame
563 {
564  OBBTreeSITFrame( EntityHandle n, EntityHandle s, int dp ) : node( n ), surf( s ), depth( dp ) {}
565  EntityHandle node;
566  EntityHandle surf;
567  int depth;
568 };
569 
570 ErrorCode OrientedBoxTreeTool::sphere_intersect_triangles( const double* center_v,
571  double radius,
572  EntityHandle tree_root,
573  std::vector< EntityHandle >& facets_out,
574  std::vector< EntityHandle >* sets_out,
575  TrvStats* accum )
576 {
577  const double radsqr = radius * radius;
578  OrientedBox b;
579  ErrorCode rval;
580  Range sets;
581  const CartVect center( center_v );
582  CartVect closest, coords[3];
583  const EntityHandle* conn;
584  int num_conn;
585 #ifndef MB_OBB_USE_VECTOR_QUERIES
586  Range tris;
587  Range::const_iterator t;
588 #else
589  std::vector< EntityHandle > tris;
590  std::vector< EntityHandle >::const_iterator t;
591 #endif
592 
593  std::vector< OBBTreeSITFrame > stack;
594  std::vector< EntityHandle > children;
595  stack.reserve( 30 );
596  stack.push_back( OBBTreeSITFrame( tree_root, 0, 0 ) );
597 
598  int max_depth = -1;
599 
600  while( !stack.empty() )
601  {
602  EntityHandle surf = stack.back().surf;
603  EntityHandle node = stack.back().node;
604  int current_depth = stack.back().depth;
605  stack.pop_back();
606 
607  // increment traversal statistics.
608  if( accum )
609  {
610  accum->increment( current_depth );
611  max_depth = std::max( max_depth, current_depth );
612  }
613 
614  if( !surf && sets_out )
615  {
616  rval = get_moab_instance()->get_entities_by_type( node, MBENTITYSET, sets );
617  if( !sets.empty() ) surf = sets.front();
618  sets.clear();
619  }
620 
621  // check if sphere intersects box
622  rval = box( node, b );
623  if( MB_SUCCESS != rval ) return rval;
624  b.closest_location_in_box( center, closest );
625  closest -= center;
626  if( ( closest % closest ) > radsqr ) continue;
627 
628  // push child boxes on stack
629  children.clear();
630  rval = instance->get_child_meshsets( node, children );
631  if( MB_SUCCESS != rval ) return rval;
632  if( !children.empty() )
633  {
634  assert( children.size() == 2 );
635  stack.push_back( OBBTreeSITFrame( children[0], surf, current_depth + 1 ) );
636  stack.push_back( OBBTreeSITFrame( children[1], surf, current_depth + 1 ) );
637  continue;
638  }
639 
640  if( accum )
641  {
642  accum->increment_leaf( current_depth );
643  }
644  // if leaf, intersect sphere with triangles
645 #ifndef MB_OBB_USE_VECTOR_QUERIES
646 #ifdef MB_OBB_USE_TYPE_QUERIES
647  rval = get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
648 #else
649  rval = get_moab_instance()->get_entities_by_handle( node, tris );
650 #endif
651  t = tris.begin();
652 #else
653  rval = get_moab_instance()->get_entities_by_handle( node, tris );
654  t = tris.lower_bound( MBTRI );
655 #endif
656  if( MB_SUCCESS != rval ) return rval;
657 
658  for( t = tris.begin(); t != tris.end(); ++t )
659  {
660 #ifndef MB_OBB_USE_VECTOR_QUERIES
661  if( TYPE_FROM_HANDLE( *t ) != MBTRI ) break;
662 #elif !defined( MB_OBB_USE_TYPE_QUERIES )
663  if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
664 #endif
665  rval = get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
666  if( MB_SUCCESS != rval ) return rval;
667  if( num_conn != 3 ) continue;
668 
669  rval = get_moab_instance()->get_coords( conn, num_conn, coords[0].array() );
670  if( MB_SUCCESS != rval ) return rval;
671 
672  GeomUtil::closest_location_on_tri( center, coords, closest );
673  closest -= center;
674  if( ( closest % closest ) <= radsqr &&
675  std::find( facets_out.begin(), facets_out.end(), *t ) == facets_out.end() )
676  {
677  facets_out.push_back( *t );
678  if( sets_out ) sets_out->push_back( surf );
679  }
680  }
681  }
682 
683  if( accum )
684  {
685  accum->end_traversal( max_depth );
686  }
687 
688  return MB_SUCCESS;
689 }
690 
691 /********************** General Ray/Tree and Ray/Triangle Intersection ***************/
692 
693 class RayIntersector : public OrientedBoxTreeTool::Op
694 {
695  private:
696  OrientedBoxTreeTool* tool;
697  const CartVect b, m;
698  const double* len;
699  const double tol;
700  Range& boxes;
701 
702  public:
703  RayIntersector( OrientedBoxTreeTool* tool_ptr,
704  const double* ray_point,
705  const double* unit_ray_dir,
706  const double* ray_length,
707  double tolerance,
708  Range& leaf_boxes )
709  : tool( tool_ptr ), b( ray_point ), m( unit_ray_dir ), len( ray_length ), tol( tolerance ), boxes( leaf_boxes )
710  {
711  }
712 
713  virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
714  virtual ErrorCode leaf( EntityHandle node );
715 };
716 
717 //#include <stdio.h>
718 // inline void dump_fragmentation( const Range& range ) {
719 // static FILE* file = fopen( "fragmentation", "w" );
720 // unsigned ranges = 0, entities = 0;
721 // for (Range::const_pair_iterator i = range.const_pair_begin(); i != range.const_pair_end(); ++i)
722 // {
723 // ++ranges;
724 // entities += i->second - i->first + 1;
725 // }
726 // fprintf( file, "%u %u\n", ranges, entities );
727 //}
728 
729 ErrorCode OrientedBoxTreeTool::ray_intersect_triangles( std::vector< double >& intersection_distances_out,
730  std::vector< EntityHandle >& intersection_facets_out,
731  const Range& boxes,
732  double /*tolerance*/,
733  const double ray_point[3],
734  const double unit_ray_dir[3],
735  const double* ray_length,
736  unsigned int* raytri_test_count )
737 {
738  ErrorCode rval;
739  intersection_distances_out.clear();
740 #ifdef MB_OBB_USE_VECTOR_QUERIES
741  std::vector< EntityHandle > tris;
742 #endif
743 
744  const CartVect point( ray_point );
745  const CartVect dir( unit_ray_dir );
746 
747  for( Range::iterator b = boxes.begin(); b != boxes.end(); ++b )
748  {
749 #ifndef MB_OBB_USE_VECTOR_QUERIES
750  Range tris;
751 #ifdef MB_OBB_USE_TYPE_QUERIES
752  rval = instance->get_entities_by_type( *b, MBTRI, tris );
753 #else
754  rval = instance->get_entities_by_handle( *b, tris );
755 #endif
756 #else
757  tris.clear();
758  rval = instance->get_entities_by_handle( *b, tris );
759 #endif
760  if( MB_SUCCESS != rval ) return rval;
761  // dump_fragmentation( tris );
762 
763 #ifndef MB_OBB_USE_VECTOR_QUERIES
764  for( Range::iterator t = tris.begin(); t != tris.end(); ++t )
765 #else
766  for( std::vector< EntityHandle >::iterator t = tris.begin(); t != tris.end(); ++t )
767 #endif
768  {
769 #ifndef MB_OBB_USE_TYPE_QUERIES
770  if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
771 #endif
772 
773  const EntityHandle* conn = NULL;
774  int len = 0;
775  rval = instance->get_connectivity( *t, conn, len, true );
776  if( MB_SUCCESS != rval ) return rval;
777 
778  CartVect coords[3];
779  rval = instance->get_coords( conn, 3, coords[0].array() );
780  if( MB_SUCCESS != rval ) return rval;
781 
782  if( raytri_test_count ) *raytri_test_count += 1;
783 
784  double td;
785  if( GeomUtil::plucker_ray_tri_intersect( coords, point, dir, td, ray_length ) )
786  {
787  intersection_distances_out.push_back( td );
788  intersection_facets_out.push_back( *t );
789  }
790  }
791  }
792 
793  return MB_SUCCESS;
794 }
795 
796 ErrorCode OrientedBoxTreeTool::ray_intersect_triangles( std::vector< double >& intersection_distances_out,
797  std::vector< EntityHandle >& intersection_facets_out,
798  EntityHandle root_set,
799  double tolerance,
800  const double ray_point[3],
801  const double unit_ray_dir[3],
802  const double* ray_length,
803  TrvStats* accum )
804 {
805  Range boxes;
806  ErrorCode rval;
807 
808  rval = ray_intersect_boxes( boxes, root_set, tolerance, ray_point, unit_ray_dir, ray_length, accum );
809  if( MB_SUCCESS != rval ) return rval;
810 
811  return ray_intersect_triangles( intersection_distances_out, intersection_facets_out, boxes, tolerance, ray_point,
812  unit_ray_dir, ray_length, accum ? &( accum->ray_tri_tests_count ) : NULL );
813 }
814 
815 ErrorCode OrientedBoxTreeTool::ray_intersect_boxes( Range& boxes_out,
816  EntityHandle root_set,
817  double tolerance,
818  const double ray_point[3],
819  const double unit_ray_dir[3],
820  const double* ray_length,
821  TrvStats* accum )
822 {
823  RayIntersector op( this, ray_point, unit_ray_dir, ray_length, tolerance, boxes_out );
824  return preorder_traverse( root_set, op, accum );
825 }
826 
827 ErrorCode RayIntersector::visit( EntityHandle node, int, bool& descend )
828 {
829  OrientedBox box;
830  ErrorCode rval = tool->box( node, box );
831  if( MB_SUCCESS != rval ) return rval;
832 
833  descend = box.intersect_ray( b, m, tol, len );
834  return MB_SUCCESS;
835 }
836 
837 ErrorCode RayIntersector::leaf( EntityHandle node )
838 {
839  boxes.insert( node );
840  return MB_SUCCESS;
841 }
842 
843 /********************** Ray/Set Intersection ****************************/
844 
845 ErrorCode OrientedBoxTreeTool::get_close_tris( CartVect int_pt,
846  double tol,
847  const EntityHandle* rootSet,
848  const EntityHandle* geomVol,
849  const Tag* senseTag,
850  std::vector< EntityHandle >& close_tris,
851  std::vector< int >& close_senses )
852 {
853  std::vector< EntityHandle > close_surfs;
854  ErrorCode rval = sphere_intersect_triangles( int_pt.array(), tol, *rootSet, close_tris, &close_surfs );
855  assert( MB_SUCCESS == rval );
856  if( MB_SUCCESS != rval ) return rval;
857 
858  // for each surface, get the surf sense wrt parent volume
859  close_senses.resize( close_surfs.size() );
860  for( unsigned i = 0; i < close_surfs.size(); ++i )
861  {
862  EntityHandle vols[2];
863  rval = get_moab_instance()->tag_get_data( *senseTag, &( close_surfs[i] ), 1, vols );
864  assert( MB_SUCCESS == rval );
865  if( MB_SUCCESS != rval ) return rval;
866  if( vols[0] == vols[1] )
867  {
868  std::cerr << "error: surf has positive and negative sense wrt same volume" << std::endl;
869  return MB_FAILURE;
870  }
871  if( *geomVol == vols[0] )
872  {
873  close_senses[i] = 1;
874  }
875  else if( *geomVol == vols[1] )
876  {
877  close_senses[i] = -1;
878  }
879  else
880  {
881  return MB_FAILURE;
882  }
883  }
884 
885  return MB_SUCCESS;
886 }
887 
888 class RayIntersectSets : public OrientedBoxTreeTool::Op
889 {
890  private:
891  // Input
892  OrientedBoxTreeTool* tool;
893  const CartVect ray_origin;
894  const CartVect ray_direction;
895  OrientedBoxTreeTool::IntersectSearchWindow& search_win; /* length to search ahead/behind of ray origin */
896  const double tol; /* used for box.intersect_ray, radius of
897  neighborhood for adjacent triangles,
898  and old mode of add_intersection */
899  OrientedBoxTreeTool::IntRegCtxt& int_reg_callback;
900 
901  // Optional Input - to screen RTIs by orientation and edge/node intersection
902  int* surfTriOrient; /* holds desired orientation of tri wrt surface */
903  int surfTriOrient_val;
904 
905  // Other Variables
906  unsigned int* raytri_test_count;
907  EntityHandle lastSet;
908  int lastSetDepth;
909 
910  public:
911  RayIntersectSets( OrientedBoxTreeTool* tool_ptr,
912  const double* ray_point,
913  const double* unit_ray_dir,
914  const double tolerance,
915  OrientedBoxTreeTool::IntersectSearchWindow& win,
916  unsigned int* ray_tri_test_count,
917  OrientedBoxTreeTool::IntRegCtxt& intRegCallback )
918  : tool( tool_ptr ), ray_origin( ray_point ), ray_direction( unit_ray_dir ), search_win( win ), tol( tolerance ),
919  int_reg_callback( intRegCallback ), surfTriOrient_val( 0 ), raytri_test_count( ray_tri_test_count ),
920  lastSet( 0 ), lastSetDepth( 0 )
921  {
922 
923  // specified orientation should be 1 or -1, indicating ray and surface
924  // normal in the same or opposite directions, respectively.
925  if( int_reg_callback.getDesiredOrient() )
926  {
927  surfTriOrient = &surfTriOrient_val;
928  }
929  else
930  {
931  surfTriOrient = NULL;
932  }
933 
934  // check the limits
935  if( search_win.first )
936  {
937  assert( 0 <= *( search_win.first ) );
938  }
939  if( search_win.second )
940  {
941  assert( 0 >= *( search_win.second ) );
942  }
943  };
944 
945  virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
946  virtual ErrorCode leaf( EntityHandle node );
947 };
948 
949 ErrorCode RayIntersectSets::visit( EntityHandle node, int depth, bool& descend )
950 {
951  OrientedBox box;
952  ErrorCode rval = tool->box( node, box );
953  assert( MB_SUCCESS == rval );
954  if( MB_SUCCESS != rval ) return rval;
955 
956  descend = box.intersect_ray( ray_origin, ray_direction, tol, search_win.first, search_win.second );
957 
958  if( lastSet && depth <= lastSetDepth ) lastSet = 0;
959 
960  if( descend && !lastSet )
961  {
962  Range tmp_sets;
963  rval = tool->get_moab_instance()->get_entities_by_type( node, MBENTITYSET, tmp_sets );
964  assert( MB_SUCCESS == rval );
965  if( MB_SUCCESS != rval ) return rval;
966 
967  if( !tmp_sets.empty() )
968  {
969  if( tmp_sets.size() > 1 ) return MB_FAILURE;
970  lastSet = *tmp_sets.begin();
971  lastSetDepth = depth;
972 
973  rval = int_reg_callback.update_orient( lastSet, surfTriOrient );
974  if( MB_SUCCESS != rval ) return rval;
975  }
976  }
977 
978  return MB_SUCCESS;
979 }
980 
981 ErrorCode RayIntersectSets::leaf( EntityHandle node )
982 {
983  assert( lastSet );
984  if( !lastSet ) // if no surface has been visited yet, something's messed up.
985  return MB_FAILURE;
986 
987 #ifndef MB_OBB_USE_VECTOR_QUERIES
988  Range tris;
989 #ifdef MB_OBB_USE_TYPE_QUERIES
990  ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
991 #else
992  ErrorCode rval = tool->get_moab_instance()->get_entities_by_handle( node, tris );
993 #endif
994 #else
995  std::vector< EntityHandle > tris;
996  ErrorCode rval = tool->get_moab_instance()->get_entities_by_handle( node, tris );
997 #endif
998  assert( MB_SUCCESS == rval );
999  if( MB_SUCCESS != rval ) return rval;
1000 
1001 #ifndef MB_OBB_USE_VECTOR_QUERIES
1002  for( Range::iterator t = tris.begin(); t != tris.end(); ++t )
1003 #else
1004  for( std::vector< EntityHandle >::iterator t = tris.begin(); t != tris.end(); ++t )
1005 #endif
1006  {
1007 #ifndef MB_OBB_USE_TYPE_QUERIES
1008  if( TYPE_FROM_HANDLE( *t ) != MBTRI ) continue;
1009 #endif
1010 
1011  const EntityHandle* conn;
1012  int num_conn;
1013  rval = tool->get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
1014  assert( MB_SUCCESS == rval );
1015  if( MB_SUCCESS != rval ) return rval;
1016 
1017  CartVect coords[3];
1018  rval = tool->get_moab_instance()->get_coords( conn, 3, coords[0].array() );
1019  assert( MB_SUCCESS == rval );
1020  if( MB_SUCCESS != rval ) return rval;
1021 
1022  if( raytri_test_count ) *raytri_test_count += 1;
1023 
1024  double int_dist;
1025  GeomUtil::intersection_type int_type = GeomUtil::NONE;
1026 
1027  if( GeomUtil::plucker_ray_tri_intersect( coords, ray_origin, ray_direction, int_dist, search_win.first,
1028  search_win.second, surfTriOrient, &int_type ) )
1029  {
1030  int_reg_callback.register_intersection( lastSet, *t, int_dist, search_win, int_type );
1031  }
1032  }
1033  return MB_SUCCESS;
1034 }
1035 
1036 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
1037  std::vector< EntityHandle >& sets_out,
1038  std::vector< EntityHandle >& facets_out,
1039  EntityHandle root_set,
1040  const double tolerance,
1041  const double ray_point[3],
1042  const double unit_ray_dir[3],
1043  IntersectSearchWindow& search_win,
1044  IntRegCtxt& int_reg_callback,
1045  TrvStats* accum )
1046 
1047 {
1048  RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
1049  accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
1050  ErrorCode rval = preorder_traverse( root_set, op, accum );
1051 
1052  distances_out = int_reg_callback.get_intersections();
1053  sets_out = int_reg_callback.get_sets();
1054  facets_out = int_reg_callback.get_facets();
1055 
1056  return rval;
1057 }
1058 
1059 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( std::vector< double >& distances_out,
1060  std::vector< EntityHandle >& sets_out,
1061  std::vector< EntityHandle >& facets_out,
1062  EntityHandle root_set,
1063  double tolerance,
1064  const double ray_point[3],
1065  const double unit_ray_dir[3],
1066  const double* ray_length,
1067  TrvStats* accum )
1068 {
1069  IntRegCtxt int_reg_ctxt;
1070 
1071  OrientedBoxTreeTool::IntersectSearchWindow search_win( ray_length, (double*)0 );
1072 
1073  RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
1074  accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_ctxt );
1075 
1076  ErrorCode rval = preorder_traverse( root_set, op, accum );
1077 
1078  if( MB_SUCCESS != rval )
1079  {
1080  return rval;
1081  }
1082 
1083  distances_out = int_reg_ctxt.get_intersections();
1084  sets_out = int_reg_ctxt.get_sets();
1085  facets_out = int_reg_ctxt.get_facets();
1086 
1087  return MB_SUCCESS;
1088 }
1089 
1090 ErrorCode OrientedBoxTreeTool::ray_intersect_sets( EntityHandle root_set,
1091  const double tolerance,
1092  const double ray_point[3],
1093  const double unit_ray_dir[3],
1094  IntersectSearchWindow& search_win,
1095  IntRegCtxt& int_reg_callback,
1096  TrvStats* accum )
1097 
1098 {
1099  RayIntersectSets op( this, ray_point, unit_ray_dir, tolerance, search_win,
1100  accum ? &( accum->ray_tri_tests_count ) : NULL, int_reg_callback );
1101  return preorder_traverse( root_set, op, accum );
1102 }
1103 
1104 /********************** Closest Point code ***************/
1105 
1106 struct OBBTreeCPFrame
1107 {
1108  OBBTreeCPFrame( double d, EntityHandle n, EntityHandle s, int dp )
1109  : dist_sqr( d ), node( n ), mset( s ), depth( dp )
1110  {
1111  }
1112  double dist_sqr;
1113  EntityHandle node;
1114  EntityHandle mset;
1115  int depth;
1116 };
1117 
1118 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point,
1119  EntityHandle root,
1120  double* point_out,
1121  EntityHandle& facet_out,
1122  EntityHandle* set_out,
1123  TrvStats* accum )
1124 {
1125  ErrorCode rval;
1126  const CartVect loc( point );
1127  double smallest_dist_sqr = std::numeric_limits< double >::max();
1128 
1129  EntityHandle current_set = 0;
1130  Range sets;
1131  std::vector< EntityHandle > children( 2 );
1132  std::vector< double > coords;
1133  std::vector< OBBTreeCPFrame > stack;
1134  int max_depth = -1;
1135 
1136  stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
1137 
1138  while( !stack.empty() )
1139  {
1140 
1141  // pop from top of stack
1142  EntityHandle node = stack.back().node;
1143  double dist_sqr = stack.back().dist_sqr;
1144  current_set = stack.back().mset;
1145  int current_depth = stack.back().depth;
1146  stack.pop_back();
1147 
1148  // If current best result is closer than the box, skip this tree node.
1149  if( dist_sqr > smallest_dist_sqr ) continue;
1150 
1151  // increment traversal statistics.
1152  if( accum )
1153  {
1154  accum->increment( current_depth );
1155  max_depth = std::max( max_depth, current_depth );
1156  }
1157 
1158  // Check if this node has a set associated with it
1159  if( set_out && !current_set )
1160  {
1161  sets.clear();
1162  rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
1163  if( MB_SUCCESS != rval ) return rval;
1164  if( !sets.empty() )
1165  {
1166  if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
1167  current_set = sets.front();
1168  }
1169  }
1170 
1171  // Get child boxes
1172  children.clear();
1173  rval = instance->get_child_meshsets( node, children );
1174  if( MB_SUCCESS != rval ) return rval;
1175 
1176  // if not a leaf node
1177  if( !children.empty() )
1178  {
1179  if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
1180 
1181  // get boxes
1182  OrientedBox box1, box2;
1183  rval = box( children[0], box1 );
1184  if( MB_SUCCESS != rval ) return rval;
1185  rval = box( children[1], box2 );
1186  if( MB_SUCCESS != rval ) return rval;
1187 
1188  // get distance from each box
1189  CartVect pt1, pt2;
1190  box1.closest_location_in_box( loc, pt1 );
1191  box2.closest_location_in_box( loc, pt2 );
1192  pt1 -= loc;
1193  pt2 -= loc;
1194  const double dsqr1 = pt1 % pt1;
1195  const double dsqr2 = pt2 % pt2;
1196 
1197  // push children on tree such that closer one is on top
1198  if( dsqr1 < dsqr2 )
1199  {
1200  stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1201  stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1202  }
1203  else
1204  {
1205  stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1206  stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1207  }
1208  }
1209  else
1210  { // LEAF NODE
1211  if( accum )
1212  {
1213  accum->increment_leaf( current_depth );
1214  }
1215 
1216  Range facets;
1217  rval = instance->get_entities_by_dimension( node, 2, facets );
1218  if( MB_SUCCESS != rval ) return rval;
1219 
1220  const EntityHandle* conn = NULL;
1221  int len = 0;
1222  CartVect tmp, diff;
1223  for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
1224  {
1225  rval = instance->get_connectivity( *i, conn, len, true );
1226  if( MB_SUCCESS != rval ) return rval;
1227 
1228  coords.resize( 3 * len );
1229  rval = instance->get_coords( conn, len, &coords[0] );
1230  if( MB_SUCCESS != rval ) return rval;
1231 
1232  if( len == 3 )
1233  GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
1234  else
1235  GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
1236 
1237  diff = tmp - loc;
1238  dist_sqr = diff % diff;
1239  if( dist_sqr < smallest_dist_sqr )
1240  {
1241  smallest_dist_sqr = dist_sqr;
1242  facet_out = *i;
1243  tmp.get( point_out );
1244  if( set_out ) *set_out = current_set;
1245  }
1246  }
1247  } // LEAF NODE
1248  }
1249 
1250  if( accum )
1251  {
1252  accum->end_traversal( max_depth );
1253  }
1254 
1255  return MB_SUCCESS;
1256 }
1257 
1258 ErrorCode OrientedBoxTreeTool::closest_to_location( const double* point,
1259  EntityHandle root,
1260  double tolerance,
1261  std::vector< EntityHandle >& facets_out,
1262  std::vector< EntityHandle >* sets_out,
1263  TrvStats* accum )
1264 {
1265  ErrorCode rval;
1266  const CartVect loc( point );
1267  double smallest_dist_sqr = std::numeric_limits< double >::max();
1268  double smallest_dist = smallest_dist_sqr;
1269 
1270  EntityHandle current_set = 0;
1271  Range sets;
1272  std::vector< EntityHandle > children( 2 );
1273  std::vector< double > coords;
1274  std::vector< OBBTreeCPFrame > stack;
1275  int max_depth = -1;
1276 
1277  stack.push_back( OBBTreeCPFrame( 0.0, root, current_set, 0 ) );
1278 
1279  while( !stack.empty() )
1280  {
1281 
1282  // pop from top of stack
1283  EntityHandle node = stack.back().node;
1284  double dist_sqr = stack.back().dist_sqr;
1285  current_set = stack.back().mset;
1286  int current_depth = stack.back().depth;
1287  stack.pop_back();
1288 
1289  // If current best result is closer than the box, skip this tree node.
1290  if( dist_sqr > smallest_dist_sqr + tolerance ) continue;
1291 
1292  // increment traversal statistics.
1293  if( accum )
1294  {
1295  accum->increment( current_depth );
1296  max_depth = std::max( max_depth, current_depth );
1297  }
1298 
1299  // Check if this node has a set associated with it
1300  if( sets_out && !current_set )
1301  {
1302  sets.clear();
1303  rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
1304  if( MB_SUCCESS != rval ) return rval;
1305  if( !sets.empty() )
1306  {
1307  if( sets.size() != 1 ) return MB_MULTIPLE_ENTITIES_FOUND;
1308  current_set = *sets.begin();
1309  }
1310  }
1311 
1312  // Get child boxes
1313  children.clear();
1314  rval = instance->get_child_meshsets( node, children );
1315  if( MB_SUCCESS != rval ) return rval;
1316 
1317  // if not a leaf node
1318  if( !children.empty() )
1319  {
1320  if( children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
1321 
1322  // get boxes
1323  OrientedBox box1, box2;
1324  rval = box( children[0], box1 );
1325  if( MB_SUCCESS != rval ) return rval;
1326  rval = box( children[1], box2 );
1327  if( MB_SUCCESS != rval ) return rval;
1328 
1329  // get distance from each box
1330  CartVect pt1, pt2;
1331  box1.closest_location_in_box( loc, pt1 );
1332  box2.closest_location_in_box( loc, pt2 );
1333  pt1 -= loc;
1334  pt2 -= loc;
1335  const double dsqr1 = pt1 % pt1;
1336  const double dsqr2 = pt2 % pt2;
1337 
1338  // push children on tree such that closer one is on top
1339  if( dsqr1 < dsqr2 )
1340  {
1341  stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1342  stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1343  }
1344  else
1345  {
1346  stack.push_back( OBBTreeCPFrame( dsqr1, children[0], current_set, current_depth + 1 ) );
1347  stack.push_back( OBBTreeCPFrame( dsqr2, children[1], current_set, current_depth + 1 ) );
1348  }
1349  }
1350  else
1351  { // LEAF NODE
1352  if( accum )
1353  {
1354  accum->increment_leaf( current_depth );
1355  }
1356 
1357  Range facets;
1358  rval = instance->get_entities_by_dimension( node, 2, facets );
1359  if( MB_SUCCESS != rval ) return rval;
1360 
1361  const EntityHandle* conn = NULL;
1362  int len = 0;
1363  CartVect tmp, diff;
1364  for( Range::iterator i = facets.begin(); i != facets.end(); ++i )
1365  {
1366  rval = instance->get_connectivity( *i, conn, len, true );
1367  if( MB_SUCCESS != rval ) return rval;
1368 
1369  coords.resize( 3 * len );
1370  rval = instance->get_coords( conn, len, &coords[0] );
1371  if( MB_SUCCESS != rval ) return rval;
1372 
1373  if( len == 3 )
1374  GeomUtil::closest_location_on_tri( loc, (CartVect*)( &coords[0] ), tmp );
1375  else
1376  GeomUtil::closest_location_on_polygon( loc, (CartVect*)( &coords[0] ), len, tmp );
1377 
1378  diff = tmp - loc;
1379  dist_sqr = diff % diff;
1380  if( dist_sqr < smallest_dist_sqr )
1381  {
1382  if( 0.5 * dist_sqr < 0.5 * smallest_dist_sqr + tolerance * ( 0.5 * tolerance - smallest_dist ) )
1383  {
1384  facets_out.clear();
1385  if( sets_out ) sets_out->clear();
1386  }
1387  smallest_dist_sqr = dist_sqr;
1388  smallest_dist = sqrt( smallest_dist_sqr );
1389  // put closest result at start of lists
1390  facets_out.push_back( *i );
1391  std::swap( facets_out.front(), facets_out.back() );
1392  if( sets_out )
1393  {
1394  sets_out->push_back( current_set );
1395  std::swap( sets_out->front(), sets_out->back() );
1396  }
1397  }
1398  else if( dist_sqr <= smallest_dist_sqr + tolerance * ( tolerance + 2 * smallest_dist ) )
1399  {
1400  facets_out.push_back( *i );
1401  if( sets_out ) sets_out->push_back( current_set );
1402  }
1403  }
1404  } // LEAF NODE
1405  }
1406 
1407  if( accum )
1408  {
1409  accum->end_traversal( max_depth );
1410  }
1411 
1412  return MB_SUCCESS;
1413 }
1414 
1415 /********************** Tree Printing Code ****************************/
1416 
1417 class TreeLayoutPrinter : public OrientedBoxTreeTool::Op
1418 {
1419  public:
1420  TreeLayoutPrinter( std::ostream& stream, Interface* instance );
1421 
1422  virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
1423  virtual ErrorCode leaf( EntityHandle node );
1424 
1425  private:
1426  Interface* instance;
1427  std::ostream& outputStream;
1428  std::vector< bool > path;
1429 };
1430 
1431 TreeLayoutPrinter::TreeLayoutPrinter( std::ostream& stream, Interface* interface )
1432  : instance( interface ), outputStream( stream )
1433 {
1434 }
1435 
1436 ErrorCode TreeLayoutPrinter::visit( EntityHandle node, int depth, bool& descend )
1437 {
1438  descend = true;
1439 
1440  if( (unsigned)depth > path.size() )
1441  {
1442  // assert(depth+1 == path.size); // preorder traversal
1443  path.push_back( true );
1444  }
1445  else
1446  {
1447  path.resize( depth );
1448  if( depth ) path.back() = false;
1449  }
1450 
1451  for( unsigned i = 0; i + 1 < path.size(); ++i )
1452  {
1453  if( path[i] )
1454  outputStream << "| ";
1455  else
1456  outputStream << " ";
1457  }
1458  if( depth )
1459  {
1460  if( path.back() )
1461  outputStream << "+---";
1462  else
1463  outputStream << "\\---";
1464  }
1465  outputStream << instance->id_from_handle( node ) << std::endl;
1466  return MB_SUCCESS;
1467 }
1468 
1469 ErrorCode TreeLayoutPrinter::leaf( EntityHandle )
1470 {
1471  return MB_SUCCESS;
1472 }
1473 
1474 class TreeNodePrinter : public OrientedBoxTreeTool::Op
1475 {
1476  public:
1477  TreeNodePrinter( std::ostream& stream,
1478  bool list_contents,
1479  bool list_box,
1480  const char* id_tag_name,
1481  OrientedBoxTreeTool* tool_ptr );
1482 
1483  virtual ErrorCode visit( EntityHandle node, int depth, bool& descend );
1484 
1485  virtual ErrorCode leaf( EntityHandle )
1486  {
1487  return MB_SUCCESS;
1488  }
1489 
1490  private:
1491  ErrorCode print_geometry( EntityHandle node );
1492  ErrorCode print_contents( EntityHandle node );
1493  ErrorCode print_counts( EntityHandle node );
1494 
1495  bool printContents;
1496  bool printGeometry;
1497  bool haveTag;
1498  Tag tag, gidTag, geomTag;
1499  Interface* instance;
1500  OrientedBoxTreeTool* tool;
1501  std::ostream& outputStream;
1502 };
1503 
1504 TreeNodePrinter::TreeNodePrinter( std::ostream& stream,
1505  bool list_contents,
1506  bool list_box,
1507  const char* id_tag_name,
1508  OrientedBoxTreeTool* tool_ptr )
1509  : printContents( list_contents ), printGeometry( list_box ), haveTag( false ), tag( 0 ), gidTag( 0 ), geomTag( 0 ),
1510  instance( tool_ptr->get_moab_instance() ), tool( tool_ptr ), outputStream( stream )
1511 {
1512  ErrorCode rval;
1513  if( id_tag_name )
1514  {
1515  rval = instance->tag_get_handle( id_tag_name, 1, MB_TYPE_INTEGER, tag );
1516  if( !rval )
1517  {
1518  std::cerr << "Could not get tag \"" << id_tag_name << "\"\n";
1519  stream << "Could not get tag \"" << id_tag_name << "\"\n";
1520  }
1521  else
1522  {
1523  haveTag = true;
1524  }
1525  }
1526 
1527  gidTag = instance->globalId_tag();
1528 
1529  rval = instance->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );
1530  if( MB_SUCCESS != rval ) geomTag = 0;
1531 }
1532 
1533 ErrorCode TreeNodePrinter::visit( EntityHandle node, int, bool& descend )
1534 {
1535  descend = true;
1536  EntityHandle setid = instance->id_from_handle( node );
1537  outputStream << setid << ":" << std::endl;
1538 
1539  Range surfs;
1540  ErrorCode r3 = MB_SUCCESS;
1541  if( geomTag )
1542  {
1543  const int two = 2;
1544  const void* tagdata[] = { &two };
1545  r3 = instance->get_entities_by_type_and_tag( node, MBENTITYSET, &geomTag, tagdata, 1, surfs );
1546 
1547  if( MB_SUCCESS == r3 && surfs.size() == 1 )
1548  {
1549  EntityHandle surf = *surfs.begin();
1550  int id;
1551  if( gidTag && MB_SUCCESS == instance->tag_get_data( gidTag, &surf, 1, &id ) )
1552  outputStream << " Surface " << id << std::endl;
1553  else
1554  outputStream << " Surface w/ unknown ID (" << surf << ")" << std::endl;
1555  }
1556  }
1557 
1558  ErrorCode r1 = printGeometry ? print_geometry( node ) : MB_SUCCESS;
1559  ErrorCode r2 = printContents ? print_contents( node ) : print_counts( node );
1560  outputStream << std::endl;
1561 
1562  if( MB_SUCCESS != r1 )
1563  return r1;
1564  else if( MB_SUCCESS != r2 )
1565  return r2;
1566  else
1567  return r3;
1568 }
1569 
1570 ErrorCode TreeNodePrinter::print_geometry( EntityHandle node )
1571 {
1572  OrientedBox box;
1573  ErrorCode rval = tool->box( node, box );
1574  if( MB_SUCCESS != rval ) return rval;
1575 
1576  CartVect length = box.dimensions();
1577 
1578  outputStream << box.center << " Radius: " << box.inner_radius() << " - " << box.outer_radius() << std::endl
1579  << '+' << box.axis( 0 ) << " : " << length[0] << std::endl
1580  << 'x' << box.axis( 1 ) << " : " << length[1] << std::endl
1581  << 'x' << box.axis( 2 ) << " : " << length[2] << std::endl;
1582  return MB_SUCCESS;
1583 }
1584 
1585 ErrorCode TreeNodePrinter::print_counts( EntityHandle node )
1586 {
1587  for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
1588  {
1589  int count = 0;
1590  ErrorCode rval = instance->get_number_entities_by_type( node, type, count );
1591  if( MB_SUCCESS != rval ) return rval;
1592  if( count > 0 ) outputStream << " " << count << " " << CN::EntityTypeName( type ) << std::endl;
1593  }
1594  return MB_SUCCESS;
1595 }
1596 
1597 ErrorCode TreeNodePrinter::print_contents( EntityHandle node )
1598 {
1599  // list contents
1600  for( EntityType type = MBVERTEX; type != MBMAXTYPE; ++type )
1601  {
1602  Range range;
1603  ErrorCode rval = instance->get_entities_by_type( node, type, range );
1604  if( MB_SUCCESS != rval ) return rval;
1605  if( range.empty() ) continue;
1606  outputStream << " " << CN::EntityTypeName( type ) << " ";
1607  std::vector< int > ids( range.size() );
1608  if( haveTag )
1609  {
1610  rval = instance->tag_get_data( tag, range, &ids[0] );
1611  std::sort( ids.begin(), ids.end() );
1612  }
1613  else
1614  {
1615  Range::iterator ri = range.begin();
1616  std::vector< int >::iterator vi = ids.begin();
1617  while( ri != range.end() )
1618  {
1619  *vi = instance->id_from_handle( *ri );
1620  ++ri;
1621  ++vi;
1622  }
1623  }
1624 
1625  unsigned i = 0;
1626  for( ;; )
1627  {
1628  unsigned beg = i, end;
1629  do
1630  {
1631  end = i++;
1632  } while( i < ids.size() && ids[end] + 1 == ids[i] );
1633  if( end == beg )
1634  outputStream << ids[end];
1635  else if( end == beg + 1 )
1636  outputStream << ids[beg] << ", " << ids[end];
1637  else
1638  outputStream << ids[beg] << "-" << ids[end];
1639 
1640  if( i == ids.size() )
1641  {
1642  outputStream << std::endl;
1643  break;
1644  }
1645  else
1646  outputStream << ", ";
1647  }
1648  }
1649 
1650  return MB_SUCCESS;
1651 }
1652 
1653 void OrientedBoxTreeTool::print( EntityHandle set, std::ostream& str, bool list, const char* tag )
1654 {
1655  TreeLayoutPrinter op1( str, instance );
1656  TreeNodePrinter op2( str, list, true, tag, this );
1657  ErrorCode r1 = preorder_traverse( set, op1 );
1658  str << std::endl;
1659  ErrorCode r2 = preorder_traverse( set, op2 );
1660  if( r1 != MB_SUCCESS || r2 != MB_SUCCESS )
1661  {
1662  std::cerr << "Errors encountered while printing tree\n";
1663  str << "Errors encountered while printing tree\n";
1664  }
1665 }
1666 
1667 /********************* Traversal Metrics Code **************************/
1668 
1669 void OrientedBoxTreeTool::TrvStats::reset()
1670 {
1671  nodes_visited_count.clear();
1672  leaves_visited_count.clear();
1673  traversals_ended_count.clear();
1674  ray_tri_tests_count = 0;
1675 }
1676 
1677 void OrientedBoxTreeTool::TrvStats::increment( unsigned depth )
1678 {
1679 
1680  while( nodes_visited_count.size() <= depth )
1681  {
1682  nodes_visited_count.push_back( 0 );
1683  leaves_visited_count.push_back( 0 );
1684  traversals_ended_count.push_back( 0 );
1685  }
1686  nodes_visited_count[depth] += 1;
1687 }
1688 
1689 void OrientedBoxTreeTool::TrvStats::increment_leaf( unsigned depth )
1690 {
1691  // assume safe depth, because increment is called first
1692  leaves_visited_count[depth] += 1;
1693 }
1694 
1695 void OrientedBoxTreeTool::TrvStats::end_traversal( unsigned depth )
1696 {
1697  // assume safe depth, because increment is always called on a given
1698  // tree level first
1699  traversals_ended_count[depth] += 1;
1700 }
1701 
1702 void OrientedBoxTreeTool::TrvStats::print( std::ostream& str ) const
1703 {
1704 
1705  const std::string h1 = "OBBTree Depth";
1706  const std::string h2 = " - NodesVisited";
1707  const std::string h3 = " - LeavesVisited";
1708  const std::string h4 = " - TraversalsEnded";
1709 
1710  str << h1 << h2 << h3 << h4 << std::endl;
1711 
1712  unsigned num_visited = 0, num_leaves = 0, num_traversals = 0;
1713  for( unsigned i = 0; i < traversals_ended_count.size(); ++i )
1714  {
1715 
1716  num_visited += nodes_visited_count[i];
1717  num_leaves += leaves_visited_count[i];
1718  num_traversals += traversals_ended_count[i];
1719 
1720  str << std::setw( h1.length() ) << i << std::setw( h2.length() ) << nodes_visited_count[i]
1721  << std::setw( h3.length() ) << leaves_visited_count[i] << std::setw( h4.length() )
1722  << traversals_ended_count[i] << std::endl;
1723  }
1724 
1725  str << std::setw( h1.length() ) << "---- Totals:" << std::setw( h2.length() ) << num_visited
1726  << std::setw( h3.length() ) << num_leaves << std::setw( h4.length() ) << num_traversals << std::endl;
1727 
1728  if( ray_tri_tests_count )
1729  {
1730  str << std::setw( h1.length() ) << "---- Total ray-tri tests: " << ray_tri_tests_count << std::endl;
1731  }
1732 }
1733 
1734 /********************** Tree Statistics Code ****************************/
1735 
1736 class StatData
1737 {
1738  public:
1739  struct Ratio
1740  {
1741  double min, max, sum, sqr;
1742  int hist[10];
1743  Ratio()
1744  : min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() ), sum( 0.0 ),
1745  sqr( 0.0 )
1746  {
1747  hist[0] = hist[1] = hist[2] = hist[3] = hist[4] = hist[5] = hist[6] = hist[7] = hist[8] = hist[9] = 0;
1748  }
1749  void accum( double v )
1750  {
1751  if( v < min ) min = v;
1752  if( v > max ) max = v;
1753  sum += v;
1754  sqr += v * v;
1755  int i = (int)( 10 * v );
1756  if( i < 0 )
1757  i = 0;
1758  else if( i > 9 )
1759  i = 9;
1760  ++hist[i];
1761  }
1762  };
1763 
1764  template < typename T >
1765  struct Stat
1766  {
1767  T min, max;
1768  double sum, sqr;
1769  Stat() : sum( 0.0 ), sqr( 0.0 )
1770  {
1771  std::numeric_limits< T > lim;
1772  min = lim.max();
1773  if( lim.is_integer )
1774  max = lim.min();
1775  else
1776  max = -lim.max();
1777  }
1778  void accum( T v )
1779  {
1780  if( v < min ) min = v;
1781  if( v > max ) max = v;
1782  sum += v;
1783  sqr += (double)v * v;
1784  }
1785  };
1786 
1787  StatData() : count( 0 ) {}
1788 
1789  Ratio volume;
1790  Ratio entities;
1791  Ratio radius;
1792  Stat< unsigned > leaf_ent;
1793  Stat< double > vol;
1794  Stat< double > area;
1795  std::vector< unsigned > leaf_depth;
1796  unsigned count;
1797 };
1798 
1799 static int measure( const CartVect& v, double& result )
1800 {
1801  const double tol = 1e-6;
1802  int dims = 0;
1803  result = 1;
1804  for( int i = 0; i < 3; ++i )
1805  if( v[i] > tol )
1806  {
1807  ++dims;
1808  result *= v[i];
1809  }
1810  return dims;
1811 }
1812 
1813 ErrorCode OrientedBoxTreeTool::recursive_stats( OrientedBoxTreeTool* tool,
1814  Interface* inst,
1815  EntityHandle set,
1816  int depth,
1817  StatData& data,
1818  unsigned& count_out,
1819  CartVect& dimensions_out )
1820 {
1821  ErrorCode rval;
1822  OrientedBox tmp_box;
1823  std::vector< EntityHandle > children( 2 );
1824  unsigned counts[2];
1825  bool isleaf;
1826 
1827  ++data.count;
1828 
1829  rval = tool->box( set, tmp_box );
1830  if( MB_SUCCESS != rval ) return rval;
1831  children.clear();
1832  rval = inst->get_child_meshsets( set, children );
1833  if( MB_SUCCESS != rval ) return rval;
1834  isleaf = children.empty();
1835  if( !isleaf && children.size() != 2 ) return MB_MULTIPLE_ENTITIES_FOUND;
1836 
1837  dimensions_out = tmp_box.dimensions();
1838  data.radius.accum( tmp_box.inner_radius() / tmp_box.outer_radius() );
1839  data.vol.accum( tmp_box.volume() );
1840  data.area.accum( tmp_box.area() );
1841 
1842  if( isleaf )
1843  {
1844  if( data.leaf_depth.size() <= (unsigned)depth ) data.leaf_depth.resize( depth + 1, 0 );
1845  ++data.leaf_depth[depth];
1846 
1847  int count;
1848  rval = inst->get_number_entities_by_handle( set, count );
1849  if( MB_SUCCESS != rval ) return rval;
1850  count_out = count;
1851  data.leaf_ent.accum( count_out );
1852  }
1853  else
1854  {
1855  for( int i = 0; i < 2; ++i )
1856  {
1857  CartVect dims;
1858  rval = recursive_stats( tool, inst, children[i], depth + 1, data, counts[i], dims );
1859  if( MB_SUCCESS != rval ) return rval;
1860  double this_measure, chld_measure;
1861  int this_dim = measure( dimensions_out, this_measure );
1862  int chld_dim = measure( dims, chld_measure );
1863  double ratio;
1864  if( chld_dim < this_dim )
1865  ratio = 0;
1866  else
1867  ratio = chld_measure / this_measure;
1868 
1869  data.volume.accum( ratio );
1870  }
1871  count_out = counts[0] + counts[1];
1872  data.entities.accum( (double)counts[0] / count_out );
1873  data.entities.accum( (double)counts[1] / count_out );
1874  }
1875  return MB_SUCCESS;
1876 }
1877 
1878 static inline double std_dev( double sqr, double sum, double count )
1879 {
1880  sum /= count;
1881  sqr /= count;
1882  return sqrt( sqr - sum * sum );
1883 }
1884 
1885 //#define WW <<std::setw(10)<<std::fixed<<
1886 #define WE << std::setw( 10 ) <<
1887 #define WW WE
1888 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set,
1889  unsigned& total_entities,
1890  double& rv,
1891  double& tot_node_volume,
1892  double& tot_to_root_volume,
1893  unsigned& tree_height,
1894  unsigned& node_count,
1895  unsigned& num_leaves )
1896 {
1897  StatData d;
1898  ErrorCode rval;
1899  unsigned i;
1900  CartVect total_dim;
1901 
1902  rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
1903  if( MB_SUCCESS != rval ) return rval;
1904 
1905  tree_height = d.leaf_depth.size();
1906  unsigned min_leaf_depth = tree_height;
1907  num_leaves = 0;
1908  unsigned max_leaf_per_depth = 0;
1909  // double sum_leaf_depth = 0, sqr_leaf_depth = 0;
1910  for( i = 0; i < d.leaf_depth.size(); ++i )
1911  {
1912  unsigned val = d.leaf_depth[i];
1913  num_leaves += val;
1914  // sum_leaf_depth += (double)val*i;
1915  // sqr_leaf_depth += (double)val*i*i;
1916  if( val && i < min_leaf_depth ) min_leaf_depth = i;
1917  if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
1918  }
1919  rv = total_dim[0] * total_dim[1] * total_dim[2];
1920  tot_node_volume = d.vol.sum;
1921  tot_to_root_volume = d.vol.sum / rv;
1922  node_count = d.count;
1923 
1924  return MB_SUCCESS;
1925 }
1926 
1927 ErrorCode OrientedBoxTreeTool::stats( EntityHandle set, std::ostream& s )
1928 {
1929  StatData d;
1930  ErrorCode rval;
1931  unsigned total_entities, i;
1932  CartVect total_dim;
1933 
1934  rval = recursive_stats( this, instance, set, 0, d, total_entities, total_dim );
1935  if( MB_SUCCESS != rval ) return rval;
1936 
1937  unsigned tree_height = d.leaf_depth.size();
1938  unsigned min_leaf_depth = tree_height, num_leaves = 0;
1939  unsigned max_leaf_per_depth = 0;
1940  double sum_leaf_depth = 0, sqr_leaf_depth = 0;
1941  for( i = 0; i < d.leaf_depth.size(); ++i )
1942  {
1943  unsigned val = d.leaf_depth[i];
1944  num_leaves += val;
1945  sum_leaf_depth += (double)val * i;
1946  sqr_leaf_depth += (double)val * i * i;
1947  if( val && i < min_leaf_depth ) min_leaf_depth = i;
1948  if( max_leaf_per_depth < val ) max_leaf_per_depth = val;
1949  }
1950  unsigned num_non_leaf = d.count - num_leaves;
1951 
1952  double rv = total_dim[0] * total_dim[1] * total_dim[2];
1953  s << "entities in tree: " << total_entities << std::endl
1954  << "root volume: " << rv << std::endl
1955  << "total node volume: " << d.vol.sum << std::endl
1956  << "total/root volume: " << d.vol.sum / rv << std::endl
1957  << "tree height: " << tree_height << std::endl
1958  << "node count: " << d.count << std::endl
1959  << "leaf count: " << num_leaves << std::endl
1960  << std::endl;
1961 
1962  double avg_leaf_depth = sum_leaf_depth / num_leaves;
1963  double rms_leaf_depth = sqrt( sqr_leaf_depth / num_leaves );
1964  double std_leaf_depth = std_dev( sqr_leaf_depth, sum_leaf_depth, num_leaves );
1965 
1966  double avg_leaf_ent = d.leaf_ent.sum / num_leaves;
1967  double rms_leaf_ent = sqrt( d.leaf_ent.sqr / num_leaves );
1968  double std_leaf_ent = std_dev( d.leaf_ent.sqr, d.leaf_ent.sum, num_leaves );
1969 
1970  unsigned num_child = 2 * num_non_leaf;
1971 
1972  double avg_vol_ratio = d.volume.sum / num_child;
1973  double rms_vol_ratio = sqrt( d.volume.sqr / num_child );
1974  double std_vol_ratio = std_dev( d.volume.sqr, d.volume.sum, num_child );
1975 
1976  double avg_ent_ratio = d.entities.sum / num_child;
1977  double rms_ent_ratio = sqrt( d.entities.sqr / num_child );
1978  double std_ent_ratio = std_dev( d.entities.sqr, d.entities.sum, num_child );
1979 
1980  double avg_rad_ratio = d.radius.sum / d.count;
1981  double rms_rad_ratio = sqrt( d.radius.sqr / d.count );
1982  double std_rad_ratio = std_dev( d.radius.sqr, d.radius.sum, d.count );
1983 
1984  double avg_vol = d.vol.sum / d.count;
1985  double rms_vol = sqrt( d.vol.sqr / d.count );
1986  double std_vol = std_dev( d.vol.sqr, d.vol.sum, d.count );
1987 
1988  double avg_area = d.area.sum / d.count;
1989  double rms_area = sqrt( d.area.sqr / d.count );
1990  double std_area = std_dev( d.area.sqr, d.area.sum, d.count );
1991 
1992  int prec = s.precision();
1993  s << " " WW "Minimum" WW "Average" WW "RMS" WW "Maximum" WW "Std.Dev." << std::endl;
1994  s << std::setprecision( 1 )
1995  << "Leaf Depth " WW min_leaf_depth WW avg_leaf_depth WW rms_leaf_depth WW d.leaf_depth.size() -
1996  1 WW std_leaf_depth
1997  << std::endl;
1998  s << std::setprecision( 0 )
1999  << "Entities/Leaf " WW d.leaf_ent.min WW avg_leaf_ent WW rms_leaf_ent WW d.leaf_ent.max WW std_leaf_ent
2000  << std::endl;
2001  s << std::setprecision( 3 )
2002  << "Child Volume Ratio " WW d.volume.min WW avg_vol_ratio WW rms_vol_ratio WW d.volume.max WW std_vol_ratio
2003  << std::endl;
2004  s << std::setprecision( 3 )
2005  << "Child Entity Ratio " WW d.entities.min WW avg_ent_ratio WW rms_ent_ratio WW d.entities.max WW std_ent_ratio
2006  << std::endl;
2007  s << std::setprecision( 3 )
2008  << "Box Radius Ratio " WW d.radius.min WW avg_rad_ratio WW rms_rad_ratio WW d.radius.max WW std_rad_ratio
2009  << std::endl;
2010  s << std::setprecision( 0 ) << "Box volume " WE d.vol.min WE avg_vol WE rms_vol WE d.vol.max WE std_vol
2011  << std::endl;
2012  s << std::setprecision( 0 ) << "Largest side area " WE d.area.min WE avg_area WE rms_area WE d.area.max WE std_area
2013  << std::endl;
2014  s << std::setprecision( prec ) << std::endl;
2015 
2016  s << "Leaf Depth Histogram (Root depth is 0)" << std::endl;
2017  double f = 60.0 / max_leaf_per_depth;
2018  for( i = min_leaf_depth; i < d.leaf_depth.size(); ++i )
2019  s << std::setw( 2 ) << i << " " << std::setw( 5 ) << d.leaf_depth[i] << " |" << std::setfill( '*' )
2020  << std::setw( (int)floor( f * d.leaf_depth[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2021  s << std::endl;
2022 
2023  s << "Child/Parent Volume Ratio Histogram" << std::endl;
2024  f = 60.0 / *( std::max_element( d.volume.hist, d.volume.hist + 10 ) );
2025  for( i = 0; i < 10u; ++i )
2026  s << "0." << i << " " << std::setw( 5 ) << d.volume.hist[i] << " |" << std::setfill( '*' )
2027  << std::setw( (int)floor( f * d.volume.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2028  s << std::endl;
2029 
2030  s << "Child/Parent Entity Count Ratio Histogram" << std::endl;
2031  f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 10 ) );
2032  for( i = 0; i < 10u; ++i )
2033  s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
2034  << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2035  s << std::endl;
2036 
2037  s << "Inner/Outer Radius Ratio Histogram (~0.70 for cube)" << std::endl;
2038  // max radius ratio for a box is about 0.7071. Move any boxes
2039  // in the .7 bucket into .6 and print .0 to .6.
2040  d.radius.hist[6] += d.radius.hist[7];
2041  f = 60.0 / *( std::max_element( d.entities.hist, d.entities.hist + 7 ) );
2042  for( i = 0; i < 7u; ++i )
2043  s << "0." << i << " " << std::setw( 5 ) << d.entities.hist[i] << " |" << std::setfill( '*' )
2044  << std::setw( (int)floor( f * d.entities.hist[i] + 0.5 ) ) << "" << std::setfill( ' ' ) << std::endl;
2045  s << std::endl;
2046 
2047  return MB_SUCCESS;
2048 }
2049 
2050 } // namespace moab