MOAB: Mesh Oriented datABase  (version 5.2.1)
obb_analysis.cpp
Go to the documentation of this file.
00001 #include <iostream>
00002 #include <fstream>
00003 #include <sstream>
00004 
00005 #include "moab/Core.hpp"
00006 #include "moab/Interface.hpp"
00007 #include "moab/OrientedBoxTreeTool.hpp"
00008 #include "moab/CartVect.hpp"
00009 #include "moab/OrientedBox.hpp"
00010 #include "moab/GeomTopoTool.hpp"
00011 
00012 #include "cgm2moab.hpp"
00013 
00014 using namespace moab;
00015 
00016 class TriCounter : public OrientedBoxTreeTool::Op
00017 {
00018 
00019   public:
00020     int count;
00021     Interface* mbi;
00022     OrientedBoxTreeTool* tool;
00023     const CartVect& pt;
00024 
00025     TriCounter( Interface* mbi_p, OrientedBoxTreeTool* tool_p, const CartVect& p )
00026         : OrientedBoxTreeTool::Op(), count( 0 ), mbi( mbi_p ), tool( tool_p ), pt( p )
00027     {
00028     }
00029 
00030     virtual ErrorCode visit( EntityHandle node, int, bool& descend )
00031     {
00032         OrientedBox box;
00033         ErrorCode rval = tool->box( node, box );
00034 
00035         descend = box.contained( pt, 1e-6 );
00036 
00037         return rval;
00038     }
00039 
00040     virtual ErrorCode leaf( EntityHandle node )
00041     {
00042 
00043         int numtris;
00044         ErrorCode rval = tool->get_moab_instance()->get_number_entities_by_type( node, MBTRI, numtris );
00045         count += numtris;
00046         return rval;
00047     }
00048 };
00049 
00050 ErrorCode obbvis_create( GeomTopoTool& gtt, std::vector< int >& volumes, int grid, std::string& filename )
00051 {
00052     OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
00053 
00054     CartVect min, max;
00055     EntityHandle vol = gtt.entity_by_id( 3, volumes.front() );
00056     double middle[3];
00057     double axis1[3], axis2[3], axis3[3];
00058     double minPt[3], maxPt[3];
00059     ErrorCode rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
00060     // compute min and max verticies
00061     for( int i = 0; i < 3; i++ )
00062     {
00063         double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
00064         minPt[i]   = middle[i] - sum;
00065         maxPt[i]   = middle[i] + sum;
00066     }
00067     CHECKERR( gtt, rval );
00068 
00069     /* Compute an axis-aligned bounding box of all the requested volumes */
00070     for( std::vector< int >::iterator i = volumes.begin() + 1; i != volumes.end(); ++i )
00071     {
00072         CartVect i_min, i_max;
00073         vol  = gtt.entity_by_id( 3, *i );
00074         rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
00075         // compute min and max verticies
00076         for( int j = 0; j < 3; j++ )
00077         {
00078             double sum = fabs( axis1[j] ) + fabs( axis2[j] ) + fabs( axis3[j] );
00079             minPt[j]   = middle[j] - sum;
00080             maxPt[j]   = middle[j] + sum;
00081         }
00082 
00083         for( int j = 0; j < 3; ++j )
00084         {
00085             min[j] = std::min( min[j], i_min[j] );
00086             max[j] = std::max( max[j], i_max[j] );
00087         }
00088     }
00089 
00090     // These vectors could be repurposed to describe an OBB without changing the loops below
00091     CartVect center( middle );
00092     CartVect v1( axis1 );
00093     CartVect v2( axis2 );
00094     CartVect v3( axis3 );
00095 
00096     /* Compute the vertices of the visualization grid.  Calculation points are at the center
00097        of each cell in this grid, so make grid+1 vertices in each direction. */
00098     int numpoints = pow( (double)( grid + 1 ), 3 );
00099     double* pgrid = new double[numpoints * 3];
00100     int idx       = 0;
00101 
00102     for( int i = 0; i < numpoints * 3; ++i )
00103         pgrid[i] = 0.0;
00104 
00105     for( int i = 0; i <= grid; ++i )
00106     {
00107         CartVect x = -v1 + ( ( v1 * 2.0 ) * ( i / (double)grid ) );
00108 
00109         for( int j = 0; j <= grid; ++j )
00110         {
00111             CartVect y = -v2 + ( ( v2 * 2.0 ) * ( j / (double)grid ) );
00112 
00113             for( int k = 0; k <= grid; ++k )
00114             {
00115                 CartVect z = -v3 + ( ( v3 * 2.0 ) * ( k / (double)grid ) );
00116 
00117                 CartVect p = center + x + y + z;
00118                 for( int d = 0; d < 3; ++d )
00119                 {
00120                     pgrid[idx++] = p[d];
00121                 }
00122             }
00123         }
00124     }
00125 
00126     /* Create a new MOAB to use for output, and build the vertex grid */
00127     Core mb2;
00128     Range r;
00129     rval = mb2.create_vertices( pgrid, numpoints, r );
00130     CHECKERR( mb2, rval );
00131 
00132     Tag lttag;
00133     rval = mb2.tag_get_handle( "LEAFTRIS", sizeof( int ), MB_TYPE_INTEGER, lttag,
00134                                MB_TAG_EXCL | MB_TAG_CREAT | MB_TAG_BYTES | MB_TAG_DENSE, 0 );
00135     CHECKERR( mb2, rval );
00136 
00137     int row  = grid + 1;
00138     int side = row * row;
00139     EntityHandle connect[8];
00140     EntityHandle hex;
00141 
00142     // offset from grid corner to grid center
00143     CartVect grid_hex_center_offset = ( v1 + v2 + v3 ) * 2 * ( 1.0 / grid );
00144 
00145     // collect all the surfaces from the requested volumes to iterate over --
00146     // this prevents checking a shared surface more than once.
00147     Range surfs;
00148     for( std::vector< int >::iterator it = volumes.begin(); it != volumes.end(); ++it )
00149     {
00150 
00151         vol = gtt.entity_by_id( 3, *it );
00152         Range it_surfs;
00153         rval = gtt.get_moab_instance()->get_child_meshsets( vol, it_surfs );
00154         CHECKERR( gtt, rval );
00155         surfs.merge( it_surfs );
00156     }
00157     std::cout << "visualizing " << surfs.size() << " surfaces." << std::endl;
00158 
00159     /* Build hexes for any point with 1 or more leaftris */
00160     for( int i = 0; i < grid; ++i )
00161     {
00162         for( int j = 0; j < grid; ++j )
00163         {
00164             for( int k = 0; k < grid; ++k )
00165             {
00166 
00167                 idx = ( side * k ) + ( row * j ) + i;
00168                 assert( idx + 1 + row + side > numpoints - 1 );
00169 
00170                 CartVect loc = CartVect( ( pgrid + ( idx * 3 ) ) ) + grid_hex_center_offset;
00171                 TriCounter tc( gtt.get_moab_instance(), &obbtool, loc );
00172 
00173                 for( Range::iterator it = surfs.begin(); it != surfs.end(); ++it )
00174                 {
00175 
00176                     EntityHandle surf_tree;
00177                     rval = gtt.get_root( *it, surf_tree );
00178                     CHECKERR( gtt, rval );
00179 
00180                     rval = obbtool.preorder_traverse( surf_tree, tc );
00181                     CHECKERR( gtt, rval );
00182                 }
00183 
00184                 if( tc.count == 0 ) continue;
00185 
00186                 connect[0] = r[idx];
00187                 connect[1] = r[idx + 1];
00188                 connect[2] = r[idx + 1 + row];
00189                 connect[3] = r[idx + row];
00190                 connect[4] = r[idx + side];
00191                 connect[5] = r[idx + 1 + side];
00192                 connect[6] = r[idx + 1 + row + side];
00193                 connect[7] = r[idx + row + side];
00194 
00195                 rval = mb2.create_element( MBHEX, connect, 8, hex );
00196                 CHECKERR( mb2, rval );
00197 
00198                 rval = mb2.tag_set_data( lttag, &hex, 1, &( tc.count ) );
00199                 CHECKERR( mb2, rval );
00200             }
00201         }
00202     }
00203 
00204     if( verbose ) { std::cout << "Writing " << filename << std::endl; }
00205     rval = mb2.write_file( filename.c_str() );
00206     CHECKERR( mb2, rval );
00207 
00208     return rval;
00209 }
00210 
00211 // stats code borrowed from OrientedBoxTreeTool.cpp
00212 static inline double std_dev( double sqr, double sum, double count )
00213 {
00214     sum /= count;
00215     sqr /= count;
00216     return sqrt( sqr - sum * sum );
00217 }
00218 
00219 class TriStats : public OrientedBoxTreeTool::Op
00220 {
00221 
00222   public:
00223     unsigned min, max, sum, leaves;
00224     double sqr;
00225 
00226     const static unsigned ten_buckets_max = 5;
00227     unsigned ten_buckets[ten_buckets_max];
00228     double ten_buckets_vol[ten_buckets_max];
00229 
00230     Interface* mbi;
00231     OrientedBoxTreeTool* tool;
00232 
00233     double tot_vol;
00234 
00235     TriStats( Interface* mbi_p, OrientedBoxTreeTool* tool_p, EntityHandle root )
00236         : OrientedBoxTreeTool::Op(), sum( 0 ), leaves( 0 ), sqr( 0 ), mbi( mbi_p ), tool( tool_p )
00237     {
00238         min = std::numeric_limits< unsigned >::max();
00239         max = std::numeric_limits< unsigned >::min();
00240 
00241         for( unsigned i = 0; i < ten_buckets_max; ++i )
00242         {
00243             ten_buckets[i]     = 0;
00244             ten_buckets_vol[i] = 0.;
00245         }
00246 
00247         OrientedBox box;
00248         ErrorCode rval = tool->box( root, box );
00249         CHECKERR( mbi, rval );
00250         tot_vol = box.volume();
00251     }
00252 
00253     virtual ErrorCode visit( EntityHandle, int, bool& descend )
00254     {
00255 
00256         descend = true;
00257         return MB_SUCCESS;
00258     }
00259 
00260     virtual ErrorCode leaf( EntityHandle node )
00261     {
00262 
00263         Range tris;
00264         ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
00265         unsigned count = tris.size();
00266 
00267         sum += count;
00268         sqr += ( count * count );
00269         if( min > count ) min = count;
00270         if( max < count ) max = count;
00271 
00272         for( unsigned i = 0; i < ten_buckets_max; ++i )
00273         {
00274             if( count > std::pow( (double)10, (int)( i + 1 ) ) )
00275             {
00276                 ten_buckets[i] += 1;
00277                 OrientedBox box;
00278                 rval = tool->box( node, box );
00279                 CHECKERR( mbi, rval );
00280                 ten_buckets_vol[i] += box.volume();
00281             }
00282         }
00283 
00284         leaves++;
00285 
00286         return rval;
00287     }
00288 
00289     std::string commafy( int num )
00290     {
00291         std::stringstream str;
00292         str << num;
00293         std::string s = str.str();
00294 
00295         int n = s.size();
00296         for( int i = n - 3; i >= 1; i -= 3 )
00297         {
00298             s.insert( i, 1, ',' );
00299             n++;
00300         }
00301 
00302         return s;
00303     }
00304 
00305     void write_results( std::ostream& out )
00306     {
00307 
00308         out << commafy( sum ) << " triangles in " << commafy( leaves ) << " leaves." << std::endl;
00309 
00310         double avg    = sum / (double)leaves;
00311         double stddev = std_dev( sqr, sum, leaves );
00312 
00313         out << "Tris per leaf: Min " << min << ", Max " << max << ", avg " << avg << ", stddev " << stddev << std::endl;
00314 
00315         for( unsigned i = 0; i < ten_buckets_max; ++i )
00316         {
00317             if( ten_buckets[i] )
00318             {
00319                 out << "Leaves exceeding " << std::pow( (double)10, (int)( i + 1 ) )
00320                     << " triangles: " << ten_buckets[i];
00321 
00322                 double frac_total_vol = ten_buckets_vol[i] / tot_vol;
00323                 double avg_ftv        = frac_total_vol / ten_buckets[i];
00324 
00325                 out << " (avg " << avg_ftv * 100.0 << "% of OBB volume)" << std::endl;
00326             }
00327         }
00328     }
00329 };
00330 
00331 /* no longer supported
00332 static std::string make_property_string( DagMC& dag, EntityHandle eh, std::vector<std::string>
00333 &properties )
00334 {
00335   ErrorCode ret;
00336   std::string propstring;
00337   for( std::vector<std::string>::iterator p = properties.begin();
00338     p != properties.end(); ++p )
00339   {
00340     if( dag.has_prop( eh, *p ) ){
00341       std::vector< std::string> vals;
00342       ret = dag.prop_values( eh, *p, vals );
00343       CHECKERR(dag,ret);
00344       propstring += *p;
00345       if( vals.size() == 1 ){
00346         propstring += "=";
00347         propstring += vals[0];
00348       }
00349       else if( vals.size() > 1 ){
00350         // this property has multiple values, list within brackets
00351         propstring += "=[";
00352         for( std::vector<std::string>::iterator i = vals.begin();
00353              i != vals.end(); ++i )
00354         {
00355             propstring += *i;
00356             propstring += ",";
00357         }
00358         // replace the last trailing comma with a close braket
00359         propstring[ propstring.length()-1 ] = ']';
00360       }
00361       propstring += ", ";
00362     }
00363   }
00364   if( propstring.length() ){
00365     propstring.resize( propstring.length() - 2 ); // drop trailing comma
00366   }
00367   return propstring;
00368 }
00369 */
00370 
00371 ErrorCode obbstat_write( GeomTopoTool& gtt, std::vector< int >& volumes, std::vector< std::string >& properties,
00372                          std::ostream& out )
00373 {
00374 
00375     ErrorCode ret                = MB_SUCCESS;
00376     OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
00377 
00378     // can assume that volume numbers are valid.
00379     for( std::vector< int >::iterator i = volumes.begin(); i != volumes.end(); ++i )
00380     {
00381         EntityHandle vol_root;
00382         EntityHandle vol = gtt.entity_by_id( 3, *i );
00383         CHECKERR( gtt, ret );
00384 
00385         if( vol == 0 )
00386         {
00387             std::cerr << "ERROR: volume " << *i << " has no entity." << std::endl;
00388             continue;
00389         }
00390 
00391         ret = gtt.get_root( vol, vol_root );
00392         CHECKERR( gtt, ret );
00393 
00394         out << "\nVolume " << *i << " " << std::flush;
00395 
00396         if( gtt.is_implicit_complement( vol ) ) out << "(implicit complement) ";
00397         out << std::endl;
00398 
00399         // get all surfaces in volume
00400         Range surfs;
00401         ret = gtt.get_moab_instance()->get_child_meshsets( vol, surfs );
00402         CHECKERR( gtt, ret );
00403 
00404         out << "   with " << surfs.size() << " surfaces" << std::endl;
00405 
00406         TriStats ts( gtt.get_moab_instance(), &obbtool, vol_root );
00407         ret = obbtool.preorder_traverse( vol_root, ts );
00408         CHECKERR( gtt, ret );
00409         ts.write_results( out );
00410 
00411         if( verbose )
00412         {
00413             out << "Surface list: " << std::flush;
00414             for( Range::iterator j = surfs.begin(); j != surfs.end(); ++j )
00415             {
00416                 out << gtt.global_id( *j );
00417                 if( j + 1 != surfs.end() ) out << ",";
00418             }
00419             out << std::endl;
00420             ret = obbtool.stats( vol_root, out );
00421             CHECKERR( gtt, ret );
00422         }
00423 
00424         out << "\n    ------------ " << std::endl;
00425     }
00426 
00427     return ret;
00428 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines