Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
AdaptiveKDTree.cpp
Go to the documentation of this file.
1 /**\file AdaptiveKDTree.cpp
2  */
3 
4 #include "moab/AdaptiveKDTree.hpp"
5 #include "moab/Interface.hpp"
6 #include "moab/GeomUtil.hpp"
7 #include "moab/Range.hpp"
8 #include "moab/ElemEvaluator.hpp"
9 #include "moab/CpuTimer.hpp"
10 #include "Internals.hpp"
11 #include "moab/Util.hpp"
12 #include <cmath>
13 
14 #include <cassert>
15 #include <algorithm>
16 #include <limits>
17 #include <iostream>
18 #include <cstdio>
19 
20 namespace moab
21 {
22 
23 const char* AdaptiveKDTree::treeName = "AKDTree";
24 
25 #define MB_AD_KD_TREE_DEFAULT_TAG_NAME
26 
27 // If defined, use single tag for both axis and location of split plane
28 #define MB_AD_KD_TREE_USE_SINGLE_TAG
29 
30 // No effect if MB_AD_KD_TREE_USE_SINGLE_TAG is not defined.
31 // If defined, store plane axis as double so tag has consistent
32 // type (doubles for both location and axis). If not defined,
33 // store struct Plane as opaque.
34 #define MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG
35 
36 AdaptiveKDTree::AdaptiveKDTree( Interface* iface )
37  : Tree( iface ), planeTag( 0 ), axisTag( 0 ), splitsPerDir( 3 ), planeSet( SUBDIVISION_SNAP ), spherical( false ),
38  radius( 1.0 )
39 {
40  boxTagName = treeName;
41 
42  ErrorCode rval = init();
43  if( MB_SUCCESS != rval ) throw rval;
44 }
45 
46 AdaptiveKDTree::AdaptiveKDTree( Interface* iface,
47  const Range& entities,
48  EntityHandle* tree_root_set,
49  FileOptions* opts )
50  : Tree( iface ), planeTag( 0 ), axisTag( 0 ), splitsPerDir( 3 ), planeSet( SUBDIVISION_SNAP ), spherical( false ),
51  radius( 1.0 )
52 {
53  boxTagName = treeName;
54 
55  ErrorCode rval;
56  if( opts )
57  {
58  rval = parse_options( *opts );
59  if( MB_SUCCESS != rval ) throw rval;
60  }
61 
62  rval = init();
63  if( MB_SUCCESS != rval ) throw rval;
64 
65  rval = build_tree( entities, tree_root_set, opts );
66  if( MB_SUCCESS != rval ) throw rval;
67 }
68 
69 AdaptiveKDTree::~AdaptiveKDTree()
70 {
71  if( !cleanUp ) return;
72 
73  if( myRoot )
74  {
75  reset_tree();
76  myRoot = 0;
77  }
78 }
79 
80 ErrorCode AdaptiveKDTree::build_tree( const Range& entities, EntityHandle* tree_root_set, FileOptions* options )
81 {
82  ErrorCode rval;
83  CpuTimer cp;
84 
85  if( options )
86  {
87  rval = parse_options( *options );
88  if( MB_SUCCESS != rval ) return rval;
89 
90  if( !options->all_seen() ) return MB_FAILURE;
91  }
92 
93  // calculate bounding box of elements
94  BoundBox box;
95  rval = box.update( *moab(), entities, spherical, radius );
96  if( MB_SUCCESS != rval ) return rval;
97 
98  // create tree root
99  EntityHandle tmp_root;
100  if( !tree_root_set ) tree_root_set = &tmp_root;
101  rval = create_root( box.bMin.array(), box.bMax.array(), *tree_root_set );
102  if( MB_SUCCESS != rval ) return rval;
103  rval = moab()->add_entities( *tree_root_set, entities );
104  if( MB_SUCCESS != rval ) return rval;
105 
106  AdaptiveKDTreeIter iter;
107  iter.initialize( this, *tree_root_set, box.bMin.array(), box.bMax.array(), AdaptiveKDTreeIter::LEFT );
108 
109  std::vector< double > tmp_data;
110  std::vector< EntityHandle > tmp_data2;
111  for( ;; )
112  {
113 
114  int pcount;
115  rval = moab()->get_number_entities_by_handle( iter.handle(), pcount );
116  if( MB_SUCCESS != rval ) break;
117 
118  const size_t p_count = pcount;
119  Range best_left, best_right, best_both;
120  Plane best_plane = { HUGE_VAL, -1 };
121  if( (int)p_count > maxPerLeaf && (int)iter.depth() < maxDepth )
122  {
123  switch( planeSet )
124  {
125  case AdaptiveKDTree::SUBDIVISION:
126  rval = best_subdivision_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
127  minWidth );
128  break;
129  case AdaptiveKDTree::SUBDIVISION_SNAP:
130  rval = best_subdivision_snap_plane( splitsPerDir, iter, best_left, best_right, best_both,
131  best_plane, tmp_data, minWidth );
132  break;
133  case AdaptiveKDTree::VERTEX_MEDIAN:
134  rval = best_vertex_median_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
135  tmp_data, minWidth );
136  break;
137  case AdaptiveKDTree::VERTEX_SAMPLE:
138  rval = best_vertex_sample_plane( splitsPerDir, iter, best_left, best_right, best_both, best_plane,
139  tmp_data, tmp_data2, minWidth );
140  break;
141  default:
142  rval = MB_FAILURE;
143  }
144 
145  if( MB_SUCCESS != rval ) return rval;
146  }
147 
148  if( best_plane.norm >= 0 )
149  {
150  best_left.merge( best_both );
151  best_right.merge( best_both );
152  rval = split_leaf( iter, best_plane, best_left, best_right );
153  if( MB_SUCCESS != rval ) return rval;
154  }
155  else
156  {
157  rval = iter.step();
158  if( MB_ENTITY_NOT_FOUND == rval )
159  {
160  rval = treeStats.compute_stats( mbImpl, myRoot );
161  treeStats.initTime = cp.time_elapsed();
162  return rval; // at end
163  }
164  else if( MB_SUCCESS != rval )
165  break;
166  }
167  }
168 
169  reset_tree();
170 
171  treeStats.reset();
172 
173  return rval;
174 }
175 
176 ErrorCode AdaptiveKDTree::parse_options( FileOptions& opts )
177 {
178  ErrorCode rval = parse_common_options( opts );
179  if( MB_SUCCESS != rval ) return rval;
180 
181  // SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
182  int tmp_int;
183  rval = opts.get_int_option( "SPLITS_PER_DIR", tmp_int );
184  if( MB_SUCCESS == rval ) splitsPerDir = tmp_int;
185 
186  // PLANE_SET: method used to decide split planes; see CandidatePlaneSet enum (below)
187  // for possible values; default = 1 (SUBDIVISION_SNAP)
188  rval = opts.get_int_option( "PLANE_SET", tmp_int );
189  if( MB_SUCCESS == rval && ( tmp_int < SUBDIVISION || tmp_int > VERTEX_SAMPLE ) )
190  return MB_FAILURE;
191  else if( MB_ENTITY_NOT_FOUND == rval )
192  planeSet = SUBDIVISION;
193  else
194  planeSet = (CandidatePlaneSet)( tmp_int );
195 
196  rval = opts.get_toggle_option( "SPHERICAL", false, spherical );
197  if( MB_SUCCESS != rval ) spherical = false;
198 
199  double tmp = 1.0;
200  rval = opts.get_real_option( "RADIUS", tmp );
201  if( MB_SUCCESS != rval )
202  radius = 1.0;
203  else
204  radius = tmp;
205 
206  return MB_SUCCESS;
207 }
208 
209 ErrorCode AdaptiveKDTree::make_tag( Interface* iface,
210  std::string name,
211  TagType storage,
212  DataType type,
213  int count,
214  void* default_val,
215  Tag& tag_handle,
216  std::vector< Tag >& created_tags )
217 {
218  ErrorCode rval =
219  iface->tag_get_handle( name.c_str(), count, type, tag_handle, MB_TAG_CREAT | storage, default_val );
220 
221  if( MB_SUCCESS == rval )
222  {
223  if( std::find( created_tags.begin(), created_tags.end(), tag_handle ) == created_tags.end() )
224  created_tags.push_back( tag_handle );
225  }
226  else
227  {
228  while( !created_tags.empty() )
229  {
230  iface->tag_delete( created_tags.back() );
231  created_tags.pop_back();
232  }
233 
234  planeTag = axisTag = (Tag)-1;
235  }
236 
237  return rval;
238 }
239 
240 ErrorCode AdaptiveKDTree::init()
241 {
242  std::vector< Tag > ctl;
243 
244 #ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
245  // create two tags, one for axis direction and one for axis coordinate
246  std::string n1( treeName ), n2( treeName );
247  n1 += "_coord";
248  n2 += "_norm";
249  ErrorCode rval = make_tag( moab(), n1, MB_TAG_DENSE, MB_TYPE_DOUBLE, 1, 0, planeTag, ctl );
250  if( MB_SUCCESS != rval ) return rval;
251  rval = make_tag( moab(), n2, MB_TAG_DENSE, MB_TYPE_INT, 1, 0, axisTag, ctl );
252  if( MB_SUCCESS != rval ) return rval;
253 
254 #elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
255  // create tag to hold two doubles, one for location and one for axis
256  std::string double_tag_name = std::string( treeName ) + std::string( "_coord_norm" );
257  ErrorCode rval = make_tag( moab(), double_tag_name, MB_TAG_DENSE, MB_TYPE_DOUBLE, 2, 0, planeTag, ctl );
258  if( MB_SUCCESS != rval ) return rval;
259 #else
260  // create opaque tag to hold struct Plane
261  ErrorCode rval = make_tag( moab(), tagname, MB_TAG_DENSE, MB_TYPE_OPAQUE, sizeof( Plane ), 0, planeTag, ctl );
262  if( MB_SUCCESS != rval ) return rval;
263 
264 #ifdef MOAB_HAVE_HDF5
265  // create a mesh tag holding the HDF5 type for a struct Plane
266  Tag type_tag;
267  std::string type_tag_name = "__hdf5_tag_type_";
268  type_tag_name += boxTagName;
269  rval = make_tag( moab(), type_tag_name, MB_TAG_MESH, MB_TYPE_OPAQUE, sizeof( hid_t ), 0, type_tag, ctl );
270  if( MB_SUCCESS != rval ) return rval;
271  // create HDF5 type object describing struct Plane
272  Plane p;
273  hid_t handle = H5Tcreate( H5T_COMPOUND, sizeof( Plane ) );
274  H5Tinsert( handle, "coord", &( p.coord ) - &p, H5T_NATIVE_DOUBLE );
275  H5Tinsert( handle, "norm", &( p.axis ) - &p, H5T_NATIVE_INT );
276  EntityHandle root = 0;
277  rval = mbImpl->tag_set_data( type_tag, &root, 1, &handle );
278  if( MB_SUCCESS != rval ) return rval;
279 #endif
280 #endif
281 
282  return rval;
283 }
284 
285 ErrorCode AdaptiveKDTree::get_split_plane( EntityHandle entity, Plane& plane )
286 {
287 #ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
288  ErrorCode r1, r2;
289  r1 = moab()->tag_get_data( planeTag, &entity, 1, &plane.coord );
290  r2 = moab()->tag_get_data( axisTag, &entity, 1, &plane.norm );
291  return MB_SUCCESS == r1 ? r2 : r1;
292 #elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
293  double values[2];
294  ErrorCode rval = moab()->tag_get_data( planeTag, &entity, 1, values );
295  plane.coord = values[0];
296  plane.norm = (int)values[1];
297  return rval;
298 #else
299  return moab()->tag_get_data( planeTag, &entity, 1, &plane );
300 #endif
301 }
302 
303 ErrorCode AdaptiveKDTree::set_split_plane( EntityHandle entity, const Plane& plane )
304 {
305 #ifndef MB_AD_KD_TREE_USE_SINGLE_TAG
306  ErrorCode r1, r2;
307  r1 = moab()->tag_set_data( planeTag, &entity, 1, &plane.coord );
308  r2 = moab()->tag_set_data( axisTag, &entity, 1, &plane.norm );
309  return MB_SUCCESS == r1 ? r2 : r1;
310 #elif defined( MB_AD_KD_TREE_USE_TWO_DOUBLE_TAG )
311  double values[2] = { plane.coord, static_cast< double >( plane.norm ) };
312  return moab()->tag_set_data( planeTag, &entity, 1, values );
313 #else
314  return moab()->tag_set_data( planeTag, &entity, 1, &plane );
315 #endif
316 }
317 
318 ErrorCode AdaptiveKDTree::get_tree_iterator( EntityHandle root, AdaptiveKDTreeIter& iter )
319 {
320  double box[6];
321  ErrorCode rval = moab()->tag_get_data( boxTag, &root, 1, box );
322  if( MB_SUCCESS != rval ) return rval;
323 
324  return get_sub_tree_iterator( root, box, box + 3, iter );
325 }
326 
327 ErrorCode AdaptiveKDTree::get_last_iterator( EntityHandle root, AdaptiveKDTreeIter& iter )
328 {
329  double box[6];
330  ErrorCode rval = moab()->tag_get_data( boxTag, &root, 1, box );
331  if( MB_SUCCESS != rval ) return rval;
332 
333  return iter.initialize( this, root, box, box + 3, AdaptiveKDTreeIter::RIGHT );
334 }
335 
336 ErrorCode AdaptiveKDTree::get_sub_tree_iterator( EntityHandle root,
337  const double min[3],
338  const double max[3],
339  AdaptiveKDTreeIter& result )
340 {
341  return result.initialize( this, root, min, max, AdaptiveKDTreeIter::LEFT );
342 }
343 
344 ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf, Plane plane, EntityHandle& left, EntityHandle& right )
345 {
346  ErrorCode rval;
347 
348  rval = moab()->create_meshset( meshsetFlags, left );
349  if( MB_SUCCESS != rval ) return rval;
350 
351  rval = moab()->create_meshset( meshsetFlags, right );
352  if( MB_SUCCESS != rval )
353  {
354  moab()->delete_entities( &left, 1 );
355  return rval;
356  }
357 
358  if( MB_SUCCESS != set_split_plane( leaf.handle(), plane ) ||
359  MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), left ) ||
360  MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), right ) ||
361  MB_SUCCESS != leaf.step_to_first_leaf( AdaptiveKDTreeIter::LEFT ) )
362  {
363  EntityHandle children[] = { left, right };
364  moab()->delete_entities( children, 2 );
365  return MB_FAILURE;
366  }
367 
368  return MB_SUCCESS;
369 }
370 
371 ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf, Plane plane )
372 {
373  EntityHandle left, right;
374  return split_leaf( leaf, plane, left, right );
375 }
376 
377 ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf,
378  Plane plane,
379  const Range& left_entities,
380  const Range& right_entities )
381 {
382  EntityHandle left, right, parent = leaf.handle();
383  ErrorCode rval = split_leaf( leaf, plane, left, right );
384  if( MB_SUCCESS != rval ) return rval;
385 
386  if( MB_SUCCESS == moab()->add_entities( left, left_entities ) &&
387  MB_SUCCESS == moab()->add_entities( right, right_entities ) &&
388  MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
389  return MB_SUCCESS;
390 
391  moab()->remove_child_meshset( parent, left );
392  moab()->remove_child_meshset( parent, right );
393  EntityHandle children[] = { left, right };
394  moab()->delete_entities( children, 2 );
395  return MB_FAILURE;
396 }
397 
398 ErrorCode AdaptiveKDTree::split_leaf( AdaptiveKDTreeIter& leaf,
399  Plane plane,
400  const std::vector< EntityHandle >& left_entities,
401  const std::vector< EntityHandle >& right_entities )
402 {
403  EntityHandle left, right, parent = leaf.handle();
404  ErrorCode rval = split_leaf( leaf, plane, left, right );
405  if( MB_SUCCESS != rval ) return rval;
406 
407  if( MB_SUCCESS == moab()->add_entities( left, &left_entities[0], left_entities.size() ) &&
408  MB_SUCCESS == moab()->add_entities( right, &right_entities[0], right_entities.size() ) &&
409  MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
410  return MB_SUCCESS;
411 
412  moab()->remove_child_meshset( parent, left );
413  moab()->remove_child_meshset( parent, right );
414  EntityHandle children[] = { left, right };
415  moab()->delete_entities( children, 2 );
416  return MB_FAILURE;
417 }
418 
419 ErrorCode AdaptiveKDTree::merge_leaf( AdaptiveKDTreeIter& iter )
420 {
421  ErrorCode rval;
422  if( iter.depth() == 1 ) // at root
423  return MB_FAILURE;
424 
425  // Move iter to parent
426 
427  AdaptiveKDTreeIter::StackObj node = iter.mStack.back();
428  iter.mStack.pop_back();
429 
430  iter.childVect.clear();
431  rval = moab()->get_child_meshsets( iter.mStack.back().entity, iter.childVect );
432  if( MB_SUCCESS != rval ) return rval;
433  Plane plane;
434  rval = get_split_plane( iter.mStack.back().entity, plane );
435  if( MB_SUCCESS != rval ) return rval;
436 
437  int child_idx = iter.childVect[0] == node.entity ? 0 : 1;
438  assert( iter.childVect[child_idx] == node.entity );
439  iter.mBox[1 - child_idx][plane.norm] = node.coord;
440 
441  // Get all entities from children and put them in parent
442  EntityHandle parent = iter.handle();
443  moab()->remove_child_meshset( parent, iter.childVect[0] );
444  moab()->remove_child_meshset( parent, iter.childVect[1] );
445  std::vector< EntityHandle > stack( iter.childVect );
446 
447  Range range;
448  while( !stack.empty() )
449  {
450  EntityHandle h = stack.back();
451  stack.pop_back();
452  range.clear();
453  rval = moab()->get_entities_by_handle( h, range );
454  if( MB_SUCCESS != rval ) return rval;
455  rval = moab()->add_entities( parent, range );
456  if( MB_SUCCESS != rval ) return rval;
457 
458  iter.childVect.clear();
459  rval = moab()->get_child_meshsets( h, iter.childVect );MB_CHK_ERR( rval );
460  if( !iter.childVect.empty() )
461  {
462  moab()->remove_child_meshset( h, iter.childVect[0] );
463  moab()->remove_child_meshset( h, iter.childVect[1] );
464  stack.push_back( iter.childVect[0] );
465  stack.push_back( iter.childVect[1] );
466  }
467 
468  rval = moab()->delete_entities( &h, 1 );
469  if( MB_SUCCESS != rval ) return rval;
470  }
471 
472  return MB_SUCCESS;
473 }
474 
475 ErrorCode AdaptiveKDTreeIter::initialize( AdaptiveKDTree* ttool,
476  EntityHandle root,
477  const double bmin[3],
478  const double bmax[3],
479  Direction direction )
480 {
481  mStack.clear();
482  treeTool = ttool;
483  mBox[BMIN][0] = bmin[0];
484  mBox[BMIN][1] = bmin[1];
485  mBox[BMIN][2] = bmin[2];
486  mBox[BMAX][0] = bmax[0];
487  mBox[BMAX][1] = bmax[1];
488  mBox[BMAX][2] = bmax[2];
489  mStack.push_back( StackObj( root, 0 ) );
490  return step_to_first_leaf( direction );
491 }
492 
493 ErrorCode AdaptiveKDTreeIter::step_to_first_leaf( Direction direction )
494 {
495  ErrorCode rval;
496  AdaptiveKDTree::Plane plane;
497  const Direction opposite = static_cast< Direction >( 1 - direction );
498 
499  for( ;; )
500  {
501  childVect.clear();
502  treeTool->treeStats.nodesVisited++; // not sure whether this is the visit or the push_back below
503  rval = treeTool->moab()->get_child_meshsets( mStack.back().entity, childVect );
504  if( MB_SUCCESS != rval ) return rval;
505  if( childVect.empty() )
506  { // leaf
507  treeTool->treeStats.leavesVisited++;
508  break;
509  }
510 
511  rval = treeTool->get_split_plane( mStack.back().entity, plane );
512  if( MB_SUCCESS != rval ) return rval;
513 
514  mStack.push_back( StackObj( childVect[direction], mBox[opposite][plane.norm] ) );
515  mBox[opposite][plane.norm] = plane.coord;
516  }
517  return MB_SUCCESS;
518 }
519 
520 ErrorCode AdaptiveKDTreeIter::step( Direction direction )
521 {
522  StackObj node, parent;
523  ErrorCode rval;
524  AdaptiveKDTree::Plane plane;
525  const Direction opposite = static_cast< Direction >( 1 - direction );
526 
527  // If stack is empty, then either this iterator is uninitialized
528  // or we reached the end of the iteration (and return
529  // MB_ENTITY_NOT_FOUND) already.
530  if( mStack.empty() ) return MB_FAILURE;
531 
532  // Pop the current node from the stack.
533  // The stack should then contain the parent of the current node.
534  // If the stack is empty after this pop, then we've reached the end.
535  node = mStack.back();
536  mStack.pop_back();
537  treeTool->treeStats.nodesVisited++;
538  if( mStack.empty() ) treeTool->treeStats.leavesVisited++;
539 
540  while( !mStack.empty() )
541  {
542  // Get data for parent entity
543  parent = mStack.back();
544  childVect.clear();
545  rval = treeTool->moab()->get_child_meshsets( parent.entity, childVect );
546  if( MB_SUCCESS != rval ) return rval;
547  rval = treeTool->get_split_plane( parent.entity, plane );
548  if( MB_SUCCESS != rval ) return rval;
549 
550  // If we're at the left child
551  if( childVect[opposite] == node.entity )
552  {
553  // change from box of left child to box of parent
554  mBox[direction][plane.norm] = node.coord;
555  // push right child on stack
556  node.entity = childVect[direction];
557  treeTool->treeStats.nodesVisited++; // changed node
558  node.coord = mBox[opposite][plane.norm];
559  mStack.push_back( node );
560  // change from box of parent to box of right child
561  mBox[opposite][plane.norm] = plane.coord;
562  // descend to left-most leaf of the right child
563  return step_to_first_leaf( opposite );
564  }
565 
566  // The current node is the right child of the parent,
567  // continue up the tree.
568  assert( childVect[direction] == node.entity );
569  mBox[opposite][plane.norm] = node.coord;
570  node = parent;
571  treeTool->treeStats.nodesVisited++;
572  mStack.pop_back();
573  }
574 
575  return MB_ENTITY_NOT_FOUND;
576 }
577 
578 ErrorCode AdaptiveKDTreeIter::get_neighbors( AdaptiveKDTree::Axis norm,
579  bool neg,
580  std::vector< AdaptiveKDTreeIter >& results,
581  double epsilon ) const
582 {
583  StackObj node, parent;
584  ErrorCode rval;
585  AdaptiveKDTree::Plane plane;
586  int child_idx;
587 
588  // Find tree node at which the specified side of the box
589  // for this node was created.
590  AdaptiveKDTreeIter iter( *this ); // temporary iterator (don't modify *this)
591  node = iter.mStack.back();
592  iter.mStack.pop_back();
593  for( ;; )
594  {
595  // reached the root - original node was on boundary (no neighbors)
596  if( iter.mStack.empty() ) return MB_SUCCESS;
597 
598  // get parent node data
599  parent = iter.mStack.back();
600  iter.childVect.clear();
601  rval = treeTool->moab()->get_child_meshsets( parent.entity, iter.childVect );
602  if( MB_SUCCESS != rval ) return rval;
603  rval = treeTool->get_split_plane( parent.entity, plane );
604  if( MB_SUCCESS != rval ) return rval;
605 
606  child_idx = iter.childVect[0] == node.entity ? 0 : 1;
607  assert( iter.childVect[child_idx] == node.entity );
608 
609  // if we found the split plane for the desired side
610  // push neighbor on stack and stop
611  if( plane.norm == norm && (int)neg == child_idx )
612  {
613  // change from box of previous child to box of parent
614  iter.mBox[1 - child_idx][plane.norm] = node.coord;
615  // push other child of parent onto stack
616  node.entity = iter.childVect[1 - child_idx];
617  node.coord = iter.mBox[child_idx][plane.norm];
618  iter.mStack.push_back( node );
619  // change from parent box to box of new child
620  iter.mBox[child_idx][plane.norm] = plane.coord;
621  break;
622  }
623 
624  // continue up the tree
625  iter.mBox[1 - child_idx][plane.norm] = node.coord;
626  node = parent;
627  iter.mStack.pop_back();
628  }
629 
630  // now move down tree, searching for adjacent boxes
631  std::vector< AdaptiveKDTreeIter > list;
632  // loop over all potential paths to neighbors (until list is empty)
633  for( ;; )
634  {
635  // follow a single path to a leaf, append any other potential
636  // paths to neighbors to 'list'
637  node = iter.mStack.back();
638  for( ;; )
639  {
640  iter.childVect.clear();
641  rval = treeTool->moab()->get_child_meshsets( node.entity, iter.childVect );
642  if( MB_SUCCESS != rval ) return rval;
643 
644  // if leaf
645  if( iter.childVect.empty() )
646  {
647  results.push_back( iter );
648  break;
649  }
650 
651  rval = treeTool->get_split_plane( node.entity, plane );
652  if( MB_SUCCESS != rval ) return rval;
653 
654  // if split parallel to side
655  if( plane.norm == norm )
656  {
657  // continue with whichever child is on the correct side of the split
658  node.entity = iter.childVect[neg];
659  node.coord = iter.mBox[1 - neg][plane.norm];
660  iter.mStack.push_back( node );
661  iter.mBox[1 - neg][plane.norm] = plane.coord;
662  }
663  // if left child is adjacent
664  else if( this->mBox[BMIN][plane.norm] - plane.coord <= epsilon )
665  {
666  // if right child is also adjacent, add to list
667  if( plane.coord - this->mBox[BMAX][plane.norm] <= epsilon )
668  {
669  list.push_back( iter );
670  list.back().mStack.push_back( StackObj( iter.childVect[1], iter.mBox[BMIN][plane.norm] ) );
671  list.back().mBox[BMIN][plane.norm] = plane.coord;
672  }
673  // continue with left child
674  node.entity = iter.childVect[0];
675  node.coord = iter.mBox[BMAX][plane.norm];
676  iter.mStack.push_back( node );
677  iter.mBox[BMAX][plane.norm] = plane.coord;
678  }
679  // right child is adjacent
680  else
681  {
682  // if left child is not adjacent, right must be or something
683  // is really messed up.
684  assert( plane.coord - this->mBox[BMAX][plane.norm] <= epsilon );
685  // continue with left child
686  node.entity = iter.childVect[1];
687  node.coord = iter.mBox[BMIN][plane.norm];
688  iter.mStack.push_back( node );
689  iter.mBox[BMIN][plane.norm] = plane.coord;
690  }
691  }
692 
693  if( list.empty() ) break;
694 
695  iter = list.back();
696  list.pop_back();
697  }
698 
699  return MB_SUCCESS;
700 }
701 
702 ErrorCode AdaptiveKDTreeIter::sibling_side( AdaptiveKDTree::Axis& axis_out, bool& neg_out ) const
703 {
704  if( mStack.size() < 2 ) // at tree root
705  return MB_ENTITY_NOT_FOUND;
706 
707  EntityHandle parent = mStack[mStack.size() - 2].entity;
708  AdaptiveKDTree::Plane plane;
709  ErrorCode rval = tool()->get_split_plane( parent, plane );
710  if( MB_SUCCESS != rval ) return MB_FAILURE;
711 
712  childVect.clear();
713  rval = tool()->moab()->get_child_meshsets( parent, childVect );
714  if( MB_SUCCESS != rval || childVect.size() != 2 ) return MB_FAILURE;
715 
716  axis_out = static_cast< AdaptiveKDTree::Axis >( plane.norm );
717  neg_out = ( childVect[1] == handle() );
718  assert( childVect[neg_out] == handle() );
719  return MB_SUCCESS;
720 }
721 
722 ErrorCode AdaptiveKDTreeIter::get_parent_split_plane( AdaptiveKDTree::Plane& plane ) const
723 {
724  if( mStack.size() < 2 ) // at tree root
725  return MB_ENTITY_NOT_FOUND;
726 
727  EntityHandle parent = mStack[mStack.size() - 2].entity;
728  return tool()->get_split_plane( parent, plane );
729 }
730 
731 bool AdaptiveKDTreeIter::is_sibling( const AdaptiveKDTreeIter& other_leaf ) const
732 {
733  const size_t s = mStack.size();
734  return ( s > 1 ) && ( s == other_leaf.mStack.size() ) &&
735  ( other_leaf.mStack[s - 2].entity == mStack[s - 2].entity ) && other_leaf.handle() != handle();
736 }
737 
738 bool AdaptiveKDTreeIter::is_sibling( EntityHandle other_leaf ) const
739 {
740  if( mStack.size() < 2 || other_leaf == handle() ) return false;
741  EntityHandle parent = mStack[mStack.size() - 2].entity;
742  childVect.clear();
743  ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
744  if( MB_SUCCESS != rval || childVect.size() != 2 )
745  {
746  assert( false );
747  return false;
748  }
749  return childVect[0] == other_leaf || childVect[1] == other_leaf;
750 }
751 
752 bool AdaptiveKDTreeIter::sibling_is_forward() const
753 {
754  if( mStack.size() < 2 ) // if root
755  return false;
756  EntityHandle parent = mStack[mStack.size() - 2].entity;
757  childVect.clear();
758  ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
759  if( MB_SUCCESS != rval || childVect.size() != 2 )
760  {
761  assert( false );
762  return false;
763  }
764  return childVect[0] == handle();
765 }
766 
767 bool AdaptiveKDTreeIter::intersect_ray( const double ray_point[3],
768  const double ray_vect[3],
769  double& t_enter,
770  double& t_exit ) const
771 {
772  treeTool->treeStats.traversalLeafObjectTests++;
773  return GeomUtil::ray_box_intersect( CartVect( box_min() ), CartVect( box_max() ), CartVect( ray_point ),
774  CartVect( ray_vect ), t_enter, t_exit );
775 }
776 
777 ErrorCode AdaptiveKDTree::intersect_children_with_elems( const Range& elems,
778  AdaptiveKDTree::Plane plane,
779  double eps,
780  CartVect box_min,
781  CartVect box_max,
782  Range& left_tris,
783  Range& right_tris,
784  Range& both_tris,
785  double& metric_value )
786 {
787  left_tris.clear();
788  right_tris.clear();
789  both_tris.clear();
790  CartVect coords[16];
791 
792  // get extents of boxes for left and right sides
793  BoundBox left_box( box_min, box_max ), right_box( box_min, box_max );
794  right_box.bMin = box_min;
795  left_box.bMax = box_max;
796  right_box.bMin[plane.norm] = left_box.bMax[plane.norm] = plane.coord;
797  const CartVect left_cen = 0.5 * ( left_box.bMax + box_min );
798  const CartVect left_dim = 0.5 * ( left_box.bMax - box_min );
799  const CartVect right_cen = 0.5 * ( box_max + right_box.bMin );
800  const CartVect right_dim = 0.5 * ( box_max - right_box.bMin );
801  const CartVect dim = box_max - box_min;
802  const double max_tol = std::max( dim[0], std::max( dim[1], dim[2] ) ) / 10;
803 
804  // test each entity
805  ErrorCode rval;
806  int count, count2;
807  const EntityHandle *conn, *conn2;
808 
809  const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
810  const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
811  const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
812  Range::iterator left_ins = left_tris.begin();
813  Range::iterator right_ins = right_tris.begin();
814  Range::iterator both_ins = both_tris.begin();
815  Range::const_iterator i;
816 
817  // vertices
818  for( i = elems.begin(); i != elem_begin; ++i )
819  {
820  tree_stats().constructLeafObjectTests++;
821  rval = moab()->get_coords( &*i, 1, coords[0].array() );
822  if( MB_SUCCESS != rval ) return rval;
823 
824  bool lo = false, ro = false;
825  if( coords[0][plane.norm] <= plane.coord ) lo = true;
826  if( coords[0][plane.norm] >= plane.coord ) ro = true;
827 
828  if( lo && ro )
829  both_ins = both_tris.insert( both_ins, *i, *i );
830  else if( lo )
831  left_ins = left_tris.insert( left_ins, *i, *i );
832  else // if (ro)
833  right_ins = right_tris.insert( right_ins, *i, *i );
834  }
835 
836  // non-polyhedron elements
837  std::vector< EntityHandle > dum_vector;
838  for( i = elem_begin; i != poly_begin; ++i )
839  {
840  tree_stats().constructLeafObjectTests++;
841  rval = moab()->get_connectivity( *i, conn, count, true, &dum_vector );
842  if( MB_SUCCESS != rval ) return rval;
843  if( count > (int)( sizeof( coords ) / sizeof( coords[0] ) ) ) return MB_FAILURE;
844  rval = moab()->get_coords( &conn[0], count, coords[0].array() );
845  if( MB_SUCCESS != rval ) return rval;
846 
847  bool lo = false, ro = false;
848  for( int j = 0; j < count; ++j )
849  {
850  if( coords[j][plane.norm] <= plane.coord ) lo = true;
851  if( coords[j][plane.norm] >= plane.coord ) ro = true;
852  }
853 
854  // Triangle must be in at least one leaf. If test against plane
855  // identified that leaf, then we're done. If triangle is on both
856  // sides of plane, do more precise test to ensure that it is really
857  // in both.
858  // BoundBox box;
859  // box.update(*moab(), *i);
860  if( lo && ro )
861  {
862  double tol = eps;
863  lo = ro = false;
864  while( !lo && !ro && tol <= max_tol )
865  {
866  tree_stats().boxElemTests += 2;
867  lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE( *i ), left_cen, left_dim + CartVect( tol ),
868  count );
869  ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE( *i ), right_cen, right_dim + CartVect( tol ),
870  count );
871 
872  tol *= 10.0;
873  }
874  }
875  if( lo && ro )
876  both_ins = both_tris.insert( both_ins, *i, *i );
877  else if( lo )
878  left_ins = left_tris.insert( left_ins, *i, *i );
879  else if( ro )
880  right_ins = right_tris.insert( right_ins, *i, *i );
881  }
882 
883  // polyhedra
884  for( i = poly_begin; i != set_begin; ++i )
885  {
886  tree_stats().constructLeafObjectTests++;
887  rval = moab()->get_connectivity( *i, conn, count, true );
888  if( MB_SUCCESS != rval ) return rval;
889 
890  // just check the bounding box of the polyhedron
891  bool lo = false, ro = false;
892  for( int j = 0; j < count; ++j )
893  {
894  rval = moab()->get_connectivity( conn[j], conn2, count2, true );
895  if( MB_SUCCESS != rval ) return rval;
896 
897  for( int k = 0; k < count2; ++k )
898  {
899  rval = moab()->get_coords( conn2 + k, 1, coords[0].array() );
900  if( MB_SUCCESS != rval ) return rval;
901  if( coords[0][plane.norm] <= plane.coord ) lo = true;
902  if( coords[0][plane.norm] >= plane.coord ) ro = true;
903  }
904  }
905 
906  if( lo && ro )
907  both_ins = both_tris.insert( both_ins, *i, *i );
908  else if( lo )
909  left_ins = left_tris.insert( left_ins, *i, *i );
910  else if( ro )
911  right_ins = right_tris.insert( right_ins, *i, *i );
912  }
913 
914  // sets
915  BoundBox tbox;
916  for( i = set_begin; i != elems.end(); ++i )
917  {
918  tree_stats().constructLeafObjectTests++;
919  rval = tbox.update( *moab(), *i, spherical, radius );
920  if( MB_SUCCESS != rval ) return rval;
921 
922  bool lo = false, ro = false;
923  if( tbox.bMin[plane.norm] <= plane.coord ) lo = true;
924  if( tbox.bMax[plane.norm] >= plane.coord ) ro = true;
925 
926  if( lo && ro )
927  both_ins = both_tris.insert( both_ins, *i, *i );
928  else if( lo )
929  left_ins = left_tris.insert( left_ins, *i, *i );
930  else // if (ro)
931  right_ins = right_tris.insert( right_ins, *i, *i );
932  }
933 
934  CartVect box_dim = box_max - box_min;
935  double area_left = left_dim[0] * left_dim[1] + left_dim[1] * left_dim[2] + left_dim[2] * left_dim[0];
936  double area_right = right_dim[0] * right_dim[1] + right_dim[1] * right_dim[2] + right_dim[2] * right_dim[0];
937  double area_both = box_dim[0] * box_dim[1] + box_dim[1] * box_dim[2] + box_dim[2] * box_dim[0];
938  metric_value = ( area_left * left_tris.size() + area_right * right_tris.size() ) / area_both + both_tris.size();
939  return MB_SUCCESS;
940 }
941 
942 ErrorCode AdaptiveKDTree::best_subdivision_plane( int num_planes,
943  const AdaptiveKDTreeIter& iter,
944  Range& best_left,
945  Range& best_right,
946  Range& best_both,
947  AdaptiveKDTree::Plane& best_plane,
948  double eps )
949 {
950  double metric_val = std::numeric_limits< unsigned >::max();
951 
952  ErrorCode r;
953  const CartVect box_min( iter.box_min() );
954  const CartVect box_max( iter.box_max() );
955  const CartVect diff( box_max - box_min );
956 
957  Range entities;
958  r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
959  if( MB_SUCCESS != r ) return r;
960  const size_t p_count = entities.size();
961 
962  for( int axis = 0; axis < 3; ++axis )
963  {
964  int plane_count = num_planes;
965  if( ( num_planes + 1 ) * eps >= diff[axis] ) plane_count = (int)( diff[axis] / eps ) - 1;
966 
967  for( int p = 1; p <= plane_count; ++p )
968  {
969  AdaptiveKDTree::Plane plane = { box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis], axis };
970  Range left, right, both;
971  double val;
972  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
973  if( MB_SUCCESS != r ) return r;
974  const size_t sdiff = p_count - both.size();
975  if( left.size() == sdiff || right.size() == sdiff ) continue;
976 
977  if( val >= metric_val ) continue;
978 
979  metric_val = val;
980  best_plane = plane;
981  best_left.swap( left );
982  best_right.swap( right );
983  best_both.swap( both );
984  }
985  }
986 
987  return MB_SUCCESS;
988 }
989 
990 ErrorCode AdaptiveKDTree::best_subdivision_snap_plane( int num_planes,
991  const AdaptiveKDTreeIter& iter,
992  Range& best_left,
993  Range& best_right,
994  Range& best_both,
995  AdaptiveKDTree::Plane& best_plane,
996  std::vector< double >& tmp_data,
997  double eps )
998 {
999  double metric_val = std::numeric_limits< unsigned >::max();
1000 
1001  ErrorCode r;
1002  // const CartVect tol(eps*diff);
1003 
1004  Range entities, vertices;
1005  r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
1006  if( MB_SUCCESS != r ) return r;
1007  const size_t p_count = entities.size();
1008  r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
1009  if( MB_SUCCESS != r ) return r;
1010 
1011  unsigned int nverts = vertices.size();
1012  tmp_data.resize( 3 * nverts );
1013  r = iter.tool()->moab()->get_coords( vertices, &tmp_data[0], &tmp_data[nverts], &tmp_data[2 * nverts] );
1014  if( MB_SUCCESS != r ) return r;
1015 
1016  // calculate bounding box of vertices
1017  // decide based on the actual box the splitting plane
1018  // do not decide based on iterator box.
1019  // it could be too big
1020  // BoundBox box;
1021  // r = box.update(*moab(), vertices);
1022  CartVect box_min;
1023  CartVect box_max;
1024  for( int dir = 0; dir < 3; dir++ )
1025  {
1026  double amin = tmp_data[dir * nverts];
1027  double amax = amin;
1028  double* p = &tmp_data[dir * nverts + 1];
1029  for( unsigned int i = 1; i < nverts; i++ )
1030  {
1031  if( *p < amin ) amin = *p;
1032  if( *p > amax ) amax = *p;
1033  p++;
1034  }
1035  box_min[dir] = amin;
1036  box_max[dir] = amax;
1037  }
1038  CartVect diff( box_max - box_min );
1039 
1040  for( int axis = 0; axis < 3; ++axis )
1041  {
1042  int plane_count = num_planes;
1043 
1044  // if num_planes results in width < eps, reset the plane count
1045  if( ( num_planes + 1 ) * eps >= diff[axis] ) plane_count = (int)( diff[axis] / eps ) - 1;
1046 
1047  for( int p = 1; p <= plane_count; ++p )
1048  {
1049 
1050  // coord of this plane on axis
1051  double coord = box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis];
1052 
1053  // find closest vertex coordinate to this plane position
1054  unsigned int istrt = axis * nverts;
1055  double closest_coord = tmp_data[istrt];
1056  for( unsigned i = 1; i < nverts; ++i )
1057  if( fabs( coord - tmp_data[istrt + i] ) < fabs( coord - closest_coord ) )
1058  closest_coord = tmp_data[istrt + i];
1059  if( closest_coord - box_min[axis] <= eps || box_max[axis] - closest_coord <= eps ) continue;
1060 
1061  // seprate elems into left/right/both, and compute separating metric
1062  AdaptiveKDTree::Plane plane = { closest_coord, axis };
1063  Range left, right, both;
1064  double val;
1065  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
1066  if( MB_SUCCESS != r ) return r;
1067  const size_t d = p_count - both.size();
1068  if( left.size() == d || right.size() == d ) continue;
1069 
1070  if( val >= metric_val ) continue;
1071 
1072  metric_val = val;
1073  best_plane = plane;
1074  best_left.swap( left );
1075  best_right.swap( right );
1076  best_both.swap( both );
1077  }
1078  }
1079 
1080  return MB_SUCCESS;
1081 }
1082 
1083 ErrorCode AdaptiveKDTree::best_vertex_median_plane( int num_planes,
1084  const AdaptiveKDTreeIter& iter,
1085  Range& best_left,
1086  Range& best_right,
1087  Range& best_both,
1088  AdaptiveKDTree::Plane& best_plane,
1089  std::vector< double >& coords,
1090  double eps )
1091 {
1092  double metric_val = std::numeric_limits< unsigned >::max();
1093 
1094  ErrorCode r;
1095  const CartVect box_min( iter.box_min() );
1096  const CartVect box_max( iter.box_max() );
1097 
1098  Range entities, vertices;
1099  r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
1100  if( MB_SUCCESS != r ) return r;
1101  const size_t p_count = entities.size();
1102  r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
1103  if( MB_SUCCESS != r ) return r;
1104 
1105  coords.resize( vertices.size() );
1106  for( int axis = 0; axis < 3; ++axis )
1107  {
1108  if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
1109 
1110  double* ptrs[] = { 0, 0, 0 };
1111  ptrs[axis] = &coords[0];
1112  r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
1113  if( MB_SUCCESS != r ) return r;
1114 
1115  std::sort( coords.begin(), coords.end() );
1116  std::vector< double >::iterator citer;
1117  citer = std::upper_bound( coords.begin(), coords.end(), box_min[axis] + eps );
1118  const size_t count = std::upper_bound( citer, coords.end(), box_max[axis] - eps ) - citer;
1119  size_t step;
1120  int np = num_planes;
1121  if( count < 2 * (size_t)num_planes )
1122  {
1123  step = 1;
1124  np = count - 1;
1125  }
1126  else
1127  {
1128  step = count / ( num_planes + 1 );
1129  }
1130 
1131  for( int p = 1; p <= np; ++p )
1132  {
1133 
1134  citer += step;
1135  AdaptiveKDTree::Plane plane = { *citer, axis };
1136  Range left, right, both;
1137  double val;
1138  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
1139  if( MB_SUCCESS != r ) return r;
1140  const size_t diff = p_count - both.size();
1141  if( left.size() == diff || right.size() == diff ) continue;
1142 
1143  if( val >= metric_val ) continue;
1144 
1145  metric_val = val;
1146  best_plane = plane;
1147  best_left.swap( left );
1148  best_right.swap( right );
1149  best_both.swap( both );
1150  }
1151  }
1152 
1153  return MB_SUCCESS;
1154 }
1155 
1156 ErrorCode AdaptiveKDTree::best_vertex_sample_plane( int num_planes,
1157  const AdaptiveKDTreeIter& iter,
1158  Range& best_left,
1159  Range& best_right,
1160  Range& best_both,
1161  AdaptiveKDTree::Plane& best_plane,
1162  std::vector< double >& coords,
1163  std::vector< EntityHandle >& indices,
1164  double eps )
1165 {
1166  const size_t random_elem_threshold = 20 * num_planes;
1167  double metric_val = std::numeric_limits< unsigned >::max();
1168 
1169  ErrorCode r;
1170  const CartVect box_min( iter.box_min() );
1171  const CartVect box_max( iter.box_max() );
1172 
1173  Range entities, vertices;
1174  r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
1175  if( MB_SUCCESS != r ) return r;
1176 
1177  // We are selecting random vertex coordinates to use for candidate split
1178  // planes. So if element list is large, begin by selecting random elements.
1179  const size_t p_count = entities.size();
1180  coords.resize( 3 * num_planes );
1181  if( p_count < random_elem_threshold )
1182  {
1183  r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
1184  if( MB_SUCCESS != r ) return r;
1185  }
1186  else
1187  {
1188  indices.resize( random_elem_threshold );
1189  const int num_rand = p_count / RAND_MAX + 1;
1190  for( size_t j = 0; j < random_elem_threshold; ++j )
1191  {
1192  size_t rnd = rand();
1193  for( int i = num_rand; i > 1; --i )
1194  rnd *= rand();
1195  rnd %= p_count;
1196  indices[j] = entities[rnd];
1197  }
1198  r = iter.tool()->moab()->get_adjacencies( &indices[0], random_elem_threshold, 0, false, vertices,
1199  Interface::UNION );
1200  if( MB_SUCCESS != r ) return r;
1201  }
1202 
1203  coords.resize( vertices.size() );
1204  for( int axis = 0; axis < 3; ++axis )
1205  {
1206  if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
1207 
1208  double* ptrs[] = { 0, 0, 0 };
1209  ptrs[axis] = &coords[0];
1210  r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
1211  if( MB_SUCCESS != r ) return r;
1212 
1213  size_t num_valid_coords = 0;
1214  for( size_t i = 0; i < coords.size(); ++i )
1215  if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) ++num_valid_coords;
1216 
1217  if( 2 * (size_t)num_planes > num_valid_coords )
1218  {
1219  indices.clear();
1220  for( size_t i = 0; i < coords.size(); ++i )
1221  if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) indices.push_back( i );
1222  }
1223  else
1224  {
1225  indices.resize( num_planes );
1226  // make sure random indices are sufficient to cover entire range
1227  const int num_rand = coords.size() / RAND_MAX + 1;
1228  for( int j = 0; j < num_planes; ++j )
1229  {
1230  size_t rnd;
1231  do
1232  {
1233  rnd = rand();
1234  for( int i = num_rand; i > 1; --i )
1235  rnd *= rand();
1236  rnd %= coords.size();
1237  } while( coords[rnd] <= box_min[axis] + eps || coords[rnd] >= box_max[axis] - eps );
1238  indices[j] = rnd;
1239  }
1240  }
1241 
1242  for( unsigned p = 0; p < indices.size(); ++p )
1243  {
1244 
1245  AdaptiveKDTree::Plane plane = { coords[indices[p]], axis };
1246  Range left, right, both;
1247  double val;
1248  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
1249  if( MB_SUCCESS != r ) return r;
1250  const size_t diff = p_count - both.size();
1251  if( left.size() == diff || right.size() == diff ) continue;
1252 
1253  if( val >= metric_val ) continue;
1254 
1255  metric_val = val;
1256  best_plane = plane;
1257  best_left.swap( left );
1258  best_right.swap( right );
1259  best_both.swap( both );
1260  }
1261  }
1262 
1263  return MB_SUCCESS;
1264 }
1265 
1266 ErrorCode AdaptiveKDTree::point_search( const double* point,
1267  EntityHandle& leaf_out,
1268  const double iter_tol,
1269  const double inside_tol,
1270  bool* multiple_leaves,
1271  EntityHandle* start_node,
1272  CartVect* params )
1273 {
1274  std::vector< EntityHandle > children;
1275  Plane plane;
1276 
1277  treeStats.numTraversals++;
1278  leaf_out = 0;
1279  BoundBox box;
1280  // kdtrees never have multiple leaves containing a pt
1281  if( multiple_leaves ) *multiple_leaves = false;
1282 
1283  EntityHandle node = ( start_node ? *start_node : myRoot );
1284 
1285  treeStats.nodesVisited++;
1286  ErrorCode rval = get_bounding_box( box, &node );
1287  if( MB_SUCCESS != rval ) return rval;
1288  if( !box.contains_point( point, iter_tol ) ) return MB_SUCCESS;
1289 
1290  rval = moab()->get_child_meshsets( node, children );
1291  if( MB_SUCCESS != rval ) return rval;
1292 
1293  while( !children.empty() )
1294  {
1295  treeStats.nodesVisited++;
1296 
1297  rval = get_split_plane( node, plane );
1298  if( MB_SUCCESS != rval ) return rval;
1299 
1300  const double d = point[plane.norm] - plane.coord;
1301  node = children[( d > 0.0 )];
1302 
1303  children.clear();
1304  rval = moab()->get_child_meshsets( node, children );
1305  if( MB_SUCCESS != rval ) return rval;
1306  }
1307 
1308  treeStats.leavesVisited++;
1309  if( myEval && params )
1310  {
1311  rval = myEval->find_containing_entity( node, point, iter_tol, inside_tol, leaf_out, params->array(),
1312  &treeStats.traversalLeafObjectTests );
1313  if( MB_SUCCESS != rval ) return rval;
1314  }
1315  else
1316  leaf_out = node;
1317 
1318  return MB_SUCCESS;
1319 }
1320 
1321 ErrorCode AdaptiveKDTree::point_search( const double* point,
1322  AdaptiveKDTreeIter& leaf_it,
1323  const double iter_tol,
1324  const double /*inside_tol*/,
1325  bool* multiple_leaves,
1326  EntityHandle* start_node )
1327 {
1328  ErrorCode rval;
1329  treeStats.numTraversals++;
1330 
1331  // kdtrees never have multiple leaves containing a pt
1332  if( multiple_leaves ) *multiple_leaves = false;
1333 
1334  leaf_it.mBox[0] = boundBox.bMin;
1335  leaf_it.mBox[1] = boundBox.bMax;
1336 
1337  // test that point is inside tree
1338  if( !boundBox.contains_point( point, iter_tol ) )
1339  {
1340  treeStats.nodesVisited++;
1341  return MB_ENTITY_NOT_FOUND;
1342  }
1343 
1344  // initialize iterator at tree root
1345  leaf_it.treeTool = this;
1346  leaf_it.mStack.clear();
1347  leaf_it.mStack.push_back( AdaptiveKDTreeIter::StackObj( ( start_node ? *start_node : myRoot ), 0 ) );
1348 
1349  // loop until we reach a leaf
1350  AdaptiveKDTree::Plane plane;
1351  for( ;; )
1352  {
1353  treeStats.nodesVisited++;
1354 
1355  // get children
1356  leaf_it.childVect.clear();
1357  rval = moab()->get_child_meshsets( leaf_it.handle(), leaf_it.childVect );
1358  if( MB_SUCCESS != rval ) return rval;
1359 
1360  // if no children, then at leaf (done)
1361  if( leaf_it.childVect.empty() )
1362  {
1363  treeStats.leavesVisited++;
1364  break;
1365  }
1366 
1367  // get split plane
1368  rval = get_split_plane( leaf_it.handle(), plane );
1369  if( MB_SUCCESS != rval ) return rval;
1370 
1371  // step iterator to appropriate child
1372  // idx: 0->left, 1->right
1373  const int idx = ( point[plane.norm] > plane.coord );
1374  leaf_it.mStack.push_back(
1375  AdaptiveKDTreeIter::StackObj( leaf_it.childVect[idx], leaf_it.mBox[1 - idx][plane.norm] ) );
1376  leaf_it.mBox[1 - idx][plane.norm] = plane.coord;
1377  }
1378 
1379  return MB_SUCCESS;
1380 }
1381 
1382 struct NodeDistance
1383 {
1384  EntityHandle handle;
1385  CartVect dist; // from_point - closest_point_on_box
1386 };
1387 
1388 ErrorCode AdaptiveKDTree::distance_search( const double from_point[3],
1389  const double distance,
1390  std::vector< EntityHandle >& result_list,
1391  const double iter_tol,
1392  const double inside_tol,
1393  std::vector< double >* result_dists,
1394  std::vector< CartVect >* result_params,
1395  EntityHandle* tree_root )
1396 {
1397  treeStats.numTraversals++;
1398  const double dist_sqr = distance * distance;
1399  const CartVect from( from_point );
1400  std::vector< NodeDistance > list,
1401  result_list_nodes; // list of subtrees to traverse, and results
1402  // pre-allocate space for default max tree depth
1403  list.reserve( maxDepth );
1404 
1405  // misc temporary values
1406  Plane plane;
1407  NodeDistance node;
1408  ErrorCode rval;
1409  std::vector< EntityHandle > children;
1410 
1411  // Get distance from input position to bounding box of tree
1412  // (zero if inside box)
1413  BoundBox box;
1414  rval = get_bounding_box( box );
1415  if( MB_SUCCESS == rval && !box.contains_point( from_point, iter_tol ) )
1416  {
1417  treeStats.nodesVisited++;
1418  return MB_SUCCESS;
1419  }
1420 
1421  // if bounding box is not available (e.g. not starting from true root)
1422  // just start with zero. Less efficient, but will work.
1423  node.dist = CartVect( 0.0 );
1424  if( MB_SUCCESS == rval )
1425  {
1426  for( int i = 0; i < 3; ++i )
1427  {
1428  if( from_point[i] < box.bMin[i] )
1429  node.dist[i] = box.bMin[i] - from_point[i];
1430  else if( from_point[i] > box.bMax[i] )
1431  node.dist[i] = from_point[i] - box.bMax[i];
1432  }
1433  if( node.dist % node.dist > dist_sqr )
1434  {
1435  treeStats.nodesVisited++;
1436  return MB_SUCCESS;
1437  }
1438  }
1439 
1440  // begin with root in list
1441  node.handle = ( tree_root ? *tree_root : myRoot );
1442  list.push_back( node );
1443 
1444  while( !list.empty() )
1445  {
1446 
1447  node = list.back();
1448  list.pop_back();
1449  treeStats.nodesVisited++;
1450 
1451  // If leaf node, test contained triangles
1452  children.clear();
1453  rval = moab()->get_child_meshsets( node.handle, children );
1454  if( children.empty() )
1455  {
1456  treeStats.leavesVisited++;
1457  if( myEval && result_params )
1458  {
1459  EntityHandle ent;
1460  CartVect params;
1461  rval = myEval->find_containing_entity( node.handle, from_point, iter_tol, inside_tol, ent,
1462  params.array(), &treeStats.traversalLeafObjectTests );
1463  if( MB_SUCCESS != rval )
1464  return rval;
1465  else if( ent )
1466  {
1467  result_list.push_back( ent );
1468  result_params->push_back( params );
1469  if( result_dists ) result_dists->push_back( 0.0 );
1470  }
1471  }
1472  else
1473  {
1474  result_list_nodes.push_back( node );
1475  continue;
1476  }
1477  }
1478 
1479  // If not leaf node, add children to working list
1480  rval = get_split_plane( node.handle, plane );
1481  if( MB_SUCCESS != rval ) return rval;
1482 
1483  const double d = from[plane.norm] - plane.coord;
1484 
1485  // right of plane?
1486  if( d > 0 )
1487  {
1488  node.handle = children[1];
1489  list.push_back( node );
1490  // if the split plane is close to the input point, add
1491  // the left child also (we'll check the exact distance
1492  /// when we pop it from the list.)
1493  if( d <= distance )
1494  {
1495  node.dist[plane.norm] = d;
1496  if( node.dist % node.dist <= dist_sqr )
1497  {
1498  node.handle = children[0];
1499  list.push_back( node );
1500  }
1501  }
1502  }
1503  // left of plane
1504  else
1505  {
1506  node.handle = children[0];
1507  list.push_back( node );
1508  // if the split plane is close to the input point, add
1509  // the right child also (we'll check the exact distance
1510  /// when we pop it from the list.)
1511  if( -d <= distance )
1512  {
1513  node.dist[plane.norm] = -d;
1514  if( node.dist % node.dist <= dist_sqr )
1515  {
1516  node.handle = children[1];
1517  list.push_back( node );
1518  }
1519  }
1520  }
1521  }
1522 
1523  if( myEval && result_params ) return MB_SUCCESS;
1524 
1525  // separate loops to avoid if test inside loop
1526 
1527  result_list.reserve( result_list_nodes.size() );
1528  for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end(); ++vit )
1529  result_list.push_back( ( *vit ).handle );
1530 
1531  if( result_dists && distance > 0.0 )
1532  {
1533  result_dists->reserve( result_list_nodes.size() );
1534  for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end();
1535  ++vit )
1536  result_dists->push_back( ( *vit ).dist.length() );
1537  }
1538 
1539  return MB_SUCCESS;
1540 }
1541 
1542 static ErrorCode closest_to_triangles( Interface* moab,
1543  const Range& tris,
1544  const CartVect& from,
1545  double& shortest_dist_sqr,
1546  CartVect& closest_pt,
1547  EntityHandle& closest_tri )
1548 {
1549  ErrorCode rval;
1550  CartVect pos, diff, verts[3];
1551  const EntityHandle* conn = NULL;
1552  int len = 0;
1553 
1554  for( Range::iterator i = tris.begin(); i != tris.end(); ++i )
1555  {
1556  rval = moab->get_connectivity( *i, conn, len );
1557  if( MB_SUCCESS != rval ) return rval;
1558 
1559  rval = moab->get_coords( conn, 3, verts[0].array() );
1560  if( MB_SUCCESS != rval ) return rval;
1561 
1562  GeomUtil::closest_location_on_tri( from, verts, pos );
1563  diff = pos - from;
1564  double dist_sqr = diff % diff;
1565  if( dist_sqr < shortest_dist_sqr )
1566  {
1567  // new closest location
1568  shortest_dist_sqr = dist_sqr;
1569  closest_pt = pos;
1570  closest_tri = *i;
1571  }
1572  }
1573 
1574  return MB_SUCCESS;
1575 }
1576 
1577 static ErrorCode closest_to_triangles( Interface* moab,
1578  EntityHandle set_handle,
1579  const CartVect& from,
1580  double& shortest_dist_sqr,
1581  CartVect& closest_pt,
1582  EntityHandle& closest_tri )
1583 {
1584  ErrorCode rval;
1585  Range tris;
1586 
1587  rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
1588  if( MB_SUCCESS != rval ) return rval;
1589 
1590  return closest_to_triangles( moab, tris, from, shortest_dist_sqr, closest_pt, closest_tri );
1591 }
1592 
1593 ErrorCode AdaptiveKDTree::find_close_triangle( EntityHandle root,
1594  const double from[3],
1595  double pt[3],
1596  EntityHandle& triangle )
1597 {
1598  ErrorCode rval;
1599  Range tris;
1600  Plane split;
1601  std::vector< EntityHandle > stack;
1602  std::vector< EntityHandle > children( 2 );
1603  stack.reserve( 30 );
1604  assert( root );
1605  stack.push_back( root );
1606 
1607  while( !stack.empty() )
1608  {
1609  EntityHandle node = stack.back();
1610  stack.pop_back();
1611 
1612  for( ;; )
1613  { // loop until we find a leaf
1614 
1615  children.clear();
1616  rval = moab()->get_child_meshsets( node, children );
1617  if( MB_SUCCESS != rval ) return rval;
1618 
1619  // loop termination criterion
1620  if( children.empty() ) break;
1621 
1622  // if not a leaf, get split plane
1623  rval = get_split_plane( node, split );
1624  if( MB_SUCCESS != rval ) return rval;
1625 
1626  // continue down the side that contains the point,
1627  // and push the other side onto the stack in case
1628  // we need to check it later.
1629  int rs = split.right_side( from );
1630  node = children[rs];
1631  stack.push_back( children[1 - rs] );
1632  }
1633 
1634  // We should now be at a leaf.
1635  // If it has some triangles, we're done.
1636  // If not, continue searching for another leaf.
1637  tris.clear();
1638  rval = moab()->get_entities_by_type( node, MBTRI, tris );
1639  if( !tris.empty() )
1640  {
1641  double dist_sqr = HUGE_VAL;
1642  CartVect point( pt );
1643  rval = closest_to_triangles( moab(), tris, CartVect( from ), dist_sqr, point, triangle );
1644  point.get( pt );
1645  return rval;
1646  }
1647  }
1648 
1649  // If we got here, then we traversed the entire tree
1650  // and all the leaves were empty.
1651  return MB_ENTITY_NOT_FOUND;
1652 }
1653 
1654 /** Find the triangles in a set that are closer to the input
1655  * position than any triangles in the 'closest_tris' list.
1656  *
1657  * closest_tris is assumed to contain a list of triangles for
1658  * which the first is the closest known triangle to the input
1659  * position and the first entry in 'closest_pts' is the closest
1660  * location on that triangle. Any other values in the lists must
1661  * be other triangles for which the closest point is within the
1662  * input tolerance of the closest closest point. This function
1663  * will update the lists as appropriate if any closer triangles
1664  * or triangles within the tolerance of the current closest location
1665  * are found. The first entry is maintained as the closest of the
1666  * list of triangles.
1667  */
1668 /*
1669  static ErrorCode closest_to_triangles( Interface* moab,
1670  EntityHandle set_handle,
1671  double tolerance,
1672  const CartVect& from,
1673  std::vector<EntityHandle>& closest_tris,
1674  std::vector<CartVect>& closest_pts )
1675  {
1676  ErrorCode rval;
1677  Range tris;
1678  CartVect pos, diff, verts[3];
1679  const EntityHandle* conn;
1680  int len;
1681  double shortest_dist_sqr = HUGE_VAL;
1682  if (!closest_pts.empty()) {
1683  diff = from - closest_pts.front();
1684  shortest_dist_sqr = diff % diff;
1685  }
1686 
1687  rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
1688  if (MB_SUCCESS != rval)
1689  return rval;
1690 
1691  for (Range::iterator i = tris.begin(); i != tris.end(); ++i) {
1692  rval = moab->get_connectivity( *i, conn, len );
1693  if (MB_SUCCESS != rval)
1694  return rval;
1695 
1696  rval = moab->get_coords( conn, 3, verts[0].array() );
1697  if (MB_SUCCESS != rval)
1698  return rval;
1699 
1700  GeomUtil::closest_location_on_tri( from, verts, pos );
1701  diff = pos - from;
1702  double dist_sqr = diff % diff;
1703  if (dist_sqr < shortest_dist_sqr) {
1704  // new closest location
1705  shortest_dist_sqr = dist_sqr;
1706 
1707  if (closest_pts.empty()) {
1708  closest_tris.push_back( *i );
1709  closest_pts.push_back( pos );
1710  }
1711  // if have a previous closest location
1712  else {
1713  // if previous closest is more than 2*tolerance away
1714  // from new closest, then nothing in the list can
1715  // be within tolerance of new closest point.
1716  diff = pos - closest_pts.front();
1717  dist_sqr = diff % diff;
1718  if (dist_sqr > 4.0 * tolerance * tolerance) {
1719  closest_tris.clear();
1720  closest_pts.clear();
1721  closest_tris.push_back( *i );
1722  closest_pts.push_back( pos );
1723  }
1724  // otherwise need to remove any triangles that are
1725  // not within tolerance of the new closest point.
1726  else {
1727  unsigned r = 0, w = 0;
1728  for (r = 0; r < closest_pts.size(); ++r) {
1729  diff = pos - closest_pts[r];
1730  if (diff % diff <= tolerance*tolerance) {
1731  closest_pts[w] = closest_pts[r];
1732  closest_tris[w] = closest_tris[r];
1733  ++w;
1734  }
1735  }
1736  closest_pts.resize( w + 1 );
1737  closest_tris.resize( w + 1 );
1738  // always put the closest one in the front
1739  if (w > 0) {
1740  closest_pts.back() = closest_pts.front();
1741  closest_tris.back() = closest_tris.front();
1742  }
1743  closest_pts.front() = pos;
1744  closest_tris.front() = *i;
1745  }
1746  }
1747  }
1748  else {
1749  // If within tolerance of old closest triangle,
1750  // add this one to the list.
1751  diff = closest_pts.front() - pos;
1752  if (diff % diff <= tolerance*tolerance) {
1753  closest_pts.push_back( pos );
1754  closest_tris.push_back( *i );
1755  }
1756  }
1757  }
1758 
1759  return MB_SUCCESS;
1760  }
1761 */
1762 
1763 ErrorCode AdaptiveKDTree::closest_triangle( EntityHandle tree_root,
1764  const double from_coords[3],
1765  double closest_point_out[3],
1766  EntityHandle& triangle_out )
1767 {
1768  ErrorCode rval;
1769  double shortest_dist_sqr = HUGE_VAL;
1770  std::vector< EntityHandle > leaves;
1771  const CartVect from( from_coords );
1772  CartVect closest_pt;
1773 
1774  // Find the leaf containing the input point
1775  // This search does not take into account any bounding box for the
1776  // tree, so it always returns one leaf.
1777  assert( tree_root );
1778  rval = find_close_triangle( tree_root, from_coords, closest_pt.array(), triangle_out );
1779  if( MB_SUCCESS != rval ) return rval;
1780 
1781  // Find any other leaves for which the bounding box is within
1782  // the same distance from the input point as the current closest
1783  // point is.
1784  CartVect diff = closest_pt - from;
1785  rval = distance_search( from_coords, sqrt( diff % diff ), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
1786  if( MB_SUCCESS != rval ) return rval;
1787 
1788  // Check any close leaves to see if they contain triangles that
1789  // are as close to or closer than the current closest triangle(s).
1790  for( unsigned i = 0; i < leaves.size(); ++i )
1791  {
1792  rval = closest_to_triangles( moab(), leaves[i], from, shortest_dist_sqr, closest_pt, triangle_out );
1793  if( MB_SUCCESS != rval ) return rval;
1794  }
1795 
1796  // pass back resulting position
1797  closest_pt.get( closest_point_out );
1798  return MB_SUCCESS;
1799 }
1800 
1801 ErrorCode AdaptiveKDTree::sphere_intersect_triangles( EntityHandle tree_root,
1802  const double center[3],
1803  double rad,
1804  std::vector< EntityHandle >& triangles )
1805 {
1806  ErrorCode rval;
1807  std::vector< EntityHandle > leaves;
1808  const CartVect from( center );
1809  CartVect closest_pt;
1810  const EntityHandle* conn;
1811  CartVect coords[3];
1812  int conn_len;
1813 
1814  // get leaves of tree that intersect sphere
1815  assert( tree_root );
1816  rval = distance_search( center, rad, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
1817  if( MB_SUCCESS != rval ) return rval;
1818 
1819  // search each leaf for triangles intersecting sphere
1820  for( unsigned i = 0; i < leaves.size(); ++i )
1821  {
1822  Range tris;
1823  rval = moab()->get_entities_by_type( leaves[i], MBTRI, tris );
1824  if( MB_SUCCESS != rval ) return rval;
1825 
1826  for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
1827  {
1828  rval = moab()->get_connectivity( *j, conn, conn_len );
1829  if( MB_SUCCESS != rval ) return rval;
1830  rval = moab()->get_coords( conn, 3, coords[0].array() );
1831  if( MB_SUCCESS != rval ) return rval;
1832  GeomUtil::closest_location_on_tri( from, coords, closest_pt );
1833  closest_pt -= from;
1834  if( ( closest_pt % closest_pt ) <= ( rad * rad ) ) triangles.push_back( *j );
1835  }
1836  }
1837 
1838  // remove duplicates from triangle list
1839  std::sort( triangles.begin(), triangles.end() );
1840  triangles.erase( std::unique( triangles.begin(), triangles.end() ), triangles.end() );
1841  return MB_SUCCESS;
1842 }
1843 
1844 struct NodeSeg
1845 {
1846  NodeSeg( EntityHandle h, double b, double e ) : handle( h ), beg( b ), end( e ) {}
1847  EntityHandle handle;
1848  double beg, end;
1849 };
1850 
1851 ErrorCode AdaptiveKDTree::ray_intersect_triangles( EntityHandle root,
1852  const double tol,
1853  const double ray_dir_in[3],
1854  const double ray_pt_in[3],
1855  std::vector< EntityHandle >& tris_out,
1856  std::vector< double >& dists_out,
1857  int max_ints,
1858  double ray_end )
1859 {
1860  ErrorCode rval;
1861  double ray_beg = 0.0;
1862  if( ray_end < 0.0 ) ray_end = HUGE_VAL;
1863 
1864  // if root has bounding box, trim ray to that box
1865  CartVect tvec( tol );
1866  BoundBox box;
1867  const CartVect ray_pt( ray_pt_in ), ray_dir( ray_dir_in );
1868  rval = get_bounding_box( box );
1869  if( MB_SUCCESS == rval )
1870  {
1871  if( !GeomUtil::segment_box_intersect( box.bMin - tvec, box.bMax + tvec, ray_pt, ray_dir, ray_beg, ray_end ) )
1872  return MB_SUCCESS; // ray misses entire tree.
1873  }
1874 
1875  Range tris;
1876  Range::iterator iter;
1877  CartVect tri_coords[3];
1878  const EntityHandle* tri_conn;
1879  int conn_len;
1880  double tri_t;
1881 
1882  Plane plane;
1883  std::vector< EntityHandle > children;
1884  std::vector< NodeSeg > list;
1885  NodeSeg seg( root, ray_beg, ray_end );
1886  list.push_back( seg );
1887 
1888  while( !list.empty() )
1889  {
1890  seg = list.back();
1891  list.pop_back();
1892 
1893  // If we are limited to a certain number of intersections
1894  // (max_ints != 0), then ray_end will contain the distance
1895  // to the furthest intersection we have so far. If the
1896  // tree node is further than that, skip it.
1897  if( seg.beg > ray_end ) continue;
1898 
1899  // Check if at a leaf
1900  children.clear();
1901  rval = moab()->get_child_meshsets( seg.handle, children );
1902  if( MB_SUCCESS != rval ) return rval;
1903  if( children.empty() )
1904  { // leaf
1905 
1906  tris.clear();
1907  rval = moab()->get_entities_by_type( seg.handle, MBTRI, tris );
1908  if( MB_SUCCESS != rval ) return rval;
1909 
1910  for( iter = tris.begin(); iter != tris.end(); ++iter )
1911  {
1912  rval = moab()->get_connectivity( *iter, tri_conn, conn_len );
1913  if( MB_SUCCESS != rval ) return rval;
1914  rval = moab()->get_coords( tri_conn, 3, tri_coords[0].array() );
1915  if( MB_SUCCESS != rval ) return rval;
1916 
1917  if( GeomUtil::ray_tri_intersect( tri_coords, ray_pt, ray_dir, tri_t, &ray_end ) )
1918  {
1919  if( !max_ints )
1920  {
1921  if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
1922  {
1923  tris_out.push_back( *iter );
1924  dists_out.push_back( tri_t );
1925  }
1926  }
1927  else if( tri_t < ray_end )
1928  {
1929  if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
1930  {
1931  if( tris_out.size() < (unsigned)max_ints )
1932  {
1933  tris_out.resize( tris_out.size() + 1 );
1934  dists_out.resize( dists_out.size() + 1 );
1935  }
1936  int w = tris_out.size() - 1;
1937  for( ; w > 0 && tri_t < dists_out[w - 1]; --w )
1938  {
1939  tris_out[w] = tris_out[w - 1];
1940  dists_out[w] = dists_out[w - 1];
1941  }
1942  tris_out[w] = *iter;
1943  dists_out[w] = tri_t;
1944  if( tris_out.size() >= (unsigned)max_ints )
1945  // when we have already reached the max intx points, we cans safely
1946  // reset ray_end, because we will accept new points only "closer"
1947  // than the last one
1948  ray_end = dists_out.back();
1949  }
1950  }
1951  }
1952  }
1953 
1954  continue;
1955  }
1956 
1957  rval = get_split_plane( seg.handle, plane );
1958  if( MB_SUCCESS != rval ) return rval;
1959 
1960  // Consider two planes that are the split plane +/- the tolerance.
1961  // Calculate the segment parameter at which the line segment intersects
1962  // the true plane, and also the difference between that value and the
1963  // intersection with either of the +/- tol planes.
1964  const double inv_dir = 1.0 / ray_dir[plane.norm]; // only do division once
1965  const double t = ( plane.coord - ray_pt[plane.norm] ) * inv_dir; // intersection with plane
1966  const double diff = tol * inv_dir; // t adjustment for +tol plane
1967  // const double t0 = t - diff; // intersection with -tol plane
1968  // const double t1 = t + diff; // intersection with +tol plane
1969 
1970  // The index of the child tree node (0 or 1) that is on the
1971  // side of the plane to which the ray direction points. That is,
1972  // if the ray direction is opposite the plane normal, the index
1973  // of the child corresponding to the side beneath the plane. If
1974  // the ray direction is the same as the plane normal, the index
1975  // of the child corresponding to the side above the plane.
1976  const int fwd_child = ( ray_dir[plane.norm] > 0.0 );
1977 
1978  // Note: we maintain seg.beg <= seg.end at all times, so assume that here.
1979 
1980  // If segment is parallel to plane
1981  if( !Util::is_finite( t ) )
1982  {
1983  if( ray_pt[plane.norm] - tol <= plane.coord ) list.push_back( NodeSeg( children[0], seg.beg, seg.end ) );
1984  if( ray_pt[plane.norm] + tol >= plane.coord ) list.push_back( NodeSeg( children[1], seg.beg, seg.end ) );
1985  }
1986  // If segment is entirely to one side of plane such that the
1987  // intersection with the split plane is past the end of the segment
1988  else if( seg.end + diff < t )
1989  {
1990  // If segment direction is opposite that of plane normal, then
1991  // being past the end of the segment means that we are to the
1992  // right (or above) the plane and what the right child (index == 1).
1993  // Otherwise we want the left child (index == 0);
1994  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
1995  }
1996  // If the segment is entirely to one side of the plane such that
1997  // the intersection with the split plane is before the start of the
1998  // segment
1999  else if( seg.beg - diff > t )
2000  {
2001  // If segment direction is opposite that of plane normal, then
2002  // being before the start of the segment means that we are to the
2003  // left (or below) the plane and what the left child (index == 0).
2004  // Otherwise we want the right child (index == 1);
2005  list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
2006  }
2007  // Otherwise we must intersect the plane.
2008  // Note: be careful not to grow the segment if t is slightly
2009  // outside the current segment, as doing so would effectively
2010  // increase the tolerance as we descend the tree.
2011  else if( t <= seg.beg )
2012  {
2013  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.beg ) );
2014  list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
2015  }
2016  else if( t >= seg.end )
2017  {
2018  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
2019  list.push_back( NodeSeg( children[fwd_child], seg.end, seg.end ) );
2020  }
2021  else
2022  {
2023  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, t ) );
2024  list.push_back( NodeSeg( children[fwd_child], t, seg.end ) );
2025  }
2026  }
2027 
2028  return MB_SUCCESS;
2029 }
2030 
2031 ErrorCode AdaptiveKDTree::compute_depth( EntityHandle root, unsigned int& min_depth, unsigned int& max_depth )
2032 {
2033  AdaptiveKDTreeIter iter;
2034  get_tree_iterator( root, iter );
2035  iter.step_to_first_leaf( AdaptiveKDTreeIter::LEFT );
2036  min_depth = max_depth = iter.depth();
2037 
2038  int num_of_elements = 0, max, min;
2039  moab()->get_number_entities_by_handle( iter.handle(), num_of_elements );
2040  max = min = num_of_elements;
2041  while( MB_SUCCESS == iter.step() )
2042  {
2043  int temp = 0;
2044  moab()->get_number_entities_by_handle( iter.handle(), temp );
2045  max = std::max( max, temp );
2046  min = std::min( min, temp );
2047  if( iter.depth() > max_depth )
2048  max_depth = iter.depth();
2049  else if( iter.depth() < min_depth )
2050  min_depth = iter.depth();
2051  }
2052  return MB_SUCCESS;
2053 }
2054 
2055 ErrorCode AdaptiveKDTree::get_info( EntityHandle root, double bmin[3], double bmax[3], unsigned int& dep )
2056 {
2057  BoundBox box;
2058  ErrorCode result = get_bounding_box( box, &root );
2059  if( MB_SUCCESS != result ) return result;
2060  box.bMin.get( bmin );
2061  box.bMax.get( bmax );
2062 
2063  unsigned min_depth;
2064  return compute_depth( root, min_depth, dep );
2065 }
2066 
2067 static std::string mem_to_string( unsigned long mem )
2068 {
2069  char unit[3] = "B";
2070  if( mem > 9 * 1024 )
2071  {
2072  mem = ( mem + 512 ) / 1024;
2073  strcpy( unit, "kB" );
2074  }
2075  if( mem > 9 * 1024 )
2076  {
2077  mem = ( mem + 512 ) / 1024;
2078  strcpy( unit, "MB" );
2079  }
2080  if( mem > 9 * 1024 )
2081  {
2082  mem = ( mem + 512 ) / 1024;
2083  strcpy( unit, "GB" );
2084  }
2085  char buffer[256];
2086  snprintf( buffer, 256, "%lu %s", mem, unit );
2087  return buffer;
2088 }
2089 
2090 template < typename T >
2091 struct SimpleStat
2092 {
2093  T min, max, sum, sqr;
2094  size_t count;
2095  SimpleStat();
2096  void add( T value );
2097  double avg() const
2098  {
2099  return (double)sum / count;
2100  }
2101  double rms() const
2102  {
2103  return sqrt( (double)sqr / count );
2104  }
2105  double dev() const
2106  {
2107  return ( count > 1
2108  ? sqrt( ( count * (double)sqr - (double)sum * (double)sum ) / ( (double)count * ( count - 1 ) ) )
2109  : 0.0 );
2110  }
2111 };
2112 
2113 template < typename T >
2114 SimpleStat< T >::SimpleStat()
2115  : min( std::numeric_limits< T >::max() ), max( std::numeric_limits< T >::min() ), sum( 0 ), sqr( 0 ), count( 0 )
2116 {
2117 }
2118 
2119 template < typename T >
2120 void SimpleStat< T >::add( T value )
2121 {
2122  if( value < min ) min = value;
2123  if( value > max ) max = value;
2124  sum += value;
2125  sqr += value * value;
2126  ++count;
2127 }
2128 
2129 ErrorCode AdaptiveKDTree::print()
2130 {
2131  Range range;
2132 
2133  Range tree_sets, elem2d, elem3d, verts, all;
2134  moab()->get_child_meshsets( myRoot, tree_sets, 0 );
2135 
2136  {
2137  std::vector<EntityHandle> elem2d_vec, elem3d_vec, verts_vec;
2138  for( Range::iterator rit = tree_sets.begin(); rit != tree_sets.end(); ++rit )
2139  {
2140  moab()->get_entities_by_dimension( *rit, 2, elem2d_vec );
2141  moab()->get_entities_by_dimension( *rit, 3, elem3d_vec );
2142  moab()->get_entities_by_type( *rit, MBVERTEX, verts_vec );
2143  }
2144  std::sort(elem2d_vec.begin(), elem2d_vec.end());
2145  std::copy( elem2d_vec.rbegin(), elem2d_vec.rend(), range_inserter( elem2d ) );
2146  std::sort(elem3d_vec.begin(), elem3d_vec.end());
2147  std::copy( elem3d_vec.rbegin(), elem3d_vec.rend(), range_inserter( elem3d ) );
2148  std::sort(verts_vec.begin(), verts_vec.end());
2149  std::copy( verts_vec.rbegin(), verts_vec.rend(), range_inserter( verts ) );
2150  }
2151 
2152  all.merge( verts );
2153  all.merge( elem2d );
2154  all.merge( elem3d );
2155  tree_sets.insert( myRoot );
2156  unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized,
2157  elem_used, elem_amortized;
2158  moab()->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0, 0,
2159  0, &set_tag_used, &set_tag_amortized );
2160  moab()->estimated_memory_use( all, &elem_used, &elem_amortized );
2161 
2162  int num_2d = 0, num_3d = 0;
2163  ;
2164  moab()->get_number_entities_by_dimension( 0, 2, num_2d );
2165  moab()->get_number_entities_by_dimension( 0, 3, num_3d );
2166 
2167  BoundBox box;
2168  ErrorCode rval = get_bounding_box( box, &myRoot );
2169  if( MB_SUCCESS != rval || box == BoundBox() ) throw rval;
2170  double diff[3] = { box.bMax[0] - box.bMin[0], box.bMax[1] - box.bMin[1], box.bMax[2] - box.bMin[2] };
2171  double tree_vol = diff[0] * diff[1] * diff[2];
2172  double tree_surf_area = 2 * ( diff[0] * diff[1] + diff[1] * diff[2] + diff[2] * diff[0] );
2173 
2174  SimpleStat< unsigned > depth, size;
2175  SimpleStat< double > vol, surf;
2176 
2177  AdaptiveKDTreeIter iter;
2178  get_tree_iterator( myRoot, iter );
2179  do
2180  {
2181  depth.add( iter.depth() );
2182 
2183  int num_leaf_elem;
2184  moab()->get_number_entities_by_handle( iter.handle(), num_leaf_elem );
2185  size.add( num_leaf_elem );
2186 
2187  const double* n = iter.box_min();
2188  const double* x = iter.box_max();
2189  double dims[3] = { x[0] - n[0], x[1] - n[1], x[2] - n[2] };
2190 
2191  double leaf_vol = dims[0] * dims[1] * dims[2];
2192  vol.add( leaf_vol );
2193 
2194  double area = 2.0 * ( dims[0] * dims[1] + dims[1] * dims[2] + dims[2] * dims[0] );
2195  surf.add( area );
2196 
2197  } while( MB_SUCCESS == iter.step() );
2198 
2199  printf( "------------------------------------------------------------------\n" );
2200  printf( "tree volume: %f\n", tree_vol );
2201  printf( "total elements: %d\n", num_2d + num_3d );
2202  printf( "number of leaves: %lu\n", (unsigned long)depth.count );
2203  printf( "number of nodes: %lu\n", (unsigned long)tree_sets.size() );
2204  printf( "volume ratio: %0.2f%%\n", 100 * ( vol.sum / tree_vol ) );
2205  printf( "surface ratio: %0.2f%%\n", 100 * ( surf.sum / tree_surf_area ) );
2206  printf( "\nmemory: used amortized\n" );
2207  printf( " ---------- ----------\n" );
2208  printf( "elements %10s %10s\n", mem_to_string( elem_used ).c_str(), mem_to_string( elem_amortized ).c_str() );
2209  printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() );
2210  printf( "sets %10s %10s\n", mem_to_string( set_store_used ).c_str(),
2211  mem_to_string( set_store_amortized ).c_str() );
2212  printf( "set tags %10s %10s\n", mem_to_string( set_tag_used ).c_str(),
2213  mem_to_string( set_tag_amortized ).c_str() );
2214  printf( "\nleaf stats: min avg rms max std.dev\n" );
2215  printf( " ---------- ---------- ---------- ---------- ----------\n" );
2216  printf( "depth %10u %10.1f %10.1f %10u %10.2f\n", depth.min, depth.avg(), depth.rms(), depth.max,
2217  depth.dev() );
2218  printf( "triangles %10u %10.1f %10.1f %10u %10.2f\n", size.min, size.avg(), size.rms(), size.max, size.dev() );
2219  printf( "volume %10.2g %10.2g %10.2g %10.2g %10.2g\n", vol.min, vol.avg(), vol.rms(), vol.max, vol.dev() );
2220  printf( "surf. area %10.2g %10.2g %10.2g %10.2g %10.2g\n", surf.min, surf.avg(), surf.rms(), surf.max,
2221  surf.dev() );
2222  printf( "------------------------------------------------------------------\n" );
2223 
2224  return MB_SUCCESS;
2225 }
2226 
2227 } // namespace moab