MOAB: Mesh Oriented datABase  (version 5.3.1)
obb_tree_tool.cpp
Go to the documentation of this file.
00001 #define IS_BUILDING_MB
00002 #include "moab/Core.hpp"
00003 #include "moab/CartVect.hpp"
00004 #include "moab/OrientedBoxTreeTool.hpp"
00005 #include "moab/OrientedBox.hpp"
00006 #include "Internals.hpp"
00007 #include "moab/Range.hpp"
00008 
00009 #include <iostream>
00010 #include <iomanip>
00011 #include <cstdio>
00012 #include <limits>
00013 #include <cstdlib>
00014 #include <ctime>
00015 #if !defined( _MSC_VER ) && !defined( __MINGW32__ )
00016 #include <unistd.h>
00017 #include <sys/stat.h>
00018 #include <fcntl.h>
00019 #endif
00020 
00021 using namespace moab;
00022 
00023 std::string clock_to_string( clock_t t );
00024 std::string mem_to_string( unsigned long mem );
00025 
00026 const int MAX_TAG_VALUE    = 20;
00027 const char* const TAG_NAME = "OBB_ID";
00028 const char* const TREE_TAG = "OBB_ROOT";
00029 const char* root_tag       = TREE_TAG;
00030 
00031 ErrorCode get_root( Interface* moab, EntityHandle& root );
00032 EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings );
00033 void delete_existing_tree( Interface* interface );
00034 void print_stats( Interface* interface );
00035 void tag_triangles( Interface* interface );
00036 void tag_vertices( Interface* interface );
00037 void write_tree_blocks( Interface* interface, const char* file );
00038 
00039 static void usage( bool err = true )
00040 {
00041     std::ostream& s = err ? std::cerr : std::cout;
00042     s << "obb_tree_tool [-s|-S] [-d <int>] [-n <int>] <input file> <output file>" << std::endl
00043       << "obb_tree_tool [-h]" << std::endl;
00044     if( !err )
00045     {
00046         OrientedBoxTreeTool::Settings st;
00047         s << "Tool to build adaptive kd-Tree from triangles" << std::endl;
00048         s << "  -s        Use range-based sets for tree nodes" << std::endl
00049           << "  -S        Use vector-based sets for tree nodes" << std::endl
00050           << "  -d <int>  Specify maximum depth for tree. Default: " << st.max_depth << std::endl
00051           << "  -n <int>  Specify maximum entities per leaf. Default: " << st.max_leaf_entities << std::endl
00052           << "  -m <real> Specify worst split ratio. Default: " << st.worst_split_ratio << std::endl
00053           << "  -M <real> Specify best split ratio. Default: " << st.best_split_ratio << std::endl
00054           << "  -t        Tag triangles will tree cell number." << std::endl
00055           << "  -T        Write tree boxes to file." << std::endl
00056           << "  -N        Specify mesh tag containing tree root. Default: \"" << TREE_TAG << '"' << std::endl
00057           << std::endl;
00058     }
00059     exit( err );
00060 }
00061 
00062 #if defined( _MSC_VER ) || defined( __MINGW32__ )
00063 static void memory_use( unsigned long long& vsize, unsigned long long& rss )
00064 {
00065     vsize = rss = 0;
00066 }
00067 #else
00068 static void memory_use( unsigned long long& vsize, unsigned long long& rss )
00069 {
00070     char buffer[512];
00071     unsigned long lvsize;
00072     long lrss;
00073     int filp  = open( "/proc/self/stat", O_RDONLY );
00074     ssize_t r = read( filp, buffer, sizeof( buffer ) - 1 );
00075     close( filp );
00076     if( r < 0 ) r = 0;
00077     lvsize = lrss = 0;
00078     buffer[r]     = '\0';
00079     sscanf( buffer,
00080             "%*d %*s %*c "      // pid command state
00081             "%*d %*d "          // ppid pgrp
00082             "%*d %*d %*d "      // session tty_nr tpgid
00083             "%*u "              // flags
00084             "%*u %*u %*u %*u "  // minflt cminflt majflt cmajflt
00085             "%*u %*u %*d %*d "  // utime stime cutime cstime
00086             "%*d %*d %*d "      // priority nice (unused)
00087             "%*d %*u "          // itrealval starttime
00088             "%lu %ld",
00089             &lvsize, &lrss );
00090     rss   = lrss * getpagesize();
00091     vsize = lvsize;
00092 }
00093 #endif
00094 
00095 static int parseint( int& i, int argc, char* argv[] )
00096 {
00097     char* end;
00098     ++i;
00099     if( i == argc )
00100     {
00101         std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
00102         usage();
00103     }
00104 
00105     int result = strtol( argv[i], &end, 0 );
00106     if( result < 0 || *end )
00107     {
00108         std::cerr << "Expected positive integer following '" << argv[i - 1] << "'" << std::endl;
00109         usage();
00110     }
00111 
00112     return result;
00113 }
00114 
00115 static double parsedouble( int& i, int argc, char* argv[] )
00116 {
00117     char* end;
00118     ++i;
00119     if( i == argc )
00120     {
00121         std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
00122         usage();
00123     }
00124 
00125     double result = strtod( argv[i], &end );
00126     if( result < 0 || *end )
00127     {
00128         std::cerr << "Expected positive real number following '" << argv[i - 1] << "'" << std::endl;
00129         usage();
00130     }
00131 
00132     return result;
00133 }
00134 
00135 int main( int argc, char* argv[] )
00136 {
00137     const char* input_file  = 0;
00138     const char* output_file = 0;
00139     const char* tree_file   = 0;
00140     OrientedBoxTreeTool::Settings settings;
00141     bool tag_tris = false;
00142     clock_t load_time, build_time, stat_time, tag_time, write_time, block_time;
00143 
00144     for( int i = 1; i < argc; ++i )
00145     {
00146         if( argv[i][0] != '-' )
00147         {
00148             if( !input_file )
00149                 input_file = argv[i];
00150             else if( !output_file )
00151                 output_file = argv[i];
00152             else
00153                 usage();
00154             continue;
00155         }
00156 
00157         if( !argv[i][1] || argv[i][2] ) usage();
00158 
00159         switch( argv[i][1] )
00160         {
00161             case 's':
00162                 settings.set_options = MESHSET_SET;
00163                 break;
00164             case 'S':
00165                 settings.set_options = MESHSET_ORDERED;
00166                 break;
00167             case 'd':
00168                 settings.max_depth = parseint( i, argc, argv );
00169                 break;
00170             case 'n':
00171                 settings.max_leaf_entities = parseint( i, argc, argv );
00172                 break;
00173             case 'm':
00174                 settings.worst_split_ratio = parsedouble( i, argc, argv );
00175                 break;
00176             case 'M':
00177                 settings.best_split_ratio = parsedouble( i, argc, argv );
00178                 break;
00179             case 't':
00180                 tag_tris = true;
00181                 break;
00182             case 'T':
00183                 if( ++i == argc ) usage();
00184                 tree_file = argv[i];
00185                 break;
00186             case 'N':
00187                 if( ++i == argc ) usage();
00188                 root_tag = argv[i];
00189                 break;
00190             case 'h':
00191                 usage( false );
00192             default:
00193                 usage();
00194         }
00195     }
00196 
00197     if( !output_file ) usage();
00198 
00199     ErrorCode rval;
00200     Core moab_core;
00201     Interface* interface = &moab_core;
00202 
00203     load_time = clock();
00204     rval      = interface->load_mesh( input_file );
00205     if( MB_SUCCESS != rval )
00206     {
00207         std::cerr << "Error reading file: " << input_file << std::endl;
00208         exit( 2 );
00209     }
00210     load_time = clock() - load_time;
00211 
00212     delete_existing_tree( interface );
00213 
00214     std::cout << "Building tree..." << std::endl;
00215     build_time = clock();
00216     build_tree( interface, settings );
00217     build_time = clock() - build_time;
00218 
00219     std::cout << "Calculating stats..." << std::endl;
00220     print_stats( interface );
00221     stat_time = clock() - build_time;
00222 
00223     if( tag_tris )
00224     {
00225         std::cout << "Tagging tree..." << std::endl;
00226         tag_triangles( interface );
00227         tag_vertices( interface );
00228     }
00229     tag_time = clock() - stat_time;
00230 
00231     std::cout << "Writing file... ";
00232     std::cout.flush();
00233     rval = interface->write_mesh( output_file );
00234     if( MB_SUCCESS != rval )
00235     {
00236         std::cerr << "Error writing file: " << output_file << std::endl;
00237         exit( 3 );
00238     }
00239     write_time = clock() - tag_time;
00240     std::cout << "Wrote " << output_file << std::endl;
00241 
00242     if( tree_file )
00243     {
00244         std::cout << "Writing tree block rep...";
00245         std::cout.flush();
00246         write_tree_blocks( interface, tree_file );
00247         std::cout << "Wrote " << tree_file << std::endl;
00248     }
00249     block_time = clock() - write_time;
00250 
00251     std::cout << "Times:  "
00252               << "    Load"
00253               << "   Build"
00254               << "   Stats"
00255               << "   Write";
00256     if( tag_tris ) std::cout << "Tag Sets";
00257     if( tree_file ) std::cout << "Block   ";
00258     std::cout << std::endl;
00259 
00260     std::cout << "        " << std::setw( 8 ) << clock_to_string( load_time ) << std::setw( 8 )
00261               << clock_to_string( build_time ) << std::setw( 8 ) << clock_to_string( stat_time ) << std::setw( 8 )
00262               << clock_to_string( write_time );
00263     if( tag_tris ) std::cout << std::setw( 8 ) << clock_to_string( tag_time );
00264     if( tree_file ) std::cout << std::setw( 8 ) << clock_to_string( block_time );
00265     std::cout << std::endl;
00266 
00267     return 0;
00268 }
00269 
00270 ErrorCode get_root( Interface* moab, EntityHandle& root )
00271 {
00272     Tag tag;
00273     ErrorCode rval;
00274 
00275     rval = moab->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, tag );
00276     if( MB_SUCCESS != rval ) return rval;
00277 
00278     const EntityHandle mesh = 0;
00279     return moab->tag_get_data( tag, &mesh, 1, &root );
00280 }
00281 
00282 void delete_existing_tree( Interface* interface )
00283 {
00284     EntityHandle root;
00285     ErrorCode rval = get_root( interface, root );
00286     if( MB_SUCCESS == rval )
00287     {
00288         OrientedBoxTreeTool tool( interface );
00289         rval = tool.delete_tree( root );
00290         if( MB_SUCCESS != rval )
00291         {
00292             std::cerr << "Failed to destroy existing trees. Aborting" << std::endl;
00293             exit( 5 );
00294         }
00295     }
00296 }
00297 
00298 EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings )
00299 {
00300     ErrorCode rval;
00301     EntityHandle root = 0;
00302     Range triangles;
00303 
00304     rval = interface->get_entities_by_type( 0, MBTRI, triangles );
00305     if( MB_SUCCESS != rval || triangles.empty() )
00306     {
00307         std::cerr << "No triangles from which to build tree." << std::endl;
00308         exit( 4 );
00309     }
00310 
00311     OrientedBoxTreeTool tool( interface );
00312     rval = tool.build( triangles, root, &settings );
00313     if( MB_SUCCESS != rval || !root )
00314     {
00315         std::cerr << "Tree construction failed." << std::endl;
00316         exit( 4 );
00317     }
00318 
00319     // store tree root
00320     Tag roottag;
00321     rval = interface->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, roottag, MB_TAG_CREAT | MB_TAG_SPARSE );
00322     if( MB_SUCCESS != rval )
00323     {
00324         std::cout << "Failed to create root tag: \"" << root_tag << '"' << std::endl;
00325         exit( 2 );
00326     }
00327     const EntityHandle mesh = 0;
00328     rval                    = interface->tag_set_data( roottag, &mesh, 1, &root );
00329     if( MB_SUCCESS != rval )
00330     {
00331         std::cout << "Failed to set root tag: \"" << root_tag << '"' << std::endl;
00332         exit( 2 );
00333     }
00334 
00335     return root;
00336 }
00337 
00338 std::string clock_to_string( clock_t t )
00339 {
00340     char unit[5] = "s";
00341     char buffer[256];
00342     double dt = t / (double)CLOCKS_PER_SEC;
00343     // if (dt > 300) {
00344     //  dt /= 60;
00345     //  strcpy( unit, "min" );
00346     //}
00347     // if (dt > 300) {
00348     //  dt /= 60;
00349     //  strcpy( unit, "hr" );
00350     //}
00351     // if (dt > 100) {
00352     //  dt /= 24;
00353     //  strcpy( unit, "days" );
00354     //}
00355     sprintf( buffer, "%0.2f%s", dt, unit );
00356     return buffer;
00357 }
00358 
00359 std::string mem_to_string( unsigned long mem )
00360 {
00361     char unit[3] = "B";
00362     if( mem > 9 * 1024 )
00363     {
00364         mem = ( mem + 512 ) / 1024;
00365         strcpy( unit, "kB" );
00366     }
00367     if( mem > 9 * 1024 )
00368     {
00369         mem = ( mem + 512 ) / 1024;
00370         strcpy( unit, "MB" );
00371     }
00372     if( mem > 9 * 1024 )
00373     {
00374         mem = ( mem + 512 ) / 1024;
00375         strcpy( unit, "GB" );
00376     }
00377     char buffer[256];
00378     sprintf( buffer, "%lu %s", mem, unit );
00379     return buffer;
00380 }
00381 
00382 template < typename T >
00383 struct SimpleStat
00384 {
00385     T min, max, sum, sqr;
00386     size_t count;
00387     SimpleStat();
00388     void add( T value );
00389     double avg() const
00390     {
00391         return (double)sum / count;
00392     }
00393     double rms() const
00394     {
00395         return sqrt( (double)sqr / count );
00396     }
00397     double dev() const
00398     {
00399         return sqrt( ( count * (double)sqr - (double)sum * (double)sum ) / ( (double)count * ( count - 1 ) ) );
00400     }
00401 };
00402 
00403 template < typename T >
00404 SimpleStat< T >::SimpleStat()
00405     : min( std::numeric_limits< T >::max() ), max( std::numeric_limits< T >::min() ), sum( 0 ), sqr( 0 ), count( 0 )
00406 {
00407 }
00408 
00409 void print_stats( Interface* interface )
00410 {
00411     EntityHandle root;
00412     Range range;
00413     ErrorCode rval;
00414     rval = get_root( interface, root );
00415     if( MB_SUCCESS != rval )
00416     {
00417         std::cerr << "Internal error: Failed to retrieve root." << std::endl;
00418         exit( 5 );
00419     }
00420     OrientedBoxTreeTool tool( interface );
00421 
00422     Range tree_sets, triangles, verts;
00423     // interface->get_child_meshsets( root, tree_sets, 0 );
00424     interface->get_entities_by_type( 0, MBENTITYSET, tree_sets );
00425     tree_sets.erase( tree_sets.begin(), Range::lower_bound( tree_sets.begin(), tree_sets.end(), root ) );
00426     interface->get_entities_by_type( 0, MBTRI, triangles );
00427     interface->get_entities_by_type( 0, MBVERTEX, verts );
00428     triangles.merge( verts );
00429     tree_sets.insert( root );
00430     unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized,
00431         tri_used, tri_amortized;
00432     interface->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0,
00433                                      0, 0, &set_tag_used, &set_tag_amortized );
00434     interface->estimated_memory_use( triangles, &tri_used, &tri_amortized );
00435 
00436     int num_tri = 0;
00437     interface->get_number_entities_by_type( 0, MBTRI, num_tri );
00438 
00439     tool.stats( root, std::cout );
00440 
00441     unsigned long long real_rss, real_vsize;
00442     memory_use( real_vsize, real_rss );
00443 
00444     printf( "------------------------------------------------------------------\n" );
00445     printf( "\nmemory:           used  amortized\n" );
00446     printf( "            ---------- ----------\n" );
00447     printf( "triangles   %10s %10s\n", mem_to_string( tri_used ).c_str(), mem_to_string( tri_amortized ).c_str() );
00448     printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() );
00449     printf( "sets        %10s %10s\n", mem_to_string( set_store_used ).c_str(),
00450             mem_to_string( set_store_amortized ).c_str() );
00451     printf( "set tags    %10s %10s\n", mem_to_string( set_tag_used ).c_str(),
00452             mem_to_string( set_tag_amortized ).c_str() );
00453     printf( "total real  %10s %10s\n", mem_to_string( real_rss ).c_str(), mem_to_string( real_vsize ).c_str() );
00454     printf( "------------------------------------------------------------------\n" );
00455 }
00456 
00457 static int hash_handle( EntityHandle handle )
00458 {
00459     EntityID h = ID_FROM_HANDLE( handle );
00460     return (int)( ( h * 13 + 7 ) % MAX_TAG_VALUE ) + 1;
00461 }
00462 
00463 class TriTagger : public OrientedBoxTreeTool::Op
00464 {
00465   private:
00466     Interface* mMB;
00467     Tag mTag;
00468     std::vector< EntityHandle > mHandles;
00469     std::vector< int > mTagData;
00470 
00471   public:
00472     TriTagger( Tag tag, Interface* moab ) : mMB( moab ), mTag( tag ) {}
00473 
00474     ErrorCode visit( EntityHandle, int, bool& descent )
00475     {
00476         descent = true;
00477         return MB_SUCCESS;
00478     }
00479 
00480     ErrorCode leaf( EntityHandle node )
00481     {
00482         mHandles.clear();
00483         mMB->get_entities_by_handle( node, mHandles );
00484         mTagData.clear();
00485         mTagData.resize( mHandles.size(), hash_handle( node ) );
00486         mMB->tag_set_data( mTag, &mHandles[0], mHandles.size(), &mTagData[0] );
00487         return MB_SUCCESS;
00488     }
00489 };
00490 
00491 void tag_triangles( Interface* moab )
00492 {
00493     EntityHandle root;
00494     ErrorCode rval = get_root( moab, root );
00495     if( MB_SUCCESS != rval )
00496     {
00497         std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
00498         exit( 5 );
00499     }
00500 
00501     Tag tag;
00502     int zero = 0;
00503     moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00504     TriTagger op( tag, moab );
00505 
00506     OrientedBoxTreeTool tool( moab );
00507     rval = tool.preorder_traverse( root, op );
00508     if( MB_SUCCESS != rval )
00509     {
00510         std::cerr << "Internal error tagging triangles" << std::endl;
00511         exit( 5 );
00512     }
00513 }
00514 
00515 class VtxTagger : public OrientedBoxTreeTool::Op
00516 {
00517   private:
00518     Interface* mMB;
00519     Tag mTag;
00520     std::vector< EntityHandle > mHandles;
00521     std::vector< EntityHandle > mConn;
00522     std::vector< int > mTagData;
00523 
00524   public:
00525     VtxTagger( Tag tag, Interface* moab ) : mMB( moab ), mTag( tag ) {}
00526 
00527     ErrorCode visit( EntityHandle, int, bool& descent )
00528     {
00529         descent = true;
00530         return MB_SUCCESS;
00531     }
00532 
00533     ErrorCode leaf( EntityHandle node )
00534     {
00535         mHandles.clear();
00536         mMB->get_entities_by_handle( node, mHandles );
00537         mConn.clear();
00538         mMB->get_connectivity( &mHandles[0], mHandles.size(), mConn );
00539         mTagData.clear();
00540         mTagData.resize( mConn.size(), hash_handle( node ) );
00541         mMB->tag_set_data( mTag, &mConn[0], mConn.size(), &mTagData[0] );
00542         return MB_SUCCESS;
00543     }
00544 };
00545 
00546 void tag_vertices( Interface* moab )
00547 {
00548     EntityHandle root;
00549     ErrorCode rval = get_root( moab, root );
00550     if( MB_SUCCESS != rval )
00551     {
00552         std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
00553         exit( 5 );
00554     }
00555 
00556     Tag tag;
00557     int zero = 0;
00558     moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00559     VtxTagger op( tag, moab );
00560 
00561     OrientedBoxTreeTool tool( moab );
00562     rval = tool.preorder_traverse( root, op );
00563     if( MB_SUCCESS != rval )
00564     {
00565         std::cerr << "Internal error tagging vertices" << std::endl;
00566         exit( 5 );
00567     }
00568 }
00569 
00570 class LeafHexer : public OrientedBoxTreeTool::Op
00571 {
00572   private:
00573     OrientedBoxTreeTool* mTool;
00574     Interface* mOut;
00575     Tag mTag;
00576     std::vector< EntityHandle > mHandles;
00577     std::vector< EntityHandle > mConn;
00578     std::vector< int > mTagData;
00579 
00580   public:
00581     LeafHexer( OrientedBoxTreeTool* tool, Interface* mb2, Tag tag ) : mTool( tool ), mOut( mb2 ), mTag( tag ) {}
00582 
00583     ErrorCode visit( EntityHandle, int, bool& descent )
00584     {
00585         descent = true;
00586         return MB_SUCCESS;
00587     }
00588 
00589     ErrorCode leaf( EntityHandle node )
00590     {
00591         OrientedBox box;
00592         ErrorCode rval = mTool->box( node, box );
00593         if( MB_SUCCESS != rval ) return rval;
00594         EntityHandle h;
00595         rval = box.make_hex( h, mOut );
00596         if( MB_SUCCESS != rval ) return rval;
00597         int i = hash_handle( node );
00598         return mOut->tag_set_data( mTag, &h, 1, &i );
00599     }
00600 };
00601 
00602 void write_tree_blocks( Interface* interface, const char* file )
00603 {
00604     EntityHandle root;
00605     ErrorCode rval = get_root( interface, root );
00606     if( MB_SUCCESS != rval )
00607     {
00608         std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
00609         exit( 5 );
00610     }
00611 
00612     Core moab2;
00613     Tag tag;
00614     int zero = 0;
00615     moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
00616 
00617     OrientedBoxTreeTool tool( interface );
00618     LeafHexer op( &tool, &moab2, tag );
00619     rval = tool.preorder_traverse( root, op );
00620     if( MB_SUCCESS != rval )
00621     {
00622         std::cerr << "Internal error: failed to construct leaf hexes" << std::endl;
00623         exit( 5 );
00624     }
00625 
00626     moab2.write_mesh( file );
00627 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines