MOAB: Mesh Oriented datABase  (version 5.3.0)
adaptive_kd_tree_tests.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/AdaptiveKDTree.hpp"
00003 #include "moab/CartVect.hpp"
00004 #include "moab/GeomUtil.hpp"
00005 #include "moab/Range.hpp"
00006 #include "TestUtil.hpp"
00007 
00008 using namespace moab;
00009 
00010 #include <iostream>
00011 #include <algorithm>
00012 #include <sstream>
00013 
00014 #ifdef MOAB_HAVE_MPI
00015 #include "moab_mpi.h"
00016 #endif
00017 
00018 /* Utility method - compare two range boxes */
00019 bool box_equal( const AdaptiveKDTreeIter& iter, double x_min, double y_min, double z_min, double x_max, double y_max,
00020                 double z_max )
00021 {
00022     return iter.box_min()[0] == x_min && iter.box_min()[1] == y_min && iter.box_min()[2] == z_min &&
00023            iter.box_max()[0] == x_max && iter.box_max()[1] == y_max && iter.box_max()[2] == z_max;
00024 }
00025 
00026 void build_triangles( Interface* moab, Range& tris, int num_vert, const double* coords, int num_tri,
00027                       const unsigned* conn )
00028 {
00029     std::vector< EntityHandle > verts( num_vert );
00030     for( int i = 0; i < num_vert; ++i )
00031         moab->create_vertex( coords + 3 * i, verts[i] );
00032 
00033     EntityHandle tri, tri_verts[3];
00034     for( int i = 0; i < num_tri; ++i )
00035     {
00036         tri_verts[0] = verts[conn[3 * i]];
00037         tri_verts[1] = verts[conn[3 * i + 1]];
00038         tri_verts[2] = verts[conn[3 * i + 2]];
00039         moab->create_element( MBTRI, tri_verts, 3, tri );
00040         tris.insert( tri );
00041     }
00042 }
00043 
00044 /* Utility method - build a 2x2x2 box composed of two triagles on a side */
00045 void build_triangle_box_small( Interface* moab, Range& tris )
00046 {
00047     const double coords[] = { -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1 };
00048 
00049     const unsigned conn[] = { 0, 1, 5, 0, 5, 4, 2, 6, 7, 2, 7, 3, 1, 2, 6, 1, 6, 5,
00050                               0, 3, 4, 3, 7, 4, 0, 1, 3, 1, 2, 3, 4, 5, 6, 4, 6, 7 };
00051 
00052     build_triangles( moab, tris, 8, coords, 12, conn );
00053 }
00054 
00055 /* build 6x6x6 box centered at origin and composed of 18 triangles per side,
00056    Each side of the box is composed of a 3x3 grid of quads, where each
00057    quad is split diagonally to form two triangles. */
00058 void build_triangle_box_large( Interface* moab, Range& tris )
00059 {
00060     const double coords[] = { // corners
00061                               -3, -3, -3, 3, -3, -3, 3, 3, -3, -3, 3, -3, -3, -3, 3, 3, -3, 3, 3, 3, 3, -3, 3, 3,
00062 
00063                               // edges
00064                               -1, -3, -3, 1, -3, -3, 3, -1, -3, 3, 1, -3, 1, 3, -3, -1, 3, -3, -3, 1, -3, -3, -1, -3,
00065 
00066                               -1, -3, 3, 1, -3, 3, 3, -1, 3, 3, 1, 3, 1, 3, 3, -1, 3, 3, -3, 1, 3, -3, -1, 3,
00067 
00068                               -3, -3, -1, -3, -3, 1, 3, -3, -1, 3, -3, 1, 3, 3, -1, 3, 3, 1, -3, 3, -1, -3, 3, 1,
00069 
00070                               // faces
00071                               -1, -3, -1, 1, -3, -1, 1, -3, 1, -1, -3, 1,
00072 
00073                               3, -1, -1, 3, 1, -1, 3, 1, 1, 3, -1, 1,
00074 
00075                               1, 3, -1, -1, 3, -1, -1, 3, 1, 1, 3, 1,
00076 
00077                               -3, 1, -1, -3, -1, -1, -3, -1, 1, -3, 1, 1,
00078 
00079                               -1, -1, -3, 1, -1, -3, 1, 1, -3, -1, 1, -3,
00080 
00081                               -1, -1, 3, 1, -1, 3, 1, 1, 3, -1, 1, 3
00082     };
00083     const unsigned conn[] = {
00084         // face 0
00085         0, 8, 24, 8, 32, 24, 8, 33, 32, 8, 9, 33, 9, 1, 33, 1, 26, 33, 24, 35, 25, 24, 32, 35, 32, 33, 35, 33, 34, 35,
00086         33, 27, 34, 33, 26, 27, 35, 4, 25, 35, 16, 4, 35, 17, 16, 35, 34, 17, 27, 17, 34, 27, 5, 17,
00087         // face 1
00088         36, 26, 1, 36, 1, 10, 36, 10, 11, 36, 11, 37, 11, 28, 37, 11, 2, 28, 36, 27, 26, 36, 39, 27, 36, 38, 39, 36, 37,
00089         38, 37, 28, 38, 28, 29, 38, 18, 5, 27, 18, 27, 39, 18, 39, 38, 18, 38, 19, 6, 19, 38, 6, 38, 29,
00090         // face 2
00091         12, 28, 2, 12, 40, 28, 12, 41, 40, 12, 13, 41, 3, 41, 13, 3, 30, 41, 43, 29, 28, 43, 28, 40, 43, 40, 41, 43, 41,
00092         42, 41, 31, 42, 41, 30, 31, 43, 6, 29, 43, 20, 6, 43, 21, 20, 43, 42, 21, 21, 42, 31, 21, 31, 7,
00093         // face 3
00094         44, 30, 3, 44, 3, 14, 44, 14, 15, 44, 15, 45, 15, 24, 45, 15, 0, 24, 44, 31, 30, 44, 47, 31, 44, 46, 47, 44, 45,
00095         46, 46, 45, 24, 46, 24, 25, 31, 22, 7, 31, 47, 22, 46, 22, 47, 46, 23, 22, 46, 4, 23, 46, 25, 4,
00096         // face 4
00097         8, 15, 0, 8, 48, 15, 8, 49, 48, 8, 9, 49, 1, 49, 9, 1, 10, 49, 51, 14, 15, 51, 15, 48, 51, 48, 49, 51, 49, 50,
00098         11, 50, 49, 11, 49, 10, 51, 3, 14, 51, 13, 3, 51, 12, 13, 51, 50, 12, 11, 12, 50, 11, 2, 12,
00099         // face 5
00100         4, 52, 16, 4, 23, 52, 22, 52, 23, 22, 55, 52, 22, 21, 55, 22, 7, 21, 17, 16, 52, 17, 52, 53, 54, 53, 52, 54, 52,
00101         55, 54, 55, 21, 54, 21, 20, 18, 5, 17, 18, 17, 53, 18, 53, 54, 18, 54, 19, 6, 19, 54, 6, 54, 20
00102     };
00103 
00104     build_triangles( moab, tris, 56, coords, 108, conn );
00105 }
00106 
00107 /* Utility method - build 2x2x2 octahedron (3D diamond)*/
00108 void build_triangle_octahedron( Interface* moab, Range& tris )
00109 {
00110     const double coords[] = { 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1 };
00111     const unsigned conn[] = { 0, 1, 4, 1, 2, 4, 2, 3, 4, 3, 0, 4, 1, 0, 5, 2, 1, 5, 3, 2, 5, 0, 3, 5 };
00112 
00113     build_triangles( moab, tris, 6, coords, 8, conn );
00114 }
00115 
00116 void test_valid_tree( AdaptiveKDTree* tool, EntityHandle root, FileOptions& options, const Range& expected_tris )
00117 {
00118     Range all_tris;
00119     ErrorCode rval;
00120     AdaptiveKDTreeIter iter;
00121     CHECK( MB_SUCCESS == tool->get_tree_iterator( root, iter ) );
00122     do
00123     {
00124         double dep;
00125         rval = options.get_real_option( "MAX_DEPTH", dep );
00126         CHECK( MB_ENTITY_NOT_FOUND == rval || iter.depth() <= tool->get_max_depth() );
00127 
00128         Range tris;
00129         CHECK( tool->moab()->get_entities_by_type( iter.handle(), MBTRI, tris ) == MB_SUCCESS );
00130         // CHECK( !tris.empty() );
00131         all_tris.merge( tris );
00132 
00133         const CartVect min( iter.box_min() ), max( iter.box_max() );
00134         const CartVect cen( 0.5 * ( min + max ) ), hdim( 0.5 * ( max - min ) );
00135 
00136         for( Range::iterator i = tris.begin(); i != tris.end(); ++i )
00137         {
00138             EntityHandle tri = *i;
00139             const EntityHandle* conn;
00140             int conn_len;
00141             CHECK( tool->moab()->get_connectivity( tri, conn, conn_len ) == MB_SUCCESS );
00142             CHECK( conn_len == 3 );
00143             CartVect coords[3];
00144             CHECK( tool->moab()->get_coords( conn, 3, coords[0].array() ) == MB_SUCCESS );
00145             CHECK( GeomUtil::box_tri_overlap( coords, cen, hdim ) );
00146         }
00147     } while( MB_SUCCESS == ( rval = iter.step() ) );
00148 
00149     CHECK( MB_ENTITY_NOT_FOUND == rval );
00150     CHECK( expected_tris == all_tris );
00151 }
00152 
00153 /* utility method - check that all tris share a vertex */
00154 void check_common_vertex( Interface* moab, const EntityHandle* tris, unsigned num_tri, CartVect point )
00155 {
00156     for( unsigned i = 0; i < num_tri; ++i )
00157         CHECK( MBTRI == moab->type_from_handle( tris[i] ) );
00158 
00159     ErrorCode rval;
00160     CartVect tri_coords[3];
00161     const EntityHandle* conn;
00162     int conn_len;
00163     rval = moab->get_connectivity( tris[0], conn, conn_len );
00164     CHECK( MB_SUCCESS == rval );
00165     rval = moab->get_coords( conn, 3, tri_coords[0].array() );
00166     CHECK( MB_SUCCESS == rval );
00167     tri_coords[0] -= point;
00168     tri_coords[1] -= point;
00169     tri_coords[2] -= point;
00170     EntityHandle vertex = 0;
00171     if( tri_coords[0] % tri_coords[0] < 1e-6 )
00172         vertex = conn[0];
00173     else if( tri_coords[1] % tri_coords[1] < 1e-6 )
00174         vertex = conn[1];
00175     else if( tri_coords[2] % tri_coords[2] < 1e-6 )
00176         vertex = conn[2];
00177     else
00178         CHECK( false );
00179     for( unsigned j = 1; j < num_tri; ++j )
00180     {
00181         rval = moab->get_connectivity( tris[j], conn, conn_len );
00182         CHECK( conn[0] == vertex || conn[1] == vertex || conn[2] == vertex );
00183     }
00184 }
00185 
00186 /* utility method - check that all tris share a vertex */
00187 void check_point_in_triangles( Interface* moab, const EntityHandle* tris, unsigned num_tris, CartVect point )
00188 {
00189     ErrorCode rval;
00190     CartVect tri_coords[3], tript;
00191     const EntityHandle* conn;
00192     int conn_len;
00193 
00194     for( unsigned i = 0; i < num_tris; ++i )
00195     {
00196         CHECK( MBTRI == moab->type_from_handle( tris[i] ) );
00197 
00198         rval = moab->get_connectivity( tris[i], conn, conn_len );
00199         CHECK( MB_SUCCESS == rval );
00200         rval = moab->get_coords( conn, 3, tri_coords[0].array() );
00201         CHECK( MB_SUCCESS == rval );
00202 
00203         GeomUtil::closest_location_on_tri( point, tri_coords, tript );
00204         tript -= point;
00205         CHECK( fabs( tript[0] ) < 1e-6 );
00206         CHECK( fabs( tript[1] ) < 1e-6 );
00207         CHECK( fabs( tript[2] ) < 1e-6 );
00208     }
00209 }
00210 
00211 /* Create the following 2D tree (no Z splits)
00212 
00213         (6) X = -3    (5) X = 0
00214         /             |
00215         /             |      [8]
00216         /             |---------------- (7) Y = 3
00217     [5] /    [6]      |
00218         /             |      [7]
00219         /             |
00220         /             |
00221   -------------------------------------- (1) Y = 0
00222                   |                |
00223         [2]       |                |
00224 (3)----------------|                |
00225 Y = -2             |       [3]      |   [4]
00226                   |                |
00227         [1]       |                |
00228                   |                |
00229                   (2) X = -1      (4) X = 4
00230  */
00231 void create_simple_2d_tree( AdaptiveKDTree& tool, EntityHandle& root, EntityHandle leaves[9] = 0 )
00232 {
00233     ErrorCode rval;
00234     AdaptiveKDTree::Plane plane;
00235 
00236     // create a single-node tree
00237     double min[3] = { -5, -4, -1 };
00238     double max[3] = { 5, 4, 1 };
00239     rval          = tool.create_root( min, max, root );
00240     CHECK( MB_SUCCESS == rval );
00241 
00242     // get iterator for tree
00243     AdaptiveKDTreeIter iter;
00244     rval = tool.get_tree_iterator( root, iter );
00245     CHECK( MB_SUCCESS == rval );
00246     CHECK( box_equal( iter, -5, -4, -1, 5, 4, 1 ) );
00247 
00248     // split plane (1)
00249     plane.norm  = AdaptiveKDTree::Y;
00250     plane.coord = 0.0;
00251     rval        = tool.split_leaf( iter, plane );
00252     CHECK( MB_SUCCESS == rval );
00253     CHECK( box_equal( iter, -5, -4, -1, 5, 0, 1 ) );
00254 
00255     // split plane (2)
00256     plane.norm  = AdaptiveKDTree::X;
00257     plane.coord = -1.0;
00258     rval        = tool.split_leaf( iter, plane );
00259     CHECK( MB_SUCCESS == rval );
00260     CHECK( box_equal( iter, -5, -4, -1, -1, 0, 1 ) );
00261 
00262     // split plane (3), leaf [1]
00263     plane.norm  = AdaptiveKDTree::Y;
00264     plane.coord = -2.0;
00265     rval        = tool.split_leaf( iter, plane );
00266     CHECK( MB_SUCCESS == rval );
00267     CHECK( box_equal( iter, -5, -4, -1, -1, -2, 1 ) );
00268     if( leaves ) leaves[1] = iter.handle();
00269 
00270     // leaf [2]
00271     rval = iter.step();
00272     CHECK( MB_SUCCESS == rval );
00273     CHECK( box_equal( iter, -5, -2, -1, -1, 0, 1 ) );
00274     if( leaves ) leaves[2] = iter.handle();
00275 
00276     // non-leaf right of split plane (2)
00277     rval = iter.step();
00278     CHECK( MB_SUCCESS == rval );
00279     CHECK( box_equal( iter, -1, -4, -1, 5, 0, 1 ) );
00280 
00281     // split plane (4) and leaf [3]
00282     plane.norm  = AdaptiveKDTree::X;
00283     plane.coord = 4.0;
00284     rval        = tool.split_leaf( iter, plane );
00285     CHECK( MB_SUCCESS == rval );
00286     CHECK( box_equal( iter, -1, -4, -1, 4, 0, 1 ) );
00287     if( leaves ) leaves[3] = iter.handle();
00288 
00289     // leaf [4]
00290     rval = iter.step();
00291     CHECK( MB_SUCCESS == rval );
00292     CHECK( box_equal( iter, 4, -4, -1, 5, 0, 1 ) );
00293     if( leaves ) leaves[4] = iter.handle();
00294 
00295     // non-leaf above split plane (1)
00296     rval = iter.step();
00297     CHECK( MB_SUCCESS == rval );
00298     CHECK( box_equal( iter, -5, 0, -1, 5, 4, 1 ) );
00299 
00300     // split plane (5)
00301     plane.norm  = AdaptiveKDTree::X;
00302     plane.coord = 0.0;
00303     rval        = tool.split_leaf( iter, plane );
00304     CHECK( MB_SUCCESS == rval );
00305     CHECK( box_equal( iter, -5, 0, -1, 0, 4, 1 ) );
00306 
00307     // split plane (6) and leaf [5]
00308     plane.norm  = AdaptiveKDTree::X;
00309     plane.coord = -3.0;
00310     rval        = tool.split_leaf( iter, plane );
00311     CHECK( MB_SUCCESS == rval );
00312     CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) );
00313     if( leaves ) leaves[5] = iter.handle();
00314 
00315     // leaf [6];
00316     rval = iter.step();
00317     CHECK( MB_SUCCESS == rval );
00318     CHECK( box_equal( iter, -3, 0, -1, 0, 4, 1 ) );
00319     if( leaves ) leaves[6] = iter.handle();
00320 
00321     // non-leaf right of split plane (5)
00322     rval = iter.step();
00323     CHECK( MB_SUCCESS == rval );
00324     CHECK( box_equal( iter, 0, 0, -1, 5, 4, 1 ) );
00325 
00326     // split plane (7) and leaf [7]
00327     plane.norm  = AdaptiveKDTree::Y;
00328     plane.coord = 3.0;
00329     rval        = tool.split_leaf( iter, plane );
00330     CHECK( MB_SUCCESS == rval );
00331     CHECK( box_equal( iter, 0, 0, -1, 5, 3, 1 ) );
00332     if( leaves ) leaves[7] = iter.handle();
00333 
00334     // leaf [8];
00335     rval = iter.step();
00336     CHECK( MB_SUCCESS == rval );
00337     CHECK( box_equal( iter, 0, 3, -1, 5, 4, 1 ) );
00338     if( leaves ) leaves[8] = iter.handle();
00339 
00340     // the end
00341     rval = iter.step();
00342     CHECK( MB_ENTITY_NOT_FOUND == rval );
00343 }
00344 
00345 void leaf_iterator_test()
00346 {
00347     ErrorCode rval;
00348     Core moab;
00349     Interface* mb = &moab;
00350     AdaptiveKDTree tool( mb );
00351 
00352     // create a single-node tree
00353     EntityHandle root;
00354     double min[3] = { -3, -2, -1 };
00355     double max[3] = { 1, 2, 3 };
00356     rval          = tool.create_root( min, max, root );
00357     CHECK( MB_SUCCESS == rval );
00358 
00359     // get iterator for tree
00360     AdaptiveKDTreeIter iter;
00361     rval = tool.get_tree_iterator( root, iter );
00362     CHECK( MB_SUCCESS == rval );
00363     CHECK( box_equal( iter, -3, -2, -1, 1, 2, 3 ) );
00364 
00365     // check that steping the iterator fails correctly.
00366     rval = iter.step();
00367     CHECK( MB_ENTITY_NOT_FOUND == rval );
00368     rval = iter.step();
00369     CHECK( MB_SUCCESS != rval );
00370 
00371     rval = tool.reset_tree();
00372     CHECK( MB_SUCCESS == rval );
00373     root = 0;
00374 
00375     // construct a simple 2D tree for remaining tests
00376     EntityHandle leaves[9];
00377     create_simple_2d_tree( tool, root, leaves );
00378 
00379     /**** Now traverse tree again, and check neighbors of each leaf *****/
00380     std::vector< AdaptiveKDTreeIter > list;
00381 
00382     // leaf [1]
00383 
00384     rval = tool.get_tree_iterator( root, iter );
00385     CHECK( MB_SUCCESS == rval );
00386     CHECK( box_equal( iter, -5, -4, -1, -1, -2, 1 ) );
00387     CHECK( iter.handle() == leaves[1] );
00388 
00389     list.clear();
00390     rval = iter.get_neighbors( AdaptiveKDTree::X, true, list );
00391     CHECK( MB_SUCCESS == rval );
00392     CHECK( list.empty() );
00393 
00394     list.clear();
00395     rval = iter.get_neighbors( AdaptiveKDTree::X, false, list );
00396     CHECK( MB_SUCCESS == rval );
00397     CHECK( list.size() == 1 );
00398     // should be leaf [3]
00399     CHECK( box_equal( list.front(), -1, -4, -1, 4, 0, 1 ) );
00400     CHECK( list.front().handle() == leaves[3] );
00401 
00402     list.clear();
00403     rval = iter.get_neighbors( AdaptiveKDTree::Y, true, list );
00404     CHECK( MB_SUCCESS == rval );
00405     CHECK( list.empty() );
00406 
00407     list.clear();
00408     rval = iter.get_neighbors( AdaptiveKDTree::Y, false, list );
00409     CHECK( MB_SUCCESS == rval );
00410     CHECK( list.size() == 1 );
00411     // should be leaf [2]
00412     CHECK( box_equal( list.front(), -5, -2, -1, -1, 0, 1 ) );
00413     CHECK( list.front().handle() == leaves[2] );
00414 
00415     list.clear();
00416     rval = iter.get_neighbors( AdaptiveKDTree::Z, true, list );
00417     CHECK( MB_SUCCESS == rval );
00418     CHECK( list.empty() );
00419 
00420     list.clear();
00421     rval = iter.get_neighbors( AdaptiveKDTree::Z, false, list );
00422     CHECK( MB_SUCCESS == rval );
00423     CHECK( list.empty() );
00424 
00425     // leaf [2]
00426 
00427     rval = iter.step();
00428     CHECK( MB_SUCCESS == rval );
00429     CHECK( box_equal( iter, -5, -2, -1, -1, 0, 1 ) );
00430     CHECK( iter.handle() == leaves[2] );
00431 
00432     list.clear();
00433     rval = iter.get_neighbors( AdaptiveKDTree::X, true, list );
00434     CHECK( MB_SUCCESS == rval );
00435     CHECK( list.empty() );
00436 
00437     list.clear();
00438     rval = iter.get_neighbors( AdaptiveKDTree::X, false, list );
00439     CHECK( MB_SUCCESS == rval );
00440     CHECK( list.size() == 1 );
00441     // should be leaf [3]
00442     CHECK( box_equal( list.front(), -1, -4, -1, 4, 0, 1 ) );
00443     CHECK( list.front().handle() == leaves[3] );
00444 
00445     list.clear();
00446     rval = iter.get_neighbors( AdaptiveKDTree::Y, true, list );
00447     CHECK( MB_SUCCESS == rval );
00448     // should be leaf [1]
00449     CHECK( list.size() == 1 );
00450     CHECK( box_equal( list.front(), -5, -4, -1, -1, -2, 1 ) );
00451     CHECK( list.front().handle() == leaves[1] );
00452 
00453     list.clear();
00454     rval = iter.get_neighbors( AdaptiveKDTree::Y, false, list );
00455     CHECK( MB_SUCCESS == rval );
00456     CHECK( list.size() == 2 );
00457     // should be leaf [5] and leaf[6]
00458     if( list[0].handle() == leaves[6] ) std::swap( list[0], list[1] );
00459     CHECK( list[0].handle() == leaves[5] );
00460     CHECK( list[1].handle() == leaves[6] );
00461     CHECK( box_equal( list[0], -5, 0, -1, -3, 4, 1 ) );
00462     CHECK( box_equal( list[1], -3, 0, -1, 0, 4, 1 ) );
00463 
00464     list.clear();
00465     rval = iter.get_neighbors( AdaptiveKDTree::Z, true, list );
00466     CHECK( MB_SUCCESS == rval );
00467     CHECK( list.empty() );
00468 
00469     list.clear();
00470     rval = iter.get_neighbors( AdaptiveKDTree::Z, false, list );
00471     CHECK( MB_SUCCESS == rval );
00472     CHECK( list.empty() );
00473 
00474     // leaf [3]
00475     rval = iter.step();
00476     CHECK( MB_SUCCESS == rval );
00477     CHECK( box_equal( iter, -1, -4, -1, 4, 0, 1 ) );
00478     CHECK( iter.handle() == leaves[3] );
00479 
00480     // leaf [4]
00481     rval = iter.step();
00482     CHECK( MB_SUCCESS == rval );
00483     CHECK( box_equal( iter, 4, -4, -1, 5, 0, 1 ) );
00484     CHECK( iter.handle() == leaves[4] );
00485 
00486     // leaf [5]
00487     rval = iter.step();
00488     CHECK( MB_SUCCESS == rval );
00489     CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) );
00490     CHECK( iter.handle() == leaves[5] );
00491 
00492     // leaf [6]
00493     rval = iter.step();
00494     CHECK( MB_SUCCESS == rval );
00495     CHECK( box_equal( iter, -3, 0, -1, 0, 4, 1 ) );
00496     CHECK( iter.handle() == leaves[6] );
00497 
00498     // leaf [7]
00499     rval = iter.step();
00500     CHECK( MB_SUCCESS == rval );
00501     CHECK( box_equal( iter, 0, 0, -1, 5, 3, 1 ) );
00502     CHECK( iter.handle() == leaves[7] );
00503 
00504     // leaf [8]
00505     rval = iter.step();
00506     CHECK( MB_SUCCESS == rval );
00507     CHECK( box_equal( iter, 0, 3, -1, 5, 4, 1 ) );
00508     CHECK( iter.handle() == leaves[8] );
00509 }
00510 
00511 void test_tree_merge_nodes()
00512 {
00513     // build simple tree for tests
00514     ErrorCode rval;
00515     Core moab;
00516     Interface* mb = &moab;
00517     AdaptiveKDTree tool( mb );
00518     EntityHandle root;
00519     EntityHandle leaves[9];
00520     create_simple_2d_tree( tool, root, leaves );
00521 
00522     // get iterator for tree
00523     AdaptiveKDTreeIter iter;
00524     rval = tool.get_tree_iterator( root, iter );
00525     CHECK( MB_SUCCESS == rval );
00526 
00527     // merge leaves 1 and 2
00528     rval = tool.merge_leaf( iter );
00529     CHECK( MB_SUCCESS == rval );
00530     CHECK( box_equal( iter, -5, -4, -1, -1, 0, 1 ) );
00531 
00532     // merge leaf 1,2 with 3,4 (implicity merges 3 and 4)
00533     rval = tool.merge_leaf( iter );
00534     CHECK( MB_SUCCESS == rval );
00535     CHECK( box_equal( iter, -5, -4, -1, 5, 0, 1 ) );
00536 
00537     // make sure iterator remains valid
00538     rval = iter.step();
00539     CHECK( MB_SUCCESS == rval );
00540     // leaf 5
00541     CHECK( box_equal( iter, -5, 0, -1, -3, 4, 1 ) );
00542     CHECK( iter.handle() == leaves[5] );
00543 }
00544 
00545 void test_build_tree_bisect_triangles()
00546 {
00547     Core moab;
00548     AdaptiveKDTree tool( &moab );
00549     Range box_tris;
00550 
00551     std::ostringstream options;
00552     options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
00553     FileOptions opts( options.str().c_str() );
00554     EntityHandle root;
00555     ErrorCode rval;
00556 
00557     moab.delete_mesh();
00558     box_tris.clear();
00559     build_triangle_box_small( &moab, box_tris );
00560     rval = tool.build_tree( box_tris, &root, &opts );
00561     CHECK( MB_SUCCESS == rval );
00562     test_valid_tree( &tool, root, opts, box_tris );
00563     tool.reset_tree();
00564 
00565     moab.delete_mesh();
00566     box_tris.clear();
00567     build_triangle_octahedron( &moab, box_tris );
00568     rval = tool.build_tree( box_tris, &root, &opts );
00569     CHECK( MB_SUCCESS == rval );
00570     test_valid_tree( &tool, root, opts, box_tris );
00571     tool.reset_tree();
00572 
00573     moab.delete_mesh();
00574     box_tris.clear();
00575     build_triangle_box_large( &moab, box_tris );
00576     rval = tool.build_tree( box_tris, &root, &opts );
00577     CHECK( MB_SUCCESS == rval );
00578     test_valid_tree( &tool, root, opts, box_tris );
00579     tool.reset_tree();
00580 
00581     options << "MAX_DEPTH=2;";
00582     opts = FileOptions( options.str().c_str() );
00583     build_triangle_box_large( &moab, box_tris );
00584     rval = tool.build_tree( box_tris, &root, &opts );
00585     CHECK( MB_SUCCESS == rval );
00586     test_valid_tree( &tool, root, opts, box_tris );
00587     tool.reset_tree();
00588 }
00589 
00590 void test_closest_triangle()
00591 {
00592     Core moab;
00593     AdaptiveKDTree tool( &moab );
00594     Range box_tris;
00595 
00596     std::ostringstream options;
00597     options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
00598     FileOptions opts( options.str().c_str() );
00599     EntityHandle root;
00600     ErrorCode rval;
00601     EntityHandle tri;
00602 
00603     moab.delete_mesh();
00604     box_tris.clear();
00605     build_triangle_box_small( &moab, box_tris );
00606     rval = tool.build_tree( box_tris, &root, &opts );
00607     CHECK( MB_SUCCESS == rval );
00608     test_valid_tree( &tool, root, opts, box_tris );
00609 
00610     // test closest to each corner of the box.
00611     for( unsigned i = 0; i < 8; ++i )
00612     {
00613         int x           = i & 1 ? 1 : -1;
00614         int y           = i & 2 ? 1 : -1;
00615         int z           = i & 4 ? 1 : -1;
00616         double coords[] = { 2.0 * x, 2.0 * y, 2.0 * z };
00617         CartVect point;
00618         rval = tool.closest_triangle( root, coords, point.array(), tri );
00619         CHECK( MB_SUCCESS == rval );
00620         // check result position
00621         CHECK( fabs( x - point[0] ) < 1e-6 );
00622         CHECK( fabs( y - point[1] ) < 1e-6 );
00623         CHECK( fabs( z - point[2] ) < 1e-6 );
00624         // check point is at triangle vertex
00625         check_common_vertex( tool.moab(), &tri, 1, point );
00626     }
00627 
00628     // test location on each face
00629     const CartVect facepts[] = { CartVect( -1.0, 0.5, 0.5 ), CartVect( 1.0, 0.5, 0.5 ),  CartVect( 0.5, -1.0, 0.5 ),
00630                                  CartVect( 0.5, 1.0, 0.5 ),  CartVect( 0.5, 0.5, -1.0 ), CartVect( 0.5, 0.5, 1.0 ) };
00631     for( unsigned i = 0; i < 6; ++i )
00632     {
00633         CartVect point;
00634         rval = tool.closest_triangle( root, facepts[i].array(), point.array(), tri );
00635         CHECK( MB_SUCCESS == rval );
00636         // check result position
00637         const CartVect diff = facepts[i] - point;
00638         CHECK( fabs( diff[0] ) < 1e-6 );
00639         CHECK( fabs( diff[1] ) < 1e-6 );
00640         CHECK( fabs( diff[2] ) < 1e-6 );
00641         // check that point is contained within triangle
00642         check_point_in_triangles( tool.moab(), &tri, 1, point );
00643     }
00644 
00645     // test a location really far from the tree
00646     {
00647         const double far[] = { 0.75, 0.75, 200 };
00648         CartVect point;
00649         rval = tool.closest_triangle( root, far, point.array(), tri );
00650         CHECK( MB_SUCCESS == rval );
00651         CHECK( fabs( point[0] - 0.75 ) < 1e-6 );
00652         CHECK( fabs( point[1] - 0.75 ) < 1e-6 );
00653         CHECK( fabs( point[2] - 1.00 ) < 1e-6 );
00654         // check that point is contained within triangle
00655         check_point_in_triangles( tool.moab(), &tri, 1, point );
00656     }
00657 
00658     // now do it all again with a lot more triangles
00659     tool.reset_tree();
00660     moab.delete_mesh();
00661     box_tris.clear();
00662 
00663     build_triangle_box_large( &moab, box_tris );
00664     rval = tool.build_tree( box_tris, &root, &opts );
00665     CHECK( MB_SUCCESS == rval );
00666     test_valid_tree( &tool, root, opts, box_tris );
00667 
00668     // test closest to each corner of the box.
00669     for( unsigned i = 0; i < 8; ++i )
00670     {
00671         int x           = i & 1 ? 1 : -1;
00672         int y           = i & 2 ? 1 : -1;
00673         int z           = i & 4 ? 1 : -1;
00674         double coords[] = { 4.0 * x, 4.0 * y, 4.0 * z };
00675         CartVect point;
00676         rval = tool.closest_triangle( root, coords, point.array(), tri );
00677         CHECK( MB_SUCCESS == rval );
00678         // check result position
00679         CHECK( fabs( 3.0 * x - point[0] ) < 1e-6 );
00680         CHECK( fabs( 3.0 * y - point[1] ) < 1e-6 );
00681         CHECK( fabs( 3.0 * z - point[2] ) < 1e-6 );
00682         // check point is at triangle vertex
00683         check_common_vertex( tool.moab(), &tri, 1, point );
00684     }
00685 
00686     // test location on each face
00687     const CartVect facepts2[] = { CartVect( -3.0, 0.5, 0.5 ), CartVect( 3.0, 0.5, 0.5 ),  CartVect( 0.5, -3.0, 0.5 ),
00688                                   CartVect( 0.5, 3.0, 0.5 ),  CartVect( 0.5, 0.5, -3.0 ), CartVect( 0.5, 0.5, 3.0 ) };
00689     for( unsigned i = 0; i < 6; ++i )
00690     {
00691         CartVect point;
00692         rval = tool.closest_triangle( root, facepts2[i].array(), point.array(), tri );
00693         CHECK( MB_SUCCESS == rval );
00694         // check result position
00695         const CartVect diff = facepts2[i] - point;
00696         CHECK( fabs( diff[0] ) < 1e-6 );
00697         CHECK( fabs( diff[1] ) < 1e-6 );
00698         CHECK( fabs( diff[2] ) < 1e-6 );
00699         // check that point is contained within triangle
00700         check_point_in_triangles( tool.moab(), &tri, 1, point );
00701     }
00702 
00703     // test a location really far from the tree
00704     {
00705         const double far[] = { 2.75, 2.75, 200 };
00706         CartVect point;
00707         rval = tool.closest_triangle( root, far, point.array(), tri );
00708         CHECK( MB_SUCCESS == rval );
00709         CHECK( fabs( point[0] - 2.75 ) < 1e-6 );
00710         CHECK( fabs( point[1] - 2.75 ) < 1e-6 );
00711         CHECK( fabs( point[2] - 3.00 ) < 1e-6 );
00712         // check that point is contained within triangle
00713         check_point_in_triangles( tool.moab(), &tri, 1, point );
00714     }
00715 }
00716 
00717 void test_sphere_intersect_triangles()
00718 {
00719     Core moab;
00720     AdaptiveKDTree tool( &moab );
00721     Range box_tris;
00722 
00723     std::ostringstream options;
00724     options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
00725     FileOptions opts( options.str().c_str() );
00726     EntityHandle root;
00727     ErrorCode rval;
00728     std::vector< EntityHandle > triangles;
00729 
00730     moab.delete_mesh();
00731     box_tris.clear();
00732     build_triangle_box_small( &moab, box_tris );
00733     rval = tool.build_tree( box_tris, &root, &opts );
00734     CHECK( MB_SUCCESS == rval );
00735     test_valid_tree( &tool, root, opts, box_tris );
00736 
00737     // test closest to each corner of the box.
00738     for( unsigned i = 0; i < 8; ++i )
00739     {
00740         double x        = i & 1 ? 1 : -1;
00741         double y        = i & 2 ? 1 : -1;
00742         double z        = i & 4 ? 1 : -1;
00743         double center[] = { x, y, z };
00744         triangles.clear();
00745         rval = tool.sphere_intersect_triangles( root, center, 0.26, triangles );
00746         CHECK( MB_SUCCESS == rval );
00747         // expect 3 to 6 triangles, depending on the corner
00748         CHECK( triangles.size() >= 3 );
00749         CHECK( triangles.size() <= 6 );
00750         // check point is at the same vertex for each triangle
00751         check_common_vertex( tool.moab(), &triangles[0], triangles.size(), CartVect( center ) );
00752     }
00753 
00754     // now do it all again with a lot more triangles
00755 
00756     tool.reset_tree();
00757     moab.delete_mesh();
00758     box_tris.clear();
00759     build_triangle_box_large( &moab, box_tris );
00760 
00761     rval = tool.build_tree( box_tris, &root, &opts );
00762     CHECK( MB_SUCCESS == rval );
00763     test_valid_tree( &tool, root, opts, box_tris );
00764 
00765     // test closest to each corner of the box.
00766     for( unsigned i = 0; i < 8; ++i )
00767     {
00768         int x           = i & 1 ? 1 : -1;
00769         int y           = i & 2 ? 1 : -1;
00770         int z           = i & 4 ? 1 : -1;
00771         double center[] = { 3.0 * x, 3.0 * y, 3.0 * z };
00772         triangles.clear();
00773         rval = tool.sphere_intersect_triangles( root, center, 0.26, triangles );
00774         CHECK( MB_SUCCESS == rval );
00775         // expect 3 to 6 triangles, depending on the corner
00776         CHECK( triangles.size() >= 3 );
00777         CHECK( triangles.size() <= 6 );
00778         // check point is at the same vertex for each triangle
00779         check_common_vertex( tool.moab(), &triangles[0], triangles.size(), CartVect( center ) );
00780     }
00781 }
00782 
00783 void test_ray_intersect_triangles()
00784 {
00785     Core moab;
00786     AdaptiveKDTree tool( &moab );
00787     Range box_tris;
00788 
00789     std::ostringstream options;
00790     options << "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;PLANE_SET=0;";
00791     FileOptions opts( options.str().c_str() );
00792     EntityHandle root;
00793     ErrorCode rval;
00794     std::vector< EntityHandle > tris;
00795     std::vector< double > dists;
00796 
00797     moab.delete_mesh();
00798     box_tris.clear();
00799     build_triangle_box_small( &moab, box_tris );
00800     rval = tool.build_tree( box_tris, &root, &opts );
00801     CHECK( MB_SUCCESS == rval );
00802     test_valid_tree( &tool, root, opts, box_tris );
00803 
00804     // test ray through box parallel to X axis
00805     CartVect dir( 1, 0, 0 );
00806     CartVect pt( -2, 0.75, 0.75 );
00807     tris.clear();
00808     dists.clear();
00809     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
00810     CHECK( MB_SUCCESS == rval );
00811     CHECK( tris.size() == 3 );
00812     CHECK( dists.size() == tris.size() );
00813     for( unsigned i = 0; i < dists.size(); ++i )
00814     {
00815         CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 );
00816         CartVect tript = pt + dists[i] * dir;
00817         check_point_in_triangles( &moab, &tris[i], 1, tript );
00818     }
00819 
00820     // test ray through box parallel to X axis, closest tri only
00821     tris.clear();
00822     dists.clear();
00823     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 1 );
00824     CHECK( MB_SUCCESS == rval );
00825     CHECK( tris.size() == 1 );
00826     CHECK( dists.size() == tris.size() );
00827     CHECK( fabs( dists[0] - 1 ) < 1e-6 );
00828     check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir );
00829 
00830     // test ray ending within box, parallel to X axis
00831     tris.clear();
00832     dists.clear();
00833     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 0, 2.0 );
00834     CHECK( MB_SUCCESS == rval );
00835     CHECK( tris.size() == 1 );
00836     CHECK( dists.size() == tris.size() );
00837     CHECK( fabs( dists[0] - 1 ) < 1e-6 );
00838     check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir );
00839 
00840     // test ray starting within box parallel to X axis
00841     tris.clear();
00842     dists.clear();
00843     pt   = CartVect( 0, .75, .75 );
00844     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
00845     CHECK( MB_SUCCESS == rval );
00846     CHECK( tris.size() == 2 );
00847     CHECK( dists.size() == tris.size() );
00848     for( unsigned i = 0; i < dists.size(); ++i )
00849     {
00850         CHECK( fabs( dists[i] - 1 ) < 1e-6 );
00851         CartVect tript = pt + dists[i] * dir;
00852         check_point_in_triangles( &moab, &tris[i], 1, tript );
00853     }
00854 
00855     // test skew ray through box
00856     dir = CartVect( 0.5 * sqrt( 2.0 ), 0.5 * sqrt( 2.0 ), 0.0 );
00857     pt  = CartVect( 0, -1.5, 0 );
00858     tris.clear();
00859     dists.clear();
00860     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
00861     CHECK( MB_SUCCESS == rval );
00862     CHECK( tris.size() == 2 );
00863     CHECK( dists.size() == tris.size() );
00864     if( dists[0] < dists[1] )
00865     {
00866         CHECK( fabs( dists[0] - 0.5 * sqrt( 2.0 ) ) < 1e-6 );
00867         CHECK( fabs( dists[1] - sqrt( 2.0 ) ) < 1e-6 );
00868     }
00869     check_point_in_triangles( &moab, &tris[0], 1, pt + dists[0] * dir );
00870     check_point_in_triangles( &moab, &tris[1], 1, pt + dists[1] * dir );
00871 
00872     // test ray through box parallel to -Y axis
00873     dir = CartVect( 0, -1, 0 );
00874     pt  = CartVect( 0, 2, 0 );
00875     tris.clear();
00876     dists.clear();
00877     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
00878     CHECK( MB_SUCCESS == rval );
00879     CHECK( tris.size() == 4 );
00880     CHECK( dists.size() == tris.size() );
00881     for( unsigned i = 0; i < dists.size(); ++i )
00882     {
00883         CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 );
00884         CartVect tript = pt + dists[i] * dir;
00885         check_point_in_triangles( &moab, &tris[i], 1, tript );
00886     }
00887 
00888     // test ray through box parallel to Z axis
00889     dir = CartVect( 0, 0, 1 );
00890     pt  = CartVect( 0, 0, -2 );
00891     tris.clear();
00892     dists.clear();
00893     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists );
00894     CHECK( MB_SUCCESS == rval );
00895     CHECK( tris.size() == 4 );
00896     CHECK( dists.size() == tris.size() );
00897     for( unsigned i = 0; i < dists.size(); ++i )
00898     {
00899         CHECK( fabs( dists[i] - 1 ) < 1e-6 || fabs( dists[i] - 3 ) < 1e-6 );
00900         CartVect tript = pt + dists[i] * dir;
00901         check_point_in_triangles( &moab, &tris[i], 1, tript );
00902     }
00903 
00904     // test ray through box parallel to Z axis, limit 2 first 2 intersections
00905     tris.clear();
00906     dists.clear();
00907     rval = tool.ray_intersect_triangles( root, 1e-6, dir.array(), pt.array(), tris, dists, 2 );
00908     CHECK( MB_SUCCESS == rval );
00909     CHECK( tris.size() == 2 );
00910     CHECK( dists.size() == tris.size() );
00911     for( unsigned i = 0; i < dists.size(); ++i )
00912     {
00913         CHECK( fabs( dists[i] - 1 ) < 1e-6 );
00914         CartVect tript = pt + dists[i] * dir;
00915         check_point_in_triangles( &moab, &tris[i], 1, tript );
00916     }
00917 }
00918 
00919 void test_leaf_volume()
00920 {
00921     Core moab;
00922     AdaptiveKDTree tool( &moab );
00923     EntityHandle root;
00924 
00925     create_simple_2d_tree( tool, root );
00926     AdaptiveKDTreeIter iter;
00927     ErrorCode rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
00928 
00929     CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 );  // 1
00930     CHECK_ERR( iter.step() );
00931     CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 );  // 2
00932     CHECK_ERR( iter.step() );
00933     CHECK_REAL_EQUAL( 40.0, iter.volume(), 1e-12 );  // 3
00934     CHECK_ERR( iter.step() );
00935     CHECK_REAL_EQUAL( 8.0, iter.volume(), 1e-12 );  // 4
00936     CHECK_ERR( iter.step() );
00937     CHECK_REAL_EQUAL( 16.0, iter.volume(), 1e-12 );  // 5
00938     CHECK_ERR( iter.step() );
00939     CHECK_REAL_EQUAL( 24.0, iter.volume(), 1e-12 );  // 6
00940     CHECK_ERR( iter.step() );
00941     CHECK_REAL_EQUAL( 30.0, iter.volume(), 1e-12 );  // 7
00942     CHECK_ERR( iter.step() );
00943     CHECK_REAL_EQUAL( 10.0, iter.volume(), 1e-12 );  // 8
00944 }
00945 
00946 void test_leaf_sibling()
00947 {
00948     ErrorCode rval;
00949     Core moab;
00950     AdaptiveKDTree tool( &moab );
00951     EntityHandle root;
00952 
00953     create_simple_2d_tree( tool, root );
00954     AdaptiveKDTreeIter iter1, iter2;
00955     rval = tool.get_tree_iterator( root, iter1 );CHECK_ERR( rval );
00956     rval = tool.get_tree_iterator( root, iter2 );CHECK_ERR( rval );
00957 
00958     // iter1 == 1, iter2 == 1
00959     CHECK( !iter1.is_sibling( iter1 ) );
00960     CHECK( !iter1.is_sibling( iter1.handle() ) );
00961     CHECK( !iter1.is_sibling( iter2 ) );
00962     CHECK( iter1.sibling_is_forward() );
00963 
00964     // iter1 == 1, iter2 == 2
00965     rval = iter2.step();CHECK_ERR( rval );
00966     CHECK( iter1.is_sibling( iter2 ) );
00967     CHECK( iter1.is_sibling( iter2.handle() ) );
00968     CHECK( iter2.is_sibling( iter1 ) );
00969     CHECK( iter2.is_sibling( iter1.handle() ) );
00970     CHECK( !iter2.sibling_is_forward() );
00971 
00972     // iter1 == 1, iter2 == 3
00973     rval = iter2.step();CHECK_ERR( rval );
00974     CHECK( !iter1.is_sibling( iter2 ) );
00975     CHECK( !iter1.is_sibling( iter2.handle() ) );
00976     CHECK( !iter2.is_sibling( iter1 ) );
00977     CHECK( !iter2.is_sibling( iter1.handle() ) );
00978     CHECK( iter2.sibling_is_forward() );
00979 
00980     // iter1 == 2, iter2 == 3
00981     rval = iter1.step();CHECK_ERR( rval );
00982     CHECK( !iter1.is_sibling( iter2 ) );
00983     CHECK( !iter1.is_sibling( iter2.handle() ) );
00984     CHECK( !iter2.is_sibling( iter1 ) );
00985     CHECK( !iter2.is_sibling( iter1.handle() ) );
00986     CHECK( !iter1.sibling_is_forward() );
00987 
00988     // iter1 == 4, iter2 == 3
00989     rval = iter1.step();CHECK_ERR( rval );
00990     CHECK_EQUAL( iter1.handle(), iter2.handle() );
00991     rval = iter1.step();CHECK_ERR( rval );
00992     CHECK( iter1.is_sibling( iter2 ) );
00993     CHECK( iter1.is_sibling( iter2.handle() ) );
00994     CHECK( iter2.is_sibling( iter1 ) );
00995     CHECK( iter2.is_sibling( iter1.handle() ) );
00996     CHECK( !iter1.sibling_is_forward() );
00997 }
00998 
00999 void test_leaf_intersects_plane()
01000 {
01001     ErrorCode rval;
01002     Core moab;
01003     AdaptiveKDTree tool( &moab );
01004 
01005     EntityHandle root;
01006     const double min[3] = { -5, -4, -1 };
01007     const double max[3] = { 1, 2, 3 };
01008     rval                = tool.create_root( min, max, root );CHECK_ERR( rval );
01009 
01010     AdaptiveKDTreeIter iter;
01011     rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
01012 
01013     AdaptiveKDTree::Plane x1 = { min[0] - 1, AdaptiveKDTree::X };
01014     CHECK( !iter.intersects( x1 ) );
01015     AdaptiveKDTree::Plane x2 = { min[0], AdaptiveKDTree::X };
01016     CHECK( iter.intersects( x2 ) );
01017     AdaptiveKDTree::Plane x3 = { min[0] + 1, AdaptiveKDTree::X };
01018     CHECK( iter.intersects( x3 ) );
01019     AdaptiveKDTree::Plane x4 = { max[0] - 1, AdaptiveKDTree::X };
01020     CHECK( iter.intersects( x4 ) );
01021     AdaptiveKDTree::Plane x5 = { max[0], AdaptiveKDTree::X };
01022     CHECK( iter.intersects( x5 ) );
01023     AdaptiveKDTree::Plane x6 = { max[0] + 1, AdaptiveKDTree::X };
01024     CHECK( !iter.intersects( x6 ) );
01025 
01026     AdaptiveKDTree::Plane y1 = { min[1] - 1, AdaptiveKDTree::Y };
01027     CHECK( !iter.intersects( y1 ) );
01028     AdaptiveKDTree::Plane y2 = { min[1], AdaptiveKDTree::Y };
01029     CHECK( iter.intersects( y2 ) );
01030     AdaptiveKDTree::Plane y3 = { min[1] + 1, AdaptiveKDTree::Y };
01031     CHECK( iter.intersects( y3 ) );
01032     AdaptiveKDTree::Plane y4 = { max[1] - 1, AdaptiveKDTree::Y };
01033     CHECK( iter.intersects( y4 ) );
01034     AdaptiveKDTree::Plane y5 = { max[1], AdaptiveKDTree::Y };
01035     CHECK( iter.intersects( y5 ) );
01036     AdaptiveKDTree::Plane y6 = { max[1] + 1, AdaptiveKDTree::Y };
01037     CHECK( !iter.intersects( y6 ) );
01038 
01039     AdaptiveKDTree::Plane z1 = { min[2] - 1, AdaptiveKDTree::Z };
01040     CHECK( !iter.intersects( z1 ) );
01041     AdaptiveKDTree::Plane z2 = { min[2], AdaptiveKDTree::Z };
01042     CHECK( iter.intersects( z2 ) );
01043     AdaptiveKDTree::Plane z3 = { min[2] + 1, AdaptiveKDTree::Z };
01044     CHECK( iter.intersects( z3 ) );
01045     AdaptiveKDTree::Plane z4 = { max[2] - 1, AdaptiveKDTree::Z };
01046     CHECK( iter.intersects( z4 ) );
01047     AdaptiveKDTree::Plane z5 = { max[2], AdaptiveKDTree::Z };
01048     CHECK( iter.intersects( z5 ) );
01049     AdaptiveKDTree::Plane z6 = { max[2] + 1, AdaptiveKDTree::Z };
01050     CHECK( !iter.intersects( z6 ) );
01051 }
01052 
01053 #define CHECK_RAY_XSECTS( PT, DIR, T_IN, T_OUT )                     \
01054     do                                                               \
01055     {                                                                \
01056         CHECK( iter.intersect_ray( ( PT ), ( DIR ), t_in, t_out ) ); \
01057         CHECK_REAL_EQUAL( ( T_IN ), t_in, 1e-6 );                    \
01058         CHECK_REAL_EQUAL( ( T_OUT ), t_out, 1e-6 );                  \
01059     } while( false )
01060 
01061 void test_leaf_intersects_ray()
01062 {
01063     ErrorCode rval;
01064     Core moab;
01065     AdaptiveKDTree tool( &moab );
01066     double t_in, t_out;
01067 
01068     EntityHandle root;
01069     const double min[3] = { -5, -4, -1 };
01070     const double max[3] = { 1, 2, 3 };
01071     rval                = tool.create_root( min, max, root );CHECK_ERR( rval );
01072 
01073     AdaptiveKDTreeIter iter;
01074     rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
01075 
01076     // start with point inside box
01077     const double pt1[]  = { 0, 0, 0 };
01078     const double dir1[] = { 0.1, 0.1, 0.1 };
01079     CHECK_RAY_XSECTS( pt1, dir1, 0, 10 );
01080     const double dir2[] = { 5, 5, 5 };
01081     CHECK_RAY_XSECTS( pt1, dir2, 0, 0.2 );
01082     const double pxdir[] = { 1, 0, 0 };
01083     CHECK_RAY_XSECTS( pt1, pxdir, 0, 1 );
01084     const double nxdir[] = { -1, 0, 0 };
01085     CHECK_RAY_XSECTS( pt1, nxdir, 0, 5 );
01086     const double pydir[] = { 0, 1, 0 };
01087     CHECK_RAY_XSECTS( pt1, pydir, 0, 2 );
01088     const double nydir[] = { 0, -1, 0 };
01089     CHECK_RAY_XSECTS( pt1, nydir, 0, 4 );
01090     const double pzdir[] = { 0, 0, 1 };
01091     CHECK_RAY_XSECTS( pt1, pzdir, 0, 3 );
01092     const double nzdir[] = { 0, 0, -1 };
01093     CHECK_RAY_XSECTS( pt1, nzdir, 0, 1 );
01094 
01095     // point below box
01096     const double pt2[] = { 0, 0, -2 };
01097     CHECK_RAY_XSECTS( pt2, dir1, 10, 10 );
01098     CHECK_RAY_XSECTS( pt2, dir2, 0.2, 0.2 );
01099     CHECK( !iter.intersect_ray( pt2, pxdir, t_in, t_out ) );
01100     CHECK( !iter.intersect_ray( pt2, nxdir, t_in, t_out ) );
01101     CHECK( !iter.intersect_ray( pt2, pydir, t_in, t_out ) );
01102     CHECK( !iter.intersect_ray( pt2, nydir, t_in, t_out ) );
01103     CHECK_RAY_XSECTS( pt2, pzdir, 1, 5 );
01104     CHECK( !iter.intersect_ray( pt2, nzdir, t_in, t_out ) );
01105 
01106     // point right of box
01107     const double pt3[] = { 3, 0, 0 };
01108     CHECK( !iter.intersect_ray( pt3, dir1, t_in, t_out ) );
01109     CHECK( !iter.intersect_ray( pt3, dir2, t_in, t_out ) );
01110     CHECK( !iter.intersect_ray( pt3, pxdir, t_in, t_out ) );
01111     CHECK_RAY_XSECTS( pt3, nxdir, 2, 8 );
01112     CHECK( !iter.intersect_ray( pt3, pydir, t_in, t_out ) );
01113     CHECK( !iter.intersect_ray( pt3, nydir, t_in, t_out ) );
01114     CHECK( !iter.intersect_ray( pt3, pzdir, t_in, t_out ) );
01115     CHECK( !iter.intersect_ray( pt3, nzdir, t_in, t_out ) );
01116 
01117     // a few more complex test cases
01118     const double dira[] = { -3, 0, 3 };
01119     CHECK_RAY_XSECTS( pt3, dira, 2. / 3., 1.0 );
01120     const double dirb[] = { -2, 0, 3.1 };
01121     CHECK( !iter.intersect_ray( pt3, dirb, t_in, t_out ) );
01122 }
01123 
01124 int main()
01125 {
01126 #ifdef MOAB_HAVE_MPI
01127     int fail = MPI_Init( 0, 0 );
01128     if( fail ) return fail;
01129 #endif
01130 
01131     int error_count = 0;
01132 
01133     error_count += RUN_TEST( leaf_iterator_test );
01134     error_count += RUN_TEST( test_build_tree_bisect_triangles );
01135     error_count += RUN_TEST( test_closest_triangle );
01136     error_count += RUN_TEST( test_sphere_intersect_triangles );
01137     error_count += RUN_TEST( test_ray_intersect_triangles );
01138     error_count += RUN_TEST( test_tree_merge_nodes );
01139     error_count += RUN_TEST( test_leaf_volume );
01140     error_count += RUN_TEST( test_leaf_sibling );
01141     error_count += RUN_TEST( test_leaf_intersects_plane );
01142     error_count += RUN_TEST( test_leaf_intersects_ray );
01143 
01144 #ifdef MOAB_HAVE_MPI
01145     fail = MPI_Finalize();
01146     if( fail ) return fail;
01147 #endif
01148 
01149     return error_count;
01150 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines