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 2008 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 BSPTree.cpp
17 : : *\author Jason Kraftcheck ([email protected])
18 : : *\date 2008-05-13
19 : : */
20 : :
21 : : #include "moab/BSPTree.hpp"
22 : : #include "moab/GeomUtil.hpp"
23 : : #include "moab/Range.hpp"
24 : : #include "Internals.hpp"
25 : : #include "moab/BSPTreePoly.hpp"
26 : : #include "moab/Util.hpp"
27 : :
28 : : #include <assert.h>
29 : : #include <string.h>
30 : : #include <algorithm>
31 : : #include <limits>
32 : :
33 : : #if defined( MOAB_HAVE_IEEEFP_H )
34 : : #include <ieeefp.h>
35 : : #endif
36 : :
37 : : #define MB_BSP_TREE_DEFAULT_TAG_NAME "BSPTree"
38 : :
39 : : namespace moab
40 : : {
41 : :
42 : 12 : static void corners_from_box( const double box_min[3], const double box_max[3], double corners[8][3] )
43 : : {
44 : 12 : const double* ranges[] = { box_min, box_max };
45 [ + + ]: 36 : for( int z = 0; z < 2; ++z )
46 : : {
47 : 24 : corners[4 * z][0] = box_min[0];
48 : 24 : corners[4 * z][1] = box_min[1];
49 : 24 : corners[4 * z][2] = ranges[z][2];
50 : :
51 : 24 : corners[4 * z + 1][0] = box_max[0];
52 : 24 : corners[4 * z + 1][1] = box_min[1];
53 : 24 : corners[4 * z + 1][2] = ranges[z][2];
54 : :
55 : 24 : corners[4 * z + 2][0] = box_max[0];
56 : 24 : corners[4 * z + 2][1] = box_max[1];
57 : 24 : corners[4 * z + 2][2] = ranges[z][2];
58 : :
59 : 24 : corners[4 * z + 3][0] = box_min[0];
60 : 24 : corners[4 * z + 3][1] = box_max[1];
61 : 24 : corners[4 * z + 3][2] = ranges[z][2];
62 : : }
63 : 12 : }
64 : :
65 : : // assume box has planar sides
66 : : // test if point is contained in box
67 : 6 : static bool point_in_box( const double corners[8][3], const double point[3] )
68 : : {
69 : 6 : const unsigned side_verts[6][3] = { { 0, 3, 1 }, { 4, 5, 7 }, { 0, 1, 4 }, { 1, 2, 5 }, { 2, 3, 6 }, { 3, 0, 7 } };
70 : : // If we assume planar sides, then the box is the intersection
71 : : // of 6 half-spaces defined by the planes of the sides.
72 [ + - ]: 6 : const CartVect pt( point );
73 [ + + ]: 36 : for( unsigned s = 0; s < 6; ++s )
74 : : {
75 [ + - ]: 32 : CartVect v0( corners[side_verts[s][0]] );
76 [ + - ]: 32 : CartVect v1( corners[side_verts[s][1]] );
77 [ + - ]: 32 : CartVect v2( corners[side_verts[s][2]] );
78 [ + - ][ + - ]: 32 : CartVect N = ( v1 - v0 ) * ( v2 - v0 );
[ + - ]
79 [ + - ][ + - ]: 32 : if( ( v0 - pt ) % N < 0 ) return false;
[ + + ]
80 : : }
81 : 6 : return true;
82 : : }
83 : :
84 : 4 : void BSPTree::Plane::set( const double pt1[3], const double pt2[3], const double pt3[3] )
85 : : {
86 : 4 : const double v1[] = { pt2[0] - pt1[0], pt2[1] - pt1[1], pt2[2] - pt1[2] };
87 : 4 : const double v2[] = { pt3[0] - pt1[0], pt3[1] - pt1[1], pt3[2] - pt1[2] };
88 : 8 : const double nrm[] = { v1[1] * v2[2] - v1[2] * v2[1], v1[2] * v2[0] - v1[0] * v2[2],
89 : 8 : v1[0] * v2[1] - v1[1] * v2[0] };
90 [ + - ]: 4 : set( nrm, pt1 );
91 : 4 : }
92 : :
93 : 15 : ErrorCode BSPTree::init_tags( const char* tagname )
94 : : {
95 [ + - ]: 15 : if( !tagname ) tagname = MB_BSP_TREE_DEFAULT_TAG_NAME;
96 : :
97 [ + - ]: 15 : std::string rootname( tagname );
98 [ + - ]: 15 : rootname += "_box";
99 : :
100 [ + - ][ + - ]: 15 : ErrorCode rval = moab()->tag_get_handle( tagname, 4, MB_TYPE_DOUBLE, planeTag, MB_TAG_CREAT | MB_TAG_DENSE );
101 [ - + ]: 15 : if( MB_SUCCESS != rval )
102 : 0 : planeTag = 0;
103 : : else
104 [ + - ][ + - ]: 15 : rval = moab()->tag_get_handle( rootname.c_str(), 24, MB_TYPE_DOUBLE, rootTag, MB_TAG_CREAT | MB_TAG_SPARSE );
105 [ - + ]: 15 : if( MB_SUCCESS != rval ) rootTag = 0;
106 : 15 : return rval;
107 : : }
108 : :
109 : 15 : BSPTree::BSPTree( Interface* mb, const char* tagname, unsigned set_flags )
110 : 15 : : mbInstance( mb ), meshSetFlags( set_flags ), cleanUpTrees( false )
111 : : {
112 [ + - ]: 15 : init_tags( tagname );
113 : 15 : }
114 : :
115 : 0 : BSPTree::BSPTree( Interface* mb, bool destroy_created_trees, const char* tagname, unsigned set_flags )
116 : 0 : : mbInstance( mb ), meshSetFlags( set_flags ), cleanUpTrees( destroy_created_trees )
117 : : {
118 [ # # ]: 0 : init_tags( tagname );
119 : 0 : }
120 : :
121 [ - + ]: 30 : BSPTree::~BSPTree()
122 : : {
123 [ + - ]: 15 : if( !cleanUpTrees ) return;
124 : :
125 [ # # ]: 0 : while( !createdTrees.empty() )
126 : : {
127 : 0 : EntityHandle tree = createdTrees.back();
128 : : // make sure this is a tree (rather than some other, stale handle)
129 : 0 : const void* data_ptr = 0;
130 : 0 : ErrorCode rval = moab()->tag_get_by_ptr( rootTag, &tree, 1, &data_ptr );
131 [ # # ]: 0 : if( MB_SUCCESS == rval ) rval = delete_tree( tree );
132 [ # # ]: 0 : if( MB_SUCCESS != rval ) createdTrees.pop_back();
133 : : }
134 : 0 : }
135 : :
136 : 101 : ErrorCode BSPTree::set_split_plane( EntityHandle node, const Plane& p )
137 : : {
138 : : // check for unit-length normal
139 : 101 : const double lensqr = p.norm[0] * p.norm[0] + p.norm[1] * p.norm[1] + p.norm[2] * p.norm[2];
140 [ + + ]: 101 : if( fabs( lensqr - 1.0 ) < std::numeric_limits< double >::epsilon() )
141 [ + - ][ + - ]: 50 : return moab()->tag_set_data( planeTag, &node, 1, &p );
142 : :
143 : 51 : const double inv_len = 1.0 / sqrt( lensqr );
144 : 51 : Plane p2( p );
145 : 51 : p2.norm[0] *= inv_len;
146 : 51 : p2.norm[1] *= inv_len;
147 : 51 : p2.norm[2] *= inv_len;
148 : 51 : p2.coeff *= inv_len;
149 : :
150 : : // check for zero-length normal
151 [ + - ][ - + ]: 51 : if( !Util::is_finite( p2.norm[0] + p2.norm[1] + p2.norm[2] + p2.coeff ) ) return MB_FAILURE;
152 : :
153 : : // store plane
154 [ + - ][ + - ]: 101 : return moab()->tag_set_data( planeTag, &node, 1, &p2 );
155 : : }
156 : :
157 : 0 : ErrorCode BSPTree::set_tree_box( EntityHandle root_handle, const double box_min[3], const double box_max[3] )
158 : : {
159 : : double corners[8][3];
160 : 0 : corners_from_box( box_min, box_max, corners );
161 [ # # ]: 0 : return set_tree_box( root_handle, corners );
162 : : }
163 : :
164 : 15 : ErrorCode BSPTree::set_tree_box( EntityHandle root_handle, const double corners[8][3] )
165 : : {
166 : 15 : return moab()->tag_set_data( rootTag, &root_handle, 1, corners );
167 : : }
168 : :
169 : 72 : ErrorCode BSPTree::get_tree_box( EntityHandle root_handle, double corners[8][3] )
170 : : {
171 : 72 : return moab()->tag_get_data( rootTag, &root_handle, 1, corners );
172 : : }
173 : :
174 : 24 : ErrorCode BSPTree::get_tree_box( EntityHandle root_handle, double corners[24] )
175 : : {
176 : 24 : return moab()->tag_get_data( rootTag, &root_handle, 1, corners );
177 : : }
178 : :
179 : 3 : ErrorCode BSPTree::create_tree( EntityHandle& root_handle )
180 : : {
181 : 3 : const double min[3] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL };
182 : 3 : const double max[3] = { HUGE_VAL, HUGE_VAL, HUGE_VAL };
183 [ + - ]: 3 : return create_tree( min, max, root_handle );
184 : : }
185 : :
186 : 15 : ErrorCode BSPTree::create_tree( const double corners[8][3], EntityHandle& root_handle )
187 : : {
188 : 15 : ErrorCode rval = moab()->create_meshset( meshSetFlags, root_handle );
189 [ - + ]: 15 : if( MB_SUCCESS != rval ) return rval;
190 : :
191 : 15 : rval = set_tree_box( root_handle, corners );
192 [ - + ]: 15 : if( MB_SUCCESS != rval )
193 : : {
194 : 0 : moab()->delete_entities( &root_handle, 1 );
195 : 0 : root_handle = 0;
196 : 0 : return rval;
197 : : }
198 : :
199 : 15 : createdTrees.push_back( root_handle );
200 : 15 : return MB_SUCCESS;
201 : : }
202 : :
203 : 12 : ErrorCode BSPTree::create_tree( const double box_min[3], const double box_max[3], EntityHandle& root_handle )
204 : : {
205 : : double corners[8][3];
206 : 12 : corners_from_box( box_min, box_max, corners );
207 [ + - ]: 12 : return create_tree( corners, root_handle );
208 : : }
209 : :
210 : 0 : ErrorCode BSPTree::delete_tree( EntityHandle root_handle )
211 : : {
212 : : ErrorCode rval;
213 : :
214 [ # # ][ # # ]: 0 : std::vector< EntityHandle > children, dead_sets, current_sets;
[ # # ]
215 [ # # ]: 0 : current_sets.push_back( root_handle );
216 [ # # ]: 0 : while( !current_sets.empty() )
217 : : {
218 [ # # ]: 0 : EntityHandle set = current_sets.back();
219 [ # # ]: 0 : current_sets.pop_back();
220 [ # # ]: 0 : dead_sets.push_back( set );
221 [ # # ][ # # ]: 0 : rval = moab()->get_child_meshsets( set, children );
222 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
223 [ # # ][ # # ]: 0 : std::copy( children.begin(), children.end(), std::back_inserter( current_sets ) );
224 : 0 : children.clear();
225 : : }
226 : :
227 [ # # ][ # # ]: 0 : rval = moab()->tag_delete_data( rootTag, &root_handle, 1 );
228 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
229 : :
230 [ # # ][ # # ]: 0 : createdTrees.erase( std::remove( createdTrees.begin(), createdTrees.end(), root_handle ), createdTrees.end() );
231 [ # # ][ # # ]: 0 : return moab()->delete_entities( &dead_sets[0], dead_sets.size() );
[ # # ]
232 : : }
233 : :
234 : 0 : ErrorCode BSPTree::find_all_trees( Range& results )
235 : : {
236 : 0 : return moab()->get_entities_by_type_and_tag( 0, MBENTITYSET, &rootTag, 0, 1, results );
237 : : }
238 : :
239 : 30 : ErrorCode BSPTree::get_tree_iterator( EntityHandle root, BSPTreeIter& iter )
240 : : {
241 : 30 : ErrorCode rval = iter.initialize( this, root );
242 [ - + ]: 30 : if( MB_SUCCESS != rval ) return rval;
243 : 30 : return iter.step_to_first_leaf( BSPTreeIter::LEFT );
244 : : }
245 : :
246 : 2 : ErrorCode BSPTree::get_tree_end_iterator( EntityHandle root, BSPTreeIter& iter )
247 : : {
248 : 2 : ErrorCode rval = iter.initialize( this, root );
249 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
250 : 2 : return iter.step_to_first_leaf( BSPTreeIter::RIGHT );
251 : : }
252 : :
253 : 101 : ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane, EntityHandle& left, EntityHandle& right )
254 : : {
255 : : ErrorCode rval;
256 : :
257 : 101 : rval = moab()->create_meshset( meshSetFlags, left );
258 [ - + ]: 101 : if( MB_SUCCESS != rval ) return rval;
259 : :
260 : 101 : rval = moab()->create_meshset( meshSetFlags, right );
261 [ - + ]: 101 : if( MB_SUCCESS != rval )
262 : : {
263 : 0 : moab()->delete_entities( &left, 1 );
264 : 0 : return rval;
265 : : }
266 : :
267 [ + - ][ + + ]: 303 : if( MB_SUCCESS != set_split_plane( leaf.handle(), plane ) ||
268 [ + - ]: 202 : MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), left ) ||
269 [ + - ][ + + ]: 303 : MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), right ) ||
270 : 101 : MB_SUCCESS != leaf.step_to_first_leaf( BSPTreeIter::LEFT ) )
271 : : {
272 : 42 : EntityHandle children[] = { left, right };
273 [ + - ][ + - ]: 42 : moab()->delete_entities( children, 2 );
274 : 42 : return MB_FAILURE;
275 : : }
276 : :
277 : 101 : return MB_SUCCESS;
278 : : }
279 : :
280 : 82 : ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane )
281 : : {
282 : : EntityHandle left, right;
283 [ + - ]: 82 : return split_leaf( leaf, plane, left, right );
284 : : }
285 : :
286 : 0 : ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane, const Range& left_entities, const Range& right_entities )
287 : : {
288 [ # # ]: 0 : EntityHandle left, right, parent = leaf.handle();
289 [ # # ]: 0 : ErrorCode rval = split_leaf( leaf, plane, left, right );
290 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
291 : :
292 [ # # ][ # # ]: 0 : if( MB_SUCCESS == moab()->add_entities( left, left_entities ) &&
[ # # ][ # # ]
293 [ # # ][ # # ]: 0 : MB_SUCCESS == moab()->add_entities( right, right_entities ) &&
[ # # ][ # # ]
294 [ # # ][ # # ]: 0 : MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
295 : 0 : return MB_SUCCESS;
296 : :
297 [ # # ][ # # ]: 0 : moab()->remove_child_meshset( parent, left );
298 [ # # ][ # # ]: 0 : moab()->remove_child_meshset( parent, right );
299 : 0 : EntityHandle children[] = { left, right };
300 [ # # ][ # # ]: 0 : moab()->delete_entities( children, 2 );
301 : 0 : return MB_FAILURE;
302 : : }
303 : :
304 : 19 : ErrorCode BSPTree::split_leaf( BSPTreeIter& leaf, Plane plane, const std::vector< EntityHandle >& left_entities,
305 : : const std::vector< EntityHandle >& right_entities )
306 : : {
307 [ + - ]: 19 : EntityHandle left, right, parent = leaf.handle();
308 [ + - ]: 19 : ErrorCode rval = split_leaf( leaf, plane, left, right );
309 [ - + ]: 19 : if( MB_SUCCESS != rval ) return rval;
310 : :
311 [ + - ][ + - ]: 57 : if( MB_SUCCESS == moab()->add_entities( left, &left_entities[0], left_entities.size() ) &&
[ + - ][ + - ]
[ + - ]
312 [ + - ][ + - ]: 38 : MB_SUCCESS == moab()->add_entities( right, &right_entities[0], right_entities.size() ) &&
[ + - ][ + - ]
[ + - ]
313 [ + - ][ + - ]: 19 : MB_SUCCESS == moab()->clear_meshset( &parent, 1 ) )
314 : 19 : return MB_SUCCESS;
315 : :
316 [ # # ][ # # ]: 0 : moab()->remove_child_meshset( parent, left );
317 [ # # ][ # # ]: 0 : moab()->remove_child_meshset( parent, right );
318 : 0 : EntityHandle children[] = { left, right };
319 [ # # ][ # # ]: 0 : moab()->delete_entities( children, 2 );
320 : 19 : return MB_FAILURE;
321 : : }
322 : :
323 : 2 : ErrorCode BSPTree::merge_leaf( BSPTreeIter& iter )
324 : : {
325 : : ErrorCode rval;
326 [ + - ][ - + ]: 2 : if( iter.depth() == 1 ) // at root
327 : 0 : return MB_FAILURE;
328 : :
329 : : // Move iter to parent
330 [ + - ]: 2 : iter.up();
331 : :
332 : : // Get sets to merge
333 [ + - ]: 2 : EntityHandle parent = iter.handle();
334 : 2 : iter.childVect.clear();
335 [ + - ][ + - ]: 2 : rval = moab()->get_child_meshsets( parent, iter.childVect );
336 [ - + ]: 2 : if( MB_SUCCESS != rval ) return rval;
337 : :
338 : : // Remove child links
339 [ + - ][ + - ]: 2 : moab()->remove_child_meshset( parent, iter.childVect[0] );
[ + - ]
340 [ + - ][ + - ]: 2 : moab()->remove_child_meshset( parent, iter.childVect[1] );
[ + - ]
341 [ + - ]: 2 : std::vector< EntityHandle > stack( iter.childVect );
342 : :
343 : : // Get all entities from children and put them in parent
344 [ + - ]: 4 : Range range;
345 [ + + ]: 6 : while( !stack.empty() )
346 : : {
347 [ + - ]: 4 : EntityHandle h = stack.back();
348 [ + - ]: 4 : stack.pop_back();
349 [ + - ]: 4 : range.clear();
350 [ + - ][ + - ]: 4 : rval = moab()->get_entities_by_handle( h, range );
351 [ - + ]: 4 : if( MB_SUCCESS != rval ) return rval;
352 [ + - ][ + - ]: 4 : rval = moab()->add_entities( parent, range );
353 [ - + ]: 4 : if( MB_SUCCESS != rval ) return rval;
354 : :
355 : 4 : iter.childVect.clear();
356 [ + - ][ + - ]: 4 : rval = moab()->get_child_meshsets( h, iter.childVect );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
357 [ - + ]: 4 : if( !iter.childVect.empty() )
358 : : {
359 [ # # ][ # # ]: 0 : moab()->remove_child_meshset( h, iter.childVect[0] );
[ # # ]
360 [ # # ][ # # ]: 0 : moab()->remove_child_meshset( h, iter.childVect[1] );
[ # # ]
361 [ # # ][ # # ]: 0 : stack.push_back( iter.childVect[0] );
362 [ # # ][ # # ]: 0 : stack.push_back( iter.childVect[1] );
363 : : }
364 : :
365 [ + - ][ + - ]: 4 : rval = moab()->delete_entities( &h, 1 );
366 [ - + ]: 4 : if( MB_SUCCESS != rval ) return rval;
367 : : }
368 : :
369 : 4 : return MB_SUCCESS;
370 : : }
371 : :
372 : 46 : ErrorCode BSPTreeIter::initialize( BSPTree* btool, EntityHandle root, const double* /*point*/ )
373 : : {
374 : 46 : treeTool = btool;
375 : 46 : mStack.clear();
376 : 46 : mStack.push_back( root );
377 : 46 : return MB_SUCCESS;
378 : : }
379 : :
380 : 119 : ErrorCode BSPTreeIter::step_to_first_leaf( Direction direction )
381 : : {
382 : : ErrorCode rval;
383 : : for( ;; )
384 : : {
385 : 191 : childVect.clear();
386 : 191 : rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
387 [ - + ]: 191 : if( MB_SUCCESS != rval ) return rval;
388 [ + + ]: 191 : if( childVect.empty() ) // leaf
389 : 119 : break;
390 : :
391 : 72 : mStack.push_back( childVect[direction] );
392 : : }
393 : 191 : return MB_SUCCESS;
394 : : }
395 : :
396 : 72 : ErrorCode BSPTreeIter::step( Direction direction )
397 : : {
398 : : EntityHandle node, parent;
399 : : ErrorCode rval;
400 : 72 : const Direction opposite = static_cast< Direction >( 1 - direction );
401 : :
402 : : // If stack is empty, then either this iterator is uninitialized
403 : : // or we reached the end of the iteration (and return
404 : : // MB_ENTITY_NOT_FOUND) already.
405 [ - + ]: 72 : if( mStack.empty() ) return MB_FAILURE;
406 : :
407 : : // Pop the current node from the stack.
408 : : // The stack should then contain the parent of the current node.
409 : : // If the stack is empty after this pop, then we've reached the end.
410 : 72 : node = mStack.back();
411 : 72 : mStack.pop_back();
412 : :
413 [ + + ]: 124 : while( !mStack.empty() )
414 : : {
415 : : // Get data for parent entity
416 : 115 : parent = mStack.back();
417 : 115 : childVect.clear();
418 : 115 : rval = tool()->moab()->get_child_meshsets( parent, childVect );
419 [ - + ]: 115 : if( MB_SUCCESS != rval ) return rval;
420 : :
421 : : // If we're at the left child
422 [ + + ]: 115 : if( childVect[opposite] == node )
423 : : {
424 : : // push right child on stack
425 : 63 : mStack.push_back( childVect[direction] );
426 : : // descend to left-most leaf of the right child
427 : 63 : return step_to_first_leaf( opposite );
428 : : }
429 : :
430 : : // The current node is the right child of the parent,
431 : : // continue up the tree.
432 [ - + ]: 52 : assert( childVect[direction] == node );
433 : 52 : node = parent;
434 : 52 : mStack.pop_back();
435 : : }
436 : :
437 : 9 : return MB_ENTITY_NOT_FOUND;
438 : : }
439 : :
440 : 0 : ErrorCode BSPTreeIter::up()
441 : : {
442 [ # # ]: 0 : if( mStack.size() < 2 ) return MB_ENTITY_NOT_FOUND;
443 : 0 : mStack.pop_back();
444 : 0 : return MB_SUCCESS;
445 : : }
446 : :
447 : 16 : ErrorCode BSPTreeIter::down( const BSPTree::Plane& /*plane*/, Direction dir )
448 : : {
449 : 16 : childVect.clear();
450 : 16 : ErrorCode rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
451 [ - + ]: 16 : if( MB_SUCCESS != rval ) return rval;
452 [ - + ]: 16 : if( childVect.empty() ) return MB_ENTITY_NOT_FOUND;
453 : :
454 : 16 : mStack.push_back( childVect[dir] );
455 : 16 : return MB_SUCCESS;
456 : : }
457 : :
458 : 35 : ErrorCode BSPTreeIter::get_parent_split_plane( BSPTree::Plane& plane ) const
459 : : {
460 [ + + ]: 35 : if( mStack.size() < 2 ) // at tree root
461 : 2 : return MB_ENTITY_NOT_FOUND;
462 : :
463 : 33 : EntityHandle parent = mStack[mStack.size() - 2];
464 : 33 : return tool()->get_split_plane( parent, plane );
465 : : }
466 : :
467 : 4 : double BSPTreeIter::volume() const
468 : : {
469 [ + - ]: 4 : BSPTreePoly polyhedron;
470 [ + - ]: 4 : ErrorCode rval = calculate_polyhedron( polyhedron );
471 [ + - ][ + - ]: 4 : return MB_SUCCESS == rval ? polyhedron.volume() : -1.0;
472 : : }
473 : :
474 : 10 : bool BSPTreeIter::is_sibling( const BSPTreeIter& other_leaf ) const
475 : : {
476 : 10 : const size_t s = mStack.size();
477 [ + - ][ + - ]: 16 : return ( s > 1 ) && ( s == other_leaf.mStack.size() ) && ( other_leaf.mStack[s - 2] == mStack[s - 2] ) &&
[ + + + + ]
478 : 16 : other_leaf.handle() != handle();
479 : : }
480 : :
481 : 9 : bool BSPTreeIter::is_sibling( EntityHandle other_leaf ) const
482 : : {
483 [ + - ][ + + ]: 9 : if( mStack.size() < 2 || other_leaf == handle() ) return false;
[ + + ]
484 : 8 : EntityHandle parent = mStack[mStack.size() - 2];
485 : 8 : childVect.clear();
486 : 8 : ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
487 [ + - ][ - + ]: 8 : if( MB_SUCCESS != rval || childVect.size() != 2 )
[ - + ]
488 : : {
489 : 0 : assert( false );
490 : : return false;
491 : : }
492 [ + + ][ + + ]: 8 : return childVect[0] == other_leaf || childVect[1] == other_leaf;
493 : : }
494 : :
495 : 5 : bool BSPTreeIter::sibling_is_forward() const
496 : : {
497 [ - + ]: 5 : if( mStack.size() < 2 ) // if root
498 : 0 : return false;
499 : 5 : EntityHandle parent = mStack[mStack.size() - 2];
500 : 5 : childVect.clear();
501 : 5 : ErrorCode rval = tool()->moab()->get_child_meshsets( parent, childVect );
502 [ + - ][ - + ]: 5 : if( MB_SUCCESS != rval || childVect.size() != 2 )
[ - + ]
503 : : {
504 : 0 : assert( false );
505 : : return false;
506 : : }
507 : 5 : return childVect[0] == handle();
508 : : }
509 : :
510 : 24 : ErrorCode BSPTreeIter::calculate_polyhedron( BSPTreePoly& poly_out ) const
511 : : {
512 : : ErrorCode rval;
513 : :
514 [ - + ]: 24 : assert( sizeof( CartVect ) == 3 * sizeof( double ) );
515 [ + - ][ + + ]: 216 : CartVect corners[8];
516 [ + - ][ + - ]: 24 : rval = treeTool->get_tree_box( mStack.front(), corners[0].array() );
[ + - ]
517 [ - + ]: 24 : if( MB_SUCCESS != rval ) return rval;
518 : :
519 [ + - ]: 24 : rval = poly_out.set( corners );
520 [ - + ]: 24 : if( MB_SUCCESS != rval ) return rval;
521 : :
522 [ + - ]: 24 : BSPTree::Plane plane;
523 : 24 : std::vector< EntityHandle >::const_iterator i = mStack.begin();
524 [ + - ]: 24 : std::vector< EntityHandle >::const_iterator here = mStack.end() - 1;
525 [ + - ][ + + ]: 123 : while( i != here )
526 : : {
527 [ + - ][ + - ]: 99 : rval = treeTool->get_split_plane( *i, plane );
528 [ - + ]: 99 : if( MB_SUCCESS != rval ) return rval;
529 : :
530 : 99 : childVect.clear();
531 [ + - ][ + - ]: 99 : rval = treeTool->moab()->get_child_meshsets( *i, childVect );
[ + - ]
532 [ - + ]: 99 : if( MB_SUCCESS != rval ) return rval;
533 [ - + ]: 99 : if( childVect.size() != 2 ) return MB_FAILURE;
534 : :
535 [ + - ]: 99 : ++i;
536 [ + - ][ + - ]: 99 : if( childVect[1] == *i ) plane.flip();
[ + + ][ + - ]
537 : :
538 [ + - ]: 99 : CartVect norm( plane.norm );
539 [ + - ]: 99 : poly_out.cut_polyhedron( norm, plane.coeff );
540 : : }
541 : :
542 : 24 : return MB_SUCCESS;
543 : : }
544 : :
545 : 19 : ErrorCode BSPTreeBoxIter::initialize( BSPTree* tool_ptr, EntityHandle root, const double* point )
546 : : {
547 : 19 : ErrorCode rval = BSPTreeIter::initialize( tool_ptr, root );
548 [ - + ]: 19 : if( MB_SUCCESS != rval ) return rval;
549 : :
550 : 19 : rval = tool()->get_tree_box( root, leafCoords );
551 [ - + ]: 19 : if( MB_SUCCESS != rval ) return rval;
552 : :
553 [ + + ][ + + ]: 19 : if( point && !point_in_box( leafCoords, point ) ) return MB_ENTITY_NOT_FOUND;
[ + + ]
554 : :
555 : 17 : stackData.resize( 1 );
556 : 17 : return MB_SUCCESS;
557 : : }
558 : :
559 : 367 : BSPTreeBoxIter::SideBits BSPTreeBoxIter::side_above_plane( const double hex_coords[8][3], const BSPTree::Plane& plane )
560 : : {
561 : 367 : unsigned result = 0;
562 [ + + ]: 3303 : for( unsigned i = 0; i < 8u; ++i )
563 : 2936 : result |= plane.above( hex_coords[i] ) << i;
564 : 367 : return (BSPTreeBoxIter::SideBits)result;
565 : : }
566 : :
567 : 210 : BSPTreeBoxIter::SideBits BSPTreeBoxIter::side_on_plane( const double hex_coords[8][3], const BSPTree::Plane& plane )
568 : : {
569 : 210 : unsigned result = 0;
570 [ + + ]: 1890 : for( unsigned i = 0; i < 8u; ++i )
571 : : {
572 : 1680 : bool on = plane.distance( hex_coords[i] ) <= BSPTree::epsilon();
573 : 1680 : result |= on << i;
574 : : }
575 : 210 : return (BSPTreeBoxIter::SideBits)result;
576 : : }
577 : :
578 : 272 : static inline void copy_coords( const double src[3], double dest[3] )
579 : : {
580 : 272 : dest[0] = src[0];
581 : 272 : dest[1] = src[1];
582 : 272 : dest[2] = src[2];
583 : 272 : }
584 : :
585 : 68 : ErrorCode BSPTreeBoxIter::face_corners( const SideBits face, const double hex_corners[8][3], double face_corners[4][3] )
586 : : {
587 [ + + + + : 68 : switch( face )
+ + - ]
588 : : {
589 : : case BSPTreeBoxIter::B0154:
590 : 12 : copy_coords( hex_corners[0], face_corners[0] );
591 : 12 : copy_coords( hex_corners[1], face_corners[1] );
592 : 12 : copy_coords( hex_corners[5], face_corners[2] );
593 : 12 : copy_coords( hex_corners[4], face_corners[3] );
594 : 12 : break;
595 : : case BSPTreeBoxIter::B1265:
596 : 14 : copy_coords( hex_corners[1], face_corners[0] );
597 : 14 : copy_coords( hex_corners[2], face_corners[1] );
598 : 14 : copy_coords( hex_corners[6], face_corners[2] );
599 : 14 : copy_coords( hex_corners[5], face_corners[3] );
600 : 14 : break;
601 : : case BSPTreeBoxIter::B2376:
602 : 12 : copy_coords( hex_corners[2], face_corners[0] );
603 : 12 : copy_coords( hex_corners[3], face_corners[1] );
604 : 12 : copy_coords( hex_corners[7], face_corners[2] );
605 : 12 : copy_coords( hex_corners[6], face_corners[3] );
606 : 12 : break;
607 : : case BSPTreeBoxIter::B3047:
608 : 14 : copy_coords( hex_corners[3], face_corners[0] );
609 : 14 : copy_coords( hex_corners[0], face_corners[1] );
610 : 14 : copy_coords( hex_corners[4], face_corners[2] );
611 : 14 : copy_coords( hex_corners[7], face_corners[3] );
612 : 14 : break;
613 : : case BSPTreeBoxIter::B3210:
614 : 8 : copy_coords( hex_corners[3], face_corners[0] );
615 : 8 : copy_coords( hex_corners[2], face_corners[1] );
616 : 8 : copy_coords( hex_corners[1], face_corners[2] );
617 : 8 : copy_coords( hex_corners[0], face_corners[3] );
618 : 8 : break;
619 : : case BSPTreeBoxIter::B4567:
620 : 8 : copy_coords( hex_corners[4], face_corners[0] );
621 : 8 : copy_coords( hex_corners[5], face_corners[1] );
622 : 8 : copy_coords( hex_corners[6], face_corners[2] );
623 : 8 : copy_coords( hex_corners[7], face_corners[3] );
624 : 8 : break;
625 : : default:
626 : 0 : return MB_FAILURE; // child is not a box
627 : : }
628 : :
629 : 68 : return MB_SUCCESS;
630 : : }
631 : :
632 : : /** \brief Clip an edge using a plane
633 : : *
634 : : * Given an edge from keep_end_coords to cut_end_coords,
635 : : * cut the edge using the passed plane, such that cut_end_coords
636 : : * is updated with a new location on the plane, and old_coords_out
637 : : * contains the original value of cut_end_coords.
638 : : */
639 : 676 : static inline void plane_cut_edge( double old_coords_out[3], const double keep_end_coords[3], double cut_end_coords[3],
640 : : const BSPTree::Plane& plane )
641 : : {
642 [ + - ][ + - ]: 676 : const CartVect start( keep_end_coords ), end( cut_end_coords );
643 [ + - ]: 676 : const CartVect norm( plane.norm );
644 [ + - ]: 676 : CartVect xsect_point;
645 : :
646 [ + - ]: 676 : const CartVect m = end - start;
647 [ + - ][ + - ]: 676 : const double t = -( norm % start + plane.coeff ) / ( norm % m );
648 [ + - ][ - + ]: 676 : assert( t > 0.0 && t < 1.0 );
649 [ + - ][ + - ]: 676 : xsect_point = start + t * m;
650 : :
651 [ + - ]: 676 : end.get( old_coords_out );
652 [ + - ]: 676 : xsect_point.get( cut_end_coords );
653 : 676 : }
654 : :
655 : : /** Given the corners of a hexahedron in corners_input and a
656 : : * plane, cut the hex with the plane, updating corners_input
657 : : * and storing the original,cut-off side of the hex in cut_face_out.
658 : : *
659 : : * The portion of the hex below the plane is retained. cut_face_out
660 : : * will contain the side of the hex that is entirely above the plane.
661 : : *\return MB_FAILURE if plane/hex intersection is not a quadrilateral.
662 : : */
663 : 211 : static ErrorCode plane_cut_box( double cut_face_out[4][3], double corners_inout[8][3], const BSPTree::Plane& plane )
664 : : {
665 [ + + + + : 211 : switch( BSPTreeBoxIter::side_above_plane( corners_inout, plane ) )
+ + + ]
666 : : {
667 : : case BSPTreeBoxIter::B0154:
668 : 39 : plane_cut_edge( cut_face_out[0], corners_inout[3], corners_inout[0], plane );
669 : 39 : plane_cut_edge( cut_face_out[1], corners_inout[2], corners_inout[1], plane );
670 : 39 : plane_cut_edge( cut_face_out[2], corners_inout[6], corners_inout[5], plane );
671 : 39 : plane_cut_edge( cut_face_out[3], corners_inout[7], corners_inout[4], plane );
672 : 39 : break;
673 : : case BSPTreeBoxIter::B1265:
674 : 44 : plane_cut_edge( cut_face_out[0], corners_inout[0], corners_inout[1], plane );
675 : 44 : plane_cut_edge( cut_face_out[1], corners_inout[3], corners_inout[2], plane );
676 : 44 : plane_cut_edge( cut_face_out[2], corners_inout[7], corners_inout[6], plane );
677 : 44 : plane_cut_edge( cut_face_out[3], corners_inout[4], corners_inout[5], plane );
678 : 44 : break;
679 : : case BSPTreeBoxIter::B2376:
680 : 38 : plane_cut_edge( cut_face_out[0], corners_inout[1], corners_inout[2], plane );
681 : 38 : plane_cut_edge( cut_face_out[1], corners_inout[0], corners_inout[3], plane );
682 : 38 : plane_cut_edge( cut_face_out[2], corners_inout[4], corners_inout[7], plane );
683 : 38 : plane_cut_edge( cut_face_out[3], corners_inout[5], corners_inout[6], plane );
684 : 38 : break;
685 : : case BSPTreeBoxIter::B3047:
686 : 41 : plane_cut_edge( cut_face_out[0], corners_inout[2], corners_inout[3], plane );
687 : 41 : plane_cut_edge( cut_face_out[1], corners_inout[1], corners_inout[0], plane );
688 : 41 : plane_cut_edge( cut_face_out[2], corners_inout[5], corners_inout[4], plane );
689 : 41 : plane_cut_edge( cut_face_out[3], corners_inout[6], corners_inout[7], plane );
690 : 41 : break;
691 : : case BSPTreeBoxIter::B3210:
692 : 3 : plane_cut_edge( cut_face_out[0], corners_inout[7], corners_inout[3], plane );
693 : 3 : plane_cut_edge( cut_face_out[1], corners_inout[6], corners_inout[2], plane );
694 : 3 : plane_cut_edge( cut_face_out[2], corners_inout[5], corners_inout[1], plane );
695 : 3 : plane_cut_edge( cut_face_out[3], corners_inout[4], corners_inout[0], plane );
696 : 3 : break;
697 : : case BSPTreeBoxIter::B4567:
698 : 4 : plane_cut_edge( cut_face_out[0], corners_inout[0], corners_inout[4], plane );
699 : 4 : plane_cut_edge( cut_face_out[1], corners_inout[1], corners_inout[5], plane );
700 : 4 : plane_cut_edge( cut_face_out[2], corners_inout[2], corners_inout[6], plane );
701 : 4 : plane_cut_edge( cut_face_out[3], corners_inout[3], corners_inout[7], plane );
702 : 4 : break;
703 : : default:
704 : 42 : return MB_FAILURE; // child is not a box
705 : : }
706 : :
707 : 169 : return MB_SUCCESS;
708 : : }
709 : :
710 : 840 : static inline void copy_coords( double dest[3], const double source[3] )
711 : : {
712 : 840 : dest[0] = source[0];
713 : 840 : dest[1] = source[1];
714 : 840 : dest[2] = source[2];
715 : 840 : }
716 : :
717 : : /** reverse of plane_cut_box */
718 : 210 : static inline ErrorCode plane_uncut_box( const double cut_face_in[4][3], double corners_inout[8][3],
719 : : const BSPTree::Plane& plane )
720 : : {
721 [ + + + + : 210 : switch( BSPTreeBoxIter::side_on_plane( corners_inout, plane ) )
+ + - ]
722 : : {
723 : : case BSPTreeBoxIter::B0154:
724 : 43 : copy_coords( corners_inout[0], cut_face_in[0] );
725 : 43 : copy_coords( corners_inout[1], cut_face_in[1] );
726 : 43 : copy_coords( corners_inout[5], cut_face_in[2] );
727 : 43 : copy_coords( corners_inout[4], cut_face_in[3] );
728 : 43 : break;
729 : : case BSPTreeBoxIter::B1265:
730 : 62 : copy_coords( corners_inout[1], cut_face_in[0] );
731 : 62 : copy_coords( corners_inout[2], cut_face_in[1] );
732 : 62 : copy_coords( corners_inout[6], cut_face_in[2] );
733 : 62 : copy_coords( corners_inout[5], cut_face_in[3] );
734 : 62 : break;
735 : : case BSPTreeBoxIter::B2376:
736 : 42 : copy_coords( corners_inout[2], cut_face_in[0] );
737 : 42 : copy_coords( corners_inout[3], cut_face_in[1] );
738 : 42 : copy_coords( corners_inout[7], cut_face_in[2] );
739 : 42 : copy_coords( corners_inout[6], cut_face_in[3] );
740 : 42 : break;
741 : : case BSPTreeBoxIter::B3047:
742 : 58 : copy_coords( corners_inout[3], cut_face_in[0] );
743 : 58 : copy_coords( corners_inout[0], cut_face_in[1] );
744 : 58 : copy_coords( corners_inout[4], cut_face_in[2] );
745 : 58 : copy_coords( corners_inout[7], cut_face_in[3] );
746 : 58 : break;
747 : : case BSPTreeBoxIter::B3210:
748 : 3 : copy_coords( corners_inout[3], cut_face_in[0] );
749 : 3 : copy_coords( corners_inout[2], cut_face_in[1] );
750 : 3 : copy_coords( corners_inout[1], cut_face_in[2] );
751 : 3 : copy_coords( corners_inout[0], cut_face_in[3] );
752 : 3 : break;
753 : : case BSPTreeBoxIter::B4567:
754 : 2 : copy_coords( corners_inout[4], cut_face_in[0] );
755 : 2 : copy_coords( corners_inout[5], cut_face_in[1] );
756 : 2 : copy_coords( corners_inout[6], cut_face_in[2] );
757 : 2 : copy_coords( corners_inout[7], cut_face_in[3] );
758 : 2 : break;
759 : : default:
760 : 0 : return MB_FAILURE; // child is not a box
761 : : }
762 : :
763 : 210 : return MB_SUCCESS;
764 : : }
765 : :
766 : 111 : ErrorCode BSPTreeBoxIter::step_to_first_leaf( Direction direction )
767 : : {
768 : : ErrorCode rval;
769 [ + - ]: 111 : BSPTree::Plane plane;
770 : : Corners clipped_corners;
771 : :
772 : : for( ;; )
773 : : {
774 : 152 : childVect.clear();
775 [ + - ][ + - ]: 152 : rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
[ + - ][ + - ]
776 [ - + ]: 152 : if( MB_SUCCESS != rval ) return rval;
777 [ + + ]: 152 : if( childVect.empty() ) // leaf
778 : 69 : break;
779 : :
780 [ + - ][ + - ]: 83 : rval = tool()->get_split_plane( mStack.back(), plane );
[ + - ]
781 [ - + ]: 83 : if( MB_SUCCESS != rval ) return rval;
782 : :
783 [ - + ][ # # ]: 83 : if( direction == RIGHT ) plane.flip();
784 [ + - ]: 83 : rval = plane_cut_box( clipped_corners.coords, leafCoords, plane );
785 [ + + ]: 83 : if( MB_SUCCESS != rval ) return rval;
786 [ + - ][ + - ]: 41 : mStack.push_back( childVect[direction] );
787 [ + - ]: 41 : stackData.push_back( clipped_corners );
788 : : }
789 : 152 : return MB_SUCCESS;
790 : : }
791 : :
792 : 186 : ErrorCode BSPTreeBoxIter::up()
793 : : {
794 : : ErrorCode rval;
795 [ + + ]: 186 : if( mStack.size() == 1 ) return MB_ENTITY_NOT_FOUND;
796 : :
797 [ + - ]: 158 : EntityHandle node = mStack.back();
798 [ + - ]: 158 : Corners clipped_face = stackData.back();
799 [ + - ]: 158 : mStack.pop_back();
800 [ + - ]: 158 : stackData.pop_back();
801 : :
802 [ + - ]: 158 : BSPTree::Plane plane;
803 [ + - ][ + - ]: 158 : rval = tool()->get_split_plane( mStack.back(), plane );
[ + - ]
804 [ - + ]: 158 : if( MB_SUCCESS != rval )
805 : : {
806 [ # # ]: 0 : mStack.push_back( node );
807 [ # # ]: 0 : stackData.push_back( clipped_face );
808 : 0 : return rval;
809 : : }
810 : :
811 [ + - ]: 158 : rval = plane_uncut_box( clipped_face.coords, leafCoords, plane );
812 [ - + ]: 158 : if( MB_SUCCESS != rval )
813 : : {
814 [ # # ]: 0 : mStack.push_back( node );
815 [ # # ]: 0 : stackData.push_back( clipped_face );
816 : 0 : return rval;
817 : : }
818 : :
819 : 186 : return MB_SUCCESS;
820 : : }
821 : :
822 : 94 : ErrorCode BSPTreeBoxIter::down( const BSPTree::Plane& plane_ref, Direction direction )
823 : : {
824 : 94 : childVect.clear();
825 [ + - ][ + - ]: 94 : ErrorCode rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
[ + - ][ + - ]
826 [ - + ]: 94 : if( MB_SUCCESS != rval ) return rval;
827 [ - + ]: 94 : if( childVect.empty() ) return MB_ENTITY_NOT_FOUND;
828 : :
829 : 94 : BSPTree::Plane plane( plane_ref );
830 [ + + ][ + - ]: 94 : if( direction == RIGHT ) plane.flip();
831 : :
832 : : Corners clipped_face;
833 [ + - ]: 94 : rval = plane_cut_box( clipped_face.coords, leafCoords, plane );
834 [ - + ]: 94 : if( MB_SUCCESS != rval ) return rval;
835 : :
836 [ + - ][ + - ]: 94 : mStack.push_back( childVect[direction] );
837 [ + - ]: 94 : stackData.push_back( clipped_face );
838 : 94 : return MB_SUCCESS;
839 : : }
840 : :
841 : 38 : ErrorCode BSPTreeBoxIter::step( Direction direction )
842 : : {
843 : : EntityHandle node, parent;
844 : : Corners clipped_face;
845 : : ErrorCode rval;
846 [ + - ]: 38 : BSPTree::Plane plane;
847 : 38 : const Direction opposite = static_cast< Direction >( 1 - direction );
848 : :
849 : : // If stack is empty, then either this iterator is uninitialized
850 : : // or we reached the end of the iteration (and return
851 : : // MB_ENTITY_NOT_FOUND) already.
852 [ - + ]: 38 : if( mStack.empty() ) return MB_FAILURE;
853 : :
854 : : // Pop the current node from the stack.
855 : : // The stack should then contain the parent of the current node.
856 : : // If the stack is empty after this pop, then we've reached the end.
857 [ + - ]: 38 : node = mStack.back();
858 [ + - ]: 38 : mStack.pop_back();
859 [ + - ]: 38 : clipped_face = stackData.back();
860 [ + - ]: 38 : stackData.pop_back();
861 : :
862 [ + + ]: 56 : while( !mStack.empty() )
863 : : {
864 : : // Get data for parent entity
865 [ + - ]: 52 : parent = mStack.back();
866 : 52 : childVect.clear();
867 [ + - ][ + - ]: 52 : rval = tool()->moab()->get_child_meshsets( parent, childVect );
[ + - ]
868 [ - + ]: 52 : if( MB_SUCCESS != rval ) return rval;
869 [ + - ][ + - ]: 52 : rval = tool()->get_split_plane( parent, plane );
870 [ - + ]: 52 : if( MB_SUCCESS != rval ) return rval;
871 [ + + ][ + - ]: 52 : if( direction == LEFT ) plane.flip();
872 : :
873 : : // If we're at the left child
874 [ + - ][ + + ]: 52 : if( childVect[opposite] == node )
875 : : {
876 : : // change from box of left child to box of parent
877 [ + - ]: 34 : plane_uncut_box( clipped_face.coords, leafCoords, plane );
878 : : // change from box of parent to box of right child
879 [ + - ]: 34 : plane.flip();
880 [ + - ]: 34 : plane_cut_box( clipped_face.coords, leafCoords, plane );
881 : : // push right child on stack
882 [ + - ][ + - ]: 34 : mStack.push_back( childVect[direction] );
883 [ + - ]: 34 : stackData.push_back( clipped_face );
884 : : // descend to left-most leaf of the right child
885 [ + - ]: 34 : return step_to_first_leaf( opposite );
886 : : }
887 : :
888 : : // The current node is the right child of the parent,
889 : : // continue up the tree.
890 [ + - ][ - + ]: 18 : assert( childVect[direction] == node );
891 [ + - ]: 18 : plane.flip();
892 [ + - ]: 18 : plane_uncut_box( clipped_face.coords, leafCoords, plane );
893 : 18 : node = parent;
894 [ + - ]: 18 : clipped_face = stackData.back();
895 [ + - ]: 18 : mStack.pop_back();
896 [ + - ]: 18 : stackData.pop_back();
897 : : }
898 : :
899 : 38 : return MB_ENTITY_NOT_FOUND;
900 : : }
901 : :
902 : 26 : ErrorCode BSPTreeBoxIter::get_box_corners( double coords[8][3] ) const
903 : : {
904 : 26 : memcpy( coords, leafCoords, 24 * sizeof( double ) );
905 : 26 : return MB_SUCCESS;
906 : : }
907 : :
908 : : // result = a - b
909 : 12 : static void subtr( double result[3], const double a[3], const double b[3] )
910 : : {
911 : 12 : result[0] = a[0] - b[0];
912 : 12 : result[1] = a[1] - b[1];
913 : 12 : result[2] = a[2] - b[2];
914 : 12 : }
915 : :
916 : : // result = a + b + c + d
917 : 24 : static void sum( double result[3], const double a[3], const double b[3], const double c[3], const double d[3] )
918 : : {
919 : 24 : result[0] = a[0] + b[0] + c[0] + d[0];
920 : 24 : result[1] = a[1] + b[1] + c[1] + d[1];
921 : 24 : result[2] = a[2] + b[2] + c[2] + d[2];
922 : 24 : }
923 : :
924 : : // result = a cross b
925 : 4 : static void cross( double result[3], const double a[3], const double b[3] )
926 : : {
927 : 4 : result[0] = a[1] * b[2] - a[2] * b[1];
928 : 4 : result[1] = a[2] * b[0] - a[0] * b[2];
929 : 4 : result[2] = a[0] * b[1] - a[1] * b[0];
930 : 4 : }
931 : :
932 : 4 : static double dot( const double a[3], const double b[3] )
933 : : {
934 : 4 : return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
935 : : }
936 : :
937 : 4 : double BSPTreeBoxIter::volume() const
938 : : {
939 : : // have planar sides, so use mid-face tripple product
940 : : double f1[3], f2[3], f3[3], f4[3], f5[3], f6[3];
941 : 4 : sum( f1, leafCoords[0], leafCoords[1], leafCoords[4], leafCoords[5] );
942 : 4 : sum( f2, leafCoords[1], leafCoords[2], leafCoords[5], leafCoords[6] );
943 : 4 : sum( f3, leafCoords[2], leafCoords[3], leafCoords[6], leafCoords[7] );
944 : 4 : sum( f4, leafCoords[0], leafCoords[3], leafCoords[4], leafCoords[7] );
945 : 4 : sum( f5, leafCoords[0], leafCoords[1], leafCoords[2], leafCoords[3] );
946 : 4 : sum( f6, leafCoords[4], leafCoords[5], leafCoords[6], leafCoords[7] );
947 : : double v13[3], v24[3], v65[3];
948 : 4 : subtr( v13, f1, f3 );
949 : 4 : subtr( v24, f2, f4 );
950 : 4 : subtr( v65, f6, f5 );
951 : : double cr[3];
952 : 4 : cross( cr, v13, v24 );
953 : 4 : return ( 1. / 64 ) * dot( cr, v65 );
954 : : }
955 : :
956 : 16 : BSPTreeBoxIter::XSect BSPTreeBoxIter::splits( const BSPTree::Plane& plane ) const
957 : : {
958 : : // test each corner relative to the plane
959 : 16 : unsigned result = 0;
960 [ + + ]: 114 : for( unsigned i = 0; i < 8u; ++i )
961 : : {
962 : 102 : double d = plane.signed_distance( leafCoords[i] );
963 : : // if corner is on plane, than intersection
964 : : // will result in a degenerate hex
965 [ + + ]: 102 : if( fabs( d ) < BSPTree::epsilon() ) return NONHEX;
966 : : // if mark vertices above plane
967 [ + + ]: 98 : if( d > 0.0 ) result |= 1 << i;
968 : : }
969 : :
970 [ + + + ]: 12 : switch( result )
971 : : {
972 : : // if all vertices or no vertices above plane,
973 : : // then plane doesn't intersect
974 : : case 0:
975 : : case 0xFF:
976 : 4 : return MISS;
977 : :
978 : : // if there are four vertices above the plane
979 : : // and they compose a single face of the hex,
980 : : // then the cut will result in two hexes
981 : : case B0154:
982 : : case B1265:
983 : : case B2376:
984 : : case B3047:
985 : : case B3210:
986 : : case B4567:
987 : 4 : return SPLIT;
988 : :
989 : : // otherwise intersects, but split would not result
990 : : // in two hexahedrons
991 : : default:
992 : 4 : return NONHEX;
993 : : }
994 : : }
995 : :
996 : 12 : bool BSPTreeBoxIter::intersects( const BSPTree::Plane& plane ) const
997 : : {
998 : : // test each corner relative to the plane
999 : 12 : unsigned count = 0;
1000 [ + + ]: 108 : for( unsigned i = 0; i < 8u; ++i )
1001 : 96 : count += plane.above( leafCoords[i] );
1002 [ + + ][ + + ]: 12 : return count > 0 && count < 8u;
1003 : : }
1004 : :
1005 : 0 : ErrorCode BSPTreeBoxIter::sibling_side( SideBits& side_out ) const
1006 : : {
1007 [ # # ]: 0 : if( mStack.size() < 2 ) // at tree root
1008 : 0 : return MB_ENTITY_NOT_FOUND;
1009 : :
1010 [ # # ]: 0 : EntityHandle parent = mStack[mStack.size() - 2];
1011 [ # # ]: 0 : BSPTree::Plane plane;
1012 [ # # ][ # # ]: 0 : ErrorCode rval = tool()->get_split_plane( parent, plane );
1013 [ # # ]: 0 : if( MB_SUCCESS != rval ) return MB_FAILURE;
1014 : :
1015 [ # # ]: 0 : side_out = side_on_plane( leafCoords, plane );
1016 : 0 : return MB_SUCCESS;
1017 : : }
1018 : :
1019 : 68 : ErrorCode BSPTreeBoxIter::get_neighbors( SideBits side, std::vector< BSPTreeBoxIter >& results, double epsilon ) const
1020 : : {
1021 : : EntityHandle tmp_handle;
1022 [ + - ]: 68 : BSPTree::Plane plane;
1023 : : ErrorCode rval;
1024 : : int n;
1025 : :
1026 : : Corners face;
1027 [ + - ]: 68 : rval = face_corners( side, leafCoords, face.coords );
1028 [ - + ]: 68 : if( MB_SUCCESS != rval ) return rval;
1029 : :
1030 : : // Move up tree until we find the split that created the specified side.
1031 : : // Push the sibling at that level onto the iterator stack as
1032 : : // all neighbors will be rooted at that node.
1033 [ + - ]: 68 : BSPTreeBoxIter iter( *this ); // temporary iterator (don't modifiy *this)
1034 : : for( ;; )
1035 : : {
1036 [ + - ]: 184 : tmp_handle = iter.handle();
1037 : :
1038 [ + - ]: 184 : rval = iter.up();
1039 [ + + ]: 184 : if( MB_SUCCESS != rval ) // reached root - no neighbors on that side
1040 [ - + ]: 28 : return ( rval == MB_ENTITY_NOT_FOUND ) ? MB_SUCCESS : rval;
1041 : :
1042 : 156 : iter.childVect.clear();
1043 [ + - ][ + - ]: 156 : rval = tool()->moab()->get_child_meshsets( iter.handle(), iter.childVect );
[ + - ][ + - ]
1044 [ - + ]: 156 : if( MB_SUCCESS != rval ) return rval;
1045 : :
1046 [ + - ][ + - ]: 156 : rval = tool()->get_split_plane( iter.handle(), plane );
[ + - ]
1047 [ - + ]: 156 : if( MB_SUCCESS != rval ) return rval;
1048 [ + - ]: 156 : SideBits s = side_above_plane( iter.leafCoords, plane );
1049 : :
1050 [ + - ][ + + ]: 156 : if( tmp_handle == iter.childVect[0] && s == side )
[ + + ][ + + ]
1051 : : {
1052 [ + - ]: 20 : rval = iter.down( plane, RIGHT );
1053 [ - + ]: 20 : if( MB_SUCCESS != rval ) return rval;
1054 : 20 : break;
1055 : : }
1056 [ + - ][ + + ]: 136 : else if( tmp_handle == iter.childVect[1] && opposite_face( s ) == side )
[ + - ][ + + ]
[ + + ]
1057 : : {
1058 [ + - ]: 20 : rval = iter.down( plane, LEFT );
1059 [ - + ]: 20 : if( MB_SUCCESS != rval ) return rval;
1060 : 20 : break;
1061 : : }
1062 : 116 : }
1063 : :
1064 : : // now move down tree, searching for adjacent boxes
1065 [ + - ]: 80 : std::vector< BSPTreeBoxIter > list;
1066 : : // loop over all potential paths to neighbors (until list is empty)
1067 : : for( ;; )
1068 : : {
1069 : : // follow a single path to a leaf, append any other potential
1070 : : // paths to neighbors to 'list'
1071 : : for( ;; )
1072 : : {
1073 [ + - ][ + - ]: 86 : rval = tool()->moab()->num_child_meshsets( iter.handle(), &n );
[ + - ][ + - ]
1074 [ - + ]: 86 : if( MB_SUCCESS != rval ) return rval;
1075 : :
1076 : : // if leaf
1077 [ + + ]: 86 : if( !n )
1078 : : {
1079 [ + - ]: 52 : results.push_back( iter );
1080 : 52 : break;
1081 : : }
1082 : :
1083 [ + - ][ + - ]: 34 : rval = tool()->get_split_plane( iter.handle(), plane );
[ + - ]
1084 [ - + ]: 34 : if( MB_SUCCESS != rval ) return rval;
1085 : :
1086 : 34 : bool some_above = false, some_below = false;
1087 [ + + ]: 170 : for( int i = 0; i < 4; ++i )
1088 : : {
1089 [ + - ]: 136 : double signed_d = plane.signed_distance( face.coords[i] );
1090 [ + + ]: 136 : if( signed_d > -epsilon ) some_above = true;
1091 [ + + ]: 136 : if( signed_d < epsilon ) some_below = true;
1092 : : }
1093 : :
1094 [ + + ][ + + ]: 34 : if( some_above && some_below )
1095 : : {
1096 [ + - ]: 12 : list.push_back( iter );
1097 [ + - ][ + - ]: 12 : list.back().down( plane, RIGHT );
1098 [ + - ]: 12 : iter.down( plane, LEFT );
1099 : : }
1100 [ + + ]: 22 : else if( some_above )
1101 : : {
1102 [ + - ]: 9 : iter.down( plane, RIGHT );
1103 : : }
1104 [ + - ]: 13 : else if( some_below )
1105 : : {
1106 [ + - ]: 13 : iter.down( plane, LEFT );
1107 : : }
1108 : : else
1109 : : {
1110 : : // tolerance issue -- epsilon to small? 2D box?
1111 : 0 : return MB_FAILURE;
1112 : : }
1113 : 34 : }
1114 : :
1115 [ + + ]: 52 : if( list.empty() ) break;
1116 : :
1117 [ + - ][ + - ]: 12 : iter = list.back();
1118 [ + - ]: 12 : list.pop_back();
1119 : 12 : }
1120 : :
1121 : 108 : return MB_SUCCESS;
1122 : : }
1123 : :
1124 : 0 : ErrorCode BSPTreeBoxIter::calculate_polyhedron( BSPTreePoly& poly_out ) const
1125 : : {
1126 : 0 : const CartVect* ptr = reinterpret_cast< const CartVect* >( leafCoords );
1127 : 0 : return poly_out.set( ptr );
1128 : : }
1129 : :
1130 : 28 : ErrorCode BSPTree::leaf_containing_point( EntityHandle tree_root, const double point[3], EntityHandle& leaf_out )
1131 : : {
1132 [ + - ]: 28 : std::vector< EntityHandle > children;
1133 [ + - ]: 28 : Plane plane;
1134 : 28 : EntityHandle node = tree_root;
1135 [ + - ][ + - ]: 28 : ErrorCode rval = moab()->get_child_meshsets( node, children );
1136 [ - + ]: 28 : if( MB_SUCCESS != rval ) return rval;
1137 [ + + ]: 135 : while( !children.empty() )
1138 : : {
1139 [ + - ]: 107 : rval = get_split_plane( node, plane );
1140 [ - + ]: 107 : if( MB_SUCCESS != rval ) return rval;
1141 : :
1142 [ + - ][ + - ]: 107 : node = children[plane.above( point )];
1143 : 107 : children.clear();
1144 [ + - ][ + - ]: 107 : rval = moab()->get_child_meshsets( node, children );
1145 [ - + ]: 107 : if( MB_SUCCESS != rval ) return rval;
1146 : : }
1147 : 28 : leaf_out = node;
1148 : 28 : return MB_SUCCESS;
1149 : : }
1150 : :
1151 : 14 : ErrorCode BSPTree::leaf_containing_point( EntityHandle root, const double point[3], BSPTreeIter& iter )
1152 : : {
1153 : : ErrorCode rval;
1154 : :
1155 : 14 : rval = iter.initialize( this, root, point );
1156 [ + + ]: 14 : if( MB_SUCCESS != rval ) return rval;
1157 : :
1158 : : for( ;; )
1159 : : {
1160 : 36 : iter.childVect.clear();
1161 [ + - ][ + - ]: 36 : rval = moab()->get_child_meshsets( iter.handle(), iter.childVect );
[ + - ]
1162 [ + - ][ + + ]: 38 : if( MB_SUCCESS != rval || iter.childVect.empty() ) return rval;
[ + + ]
1163 : :
1164 [ + - ]: 24 : Plane plane;
1165 [ + - ][ + - ]: 24 : rval = get_split_plane( iter.handle(), plane );
1166 [ - + ]: 24 : if( MB_SUCCESS != rval ) return rval;
1167 : :
1168 [ + - ][ + - ]: 24 : rval = iter.down( plane, ( BSPTreeIter::Direction )( plane.above( point ) ) );
1169 [ - + ]: 24 : if( MB_SUCCESS != rval ) return rval;
1170 : 24 : }
1171 : : }
1172 : :
1173 : : template < typename PlaneIter >
1174 : 159 : static inline bool ray_intersect_halfspaces( const CartVect& ray_pt, const CartVect& ray_dir, const PlaneIter& begin,
1175 : : const PlaneIter& end, double& t_enter, double& t_exit )
1176 : : {
1177 : 159 : const double epsilon = 1e-12;
1178 : :
1179 : : // begin with inifinite ray
1180 : 159 : t_enter = 0.0;
1181 : 159 : t_exit = std::numeric_limits< double >::infinity();
1182 : :
1183 : : // cut ray such that we keep only the portion contained
1184 : : // in each halfspace
1185 [ + - ][ + - ]: 777 : for( PlaneIter i = begin; i != end; ++i )
[ + - ][ + + ]
[ + - ][ + - ]
[ + - ][ + + ]
1186 : : {
1187 [ + - ][ + - ]: 642 : CartVect norm( i->norm );
[ + - ][ + - ]
1188 [ + - ][ + - ]: 642 : double coeff = i->coeff;
1189 [ + - ][ + - ]: 642 : double den = norm % ray_dir;
1190 [ + + ][ + + ]: 642 : if( fabs( den ) < epsilon )
1191 : : { // ray is parallel to plane
1192 [ + - ][ + - ]: 276 : if( i->above( ray_pt.array() ) ) return false; // ray entirely outside half-space
[ + - ][ - + ]
[ + - ][ + - ]
[ + - ][ + + ]
1193 : : }
1194 : : else
1195 : : {
1196 [ + - ][ + - ]: 366 : double t_xsect = ( -coeff - ( norm % ray_pt ) ) / den;
1197 : : // keep portion of ray/segment below plane
1198 [ + + ][ + + ]: 366 : if( den > 0 )
1199 : : {
1200 [ + - ][ + + ]: 183 : if( t_xsect < t_exit ) t_exit = t_xsect;
1201 : : }
1202 : : else
1203 : : {
1204 [ - + ][ + + ]: 183 : if( t_xsect > t_enter ) t_enter = t_xsect;
1205 : : }
1206 : : }
1207 : : }
1208 : :
1209 : 159 : return t_exit >= t_enter;
1210 : : }
1211 : :
1212 : : class BoxPlaneIter
1213 : : {
1214 : : int faceNum;
1215 : : BSPTree::Plane facePlanes[6];
1216 : :
1217 : : public:
1218 : : BoxPlaneIter( const double coords[8][3] );
1219 [ + + ]: 848 : BoxPlaneIter() : faceNum( 6 ) {} // initialize to 'end'
1220 : 1428 : const BSPTree::Plane* operator->() const
1221 : : {
1222 : 1428 : return facePlanes + faceNum;
1223 : : }
1224 : : bool operator==( const BoxPlaneIter& other ) const
1225 : : {
1226 : : return faceNum == other.faceNum;
1227 : : }
1228 : 670 : bool operator!=( const BoxPlaneIter& other ) const
1229 : : {
1230 : 670 : return faceNum != other.faceNum;
1231 : : }
1232 : 564 : BoxPlaneIter& operator++()
1233 : : {
1234 : 564 : ++faceNum;
1235 : 564 : return *this;
1236 : : }
1237 : : };
1238 : :
1239 : : static const int box_face_corners[6][4] = { { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 3, 7, 6 },
1240 : : { 3, 0, 4, 7 }, { 3, 2, 1, 0 }, { 4, 5, 6, 7 } };
1241 : :
1242 [ + + ]: 848 : BoxPlaneIter::BoxPlaneIter( const double coords[8][3] ) : faceNum( 0 )
1243 : : {
1244 : : // NOTE: In the case of a BSP tree, all sides of the
1245 : : // leaf will planar.
1246 : : assert( sizeof( CartVect ) == sizeof( coords[0] ) );
1247 : 106 : const CartVect* corners = reinterpret_cast< const CartVect* >( coords );
1248 [ + + ]: 742 : for( int i = 0; i < 6; ++i )
1249 : : {
1250 : 636 : const int* indices = box_face_corners[i];
1251 [ + - ]: 636 : CartVect v1 = corners[indices[1]] - corners[indices[0]];
1252 [ + - ]: 636 : CartVect v2 = corners[indices[3]] - corners[indices[0]];
1253 [ + - ]: 636 : CartVect n = v1 * v2;
1254 [ + - ][ + - ]: 636 : facePlanes[i] = BSPTree::Plane( n.array(), -( n % corners[indices[2]] ) );
[ + - ]
1255 : : }
1256 : 106 : }
1257 : :
1258 : 53 : bool BSPTreeBoxIter::intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter,
1259 : : double& t_exit ) const
1260 : : {
1261 [ + - ][ + - ]: 53 : BoxPlaneIter iter( this->leafCoords ), end;
1262 [ + - ][ + - ]: 53 : return ray_intersect_halfspaces( CartVect( ray_point ), CartVect( ray_vect ), iter, end, t_enter, t_exit );
[ + - ]
1263 : : }
1264 : :
1265 : 424 : class BSPTreePlaneIter
1266 : : {
1267 : : BSPTree* toolPtr;
1268 : : const EntityHandle* const pathToRoot;
1269 : : int pathPos;
1270 : : BSPTree::Plane tmpPlane;
1271 : : std::vector< EntityHandle > tmpChildren;
1272 : :
1273 : : public:
1274 : 53 : BSPTreePlaneIter( BSPTree* tool, const EntityHandle* path, int path_len )
1275 : 53 : : toolPtr( tool ), pathToRoot( path ), pathPos( path_len - 1 )
1276 : : {
1277 [ + - ]: 53 : operator++();
1278 : 53 : }
1279 : 53 : BSPTreePlaneIter() // initialize to 'end'
1280 : 53 : : toolPtr( 0 ), pathToRoot( 0 ), pathPos( -1 )
1281 : : {
1282 : 53 : }
1283 : :
1284 : 132 : const BSPTree::Plane* operator->() const
1285 : : {
1286 : 132 : return &tmpPlane;
1287 : : }
1288 : : bool operator==( const BSPTreePlaneIter& other ) const
1289 : : {
1290 : : return pathPos == other.pathPos;
1291 : : }
1292 : 107 : bool operator!=( const BSPTreePlaneIter& other ) const
1293 : : {
1294 : 107 : return pathPos != other.pathPos;
1295 : : }
1296 : : BSPTreePlaneIter& operator++();
1297 : : };
1298 : :
1299 : 107 : BSPTreePlaneIter& BSPTreePlaneIter::operator++()
1300 : : {
1301 [ + + ]: 107 : if( --pathPos < 0 ) return *this;
1302 : :
1303 : 54 : EntityHandle prev = pathToRoot[pathPos + 1];
1304 : 54 : EntityHandle curr = pathToRoot[pathPos];
1305 : :
1306 : 54 : ErrorCode rval = toolPtr->get_split_plane( curr, tmpPlane );
1307 [ - + ]: 54 : if( MB_SUCCESS != rval )
1308 : : {
1309 : 0 : assert( false );
1310 : : pathPos = 0;
1311 : : return *this;
1312 : : }
1313 : :
1314 : 54 : tmpChildren.clear();
1315 : 54 : rval = toolPtr->moab()->get_child_meshsets( curr, tmpChildren );
1316 [ + - ][ - + ]: 54 : if( MB_SUCCESS != rval || tmpChildren.size() != 2 )
[ - + ]
1317 : : {
1318 : 0 : assert( false );
1319 : : pathPos = 0;
1320 : : return *this;
1321 : : }
1322 : :
1323 [ + + ]: 54 : if( tmpChildren[1] == prev ) tmpPlane.flip();
1324 : 54 : return *this;
1325 : : }
1326 : :
1327 : 53 : bool BSPTreeIter::intersect_ray( const double ray_point[3], const double ray_vect[3], double& t_enter,
1328 : : double& t_exit ) const
1329 : : {
1330 : : // intersect with half-spaces defining tree
1331 [ + - ][ + - ]: 106 : BSPTreePlaneIter iter1( tool(), &mStack[0], mStack.size() ), end1;
[ + - ][ + - ]
1332 [ + - ][ + - ]: 53 : if( !ray_intersect_halfspaces( CartVect( ray_point ), CartVect( ray_vect ), iter1, end1, t_enter, t_exit ) )
[ + - ][ - + ]
1333 : 0 : return false;
1334 : :
1335 : : // itersect with box bounding entire tree
1336 : : double corners[8][3];
1337 [ + - ][ + - ]: 53 : ErrorCode rval = tool()->get_tree_box( mStack.front(), corners );
[ + - ]
1338 [ - + ]: 53 : if( MB_SUCCESS != rval )
1339 : : {
1340 : 0 : assert( false );
1341 : : return false;
1342 : : }
1343 : :
1344 [ + - ][ + - ]: 53 : BoxPlaneIter iter2( corners ), end2;
1345 : : double t2_enter, t2_exit;
1346 [ + - ][ + - ]: 53 : if( !ray_intersect_halfspaces( CartVect( ray_point ), CartVect( ray_vect ), iter2, end2, t2_enter, t2_exit ) )
[ + - ][ + + ]
1347 : 30 : return false;
1348 : :
1349 : : // if intersect both box and halfspaces, check that
1350 : : // two intersections overlap
1351 [ + + ]: 23 : if( t_enter < t2_enter ) t_enter = t2_enter;
1352 [ + + ]: 23 : if( t_exit > t2_exit ) t_exit = t2_exit;
1353 : 76 : return t_enter <= t_exit;
1354 : : }
1355 : :
1356 [ + - ][ + - ]: 4 : } // namespace moab
|
