MOAB: Mesh Oriented datABase  (version 5.3.0)
kd_tree_test.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/AdaptiveKDTree.hpp"
00003 #include "moab/Range.hpp"
00004 #include "moab/CartVect.hpp"
00005 
00006 #ifdef MOAB_HAVE_MPI
00007 #include "moab_mpi.h"
00008 #endif
00009 
00010 #include <cmath>
00011 #include <cassert>
00012 #include <cfloat>
00013 #include <cstdio>
00014 
00015 #include "TestUtil.hpp"
00016 
00017 using namespace moab;
00018 
00019 const unsigned INTERVALS = 4;
00020 const unsigned DEPTH     = 7;  // 3*log2(INTERVALS)+1
00021 const char* TAG_NAME     = "TEST_DATA";
00022 
00023 EntityHandle create_tree( AdaptiveKDTree& tool, unsigned depth, int intervals, Tag* tag_handle = 0 );
00024 void validate_tree( AdaptiveKDTree& tool, EntityHandle root, int depth, double intervals );
00025 
00026 void test_tree_create();
00027 void test_leaf_merge();
00028 void test_tree_readwrite();
00029 void test_tree_delete();
00030 void test_iterator_back();
00031 void test_point_search();
00032 
00033 int main( int argc, char** argv )
00034 {
00035 #ifdef MOAB_HAVE_MPI
00036     int fail = MPI_Init( &argc, &argv );
00037     if( fail ) return fail;
00038 #else
00039     // silence the warning of parameters not used, in serial; there should be a smarter way :(
00040     argv[0] = argv[argc - argc];
00041 #endif
00042 
00043     int err = RUN_TEST( test_tree_create );
00044     if( err )  // can't run other tests if can't create tree
00045         return 1;
00046     err += RUN_TEST( test_leaf_merge );
00047 #ifdef MOAB_HAVE_HDF5
00048     err += RUN_TEST( test_tree_readwrite );
00049 #endif
00050     err += RUN_TEST( test_tree_delete );
00051     err += RUN_TEST( test_iterator_back );
00052     err += RUN_TEST( test_point_search );
00053 
00054 #ifdef MOAB_HAVE_MPI
00055     fail = MPI_Finalize();
00056     if( fail ) return fail;
00057 #endif
00058 
00059     return err;
00060 }
00061 
00062 EntityHandle create_tree( AdaptiveKDTree& tool, unsigned depth, int intervals, Tag* tag_handle )
00063 {
00064     // Create tree root
00065     ErrorCode err;
00066     EntityHandle root, leaf;
00067     err = tool.create_root( CartVect( 0.0 ).array(), CartVect( intervals ).array(), root );
00068     assert( !err );
00069 
00070     // Use iterator to create tree to fixed depth of DEPTH
00071     AdaptiveKDTreeIter iter;
00072     err = tool.get_tree_iterator( root, iter );CHECK_ERR( err );
00073     while( err == MB_SUCCESS )
00074     {
00075         if( iter.depth() < depth )
00076         {
00077             // bisect leaves along alternating axes
00078             AdaptiveKDTree::Plane split;
00079             split.norm  = iter.depth() % 3;  // alternate split axes;
00080             split.coord = 0.5 * ( iter.box_min()[split.norm] + iter.box_max()[split.norm] );
00081             err         = tool.split_leaf( iter, split );  // advances iter to first new leaf
00082             CHECK_ERR( err );
00083         }
00084         // if current leaf is at desired depth, advance to next one
00085         else
00086         {
00087             err = iter.step();
00088         }
00089     }
00090     CHECK( MB_ENTITY_NOT_FOUND == err );
00091 
00092     if( !tag_handle ) return root;
00093 
00094     // define a tag to use to store integer values on tree leaves
00095     err = tool.moab()->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, *tag_handle, MB_TAG_DENSE | MB_TAG_EXCL );CHECK_ERR( err );
00096 
00097     // iterate over tree setting data
00098     int counter = 0;
00099     for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
00100     {
00101         // store integer value on leaf
00102         ++counter;
00103         leaf = iter.handle();
00104         err  = tool.moab()->tag_set_data( *tag_handle, &leaf, 1, &counter );CHECK_ERR( err );
00105     }
00106 
00107     return root;
00108 }
00109 
00110 void validate_tree( AdaptiveKDTree& tool, EntityHandle root, unsigned depth, int intervals, Tag data )
00111 {
00112     ErrorCode err;
00113     const double VOL = 1.0;  // all leaves should be 1x1x1 boxes
00114     int val;
00115 
00116     // iterate over tree, verifying leaves
00117     AdaptiveKDTreeIter iter;
00118     int counter = 0;
00119     for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
00120     {
00121         // store integer value on leaf
00122         ++counter;
00123         EntityHandle leaf = iter.handle();
00124         CHECK( leaf != 0 );
00125         CHECK_EQUAL( MBENTITYSET, TYPE_FROM_HANDLE( leaf ) );
00126 
00127         // check size of leaf
00128         const double* min = iter.box_min();
00129         const double* max = iter.box_max();
00130         double dims[]     = { max[0] - min[0], max[1] - min[1], max[2] - min[2] };
00131         double volume     = dims[0] * dims[1] * dims[2];
00132         CHECK_REAL_EQUAL( VOL, volume, DBL_EPSILON );
00133 
00134         // check depth of leaf
00135         CHECK_EQUAL( depth, iter.depth() );
00136 
00137         // check tag value on leaf
00138         err = tool.moab()->tag_get_data( data, &leaf, 1, &val );CHECK_ERR( err );
00139         CHECK_EQUAL( counter, val );
00140     }
00141     // check number of leaves
00142     const int num_leaves = intervals * intervals * intervals;
00143     CHECK_EQUAL( num_leaves, counter );
00144 }
00145 
00146 void test_tree_create()
00147 {
00148     Tag tag;
00149     Core mb;
00150     AdaptiveKDTree tool( &mb );
00151     const EntityHandle root = create_tree( tool, DEPTH, INTERVALS, &tag );
00152     validate_tree( tool, root, DEPTH, INTERVALS, tag );
00153 }
00154 
00155 void test_leaf_merge()
00156 {
00157     ErrorCode err;
00158     Core mb;
00159     AdaptiveKDTree tool( &mb );
00160     Tag data;
00161     const EntityHandle root = create_tree( tool, DEPTH, INTERVALS, &data );
00162 
00163     // reduce tree depth to DEPTH-1 by merging adjacent leaf pairs,
00164     // make new "leaf" have smaller of two data values on original pair
00165     AdaptiveKDTreeIter iter;
00166     for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
00167     {
00168         // get data for first leaf
00169         int data1;
00170         EntityHandle leaf = iter.handle();
00171         err               = mb.tag_get_data( data, &leaf, 1, &data1 );CHECK_ERR( err );
00172         // tree traversal is always such that two leaves with same parent are consective
00173         err = iter.step();CHECK_ERR( err );
00174         // get data for sibling
00175         int data2;
00176         leaf = iter.handle();
00177         err  = mb.tag_get_data( data, &leaf, 1, &data2 );CHECK_ERR( err );
00178         // as we stored increasing values, these had better be increasing
00179         CHECK_EQUAL( 1, data2 - data1 );
00180         // merge leaf pair (iter can be at either one)
00181         err = tool.merge_leaf( iter );  // changes iter to be new "merged" leaf
00182         CHECK_ERR( err );
00183         // store smaller of two values on new leaf
00184         leaf = iter.handle();
00185         err  = mb.tag_set_data( data, &leaf, 1, &data1 );CHECK_ERR( err );
00186     }
00187 
00188     // Iterate over tree, verifying leaves and checking data
00189     // Initial leaves had volume of 1 : merged pairs of leaves so volume should now be 2.
00190     // Initial leaves were enumerated in order : merged pairs so new leaves should
00191     //   have data incrementing in steps of 2.
00192     int counter = 1;
00193     for( err = tool.get_tree_iterator( root, iter ); !err; err = iter.step() )
00194     {
00195         // store integer value on leaf
00196         int data1;
00197         EntityHandle leaf = iter.handle();
00198         err               = mb.tag_get_data( data, &leaf, 1, &data1 );CHECK_ERR( err );
00199         CHECK_EQUAL( counter, data1 );
00200         counter += 2;
00201 
00202         // check size of leaf
00203         const double* min = iter.box_min();
00204         const double* max = iter.box_max();
00205         double dims[]     = { max[0] - min[0], max[1] - min[1], max[2] - min[2] };
00206         double volume     = dims[0] * dims[1] * dims[2];
00207         CHECK_REAL_EQUAL( 2.0, volume, DBL_EPSILON );
00208 
00209         // check depth of leaf
00210         CHECK_EQUAL( DEPTH - 1, iter.depth() );
00211     }
00212 }
00213 
00214 void test_tree_readwrite()
00215 {
00216     ErrorCode err;
00217     Tag tag;
00218     Core mb;
00219     AdaptiveKDTree tool( &mb );
00220     create_tree( tool, DEPTH, INTERVALS, &tag );
00221 
00222     // write to file
00223     err = mb.write_file( "tree.h5m" );CHECK_ERR( err );
00224 
00225     // clear everything
00226     mb.delete_mesh();
00227 
00228     // read tree from file
00229     err = mb.load_file( "tree.h5m" );
00230     remove( "tree.h5m" );CHECK_ERR( err );
00231 
00232     // get tag handle by name, because the handle may have changed
00233     err = mb.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag );CHECK_ERR( err );
00234 
00235     // get root handle for tree
00236     Range range;
00237     err = tool.find_all_trees( range );
00238     assert( !err );
00239     assert( range.size() == 1 );
00240     EntityHandle root = range.front();  // first (only) handle
00241 
00242     validate_tree( tool, root, DEPTH, INTERVALS, tag );
00243 }
00244 
00245 void test_tree_delete()
00246 {
00247     ErrorCode err;
00248     Core mb;
00249     AdaptiveKDTree tool( &mb );
00250     Tag data;
00251     create_tree( tool, DEPTH, INTERVALS, &data );
00252 
00253     err = tool.reset_tree();CHECK_ERR( err );
00254 
00255     Range ents;
00256     err = mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &data, 0, 1, ents );CHECK_ERR( err );
00257     CHECK( ents.empty() );
00258 }
00259 
00260 void test_iterator_back()
00261 {
00262     Core mb;
00263     AdaptiveKDTree tool( &mb );
00264     const EntityHandle root = create_tree( tool, DEPTH, INTERVALS );
00265 
00266     AdaptiveKDTreeIter iter;
00267     ErrorCode rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
00268 
00269     CartVect min( iter.box_min() );
00270     CartVect max( iter.box_max() );
00271     EntityHandle leaf = iter.handle();
00272 
00273     // going back from first location should fail.
00274     rval = iter.back();
00275     CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
00276     rval = tool.get_tree_iterator( root, iter );CHECK_ERR( rval );
00277 
00278     // make sure iterator is valid
00279     CHECK_REAL_EQUAL( min[0], iter.box_min()[0], DBL_EPSILON );
00280     CHECK_REAL_EQUAL( min[1], iter.box_min()[1], DBL_EPSILON );
00281     CHECK_REAL_EQUAL( min[2], iter.box_min()[2], DBL_EPSILON );
00282     CHECK_REAL_EQUAL( max[0], iter.box_max()[0], DBL_EPSILON );
00283     CHECK_REAL_EQUAL( max[1], iter.box_max()[1], DBL_EPSILON );
00284     CHECK_REAL_EQUAL( max[2], iter.box_max()[2], DBL_EPSILON );
00285     CHECK_EQUAL( leaf, iter.handle() );
00286 
00287     while( MB_SUCCESS == iter.step() )
00288     {
00289         // Get values at current iterator location
00290         CartVect next_min( iter.box_min() );
00291         CartVect next_max( iter.box_max() );
00292         EntityHandle next_leaf = iter.handle();
00293 
00294         // step back to previous location
00295         rval = iter.back();CHECK_ERR( rval );
00296 
00297         // check expected values for previous location
00298         CHECK_REAL_EQUAL( min[0], iter.box_min()[0], DBL_EPSILON );
00299         CHECK_REAL_EQUAL( min[1], iter.box_min()[1], DBL_EPSILON );
00300         CHECK_REAL_EQUAL( min[2], iter.box_min()[2], DBL_EPSILON );
00301         CHECK_REAL_EQUAL( max[0], iter.box_max()[0], DBL_EPSILON );
00302         CHECK_REAL_EQUAL( max[1], iter.box_max()[1], DBL_EPSILON );
00303         CHECK_REAL_EQUAL( max[2], iter.box_max()[2], DBL_EPSILON );
00304         CHECK_EQUAL( leaf, iter.handle() );
00305 
00306         // advance iterator to 'current' location
00307         rval = iter.step();CHECK_ERR( rval );
00308 
00309         // check that iterator values are correct
00310         CHECK_REAL_EQUAL( next_min[0], iter.box_min()[0], DBL_EPSILON );
00311         CHECK_REAL_EQUAL( next_min[1], iter.box_min()[1], DBL_EPSILON );
00312         CHECK_REAL_EQUAL( next_min[2], iter.box_min()[2], DBL_EPSILON );
00313         CHECK_REAL_EQUAL( next_max[0], iter.box_max()[0], DBL_EPSILON );
00314         CHECK_REAL_EQUAL( next_max[1], iter.box_max()[1], DBL_EPSILON );
00315         CHECK_REAL_EQUAL( next_max[2], iter.box_max()[2], DBL_EPSILON );
00316         CHECK_EQUAL( next_leaf, iter.handle() );
00317 
00318         // store values for next iteration
00319         min  = next_min;
00320         max  = next_max;
00321         leaf = next_leaf;
00322     }
00323 }
00324 
00325 void test_point_search()
00326 {
00327     Core mb;
00328     AdaptiveKDTree tool( &mb );
00329     const EntityHandle root = create_tree( tool, DEPTH, INTERVALS );
00330 
00331     ErrorCode rval;
00332     EntityHandle leaf;
00333     AdaptiveKDTreeIter iter, iter2;
00334 
00335     // points to test
00336     CartVect left( 0.5 );
00337     CartVect right( CartVect( INTERVALS ) - left );
00338 
00339     // compare leaf search to iterator search
00340     rval = tool.point_search( left.array(), leaf );CHECK_ERR( rval );
00341     rval = tool.point_search( left.array(), iter );CHECK_ERR( rval );
00342     CHECK_EQUAL( leaf, iter.handle() );
00343 
00344     // iterator should be at 'first' leaf
00345     rval = tool.get_tree_iterator( root, iter2 );CHECK_ERR( rval );
00346     for( ;; )
00347     {
00348         CHECK_EQUAL( iter.handle(), iter2.handle() );
00349         CHECK_EQUAL( iter.depth(), iter2.depth() );
00350         CHECK_REAL_EQUAL( iter.box_min()[0], iter2.box_min()[0], DBL_EPSILON );
00351         CHECK_REAL_EQUAL( iter.box_min()[1], iter2.box_min()[1], DBL_EPSILON );
00352         CHECK_REAL_EQUAL( iter.box_min()[2], iter2.box_min()[2], DBL_EPSILON );
00353         CHECK_REAL_EQUAL( iter.box_max()[0], iter2.box_max()[0], DBL_EPSILON );
00354         CHECK_REAL_EQUAL( iter.box_max()[1], iter2.box_max()[1], DBL_EPSILON );
00355         CHECK_REAL_EQUAL( iter.box_max()[2], iter2.box_max()[2], DBL_EPSILON );
00356 
00357         rval = iter2.step();
00358         if( MB_ENTITY_NOT_FOUND == rval ) break;CHECK_ERR( rval );
00359         rval = iter.step();CHECK_ERR( rval );
00360     }
00361 
00362     // compare leaf search to iterator search
00363     rval = tool.point_search( right.array(), leaf, 0.0, 0.0, NULL, const_cast< EntityHandle* >( &root ) );CHECK_ERR( rval );
00364     rval = tool.point_search( right.array(), iter, 0.0, 0.0, NULL, const_cast< EntityHandle* >( &root ) );CHECK_ERR( rval );
00365     assert( iter.handle() == leaf );
00366 
00367     // iterator should be at 'last' leaf
00368     rval = tool.get_last_iterator( root, iter2 );CHECK_ERR( rval );
00369     for( ;; )
00370     {
00371         CHECK_EQUAL( iter.handle(), iter2.handle() );
00372         CHECK_EQUAL( iter.depth(), iter2.depth() );
00373         CHECK_REAL_EQUAL( iter.box_min()[0], iter2.box_min()[0], DBL_EPSILON );
00374         CHECK_REAL_EQUAL( iter.box_min()[1], iter2.box_min()[1], DBL_EPSILON );
00375         CHECK_REAL_EQUAL( iter.box_min()[2], iter2.box_min()[2], DBL_EPSILON );
00376         CHECK_REAL_EQUAL( iter.box_max()[0], iter2.box_max()[0], DBL_EPSILON );
00377         CHECK_REAL_EQUAL( iter.box_max()[1], iter2.box_max()[1], DBL_EPSILON );
00378         CHECK_REAL_EQUAL( iter.box_max()[2], iter2.box_max()[2], DBL_EPSILON );
00379 
00380         rval = iter2.back();
00381         if( MB_ENTITY_NOT_FOUND == rval ) break;CHECK_ERR( rval );
00382         rval = iter.back();CHECK_ERR( rval );
00383     }
00384 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines