 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
1.9.1.