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