MOAB: Mesh Oriented datABase  (version 5.2.1)
surfplot.cpp
Go to the documentation of this file.
00001 #include "moab/Core.hpp"
00002 #include "moab/Range.hpp"
00003 #include "MBTagConventions.hpp"
00004 #include <iostream>
00005 #include <fstream>
00006 #include <limits>
00007 #include <cstdlib>
00008 #include <math.h>
00009 
00010 /* Exit values */
00011 #define USAGE_ERROR       1
00012 #define READ_ERROR        2
00013 #define WRITE_ERROR       3
00014 #define SURFACE_NOT_FOUND 4
00015 #define OTHER_ERROR       5
00016 
00017 static void usage_error( const char* name )
00018 {
00019     std::cerr << "Usage: " << name << " [-g|-p] <Surface_ID> <input_file>" << std::endl
00020               << "\t-g           -  Write GNU Plot data file (default)." << std::endl
00021               << "\t-p           -  Write encapsulated postscript file" << std::endl
00022               << "\t-s           -  Write an SVG file" << std::endl
00023               << "\t<Surface_ID> -  ID of surface containing mesh to export (0 for entire file)." << std::endl
00024               << "\t<input_file> -  Mesh file to read." << std::endl
00025               << std::endl
00026               << "  This utility plots the mesh of a single geometric surface "
00027               << "projected to a plane.  The output file is written to stdout." << std::endl;
00028 
00029     exit( USAGE_ERROR );
00030 }
00031 
00032 struct CartVect3D
00033 {
00034 
00035     double x, y, z;
00036 
00037     CartVect3D() : x( 0.0 ), y( 0.0 ), z( 0.0 ) {}
00038 
00039     CartVect3D( double x_, double y_, double z_ ) : x( x_ ), y( y_ ), z( z_ ) {}
00040 
00041     CartVect3D& operator+=( const CartVect3D& o )
00042     {
00043         x += o.x;
00044         y += o.y;
00045         z += o.z;
00046         return *this;
00047     }
00048 
00049     CartVect3D& operator-=( const CartVect3D& o )
00050     {
00051         x -= o.x;
00052         y -= o.y;
00053         z -= o.z;
00054         return *this;
00055     }
00056 
00057     CartVect3D& operator*=( const CartVect3D& );
00058 
00059     CartVect3D& operator+=( double v )
00060     {
00061         x += v;
00062         y += v;
00063         z += v;
00064         return *this;
00065     }
00066 
00067     CartVect3D& operator-=( double v )
00068     {
00069         x -= v;
00070         y -= v;
00071         z -= v;
00072         return *this;
00073     }
00074 
00075     CartVect3D& operator*=( double v )
00076     {
00077         x *= v;
00078         y *= v;
00079         z *= v;
00080         return *this;
00081     }
00082 
00083     CartVect3D& operator/=( double v )
00084     {
00085         x /= v;
00086         y /= v;
00087         z /= v;
00088         return *this;
00089     }
00090 
00091     double len() const
00092     {
00093         return sqrt( x * x + y * y + z * z );
00094     }
00095 };
00096 
00097 // static CartVect3D operator-( const CartVect3D &a )
00098 //  { return CartVect3D(-a.z, -a.y, -a.z);  }
00099 
00100 // static CartVect3D operator+( const CartVect3D &a, const CartVect3D &b )
00101 //  { return CartVect3D(a.x+b.x, a.y+b.y, a.z+b.z); }
00102 
00103 static CartVect3D operator-( const CartVect3D& a, const CartVect3D& b )
00104 {
00105     return CartVect3D( a.x - b.x, a.y - b.y, a.z - b.z );
00106 }
00107 
00108 static double operator%( const CartVect3D& a, const CartVect3D& b )
00109 {
00110     return a.x * b.x + a.y * b.y + a.z * b.z;
00111 }
00112 
00113 static CartVect3D operator*( const CartVect3D& a, const CartVect3D& b )
00114 {
00115     CartVect3D result;
00116     result.x = a.y * b.z - a.z * b.y;
00117     result.y = a.z * b.x - a.x * b.z;
00118     result.z = a.x * b.y - a.y * b.x;
00119     return result;
00120 }
00121 
00122 CartVect3D& CartVect3D::operator*=( const CartVect3D& o )
00123 {
00124     *this = *this * o;
00125     return *this;
00126 }
00127 
00128 static void find_rotation( CartVect3D plane_normal, double matrix[3][3] )
00129 {
00130     // normalize
00131     plane_normal /= plane_normal.len();
00132     if( fabs( plane_normal.x ) < 0.1 ) plane_normal.x = 0.0;
00133     if( fabs( plane_normal.y ) < 0.1 ) plane_normal.y = 0.0;
00134     if( fabs( plane_normal.z ) < 0.1 ) plane_normal.z = 0.0;
00135 
00136     // calculate vector to rotate about
00137     const CartVect3D Z( 0, 0, 1 );
00138     CartVect3D vector = plane_normal * Z;
00139     const double len  = vector.len();
00140 
00141     // If vector is zero, no rotation
00142     if( len < 1e-2 )
00143     {
00144         matrix[0][0] = matrix[1][1] = matrix[2][2] = 1.0;
00145         matrix[0][1] = matrix[1][0] = 0.0;
00146         matrix[0][2] = matrix[2][0] = 0.0;
00147         matrix[1][2] = matrix[2][1] = 0.0;
00148         return;
00149     }
00150     vector /= len;
00151 
00152     const double cosine = plane_normal % Z;
00153     const double sine   = sqrt( 1 - cosine * cosine );
00154 
00155     std::cerr << "Rotation: " << acos( cosine ) << " [" << vector.x << ' ' << vector.y << ' ' << vector.z << ']'
00156               << std::endl;
00157 
00158     const double x = vector.x;
00159     const double y = vector.y;
00160     const double z = vector.z;
00161     const double c = cosine;
00162     const double s = sine;
00163     const double o = 1.0 - cosine;
00164 
00165     matrix[0][0] = c + x * x * o;
00166     matrix[0][1] = -z * s + x * y * o;
00167     matrix[0][2] = y * s + x * z * o;
00168 
00169     matrix[1][0] = z * s + x * z * o;
00170     matrix[1][1] = c + y * y * o;
00171     matrix[1][2] = -x * s + y * z * o;
00172 
00173     matrix[2][0] = -y * s + x * z * o;
00174     matrix[2][1] = x * s + y * z * o;
00175     matrix[2][2] = c + z * z * o;
00176 }
00177 
00178 static void transform_point( CartVect3D& point, double matrix[3][3] )
00179 {
00180     const double x = point.x;
00181     const double y = point.y;
00182     const double z = point.z;
00183 
00184     point.x = x * matrix[0][0] + y * matrix[0][1] + z * matrix[0][2];
00185     point.y = x * matrix[1][0] + y * matrix[1][1] + z * matrix[1][2];
00186     point.z = x * matrix[2][0] + y * matrix[2][1] + z * matrix[2][2];
00187 }
00188 
00189 static void write_gnuplot( std::ostream& stream, const std::vector< CartVect3D >& list );
00190 
00191 static void write_svg( std::ostream& stream, const std::vector< CartVect3D >& list );
00192 
00193 static void write_eps( std::ostream& stream, const std::vector< CartVect3D >& list, int surface_id );
00194 
00195 enum FileType
00196 {
00197     POSTSCRIPT,
00198     GNUPLOT,
00199     SVG
00200 };
00201 
00202 using namespace moab;
00203 
00204 int main( int argc, char* argv[] )
00205 {
00206     Interface* moab = new Core();
00207     ErrorCode result;
00208     std::vector< CartVect3D >::iterator iter;
00209     FileType type = GNUPLOT;
00210 
00211     int idx = 1;
00212     if( argc == 4 )
00213     {
00214         if( !strcmp( argv[idx], "-p" ) )
00215             type = POSTSCRIPT;
00216         else if( !strcmp( argv[idx], "-g" ) )
00217             type = GNUPLOT;
00218         else if( !strcmp( argv[idx], "-s" ) )
00219             type = SVG;
00220         else
00221             usage_error( argv[0] );
00222         ++idx;
00223     }
00224 
00225     // scan CL args
00226     int surface_id;
00227     if( argc - idx != 2 ) usage_error( argv[0] );
00228     char* endptr;
00229     surface_id = strtol( argv[idx], &endptr, 0 );
00230     if( !endptr || *endptr ) usage_error( argv[0] );
00231     ++idx;
00232 
00233     // Load mesh
00234     result = moab->load_mesh( argv[idx] );
00235     if( MB_SUCCESS != result )
00236     {
00237         if( MB_FILE_DOES_NOT_EXIST == result )
00238             std::cerr << argv[idx] << " : open failed.\n";
00239         else
00240             std::cerr << argv[idx] << " : error reading file.\n";
00241         return READ_ERROR;
00242     }
00243 
00244     // Get tag handles
00245     EntityHandle surface;
00246     const int dimension = 2;  // surface
00247     if( surface_id )
00248     {
00249         Tag tags[2];
00250         result = moab->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tags[0] );
00251         if( MB_SUCCESS != result )
00252         {
00253             std::cerr << "No geometry tag.\n";
00254             return OTHER_ERROR;
00255         }
00256         tags[1] = moab->globalId_tag();
00257 
00258         // Find entityset for surface.
00259         const void* tag_values[] = { &dimension, &surface_id };
00260         Range surfaces;
00261         moab->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, tag_values, 2, surfaces );
00262         if( surfaces.size() != 1 )
00263         {
00264             std::cerr << "Found " << surfaces.size() << " surfaces with ID " << surface_id << std::endl;
00265             return SURFACE_NOT_FOUND;
00266         }
00267         surface = *surfaces.begin();
00268     }
00269     else
00270     {
00271         surface = 0;
00272     }
00273 
00274     // Get surface mesh
00275     Range elements;
00276     result = moab->get_entities_by_dimension( surface, dimension, elements );
00277     if( MB_SUCCESS != result )
00278     {
00279         std::cerr << "Internal error\n";
00280         return OTHER_ERROR;
00281     }
00282 
00283     // Calculate average corner normal in surface mesh
00284     CartVect3D normal( 0, 0, 0 );
00285     std::vector< EntityHandle > vertices;
00286     std::vector< CartVect3D > coords;
00287     for( Range::iterator i = elements.begin(); i != elements.end(); ++i )
00288     {
00289         vertices.clear();
00290         result = moab->get_connectivity( &*i, 1, vertices, true );
00291         if( MB_SUCCESS != result )
00292         {
00293             std::cerr << "Internal error\n";
00294             return OTHER_ERROR;
00295         }
00296         coords.clear();
00297         coords.resize( vertices.size() );
00298         result = moab->get_coords( &vertices[0], vertices.size(), reinterpret_cast< double* >( &coords[0] ) );
00299         if( MB_SUCCESS != result )
00300         {
00301             std::cerr << "Internal error\n";
00302             return OTHER_ERROR;
00303         }
00304 
00305         for( size_t j = 0; j < coords.size(); ++j )
00306         {
00307             CartVect3D v1 = coords[( j + 1 ) % coords.size()] - coords[j];
00308             CartVect3D v2 = coords[( j + 1 ) % coords.size()] - coords[( j + 2 ) % coords.size()];
00309             normal += ( v1 * v2 );
00310         }
00311     }
00312     normal /= normal.len();
00313 
00314     // Get edges from elements
00315     Range edge_range;
00316     result = moab->get_adjacencies( elements, 1, true, edge_range, Interface::UNION );
00317     if( MB_SUCCESS != result )
00318     {
00319         std::cerr << "Internal error\n";
00320         return OTHER_ERROR;
00321     }
00322 
00323     // Get vertex coordinates for each edge
00324     std::vector< EntityHandle > edges( edge_range.size() );
00325     std::copy( edge_range.begin(), edge_range.end(), edges.begin() );
00326     vertices.clear();
00327     result = moab->get_connectivity( &edges[0], edges.size(), vertices, true );
00328     if( MB_SUCCESS != result )
00329     {
00330         std::cerr << "Internal error\n";
00331         return OTHER_ERROR;
00332     }
00333     coords.clear();
00334     coords.resize( vertices.size() );
00335     result = moab->get_coords( &vertices[0], vertices.size(), reinterpret_cast< double* >( &coords[0] ) );
00336     if( MB_SUCCESS != result )
00337     {
00338         std::cerr << "Internal error\n";
00339         return OTHER_ERROR;
00340     }
00341 
00342     // Rotate points such that the projection into the view plane
00343     // can be accomplished by discarding the 'z' coordinate of each
00344     // point.
00345 
00346     std::cerr << "Plane normal: [" << normal.x << ' ' << normal.y << ' ' << normal.z << ']' << std::endl;
00347     double transform[3][3];
00348     find_rotation( normal, transform );
00349 
00350     for( iter = coords.begin(); iter != coords.end(); ++iter )
00351         transform_point( *iter, transform );
00352 
00353     // Write the file.
00354 
00355     switch( type )
00356     {
00357         case POSTSCRIPT:
00358             write_eps( std::cout, coords, surface_id );
00359             break;
00360         case SVG:
00361             write_svg( std::cout, coords );
00362             break;
00363         default:
00364             write_gnuplot( std::cout, coords );
00365             break;
00366     }
00367     return 0;
00368 }
00369 
00370 void write_gnuplot( std::ostream& stream, const std::vector< CartVect3D >& coords )
00371 {
00372     std::vector< CartVect3D >::const_iterator iter;
00373 
00374     stream << std::endl;
00375     for( iter = coords.begin(); iter != coords.end(); ++iter )
00376     {
00377         stream << iter->x << ' ' << iter->y << std::endl;
00378         ++iter;
00379         if( iter == coords.end() )
00380         {
00381             stream << std::endl;
00382             break;
00383         }
00384         stream << iter->x << ' ' << iter->y << std::endl;
00385         stream << std::endl;
00386     }
00387     std::cerr << "Display with gnuplot command \"plot with lines\"\n";
00388 }
00389 
00390 static void box_max( CartVect3D& cur_max, const CartVect3D& pt )
00391 {
00392     if( pt.x > cur_max.x ) cur_max.x = pt.x;
00393     if( pt.y > cur_max.y ) cur_max.y = pt.y;
00394     // if (pt.z > cur_max.z)
00395     //  cur_max.z = pt.z;
00396 }
00397 
00398 static void box_min( CartVect3D& cur_min, const CartVect3D& pt )
00399 {
00400     if( pt.x < cur_min.x ) cur_min.x = pt.x;
00401     if( pt.y < cur_min.y ) cur_min.y = pt.y;
00402     // if (pt.z > cur_min.z)
00403     //  cur_min.z = pt.z;
00404 }
00405 
00406 void write_eps( std::ostream& s, const std::vector< CartVect3D >& coords, int id )
00407 {
00408     // Coordinate range to use within EPS file
00409     const int X_MAX = 540;  // 540 pts / 72 pts/inch = 7.5 inches
00410     const int Y_MAX = 720;  // 720 pts / 72 pts/inch = 10 inches
00411 
00412     std::vector< CartVect3D >::const_iterator iter;
00413 
00414     // Get bounding box
00415     const double D_MAX = std::numeric_limits< double >::max();
00416     CartVect3D min( D_MAX, D_MAX, 0 );
00417     CartVect3D max( -D_MAX, -D_MAX, 0 );
00418     for( iter = coords.begin(); iter != coords.end(); ++iter )
00419     {
00420         box_max( max, *iter );
00421         box_min( min, *iter );
00422     }
00423 
00424     // Calculate translation to page coordinates
00425     CartVect3D offset = CartVect3D( 0, 0, 0 ) - min;
00426     CartVect3D scale  = max - min;
00427     scale.x           = X_MAX / scale.x;
00428     scale.y           = Y_MAX / scale.y;
00429     if( scale.x > scale.y )  // keep proportions
00430         scale.x = scale.y;
00431     else
00432         scale.y = scale.x;
00433 
00434     // std::cerr << "Min: " << min.x << ' ' << min.y <<
00435     //           "  Max: " << max.x << ' ' << max.y << std::endl
00436     //          << "Offset: " << offset.x << ' ' << offset.y <<
00437     //           "  Scale: " << scale.x << ' ' << scale.y << std::endl;
00438 
00439     // Write header stuff
00440     s << "%!PS-Adobe-2.0 EPSF-2.0" << std::endl;
00441     s << "%%Creator: MOAB surfplot" << std::endl;
00442     s << "%%Title: Surface " << id << std::endl;
00443     s << "%%DocumentData: Clean7Bit" << std::endl;
00444     s << "%%Origin: 0 0" << std::endl;
00445     int max_x = (int)( ( max.x + offset.x ) * scale.x );
00446     int max_y = (int)( ( max.y + offset.y ) * scale.y );
00447     s << "%%BoundingBox: 0 0 " << max_x << ' ' << max_y << std::endl;
00448     s << "%%Pages: 1" << std::endl;
00449 
00450     s << "%%BeginProlog" << std::endl;
00451     s << "save" << std::endl;
00452     s << "countdictstack" << std::endl;
00453     s << "mark" << std::endl;
00454     s << "newpath" << std::endl;
00455     s << "/showpage {} def" << std::endl;
00456     s << "/setpagedevice {pop} def" << std::endl;
00457     s << "%%EndProlog" << std::endl;
00458 
00459     s << "%%Page: 1 1" << std::endl;
00460     s << "1 setlinewidth" << std::endl;
00461     s << "0.0 setgray" << std::endl;
00462 
00463     for( iter = coords.begin(); iter != coords.end(); ++iter )
00464     {
00465         double x1 = ( iter->x + offset.x ) * scale.x;
00466         double y1 = ( iter->y + offset.y ) * scale.y;
00467         if( ++iter == coords.end() ) break;
00468         double x2 = ( iter->x + offset.x ) * scale.x;
00469         double y2 = ( iter->y + offset.y ) * scale.y;
00470 
00471         s << "newpath" << std::endl;
00472         s << x1 << ' ' << y1 << " moveto" << std::endl;
00473         s << x2 << ' ' << y2 << " lineto" << std::endl;
00474         s << "stroke" << std::endl;
00475     }
00476 
00477     s << "%%Trailer" << std::endl;
00478     s << "cleartomark" << std::endl;
00479     s << "countdictstack" << std::endl;
00480     s << "exch sub { end } repeat" << std::endl;
00481     s << "restore" << std::endl;
00482     s << "%%EOF" << std::endl;
00483 }
00484 
00485 void write_svg( std::ostream& file, const std::vector< CartVect3D >& coords )
00486 {
00487 
00488     std::vector< CartVect3D >::const_iterator iter;
00489 
00490     // Get bounding box
00491     const double D_MAX = std::numeric_limits< double >::max();
00492     CartVect3D min( D_MAX, D_MAX, 0 );
00493     CartVect3D max( -D_MAX, -D_MAX, 0 );
00494     for( iter = coords.begin(); iter != coords.end(); ++iter )
00495     {
00496         box_max( max, *iter );
00497         box_min( min, *iter );
00498     }
00499     CartVect3D size = max - min;
00500 
00501     // scale to 640 pixels on a side
00502     double scale = 640.0 / ( size.x > size.y ? size.x : size.y );
00503     size *= scale;
00504 
00505     file << "<?xml version=\"1.0\" standalone=\"no\"?>" << std::endl;
00506     file << "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" "
00507          << "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << std::endl;
00508     file << std::endl;
00509     file << "<svg width=\"" << (int)size.x << "\" height=\"" << (int)size.y
00510          << "\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" << std::endl;
00511 
00512     int left = (int)( min.x * scale );
00513     int top  = (int)( min.y * scale );
00514     iter     = coords.begin();
00515     while( iter != coords.end() )
00516     {
00517         file << "<line "
00518              << "x1=\"" << (int)( scale * iter->x ) - left << "\" "
00519              << "y1=\"" << (int)( scale * iter->y ) - top << "\" ";
00520         ++iter;
00521         file << "x2=\"" << (int)( scale * iter->x ) - left << "\" "
00522              << "y2=\"" << (int)( scale * iter->y ) - top << "\" "
00523              << " style=\"stroke:rgb(99,99,99);stroke-width:2\""
00524              << "/>" << std::endl;
00525         ++iter;
00526     }
00527 
00528     // Write footer
00529     file << "</svg>" << std::endl;
00530 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines