obb_analysis.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <sstream>
 
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/OrientedBoxTreeTool.hpp"
#include "moab/CartVect.hpp"
#include "moab/OrientedBox.hpp"
#include "moab/GeomTopoTool.hpp"
 
#include "cgm2moab.hpp"
 
using namespace moab;
 
class TriCounter : public OrientedBoxTreeTool::Op
{
 
  public:
    int count;
    Interface* mbi;
    OrientedBoxTreeTool* tool;
    const CartVect& pt;
 
    TriCounter( Interface* mbi_p, OrientedBoxTreeTool* tool_p, const CartVect& p )
        : OrientedBoxTreeTool::Op(), count( 0 ), mbi( mbi_p ), tool( tool_p ), pt( p )
    {
    }
 
    virtual ErrorCode visit( EntityHandle node, int, bool& descend )
    {
        OrientedBox box;
        ErrorCode rval = tool->box( node, box );
 
        descend = box.contained( pt, 1e-6 );
 
        return rval;
    }
 
    virtual ErrorCode leaf( EntityHandle node )
    {
 
        int numtris;
        ErrorCode rval = tool->get_moab_instance()->get_number_entities_by_type( node, MBTRI, numtris );
        count += numtris;
        return rval;
    }
};
 
ErrorCode obbvis_create( GeomTopoTool& gtt, std::vector< int >& volumes, int grid, std::string& filename )
{
    OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
 
    CartVect min, max;
    EntityHandle vol = gtt.entity_by_id( 3, volumes.front() );
    double middle[3];
    double axis1[3], axis2[3], axis3[3];
    double minPt[3], maxPt[3];
    ErrorCode rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
    // compute min and max verticies
    for( int i = 0; i < 3; i++ )
    {
        double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
        minPt[i]   = middle[i] - sum;
        maxPt[i]   = middle[i] + sum;
    }
    CHECKERR( gtt, rval );
 
    /* Compute an axis-aligned bounding box of all the requested volumes */
    for( std::vector< int >::iterator i = volumes.begin() + 1; i != volumes.end(); ++i )
    {
        CartVect i_min, i_max;
        vol  = gtt.entity_by_id( 3, *i );
        rval = gtt.get_obb( vol, middle, axis1, axis2, axis3 );
        // compute min and max verticies
        for( int j = 0; j < 3; j++ )
        {
            double sum = fabs( axis1[j] ) + fabs( axis2[j] ) + fabs( axis3[j] );
            minPt[j]   = middle[j] - sum;
            maxPt[j]   = middle[j] + sum;
        }
 
        for( int j = 0; j < 3; ++j )
        {
            min[j] = std::min( min[j], i_min[j] );
            max[j] = std::max( max[j], i_max[j] );
        }
    }
 
    // These vectors could be repurposed to describe an OBB without changing the loops below
    CartVect center( middle );
    CartVect v1( axis1 );
    CartVect v2( axis2 );
    CartVect v3( axis3 );
 
    /* Compute the vertices of the visualization grid.  Calculation points are at the center
       of each cell in this grid, so make grid+1 vertices in each direction. */
    int numpoints = pow( (double)( grid + 1 ), 3 );
    double* pgrid = new double[numpoints * 3];
    int idx       = 0;
 
    for( int i = 0; i < numpoints * 3; ++i )
        pgrid[i] = 0.0;
 
    for( int i = 0; i <= grid; ++i )
    {
        CartVect x = -v1 + ( ( v1 * 2.0 ) * ( i / (double)grid ) );
 
        for( int j = 0; j <= grid; ++j )
        {
            CartVect y = -v2 + ( ( v2 * 2.0 ) * ( j / (double)grid ) );
 
            for( int k = 0; k <= grid; ++k )
            {
                CartVect z = -v3 + ( ( v3 * 2.0 ) * ( k / (double)grid ) );
 
                CartVect p = center + x + y + z;
                for( int d = 0; d < 3; ++d )
                {
                    pgrid[idx++] = p[d];
                }
            }
        }
    }
 
    /* Create a new MOAB to use for output, and build the vertex grid */
    Core mb2;
    Range r;
    rval = mb2.create_vertices( pgrid, numpoints, r );
    CHECKERR( mb2, rval );
 
    Tag lttag;
    rval = mb2.tag_get_handle( "LEAFTRIS", sizeof( int ), MB_TYPE_INTEGER, lttag,
                               MB_TAG_EXCL | MB_TAG_CREAT | MB_TAG_BYTES | MB_TAG_DENSE, 0 );
    CHECKERR( mb2, rval );
 
    int row  = grid + 1;
    int side = row * row;
    EntityHandle connect[8];
    EntityHandle hex;
 
    // offset from grid corner to grid center
    CartVect grid_hex_center_offset = ( v1 + v2 + v3 ) * 2 * ( 1.0 / grid );
 
    // collect all the surfaces from the requested volumes to iterate over --
    // this prevents checking a shared surface more than once.
    Range surfs;
    for( std::vector< int >::iterator it = volumes.begin(); it != volumes.end(); ++it )
    {
 
        vol = gtt.entity_by_id( 3, *it );
        Range it_surfs;
        rval = gtt.get_moab_instance()->get_child_meshsets( vol, it_surfs );
        CHECKERR( gtt, rval );
        surfs.merge( it_surfs );
    }
    std::cout << "visualizing " << surfs.size() << " surfaces." << std::endl;
 
    /* Build hexes for any point with 1 or more leaftris */
    for( int i = 0; i < grid; ++i )
    {
        for( int j = 0; j < grid; ++j )
        {
            for( int k = 0; k < grid; ++k )
            {
 
                idx = ( side * k ) + ( row * j ) + i;
                assert( idx + 1 + row + side > numpoints - 1 );
 
                CartVect loc = CartVect( ( pgrid + ( idx * 3 ) ) ) + grid_hex_center_offset;
                TriCounter tc( gtt.get_moab_instance(), &obbtool, loc );
 
                for( Range::iterator it = surfs.begin(); it != surfs.end(); ++it )
                {
 
                    EntityHandle surf_tree;
                    rval = gtt.get_root( *it, surf_tree );
                    CHECKERR( gtt, rval );
 
                    rval = obbtool.preorder_traverse( surf_tree, tc );
                    CHECKERR( gtt, rval );
                }
 
                if( tc.count == 0 ) continue;
 
                connect[0] = r[idx];
                connect[1] = r[idx + 1];
                connect[2] = r[idx + 1 + row];
                connect[3] = r[idx + row];
                connect[4] = r[idx + side];
                connect[5] = r[idx + 1 + side];
                connect[6] = r[idx + 1 + row + side];
                connect[7] = r[idx + row + side];
 
                rval = mb2.create_element( MBHEX, connect, 8, hex );
                CHECKERR( mb2, rval );
 
                rval = mb2.tag_set_data( lttag, &hex, 1, &( tc.count ) );
                CHECKERR( mb2, rval );
            }
        }
    }
 
    if( verbose )
    {
        std::cout << "Writing " << filename << std::endl;
    }
    rval = mb2.write_file( filename.c_str() );
    CHECKERR( mb2, rval );
 
    return rval;
}
 
// stats code borrowed from OrientedBoxTreeTool.cpp
static inline double std_dev( double sqr, double sum, double count )
{
    sum /= count;
    sqr /= count;
    return sqrt( sqr - sum * sum );
}
 
class TriStats : public OrientedBoxTreeTool::Op
{
 
  public:
    unsigned min, max, sum, leaves;
    double sqr;
 
    const static unsigned ten_buckets_max = 5;
    unsigned ten_buckets[ten_buckets_max];
    double ten_buckets_vol[ten_buckets_max];
 
    Interface* mbi;
    OrientedBoxTreeTool* tool;
 
    double tot_vol;
 
    TriStats( Interface* mbi_p, OrientedBoxTreeTool* tool_p, EntityHandle root )
        : OrientedBoxTreeTool::Op(), sum( 0 ), leaves( 0 ), sqr( 0 ), mbi( mbi_p ), tool( tool_p )
    {
        min = std::numeric_limits< unsigned >::max();
        max = std::numeric_limits< unsigned >::min();
 
        for( unsigned i = 0; i < ten_buckets_max; ++i )
        {
            ten_buckets[i]     = 0;
            ten_buckets_vol[i] = 0.;
        }
 
        OrientedBox box;
        ErrorCode rval = tool->box( root, box );
        CHECKERR( mbi, rval );
        tot_vol = box.volume();
    }
 
    virtual ErrorCode visit( EntityHandle, int, bool& descend )
    {
 
        descend = true;
        return MB_SUCCESS;
    }
 
    virtual ErrorCode leaf( EntityHandle node )
    {
 
        Range tris;
        ErrorCode rval = tool->get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
        unsigned count = tris.size();
 
        sum += count;
        sqr += ( count * count );
        if( min > count ) min = count;
        if( max < count ) max = count;
 
        for( unsigned i = 0; i < ten_buckets_max; ++i )
        {
            if( count > std::pow( (double)10, (int)( i + 1 ) ) )
            {
                ten_buckets[i] += 1;
                OrientedBox box;
                rval = tool->box( node, box );
                CHECKERR( mbi, rval );
                ten_buckets_vol[i] += box.volume();
            }
        }
 
        leaves++;
 
        return rval;
    }
 
