Mesh Oriented datABase  (version 5.6.0)
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  MB_CHK_ERR( moab()->get_child_meshsets( h, iter.childVect ) );
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( eps > 0.0 && std::isfinite(diff[axis]) && ( num_planes + 1 ) * eps >= diff[axis] )
1046  plane_count = (int)( diff[axis] / eps ) - 1;
1047 
1048  for( int p = 1; p <= plane_count; ++p )
1049  {
1050 
1051  // coord of this plane on axis
1052  double coord = box_min[axis] + ( p / ( 1.0 + plane_count ) ) * diff[axis];
1053 
1054  // find closest vertex coordinate to this plane position
1055  unsigned int istrt = axis * nverts;
1056  double closest_coord = tmp_data[istrt];
1057  for( unsigned i = 1; i < nverts; ++i )
1058  if( fabs( coord - tmp_data[istrt + i] ) < fabs( coord - closest_coord ) )
1059  closest_coord = tmp_data[istrt + i];
1060  if( closest_coord - box_min[axis] <= eps || box_max[axis] - closest_coord <= eps ) continue;
1061 
1062  // seprate elems into left/right/both, and compute separating metric
1063  AdaptiveKDTree::Plane plane = { closest_coord, axis };
1064  Range left, right, both;
1065  double val;
1066  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
1067  if( MB_SUCCESS != r ) return r;
1068  const size_t d = p_count - both.size();
1069  if( left.size() == d || right.size() == d ) continue;
1070 
1071  if( val >= metric_val ) continue;
1072 
1073  metric_val = val;
1074  best_plane = plane;
1075  best_left.swap( left );
1076  best_right.swap( right );
1077  best_both.swap( both );
1078  }
1079  }
1080 
1081  return MB_SUCCESS;
1082 }
1083 
1084 ErrorCode AdaptiveKDTree::best_vertex_median_plane( int num_planes,
1085  const AdaptiveKDTreeIter& iter,
1086  Range& best_left,
1087  Range& best_right,
1088  Range& best_both,
1089  AdaptiveKDTree::Plane& best_plane,
1090  std::vector< double >& coords,
1091  double eps )
1092 {
1093  double metric_val = std::numeric_limits< unsigned >::max();
1094 
1095  ErrorCode r;
1096  const CartVect box_min( iter.box_min() );
1097  const CartVect box_max( iter.box_max() );
1098 
1099  Range entities, vertices;
1100  r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
1101  if( MB_SUCCESS != r ) return r;
1102  const size_t p_count = entities.size();
1103  r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
1104  if( MB_SUCCESS != r ) return r;
1105 
1106  coords.resize( vertices.size() );
1107  for( int axis = 0; axis < 3; ++axis )
1108  {
1109  if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
1110 
1111  double* ptrs[] = { 0, 0, 0 };
1112  ptrs[axis] = &coords[0];
1113  r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
1114  if( MB_SUCCESS != r ) return r;
1115 
1116  std::sort( coords.begin(), coords.end() );
1117  std::vector< double >::iterator citer;
1118  citer = std::upper_bound( coords.begin(), coords.end(), box_min[axis] + eps );
1119  const size_t count = std::upper_bound( citer, coords.end(), box_max[axis] - eps ) - citer;
1120  size_t step;
1121  int np = num_planes;
1122  if( count < 2 * (size_t)num_planes )
1123  {
1124  step = 1;
1125  np = count - 1;
1126  }
1127  else
1128  {
1129  step = count / ( num_planes + 1 );
1130  }
1131 
1132  for( int p = 1; p <= np; ++p )
1133  {
1134 
1135  citer += step;
1136  AdaptiveKDTree::Plane plane = { *citer, axis };
1137  Range left, right, both;
1138  double val;
1139  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
1140  if( MB_SUCCESS != r ) return r;
1141  const size_t diff = p_count - both.size();
1142  if( left.size() == diff || right.size() == diff ) continue;
1143 
1144  if( val >= metric_val ) continue;
1145 
1146  metric_val = val;
1147  best_plane = plane;
1148  best_left.swap( left );
1149  best_right.swap( right );
1150  best_both.swap( both );
1151  }
1152  }
1153 
1154  return MB_SUCCESS;
1155 }
1156 
1157 ErrorCode AdaptiveKDTree::best_vertex_sample_plane( int num_planes,
1158  const AdaptiveKDTreeIter& iter,
1159  Range& best_left,
1160  Range& best_right,
1161  Range& best_both,
1162  AdaptiveKDTree::Plane& best_plane,
1163  std::vector< double >& coords,
1164  std::vector< EntityHandle >& indices,
1165  double eps )
1166 {
1167  const size_t random_elem_threshold = 20 * num_planes;
1168  double metric_val = std::numeric_limits< unsigned >::max();
1169 
1170  ErrorCode r;
1171  const CartVect box_min( iter.box_min() );
1172  const CartVect box_max( iter.box_max() );
1173 
1174  Range entities, vertices;
1175  r = iter.tool()->moab()->get_entities_by_handle( iter.handle(), entities );
1176  if( MB_SUCCESS != r ) return r;
1177 
1178  // We are selecting random vertex coordinates to use for candidate split
1179  // planes. So if element list is large, begin by selecting random elements.
1180  const size_t p_count = entities.size();
1181  coords.resize( 3 * num_planes );
1182  if( p_count < random_elem_threshold )
1183  {
1184  r = iter.tool()->moab()->get_adjacencies( entities, 0, false, vertices, Interface::UNION );
1185  if( MB_SUCCESS != r ) return r;
1186  }
1187  else
1188  {
1189  indices.resize( random_elem_threshold );
1190  const int num_rand = p_count / RAND_MAX + 1;
1191  for( size_t j = 0; j < random_elem_threshold; ++j )
1192  {
1193  size_t rnd = rand();
1194  for( int i = num_rand; i > 1; --i )
1195  rnd *= rand();
1196  rnd %= p_count;
1197  indices[j] = entities[rnd];
1198  }
1199  r = iter.tool()->moab()->get_adjacencies( &indices[0], random_elem_threshold, 0, false, vertices,
1200  Interface::UNION );
1201  if( MB_SUCCESS != r ) return r;
1202  }
1203 
1204  coords.resize( vertices.size() );
1205  for( int axis = 0; axis < 3; ++axis )
1206  {
1207  if( box_max[axis] - box_min[axis] <= 2 * eps ) continue;
1208 
1209  double* ptrs[] = { 0, 0, 0 };
1210  ptrs[axis] = &coords[0];
1211  r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
1212  if( MB_SUCCESS != r ) return r;
1213 
1214  size_t num_valid_coords = 0;
1215  for( size_t i = 0; i < coords.size(); ++i )
1216  if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) ++num_valid_coords;
1217 
1218  if( 2 * (size_t)num_planes > num_valid_coords )
1219  {
1220  indices.clear();
1221  for( size_t i = 0; i < coords.size(); ++i )
1222  if( coords[i] > box_min[axis] + eps && coords[i] < box_max[axis] - eps ) indices.push_back( i );
1223  }
1224  else
1225  {
1226  indices.resize( num_planes );
1227  // make sure random indices are sufficient to cover entire range
1228  const int num_rand = coords.size() / RAND_MAX + 1;
1229  for( int j = 0; j < num_planes; ++j )
1230  {
1231  size_t rnd;
1232  do
1233  {
1234  rnd = rand();
1235  for( int i = num_rand; i > 1; --i )
1236  rnd *= rand();
1237  rnd %= coords.size();
1238  } while( coords[rnd] <= box_min[axis] + eps || coords[rnd] >= box_max[axis] - eps );
1239  indices[j] = rnd;
1240  }
1241  }
1242 
1243  for( unsigned p = 0; p < indices.size(); ++p )
1244  {
1245 
1246  AdaptiveKDTree::Plane plane = { coords[indices[p]], axis };
1247  Range left, right, both;
1248  double val;
1249  r = intersect_children_with_elems( entities, plane, eps, box_min, box_max, left, right, both, val );
1250  if( MB_SUCCESS != r ) return r;
1251  const size_t diff = p_count - both.size();
1252  if( left.size() == diff || right.size() == diff ) continue;
1253 
1254  if( val >= metric_val ) continue;
1255 
1256  metric_val = val;
1257  best_plane = plane;
1258  best_left.swap( left );
1259  best_right.swap( right );
1260  best_both.swap( both );
1261  }
1262  }
1263 
1264  return MB_SUCCESS;
1265 }
1266 
1267 ErrorCode AdaptiveKDTree::point_search( const double* point,
1268  EntityHandle& leaf_out,
1269  const double iter_tol,
1270  const double inside_tol,
1271  bool* multiple_leaves,
1272  EntityHandle* start_node,
1273  CartVect* params )
1274 {
1275  std::vector< EntityHandle > children;
1276  Plane plane;
1277 
1278  treeStats.numTraversals++;
1279  leaf_out = 0;
1280  BoundBox box;
1281  // kdtrees never have multiple leaves containing a pt
1282  if( multiple_leaves ) *multiple_leaves = false;
1283 
1284  EntityHandle node = ( start_node ? *start_node : myRoot );
1285 
1286  treeStats.nodesVisited++;
1287  ErrorCode rval = get_bounding_box( box, &node );
1288  if( MB_SUCCESS != rval ) return rval;
1289  if( !box.contains_point( point, iter_tol ) ) return MB_SUCCESS;
1290 
1291  rval = moab()->get_child_meshsets( node, children );
1292  if( MB_SUCCESS != rval ) return rval;
1293 
1294  while( !children.empty() )
1295  {
1296  treeStats.nodesVisited++;
1297 
1298  rval = get_split_plane( node, plane );
1299  if( MB_SUCCESS != rval ) return rval;
1300 
1301  const double d = point[plane.norm] - plane.coord;
1302  node = children[( d > 0.0 )];
1303 
1304  children.clear();
1305  rval = moab()->get_child_meshsets( node, children );
1306  if( MB_SUCCESS != rval ) return rval;
1307  }
1308 
1309  treeStats.leavesVisited++;
1310  if( myEval && params )
1311  {
1312  rval = myEval->find_containing_entity( node, point, iter_tol, inside_tol, leaf_out, params->array(),
1313  &treeStats.traversalLeafObjectTests );
1314  if( MB_SUCCESS != rval ) return rval;
1315  }
1316  else
1317  leaf_out = node;
1318 
1319  return MB_SUCCESS;
1320 }
1321 
1322 ErrorCode AdaptiveKDTree::point_search( const double* point,
1323  AdaptiveKDTreeIter& leaf_it,
1324  const double iter_tol,
1325  const double /*inside_tol*/,
1326  bool* multiple_leaves,
1327  EntityHandle* start_node )
1328 {
1329  ErrorCode rval;
1330  treeStats.numTraversals++;
1331 
1332  // kdtrees never have multiple leaves containing a pt
1333  if( multiple_leaves ) *multiple_leaves = false;
1334 
1335  leaf_it.mBox[0] = boundBox.bMin;
1336  leaf_it.mBox[1] = boundBox.bMax;
1337 
1338  // test that point is inside tree
1339  if( !boundBox.contains_point( point, iter_tol ) )
1340  {
1341  treeStats.nodesVisited++;
1342  return MB_ENTITY_NOT_FOUND;
1343  }
1344 
1345  // initialize iterator at tree root
1346  leaf_it.treeTool = this;
1347  leaf_it.mStack.clear();
1348  leaf_it.mStack.push_back( AdaptiveKDTreeIter::StackObj( ( start_node ? *start_node : myRoot ), 0 ) );
1349 
1350  // loop until we reach a leaf
1351  AdaptiveKDTree::Plane plane;
1352  for( ;; )
1353  {
1354  treeStats.nodesVisited++;
1355 
1356  // get children
1357  leaf_it.childVect.clear();
1358  rval = moab()->get_child_meshsets( leaf_it.handle(), leaf_it.childVect );
1359  if( MB_SUCCESS != rval ) return rval;
1360 
1361  // if no children, then at leaf (done)
1362  if( leaf_it.childVect.empty() )
1363  {
1364  treeStats.leavesVisited++;
1365  break;
1366  }
1367 
1368  // get split plane
1369  rval = get_split_plane( leaf_it.handle(), plane );
1370  if( MB_SUCCESS != rval ) return rval;
1371 
1372  // step iterator to appropriate child
1373  // idx: 0->left, 1->right
1374  const int idx = ( point[plane.norm] > plane.coord );
1375  leaf_it.mStack.push_back(
1376  AdaptiveKDTreeIter::StackObj( leaf_it.childVect[idx], leaf_it.mBox[1 - idx][plane.norm] ) );
1377  leaf_it.mBox[1 - idx][plane.norm] = plane.coord;
1378  }
1379 
1380  return MB_SUCCESS;
1381 }
1382 
1383 struct NodeDistance
1384 {
1385  EntityHandle handle;
1386  CartVect dist; // from_point - closest_point_on_box
1387 };
1388 
1389 ErrorCode AdaptiveKDTree::distance_search( const double from_point[3],
1390  const double distance,
1391  std::vector< EntityHandle >& result_list,
1392  const double iter_tol,
1393  const double inside_tol,
1394  std::vector< double >* result_dists,
1395  std::vector< CartVect >* result_params,
1396  EntityHandle* tree_root )
1397 {
1398  treeStats.numTraversals++;
1399  const double dist_sqr = distance * distance;
1400  const CartVect from( from_point );
1401  std::vector< NodeDistance > list,
1402  result_list_nodes; // list of subtrees to traverse, and results
1403  // pre-allocate space for default max tree depth
1404  list.reserve( maxDepth );
1405 
1406  // misc temporary values
1407  Plane plane;
1408  NodeDistance node;
1409  ErrorCode rval;
1410  std::vector< EntityHandle > children;
1411 
1412  // Get distance from input position to bounding box of tree
1413  // (zero if inside box)
1414  BoundBox box;
1415  rval = get_bounding_box( box );
1416  if( MB_SUCCESS == rval && !box.contains_point( from_point, iter_tol ) )
1417  {
1418  treeStats.nodesVisited++;
1419  return MB_SUCCESS;
1420  }
1421 
1422  // if bounding box is not available (e.g. not starting from true root)
1423  // just start with zero. Less efficient, but will work.
1424  node.dist = CartVect( 0.0 );
1425  if( MB_SUCCESS == rval )
1426  {
1427  for( int i = 0; i < 3; ++i )
1428  {
1429  if( from_point[i] < box.bMin[i] )
1430  node.dist[i] = box.bMin[i] - from_point[i];
1431  else if( from_point[i] > box.bMax[i] )
1432  node.dist[i] = from_point[i] - box.bMax[i];
1433  }
1434  if( node.dist % node.dist > dist_sqr )
1435  {
1436  treeStats.nodesVisited++;
1437  return MB_SUCCESS;
1438  }
1439  }
1440 
1441  // begin with root in list
1442  node.handle = ( tree_root ? *tree_root : myRoot );
1443  list.push_back( node );
1444 
1445  while( !list.empty() )
1446  {
1447 
1448  node = list.back();
1449  list.pop_back();
1450  treeStats.nodesVisited++;
1451 
1452  // If leaf node, test contained triangles
1453  children.clear();
1454  rval = moab()->get_child_meshsets( node.handle, children );
1455  if( children.empty() )
1456  {
1457  treeStats.leavesVisited++;
1458  if( myEval && result_params )
1459  {
1460  EntityHandle ent;
1461  CartVect params;
1462  rval = myEval->find_containing_entity( node.handle, from_point, iter_tol, inside_tol, ent,
1463  params.array(), &treeStats.traversalLeafObjectTests );
1464  if( MB_SUCCESS != rval )
1465  return rval;
1466  else if( ent )
1467  {
1468  result_list.push_back( ent );
1469  result_params->push_back( params );
1470  if( result_dists ) result_dists->push_back( 0.0 );
1471  }
1472  }
1473  else
1474  {
1475  result_list_nodes.push_back( node );
1476  continue;
1477  }
1478  }
1479 
1480  // If not leaf node, add children to working list
1481  rval = get_split_plane( node.handle, plane );
1482  if( MB_SUCCESS != rval ) return rval;
1483 
1484  const double d = from[plane.norm] - plane.coord;
1485 
1486  // right of plane?
1487  if( d > 0 )
1488  {
1489  node.handle = children[1];
1490  list.push_back( node );
1491  // if the split plane is close to the input point, add
1492  // the left child also (we'll check the exact distance
1493  /// when we pop it from the list.)
1494  if( d <= distance )
1495  {
1496  node.dist[plane.norm] = d;
1497  if( node.dist % node.dist <= dist_sqr )
1498  {
1499  node.handle = children[0];
1500  list.push_back( node );
1501  }
1502  }
1503  }
1504  // left of plane
1505  else
1506  {
1507  node.handle = children[0];
1508  list.push_back( node );
1509  // if the split plane is close to the input point, add
1510  // the right child also (we'll check the exact distance
1511  /// when we pop it from the list.)
1512  if( -d <= distance )
1513  {
1514  node.dist[plane.norm] = -d;
1515  if( node.dist % node.dist <= dist_sqr )
1516  {
1517  node.handle = children[1];
1518  list.push_back( node );
1519  }
1520  }
1521  }
1522  }
1523 
1524  if( myEval && result_params ) return MB_SUCCESS;
1525 
1526  // separate loops to avoid if test inside loop
1527 
1528  result_list.reserve( result_list_nodes.size() );
1529  for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end(); ++vit )
1530  result_list.push_back( ( *vit ).handle );
1531 
1532  if( result_dists && distance > 0.0 )
1533  {
1534  result_dists->reserve( result_list_nodes.size() );
1535  for( std::vector< NodeDistance >::iterator vit = result_list_nodes.begin(); vit != result_list_nodes.end();
1536  ++vit )
1537  result_dists->push_back( ( *vit ).dist.length() );
1538  }
1539 
1540  return MB_SUCCESS;
1541 }
1542 
1543 static ErrorCode closest_to_triangles( Interface* moab,
1544  const Range& tris,
1545  const CartVect& from,
1546  double& shortest_dist_sqr,
1547  CartVect& closest_pt,
1548  EntityHandle& closest_tri )
1549 {
1550  ErrorCode rval;
1551  CartVect pos, diff, verts[3];
1552  const EntityHandle* conn = NULL;
1553  int len = 0;
1554 
1555  for( Range::iterator i = tris.begin(); i != tris.end(); ++i )
1556  {
1557  rval = moab->get_connectivity( *i, conn, len );
1558  if( MB_SUCCESS != rval ) return rval;
1559 
1560  rval = moab->get_coords( conn, 3, verts[0].array() );
1561  if( MB_SUCCESS != rval ) return rval;
1562 
1563  GeomUtil::closest_location_on_tri( from, verts, pos );
1564  diff = pos - from;
1565  double dist_sqr = diff % diff;
1566  if( dist_sqr < shortest_dist_sqr )
1567  {
1568  // new closest location
1569  shortest_dist_sqr = dist_sqr;
1570  closest_pt = pos;
1571  closest_tri = *i;
1572  }
1573  }
1574 
1575  return MB_SUCCESS;
1576 }
1577 
1578 static ErrorCode closest_to_triangles( Interface* moab,
1579  EntityHandle set_handle,
1580  const CartVect& from,
1581  double& shortest_dist_sqr,
1582  CartVect& closest_pt,
1583  EntityHandle& closest_tri )
1584 {
1585  ErrorCode rval;
1586  Range tris;
1587 
1588  rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
1589  if( MB_SUCCESS != rval ) return rval;
1590 
1591  return closest_to_triangles( moab, tris, from, shortest_dist_sqr, closest_pt, closest_tri );
1592 }
1593 
1594 ErrorCode AdaptiveKDTree::find_close_triangle( EntityHandle root,
1595  const double from[3],
1596  double pt[3],
1597  EntityHandle& triangle )
1598 {
1599  ErrorCode rval;
1600  Range tris;
1601  Plane split;
1602  std::vector< EntityHandle > stack;
1603  std::vector< EntityHandle > children( 2 );
1604  stack.reserve( 30 );
1605  assert( root );
1606  stack.push_back( root );
1607 
1608  while( !stack.empty() )
1609  {
1610  EntityHandle node = stack.back();
1611  stack.pop_back();
1612 
1613  for( ;; )
1614  { // loop until we find a leaf
1615 
1616  children.clear();
1617  rval = moab()->get_child_meshsets( node, children );
1618  if( MB_SUCCESS != rval ) return rval;
1619 
1620  // loop termination criterion
1621  if( children.empty() ) break;
1622 
1623  // if not a leaf, get split plane
1624  rval = get_split_plane( node, split );
1625  if( MB_SUCCESS != rval ) return rval;
1626 
1627  // continue down the side that contains the point,
1628  // and push the other side onto the stack in case
1629  // we need to check it later.
1630  int rs = split.right_side( from );
1631  node = children[rs];
1632  stack.push_back( children[1 - rs] );
1633  }
1634 
1635  // We should now be at a leaf.
1636  // If it has some triangles, we're done.
1637  // If not, continue searching for another leaf.
1638  tris.clear();
1639  rval = moab()->get_entities_by_type( node, MBTRI, tris );
1640  if( !tris.empty() )
1641  {
1642  double dist_sqr = HUGE_VAL;
1643  CartVect point( pt );
1644  rval = closest_to_triangles( moab(), tris, CartVect( from ), dist_sqr, point, triangle );
1645  point.get( pt );
1646  return rval;
1647  }
1648  }
1649 
1650  // If we got here, then we traversed the entire tree
1651  // and all the leaves were empty.
1652  return MB_ENTITY_NOT_FOUND;
1653 }
1654 
1655 /** Find the triangles in a set that are closer to the input
1656  * position than any triangles in the 'closest_tris' list.
1657  *
1658  * closest_tris is assumed to contain a list of triangles for
1659  * which the first is the closest known triangle to the input
1660  * position and the first entry in 'closest_pts' is the closest
1661  * location on that triangle. Any other values in the lists must
1662  * be other triangles for which the closest point is within the
1663  * input tolerance of the closest closest point. This function
1664  * will update the lists as appropriate if any closer triangles
1665  * or triangles within the tolerance of the current closest location
1666  * are found. The first entry is maintained as the closest of the
1667  * list of triangles.
1668  */
1669 /*
1670  static ErrorCode closest_to_triangles( Interface* moab,
1671  EntityHandle set_handle,
1672  double tolerance,
1673  const CartVect& from,
1674  std::vector<EntityHandle>& closest_tris,
1675  std::vector<CartVect>& closest_pts )
1676  {
1677  ErrorCode rval;
1678  Range tris;
1679  CartVect pos, diff, verts[3];
1680  const EntityHandle* conn;
1681  int len;
1682  double shortest_dist_sqr = HUGE_VAL;
1683  if (!closest_pts.empty()) {
1684  diff = from - closest_pts.front();
1685  shortest_dist_sqr = diff % diff;
1686  }
1687 
1688  rval = moab->get_entities_by_type( set_handle, MBTRI, tris );
1689  if (MB_SUCCESS != rval)
1690  return rval;
1691 
1692  for (Range::iterator i = tris.begin(); i != tris.end(); ++i) {
1693  rval = moab->get_connectivity( *i, conn, len );
1694  if (MB_SUCCESS != rval)
1695  return rval;
1696 
1697  rval = moab->get_coords( conn, 3, verts[0].array() );
1698  if (MB_SUCCESS != rval)
1699  return rval;
1700 
1701  GeomUtil::closest_location_on_tri( from, verts, pos );
1702  diff = pos - from;
1703  double dist_sqr = diff % diff;
1704  if (dist_sqr < shortest_dist_sqr) {
1705  // new closest location
1706  shortest_dist_sqr = dist_sqr;
1707 
1708  if (closest_pts.empty()) {
1709  closest_tris.push_back( *i );
1710  closest_pts.push_back( pos );
1711  }
1712  // if have a previous closest location
1713  else {
1714  // if previous closest is more than 2*tolerance away
1715  // from new closest, then nothing in the list can
1716  // be within tolerance of new closest point.
1717  diff = pos - closest_pts.front();
1718  dist_sqr = diff % diff;
1719  if (dist_sqr > 4.0 * tolerance * tolerance) {
1720  closest_tris.clear();
1721  closest_pts.clear();
1722  closest_tris.push_back( *i );
1723  closest_pts.push_back( pos );
1724  }
1725  // otherwise need to remove any triangles that are
1726  // not within tolerance of the new closest point.
1727  else {
1728  unsigned r = 0, w = 0;
1729  for (r = 0; r < closest_pts.size(); ++r) {
1730  diff = pos - closest_pts[r];
1731  if (diff % diff <= tolerance*tolerance) {
1732  closest_pts[w] = closest_pts[r];
1733  closest_tris[w] = closest_tris[r];
1734  ++w;
1735  }
1736  }
1737  closest_pts.resize( w + 1 );
1738  closest_tris.resize( w + 1 );
1739  // always put the closest one in the front
1740  if (w > 0) {
1741  closest_pts.back() = closest_pts.front();
1742  closest_tris.back() = closest_tris.front();
1743  }
1744  closest_pts.front() = pos;
1745  closest_tris.front() = *i;
1746  }
1747  }
1748  }
1749  else {
1750  // If within tolerance of old closest triangle,
1751  // add this one to the list.
1752  diff = closest_pts.front() - pos;
1753  if (diff % diff <= tolerance*tolerance) {
1754  closest_pts.push_back( pos );
1755  closest_tris.push_back( *i );
1756  }
1757  }
1758  }
1759 
1760  return MB_SUCCESS;
1761  }
1762 */
1763 
1764 ErrorCode AdaptiveKDTree::closest_triangle( EntityHandle tree_root,
1765  const double from_coords[3],
1766  double closest_point_out[3],
1767  EntityHandle& triangle_out )
1768 {
1769  ErrorCode rval;
1770  double shortest_dist_sqr = HUGE_VAL;
1771  std::vector< EntityHandle > leaves;
1772  const CartVect from( from_coords );
1773  CartVect closest_pt;
1774 
1775  // Find the leaf containing the input point
1776  // This search does not take into account any bounding box for the
1777  // tree, so it always returns one leaf.
1778  assert( tree_root );
1779  rval = find_close_triangle( tree_root, from_coords, closest_pt.array(), triangle_out );
1780  if( MB_SUCCESS != rval ) return rval;
1781 
1782  // Find any other leaves for which the bounding box is within
1783  // the same distance from the input point as the current closest
1784  // point is.
1785  CartVect diff = closest_pt - from;
1786  rval = distance_search( from_coords, sqrt( diff % diff ), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
1787  if( MB_SUCCESS != rval ) return rval;
1788 
1789  // Check any close leaves to see if they contain triangles that
1790  // are as close to or closer than the current closest triangle(s).
1791  for( unsigned i = 0; i < leaves.size(); ++i )
1792  {
1793  rval = closest_to_triangles( moab(), leaves[i], from, shortest_dist_sqr, closest_pt, triangle_out );
1794  if( MB_SUCCESS != rval ) return rval;
1795  }
1796 
1797  // pass back resulting position
1798  closest_pt.get( closest_point_out );
1799  return MB_SUCCESS;
1800 }
1801 
1802 ErrorCode AdaptiveKDTree::sphere_intersect_triangles( EntityHandle tree_root,
1803  const double center[3],
1804  double rad,
1805  std::vector< EntityHandle >& triangles )
1806 {
1807  ErrorCode rval;
1808  std::vector< EntityHandle > leaves;
1809  const CartVect from( center );
1810  CartVect closest_pt;
1811  const EntityHandle* conn;
1812  CartVect coords[3];
1813  int conn_len;
1814 
1815  // get leaves of tree that intersect sphere
1816  assert( tree_root );
1817  rval = distance_search( center, rad, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root );
1818  if( MB_SUCCESS != rval ) return rval;
1819 
1820  // search each leaf for triangles intersecting sphere
1821  for( unsigned i = 0; i < leaves.size(); ++i )
1822  {
1823  Range tris;
1824  rval = moab()->get_entities_by_type( leaves[i], MBTRI, tris );
1825  if( MB_SUCCESS != rval ) return rval;
1826 
1827  for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
1828  {
1829  rval = moab()->get_connectivity( *j, conn, conn_len );
1830  if( MB_SUCCESS != rval ) return rval;
1831  rval = moab()->get_coords( conn, 3, coords[0].array() );
1832  if( MB_SUCCESS != rval ) return rval;
1833  GeomUtil::closest_location_on_tri( from, coords, closest_pt );
1834  closest_pt -= from;
1835  if( ( closest_pt % closest_pt ) <= ( rad * rad ) ) triangles.push_back( *j );
1836  }
1837  }
1838 
1839  // remove duplicates from triangle list
1840  std::sort( triangles.begin(), triangles.end() );
1841  triangles.erase( std::unique( triangles.begin(), triangles.end() ), triangles.end() );
1842  return MB_SUCCESS;
1843 }
1844 
1845 struct NodeSeg
1846 {
1847  NodeSeg( EntityHandle h, double b, double e ) : handle( h ), beg( b ), end( e ) {}
1848  EntityHandle handle;
1849  double beg, end;
1850 };
1851 
1852 ErrorCode AdaptiveKDTree::ray_intersect_triangles( EntityHandle root,
1853  const double tol,
1854  const double ray_dir_in[3],
1855  const double ray_pt_in[3],
1856  std::vector< EntityHandle >& tris_out,
1857  std::vector< double >& dists_out,
1858  int max_ints,
1859  double ray_end )
1860 {
1861  ErrorCode rval;
1862  double ray_beg = 0.0;
1863  if( ray_end < 0.0 ) ray_end = HUGE_VAL;
1864 
1865  // if root has bounding box, trim ray to that box
1866  CartVect tvec( tol );
1867  BoundBox box;
1868  const CartVect ray_pt( ray_pt_in ), ray_dir( ray_dir_in );
1869  rval = get_bounding_box( box );
1870  if( MB_SUCCESS == rval )
1871  {
1872  if( !GeomUtil::segment_box_intersect( box.bMin - tvec, box.bMax + tvec, ray_pt, ray_dir, ray_beg, ray_end ) )
1873  return MB_SUCCESS; // ray misses entire tree.
1874  }
1875 
1876  Range tris;
1877  Range::iterator iter;
1878  CartVect tri_coords[3];
1879  const EntityHandle* tri_conn;
1880  int conn_len;
1881  double tri_t;
1882 
1883  Plane plane;
1884  std::vector< EntityHandle > children;
1885  std::vector< NodeSeg > list;
1886  NodeSeg seg( root, ray_beg, ray_end );
1887  list.push_back( seg );
1888 
1889  while( !list.empty() )
1890  {
1891  seg = list.back();
1892  list.pop_back();
1893 
1894  // If we are limited to a certain number of intersections
1895  // (max_ints != 0), then ray_end will contain the distance
1896  // to the furthest intersection we have so far. If the
1897  // tree node is further than that, skip it.
1898  if( seg.beg > ray_end ) continue;
1899 
1900  // Check if at a leaf
1901  children.clear();
1902  rval = moab()->get_child_meshsets( seg.handle, children );
1903  if( MB_SUCCESS != rval ) return rval;
1904  if( children.empty() )
1905  { // leaf
1906 
1907  tris.clear();
1908  rval = moab()->get_entities_by_type( seg.handle, MBTRI, tris );
1909  if( MB_SUCCESS != rval ) return rval;
1910 
1911  for( iter = tris.begin(); iter != tris.end(); ++iter )
1912  {
1913  rval = moab()->get_connectivity( *iter, tri_conn, conn_len );
1914  if( MB_SUCCESS != rval ) return rval;
1915  rval = moab()->get_coords( tri_conn, 3, tri_coords[0].array() );
1916  if( MB_SUCCESS != rval ) return rval;
1917 
1918  if( GeomUtil::ray_tri_intersect( tri_coords, ray_pt, ray_dir, tri_t, &ray_end ) )
1919  {
1920  if( !max_ints )
1921  {
1922  if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
1923  {
1924  tris_out.push_back( *iter );
1925  dists_out.push_back( tri_t );
1926  }
1927  }
1928  else if( tri_t < ray_end )
1929  {
1930  if( std::find( tris_out.begin(), tris_out.end(), *iter ) == tris_out.end() )
1931  {
1932  if( tris_out.size() < (unsigned)max_ints )
1933  {
1934  tris_out.resize( tris_out.size() + 1 );
1935  dists_out.resize( dists_out.size() + 1 );
1936  }
1937  int w = tris_out.size() - 1;
1938  for( ; w > 0 && tri_t < dists_out[w - 1]; --w )
1939  {
1940  tris_out[w] = tris_out[w - 1];
1941  dists_out[w] = dists_out[w - 1];
1942  }
1943  tris_out[w] = *iter;
1944  dists_out[w] = tri_t;
1945  if( tris_out.size() >= (unsigned)max_ints )
1946  // when we have already reached the max intx points, we cans safely
1947  // reset ray_end, because we will accept new points only "closer"
1948  // than the last one
1949  ray_end = dists_out.back();
1950  }
1951  }
1952  }
1953  }
1954 
1955  continue;
1956  }
1957 
1958  rval = get_split_plane( seg.handle, plane );
1959  if( MB_SUCCESS != rval ) return rval;
1960 
1961  // Consider two planes that are the split plane +/- the tolerance.
1962  // Calculate the segment parameter at which the line segment intersects
1963  // the true plane, and also the difference between that value and the
1964  // intersection with either of the +/- tol planes.
1965  const double inv_dir = 1.0 / ray_dir[plane.norm]; // only do division once
1966  const double t = ( plane.coord - ray_pt[plane.norm] ) * inv_dir; // intersection with plane
1967  const double diff = tol * inv_dir; // t adjustment for +tol plane
1968  // const double t0 = t - diff; // intersection with -tol plane
1969  // const double t1 = t + diff; // intersection with +tol plane
1970 
1971  // The index of the child tree node (0 or 1) that is on the
1972  // side of the plane to which the ray direction points. That is,
1973  // if the ray direction is opposite the plane normal, the index
1974  // of the child corresponding to the side beneath the plane. If
1975  // the ray direction is the same as the plane normal, the index
1976  // of the child corresponding to the side above the plane.
1977  const int fwd_child = ( ray_dir[plane.norm] > 0.0 );
1978 
1979  // Note: we maintain seg.beg <= seg.end at all times, so assume that here.
1980 
1981  // If segment is parallel to plane
1982  if( !Util::is_finite( t ) )
1983  {
1984  if( ray_pt[plane.norm] - tol <= plane.coord ) list.push_back( NodeSeg( children[0], seg.beg, seg.end ) );
1985  if( ray_pt[plane.norm] + tol >= plane.coord ) list.push_back( NodeSeg( children[1], seg.beg, seg.end ) );
1986  }
1987  // If segment is entirely to one side of plane such that the
1988  // intersection with the split plane is past the end of the segment
1989  else if( seg.end + diff < t )
1990  {
1991  // If segment direction is opposite that of plane normal, then
1992  // being past the end of the segment means that we are to the
1993  // right (or above) the plane and what the right child (index == 1).
1994  // Otherwise we want the left child (index == 0);
1995  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
1996  }
1997  // If the segment is entirely to one side of the plane such that
1998  // the intersection with the split plane is before the start of the
1999  // segment
2000  else if( seg.beg - diff > t )
2001  {
2002  // If segment direction is opposite that of plane normal, then
2003  // being before the start of the segment means that we are to the
2004  // left (or below) the plane and what the left child (index == 0).
2005  // Otherwise we want the right child (index == 1);
2006  list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
2007  }
2008  // Otherwise we must intersect the plane.
2009  // Note: be careful not to grow the segment if t is slightly
2010  // outside the current segment, as doing so would effectively
2011  // increase the tolerance as we descend the tree.
2012  else if( t <= seg.beg )
2013  {
2014  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.beg ) );
2015  list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) );
2016  }
2017  else if( t >= seg.end )
2018  {
2019  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, seg.end ) );
2020  list.push_back( NodeSeg( children[fwd_child], seg.end, seg.end ) );
2021  }
2022  else
2023  {
2024  list.push_back( NodeSeg( children[1 - fwd_child], seg.beg, t ) );
2025  list.push_back( NodeSeg( children[fwd_child], t, seg.end ) );
2026  }
2027  }
2028 
2029  return MB_SUCCESS;
2030 }
2031 
2032 ErrorCode AdaptiveKDTree::compute_depth( EntityHandle root, unsigned int& min_depth, unsigned int& max_depth )
2033 {
2034  AdaptiveKDTreeIter iter;
2035  get_tree_iterator( root, iter );
2036  iter.step_to_first_leaf( AdaptiveKDTreeIter::LEFT );
2037  min_depth = max_depth = iter.depth();
2038 
2039  int num_of_elements = 0, max, min;
2040  moab()->get_number_entities_by_handle( iter.handle(), num_of_elements );
2041  max = min = num_of_elements;
2042  while( MB_SUCCESS == iter.step() )
2043  {
2044  int temp = 0;
2045  moab()->get_number_entities_by_handle( iter.handle(), temp );
2046  max = std::max( max, temp );
2047  min = std::min( min, temp );
2048  if( iter.depth() > max_depth )
2049  max_depth = iter.depth();
2050  else if( iter.depth() < min_depth )
2051  min_depth = iter.depth();
2052  }
2053  return MB_SUCCESS;
2054 }
2055 
2056 ErrorCode AdaptiveKDTree::get_info( EntityHandle root, double bmin[3], double bmax[3], unsigned int& dep )
2057 {
2058  BoundBox box;
2059  ErrorCode result = get_bounding_box( box, &root );
2060  if( MB_SUCCESS != result ) return result;
2061  box.bMin.get( bmin );
2062  box.bMax.get( bmax );
2063 
2064  unsigned min_depth;
2065  return compute_depth( root, min_depth, dep );
2066 }
2067 
2068 static std::string mem_to_string( unsigned long mem )
2069 {
2070  char unit[3] = "B";
2071  if( mem > 9 * 1024 )
2072  {
2073  mem = ( mem + 512 ) / 1024;
2074  strcpy( unit, "kB" );
2075  }
2076  if( mem > 9 * 1024 )
2077  {
2078  mem = ( mem + 512 ) / 1024;
2079  strcpy( unit, "MB" );
2080  }
2081  if( mem > 9 * 1024 )
2082  {
2083  mem = ( mem + 512 ) / 1024;
2084  strcpy( unit, "GB" );
2085  }
2086  char buffer[256];
2087  snprintf( buffer, 256, "%lu %s", mem, unit );
2088  return buffer;
2089 }
2090 
2091 template < typename T >
2092 struct SimpleStat
2093 {
2094  T min, max, sum, sqr;
2095  size_t count;
2096  SimpleStat();
2097  void add( T value );
2098  double avg() const
2099  {
2100  return (double)sum / count;
2101  }
2102  double rms() const
2103  {
2104  return sqrt( (double)sqr / count );
2105  }
2106  double dev() const
2107  {
2108  return ( count > 1
2109  ? sqrt( ( count * (double)sqr - (double)sum * (double)sum ) / ( (double)count * ( count - 1 ) ) )
2110  : 0.0 );
2111  }
2112 };
2113 
2114 template < typename T >
2115 SimpleStat< T >::SimpleStat()
2116  : min( std::numeric_limits< T >::max() ), max( std::numeric_limits< T >::min() ), sum( 0 ), sqr( 0 ), count( 0 )
2117 {
2118 }
2119 
2120 template < typename T >
2121 void SimpleStat< T >::add( T value )
2122 {
2123  if( value < min ) min = value;
2124  if( value > max ) max = value;
2125  sum += value;
2126  sqr += value * value;
2127  ++count;
2128 }
2129 
2130 ErrorCode AdaptiveKDTree::print()
2131 {
2132  Range range;
2133 
2134  Range tree_sets, elem2d, elem3d, verts, all;
2135  moab()->get_child_meshsets( myRoot, tree_sets, 0 );
2136 
2137  {
2138  std::vector< EntityHandle > elem2d_vec, elem3d_vec, verts_vec;
2139  for( Range::iterator rit = tree_sets.begin(); rit != tree_sets.end(); ++rit )
2140  {
2141  moab()->get_entities_by_dimension( *rit, 2, elem2d_vec );
2142  moab()->get_entities_by_dimension( *rit, 3, elem3d_vec );
2143  moab()->get_entities_by_type( *rit, MBVERTEX, verts_vec );
2144  }
2145  std::sort( elem2d_vec.begin(), elem2d_vec.end() );
2146  std::copy( elem2d_vec.rbegin(), elem2d_vec.rend(), range_inserter( elem2d ) );
2147  std::sort( elem3d_vec.begin(), elem3d_vec.end() );
2148  std::copy( elem3d_vec.rbegin(), elem3d_vec.rend(), range_inserter( elem3d ) );
2149  std::sort( verts_vec.begin(), verts_vec.end() );
2150  std::copy( verts_vec.rbegin(), verts_vec.rend(), range_inserter( verts ) );
2151  }
2152 
2153  all.merge( verts );
2154  all.merge( elem2d );
2155  all.merge( elem3d );
2156  tree_sets.insert( myRoot );
2157  unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized,
2158  elem_used, elem_amortized;
2159  moab()->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0, 0,
2160  0, &set_tag_used, &set_tag_amortized );
2161  moab()->estimated_memory_use( all, &elem_used, &elem_amortized );
2162 
2163  int num_2d = 0, num_3d = 0;
2164  ;
2165  moab()->get_number_entities_by_dimension( 0, 2, num_2d );
2166  moab()->get_number_entities_by_dimension( 0, 3, num_3d );
2167 
2168  BoundBox box;
2169  ErrorCode rval = get_bounding_box( box, &myRoot );
2170  if( MB_SUCCESS != rval || box == BoundBox() ) throw rval;
2171  double diff[3] = { box.bMax[0] - box.bMin[0], box.bMax[1] - box.bMin[1], box.bMax[2] - box.bMin[2] };
2172  double tree_vol = diff[0] * diff[1] * diff[2];
2173  double tree_surf_area = 2 * ( diff[0] * diff[1] + diff[1] * diff[2] + diff[2] * diff[0] );
2174 
2175  SimpleStat< unsigned > depth, size;
2176  SimpleStat< double > vol, surf;
2177 
2178  AdaptiveKDTreeIter iter;
2179  get_tree_iterator( myRoot, iter );
2180  do
2181  {
2182  depth.add( iter.depth() );
2183 
2184  int num_leaf_elem;
2185  moab()->get_number_entities_by_handle( iter.handle(), num_leaf_elem );
2186  size.add( num_leaf_elem );
2187 
2188  const double* n = iter.box_min();
2189  const double* x = iter.box_max();
2190  double dims[3] = { x[0] - n[0], x[1] - n[1], x[2] - n[2] };
2191 
2192  double leaf_vol = dims[0] * dims[1] * dims[2];
2193  vol.add( leaf_vol );
2194 
2195  double area = 2.0 * ( dims[0] * dims[1] + dims[1] * dims[2] + dims[2] * dims[0] );
2196  surf.add( area );
2197 
2198  } while( MB_SUCCESS == iter.step() );
2199 
2200  printf( "------------------------------------------------------------------\n" );
2201  printf( "tree volume: %f\n", tree_vol );
2202  printf( "total elements: %d\n", num_2d + num_3d );
2203  printf( "number of leaves: %lu\n", (unsigned long)depth.count );
2204  printf( "number of nodes: %lu\n", (unsigned long)tree_sets.size() );
2205  printf( "volume ratio: %0.2f%%\n", 100 * ( vol.sum / tree_vol ) );
2206  printf( "surface ratio: %0.2f%%\n", 100 * ( surf.sum / tree_surf_area ) );
2207  printf( "\nmemory: used amortized\n" );
2208  printf( " ---------- ----------\n" );
2209  printf( "elements %10s %10s\n", mem_to_string( elem_used ).c_str(), mem_to_string( elem_amortized ).c_str() );
2210  printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() );
2211  printf( "sets %10s %10s\n", mem_to_string( set_store_used ).c_str(),
2212  mem_to_string( set_store_amortized ).c_str() );
2213  printf( "set tags %10s %10s\n", mem_to_string( set_tag_used ).c_str(),
2214  mem_to_string( set_tag_amortized ).c_str() );
2215  printf( "\nleaf stats: min avg rms max std.dev\n" );
2216  printf( " ---------- ---------- ---------- ---------- ----------\n" );
2217  printf( "depth %10u %10.1f %10.1f %10u %10.2f\n", depth.min, depth.avg(), depth.rms(), depth.max,
2218  depth.dev() );
2219  printf( "triangles %10u %10.1f %10.1f %10u %10.2f\n", size.min, size.avg(), size.rms(), size.max, size.dev() );
2220  printf( "volume %10.2g %10.2g %10.2g %10.2g %10.2g\n", vol.min, vol.avg(), vol.rms(), vol.max, vol.dev() );
2221  printf( "surf. area %10.2g %10.2g %10.2g %10.2g %10.2g\n", surf.min, surf.avg(), surf.rms(), surf.max,
2222  surf.dev() );
2223  printf( "------------------------------------------------------------------\n" );
2224 
2225  return MB_SUCCESS;
2226 }
2227 
2228 } // namespace moab