An Implementation of the Marching Cubes Algorithm  1.0
surfbuilder.cpp
Go to the documentation of this file.
00001 
00026 #include "surfbuilder.hpp"    // mc::SurfBuilder
00027 
00028 #include "mathmacros.hpp"     // SIGN
00029 
00030 #include <sstream>            // std::stringstream
00031 
00032 
00049 namespace mc {
00050 
00051 
00052   using common::t3DVector ;
00053   using common::ExceptionObject ;
00054 
00055 
00056   // -----------------------------------------------------------------
00057   //
00058   // Public methods
00059   //
00060   // -----------------------------------------------------------------
00061 
00069   void 
00070   SurfBuilder::run() 
00071   {
00072     /*
00073      * Create a  table to  store all possible  ways of  traversing the
00074      * edges of a grid cube  to create the triangles of the simplicial
00075      * surface.
00076      */
00077     make_lookup_table() ;
00078 
00079     /*
00080      * For each cube belonging to the grid and containing the vertices
00081      * (i,j,k),   (i+1,j,k),    (i,j+1,k),   (i+1,j+1,k),   (i,j,k+1),
00082      * (i+1,j,k+1), (i,j+1,k+1),  (i+1,j+1,k+1), obtain the  values of
00083      * the unknown implicit function at  the vertices of the cube, and
00084      * then compute the intersection  points, if any, between the cube
00085      * edges and the surface defined by the unknown implicit function.
00086      */
00087 
00088     /*
00089      * Get the number of grid points in each dimension of the grid. 
00090      */
00091     unsigned sx = _grid->get_size_x() - 1 ;
00092     unsigned sy = _grid->get_size_y() - 1 ;
00093     unsigned sz = _grid->get_size_z() - 1 ;
00094 
00095     /*
00096      * Create an array to store the function values.
00097      */
00098     double fv[ 8 ] ;
00099 
00100     /*
00101      * Visit each cube of the grid and compute the intersection points
00102      * between  its  edges and  the  surface  defined  by the  unknown
00103      * implicit function.
00104      */
00105     for ( unsigned k = 0 ; k < sz ; k++ ) {
00106       for ( unsigned j = 0 ; j < sy ; j++ ) {
00107         for ( unsigned i = 0 ; i < sx ; i++ ) {
00108           /*
00109            * Get the value of the implicit function at the vertices of
00110            * the cube whose origin vertex is given by the triple index
00111            * (i,j,k).
00112            */
00113           fv[ 0 ] = _grid->get_value(     i ,     j ,     k ) ;
00114           fv[ 1 ] = _grid->get_value(     i ,     j , k + 1 ) ;
00115           fv[ 2 ] = _grid->get_value(     i , j + 1 ,     k ) ;
00116           fv[ 3 ] = _grid->get_value(     i , j + 1 , k + 1 ) ;
00117           fv[ 4 ] = _grid->get_value( i + 1 ,     j ,     k ) ;
00118           fv[ 5 ] = _grid->get_value( i + 1 ,     j , k + 1 ) ;
00119           fv[ 6 ] = _grid->get_value( i + 1 , j + 1 ,     k ) ;
00120           fv[ 7 ] = _grid->get_value( i + 1 , j + 1 , k + 1 ) ;
00121 
00122           /*
00123            * Determine an entry of the  lookup table by using the sign
00124            * of the values of the implicit function at the vertices of
00125            * the cube.
00126            */
00127           unsigned n , index = 0 ;
00128           for ( n = 0 ; n < 8 ; n++ ) {
00129             if ( fv[ n ] > 0 ) {
00130               index += unsigned( 1 << n ) ;
00131             }
00132           }
00133 
00134           /*
00135            * Traverse the faces of the cube according to the traversal
00136            * stored in  the lookup table  at the table  entry computed
00137            * above.
00138            */
00139           for (
00140                TABLEENTRY::const_iterator tit = _cubetable[ index ].begin() ;
00141                tit != _cubetable[ index ].end() ;
00142                ++tit
00143               ) 
00144           {
00145             /*
00146              * Get the  edges intersected  by the implicit  surface in
00147              * the order  they appear in the traversals  stored in the
00148              * table.
00149              */
00150             for (
00151                  TRAVERSAL::const_iterator eit = (*tit).begin() ;
00152                  eit != (*tit).end() ;
00153                  eit++
00154                  ) 
00155             {
00156               /*
00157                * Get the edge index.
00158                */
00159               EDGES ie = (EDGES) *eit ;
00160 
00161               /*
00162                * Get the indices of the two vertices incident to the edge.
00163                */
00164               VERTS v1 = get_1st_vertex( ie ) ;
00165               VERTS v2 = get_2nd_vertex( ie ) ;
00166 
00167               /*
00168                * Compute the  intersection point between  the edge and
00169                * the surface defined by the unknown implicit function.
00170                */
00171               t3DVector pt ;
00172               try {
00173                 compute_intersection(
00174                                      i ,
00175                                      j ,
00176                                      k ,
00177                                      v1 ,
00178                                      v2 , 
00179                                      fv[ v1 ] ,
00180                                      fv[ v2 ] ,
00181                                      pt
00182                                     ) ;
00183               }
00184               catch ( const ExceptionObject& xpt ) {
00185                 treat_exception( xpt ) ;
00186                 exit( 0 ) ;
00187               }
00188 
00189               // DO SOMETHING WITH THE POINT JUST COMPUTED
00190             }
00191 
00192             // DO SOMETHING WITH ALL POINTS OF INTERSECTION (THEY DEFINE ONE PATCH)
00193           }
00194         }
00195       }
00196     }
00197 
00198     return ;
00199   }
00200 
00201 
00212   void
00213   SurfBuilder::make_lookup_table()
00214   {
00215     /*
00216      * For each  possible way of assigning  one of two  values for the
00217      * eight vertices  of a cube,  store traversal information  into a
00218      * table entry.  Each traversal contains a sequence  of cube edges
00219      * intersected  by the  surface  defined by  the unknown  implicit
00220      * function.
00221      */
00222     for ( unsigned i = 0 ; i < 256 ; i++ ) {
00223       /*
00224        * Initialize an array to indicate  whether a grid cube edge has
00225        * been visited during one traversal.
00226        */
00227       bool done[ 12 ] ;
00228       for ( unsigned e = 0 ; e < 12 ; e++ ) {
00229         done[ e ] = false ;
00230       }
00231 
00232       /*
00233        * Compute  an  array with  powers  of  2  corresponding to  the
00234        * indices of the vertices of  the same cube. The vertex indices
00235        * are 0, 1, ..., 7, and the respective power of 2 are 2^0, 2^1,
00236        * ..., 2^7.
00237        */
00238       bool pos[ 8 ] ;
00239       for ( unsigned j = 0 ; j < 8 ; j++ ) {
00240         pos[ j ] = compute_sign( i , j ) ;
00241       }
00242 
00243       /*
00244        * Compute a traversal  of the cube that determines  the way the
00245        * points inside the cube must  be connected to form one or more
00246        * triangles.
00247        */
00248       for ( unsigned e = 0 ; e < 12 ; e++ ) {
00249         /*
00250          * For each edge of the cube that has not been visited yet and
00251          * whose  vertices   have  opposite  signs,   initiate  a  new
00252          * traversal  and insert  the  traversal into  a lookup  table
00253          * entry.
00254          */
00255         EDGES ie = ( EDGES ) e ;
00256 
00257         if ( !done[ ie ] && 
00258              ( pos[ get_1st_vertex( ie ) ] != pos[ get_2nd_vertex( ie ) ] ) 
00259            ) 
00260         {
00261           EDGES stae = ie ;
00262           EDGES edge = ie ;
00263 
00264           /*
00265            * Obtain the cube face on  the right side of the edge "ie",
00266            * where  "side"  is  defined   with  respect  to  the  edge
00267            * orientation given  from the vertex with  positive sign to
00268            * the vertex with negative sign.
00269            */
00270           FACES face = pos[ get_1st_vertex( ie ) ] ? 
00271                        right_face( ie ) : left_face( ie ) ;
00272 
00273           /*
00274            * Create a traversal.
00275            */
00276           TRAVERSAL curr ;
00277 
00278           /*
00279            * Compute  a traversal along  the faces  of the  cube. This
00280            * traversal  consists  of  edges  whose  the  sign  of  the
00281            * implicit  function  at  the  two  incident  vertices  are
00282            * opposite.
00283            */
00284           do {
00285             /*
00286              * Pick the next edge in the current face.
00287              */
00288             edge = next_cw_edge( edge , face ) ;
00289 
00290             /*
00291              * Mark the new current edge as visited.
00292              */
00293             done[ edge ] = true ;
00294 
00295             /*
00296              * If the sign of the implicit function at the vertices of
00297              * the current  edge are opposite, insert the  edge in the
00298              * current traversal and change face to find the remaining
00299              * edges.
00300              */
00301             if ( pos[ get_1st_vertex( edge ) ] != pos[ get_2nd_vertex( edge ) ] ) {
00302               curr.push_back( edge ) ;
00303               face = change_face(edge , face ) ;
00304             }
00305 
00306           } while ( stae != edge ) ;
00307 
00308           /*
00309            * Insert the new traversal into the lookup table.
00310            */
00311           _cubetable[ i ].push_back( curr ) ;
00312         }
00313       }
00314     }
00315 
00316     return ;
00317   }
00318 
00319 
00329   SurfBuilder::VERTS 
00330   SurfBuilder::get_1st_vertex( EDGES e ) const
00331   {
00332     VERTS index ;
00333 
00334     if ( e == LB ) {
00335       index = LBN ;
00336     }
00337     else if ( e == LT ) {
00338       index = LTN ;
00339     }
00340     else if ( e == LN ) {
00341       index = LBN ;
00342     }
00343     else if ( e == LF ) {
00344       index = LBF ;
00345     }
00346     else if ( e == RB ) {
00347       index = RBN ;
00348     }
00349     else if ( e == RT ) {
00350       index = RTN ;
00351     }
00352     else if ( e == RN ) {
00353       index = RBN ;
00354     }
00355     else if ( e == RF ) {
00356       index = RBF ;
00357     }
00358     else if ( e == BN ) {
00359       index = LBN ;
00360     }
00361     else if ( e == BF ) {
00362       index = LBF ;
00363     }
00364     else if ( e == TN ) {
00365       index = LTN ;
00366     }
00367     else { // if ( e == TF ) {
00368       index = LTF ;
00369     }
00370 
00371     return index ;
00372   }
00373 
00374 
00384   SurfBuilder::VERTS
00385   SurfBuilder::get_2nd_vertex( EDGES e ) const
00386   {
00387     VERTS index ;
00388 
00389     if ( e == LB ) {
00390       index = LBF ;
00391     }
00392     else if ( e == LT ) {
00393       index = LTF ;
00394     }
00395     else if ( e == LN ) {
00396       index = LTN ;
00397     }
00398     else if ( e == LF ) {
00399       index = LTF ;
00400     }
00401     else if ( e == RB ) {
00402       index = RBF ;
00403     }
00404     else if ( e == RT ) {
00405       index = RTF ;
00406     }
00407     else if ( e == RN ) {
00408       index = RTN ;
00409     }
00410     else if ( e == RF ) {
00411       index = RTF ;
00412     }
00413     else if ( e == BN ) {
00414       index = RBN ;
00415     }
00416     else if ( e == BF ) {
00417       index = RBF ;
00418     }
00419     else if ( e == TN ) {
00420       index = RTN ;
00421     }
00422     else { // if ( e == TF ) {
00423       index = RTF ;
00424     }
00425 
00426     return index ;
00427   }
00428 
00429 
00441   SurfBuilder::FACES
00442   SurfBuilder::left_face( EDGES e ) const
00443   {
00444     FACES index ;
00445 
00446     if ( e == LB ) {
00447       index = B ;
00448     }
00449     else if ( e == LT ) {
00450       index = L ;
00451     }
00452     else if ( e == LN ) {
00453       index = L ;
00454     }
00455     else if ( e == LF ) {
00456       index = F ;
00457     }
00458     else if ( e == RB ) {
00459       index = R ;
00460     }
00461     else if ( e == RT ) {
00462       index = T ;
00463     }
00464     else if ( e == RN ) {
00465       index = N ;
00466     }
00467     else if ( e == RF ) {
00468       index = R ;
00469     }
00470     else if ( e == BN ) {
00471       index = N ;
00472     }
00473     else if ( e == BF ) {
00474       index = B ;
00475     }
00476     else if ( e == TN ) {
00477       index = T ;
00478     }
00479     else { // if ( e == TF ) {
00480       index = F ;
00481     }
00482 
00483     return index ;
00484   }
00485 
00486 
00498   SurfBuilder::FACES 
00499   SurfBuilder::right_face( EDGES e ) const
00500   {
00501     FACES index ;
00502 
00503     if ( e == LB ) {
00504       index = L ;
00505     }
00506     else if ( e == LT ) {
00507       index = T ;
00508     }
00509     else if ( e == LN ) {
00510       index = N ;
00511     }
00512     else if ( e == LF ) {
00513       index = L ;
00514     }
00515     else if ( e == RB ) {
00516       index = B ;
00517     }
00518     else if ( e == RT ) {
00519       index = R ;
00520     }
00521     else if ( e == RN ) {
00522       index = R ;
00523     }
00524     else if ( e == RF ) {
00525       index = F ;
00526     }
00527     else if ( e == BN ) {
00528       index = B ;
00529     }
00530     else if ( e == BF ) {
00531       index = F ;
00532     }
00533     else if ( e == TN ) {
00534       index = N ;
00535     }
00536     else { // if ( e == TF ) {
00537       index = T ;
00538     }
00539 
00540     return index ;
00541   }
00542 
00543 
00559   SurfBuilder::EDGES
00560   SurfBuilder::next_cw_edge( EDGES e , FACES f ) const
00561   {
00562     EDGES index ;
00563 
00564     if ( e == LB ) {
00565       index = ( f == L ) ? LF : BN ;
00566     }
00567     else if ( e == LT ) {
00568       index = ( f == L ) ? LN : TF ;
00569     }
00570     else if ( e == LN ) {
00571       index = ( f == L ) ? LB : TN ;
00572     }
00573     else if ( e == LF ) {
00574       index = ( f == L ) ? LT : BF ;
00575     }
00576     else if ( e == RB ) {
00577       index = ( f == R ) ? RN : BF ;
00578     }
00579     else if ( e == RT ) {
00580       index = ( f == R ) ? RF : TN ;
00581     }
00582     else if ( e == RN ) {
00583       index = ( f == R ) ? RT : BN ;
00584     }
00585     else if ( e == RF ) {
00586       index = ( f == R ) ? RB : TF ;
00587     }
00588     else if ( e == BN ) {
00589       index = ( f == B ) ? RB : LN ;
00590     }
00591     else if ( e == BF ) {
00592       index = ( f == B ) ? LB : RF ;
00593     }
00594     else if ( e == TN ) {
00595       index = ( f == T ) ? LT : RN ;
00596     }
00597     else { // if ( e == TF ) {
00598       index = ( f == T ) ? RT : LF ;
00599     }
00600 
00601     return index ;
00602   }
00603 
00604 
00628   void
00629   SurfBuilder::compute_intersection(
00630                                     unsigned i ,
00631                                     unsigned j ,
00632                                     unsigned k ,
00633                                     VERTS v1 ,
00634                                     VERTS v2 ,
00635                                     double f1 ,
00636                                     double f2 ,
00637                                     t3DVector& pt
00638                                    )
00639     const throw( ExceptionObject )
00640   {
00641     if ( SIGN( f1 ) == SIGN( f2 ) ) {
00642       std::stringstream ss ( std::stringstream::in | 
00643                              std::stringstream::out ) ;
00644       ss << "Trying to compute intersection point on non-interesecting edge" ;
00645       throw ExceptionObject( __FILE__ , __LINE__ , ss.str() ) ;
00646     }
00647 
00648     unsigned i1 , j1 , k1 , i2 , j2 , k2 ;
00649 
00650     /*
00651      * Initialize indices of "v1" and "v2".
00652      */
00653     i1 = i2 = i ;
00654     j1 = j2 = j ;
00655     k1 = k2 = k ;
00656 
00657     /*
00658      * Compute the indices of "v1" and "v2" with respect to the grid.
00659      */
00660     if ( ( v1 == RBN ) || ( v1 == RBF ) || ( v1 == RTN ) || ( v1 == RTF ) ) {
00661       i1 += 1 ;
00662     }
00663 
00664     if ( ( v1 == LTN ) || ( v1 == LTF ) || ( v1 == RTN ) || ( v1 == RTF ) ) {
00665       j1 += 1 ;
00666     }
00667 
00668     if ( ( v1 == LBF ) || ( v1 == LTF ) || ( v1 == RBF ) || ( v1 == RTF ) ) {
00669       k1 += 1 ;
00670     }
00671 
00672     if ( ( v2 == RBN ) || ( v2 == RBF ) || ( v2 == RTN ) || ( v2 == RTF ) ) {
00673       i2 += 1 ;
00674     }
00675 
00676     if ( ( v2 == LTN ) || ( v2 == LTF ) || ( v2 == RTN ) || ( v2 == RTF ) ) {
00677       j2 += 1 ;
00678     }
00679 
00680     if ( ( v2 == LBF ) || ( v2 == LTF ) || ( v2 == RBF ) || ( v2 == RTF ) ) {
00681       k2 += 1 ;
00682     }
00683 
00684     /*
00685      * DO SOMETHING TO AVOID COMPUTING THE SAME VERTEX TWICE!
00686      */
00687 
00688     /*
00689      * Obtain the coordinates of vertices "v1" and "v2".
00690      */
00691     t3DVector go = _grid->get_origin() ;
00692 
00693     double spacx = _grid->get_spacing_x() ;
00694     double spacy = _grid->get_spacing_y() ;
00695     double spacz = _grid->get_spacing_z() ;
00696  
00697     t3DVector p1 = go + t3DVector( spacx * i1 , spacy * j1 , spacz * k1 ) ;
00698     t3DVector p2 = go + t3DVector( spacx * i2 , spacy * j2 , spacz * k2 ) ;
00699 
00700     /*
00701      * Compute position  of intersection  point bewteen the  cube edge
00702      * and the surface defined  by the unknown implicit function using
00703      * affine interpolation.
00704      */
00705     const double coeff = f2 / ( f2 - f1 ) ;
00706     pt = ( coeff * p1 ) + ( ( 1 - coeff ) * p2 ) ;
00707 
00708     // DO SOMETHING WITH THE INTERSECTION POINT
00709 
00710     return ;
00711   }
00712 
00713 }         // end of namespace mc
00714  //end of group class.