Branch data Line data Source code
1 : : #include "moab/BVHTree.hpp"
2 : : #include "moab/Interface.hpp"
3 : : #include "moab/ElemEvaluator.hpp"
4 : : #include "moab/ReadUtilIface.hpp"
5 : : #include "moab/CpuTimer.hpp"
6 : :
7 : : namespace moab
8 : : {
9 : : const char* BVHTree::treeName = "BVHTree";
10 : :
11 : 1 : ErrorCode BVHTree::build_tree( const Range& entities, EntityHandle* tree_root_set, FileOptions* options )
12 : : {
13 : : ErrorCode rval;
14 [ + - ]: 1 : CpuTimer cp;
15 : :
16 [ - + ]: 1 : if( options )
17 : : {
18 [ # # ]: 0 : rval = parse_options( *options );
19 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
20 : :
21 [ # # ][ # # ]: 0 : if( !options->all_seen() ) return MB_FAILURE;
22 : : }
23 : :
24 : : // calculate bounding box of elements
25 [ + - ]: 1 : BoundBox box;
26 [ + - ][ + - ]: 1 : rval = box.update( *moab(), entities );
27 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
28 : :
29 : : // create tree root
30 : : EntityHandle tmp_root;
31 [ + - ]: 1 : if( !tree_root_set ) tree_root_set = &tmp_root;
32 [ + - ][ + - ]: 1 : rval = create_root( box.bMin.array(), box.bMax.array(), *tree_root_set );
[ + - ]
33 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
34 [ + - ]: 1 : rval = mbImpl->add_entities( *tree_root_set, entities );
35 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
36 : :
37 : : // a fully balanced tree will have 2*_entities.size()
38 : : // which is one doubling away..
39 [ + - ]: 2 : std::vector< Node > tree_nodes;
40 [ + - ][ + - ]: 1 : tree_nodes.reserve( entities.size() / maxPerLeaf );
41 [ + - ]: 2 : std::vector< HandleData > handle_data_vec;
42 [ + - ]: 1 : rval = construct_element_vec( handle_data_vec, entities, boundBox );
43 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
44 : :
45 : : #ifndef NDEBUG
46 [ + - ][ + - ]: 730 : for( std::vector< HandleData >::const_iterator i = handle_data_vec.begin(); i != handle_data_vec.end(); ++i )
[ + - ][ + + ]
47 : : {
48 [ + - ][ + - ]: 729 : if( !boundBox.intersects_box( i->myBox, 0 ) )
[ - + ]
49 : : {
50 [ # # ][ # # ]: 0 : std::cerr << "BB:" << boundBox << "EB:" << i->myBox << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
51 : 0 : return MB_FAILURE;
52 : : }
53 : : }
54 : : #endif
55 : : // We only build nonempty trees
56 [ + - ]: 1 : if( !handle_data_vec.empty() )
57 : : {
58 : : // initially all bits are set
59 [ + - ][ + - ]: 1 : tree_nodes.push_back( Node() );
60 [ + - ]: 1 : const int depth = local_build_tree( tree_nodes, handle_data_vec.begin(), handle_data_vec.end(), 0, boundBox );
61 : : #ifndef NDEBUG
62 [ + - ]: 1 : std::set< EntityHandle > entity_handles;
63 [ + - ][ + - ]: 378 : for( std::vector< Node >::iterator n = tree_nodes.begin(); n != tree_nodes.end(); ++n )
[ + + ]
64 : : {
65 [ + - ][ + - ]: 1106 : for( HandleDataVec::const_iterator j = n->entities.begin(); j != n->entities.end(); ++j )
[ + - ][ + - ]
[ + - ][ + + ]
66 : : {
67 [ + - ][ + - ]: 729 : entity_handles.insert( j->myHandle );
68 : : }
69 : : }
70 [ + - ][ - + ]: 1 : if( entity_handles.size() != entities.size() ) { std::cout << "Entity Handle Size Mismatch!" << std::endl; }
[ # # ][ # # ]
71 [ + - ][ + - ]: 730 : for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
[ + - ][ + - ]
[ + + ]
72 : : {
73 [ + - ][ + - ]: 729 : if( entity_handles.find( *i ) == entity_handles.end() )
[ + - ][ - + ]
74 [ # # ][ # # ]: 0 : std::cout << "Tree is missing an entity! " << std::endl;
75 : : }
76 : : #endif
77 [ + - ]: 1 : treeDepth = std::max( depth, treeDepth );
78 : : }
79 : :
80 : : // convert vector of Node's to entity sets and vector of TreeNode's
81 [ + - ]: 1 : rval = convert_tree( tree_nodes );
82 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
83 : :
84 [ + - ]: 1 : treeStats.reset();
85 [ + - ]: 1 : rval = treeStats.compute_stats( mbImpl, startSetHandle );
86 [ + - ]: 1 : treeStats.initTime = cp.time_elapsed();
87 : :
88 : 2 : return rval;
89 : : }
90 : :
91 : 1 : ErrorCode BVHTree::convert_tree( std::vector< Node >& tree_nodes )
92 : : {
93 : : // first construct the proper number of entity sets
94 : : ReadUtilIface* read_util;
95 [ + - ]: 1 : ErrorCode rval = mbImpl->query_interface( read_util );
96 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
97 : :
98 : : { // isolate potentially-large std::vector so it gets deleted earlier
99 [ + - ]: 1 : std::vector< unsigned int > tmp_flags( tree_nodes.size(), meshsetFlags );
100 [ + - ][ + - ]: 1 : rval = read_util->create_entity_sets( tree_nodes.size(), &tmp_flags[0], 0, startSetHandle );
101 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
102 [ + - ]: 1 : rval = mbImpl->release_interface( read_util );
103 [ - + ][ + - ]: 1 : if( MB_SUCCESS != rval ) return rval;
104 : : }
105 : :
106 : : // populate the sets and the TreeNode vector
107 : 1 : EntityHandle set_handle = startSetHandle;
108 : 1 : std::vector< Node >::iterator it;
109 [ + - ]: 1 : myTree.reserve( tree_nodes.size() );
110 [ + - ][ + - ]: 378 : for( it = tree_nodes.begin(); it != tree_nodes.end(); ++it, set_handle++ )
[ + + ]
111 : : {
112 [ + - ][ + + ]: 377 : if( it != tree_nodes.begin() && !it->entities.empty() )
[ + - ][ + + ]
[ + - ]
[ + + # # ]
113 : : {
114 [ + - ]: 189 : Range range;
115 [ + - ][ + - ]: 918 : for( HandleDataVec::iterator hit = it->entities.begin(); hit != it->entities.end(); ++hit )
[ + - ][ + - ]
[ + + ]
116 [ + - ][ + - ]: 729 : range.insert( hit->myHandle );
117 [ + - ]: 189 : rval = mbImpl->add_entities( set_handle, range );
118 [ - + ][ + - ]: 189 : if( MB_SUCCESS != rval ) return rval;
119 : : }
120 [ + - ][ + - ]: 377 : myTree.push_back( TreeNode( it->dim, it->child, it->Lmax, it->Rmin, it->box ) );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
121 : :
122 [ + - ][ + + ]: 377 : if( it->dim != 3 )
123 : : {
124 [ + - ][ + - ]: 188 : rval = mbImpl->add_child_meshset( set_handle, startSetHandle + it->child );
125 [ - + ]: 188 : if( MB_SUCCESS != rval ) return rval;
126 [ + - ][ + - ]: 188 : rval = mbImpl->add_child_meshset( set_handle, startSetHandle + it->child + 1 );
127 [ - + ]: 188 : if( MB_SUCCESS != rval ) return rval;
128 : : }
129 : : }
130 : :
131 : 1 : return MB_SUCCESS;
132 : : }
133 : :
134 : 1 : ErrorCode BVHTree::parse_options( FileOptions& opts )
135 : : {
136 [ + - ]: 1 : ErrorCode rval = parse_common_options( opts );
137 [ - + ]: 1 : if( MB_SUCCESS != rval ) return rval;
138 : :
139 : : // SPLITS_PER_DIR: number of candidate splits considered per direction; default = 3
140 : : int tmp_int;
141 [ + - ]: 1 : rval = opts.get_int_option( "SPLITS_PER_DIR", tmp_int );
142 [ - + ]: 1 : if( MB_SUCCESS == rval ) splitsPerDir = tmp_int;
143 : :
144 : 1 : return MB_SUCCESS;
145 : : }
146 : :
147 : 188 : void BVHTree::establish_buckets( HandleDataVec::const_iterator begin, HandleDataVec::const_iterator end,
148 : : const BoundBox& interval, std::vector< std::vector< Bucket > >& buckets ) const
149 : : {
150 : : // put each element into its bucket
151 [ + - ][ + - ]: 6686 : for( HandleDataVec::const_iterator i = begin; i != end; ++i )
[ + + ]
152 : : {
153 [ + - ]: 6498 : const BoundBox& box = i->myBox;
154 [ + + ]: 25992 : for( unsigned int dim = 0; dim < 3; ++dim )
155 : : {
156 [ + - ]: 19494 : const unsigned int index = Bucket::bucket_index( splitsPerDir, box, interval, dim );
157 [ + - ][ - + ]: 19494 : assert( index < buckets[dim].size() );
158 [ + - ][ + - ]: 19494 : Bucket& bucket = buckets[dim][index];
159 [ + + ]: 19494 : if( bucket.mySize > 0 )
160 [ + - ]: 18085 : bucket.boundingBox.update( box );
161 : : else
162 [ + - ]: 1409 : bucket.boundingBox = box;
163 : 19494 : bucket.mySize++;
164 : : }
165 : : }
166 : :
167 : : #ifndef NDEBUG
168 [ + - ]: 188 : BoundBox elt_union = begin->myBox;
169 [ + - ][ + - ]: 6686 : for( HandleDataVec::const_iterator i = begin; i != end; ++i )
[ + + ]
170 : : {
171 [ + - ]: 6498 : const BoundBox& box = i->myBox;
172 [ + - ]: 6498 : elt_union.update( box );
173 [ + + ]: 25992 : for( unsigned int dim = 0; dim < 3; ++dim )
174 : : {
175 [ + - ]: 19494 : const unsigned int index = Bucket::bucket_index( splitsPerDir, box, interval, dim );
176 [ + - ][ + - ]: 19494 : Bucket& bucket = buckets[dim][index];
177 [ + - ][ - + ]: 19494 : if( !bucket.boundingBox.intersects_box( box ) ) std::cerr << "Buckets not covering elements!" << std::endl;
[ # # ][ # # ]
178 : : }
179 : : }
180 [ + - ][ - + ]: 188 : if( !elt_union.intersects_box( interval ) )
181 : : {
182 [ # # ][ # # ]: 0 : std::cout << "element union: " << std::endl << elt_union;
[ # # ]
183 [ # # ][ # # ]: 0 : std::cout << "intervals: " << std::endl << interval;
[ # # ]
184 [ # # ][ # # ]: 0 : std::cout << "union of elts does not contain original box!" << std::endl;
185 : : }
186 [ + - ][ - + ]: 188 : if( !interval.intersects_box( elt_union ) )
187 : : {
188 [ # # ][ # # ]: 0 : std::cout << "original box does not contain union of elts" << std::endl;
189 [ # # ][ # # ]: 0 : std::cout << interval << std::endl << elt_union << std::endl;
[ # # ][ # # ]
190 : : }
191 [ + + ]: 752 : for( unsigned int d = 0; d < 3; ++d )
192 : : {
193 [ + - ]: 564 : std::vector< unsigned int > nonempty;
194 [ + - ]: 564 : const std::vector< Bucket >& buckets_ = buckets[d];
195 : 564 : unsigned int j = 0;
196 [ + - ][ + - ]: 2820 : for( std::vector< Bucket >::const_iterator i = buckets_.begin(); i != buckets_.end(); ++i, ++j )
[ + + ]
197 : : {
198 [ + - ][ + + ]: 2256 : if( i->mySize > 0 ) { nonempty.push_back( j ); }
[ + - ]
199 : : }
200 [ + - ][ + - ]: 1128 : BoundBox test_box = buckets_[nonempty.front()].boundingBox;
201 [ + + ]: 1973 : for( unsigned int i = 0; i < nonempty.size(); ++i )
202 [ + - ][ + - ]: 1409 : test_box.update( buckets_[nonempty[i]].boundingBox );
[ + - ]
203 : :
204 [ + - ][ - + ]: 564 : if( !test_box.intersects_box( interval ) )
205 [ # # ][ # # ]: 0 : std::cout << "union of buckets in dimension: " << d << "does not contain original box!" << std::endl;
[ # # ][ # # ]
206 [ + - ][ - + ]: 564 : if( !interval.intersects_box( test_box ) )
207 : : {
208 [ # # ]: 0 : std::cout << "original box does "
209 [ # # ]: 0 : << "not contain union of buckets"
210 [ # # ][ # # ]: 0 : << "in dimension: " << d << std::endl;
[ # # ]
211 [ # # ][ # # ]: 0 : std::cout << interval << std::endl << test_box << std::endl;
[ # # ][ # # ]
212 : : }
213 : 752 : }
214 : : #endif
215 : 188 : }
216 : :
217 : 188 : void BVHTree::initialize_splits( std::vector< std::vector< SplitData > >& splits,
218 : : const std::vector< std::vector< Bucket > >& buckets, const SplitData& data ) const
219 : : {
220 [ + + ]: 752 : for( unsigned int d = 0; d < 3; ++d )
221 : : {
222 [ + - ]: 564 : std::vector< SplitData >::iterator splits_begin = splits[d].begin();
223 [ + - ]: 564 : std::vector< SplitData >::iterator splits_end = splits[d].end();
224 [ + - ]: 564 : std::vector< Bucket >::const_iterator left_begin = buckets[d].begin();
225 [ + - ]: 564 : std::vector< Bucket >::const_iterator _end = buckets[d].end();
226 [ + - ][ + - ]: 564 : std::vector< Bucket >::const_iterator left_end = buckets[d].begin() + 1;
227 [ + - ][ + - ]: 2256 : for( std::vector< SplitData >::iterator s = splits_begin; s != splits_end; ++s, ++left_end )
[ + - ][ + + ]
228 : : {
229 [ + - ][ + - ]: 1692 : s->nl = set_interval( s->leftBox, left_begin, left_end );
[ + - ]
230 [ + - ][ + + ]: 1692 : if( s->nl == 0 )
231 : : {
232 [ + - ][ + - ]: 107 : s->leftBox = data.boundingBox;
233 [ + - ][ + - ]: 107 : s->leftBox.bMax[d] = s->leftBox.bMin[d];
[ + - ][ + - ]
234 : : }
235 [ + - ][ + - ]: 1692 : s->nr = set_interval( s->rightBox, left_end, _end );
[ + - ]
236 [ + - ][ + + ]: 1692 : if( s->nr == 0 )
237 : : {
238 [ + - ][ + - ]: 466 : s->rightBox = data.boundingBox;
239 [ + - ][ + - ]: 466 : s->rightBox.bMin[d] = s->rightBox.bMax[d];
[ + - ][ + - ]
240 : : }
241 [ + - ][ + - ]: 1692 : s->Lmax = s->leftBox.bMax[d];
[ + - ]
242 [ + - ][ + - ]: 1692 : s->Rmin = s->rightBox.bMin[d];
[ + - ]
243 [ + - ][ + - ]: 1692 : s->split = std::distance( splits_begin, s );
244 [ + - ]: 1692 : s->dim = d;
245 : : }
246 : : #ifndef NDEBUG
247 [ + - ][ + - ]: 2256 : for( std::vector< SplitData >::iterator s = splits_begin; s != splits_end; ++s )
[ + + ]
248 : : {
249 [ + - ]: 1692 : BoundBox test_box = s->leftBox;
250 [ + - ][ + - ]: 1692 : test_box.update( s->rightBox );
251 [ + - ][ - + ]: 1692 : if( !data.boundingBox.intersects_box( test_box ) )
252 : : {
253 [ # # ][ # # ]: 0 : std::cout << "nr: " << s->nr << std::endl;
[ # # ][ # # ]
254 [ # # ][ # # ]: 0 : std::cout << "Test box: " << std::endl << test_box;
[ # # ]
255 [ # # ][ # # ]: 0 : std::cout << "Left box: " << std::endl << s->leftBox;
[ # # ][ # # ]
256 [ # # ][ # # ]: 0 : std::cout << "Right box: " << std::endl << s->rightBox;
[ # # ][ # # ]
257 [ # # ][ # # ]: 0 : std::cout << "Interval: " << std::endl << data.boundingBox;
[ # # ]
258 [ # # ][ # # ]: 0 : std::cout << "Split boxes larger than bb" << std::endl;
259 : : }
260 [ + - ][ - + ]: 1692 : if( !test_box.intersects_box( data.boundingBox ) )
261 [ # # ][ # # ]: 0 : { std::cout << "bb larger than union of split boxes" << std::endl; }
262 : 1692 : }
263 : : #endif
264 : : }
265 : 188 : }
266 : :
267 : 92 : void BVHTree::median_order( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const
268 : : {
269 : 92 : int dim = data.dim;
270 [ + - ][ + - ]: 1526 : for( HandleDataVec::iterator i = begin; i != end; ++i )
[ + + ]
271 : : {
272 [ + - ][ + - ]: 1434 : i->myDim = 0.5 * ( i->myBox.bMin[dim], i->myBox.bMax[dim] );
[ + - ][ + - ]
[ + - ]
273 : : }
274 [ + - ]: 92 : std::sort( begin, end, BVHTree::HandleData_comparator() );
275 [ + - ]: 92 : const unsigned int total = std::distance( begin, end );
276 [ + - ]: 92 : HandleDataVec::iterator middle = begin + ( total / 2 );
277 [ + - ]: 92 : double middle_center = middle->myDim;
278 [ + - ][ + - ]: 92 : middle_center += ( ++middle )->myDim;
279 : 92 : middle_center *= 0.5;
280 : 92 : data.split = middle_center;
281 [ + - ]: 92 : data.nl = std::distance( begin, middle ) + 1;
282 : 92 : data.nr = total - data.nl;
283 [ + - ]: 92 : ++middle;
284 [ + - ][ + - ]: 92 : data.leftBox = begin->myBox;
285 [ + - ][ + - ]: 92 : data.rightBox = middle->myBox;
286 [ + - ][ + - ]: 964 : for( HandleDataVec::iterator i = begin; i != middle; ++i )
[ + + ]
287 : : {
288 [ + - ]: 872 : i->myDim = 0;
289 [ + - ][ + - ]: 872 : data.leftBox.update( i->myBox );
290 : : }
291 [ + - ][ + - ]: 654 : for( HandleDataVec::iterator i = middle; i != end; ++i )
[ + + ]
292 : : {
293 [ + - ]: 562 : i->myDim = 1;
294 [ + - ][ + - ]: 562 : data.rightBox.update( i->myBox );
295 : : }
296 [ + - ]: 92 : data.Rmin = data.rightBox.bMin[data.dim];
297 [ + - ]: 92 : data.Lmax = data.leftBox.bMax[data.dim];
298 : : #ifndef NDEBUG
299 : 92 : BoundBox test_box( data.rightBox );
300 [ + - ][ - + ]: 92 : if( !data.boundingBox.intersects_box( test_box ) )
301 : : {
302 [ # # ][ # # ]: 0 : std::cerr << "MEDIAN: BB Does not contain splits" << std::endl;
303 [ # # ][ # # ]: 0 : std::cerr << "test_box: " << test_box << std::endl;
[ # # ]
304 [ # # ][ # # ]: 0 : std::cerr << "data.boundingBox: " << data.boundingBox << std::endl;
[ # # ]
305 : 92 : }
306 : : #endif
307 : 92 : }
308 : :
309 : 188 : void BVHTree::find_split( HandleDataVec::iterator& begin, HandleDataVec::iterator& end, SplitData& data ) const
310 : : {
311 [ + - ][ + - ]: 188 : std::vector< std::vector< Bucket > > buckets( 3, std::vector< Bucket >( splitsPerDir + 1 ) );
312 [ + - ][ + - ]: 376 : std::vector< std::vector< SplitData > > splits( 3, std::vector< SplitData >( splitsPerDir, data ) );
313 : :
314 : 188 : const BoundBox interval = data.boundingBox;
315 [ + - ][ + - ]: 188 : establish_buckets( begin, end, interval, buckets );
[ + - ]
316 [ + - ]: 188 : initialize_splits( splits, buckets, data );
317 [ + - ]: 188 : choose_best_split( splits, data );
318 [ + + ][ - + ]: 188 : const bool use_median = ( 0 == data.nl ) || ( data.nr == 0 );
319 [ + + ]: 188 : if( !use_median )
320 [ + - ]: 96 : order_elements( begin, end, data );
321 : : else
322 [ + - ]: 92 : median_order( begin, end, data );
323 : :
324 : : #ifndef NDEBUG
325 : 188 : bool seen_one = false, issue = false;
326 : 188 : bool first_left = true, first_right = true;
327 : 188 : unsigned int count_left = 0, count_right = 0;
328 [ + - ][ + - ]: 376 : BoundBox left_box, right_box;
329 [ + - ][ + - ]: 6686 : for( HandleDataVec::iterator i = begin; i != end; ++i )
[ + + ]
330 : : {
331 [ + - ]: 6498 : int order = i->myDim;
332 [ + + ][ - + ]: 6498 : if( order != 0 && order != 1 )
333 : : {
334 [ # # ]: 0 : std::cerr << "Invalid order element !";
335 [ # # ][ # # ]: 0 : std::cerr << order << std::endl;
336 : 0 : std::exit( -1 );
337 : : }
338 [ + + ]: 6498 : if( order == 1 )
339 : : {
340 : 4315 : seen_one = 1;
341 : 4315 : count_right++;
342 [ + + ]: 4315 : if( first_right )
343 : : {
344 [ + - ][ + - ]: 188 : right_box = i->myBox;
345 : 188 : first_right = false;
346 : : }
347 : : else
348 : : {
349 [ + - ][ + - ]: 4127 : right_box.update( i->myBox );
350 : : }
351 [ + - ][ + - ]: 4315 : if( !right_box.intersects_box( i->myBox ) )
[ - + ]
352 : : {
353 [ # # ][ # # ]: 0 : if( !issue ) { std::cerr << "Bounding right box issue!" << std::endl; }
[ # # ]
354 : 0 : issue = true;
355 : : }
356 : : }
357 [ + + ]: 6498 : if( order == 0 )
358 : : {
359 : 2183 : count_left++;
360 [ + + ]: 2183 : if( first_left )
361 : : {
362 [ + - ][ + - ]: 188 : left_box = i->myBox;
363 : 188 : first_left = false;
364 : : }
365 : : else
366 : : {
367 [ + - ][ + - ]: 1995 : left_box.update( i->myBox );
368 : : }
369 [ + - ][ + - ]: 2183 : if( !data.leftBox.intersects_box( i->myBox ) )
[ - + ]
370 : : {
371 [ # # ][ # # ]: 0 : if( !issue ) { std::cerr << "Bounding left box issue!" << std::endl; }
[ # # ]
372 : 0 : issue = true;
373 : : }
374 [ - + ]: 2183 : if( seen_one )
375 : : {
376 [ # # ][ # # ]: 0 : std::cerr << "Invalid ordering!" << std::endl;
377 [ # # ][ # # ]: 0 : std::cout << ( i - 1 )->myDim << order << std::endl;
[ # # ][ # # ]
[ # # ]
378 : 0 : exit( -1 );
379 : : }
380 : : }
381 : : }
382 [ + - ][ - + ]: 188 : if( !left_box.intersects_box( data.leftBox ) ) std::cout << "left elts do not contain left box" << std::endl;
[ # # ][ # # ]
383 [ + - ][ - + ]: 188 : if( !data.leftBox.intersects_box( left_box ) ) std::cout << "left box does not contain left elts" << std::endl;
[ # # ][ # # ]
384 [ + - ][ - + ]: 188 : if( !right_box.intersects_box( data.rightBox ) ) std::cout << "right elts do not contain right box" << std::endl;
[ # # ][ # # ]
385 [ + - ][ - + ]: 188 : if( !data.rightBox.intersects_box( right_box ) ) std::cout << "right box do not contain right elts" << std::endl;
[ # # ][ # # ]
386 : :
387 [ + - ][ - + ]: 188 : if( count_left != data.nl || count_right != data.nr )
388 : : {
389 [ # # ][ # # ]: 0 : std::cerr << "counts are off!" << std::endl;
390 [ # # ][ # # ]: 0 : std::cerr << "total: " << std::distance( begin, end ) << std::endl;
[ # # ][ # # ]
391 [ # # ][ # # ]: 0 : std::cerr << "Dim: " << data.dim << std::endl;
[ # # ]
392 [ # # ][ # # ]: 0 : std::cerr << data.Lmax << " , " << data.Rmin << std::endl;
[ # # ][ # # ]
393 [ # # ][ # # ]: 0 : std::cerr << "Right box: " << std::endl << data.rightBox << "Left box: " << std::endl << data.leftBox;
[ # # ][ # # ]
[ # # ][ # # ]
394 [ # # ][ # # ]: 0 : std::cerr << "supposed to be: " << data.nl << " " << data.nr << std::endl;
[ # # ][ # # ]
[ # # ]
395 [ # # ][ # # ]: 0 : std::cerr << "accountant says: " << count_left << " " << count_right << std::endl;
[ # # ][ # # ]
[ # # ]
396 : 0 : std::exit( -1 );
397 : 188 : }
398 : : #endif
399 : 188 : }
400 : :
401 : 377 : int BVHTree::local_build_tree( std::vector< Node >& tree_nodes, HandleDataVec::iterator begin,
402 : : HandleDataVec::iterator end, const int index, const BoundBox& box, const int depth )
403 : : {
404 : : #ifndef NDEBUG
405 [ + - ][ + - ]: 7604 : for( HandleDataVec::const_iterator i = begin; i != end; ++i )
[ + - ][ + + ]
406 : : {
407 [ + - ][ + - ]: 7227 : if( !box.intersects_box( i->myBox, 0 ) )
[ - + ]
408 : : {
409 [ # # ][ # # ]: 0 : std::cerr << "depth: " << depth << std::endl;
[ # # ]
410 [ # # ][ # # ]: 0 : std::cerr << "BB:" << box << "EB:" << i->myBox << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
411 : 0 : std::exit( -1 );
412 : : }
413 : : }
414 : : #endif
415 : :
416 : 377 : const unsigned int total_num_elements = std::distance( begin, end );
417 : 377 : tree_nodes[index].box = box;
418 : :
419 : : // logic for splitting conditions
420 [ + + ][ + - ]: 377 : if( (int)total_num_elements > maxPerLeaf && depth < maxDepth )
421 : : {
422 [ + - ]: 188 : SplitData data;
423 [ + - ]: 188 : data.boundingBox = box;
424 [ + - ]: 188 : find_split( begin, end, data );
425 : : // assign data to node
426 [ + - ]: 188 : tree_nodes[index].Lmax = data.Lmax;
427 [ + - ]: 188 : tree_nodes[index].Rmin = data.Rmin;
428 [ + - ]: 188 : tree_nodes[index].dim = data.dim;
429 [ + - ]: 188 : tree_nodes[index].child = tree_nodes.size();
430 : : // insert left, right children;
431 [ + - ][ + - ]: 188 : tree_nodes.push_back( Node() );
432 [ + - ][ + - ]: 188 : tree_nodes.push_back( Node() );
433 : : const int left_depth =
434 [ + - ][ + - ]: 188 : local_build_tree( tree_nodes, begin, begin + data.nl, tree_nodes[index].child, data.leftBox, depth + 1 );
[ + - ]
435 : : const int right_depth =
436 [ + - ][ + - ]: 188 : local_build_tree( tree_nodes, begin + data.nl, end, tree_nodes[index].child + 1, data.rightBox, depth + 1 );
[ + - ]
437 [ + - ]: 188 : return std::max( left_depth, right_depth );
438 : : }
439 : :
440 : 189 : tree_nodes[index].dim = 3;
441 : 189 : std::copy( begin, end, std::back_inserter( tree_nodes[index].entities ) );
442 : 377 : return depth;
443 : : }
444 : :
445 : 0 : ErrorCode BVHTree::find_point( const std::vector< double >& point, const unsigned int& index, const double iter_tol,
446 : : const double inside_tol, std::pair< EntityHandle, CartVect >& result )
447 : : {
448 [ # # ]: 0 : if( index == 0 ) treeStats.numTraversals++;
449 [ # # ]: 0 : const TreeNode& node = myTree[index];
450 : 0 : treeStats.nodesVisited++;
451 [ # # ]: 0 : CartVect params;
452 : : int is_inside;
453 : 0 : ErrorCode rval = MB_SUCCESS;
454 [ # # ]: 0 : if( node.dim == 3 )
455 : : {
456 : 0 : treeStats.leavesVisited++;
457 [ # # ]: 0 : Range entities;
458 [ # # ]: 0 : rval = mbImpl->get_entities_by_handle( startSetHandle + index, entities );
459 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
460 : :
461 [ # # ][ # # ]: 0 : for( Range::iterator i = entities.begin(); i != entities.end(); ++i )
[ # # ][ # # ]
[ # # ]
462 : : {
463 : 0 : treeStats.traversalLeafObjectTests++;
464 [ # # ][ # # ]: 0 : myEval->set_ent_handle( *i );
465 [ # # ][ # # ]: 0 : myEval->reverse_eval( &point[0], iter_tol, inside_tol, params.array(), &is_inside );
[ # # ]
466 [ # # ]: 0 : if( is_inside )
467 : : {
468 [ # # ]: 0 : result.first = *i;
469 : 0 : result.second = params;
470 : 0 : return MB_SUCCESS;
471 : : }
472 : : }
473 : 0 : result.first = 0;
474 : 0 : return MB_SUCCESS;
475 : : }
476 : : // the extra tol here considers the case where
477 : : // 0 < Rmin - Lmax < 2tol
478 [ # # ]: 0 : std::vector< EntityHandle > children;
479 [ # # ]: 0 : rval = mbImpl->get_child_meshsets( startSetHandle + index, children );
480 [ # # ][ # # ]: 0 : if( MB_SUCCESS != rval || children.size() != 2 ) return rval;
[ # # ]
481 : :
482 [ # # ]: 0 : if( ( node.Lmax + iter_tol ) < ( node.Rmin - iter_tol ) )
483 : : {
484 [ # # ][ # # ]: 0 : if( point[node.dim] <= ( node.Lmax + iter_tol ) )
485 [ # # ][ # # ]: 0 : return find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
486 [ # # ][ # # ]: 0 : else if( point[node.dim] >= ( node.Rmin - iter_tol ) )
487 [ # # ][ # # ]: 0 : return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
488 [ # # ][ # # ]: 0 : result = std::make_pair( 0, CartVect( &point[0] ) ); // point lies in empty space.
[ # # ][ # # ]
489 : 0 : return MB_SUCCESS;
490 : : }
491 : :
492 : : // Boxes overlap
493 : : // left of Rmin, you must be on the left
494 : : // we can't be sure about the boundaries since the boxes overlap
495 : : // this was a typo in the paper which caused pain.
496 [ # # ][ # # ]: 0 : if( point[node.dim] < ( node.Rmin - iter_tol ) )
497 [ # # ][ # # ]: 0 : return find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
498 : : // if you are on the right Lmax, you must be on the right
499 [ # # ][ # # ]: 0 : else if( point[node.dim] > ( node.Lmax + iter_tol ) )
500 [ # # ][ # # ]: 0 : return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result );
501 : :
502 : : /* pg5 of paper
503 : : * However, instead of always traversing either subtree
504 : : * first (e.g. left always before right), we first traverse
505 : : * the subtree whose bounding plane has the larger distance to the
506 : : * sought point. This results in less overall traversal, and the correct
507 : : * cell is identified more quickly.
508 : : */
509 : : // So far all testing confirms that this 'heuristic' is
510 : : // significantly slower.
511 : : // I conjecture this is because it gets improperly
512 : : // branch predicted..
513 : : // bool dir = (point[ node.dim] - node.Rmin) <=
514 : : // (node.Lmax - point[ node.dim]);
515 [ # # ][ # # ]: 0 : find_point( point, children[0] - startSetHandle, iter_tol, inside_tol, result );
516 [ # # ][ # # ]: 0 : if( result.first == 0 ) { return find_point( point, children[1] - startSetHandle, iter_tol, inside_tol, result ); }
[ # # ]
517 : 0 : return MB_SUCCESS;
518 : : }
519 : :
520 : 0 : EntityHandle BVHTree::bruteforce_find( const double* point, const double iter_tol, const double inside_tol )
521 : : {
522 : 0 : treeStats.numTraversals++;
523 [ # # ]: 0 : CartVect params;
524 [ # # ]: 0 : for( unsigned int i = 0; i < myTree.size(); i++ )
525 : : {
526 [ # # ][ # # ]: 0 : if( myTree[i].dim != 3 || !myTree[i].box.contains_point( point, iter_tol ) ) continue;
[ # # ][ # # ]
[ # # ][ # # ]
527 [ # # ]: 0 : if( myEval )
528 : : {
529 : 0 : EntityHandle entity = 0;
530 : 0 : treeStats.leavesVisited++;
531 : : ErrorCode rval = myEval->find_containing_entity( startSetHandle + i, point, iter_tol, inside_tol, entity,
532 [ # # ][ # # ]: 0 : params.array(), &treeStats.traversalLeafObjectTests );
533 [ # # ]: 0 : if( entity )
534 : 0 : return entity;
535 [ # # ]: 0 : else if( MB_SUCCESS != rval )
536 : 0 : return 0;
537 : : }
538 : : else
539 : 0 : return startSetHandle + i;
540 : : }
541 : 0 : return 0;
542 : : }
543 : :
544 : 32854 : ErrorCode BVHTree::get_bounding_box( BoundBox& box, EntityHandle* tree_node ) const
545 : : {
546 [ + + ][ + + ]: 32854 : if( !tree_node || *tree_node == startSetHandle )
547 : : {
548 : 2001 : box = boundBox;
549 : 2001 : return MB_SUCCESS;
550 : : }
551 [ + - ][ + - ]: 61706 : else if( ( tree_node && !startSetHandle ) || *tree_node < startSetHandle ||
[ + - - + ]
[ - + ]
552 : 30853 : *tree_node - startSetHandle > myTree.size() )
553 : 0 : return MB_FAILURE;
554 : :
555 : 30853 : box = myTree[*tree_node - startSetHandle].box;
556 : 30853 : return MB_SUCCESS;
557 : : }
558 : :
559 : 2000 : ErrorCode BVHTree::point_search( const double* point, EntityHandle& leaf_out, const double iter_tol,
560 : : const double inside_tol, bool* multiple_leaves, EntityHandle* start_node,
561 : : CartVect* params )
562 : : {
563 : 2000 : treeStats.numTraversals++;
564 : :
565 [ - + ]: 2000 : EntityHandle this_set = ( start_node ? *start_node : startSetHandle );
566 : : // convoluted check because the root is different from startSetHandle
567 [ + - ][ + - ]: 2000 : if( this_set != myRoot && ( this_set < startSetHandle || this_set >= startSetHandle + myTree.size() ) )
[ - + ][ - + ]
568 : 0 : return MB_FAILURE;
569 [ - + ]: 2000 : else if( this_set == myRoot )
570 : 0 : this_set = startSetHandle;
571 : :
572 [ + - ]: 2000 : std::vector< EntityHandle > candidates,
573 [ + - ]: 4000 : result_list; // list of subtrees to traverse, and results
574 [ + - ]: 2000 : candidates.push_back( this_set - startSetHandle );
575 : :
576 [ + - ]: 4000 : BoundBox box;
577 [ + + ]: 33853 : while( !candidates.empty() )
578 : : {
579 [ + - ]: 32853 : EntityHandle ind = candidates.back();
580 : 32853 : treeStats.nodesVisited++;
581 [ + - ][ + + ]: 32853 : if( myTree[ind].dim == 3 ) treeStats.leavesVisited++;
582 : 32853 : this_set = startSetHandle + ind;
583 [ + - ]: 32853 : candidates.pop_back();
584 : :
585 : : // test box of this node
586 [ + - ]: 32853 : ErrorCode rval = get_bounding_box( box, &this_set );
587 [ - + ]: 32853 : if( MB_SUCCESS != rval ) return rval;
588 [ + - ][ + + ]: 32853 : if( !box.contains_point( point, iter_tol ) ) continue;
589 : :
590 : : // else if not a leaf, test children & put on list
591 [ + - ][ + + ]: 20827 : else if( myTree[ind].dim != 3 )
592 : : {
593 [ + - ][ + - ]: 18354 : candidates.push_back( myTree[ind].child );
594 [ + - ][ + - ]: 18354 : candidates.push_back( myTree[ind].child + 1 );
595 : 18354 : continue;
596 : : }
597 [ + - ][ + - ]: 2473 : else if( myTree[ind].dim == 3 && myEval && params )
[ + + ][ + - ]
[ + + ]
598 : : {
599 : : rval = myEval->find_containing_entity( startSetHandle + ind, point, iter_tol, inside_tol, leaf_out,
600 [ + - ][ + - ]: 1142 : params->array(), &treeStats.traversalLeafObjectTests );
601 [ + + ][ - + ]: 1142 : if( leaf_out || MB_SUCCESS != rval ) return rval;
602 : : }
603 : : else
604 : : {
605 : : // leaf node within distance; return in list
606 [ + - ]: 1331 : result_list.push_back( this_set );
607 : : }
608 : : }
609 : :
610 [ + - ][ + - ]: 1000 : if( !result_list.empty() ) leaf_out = result_list[0];
611 [ - + ][ # # ]: 1000 : if( multiple_leaves && result_list.size() > 1 ) *multiple_leaves = true;
[ - + ]
612 : 3000 : return MB_SUCCESS;
613 : : }
614 : :
615 : 0 : ErrorCode BVHTree::distance_search( const double from_point[3], const double distance,
616 : : std::vector< EntityHandle >& result_list, const double iter_tol,
617 : : const double inside_tol, std::vector< double >* result_dists,
618 : : std::vector< CartVect >* result_params, EntityHandle* tree_root )
619 : : {
620 : : // non-NULL root should be in tree
621 : : // convoluted check because the root is different from startSetHandle
622 [ # # ]: 0 : EntityHandle this_set = ( tree_root ? *tree_root : startSetHandle );
623 [ # # ][ # # ]: 0 : if( this_set != myRoot && ( this_set < startSetHandle || this_set >= startSetHandle + myTree.size() ) )
[ # # ][ # # ]
624 : 0 : return MB_FAILURE;
625 [ # # ]: 0 : else if( this_set == myRoot )
626 : 0 : this_set = startSetHandle;
627 : :
628 : 0 : treeStats.numTraversals++;
629 : :
630 : 0 : const double dist_sqr = distance * distance;
631 [ # # ]: 0 : const CartVect from( from_point );
632 [ # # ]: 0 : std::vector< EntityHandle > candidates; // list of subtrees to traverse
633 : : // pre-allocate space for default max tree depth
634 [ # # ]: 0 : candidates.reserve( maxDepth );
635 : :
636 : : // misc temporary values
637 : : ErrorCode rval;
638 [ # # ]: 0 : BoundBox box;
639 : :
640 [ # # ]: 0 : candidates.push_back( this_set - startSetHandle );
641 : :
642 [ # # ]: 0 : while( !candidates.empty() )
643 : : {
644 : :
645 [ # # ]: 0 : EntityHandle ind = candidates.back();
646 : 0 : this_set = startSetHandle + ind;
647 [ # # ]: 0 : candidates.pop_back();
648 : 0 : treeStats.nodesVisited++;
649 [ # # ][ # # ]: 0 : if( myTree[ind].dim == 3 ) treeStats.leavesVisited++;
650 : :
651 : : // test box of this node
652 [ # # ]: 0 : rval = get_bounding_box( box, &this_set );
653 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
654 [ # # ]: 0 : double d_sqr = box.distance_squared( from_point );
655 : :
656 : : // if greater than cutoff, continue
657 [ # # ]: 0 : if( d_sqr > dist_sqr ) continue;
658 : :
659 : : // else if not a leaf, test children & put on list
660 [ # # ][ # # ]: 0 : else if( myTree[ind].dim != 3 )
661 : : {
662 [ # # ][ # # ]: 0 : candidates.push_back( myTree[ind].child );
663 [ # # ][ # # ]: 0 : candidates.push_back( myTree[ind].child + 1 );
664 : 0 : continue;
665 : : }
666 : :
667 [ # # ][ # # ]: 0 : if( myEval && result_params )
668 : : {
669 : : EntityHandle ent;
670 [ # # ]: 0 : CartVect params;
671 : : rval = myEval->find_containing_entity( startSetHandle + ind, from_point, iter_tol, inside_tol, ent,
672 [ # # ][ # # ]: 0 : params.array(), &treeStats.traversalLeafObjectTests );
673 [ # # ]: 0 : if( MB_SUCCESS != rval )
674 : 0 : return rval;
675 [ # # ]: 0 : else if( ent )
676 : : {
677 [ # # ]: 0 : result_list.push_back( ent );
678 [ # # ]: 0 : result_params->push_back( params );
679 [ # # ][ # # ]: 0 : if( result_dists ) result_dists->push_back( 0.0 );
680 : 0 : }
681 : : }
682 : : else
683 : : {
684 : : // leaf node within distance; return in list
685 [ # # ]: 0 : result_list.push_back( this_set );
686 [ # # ][ # # ]: 0 : if( result_dists ) result_dists->push_back( sqrt( d_sqr ) );
687 : : }
688 : : }
689 : :
690 : 0 : return MB_SUCCESS;
691 : : }
692 : :
693 : 0 : ErrorCode BVHTree::print_nodes( std::vector< Node >& nodes )
694 : : {
695 : : int i;
696 : 0 : std::vector< Node >::iterator it;
697 [ # # ][ # # ]: 0 : for( it = nodes.begin(), i = 0; it != nodes.end(); ++it, i++ )
[ # # ]
698 : : {
699 [ # # ][ # # ]: 0 : std::cout << "Node " << i << ": dim = " << it->dim << ", child = " << it->child << ", Lmax/Rmin = " << it->Lmax
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
700 [ # # ][ # # ]: 0 : << "/" << it->Rmin << ", box = " << it->box << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
701 : : }
702 : 0 : return MB_SUCCESS;
703 : : }
704 : :
705 : 0 : ErrorCode BVHTree::print()
706 : : {
707 : : int i;
708 : 0 : std::vector< TreeNode >::iterator it;
709 [ # # ][ # # ]: 0 : for( it = myTree.begin(), i = 0; it != myTree.end(); ++it, i++ )
[ # # ]
710 : : {
711 [ # # ][ # # ]: 0 : std::cout << "Node " << i << ": dim = " << it->dim << ", child = " << it->child << ", Lmax/Rmin = " << it->Lmax
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
712 [ # # ][ # # ]: 0 : << "/" << it->Rmin << ", box = " << it->box << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
713 : : }
714 : 0 : return MB_SUCCESS;
715 : : }
716 : :
717 [ + - ][ + - ]: 4 : } // namespace moab
|
