1 #include "moab/Core.hpp"
2 #include "TestUtil.hpp"
3 #include "moab/BSPTree.hpp"
4 #include "moab/CartVect.hpp"
5 #include "moab/BSPTreePoly.hpp"
6 #include "moab/Range.hpp"
7 #include <algorithm>
8 
9 using namespace moab;
10 
11 void test_set_plane();
12 void test_iterator();
13 void test_box_iterator();
14 void test_tree_create();
15 void test_box_tree_create();
16 void test_leaf_containing_point_bounded_tree();
17 void test_leaf_containing_point_unbounded_tree();
18 void test_merge_leaf();
19 void test_box_iter_neighbors();
20 void test_leaf_sibling();
21 void test_leaf_volume( bool box );
test_leaf_volume_box()22 void test_leaf_volume_box() { test_leaf_volume(true); }
test_leaf_volume_gen()23 void test_leaf_volume_gen() { test_leaf_volume(false); }
24 void test_leaf_splits_intersects();
25 void test_leaf_intersects_ray_common( bool box );
test_box_leaf_intersects_ray()26 void test_box_leaf_intersects_ray() { test_leaf_intersects_ray_common(true); }
test_gen_leaf_intersects_ray()27 void test_gen_leaf_intersects_ray() { test_leaf_intersects_ray_common(false); }
28 void test_leaf_polyhedron();
29 
main()30 int main()
31 {
32   int failures = 0;
33 
34   failures += RUN_TEST( test_set_plane );
35   failures += RUN_TEST( test_iterator );
36   failures += RUN_TEST( test_box_iterator );
37   failures += RUN_TEST( test_tree_create );
38   failures += RUN_TEST( test_box_tree_create );
39   failures += RUN_TEST( test_leaf_containing_point_bounded_tree );
40   failures += RUN_TEST( test_leaf_containing_point_unbounded_tree );
41   failures += RUN_TEST( test_merge_leaf );
42   failures += RUN_TEST( test_box_iter_neighbors );
43   failures += RUN_TEST( test_leaf_sibling );
44   failures += RUN_TEST( test_leaf_volume_box );
45   failures += RUN_TEST( test_leaf_volume_gen );
46   failures += RUN_TEST( test_leaf_splits_intersects );
47   failures += RUN_TEST( test_box_leaf_intersects_ray );
48   failures += RUN_TEST( test_gen_leaf_intersects_ray );
49   failures += RUN_TEST( test_leaf_polyhedron );
50 
51   return failures;
52 }
53 
54 
55 // Make CHECK_EQUAL macro work for planes
check_equal(const BSPTree::Plane & p1,const BSPTree::Plane & p2,const char * exp1,const char * exp2,int line,const char * file)56 void check_equal( const BSPTree::Plane& p1,
57                   const BSPTree::Plane& p2,
58                   const char* exp1,
59                   const char* exp2,
60                   int line,
61                   const char* file )
62 {
63   if (fabs(p1.norm[0]-p2.norm[0]) > 1e-6 ||
64       fabs(p1.norm[1]-p2.norm[1]) > 1e-6 ||
65       fabs(p1.norm[2]-p2.norm[2]) > 1e-6 ||
66       fabs(p1.coeff  -p2.coeff  ) > 1e-6) {
67     printf( "Equality Test Failed: %s == %s\n", exp1, exp2 );
68     printf( "  at line %d of '%s'\n", line, file );
69     printf( "  Expected: %f x + %f y + %f z + %f = 0\n",
70             p1.norm[0], p1.norm[1], p1.norm[2], p1.coeff );
71     printf( "  Actual  : %f x + %f y + %f z + %f = 0\n",
72             p2.norm[0], p2.norm[1], p2.norm[2], p2.coeff );
73     printf( "\n" );
74     FLAG_ERROR;
75   }
76 }
77 
test_set_plane()78 void test_set_plane()
79 {
80   BSPTree::Plane p;
81   const double points[3][3] = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
82   p.set( points[0], points[1], points[2] );
83   CHECK_REAL_EQUAL( 0.0, p.distance( points[0] ), 1e-10 );
84   CHECK_REAL_EQUAL( 0.0, p.distance( points[1] ), 1e-10 );
85   CHECK_REAL_EQUAL( 0.0, p.distance( points[2] ), 1e-10 );
86 }
87 
test_iterator()88 void test_iterator()
89 {
90   Core moab;
91   BSPTree tool( &moab );
92   ErrorCode rval;
93   EntityHandle root;
94   BSPTreeIter iter;
95 
96     // create a depth-1 tree
97   rval = tool.create_tree( root );
98   CHECK_ERR(rval);
99   rval = tool.get_tree_iterator( root, iter );
100   CHECK_ERR(rval);
101 
102     // check initial state of iterator
103   CHECK_EQUAL( &tool, iter.tool() );
104   CHECK_EQUAL(  root, iter.handle() );
105   CHECK_EQUAL(  1u,   iter.depth() );
106 
107     // should fail if at root
108   BSPTree::Plane p;
109   rval = iter.get_parent_split_plane( p );
110   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
111 
112     // check that step past end returns correct value
113   rval = iter.step();
114   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
115 
116     // reset iterator
117   rval = tool.get_tree_iterator( root, iter );
118   CHECK_ERR(rval);
119 
120     // check that step past start returns correct value
121   rval = iter.back();
122   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
123 
124     // reset iterator
125   rval = tool.get_tree_iterator( root, iter );
126   CHECK_ERR(rval);
127 
128     // insert a single split plane
129   rval = tool.split_leaf( iter, BSPTree::Plane(2,0,0,0) );
130   CHECK_ERR(rval);
131 
132     // check initial location
133   CHECK_EQUAL( 2u, iter.depth() );
134   CHECK( iter.handle() != root );
135 
136     // create new iterators at left and right ends
137   BSPTreeIter left, right;
138   rval = tool.get_tree_iterator( root, left );
139   CHECK_ERR(rval);
140   rval = tool.get_tree_end_iterator( root, right );
141   CHECK_ERR(rval);
142 
143     // compare post-split iterator to left
144   CHECK_EQUAL( left.depth(), iter.depth() );
145   CHECK_EQUAL( left.handle(), iter.handle() );
146 
147     // step to other leaf
148   rval = iter.step();
149   CHECK_ERR(rval);
150 
151     // check location
152   CHECK_EQUAL( 2u, iter.depth() );
153   CHECK( iter.handle() != root );
154 
155     // compare stepped iterator to right
156   CHECK_EQUAL( right.depth(), iter.depth() );
157   CHECK_EQUAL( right.handle(), iter.handle() );
158 
159     // step to back to first leaf
160   rval = iter.back();
161   CHECK_ERR(rval);
162 
163     // compare to left
164   CHECK_EQUAL( left.depth(), iter.depth() );
165   CHECK_EQUAL( left.handle(), iter.handle() );
166 
167     // check plane
168     // should have unit normal
169   left.get_parent_split_plane( p );
170   CHECK_EQUAL( BSPTree::Plane(1,0,0,0), p );
171   p.norm[0] = 11;
172   right.get_parent_split_plane( p );
173   CHECK_EQUAL( BSPTree::Plane(1,0,0,0), p );
174 
175     // check step past ends
176   rval = left.back();
177   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
178   rval = right.step();
179   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
180 }
181 
compare_hexes(const double expected[8][3],const double actual[8][3],double epsilon)182 bool compare_hexes( const double expected[8][3],
183                     const double actual[8][3],
184                     double epsilon )
185 {
186     // for each of three possible rotations
187   const int rotation_maps[3][8] = { { 0, 1, 2, 3, 4, 5, 6, 7 },
188                                     { 3, 2, 6, 7, 0, 1, 5, 4 },
189                                     { 7, 3, 2, 6, 4, 0, 1, 5 } };
190   for (int i = 0; i < 3; ++i) {
191       // compare for rotationts about axis from face 4 to 5
192     for (int j = 0; j < 4; ++j) {
193       bool match = true;
194         // for each face vertex
195       for (int k = 0; k < 4 && match; ++k) {
196           // for each coordinate
197         for (int d = 0; d < 3; ++d) {
198           if (fabs(expected[(j+k)%4  ][d] - actual[rotation_maps[i][k  ]][d]) > epsilon
199            || fabs(expected[(j+k)%4+4][d] - actual[rotation_maps[i][k+4]][d]) > epsilon) {
200             match = false;
201             break;
202           }
203         }
204       }
205 
206       if (match)
207         return true;
208     }
209   }
210 
211   printf("Hex vertex copmarison failed.\n");
212   printf("           Exected                         Actual\n");
213   for (int v = 0; v < 8; ++v) {
214     printf("% 9f % 9f % 9f   % 9f % 9f % 9f\n",
215       expected[v][0], expected[v][1], expected[v][2],
216       actual[v][0], actual[v][1], actual[v][2]);
217   }
218 
219 
220   return false;
221 }
222 
aabox_corners(const double min[3],const double max[3],double corners[8][3])223 static void aabox_corners( const double min[3],
224                            const double max[3],
225                            double corners[8][3] )
226 {
227   const double expt[8][3] = { { min[0], min[1], min[2] },
228                               { max[0], min[1], min[2] },
229                               { max[0], max[1], min[2] },
230                               { min[0], max[1], min[2] },
231                               { min[0], min[1], max[2] },
232                               { max[0], min[1], max[2] },
233                               { max[0], max[1], max[2] },
234                               { min[0], max[1], max[2] } };
235   memcpy( corners, expt, sizeof(expt) );
236 }
237 
aabox_corners(double min_x,double min_y,double min_z,double max_x,double max_y,double max_z,double corners[8][3])238 static void aabox_corners( double min_x, double min_y, double min_z,
239                            double max_x, double max_y, double max_z,
240                            double corners[8][3] )
241 {
242   double min[] = { min_x, min_y, min_z };
243   double max[] = { max_x, max_y, max_z };
244   aabox_corners( min, max, corners );
245 }
246 
247 
test_box_iterator()248 void test_box_iterator()
249 {
250   Core moab;
251   BSPTree tool( &moab );
252   ErrorCode rval;
253   EntityHandle root;
254   BSPTreeBoxIter iter;
255   const double min[3] = { -1, -2, -3 };
256   const double max[3] = {  3,  2,  1 };
257 
258     // create a depth-1 tree
259   rval = tool.create_tree( min, max, root );
260   CHECK_ERR(rval);
261   rval = tool.get_tree_iterator( root, iter );
262   CHECK_ERR(rval);
263 
264     // check initial state of iterator
265   CHECK_EQUAL( &tool, iter.tool() );
266   CHECK_EQUAL(  root, iter.handle() );
267   CHECK_EQUAL(  1u,   iter.depth() );
268 
269     // check initial corner values
270   double corners[8][3], expt[8][3];
271   aabox_corners( min, max, expt );
272   rval = iter.get_box_corners( corners );
273   CHECK_ERR( rval );
274   CHECK( compare_hexes( expt, corners, 1e-6 ) );
275 
276     // should fail if at root
277   BSPTree::Plane p;
278   rval = iter.get_parent_split_plane( p );
279   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
280 
281     // check that step past end returns correct value
282   rval = iter.step();
283   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
284 
285     // reset iterator
286   rval = tool.get_tree_iterator( root, iter );
287   CHECK_ERR(rval);
288 
289     // check that step past start returns correct value
290   rval = iter.back();
291   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
292 
293     // reset iterator
294   rval = tool.get_tree_iterator( root, iter );
295   CHECK_ERR(rval);
296 
297     // insert a single split plane
298   rval = tool.split_leaf( iter, BSPTree::Plane(2,0,0,0) );
299   CHECK_ERR(rval);
300 
301     // check initial location
302   CHECK_EQUAL( 2u, iter.depth() );
303   CHECK( iter.handle() != root );
304 
305     // create new iterators at left and right ends
306   BSPTreeIter left, right;
307   rval = tool.get_tree_iterator( root, left );
308   CHECK_ERR(rval);
309   rval = tool.get_tree_end_iterator( root, right );
310   CHECK_ERR(rval);
311 
312     // compare post-split iterator to left
313   CHECK_EQUAL( left.depth(), iter.depth() );
314   CHECK_EQUAL( left.handle(), iter.handle() );
315 
316     // check box
317   aabox_corners( min, max, expt );
318   for (int i = 0; i < 8; ++i)
319     if (expt[i][0] > 0)
320       expt[i][0] = 0;
321   rval = iter.get_box_corners( corners );
322   CHECK_ERR( rval );
323   CHECK( compare_hexes( expt, corners, 1e-6 ) );
324 
325     // step to other leaf
326   rval = iter.step();
327   CHECK_ERR(rval);
328 
329     // check location
330   CHECK_EQUAL( 2u, iter.depth() );
331   CHECK( iter.handle() != root );
332 
333     // compare stepped iterator to right
334   CHECK_EQUAL( right.depth(), iter.depth() );
335   CHECK_EQUAL( right.handle(), iter.handle() );
336 
337     // check box
338   aabox_corners( min, max, expt );
339   for (int i = 0; i < 8; ++i)
340     if (expt[i][0] < 0)
341       expt[i][0] = 0;
342   rval = iter.get_box_corners( corners );
343   CHECK_ERR( rval );
344   CHECK( compare_hexes( expt, corners, 1e-6 ) );
345 
346     // step to back to first leaf
347   rval = iter.back();
348   CHECK_ERR(rval);
349 
350     // compare to left
351   CHECK_EQUAL( left.depth(), iter.depth() );
352   CHECK_EQUAL( left.handle(), iter.handle() );
353 
354     // check box
355   aabox_corners( min, max, expt );
356   for (int i = 0; i < 8; ++i)
357     if (expt[i][0] > 0)
358       expt[i][0] = 0;
359   rval = iter.get_box_corners( corners );
360   CHECK_ERR( rval );
361   CHECK( compare_hexes( expt, corners, 1e-6 ) );
362 
363     // check plane
364     // should have unit normal
365   left.get_parent_split_plane( p );
366   CHECK_EQUAL( BSPTree::Plane(1,0,0,0), p );
367   p.norm[0] = 11;
368   right.get_parent_split_plane( p );
369   CHECK_EQUAL( BSPTree::Plane(1,0,0,0), p );
370 
371     // check step past ends
372   rval = left.back();
373   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
374   rval = right.step();
375   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
376 }
377 
test_tree_create()378 void test_tree_create()
379 {
380   Core moab;
381   BSPTree tool( &moab );
382   ErrorCode rval;
383   EntityHandle root;
384   BSPTreeIter iter;
385   BSPTree::Plane p;
386 
387     // create a depth-1 tree
388   rval = tool.create_tree( root );
389   CHECK_ERR(rval);
390   rval = tool.get_tree_iterator( root, iter );
391   CHECK_ERR(rval);
392 
393     // check initial state of iterator
394   CHECK_EQUAL( &tool, iter.tool() );
395   CHECK_EQUAL(  root, iter.handle() );
396   CHECK_EQUAL(  1u,   iter.depth() );
397 
398     // split with z=0
399   rval = tool.split_leaf( iter, BSPTree::Plane(0,0,0.5,0) );
400   CHECK_ERR(rval);
401   CHECK_EQUAL( 2u, iter.depth() );
402 
403     // check that parent split plane is correct
404   rval = iter.get_parent_split_plane( p );
405   CHECK_ERR(rval);
406   CHECK_EQUAL( BSPTree::Plane(0,0,1,0), p );
407 
408     // split lower leaf with diagonal plane
409   rval = tool.split_leaf( iter, BSPTree::Plane(1,1,0,0) );
410   CHECK_ERR(rval);
411   CHECK_EQUAL( 3u, iter.depth() );
412 
413   const double r = 0.5*sqrt(2.0);
414 
415     // check that parent split plane is correct
416   rval = iter.get_parent_split_plane( p );
417   CHECK_ERR(rval);
418   CHECK_EQUAL( BSPTree::Plane(r,r,0,0), p );
419 
420     // step to upper leaf
421   rval = iter.step();
422   CHECK_ERR(rval);
423 
424     // split upper leaf with diagonal plane
425   rval = tool.split_leaf( iter, BSPTree::Plane(-1,1,0,0) );
426   CHECK_ERR(rval);
427   CHECK_EQUAL( 4u, iter.depth() );
428 
429     // check that parent split plane is correct
430   rval = iter.get_parent_split_plane( p );
431   CHECK_ERR(rval);
432   CHECK_EQUAL( BSPTree::Plane(-r,r,0,0), p );
433 
434     // iterate over four leaves
435 
436     // first leaf
437   rval = tool.get_tree_iterator( root, iter );
438   CHECK_ERR(rval);
439   CHECK_EQUAL( 3u, iter.depth() );
440   rval = iter.get_parent_split_plane( p );
441   CHECK_ERR(rval);
442   CHECK_EQUAL( BSPTree::Plane(r,r,0,0), p );
443 
444     // second leaf
445   rval = iter.step();
446   CHECK_ERR(rval);
447   CHECK_EQUAL( 4u, iter.depth() );
448   rval = iter.get_parent_split_plane( p );
449   CHECK_ERR(rval);
450   CHECK_EQUAL( BSPTree::Plane(-r,r,0,0), p );
451 
452     // third leaf
453   rval = iter.step();
454   CHECK_ERR(rval);
455   CHECK_EQUAL( 4u, iter.depth() );
456   rval = iter.get_parent_split_plane( p );
457   CHECK_ERR(rval);
458   CHECK_EQUAL( BSPTree::Plane(-r,r,0,0), p );
459 
460     // fourth leaf
461   rval = iter.step();
462   CHECK_ERR(rval);
463   CHECK_EQUAL( 2u, iter.depth() );
464   rval = iter.get_parent_split_plane( p );
465   CHECK_ERR(rval);
466   CHECK_EQUAL( BSPTree::Plane(0,0,1,0), p );
467 
468     // no more leaves
469   rval = iter.step();
470   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
471 }
472 
test_box_tree_create()473 void test_box_tree_create()
474 {
475   Core moab;
476   BSPTree tool( &moab );
477   ErrorCode rval;
478   EntityHandle root;
479   BSPTreeBoxIter iter;
480   BSPTree::Plane p;
481   const double min[3] = { -2, -2, -2 };
482   const double max[3] = {  2,  2,  2 };
483 
484     // create a depth-1 tree
485   rval = tool.create_tree( min, max, root );
486   CHECK_ERR(rval);
487   rval = tool.get_tree_iterator( root, iter );
488   CHECK_ERR(rval);
489 
490     // check initial state of iterator
491   CHECK_EQUAL( &tool, iter.tool() );
492   CHECK_EQUAL(  root, iter.handle() );
493   CHECK_EQUAL(  1u,   iter.depth() );
494 
495     // check initial corner values
496   double corners[8][3], expt[8][3];
497   aabox_corners( min, max, expt );
498   rval = iter.get_box_corners( corners );
499   CHECK_ERR( rval );
500   CHECK( compare_hexes( expt, corners, 1e-6 ) );
501 
502 
503     // Try splits that should fail because they
504     // do not intersect the leaf at all
505   rval = tool.split_leaf( iter, BSPTree::Plane(0,0,1, 4) );
506   CHECK_EQUAL( MB_FAILURE, rval );
507   rval = tool.split_leaf( iter, BSPTree::Plane(0,0,1,-4) );
508   CHECK_EQUAL( MB_FAILURE, rval );
509   rval = tool.split_leaf( iter, BSPTree::Plane(0,1,0, 4) );
510   CHECK_EQUAL( MB_FAILURE, rval );
511   rval = tool.split_leaf( iter, BSPTree::Plane(0,1,0,-4) );
512   CHECK_EQUAL( MB_FAILURE, rval );
513   rval = tool.split_leaf( iter, BSPTree::Plane(1,0,0, 4) );
514   CHECK_EQUAL( MB_FAILURE, rval );
515   rval = tool.split_leaf( iter, BSPTree::Plane(1,0,0,-4) );
516   CHECK_EQUAL( MB_FAILURE, rval );
517   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1,-1, 7) );
518   CHECK_EQUAL( MB_FAILURE, rval );
519   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1,-1,-7) );
520   CHECK_EQUAL( MB_FAILURE, rval );
521   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1,-1, 7) );
522   CHECK_EQUAL( MB_FAILURE, rval );
523   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1,-1,-7) );
524   CHECK_EQUAL( MB_FAILURE, rval );
525   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1,-1, 7) );
526   CHECK_EQUAL( MB_FAILURE, rval );
527   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1,-1,-7) );
528   CHECK_EQUAL( MB_FAILURE, rval );
529   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1,-1, 7) );
530   CHECK_EQUAL( MB_FAILURE, rval );
531   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1,-1,-7) );
532   CHECK_EQUAL( MB_FAILURE, rval );
533 
534 
535     // Try a split that should fail because the
536     // resulting leaf would not be a hexahedron.
537     // Clip each corner twice using planes with opposing normals
538   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1,-1,-4) );
539   CHECK_EQUAL( MB_FAILURE, rval );
540   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1, 1, 4) );
541   CHECK_EQUAL( MB_FAILURE, rval );
542   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1,-1,-4) );
543   CHECK_EQUAL( MB_FAILURE, rval );
544   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1, 1, 4) );
545   CHECK_EQUAL( MB_FAILURE, rval );
546   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1,-1,-4) );
547   CHECK_EQUAL( MB_FAILURE, rval );
548   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1, 1, 4) );
549   CHECK_EQUAL( MB_FAILURE, rval );
550   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1,-1,-4) );
551   CHECK_EQUAL( MB_FAILURE, rval );
552   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1, 1, 4) );
553   CHECK_EQUAL( MB_FAILURE, rval );
554   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1, 1,-4) );
555   CHECK_EQUAL( MB_FAILURE, rval );
556   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1,-1, 4) );
557   CHECK_EQUAL( MB_FAILURE, rval );
558   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1, 1,-4) );
559   CHECK_EQUAL( MB_FAILURE, rval );
560   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1,-1, 4) );
561   CHECK_EQUAL( MB_FAILURE, rval );
562   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1, 1,-4) );
563   CHECK_EQUAL( MB_FAILURE, rval );
564   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1,-1, 4) );
565   CHECK_EQUAL( MB_FAILURE, rval );
566   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1, 1,-4) );
567   CHECK_EQUAL( MB_FAILURE, rval );
568   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1,-1, 4) );
569   CHECK_EQUAL( MB_FAILURE, rval );
570     // Clip each edge
571   rval = tool.split_leaf( iter, BSPTree::Plane(-1,-1, 0,-2) );
572   CHECK_EQUAL( MB_FAILURE, rval );
573   rval = tool.split_leaf( iter, BSPTree::Plane( 1,-1, 0,-2) );
574   CHECK_EQUAL( MB_FAILURE, rval );
575   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 1, 0,-2) );
576   CHECK_EQUAL( MB_FAILURE, rval );
577   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 1, 0,-2) );
578   CHECK_EQUAL( MB_FAILURE, rval );
579   rval = tool.split_leaf( iter, BSPTree::Plane( 0,-1,-1,-2) );
580   CHECK_EQUAL( MB_FAILURE, rval );
581   rval = tool.split_leaf( iter, BSPTree::Plane( 0, 1,-1,-2) );
582   CHECK_EQUAL( MB_FAILURE, rval );
583   rval = tool.split_leaf( iter, BSPTree::Plane( 0, 1, 1,-2) );
584   CHECK_EQUAL( MB_FAILURE, rval );
585   rval = tool.split_leaf( iter, BSPTree::Plane( 0,-1, 1,-2) );
586   CHECK_EQUAL( MB_FAILURE, rval );
587   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 0,-1,-2) );
588   CHECK_EQUAL( MB_FAILURE, rval );
589   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 0,-1,-2) );
590   CHECK_EQUAL( MB_FAILURE, rval );
591   rval = tool.split_leaf( iter, BSPTree::Plane( 1, 0, 1,-2) );
592   CHECK_EQUAL( MB_FAILURE, rval );
593   rval = tool.split_leaf( iter, BSPTree::Plane(-1, 0, 1,-2) );
594   CHECK_EQUAL( MB_FAILURE, rval );
595 
596     // verify that iterator is still valid after many failed splits
597   CHECK_EQUAL( &tool, iter.tool() );
598   CHECK_EQUAL(  root, iter.handle() );
599   CHECK_EQUAL(  1u,   iter.depth() );
600   rval = iter.get_box_corners( corners );
601   CHECK_ERR( rval );
602   CHECK( compare_hexes( expt, corners, 1e-6 ) );
603 
604 
605     // split with z=0
606   rval = tool.split_leaf( iter, BSPTree::Plane(0,0,0.5,0) );
607   CHECK_ERR(rval);
608   CHECK_EQUAL( 2u, iter.depth() );
609 
610     // check that parent split plane is correct
611   rval = iter.get_parent_split_plane( p );
612   CHECK_ERR(rval);
613   CHECK_EQUAL( BSPTree::Plane(0,0,1,0), p );
614 
615     // check that box corners are correct
616   aabox_corners( min, max, expt );
617   for (unsigned i = 0; i < 8; ++i)
618     if (expt[i][2] > 0.0)
619       expt[i][2] = 0.0;
620   rval = iter.get_box_corners( corners );
621   CHECK_ERR( rval );
622   CHECK( compare_hexes( expt, corners, 1e-6 ) );
623 
624 
625     // split with z=-1 and normal in opposite direction
626   rval = tool.split_leaf( iter, BSPTree::Plane(0,0,-2,-2) );
627   CHECK_ERR(rval);
628   CHECK_EQUAL( 3u, iter.depth() );
629   for (unsigned i = 0; i < 8; ++i)
630     if (expt[i][2] < -1.0)
631       expt[i][2] = -1.0;
632   rval = iter.get_box_corners( corners );
633   CHECK_ERR( rval );
634   CHECK( compare_hexes( expt, corners, 1e-6 ) );
635 
636     // step to next leaf (z < -1)
637   rval = iter.step();
638   CHECK_ERR(rval);
639   CHECK_EQUAL( 3u, iter.depth() );
640   aabox_corners( min, max, expt );
641   for (unsigned i = 0; i < 8; ++i)
642     if (expt[i][2] > -1.0)
643       expt[i][2] = -1.0;
644   rval = iter.get_box_corners( corners );
645   CHECK_ERR( rval );
646   CHECK( compare_hexes( expt, corners, 1e-6 ) );
647 
648 
649     // split at x = -1
650   rval = tool.split_leaf( iter, BSPTree::Plane(-0.1,0,0,-0.1) );
651   CHECK_ERR(rval);
652   CHECK_EQUAL( 4u, iter.depth() );
653 
654     // check that parent split plane is correct
655   rval = iter.get_parent_split_plane( p );
656   CHECK_ERR(rval);
657   CHECK_EQUAL( BSPTree::Plane(-1,0,0,-1), p );
658 
659     // check that leaf box is correct
660   aabox_corners( -1, -2, -2, 2, 2, -1, expt );
661   rval = iter.get_box_corners( corners );
662   CHECK_ERR( rval );
663   CHECK( compare_hexes( expt, corners, 1e-6 ) );
664 
665 
666     // split at x = 1
667   rval = tool.split_leaf( iter, BSPTree::Plane(5,0,0,-5) );
668   CHECK_ERR(rval);
669   CHECK_EQUAL( 5u, iter.depth() );
670 
671     // check that parent split plane is correct
672   rval = iter.get_parent_split_plane( p );
673   CHECK_ERR(rval);
674   CHECK_EQUAL( BSPTree::Plane(1,0,0,-1), p );
675 
676     // check that leaf box is correct
677   aabox_corners( -1, -2, -2, 1, 2, -1, expt );
678   rval = iter.get_box_corners( corners );
679   CHECK_ERR( rval );
680   CHECK( compare_hexes( expt, corners, 1e-6 ) );
681 
682 
683     // split at y = -1
684   rval = tool.split_leaf( iter, BSPTree::Plane(0,-1,0,-1) );
685   CHECK_ERR(rval);
686   CHECK_EQUAL( 6u, iter.depth() );
687 
688     // check that parent split plane is correct
689   rval = iter.get_parent_split_plane( p );
690   CHECK_ERR(rval);
691   CHECK_EQUAL( BSPTree::Plane(0,-1,0,-1), p );
692 
693     // check that leaf box is correct
694   aabox_corners( -1, -1, -2, 1, 2, -1, expt );
695   rval = iter.get_box_corners( corners );
696   CHECK_ERR( rval );
697   CHECK( compare_hexes( expt, corners, 1e-6 ) );
698 
699 
700     // split at y = 1
701   rval = tool.split_leaf( iter, BSPTree::Plane(0,1,0,-1) );
702   CHECK_ERR(rval);
703   CHECK_EQUAL( 7u, iter.depth() );
704 
705     // check that parent split plane is correct
706   rval = iter.get_parent_split_plane( p );
707   CHECK_ERR(rval);
708   CHECK_EQUAL( BSPTree::Plane(0,1,0,-1), p );
709 
710     // check that leaf box is correct
711   aabox_corners( -1, -1, -2, 1, 1, -1, expt );
712   rval = iter.get_box_corners( corners );
713   CHECK_ERR( rval );
714   CHECK( compare_hexes( expt, corners, 1e-6 ) );
715 
716 
717 
718     // iterate over tree, verifying
719   rval = tool.get_tree_iterator( root, iter );
720   CHECK_ERR(rval);
721   CHECK_EQUAL( 3u, iter.depth() );
722   rval = iter.get_parent_split_plane( p );
723   CHECK_ERR(rval);
724   CHECK_EQUAL( BSPTree::Plane(0,0,-1,-1), p );
725   aabox_corners( -2, -2, -1, 2, 2, 0, expt );
726   rval = iter.get_box_corners( corners );
727   CHECK_ERR( rval );
728   CHECK( compare_hexes( expt, corners, 1e-6 ) );
729 
730   rval = iter.step();
731   CHECK_ERR(rval);
732   CHECK_EQUAL( 7u, iter.depth() );
733   rval = iter.get_parent_split_plane( p );
734   CHECK_ERR(rval);
735   CHECK_EQUAL( BSPTree::Plane(0,1,0,-1), p );
736   aabox_corners( -1, -1, -2, 1, 1, -1, expt );
737   rval = iter.get_box_corners( corners );
738   CHECK_ERR( rval );
739   CHECK( compare_hexes( expt, corners, 1e-6 ) );
740 
741   rval = iter.step();
742   CHECK_ERR(rval);
743   CHECK_EQUAL( 7u, iter.depth() );
744   rval = iter.get_parent_split_plane( p );
745   CHECK_ERR(rval);
746   CHECK_EQUAL( BSPTree::Plane(0,1,0,-1), p );
747   aabox_corners( -1, 1, -2, 1, 2, -1, expt );
748   rval = iter.get_box_corners( corners );
749   CHECK_ERR( rval );
750   CHECK( compare_hexes( expt, corners, 1e-6 ) );
751 
752   rval = iter.step();
753   CHECK_ERR(rval);
754   CHECK_EQUAL( 6u, iter.depth() );
755   rval = iter.get_parent_split_plane( p );
756   CHECK_ERR(rval);
757   CHECK_EQUAL( BSPTree::Plane(0,-1,0,-1), p );
758   aabox_corners( -1, -2, -2, 1, -1, -1, expt );
759   rval = iter.get_box_corners( corners );
760   CHECK_ERR( rval );
761   CHECK( compare_hexes( expt, corners, 1e-6 ) );
762 
763   rval = iter.step();
764   CHECK_ERR(rval);
765   CHECK_EQUAL( 5u, iter.depth() );
766   rval = iter.get_parent_split_plane( p );
767   CHECK_ERR(rval);
768   CHECK_EQUAL( BSPTree::Plane(1,0,0,-1), p );
769   aabox_corners( 1, -2, -2, 2, 2, -1, expt );
770   rval = iter.get_box_corners( corners );
771   CHECK_ERR( rval );
772   CHECK( compare_hexes( expt, corners, 1e-6 ) );
773 
774   rval = iter.step();
775   CHECK_ERR(rval);
776   CHECK_EQUAL( 4u, iter.depth() );
777   rval = iter.get_parent_split_plane( p );
778   CHECK_ERR(rval);
779   CHECK_EQUAL( BSPTree::Plane(-1,0,0,-1), p );
780   aabox_corners( -2, -2, -2, -1, 2, -1, expt );
781   rval = iter.get_box_corners( corners );
782   CHECK_ERR( rval );
783   CHECK( compare_hexes( expt, corners, 1e-6 ) );
784 
785   rval = iter.step();
786   CHECK_ERR(rval);
787   CHECK_EQUAL( 2u, iter.depth() );
788   rval = iter.get_parent_split_plane( p );
789   CHECK_ERR(rval);
790   CHECK_EQUAL( BSPTree::Plane(0,0,1,0), p );
791   aabox_corners( -2, -2, 0, 2, 2, 2, expt );
792   rval = iter.get_box_corners( corners );
793   CHECK_ERR( rval );
794   CHECK( compare_hexes( expt, corners, 1e-6 ) );
795 
796     // no more leaves
797   rval = iter.step();
798   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
799 }
800 
test_leaf_containing_point_bounded_tree()801 void test_leaf_containing_point_bounded_tree()
802 {
803   Core moab;
804   BSPTree tool( &moab );
805   ErrorCode rval;
806   EntityHandle root;
807   BSPTreeIter iter;
808   BSPTreeBoxIter b_iter;
809   BSPTree::Plane p;
810   EntityHandle h;
811   double expected[8][3], corners[8][3];
812 
813 
814 /*  Build Tree
815 
816   1.0 +---------+--------------+
817       |    A    |              |
818       |         |              |
819   0.7 +---------+      C       |
820       |         |              |
821       |         |              |
822       |    B    |              |
823       |         +--------------+ 0.3
824       |         |      D       |
825       |         |              |
826   0.0 +---------+--------------+
827       0.0       0.4            1.0  */
828 
829   const BSPTree::Plane  X_split(1.0, 0.0, 0.0,-0.4);
830   const BSPTree::Plane AB_split(0.0,-1.0, 0.0, 0.7);
831   const BSPTree::Plane CD_split(0.0,-1.0, 0.0, 0.3);
832 
833 
834   const double min[3] = { 0, 0, 0 };
835   const double max[3] = { 1, 1, 1 };
836   rval = tool.create_tree( min, max, root );
837   CHECK_ERR(rval);
838   rval = tool.get_tree_iterator( root, iter );
839   CHECK_ERR(rval);
840 
841   rval = tool.split_leaf( iter, X_split );
842   CHECK_ERR(rval);
843 
844   rval = tool.split_leaf( iter, AB_split );
845   CHECK_ERR(rval);
846   const EntityHandle A = iter.handle();
847 
848   rval = iter.step();
849   CHECK_ERR(rval);
850   const EntityHandle B = iter.handle();
851 
852   rval = iter.step();
853   CHECK_ERR(rval);
854   rval = tool.split_leaf( iter, CD_split );
855   CHECK_ERR(rval);
856   const EntityHandle C = iter.handle();
857 
858   rval = iter.step();
859   CHECK_ERR(rval);
860   const EntityHandle D = iter.handle();
861 
862 
863     // Test queries inside tree
864 
865   const double A_point[] = { 0.1, 0.8, 0.5 };
866   rval = tool.leaf_containing_point( root, A_point, h );
867   CHECK_ERR(rval);
868   CHECK_EQUAL( A, h );
869   rval = tool.leaf_containing_point( root, A_point, iter );
870   CHECK_ERR(rval);
871   CHECK_EQUAL( A, iter.handle() );
872   rval = iter.get_parent_split_plane( p );
873   CHECK_ERR(rval);
874   CHECK_EQUAL( AB_split, p );
875   rval = tool.leaf_containing_point( root, A_point, b_iter );
876   CHECK_ERR(rval);
877   CHECK_EQUAL( A, b_iter.handle() );
878   aabox_corners( 0.0, 0.7, 0.0, 0.4, 1.0, 1.0, expected );
879   rval = b_iter.get_box_corners( corners );
880   CHECK_ERR( rval );
881   CHECK( compare_hexes( expected, corners, 1e-6 ) );
882 
883   const double B_point[] = { 0.3, 0.1, 0.6 };
884   rval = tool.leaf_containing_point( root, B_point, h );
885   CHECK_ERR(rval);
886   CHECK_EQUAL( B, h );
887   rval = tool.leaf_containing_point( root, B_point, iter );
888   CHECK_ERR(rval);
889   CHECK_EQUAL( B, iter.handle() );
890   rval = iter.get_parent_split_plane( p );
891   CHECK_ERR(rval);
892   CHECK_EQUAL( AB_split, p );
893   rval = tool.leaf_containing_point( root, B_point, b_iter );
894   CHECK_ERR(rval);
895   CHECK_EQUAL( B, b_iter.handle() );
896   aabox_corners( 0.0, 0.0, 0.0, 0.4, 0.7, 1.0, expected );
897   rval = b_iter.get_box_corners( corners );
898   CHECK_ERR( rval );
899   CHECK( compare_hexes( expected, corners, 1e-6 ) );
900 
901   const double C_point[] = { 0.9, 0.9, 0.1 };
902   rval = tool.leaf_containing_point( root, C_point, h );
903   CHECK_ERR(rval);
904   CHECK_EQUAL( C, h );
905   rval = tool.leaf_containing_point( root, C_point, iter );
906   CHECK_ERR(rval);
907   CHECK_EQUAL( C, iter.handle() );
908   rval = iter.get_parent_split_plane( p );
909   CHECK_ERR(rval);
910   CHECK_EQUAL( CD_split, p );
911   rval = tool.leaf_containing_point( root, C_point, b_iter );
912   CHECK_ERR(rval);
913   CHECK_EQUAL( C, b_iter.handle() );
914   aabox_corners( 0.4, 0.3, 0.0, 1.0, 1.0, 1.0, expected );
915   rval = b_iter.get_box_corners( corners );
916   CHECK_ERR( rval );
917   CHECK( compare_hexes( expected, corners, 1e-6 ) );
918 
919   const double D_point[] = { 0.5, 0.2, 0.9 };
920   rval = tool.leaf_containing_point( root, D_point, h );
921   CHECK_ERR(rval);
922   CHECK_EQUAL( D, h );
923   rval = tool.leaf_containing_point( root, D_point, iter );
924   CHECK_ERR(rval);
925   CHECK_EQUAL( D, iter.handle() );
926   rval = iter.get_parent_split_plane( p );
927   CHECK_ERR(rval);
928   CHECK_EQUAL( CD_split, p );
929   rval = tool.leaf_containing_point( root, D_point, b_iter );
930   CHECK_ERR(rval);
931   CHECK_EQUAL( D, b_iter.handle() );
932   aabox_corners( 0.4, 0.0, 0.0, 1.0, 0.3, 1.0, expected );
933   rval = b_iter.get_box_corners( corners );
934   CHECK_ERR( rval );
935   CHECK( compare_hexes( expected, corners, 1e-6 ) );
936 
937 
938     // Try a couple points outside of the tree
939 
940   const double above_pt[] = { 0.5, 0.5, 2.0 };
941   rval = tool.leaf_containing_point( root, above_pt, b_iter );
942   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
943 
944   const double left_pt[] = { -1.0, 0.5, 0.5 };
945   rval = tool.leaf_containing_point( root, left_pt, b_iter );
946   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
947 }
948 
test_leaf_containing_point_unbounded_tree()949 void test_leaf_containing_point_unbounded_tree()
950 {
951   Core moab;
952   BSPTree tool( &moab );
953   ErrorCode rval;
954   EntityHandle root;
955   BSPTreeIter iter;
956   BSPTree::Plane p;
957   EntityHandle h;
958 
959 
960 /*  Build Tree
961 
962     \      |
963      \  C  o (0,4,0)
964       \    |
965        \   |
966         \  |
967      B   \ |         D
968           \|
969    ________o (0,0,0)
970             \
971              \
972        A      \
973                o (1,-2,0)
974                 \
975  */
976   const BSPTree::Plane  X_split( 2.0, 1.0, 0.0, 0.0);
977   const BSPTree::Plane AB_split( 0.0, 1.0, 0.0, 0.0);
978   const BSPTree::Plane CD_split( 1.0, 0.0, 0.0, 0.0);
979 
980 
981   rval = tool.create_tree( root );
982   CHECK_ERR(rval);
983   rval = tool.get_tree_iterator( root, iter );
984   CHECK_ERR(rval);
985 
986   rval = tool.split_leaf( iter, X_split );
987   CHECK_ERR(rval);
988 
989   rval = tool.split_leaf( iter, AB_split );
990   CHECK_ERR(rval);
991   const EntityHandle A = iter.handle();
992 
993   rval = iter.step();
994   CHECK_ERR(rval);
995   const EntityHandle B = iter.handle();
996 
997   rval = iter.step();
998   CHECK_ERR(rval);
999   rval = tool.split_leaf( iter, CD_split );
1000   CHECK_ERR(rval);
1001   const EntityHandle C = iter.handle();
1002 
1003   rval = iter.step();
1004   CHECK_ERR(rval);
1005   const EntityHandle D = iter.handle();
1006 
1007 
1008     // Test queries inside tree
1009 
1010   const double A_point[] = { -1000, -1000, -1000 };
1011   rval = tool.leaf_containing_point( root, A_point, h );
1012   CHECK_ERR(rval);
1013   CHECK_EQUAL( A, h );
1014   rval = tool.leaf_containing_point( root, A_point, iter );
1015   CHECK_ERR(rval);
1016   CHECK_EQUAL( A, iter.handle() );
1017   rval = iter.get_parent_split_plane( p );
1018   CHECK_ERR(rval);
1019   CHECK_EQUAL( AB_split, p );
1020 
1021   const double B_point[] = { -3, 1, 100 };
1022   rval = tool.leaf_containing_point( root, B_point, h );
1023   CHECK_ERR(rval);
1024   CHECK_EQUAL( B, h );
1025   rval = tool.leaf_containing_point( root, B_point, iter );
1026   CHECK_ERR(rval);
1027   CHECK_EQUAL( B, iter.handle() );
1028   rval = iter.get_parent_split_plane( p );
1029   CHECK_ERR(rval);
1030   CHECK_EQUAL( AB_split, p );
1031 
1032   const double C_point[] = { -10, 500, 0 };
1033   rval = tool.leaf_containing_point( root, C_point, h );
1034   CHECK_ERR(rval);
1035   CHECK_EQUAL( C, h );
1036   rval = tool.leaf_containing_point( root, C_point, iter );
1037   CHECK_ERR(rval);
1038   CHECK_EQUAL( C, iter.handle() );
1039   rval = iter.get_parent_split_plane( p );
1040   CHECK_ERR(rval);
1041   CHECK_EQUAL( CD_split, p );
1042 
1043   const double D_point[] = { 10, 500, 0 };
1044   rval = tool.leaf_containing_point( root, D_point, h );
1045   CHECK_ERR(rval);
1046   CHECK_EQUAL( D, h );
1047   rval = tool.leaf_containing_point( root, D_point, iter );
1048   CHECK_ERR(rval);
1049   CHECK_EQUAL( D, iter.handle() );
1050   rval = iter.get_parent_split_plane( p );
1051   CHECK_ERR(rval);
1052   CHECK_EQUAL( CD_split, p );
1053 }
1054 
test_merge_leaf()1055 void test_merge_leaf()
1056 {
1057   Core moab;
1058   BSPTree tool( &moab );
1059   ErrorCode rval;
1060   EntityHandle root;
1061   BSPTreeBoxIter iter;
1062   BSPTree::Plane p;
1063   double expected[8][3], corners[8][3];
1064 
1065 
1066 /*  Build Tree
1067 
1068   1.0 +---------+--------------+
1069       |    A    |              |
1070       |         |              |
1071   0.7 +---------+      C       |
1072       |         |              |
1073       |         |              |
1074       |    B    |              |
1075       |         +--------------+ 0.3
1076       |         |      D       |
1077       |         |              |
1078   0.0 +---------+--------------+
1079       0.0       0.4            1.0  */
1080 
1081   const BSPTree::Plane  X_split(1.0, 0.0, 0.0,-0.4);
1082   const BSPTree::Plane AB_split(0.0,-1.0, 0.0, 0.7);
1083   const BSPTree::Plane CD_split(0.0,-1.0, 0.0, 0.3);
1084 
1085   const double min[3] = { 0, 0, 0 };
1086   const double max[3] = { 1, 1, 1 };
1087   rval = tool.create_tree( min, max, root );
1088   CHECK_ERR(rval);
1089   rval = tool.get_tree_iterator( root, iter );
1090   CHECK_ERR(rval);
1091   rval = tool.split_leaf( iter, X_split );
1092   CHECK_ERR(rval);
1093   const EntityHandle AB = iter.handle();
1094   rval = tool.split_leaf( iter, AB_split );
1095   CHECK_ERR(rval);
1096   rval = iter.step();
1097   CHECK_ERR(rval);
1098   rval = iter.step();
1099   CHECK_ERR(rval);
1100   const EntityHandle CD = iter.handle();
1101   rval = tool.split_leaf( iter, CD_split );
1102   CHECK_ERR(rval);
1103   rval = iter.step();
1104   CHECK_ERR(rval);
1105 
1106   rval = tool.get_tree_iterator( root, iter );
1107   CHECK_ERR(rval);
1108   rval = tool.merge_leaf( iter );
1109   CHECK_ERR(rval);
1110   CHECK_EQUAL( AB, iter.handle() );
1111   CHECK_EQUAL( 2u, iter.depth() );
1112   rval = iter.get_parent_split_plane( p );
1113   CHECK_ERR(rval);
1114   CHECK_EQUAL( X_split, p );
1115   aabox_corners( 0.0, 0.0, 0.0, 0.4, 1.0, 1.0, expected );
1116   rval = iter.get_box_corners( corners );
1117   CHECK_ERR( rval );
1118   CHECK( compare_hexes( expected, corners, 1e-6 ) );
1119 
1120   rval = iter.step();
1121   CHECK_ERR(rval);
1122   rval = iter.step();
1123   CHECK_ERR(rval);
1124   rval = tool.merge_leaf( iter );
1125   CHECK_ERR(rval);
1126   CHECK_EQUAL( CD, iter.handle() );
1127   CHECK_EQUAL( 2u, iter.depth() );
1128   rval = iter.get_parent_split_plane( p );
1129   CHECK_ERR(rval);
1130   CHECK_EQUAL( X_split, p );
1131   aabox_corners( 0.4, 0.0, 0.0, 1.0, 1.0, 1.0, expected );
1132   rval = iter.get_box_corners( corners );
1133   CHECK_ERR( rval );
1134   CHECK( compare_hexes( expected, corners, 1e-6 ) );
1135 
1136   rval = iter.step();
1137   CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
1138 }
1139 
1140 static std::vector<int>
neighbors(const BSPTreeBoxIter & iter,const EntityHandle leaves[8],BSPTreeBoxIter::SideBits side,double epsilon)1141 neighbors( const BSPTreeBoxIter& iter,
1142            const EntityHandle leaves[8],
1143            BSPTreeBoxIter::SideBits side,
1144            double epsilon )
1145 {
1146   std::vector<BSPTreeBoxIter> list;
1147   ErrorCode rval = iter.get_neighbors( side, list, epsilon );
1148   CHECK_ERR(rval);
1149 
1150   std::vector<int> results;
1151   for (size_t i = 0; i < list.size(); ++i)
1152     results.push_back( std::find( leaves, leaves+8, list[i].handle() ) - leaves );
1153   std::sort( results.begin(), results.end() );
1154   return results;
1155 }
1156 
test_box_iter_neighbors()1157 void test_box_iter_neighbors()
1158 {
1159   Core moab;
1160   BSPTree tool( &moab );
1161   ErrorCode rval;
1162   EntityHandle root;
1163   BSPTreeBoxIter iter;
1164   BSPTree::Plane p;
1165 
1166 
1167 /*  Build Tree */
1168 
1169 
1170   const double corners[8][3] = { { 0, 0, 0 },
1171                                  { 8, 0, 0 },
1172                                  { 8, 2, 0 },
1173                                  { 0, 2, 0 },
1174                                  { 0, 0, 1 },
1175                                  { 8, 0, 1 },
1176                                  { 8, 2, 1 },
1177                                  { 0, 2, 1 } };
1178   rval = tool.create_tree( corners, root );
1179   CHECK_ERR(rval);
1180   rval = tool.get_tree_iterator( root, iter );
1181   CHECK_ERR(rval);
1182   EntityHandle leaves[8];
1183 
1184   /* +----------------------------------------+
1185      |                                        |
1186      |                   0*                   |
1187      |                                        |
1188      +----------------------------------------+ */
1189 
1190   p = BSPTree::Plane( 1, 0, 0, -4 );
1191   rval = tool.split_leaf( iter, p );
1192   CHECK_ERR(rval);
1193 
1194   /* +-------------------+--------------------+
1195      |                   |                    |
1196      |         0*        |         1          |
1197      |                   |                    |
1198      +----------------------------------------+ */
1199 
1200   p = BSPTree::Plane( -1, 0, 0, 2 );
1201   rval = tool.split_leaf( iter, p );
1202   CHECK_ERR(rval);
1203 
1204   /* +---------+---------+--------------------+
1205      |         |         |                    |
1206      |    1    |    0*   |         2          |
1207      |         |         |                    |
1208      +----------------------------------------+ */
1209 
1210   p = BSPTree::Plane( 0, 1, 0, -1 );
1211   rval = tool.split_leaf( iter, p );
1212   CHECK_ERR(rval);
1213 
1214   /* +---------+---------+--------------------+
1215      |         |    1    |                    |
1216      |    2    +---------+         3          |
1217      |         |    0*   |                    |
1218      +----------------------------------------+ */
1219 
1220   leaves[0] = iter.handle();
1221   rval = iter.step(); CHECK_ERR(rval);
1222   leaves[1] = iter.handle();
1223   rval = iter.step(); CHECK_ERR(rval);
1224 
1225   /* +---------+---------+--------------------+
1226      |         |    1    |                    |
1227      |    2*   +---------+         3          |
1228      |         |    0    |                    |
1229      +----------------------------------------+ */
1230 
1231   p = BSPTree::Plane( 0, -1, 0, 1 );
1232   rval = tool.split_leaf( iter, p );
1233   CHECK_ERR(rval);
1234 
1235   /* +---------+---------+--------------------+
1236      |    2*   |    1    |                    |
1237      +---------+---------+         4          |
1238      |    3    |    0    |                    |
1239      +----------------------------------------+ */
1240 
1241   leaves[2] = iter.handle();
1242   rval = iter.step(); CHECK_ERR(rval);
1243   leaves[3] = iter.handle();
1244   rval = iter.step(); CHECK_ERR(rval);
1245 
1246   /* +---------+---------+--------------------+
1247      |    2    |    1    |                    |
1248      +---------+---------+         4*         |
1249      |    3    |    0    |                    |
1250      +----------------------------------------+ */
1251 
1252   p = BSPTree::Plane( 0, 1, 0, -1 );
1253   rval = tool.split_leaf( iter, p );
1254   CHECK_ERR(rval);
1255 
1256   /* +---------+---------+--------------------+
1257      |    2    |    1    |         5          |
1258      +---------+---------+--------------------+
1259      |    3    |    0    |         4*         |
1260      +----------------------------------------+ */
1261 
1262   p = BSPTree::Plane( 1, 0, 0, -6 );
1263   rval = tool.split_leaf( iter, p );
1264   CHECK_ERR(rval);
1265 
1266   /* +---------+---------+--------------------+
1267      |    2    |    1    |          6         |
1268      +---------+---------+----------+---------+
1269      |    3    |    0    |    4*    |    5    |
1270      +------------------------------+---------+ */
1271 
1272   leaves[4] = iter.handle();
1273   rval = iter.step(); CHECK_ERR(rval);
1274   leaves[5] = iter.handle();
1275   rval = iter.step(); CHECK_ERR(rval);
1276 
1277   /* +---------+---------+--------------------+
1278      |    2    |    1    |          6*        |
1279      +---------+---------+----------+---------+
1280      |    3    |    0    |     4    |    5    |
1281      +------------------------------+---------+ */
1282 
1283   p = BSPTree::Plane( -1, 0, 0, 6 );
1284   rval = tool.split_leaf( iter, p );
1285   CHECK_ERR(rval);
1286 
1287   /* +---------+---------+--------------------+
1288      |    2    |    1    |     7    |    6    |
1289      +---------+---------+----------+---------+
1290      |    3    |    0    |     4    |    5    |
1291      +------------------------------+---------+ */
1292 
1293   leaves[6] = iter.handle();
1294   rval = iter.step(); CHECK_ERR(rval);
1295   leaves[7] = iter.handle();
1296 
1297 
1298     /* check all neighbors of each leaf */
1299 
1300   rval = tool.get_tree_iterator( root, iter );
1301   CHECK_ERR(rval);
1302 
1303   // NOTE:  Several values in the expected result list are
1304   //        commented out in the tests below.  When the tolerance
1305   //        is greater than zero, the search algorithm may or may
1306   //        not return leaves that are not face-adjacent (e.g. adjacent
1307   //        only along edges or at corners.)  The determining factor
1308   //        for whether or not such a neighbor is returned is which
1309   //        sub-tree of the split that defined the source leaf side
1310   //        the neighbor is on.  The algorithm will not search the subtree
1311   //        of the split that created the side and that contains the
1312   //        source leaf.
1313 
1314 
1315     // check neighbors of leaf 0
1316   std::vector<int> expected;
1317   CHECK_EQUAL( leaves[0], iter.handle() );
1318   expected.clear();
1319   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, 1e-6 ) );
1320   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1321   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1322   expected.push_back( 3 );
1323   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, -1e-6 ) );
1324   expected.insert( expected.begin(), 2 );
1325   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047,  1e-6 ) );
1326   expected.clear(); expected.push_back( 1 );
1327   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, -1e-6 ) );
1328   // See NOTE //  expected.push_back( 2 ); expected.push_back( 7 );
1329   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376,  1e-6 ) );
1330   expected.clear(); expected.push_back( 4 );
1331   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, -1e-6 ) );
1332   expected.push_back( 7 );
1333   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265,  1e-6 ) );
1334 
1335     // check neighbors of leaf 1
1336   CHECK_ERR( iter.step() );
1337   CHECK_EQUAL( leaves[1], iter.handle() );
1338   expected.clear();
1339   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, 1e-6 ) );
1340   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1341   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1342   expected.push_back( 2 );
1343   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, -1e-6 ) );
1344   expected.push_back( 3 );
1345   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047,  1e-6 ) );
1346   expected.clear(); expected.push_back( 0 );
1347   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, -1e-6 ) );
1348   // See NOTE //  expected.push_back( 3 ); expected.push_back( 4 );
1349   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154,  1e-6 ) );
1350   expected.clear(); expected.push_back( 7 );
1351   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, -1e-6 ) );
1352   expected.insert( expected.begin(), 4 );
1353   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265,  1e-6 ) );
1354 
1355   /* +---------+---------+--------------------+
1356      |    2    |    1    |     7    |    6    |
1357      +---------+---------+----------+---------+
1358      |    3    |    0    |     4    |    5    |
1359      +------------------------------+---------+ */
1360 
1361     // check neighbors of leaf 2
1362   CHECK_ERR( iter.step() );
1363   CHECK_EQUAL( leaves[2], iter.handle() );
1364   expected.clear();
1365   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, 1e-6 ) );
1366   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, 1e-6 ) );
1367   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1368   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1369   expected.push_back( 3 );
1370   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, -1e-6 ) );
1371   // See NOTE // expected.insert( expected.begin(), 0 );
1372   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154,  1e-6 ) );
1373   expected.clear(); expected.push_back( 1 );
1374   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, -1e-6 ) );
1375   expected.insert( expected.begin(), 0 );
1376   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265,  1e-6 ) );
1377 
1378     // check neighbors of leaf 3
1379   CHECK_ERR( iter.step() );
1380   CHECK_EQUAL( leaves[3], iter.handle() );
1381   expected.clear();
1382   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, 1e-6 ) );
1383   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, 1e-6 ) );
1384   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1385   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1386   expected.push_back( 2 );
1387   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, -1e-6 ) );
1388   // See NOTE // expected.insert( expected.begin(), 1 );
1389   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376,  1e-6 ) );
1390   expected.clear(); expected.push_back( 0 );
1391   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, -1e-6 ) );
1392   expected.push_back( 1 );
1393   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265,  1e-6 ) );
1394 
1395   /* +---------+---------+--------------------+
1396      |    2    |    1    |     7    |    6    |
1397      +---------+---------+----------+---------+
1398      |    3    |    0    |     4    |    5    |
1399      +------------------------------+---------+ */
1400 
1401     // check neighbors of leaf 4
1402   CHECK_ERR( iter.step() );
1403   CHECK_EQUAL( leaves[4], iter.handle() );
1404   expected.clear();
1405   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, 1e-6 ) );
1406   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1407   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1408   expected.push_back( 0 );
1409   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, -1e-6 ) );
1410   expected.push_back( 1 );
1411   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047,  1e-6 ) );
1412   expected.clear(); expected.push_back( 7 );
1413   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, -1e-6 ) );
1414   expected.insert( expected.begin(), 6 );
1415   // See NOTE // expected.insert( expected.begin(), 1 );
1416   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376,  1e-6 ) );
1417   expected.clear(); expected.push_back( 5 );
1418   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, -1e-6 ) );
1419   // See NOTE // expected.push_back( 6 );
1420   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265,  1e-6 ) );
1421 
1422     // check neighbors of leaf 5
1423   CHECK_ERR( iter.step() );
1424   CHECK_EQUAL( leaves[5], iter.handle() );
1425   expected.clear();
1426   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, 1e-6 ) );
1427   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, 1e-6 ) );
1428   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1429   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1430   expected.push_back( 4 );
1431   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, -1e-6 ) );
1432   // See NOTE // expected.push_back( 7 );
1433   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047,  1e-6 ) );
1434   expected.clear(); expected.push_back( 6 );
1435   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, -1e-6 ) );
1436   expected.push_back( 7 );
1437   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376,  1e-6 ) );
1438 
1439   /* +---------+---------+--------------------+
1440      |    2    |    1    |     7    |    6    |
1441      +---------+---------+----------+---------+
1442      |    3    |    0    |     4    |    5    |
1443      +------------------------------+---------+ */
1444 
1445     // check neighbors of leaf 6
1446   CHECK_ERR( iter.step() );
1447   CHECK_EQUAL( leaves[6], iter.handle() );
1448   expected.clear();
1449   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, 1e-6 ) );
1450   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, 1e-6 ) );
1451   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1452   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1453   expected.push_back( 7 );
1454   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, -1e-6 ) );
1455   // See NOTE // expected.insert( expected.begin(), 4 );
1456   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047,  1e-6 ) );
1457   expected.clear(); expected.push_back( 5 );
1458   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, -1e-6 ) );
1459   expected.insert( expected.begin(), 4 );
1460   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154,  1e-6 ) );
1461 
1462     // check neighbors of leaf 7
1463   CHECK_ERR( iter.step() );
1464   CHECK_EQUAL( leaves[7], iter.handle() );
1465   expected.clear();
1466   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B2376, 1e-6 ) );
1467   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3210, 1e-6 ) );
1468   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B4567, 1e-6 ) );
1469   expected.push_back( 1 );
1470   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047, -1e-6 ) );
1471   expected.insert( expected.begin(), 0 );
1472   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B3047,  1e-6 ) );
1473   expected.clear(); expected.push_back( 4 );
1474   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154, -1e-6 ) );
1475   // See NOTE // expected.insert( expected.begin(), 0 );
1476   expected.push_back( 5 );
1477   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B0154,  1e-6 ) );
1478   expected.clear(); expected.push_back( 6 );
1479   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265, -1e-6 ) );
1480   // See NOTE // expected.insert( expected.begin(), 5 );
1481   CHECK_EQUAL( expected, neighbors( iter, leaves, BSPTreeBoxIter::B1265,  1e-6 ) );
1482 }
1483 
1484 
test_leaf_sibling()1485 void test_leaf_sibling()
1486 {
1487   Core moab;
1488   BSPTree tool( &moab );
1489   ErrorCode rval;
1490   EntityHandle root;
1491   BSPTreeIter iter;
1492 
1493 
1494 /*  Build Tree
1495 
1496   1.0 +---------+--------------+
1497       |    A    |              |
1498       |         |              |
1499   0.7 +---------+      C       |
1500       |         |              |
1501       |         |              |
1502       |    B    |              |
1503       |         +--------------+ 0.3
1504       |         |      D       |
1505       |         |              |
1506   0.0 +---------+--------------+
1507       0.0       0.4            1.0  */
1508 
1509   const BSPTree::Plane  X_split(1.0, 0.0, 0.0,-0.4);
1510   const BSPTree::Plane AB_split(0.0,-1.0, 0.0, 0.7);
1511   const BSPTree::Plane CD_split(0.0,-1.0, 0.0, 0.3);
1512 
1513   const double min[3] = { 0, 0, 0 };
1514   const double max[3] = { 1, 1, 1 };
1515   rval = tool.create_tree( min, max, root );
1516   CHECK_ERR(rval);
1517   rval = tool.get_tree_iterator( root, iter );
1518   CHECK_ERR(rval);
1519   rval = tool.split_leaf( iter, X_split );
1520   CHECK_ERR(rval);
1521   rval = tool.split_leaf( iter, AB_split );
1522   CHECK_ERR(rval);
1523   rval = iter.step();
1524   CHECK_ERR(rval);
1525   rval = iter.step();
1526   CHECK_ERR(rval);
1527   rval = tool.split_leaf( iter, CD_split );
1528   CHECK_ERR(rval);
1529 
1530     // create two iterators for testing
1531   BSPTreeIter iter1, iter2;
1532   rval = tool.get_tree_iterator( root, iter1 );
1533   CHECK_ERR(rval);
1534   rval = tool.get_tree_iterator( root, iter2 );
1535   CHECK_ERR(rval);
1536 
1537 
1538     // iter1 == A, iter2 == A
1539   CHECK( !iter1.is_sibling( iter1 ) );
1540   CHECK( !iter1.is_sibling( iter1.handle() ) );
1541   CHECK( !iter1.is_sibling( iter2 ) );
1542   CHECK(  iter1.sibling_is_forward() );
1543 
1544     // iter1 == A, iter2 == B
1545   rval = iter2.step();
1546   CHECK_ERR(rval);
1547   CHECK(  iter1.is_sibling( iter2 ) );
1548   CHECK(  iter1.is_sibling( iter2.handle() ) );
1549   CHECK(  iter2.is_sibling( iter1 ) );
1550   CHECK(  iter2.is_sibling( iter1.handle() ) );
1551   CHECK( !iter2.sibling_is_forward() );
1552 
1553     // iter1 == A, iter2 == C
1554   rval = iter2.step();
1555   CHECK_ERR(rval);
1556   CHECK( !iter1.is_sibling( iter2 ) );
1557   CHECK( !iter1.is_sibling( iter2.handle() ) );
1558   CHECK( !iter2.is_sibling( iter1 ) );
1559   CHECK( !iter2.is_sibling( iter1.handle() ) );
1560   CHECK(  iter2.sibling_is_forward() );
1561 
1562     // iter1 == B, iter2 == C
1563   rval = iter1.step();
1564   CHECK_ERR(rval);
1565   CHECK( !iter1.is_sibling( iter2 ) );
1566   CHECK( !iter1.is_sibling( iter2.handle() ) );
1567   CHECK( !iter2.is_sibling( iter1 ) );
1568   CHECK( !iter2.is_sibling( iter1.handle() ) );
1569   CHECK( !iter1.sibling_is_forward() );
1570 
1571     // iter1 == D, iter2 == C
1572   rval = iter1.step();
1573   CHECK_ERR(rval);
1574   CHECK_EQUAL( iter1.handle(), iter2.handle() );
1575   rval = iter1.step();
1576   CHECK_ERR(rval);
1577   CHECK(  iter1.is_sibling( iter2 ) );
1578   CHECK(  iter1.is_sibling( iter2.handle() ) );
1579   CHECK(  iter2.is_sibling( iter1 ) );
1580   CHECK(  iter2.is_sibling( iter1.handle() ) );
1581   CHECK( !iter1.sibling_is_forward() );
1582 }
1583 
test_leaf_volume(bool box)1584 void test_leaf_volume( bool box )
1585 {
1586   Core moab;
1587   BSPTree tool( &moab );
1588   ErrorCode rval;
1589   EntityHandle root;
1590   BSPTreeBoxIter b_iter;
1591   BSPTreeIter g_iter;
1592   BSPTreeIter& iter = box ? b_iter : g_iter;
1593 
1594 
1595 /*  Build Tree
1596 
1597   1.0 +---------+--------------+
1598       |    A    |              |
1599       |         |              |
1600   0.7 +---------+      C       |
1601       |         |              |
1602       |         |              |
1603       |    B    |              |
1604       |         +--------------+ 0.3
1605       |         |      D       |
1606       |         |              |
1607   0.0 +---------+--------------+
1608       0.0       0.4            1.0  */
1609 
1610   const BSPTree::Plane  X_split(1.0, 0.0, 0.0,-0.4);
1611   const BSPTree::Plane AB_split(0.0,-1.0, 0.0, 0.7);
1612   const BSPTree::Plane CD_split(0.0,-1.0, 0.0, 0.3);
1613 
1614   const double min[3] = { 0, 0, 0 };
1615   const double max[3] = { 1, 1, 1 };
1616   rval = tool.create_tree( min, max, root );
1617   CHECK_ERR(rval);
1618   rval = tool.get_tree_iterator( root, iter );
1619   CHECK_ERR(rval);
1620   rval = tool.split_leaf( iter, X_split );
1621   CHECK_ERR(rval);
1622   rval = tool.split_leaf( iter, AB_split );
1623   CHECK_ERR(rval);
1624   rval = iter.step();
1625   CHECK_ERR(rval);
1626   rval = iter.step();
1627   CHECK_ERR(rval);
1628   rval = tool.split_leaf( iter, CD_split );
1629   CHECK_ERR(rval);
1630 
1631   // reset iterator
1632   rval = tool.get_tree_iterator( root, iter );
1633   CHECK_ERR(rval);
1634 
1635   // check leaf volumes
1636   CHECK_REAL_EQUAL( 0.12, iter.volume(), 1e-12 ); // A
1637   CHECK_ERR(iter.step());
1638   CHECK_REAL_EQUAL( 0.28, iter.volume(), 1e-12 ); // B
1639   CHECK_ERR(iter.step());
1640   CHECK_REAL_EQUAL( 0.42, iter.volume(), 1e-12 ); // C
1641   CHECK_ERR(iter.step());
1642   CHECK_REAL_EQUAL( 0.18, iter.volume(), 1e-12 ); // D
1643 }
1644 
test_leaf_splits_intersects()1645 void test_leaf_splits_intersects()
1646 {
1647   Core moab;
1648   BSPTree tool( &moab );
1649   ErrorCode rval;
1650   EntityHandle root;
1651   BSPTreeBoxIter iter;
1652   BSPTree::Plane p;
1653 
1654   const double min[3] = { 0, 0, 0 };
1655   const double max[3] = { 1, 2, 3 };
1656   rval = tool.create_tree( min, max, root );
1657   CHECK_ERR(rval);
1658   rval = tool.get_tree_iterator( root, iter );
1659   CHECK_ERR(rval);
1660 
1661   p.set( 1, 0, 0, 1 ); // x == -1
1662   CHECK_EQUAL( BSPTreeBoxIter::MISS, iter.splits( p ) );
1663   CHECK( !iter.intersects( p ) );
1664   p.flip();
1665   CHECK_EQUAL( BSPTreeBoxIter::MISS, iter.splits( p ) );
1666   CHECK( !iter.intersects( p ) );
1667 
1668   p.set( 1, 0, 0, 0 ); // x == 0
1669   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1670   p.flip();
1671   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1672 
1673   p.set( 1, 0, 0, -0.5 ); // x == 0.5
1674   CHECK_EQUAL( BSPTreeBoxIter::SPLIT, iter.splits( p ) );
1675   CHECK( iter.intersects( p ) );
1676   p.flip();
1677   CHECK_EQUAL( BSPTreeBoxIter::SPLIT, iter.splits( p ) );
1678   CHECK( iter.intersects( p ) );
1679 
1680   p.set( 1, 0, 0, -1 ); // x == 1
1681   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1682   p.flip();
1683   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1684 
1685   p.set( 1, 0, 0, -2 ); // x == 2
1686   CHECK_EQUAL( BSPTreeBoxIter::MISS, iter.splits( p ) );
1687   CHECK( !iter.intersects( p ) );
1688   p.flip();
1689   CHECK_EQUAL( BSPTreeBoxIter::MISS, iter.splits( p ) );
1690   CHECK( !iter.intersects( p ) );
1691 
1692   double pt1[3] = { 0, 0, 1.5 };
1693   double pt2[3] = { 1, 0, 1.5 };
1694   double pt3[3] = { 0, 1, 3.0 };
1695   p.set( pt1, pt2, pt3 );
1696   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1697   CHECK( iter.intersects( p ) );
1698   p.flip();
1699   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1700   CHECK( iter.intersects( p ) );
1701 
1702   double pt4[3] = { 0, 2, 2.8 };
1703   p.set( pt1, pt2, pt4 );
1704   CHECK_EQUAL( BSPTreeBoxIter::SPLIT, iter.splits( p ) );
1705   CHECK( iter.intersects( p ) );
1706   p.flip();
1707   CHECK_EQUAL( BSPTreeBoxIter::SPLIT, iter.splits( p ) );
1708   CHECK( iter.intersects( p ) );
1709 
1710   double pta[3] = { 0.8, 0, 0 };
1711   double ptb[3] = { 0, 0.2, 3 };
1712   double ptc[3] = { 0.8, 2, 3 };
1713   p.set( pta, ptb, ptc );
1714   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1715   CHECK( iter.intersects( p ) );
1716   p.flip();
1717   CHECK_EQUAL( BSPTreeBoxIter::NONHEX, iter.splits( p ) );
1718   CHECK( iter.intersects( p ) );
1719 }
1720 
1721 #define CHECK_RAY_XSECTS( PT, DIR, T_IN, T_OUT ) do { \
1722   CHECK(iter.intersect_ray( (PT), (DIR), t_in, t_out )); \
1723   CHECK_REAL_EQUAL( (T_IN), t_in, 1e-6 ); \
1724   CHECK_REAL_EQUAL( (T_OUT), t_out, 1e-6 ); \
1725   } while(false)
1726 
test_leaf_intersects_ray_common(bool box)1727 void test_leaf_intersects_ray_common( bool box )
1728 {
1729   ErrorCode rval;
1730   Core moab;
1731   BSPTree tool( &moab );
1732   double t_in, t_out;
1733 
1734   /** Start with only root box for initial testing **/
1735 
1736   /*  (1,5,-3)   (2,5,-3)
1737             o----o
1738            /:   / \
1739           / :  /    \
1740  (1,5,-1)o----o       \
1741          |  :  \        \
1742          |  :  Y \        \
1743          |  :  ^   \        \
1744          |  :  |     \        \
1745          |  :  +-->X   \        \ (6,1,-3)
1746          |  o./..........\.......o
1747          | . L             \    /
1748          |. Z                \ /
1749          o--------------------o
1750   (1,1,-1)                    (6,1,-1)
1751   */
1752   EntityHandle root;
1753   const double corners[][3] = { { 1, 1, -3 },
1754                                 { 6, 1, -3 },
1755                                 { 2, 5, -3 },
1756                                 { 1, 5, -3 },
1757                                 { 1, 1, -1 },
1758                                 { 6, 1, -1 },
1759                                 { 2, 5, -1 },
1760                                 { 1, 5, -1 } };
1761   rval = tool.create_tree( corners, root );
1762   CHECK_ERR(rval);
1763 
1764   BSPTreeIter gen_iter;
1765   BSPTreeBoxIter box_iter;
1766   BSPTreeIter& iter = box ? static_cast<BSPTreeIter&>(box_iter) : gen_iter;
1767   rval = tool.get_tree_iterator( root, iter );
1768   CHECK_ERR(rval);
1769 
1770     // start with point inside box
1771   const double pt1[] = { 3.5, 3, -2 };
1772   const double dir1[] = { 0.1, 0.1, 0.1 };
1773   CHECK_RAY_XSECTS( pt1, dir1, 0, 2.5 );
1774   const double dir2[] = { 5, 5, 5 };
1775   CHECK_RAY_XSECTS( pt1, dir2, 0, 0.05 );
1776   const double pxdir[] = { 1, 0, 0 };
1777   CHECK_RAY_XSECTS( pt1, pxdir, 0, 0.5 );
1778   const double nxdir[] = { -1, 0, 0 };
1779   CHECK_RAY_XSECTS( pt1, nxdir, 0, 2.5 );
1780   const double pydir[] = { 0, 1, 0 };
1781   CHECK_RAY_XSECTS( pt1, pydir, 0, 0.5 );
1782   const double nydir[] = { 0, -1, 0 };
1783   CHECK_RAY_XSECTS( pt1, nydir, 0, 2 );
1784   const double pzdir[] = { 0, 0, 1 };
1785   CHECK_RAY_XSECTS( pt1, pzdir, 0, 1 );
1786   const double nzdir[] = { 0, 0, -1 };
1787   CHECK_RAY_XSECTS( pt1, nzdir, 0, 1 );
1788 
1789     // point below box
1790   const double pt2[] = { 3.5, 3, -4 };
1791   CHECK(!iter.intersect_ray( pt2, dir1, t_in, t_out ));
1792   CHECK(!iter.intersect_ray( pt2, dir2, t_in, t_out ));
1793   CHECK(!iter.intersect_ray( pt2, pxdir, t_in, t_out ));
1794   CHECK(!iter.intersect_ray( pt2, nxdir, t_in, t_out ));
1795   CHECK(!iter.intersect_ray( pt2, pydir, t_in, t_out ));
1796   CHECK(!iter.intersect_ray( pt2, nydir, t_in, t_out ));
1797   CHECK_RAY_XSECTS( pt2, pzdir, 1, 3 );
1798   CHECK(!iter.intersect_ray( pt2, nzdir, t_in, t_out ));
1799 
1800     // point right of box
1801   const double pt3[] = { 7, 3, -2 };
1802   CHECK(!iter.intersect_ray( pt3, dir1, t_in, t_out ));
1803   CHECK(!iter.intersect_ray( pt3, dir2, t_in, t_out ));
1804   CHECK(!iter.intersect_ray( pt3, pxdir, t_in, t_out ));
1805   CHECK_RAY_XSECTS( pt3, nxdir, 3, 6 );
1806   CHECK(!iter.intersect_ray( pt3, pydir, t_in, t_out ));
1807   CHECK(!iter.intersect_ray( pt3, nydir, t_in, t_out ));
1808   CHECK(!iter.intersect_ray( pt3, pzdir, t_in, t_out ));
1809   CHECK(!iter.intersect_ray( pt3, nzdir, t_in, t_out ));
1810 
1811     // a few more complex test cases
1812   const double dira[] = { -2, -2, 0 };
1813   CHECK_RAY_XSECTS( pt3, dira, 0.75, 1.0 );
1814   const double dirb[] = { -1, -2.1, 0 };
1815   CHECK(!iter.intersect_ray( pt3, dirb, t_in, t_out ));
1816 
1817 
1818   /** Now split twice and test the bottom right corne **/
1819 
1820   BSPTree::Plane Y3( BSPTree::Y, 3.0 );
1821   rval = tool.split_leaf( iter, Y3 );
1822   CHECK_ERR(rval);
1823   BSPTree::Plane X2( BSPTree::X, 2.0 );
1824   rval = tool.split_leaf( iter, X2 );
1825   CHECK_ERR(rval);
1826   rval = iter.step();
1827   CHECK_ERR(rval);
1828 
1829 
1830   /*
1831              (2,3,-3)
1832                  o--------o (4,3,-3)
1833                 /:       /  \
1834                / :      /     \
1835       (2,3,-1)o--------o(4,3,-1)\ (6,1,-3)
1836               |  o.......\.......o
1837               | .          \    /
1838               |.             \ /
1839               o---------------o
1840          (2,1,-1)             (6,1,-1)
1841   */
1842 
1843 
1844     // start with point inside box
1845   const double pt4[] = { 4, 2, -2 };
1846   CHECK_RAY_XSECTS( pt4, dir1, 0, 5 );
1847   CHECK_RAY_XSECTS( pt4, dir2, 0, 0.1 );
1848   CHECK_RAY_XSECTS( pt4, pxdir, 0, 1 );
1849   CHECK_RAY_XSECTS( pt4, nxdir, 0, 2 );
1850   CHECK_RAY_XSECTS( pt4, pydir, 0, 1 );
1851   CHECK_RAY_XSECTS( pt4, nydir, 0, 1 );
1852   CHECK_RAY_XSECTS( pt4, pzdir, 0, 1 );
1853   CHECK_RAY_XSECTS( pt4, nzdir, 0, 1 );
1854 
1855     // point below box
1856   const double pt5[] = { 4, 2, -4 };
1857   CHECK(!iter.intersect_ray( pt5, dir1, t_in, t_out ));
1858   CHECK(!iter.intersect_ray( pt5, dir2, t_in, t_out ));
1859   CHECK(!iter.intersect_ray( pt5, pxdir, t_in, t_out ));
1860   CHECK(!iter.intersect_ray( pt5, nxdir, t_in, t_out ));
1861   CHECK(!iter.intersect_ray( pt5, pydir, t_in, t_out ));
1862   CHECK(!iter.intersect_ray( pt5, nydir, t_in, t_out ));
1863   CHECK_RAY_XSECTS( pt5, pzdir, 1, 3 );
1864   CHECK(!iter.intersect_ray( pt5, nzdir, t_in, t_out ));
1865 
1866     // point right of box
1867   const double pt6[] = { 7, 2, -2 };
1868   CHECK(!iter.intersect_ray( pt6, dir1, t_in, t_out ));
1869   CHECK(!iter.intersect_ray( pt6, dir2, t_in, t_out ));
1870   CHECK(!iter.intersect_ray( pt6, pxdir, t_in, t_out ));
1871   CHECK_RAY_XSECTS( pt6, nxdir, 2, 5 );
1872   CHECK(!iter.intersect_ray( pt6, pydir, t_in, t_out ));
1873   CHECK(!iter.intersect_ray( pt6, nydir, t_in, t_out ));
1874   CHECK(!iter.intersect_ray( pt6, pzdir, t_in, t_out ));
1875   CHECK(!iter.intersect_ray( pt6, nzdir, t_in, t_out ));
1876 
1877     // a few more complex test cases
1878   const double dird[] = { -2, -2, 0 };
1879   CHECK_RAY_XSECTS( pt6, dird, 0.5, 0.5 );
1880   const double dire[] = { -3, -2, 0 };
1881   CHECK_RAY_XSECTS( pt6, dire, 0.4, 0.5 );
1882   const double dirf[] = { -2, -2.1, 0 };
1883   CHECK(!iter.intersect_ray( pt6, dirf, t_in, t_out ));
1884 }
1885 
box(const double pts[],int num_pts,double minpt[3],double maxpt[3])1886 static void box( const double pts[], int num_pts, double minpt[3], double maxpt[3] )
1887 {
1888   minpt[0] = maxpt[0] = pts[0];
1889   minpt[1] = maxpt[1] = pts[1];
1890   minpt[2] = maxpt[2] = pts[2];
1891   for (int i = 1; i < num_pts; ++i) {
1892     for (int d = 0; d < 3; ++d) {
1893       if (pts[3*i+d] < minpt[d])
1894         minpt[d] = pts[3*i+d];
1895       if (pts[3*i+d] > maxpt[d])
1896         maxpt[d] = pts[3*i+d];
1897     }
1898   }
1899 }
1900 
build_tree(const double points[],int num_points,BSPTree & tool)1901 static EntityHandle build_tree( const double points[], int num_points, BSPTree& tool )
1902 {
1903     // create points
1904   ErrorCode rval;
1905   std::vector<EntityHandle> pts(num_points);
1906   for (int i = 0; i < num_points; ++i) {
1907     rval = tool.moab()->create_vertex( points + 3*i, pts[i] );
1908     CHECK_ERR(rval);
1909   }
1910 
1911     // calculate bounding box of tree
1912   double minpt[3], maxpt[3];
1913   box( points, num_points, minpt, maxpt );
1914 
1915     // create initial (1-node) tree
1916   EntityHandle root;
1917   rval = tool.create_tree( minpt, maxpt, root );
1918   CHECK_ERR(rval);
1919 
1920   BSPTreeIter iter;
1921   rval = tool.get_tree_iterator( root, iter );
1922   CHECK_ERR(rval);
1923 
1924   rval = tool.moab()->add_entities( root, &pts[0], pts.size() );
1925   CHECK_ERR(rval);
1926 
1927     // build tree
1928   std::vector<EntityHandle> left_pts, right_pts;
1929   std::vector<CartVect> coords(num_points), tmp_coords;
1930   for (; MB_SUCCESS == rval; rval = iter.step()) {
1931     pts.clear();
1932     rval = tool.moab()->get_entities_by_handle( iter.handle(), pts );
1933     CHECK_ERR(rval);
1934     while (pts.size() > 1) {
1935 
1936       coords.resize(pts.size());
1937       rval = tool.moab()->get_coords( &pts[0], pts.size(), coords[0].array() );
1938       CHECK_ERR(rval);
1939 
1940         // find two points far apart apart
1941       std::vector<CartVect>* ptr;
1942       if (coords.size() < 10)
1943         ptr = &coords;
1944       else {
1945         tmp_coords.resize(16);
1946         CartVect pn, px;
1947         box( coords[0].array(), coords.size(), pn.array(), px.array() );
1948         tmp_coords[8] = pn;
1949         tmp_coords[9] = CartVect( px[0], pn[1], pn[2] );
1950         tmp_coords[10] = CartVect( px[0], px[1], pn[2] );
1951         tmp_coords[11] = CartVect( pn[0], px[1], pn[2] );
1952         tmp_coords[12] = CartVect( pn[0], pn[1], px[2] );
1953         tmp_coords[13] = CartVect( px[0], pn[1], px[2] );
1954         tmp_coords[14] = px;
1955         tmp_coords[15] = CartVect( pn[0], px[1], px[2] );
1956         for (int i = 0; i < 8; ++i) {
1957           tmp_coords[i] = coords[0];
1958           for (size_t j = 1; j < coords.size(); ++j)
1959             if ((coords[j]-tmp_coords[i+8]).length_squared() <
1960                 (tmp_coords[i] - tmp_coords[i+8]).length_squared())
1961               tmp_coords[i] = coords[j];
1962         }
1963         tmp_coords.resize(8);
1964         ptr = &tmp_coords;
1965       }
1966 
1967       size_t pt1 = 0, pt2 = 0;
1968       double lsqr = -1;
1969       for (size_t i = 0; i < ptr->size(); ++i) {
1970         for (size_t j = 0; j < ptr->size();++j) {
1971           double ls = ((*ptr)[i] - (*ptr)[j]).length_squared();
1972           if (ls > lsqr) {
1973             lsqr = ls;
1974             pt1 = i;
1975             pt2 = j;
1976           }
1977         }
1978       }
1979 
1980         // if all points are coincident
1981       if (lsqr <= 1e-12)
1982         break;
1983 
1984         // define normal orthogonal to line through two points
1985       CartVect norm = (*ptr)[pt1] - (*ptr)[pt2];
1986       norm.normalize();
1987 
1988         // find mean position for plane
1989       double coeff = 0.0;
1990       for (size_t i = 0; i < coords.size(); ++i)
1991         coeff -= norm % coords[i];
1992       coeff /= coords.size();
1993 
1994         // left/right sort points
1995       left_pts.clear();
1996       right_pts.clear();
1997       for (size_t i = 0; i < coords.size(); ++i) {
1998         double d = -(norm % coords[i]);
1999         if (d >= coeff)
2000           left_pts.push_back( pts[i] );
2001         else
2002           right_pts.push_back( pts[i] );
2003       }
2004 
2005       rval = tool.split_leaf( iter, BSPTree::Plane( norm.array(), coeff ), left_pts, right_pts );
2006       CHECK_ERR(rval);
2007       CHECK( !left_pts.empty() && !right_pts.empty() );
2008       pts.swap( left_pts );
2009     }
2010 
2011 //    printf("Leaf %d contains %d vertices: ", (int)ID_FROM_HANDLE(iter.handle()),
2012 //                                             (int)(pts.size()) );
2013 //    for (size_t i = 0; i < pts.size(); ++i)
2014 //      printf( "%d, ", (int)ID_FROM_HANDLE(pts[i]));
2015 //    printf("\n");
2016   }
2017 
2018   CHECK(rval == MB_ENTITY_NOT_FOUND);
2019 
2020 
2021     // verify that tree is constructed correctly
2022   for (int i = 0; i < num_points; ++i) {
2023     CartVect pt( points + 3*i );
2024     EntityHandle leaf;
2025     rval = tool.leaf_containing_point( root, pt.array(), leaf );
2026     CHECK_ERR(rval);
2027     Range ents;
2028     rval = tool.moab()->get_entities_by_handle( leaf, ents );
2029     CHECK_ERR(rval);
2030     bool found = false;
2031     for (Range::iterator j = ents.begin(); j != ents.end(); ++j) {
2032       CartVect ent_coords;
2033       rval = tool.moab()->get_coords( &*j, 1, ent_coords.array() );
2034       CHECK_ERR(rval);
2035       if ((pt - ent_coords).length_squared() < 1e-6)
2036         found = true;
2037     }
2038     CHECK(found);
2039   }
2040 
2041   return root;
2042 }
2043 
2044 
test_leaf_polyhedron()2045 void test_leaf_polyhedron()
2046 {
2047     // array of 20 points used to construct tree
2048   static const double points[] = {
2049        7, 6, 3,
2050        5, 3, 5,
2051        9, 2, 6,
2052        7, 2, 1,
2053        3, 9, 0,
2054        6, 0, 6,
2055        1, 6, 2,
2056        9, 7, 8,
2057        2, 0, 2,
2058        5, 7, 3,
2059        2, 2, 9,
2060        7, 9, 8,
2061        1, 6, 3,
2062        3, 9, 2,
2063        4, 9, 1,
2064        4, 8, 7,
2065        3, 0, 5,
2066        0, 1, 6,
2067        2, 3, 6,
2068        1, 6, 0};
2069   const int num_pts = sizeof(points)/(3*sizeof(double));
2070 
2071   ErrorCode rval;
2072   Core moab;
2073   BSPTree tool( &moab );
2074   EntityHandle root = build_tree( points, num_pts, tool );
2075 
2076   BSPTreeIter iter;
2077   rval = tool.get_tree_iterator( root, iter );
2078   CHECK_ERR(rval);
2079 
2080   std::vector<EntityHandle> pts;
2081   std::vector<CartVect> coords;
2082 
2083   for (; rval == MB_SUCCESS; rval = iter.step()) {
2084     BSPTreePoly poly;
2085     rval = iter.calculate_polyhedron( poly );
2086     CHECK_ERR(rval);
2087 
2088     CHECK( poly.is_valid() );
2089     CHECK( poly.volume() > 0.0 );
2090 
2091     pts.clear();
2092     rval = tool.moab()->get_entities_by_handle( iter.handle(), pts );
2093     CHECK_ERR(rval);
2094     CHECK( !pts.empty() );
2095     coords.resize(pts.size());
2096     rval = tool.moab()->get_coords( &pts[0], pts.size(), coords[0].array() );
2097     CHECK_ERR(rval);
2098 
2099     for (size_t i = 0; i < pts.size(); ++i)
2100       CHECK( poly.is_point_contained( coords[i] ) );
2101   }
2102 }
2103