obb_tree_tool.cpp

Go to the documentation of this file.

#define IS_BUILDING_MB
#include "moab/Core.hpp"
#include "moab/CartVect.hpp"
#include "moab/OrientedBoxTreeTool.hpp"
#include "moab/OrientedBox.hpp"
#include "Internals.hpp"
#include "moab/Range.hpp"
 
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <limits>
#include <cstdlib>
#include <ctime>
#if !defined( _MSC_VER ) && !defined( __MINGW32__ )
#include <unistd.h>
#include <sys/stat.h>
#include <fcntl.h>
#endif
 
using namespace moab;
 
std::string clock_to_string( clock_t t );
std::string mem_to_string( unsigned long mem );
 
const int MAX_TAG_VALUE = 20;
const char* const TAG_NAME = "OBB_ID";
const char* const TREE_TAG = "OBB_ROOT";
const char* root_tag = TREE_TAG;
 
ErrorCode get_root( Interface* moab, EntityHandle& root );
EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings );
void delete_existing_tree( Interface* interface );
void print_stats( Interface* interface );
void tag_triangles( Interface* interface );
void tag_vertices( Interface* interface );
void write_tree_blocks( Interface* interface, const char* file );
 
static void usage( bool err = true )
{
 std::ostream& s = err ? std::cerr : std::cout;
 s << "obb_tree_tool [-s|-S] [-d <int>] [-n <int>] <input file> <output file>" << std::endl
 << "obb_tree_tool [-h]" << std::endl;
 if( !err )
 {
 OrientedBoxTreeTool::Settings st;
 s << "Tool to build adaptive kd-Tree from triangles" << std::endl;
 s << " -s Use range-based sets for tree nodes" << std::endl
 << " -S Use vector-based sets for tree nodes" << std::endl
 << " -d <int> Specify maximum depth for tree. Default: " << st.max_depth << std::endl
 << " -n <int> Specify maximum entities per leaf. Default: " << st.max_leaf_entities << std::endl
 << " -m <real> Specify worst split ratio. Default: " << st.worst_split_ratio << std::endl
 << " -M <real> Specify best split ratio. Default: " << st.best_split_ratio << std::endl
 << " -t Tag triangles will tree cell number." << std::endl
 << " -T Write tree boxes to file." << std::endl
 << " -N Specify mesh tag containing tree root. Default: \"" << TREE_TAG << '"' << std::endl
 << std::endl;
 }
 exit( err );
}
 
#if defined( _MSC_VER ) || defined( __MINGW32__ )
static void memory_use( unsigned long long& vsize, unsigned long long& rss )
{
 vsize = rss = 0;
}
#else
static void memory_use( unsigned long long& vsize, unsigned long long& rss )
{
 char buffer[512];
 unsigned long lvsize;
 long lrss;
 int filp = open( "/proc/self/stat", O_RDONLY );
 ssize_t r = read( filp, buffer, sizeof( buffer ) - 1 );
 close( filp );
 if( r < 0 ) r = 0;
 lvsize = lrss = 0;
 buffer[r] = '\0';
 sscanf( buffer,
 "%*d %*s %*c " // pid command state
 "%*d %*d " // ppid pgrp
 "%*d %*d %*d " // session tty_nr tpgid
 "%*u " // flags
 "%*u %*u %*u %*u " // minflt cminflt majflt cmajflt
 "%*u %*u %*d %*d " // utime stime cutime cstime
 "%*d %*d %*d " // priority nice (unused)
 "%*d %*u " // itrealval starttime
 "%lu %ld",
 &lvsize, &lrss );
 rss = lrss * getpagesize();
 vsize = lvsize;
}
#endif
 
static int parseint( int& i, int argc, char* argv[] )
{
 char* end;
 ++i;
 if( i == argc )
 {
 std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 int result = strtol( argv[i], &end, 0 );
 if( result < 0 || *end )
 {
 std::cerr << "Expected positive integer following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 return result;
}
 
static double parsedouble( int& i, int argc, char* argv[] )
{
 char* end;
 ++i;
 if( i == argc )
 {
 std::cerr << "Expected value following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 double result = strtod( argv[i], &end );
 if( result < 0 || *end )
 {
 std::cerr << "Expected positive real number following '" << argv[i - 1] << "'" << std::endl;
 usage();
 }
 
 return result;
}
 
int main( int argc, char* argv[] )
{
 const char* input_file = 0;
 const char* output_file = 0;
 const char* tree_file = 0;
 OrientedBoxTreeTool::Settings settings;
 bool tag_tris = false;
 clock_t load_time, build_time, stat_time, tag_time, write_time, block_time;
 
 for( int i = 1; i < argc; ++i )
 {
 if( argv[i][0] != '-' )
 {
 if( !input_file )
 input_file = argv[i];
 else if( !output_file )
 output_file = argv[i];
 else
 usage();
 continue;
 }
 
 if( !argv[i][1] || argv[i][2] ) usage();
 
 switch( argv[i][1] )
 {
 case 's':
 settings.set_options = MESHSET_SET;
 break;
 case 'S':
 settings.set_options = MESHSET_ORDERED;
 break;
 case 'd':
 settings.max_depth = parseint( i, argc, argv );
 break;
 case 'n':
 settings.max_leaf_entities = parseint( i, argc, argv );
 break;
 case 'm':
 settings.worst_split_ratio = parsedouble( i, argc, argv );
 break;
 case 'M':
 settings.best_split_ratio = parsedouble( i, argc, argv );
 break;
 case 't':
 tag_tris = true;
 break;
 case 'T':
 if( ++i == argc ) usage();
 tree_file = argv[i];
 break;
 case 'N':
 if( ++i == argc ) usage();
 root_tag = argv[i];
 break;
 case 'h':
 usage( false );
 break;
 default:
 usage();
 }
 }
 
 if( !output_file ) usage();
 
 ErrorCode rval;
 Core moab_core;
 Interface* interface = &moab_core;
 
 load_time = clock();
 rval = interface->load_mesh( input_file );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Error reading file: " << input_file << std::endl;
 exit( 2 );
 }
 load_time = clock() - load_time;
 
 delete_existing_tree( interface );
 
 std::cout << "Building tree..." << std::endl;
 build_time = clock();
 build_tree( interface, settings );
 build_time = clock() - build_time;
 
 std::cout << "Calculating stats..." << std::endl;
 print_stats( interface );
 stat_time = clock() - build_time;
 
 if( tag_tris )
 {
 std::cout << "Tagging tree..." << std::endl;
 tag_triangles( interface );
 tag_vertices( interface );
 }
 tag_time = clock() - stat_time;
 
 std::cout << "Writing file... ";
 std::cout.flush();
 rval = interface->write_mesh( output_file );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Error writing file: " << output_file << std::endl;
 exit( 3 );
 }
 write_time = clock() - tag_time;
 std::cout << "Wrote " << output_file << std::endl;
 
 if( tree_file )
 {
 std::cout << "Writing tree block rep...";
 std::cout.flush();
 write_tree_blocks( interface, tree_file );
 std::cout << "Wrote " << tree_file << std::endl;
 }
 block_time = clock() - write_time;
 
 std::cout << "Times: "
 << " Load"
 << " Build"
 << " Stats"
 << " Write";
 if( tag_tris ) std::cout << "Tag Sets";
 if( tree_file ) std::cout << "Block ";
 std::cout << std::endl;
 
 std::cout << " " << std::setw( 8 ) << clock_to_string( load_time ) << std::setw( 8 )
 << clock_to_string( build_time ) << std::setw( 8 ) << clock_to_string( stat_time ) << std::setw( 8 )
 << clock_to_string( write_time );
 if( tag_tris ) std::cout << std::setw( 8 ) << clock_to_string( tag_time );
 if( tree_file ) std::cout << std::setw( 8 ) << clock_to_string( block_time );
 std::cout << std::endl;
 
 return 0;
}
 
ErrorCode get_root( Interface* moab, EntityHandle& root )
{
 Tag tag;
 ErrorCode rval;
 
 rval = moab->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, tag );
 if( MB_SUCCESS != rval ) return rval;
 
 const EntityHandle mesh = 0;
 return moab->tag_get_data( tag, &mesh, 1, &root );
}
 
void delete_existing_tree( Interface* interface )
{
 EntityHandle root;
 ErrorCode rval = get_root( interface, root );
 if( MB_SUCCESS == rval )
 {
 OrientedBoxTreeTool tool( interface );
 rval = tool.delete_tree( root );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Failed to destroy existing trees. Aborting" << std::endl;
 exit( 5 );
 }
 }
}
 
EntityHandle build_tree( Interface* interface, OrientedBoxTreeTool::Settings settings )
{
 ErrorCode rval;
 EntityHandle root = 0;
 Range triangles;
 
 rval = interface->get_entities_by_type( 0, MBTRI, triangles );
 if( MB_SUCCESS != rval || triangles.empty() )
 {
 std::cerr << "No triangles from which to build tree." << std::endl;
 exit( 4 );
 }
 
 OrientedBoxTreeTool tool( interface );
 rval = tool.build( triangles, root, &settings );
 if( MB_SUCCESS != rval || !root )
 {
 std::cerr << "Tree construction failed." << std::endl;
 exit( 4 );
 }
 
 // store tree root
 Tag roottag;
 rval = interface->tag_get_handle( root_tag, 1, MB_TYPE_HANDLE, roottag, MB_TAG_CREAT | MB_TAG_SPARSE );
 if( MB_SUCCESS != rval )
 {
 std::cout << "Failed to create root tag: \"" << root_tag << '"' << std::endl;
 exit( 2 );
 }
 const EntityHandle mesh = 0;
 rval = interface->tag_set_data( roottag, &mesh, 1, &root );
 if( MB_SUCCESS != rval )
 {
 std::cout << "Failed to set root tag: \"" << root_tag << '"' << std::endl;
 exit( 2 );
 }
 
 return root;
}
 
std::string clock_to_string( clock_t t )
{
 char unit[5] = "s";
 char buffer[256];
 double dt = t / (double)CLOCKS_PER_SEC;
 // if (dt > 300) {
 // dt /= 60;
 // strcpy( unit, "min" );
 //}
 // if (dt > 300) {
 // dt /= 60;
 // strcpy( unit, "hr" );
 //}
 // if (dt > 100) {
 // dt /= 24;
 // strcpy( unit, "days" );
 //}
 sprintf( buffer, "%0.2f%s", dt, unit );
 return buffer;
}
 
std::string mem_to_string( unsigned long mem )
{
 char unit[3] = "B";
 if( mem > 9 * 1024 )
 {
 mem = ( mem + 512 ) / 1024;
 strcpy( unit, "kB" );
 }
 if( mem > 9 * 1024 )
 {
 mem = ( mem + 512 ) / 1024;
 strcpy( unit, "MB" );
 }
 if( mem > 9 * 1024 )
 {
 mem = ( mem + 512 ) / 1024;
 strcpy( unit, "GB" );
 }
 char buffer[256];
 sprintf( buffer, "%lu %s", mem, unit );
 return buffer;
}
 
template < typename T >
struct SimpleStat
{
 T min, max, sum, sqr;
 size_t count;
 SimpleStat();
 void add( T value );
 double avg() const
 {
 return (double)sum / count;
 }
 double rms() const
 {
 return sqrt( (double)sqr / count );
 }
 double dev() const
 {
 return sqrt( ( count * (double)sqr - (double)sum * (double)sum ) / ( (double)count * ( count - 1 ) ) );
 }
};
 
template < typename T >
SimpleStat< T >::SimpleStat()
 : min( std::numeric_limits< T >::max() ), max( std::numeric_limits< T >::min() ), sum( 0 ), sqr( 0 ), count( 0 )
{
}
 
void print_stats( Interface* interface )
{
 EntityHandle root;
 Range range;
 ErrorCode rval;
 rval = get_root( interface, root );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error: Failed to retrieve root." << std::endl;
 exit( 5 );
 }
 OrientedBoxTreeTool tool( interface );
 
 Range tree_sets, triangles, verts;
 // interface->get_child_meshsets( root, tree_sets, 0 );
 interface->get_entities_by_type( 0, MBENTITYSET, tree_sets );
 tree_sets.erase( tree_sets.begin(), Range::lower_bound( tree_sets.begin(), tree_sets.end(), root ) );
 interface->get_entities_by_type( 0, MBTRI, triangles );
 interface->get_entities_by_type( 0, MBVERTEX, verts );
 triangles.merge( verts );
 tree_sets.insert( root );
 unsigned long long set_used, set_amortized, set_store_used, set_store_amortized, set_tag_used, set_tag_amortized,
 tri_used, tri_amortized;
 interface->estimated_memory_use( tree_sets, &set_used, &set_amortized, &set_store_used, &set_store_amortized, 0, 0,
 0, 0, &set_tag_used, &set_tag_amortized );
 interface->estimated_memory_use( triangles, &tri_used, &tri_amortized );
 
 int num_tri = 0;
 interface->get_number_entities_by_type( 0, MBTRI, num_tri );
 
 tool.stats( root, std::cout );
 
 unsigned long long real_rss, real_vsize;
 memory_use( real_vsize, real_rss );
 
 printf( "------------------------------------------------------------------\n" );
 printf( "\nmemory: used amortized\n" );
 printf( " ---------- ----------\n" );
 printf( "triangles %10s %10s\n", mem_to_string( tri_used ).c_str(), mem_to_string( tri_amortized ).c_str() );
 printf( "sets (total)%10s %10s\n", mem_to_string( set_used ).c_str(), mem_to_string( set_amortized ).c_str() );
 printf( "sets %10s %10s\n", mem_to_string( set_store_used ).c_str(),
 mem_to_string( set_store_amortized ).c_str() );
 printf( "set tags %10s %10s\n", mem_to_string( set_tag_used ).c_str(),
 mem_to_string( set_tag_amortized ).c_str() );
 printf( "total real %10s %10s\n", mem_to_string( real_rss ).c_str(), mem_to_string( real_vsize ).c_str() );
 printf( "------------------------------------------------------------------\n" );
}
 
