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