An Implementation of the Marching Cubes Algorithm
1.0
|
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.