    std::string commafy( int num )
    {
        std::stringstream str;
        str << num;
        std::string s = str.str();
 
        int n = s.size();
        for( int i = n - 3; i >= 1; i -= 3 )
        {
            s.insert( i, 1, ',' );
            n++;
        }
 
        return s;
    }
 
    void write_results( std::ostream& out )
    {
 
        out << commafy( sum ) << " triangles in " << commafy( leaves ) << " leaves." << std::endl;
 
        double avg    = sum / (double)leaves;
        double stddev = std_dev( sqr, sum, leaves );
 
        out << "Tris per leaf: Min " << min << ", Max " << max << ", avg " << avg << ", stddev " << stddev << std::endl;
 
        for( unsigned i = 0; i < ten_buckets_max; ++i )
        {
            if( ten_buckets[i] )
            {
                out << "Leaves exceeding " << std::pow( (double)10, (int)( i + 1 ) )
                    << " triangles: " << ten_buckets[i];
 
                double frac_total_vol = ten_buckets_vol[i] / tot_vol;
                double avg_ftv        = frac_total_vol / ten_buckets[i];
 
                out << " (avg " << avg_ftv * 100.0 << "% of OBB volume)" << std::endl;
            }
        }
    }
};
 
/* no longer supported
static std::string make_property_string( DagMC& dag, EntityHandle eh, std::vector<std::string>
&properties )
{
  ErrorCode ret;
  std::string propstring;
  for( std::vector<std::string>::iterator p = properties.begin();
    p != properties.end(); ++p )
  {
    if( dag.has_prop( eh, *p ) ){
      std::vector< std::string> vals;
      ret = dag.prop_values( eh, *p, vals );
      CHECKERR(dag,ret);
      propstring += *p;
      if( vals.size() == 1 ){
        propstring += "=";
        propstring += vals[0];
      }
      else if( vals.size() > 1 ){
        // this property has multiple values, list within brackets
        propstring += "=[";
        for( std::vector<std::string>::iterator i = vals.begin();
             i != vals.end(); ++i )
        {
            propstring += *i;
            propstring += ",";
        }
        // replace the last trailing comma with a close braket
        propstring[ propstring.length()-1 ] = ']';
      }
      propstring += ", ";
    }
  }
  if( propstring.length() ){
    propstring.resize( propstring.length() - 2 ); // drop trailing comma
  }
  return propstring;
}
*/
 
ErrorCode obbstat_write( GeomTopoTool& gtt,
                         std::vector< int >& volumes,
                         std::vector< std::string >& /*properties*/,
                         std::ostream& out )
{
 
    ErrorCode ret                = MB_SUCCESS;
    OrientedBoxTreeTool& obbtool = *gtt.obb_tree();
 
    // can assume that volume numbers are valid.
    for( std::vector< int >::iterator i = volumes.begin(); i != volumes.end(); ++i )
    {
        EntityHandle vol_root;
        EntityHandle vol = gtt.entity_by_id( 3, *i );
        CHECKERR( gtt, ret );
 
        if( vol == 0 )
        {
            std::cerr << "ERROR: volume " << *i << " has no entity." << std::endl;
            continue;
        }
 
        ret = gtt.get_root( vol, vol_root );
        CHECKERR( gtt, ret );
 
        out << "\nVolume " << *i << " " << std::flush;
 
        if( gtt.is_implicit_complement( vol ) ) out << "(implicit complement) ";
        out << std::endl;
 
        // get all surfaces in volume
        Range surfs;
        ret = gtt.get_moab_instance()->get_child_meshsets( vol, surfs );
        CHECKERR( gtt, ret );
 
        out << "   with " << surfs.size() << " surfaces" << std::endl;
 
        TriStats ts( gtt.get_moab_instance(), &obbtool, vol_root );
        ret = obbtool.preorder_traverse( vol_root, ts );
        CHECKERR( gtt, ret );
        ts.write_results( out );
 
        if( verbose )
        {
            out << "Surface list: " << std::flush;
            for( Range::iterator j = surfs.begin(); j != surfs.end(); ++j )
            {
                out << gtt.global_id( *j );
                if( j + 1 != surfs.end() ) out << ",";
            }
            out << std::endl;
            ret = obbtool.stats( vol_root, out );
            CHECKERR( gtt, ret );
        }
 
        out << "\n    ------------ " << std::endl;
    }
 
    return ret;
}