static int hash_handle( EntityHandle handle )
{
 EntityID h = ID_FROM_HANDLE( handle );
 return (int)( ( h * 13 + 7 ) % MAX_TAG_VALUE ) + 1;
}
 
class TriTagger : public OrientedBoxTreeTool::Op
{
 private:
 Interface* mMB;
 Tag mTag;
 std::vector< EntityHandle > mHandles;
 std::vector< int > mTagData;
 
 public:
 TriTagger( Tag tag, Interface* moab ) : mMB( moab ), mTag( tag ) {}
 
 ErrorCode visit( EntityHandle, int, bool& descent )
 {
 descent = true;
 return MB_SUCCESS;
 }
 
 ErrorCode leaf( EntityHandle node )
 {
 mHandles.clear();
 mMB->get_entities_by_handle( node, mHandles );
 mTagData.clear();
 mTagData.resize( mHandles.size(), hash_handle( node ) );
 mMB->tag_set_data( mTag, &mHandles[0], mHandles.size(), &mTagData[0] );
 return MB_SUCCESS;
 }
};
 
void tag_triangles( Interface* moab )
{
 EntityHandle root;
 ErrorCode rval = get_root( moab, root );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
 exit( 5 );
 }
 
 Tag tag;
 int zero = 0;
 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
 TriTagger op( tag, moab );
 
 OrientedBoxTreeTool tool( moab );
 rval = tool.preorder_traverse( root, op );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error tagging triangles" << std::endl;
 exit( 5 );
 }
}
 
class VtxTagger : public OrientedBoxTreeTool::Op
{
 private:
 Interface* mMB;
 Tag mTag;
 std::vector< EntityHandle > mHandles;
 std::vector< EntityHandle > mConn;
 std::vector< int > mTagData;
 
 public:
 VtxTagger( Tag tag, Interface* moab ) : mMB( moab ), mTag( tag ) {}
 
 ErrorCode visit( EntityHandle, int, bool& descent )
 {
 descent = true;
 return MB_SUCCESS;
 }
 
 ErrorCode leaf( EntityHandle node )
 {
 mHandles.clear();
 mMB->get_entities_by_handle( node, mHandles );
 mConn.clear();
 mMB->get_connectivity( &mHandles[0], mHandles.size(), mConn );
 mTagData.clear();
 mTagData.resize( mConn.size(), hash_handle( node ) );
 mMB->tag_set_data( mTag, &mConn[0], mConn.size(), &mTagData[0] );
 return MB_SUCCESS;
 }
};
 
void tag_vertices( Interface* moab )
{
 EntityHandle root;
 ErrorCode rval = get_root( moab, root );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
 exit( 5 );
 }
 
 Tag tag;
 int zero = 0;
 moab->tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
 VtxTagger op( tag, moab );
 
 OrientedBoxTreeTool tool( moab );
 rval = tool.preorder_traverse( root, op );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error tagging vertices" << std::endl;
 exit( 5 );
 }
}
 
class LeafHexer : public OrientedBoxTreeTool::Op
{
 private:
 OrientedBoxTreeTool* mTool;
 Interface* mOut;
 Tag mTag;
 std::vector< EntityHandle > mHandles;
 std::vector< EntityHandle > mConn;
 std::vector< int > mTagData;
 
 public:
 LeafHexer( OrientedBoxTreeTool* tool, Interface* mb2, Tag tag ) : mTool( tool ), mOut( mb2 ), mTag( tag ) {}
 
 ErrorCode visit( EntityHandle, int, bool& descent )
 {
 descent = true;
 return MB_SUCCESS;
 }
 
 ErrorCode leaf( EntityHandle node )
 {
 OrientedBox box;
 ErrorCode rval = mTool->box( node, box );
 if( MB_SUCCESS != rval ) return rval;
 EntityHandle h;
 rval = box.make_hex( h, mOut );
 if( MB_SUCCESS != rval ) return rval;
 int i = hash_handle( node );
 return mOut->tag_set_data( mTag, &h, 1, &i );
 }
};
 
void write_tree_blocks( Interface* interface, const char* file )
{
 EntityHandle root;
 ErrorCode rval = get_root( interface, root );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error: Failed to retrieve tree." << std::endl;
 exit( 5 );
 }
 
 Core moab2;
 Tag tag;
 int zero = 0;
 moab2.tag_get_handle( TAG_NAME, 1, MB_TYPE_INTEGER, tag, MB_TAG_DENSE | MB_TAG_CREAT, &zero );
 
 OrientedBoxTreeTool tool( interface );
 LeafHexer op( &tool, &moab2, tag );
 rval = tool.preorder_traverse( root, op );
 if( MB_SUCCESS != rval )
 {
 std::cerr << "Internal error: failed to construct leaf hexes" << std::endl;
 exit( 5 );
 }
 
 moab2.write_mesh( file );
}