shape_reconstruction: marching

Go to the documentation of this file.
00001 /****************************************************************************
00002 * VCGLib                                                            o o     *
00003 * Visual and Computer Graphics Library                            o     o   *
00004 *                                                                _   O  _   *
00005 * Copyright(C) 2004                                                \/)\/    *
00006 * Visual Computing Lab                                            /\/|      *
00007 * ISTI - Italian National Research Council                           |      *
00008 *                                                                    \      *
00009 * All rights reserved.                                                      *
00010 *                                                                           *
00011 * This program is free software; you can redistribute it and/or modify      *
00012 * it under the terms of the GNU General Public License as published by      *
00013 * the Free Software Foundation; either version 2 of the License, or         *
00014 * (at your option) any later version.                                       *
00015 *                                                                           *
00016 * This program is distributed in the hope that it will be useful,           *
00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
00019 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
00020 * for more details.                                                         *
00021 *                                                                           *
00022 ****************************************************************************/
00023 /***************************************************************************/
00024 
00025 #ifndef __VCG_MARCHING_CUBES
00026 #define __VCG_MARCHING_CUBES
00027 
00028 #include "mc_lookup_table.h"
00029 
00030 namespace vcg
00031 {
00032     namespace tri
00033     {
00034         // Doxygen documentation
00037 
00038         /*
00039         * Cube description:
00040         *         3 ________ 2           _____2__
00041         *         /|       /|         / |       /|
00042         *       /  |     /  |      11/  3   10/  |
00043         *   7 /_______ /    |      /__6_|__ /    |1
00044         *    |     |  |6    |     |     |  |     |
00045         *    |    0|__|_____|1    |     |__|__0__|
00046         *    |    /   |    /      7   8/   5    /
00047         *    |  /     |  /        |  /     |  /9
00048         *    |/_______|/          |/___4___|/
00049         *   4          5
00050         */
00051 
00053 
00061         template<class TRIMESH_TYPE, class WALKER_TYPE>
00062         class MarchingCubes
00063         {
00064         public:
00065             enum Dimension               {X, Y, Z};
00066             typedef typename vcg::tri::Allocator< TRIMESH_TYPE > AllocatorType;
00067             typedef typename TRIMESH_TYPE::ScalarType                   ScalarType;
00068             typedef typename TRIMESH_TYPE::VertexType                   VertexType;
00069             typedef typename TRIMESH_TYPE::VertexPointer        VertexPointer;
00070             typedef typename TRIMESH_TYPE::VertexIterator       VertexIterator;
00071             typedef typename TRIMESH_TYPE::FaceType                             FaceType;
00072             typedef typename TRIMESH_TYPE::FacePointer          FacePointer;
00073             typedef typename TRIMESH_TYPE::FaceIterator         FaceIterator;
00074             typedef typename TRIMESH_TYPE::CoordType                    CoordType;
00075             typedef typename TRIMESH_TYPE::CoordType*                   CoordPointer;
00076 
00082             MarchingCubes(TRIMESH_TYPE &mesh, WALKER_TYPE &walker)
00083             {
00084                 _mesh           = &mesh;
00085                 _walker = &walker;
00086             };
00087 
00092             void Initialize()
00093             {
00094                 _mesh->Clear();
00095             }; // end of Initialize()
00096 
00101             void Finalize()
00102             {
00103                 _mesh           = NULL;
00104                 _walker = NULL;
00105             }; // end of Finalize()
00106 
00113             void ProcessCell(const vcg::Point3i &min, const vcg::Point3i &max)
00114             {
00115                 _case = _subconfig = _config = -1;
00116                 assert(min[0]<max[0] && min[1]<max[1] && min[2]<max[2]);
00117                 _corners[0].X()=min.X();                _corners[0].Y()=min.Y();                _corners[0].Z()=min.Z();
00118                 _corners[1].X()=max.X();                _corners[1].Y()=min.Y();                _corners[1].Z()=min.Z();
00119                 _corners[2].X()=max.X();                _corners[2].Y()=max.Y();                _corners[2].Z()=min.Z();
00120                 _corners[3].X()=min.X();                _corners[3].Y()=max.Y();                _corners[3].Z()=min.Z();
00121                 _corners[4].X()=min.X();                _corners[4].Y()=min.Y();                _corners[4].Z()=max.Z();
00122                 _corners[5].X()=max.X();                _corners[5].Y()=min.Y();                _corners[5].Z()=max.Z();
00123                 _corners[6].X()=max.X();                _corners[6].Y()=max.Y();                _corners[6].Z()=max.Z();
00124                 _corners[7].X()=min.X();                _corners[7].Y()=max.Y();                _corners[7].Z()=max.Z();
00125 
00126                 for (int i=0; i<8; i++)
00127                     _field[i] = _walker->V( _corners[i].X(), _corners[i].Y(), _corners[i].Z() );
00128 
00129                 unsigned char cubetype = 0;
00130                 for (int i=0; i<8; i++)
00131                     if (_field[i]>0) cubetype += 1<<i;
00132 
00133                 _case                           = MCLookUpTable::Cases(cubetype, 0); //_case                            = cases[cubetype][0];
00134                 _config                 = MCLookUpTable::Cases(cubetype, 1); //_config                  = cases[cubetype][1];
00135                 _subconfig      = 0;
00136 
00137                 VertexPointer v12 = NULL;
00138 
00139                 switch( _case )
00140                 {
00141                 case  0 : { break ; }
00142                 case  1 : { AddTriangles( MCLookUpTable::Tiling1(_config), 1 ); break; } //case  1 : { AddTriangles( tiling1[_config], 1 ); break; }
00143                 case  2 : { AddTriangles( MCLookUpTable::Tiling2(_config), 2 ); break; } //case  2 : { AddTriangles( tiling2[_config], 2 ); break; }
00144                 case  3 :
00145                     {
00146                         //if( TestFace( test3[_config]) ) AddTriangles( tiling3_2[_config], 4 ) ; // 3.2
00147                         if( TestFace( MCLookUpTable::Test3(_config)) )
00148                             AddTriangles( MCLookUpTable::Tiling3_2(_config), 4 ) ; // 3.2
00149                         else
00150                             AddTriangles( MCLookUpTable::Tiling3_1(_config), 2 ) ; // 3.1
00151                         break ;
00152                     }
00153                 case  4 :
00154                     {
00155                         //if( TestInterior( test4[_config]) ) AddTriangles( tiling4_1[_config], 2 ) ; // 4.1.1
00156                         if( TestInterior( MCLookUpTable::Test4(_config) ) )
00157                             AddTriangles( MCLookUpTable::Tiling4_1(_config), 2 ) ; // 4.1.1
00158                         else
00159                             AddTriangles( MCLookUpTable::Tiling4_2(_config), 6 ) ; // 4.1.2
00160                         break ;
00161                     }
00162                 case  5 : { AddTriangles( MCLookUpTable::Tiling5(_config), 3 ); break; }
00163                 case  6 :
00164                     {
00165                         //if( TestFace( test6[_config][0]) )
00166                         if( TestFace( MCLookUpTable::Test6(_config, 0)) )
00167                             AddTriangles( MCLookUpTable::Tiling6_2(_config), 5 ) ; // 6.2
00168                         else
00169                         {
00170                             if( TestInterior( MCLookUpTable::Test6(_config, 1)) )
00171                                 AddTriangles( MCLookUpTable::Tiling6_1_1(_config), 3 ) ; // 6.1.1
00172                             else
00173                                 AddTriangles( MCLookUpTable::Tiling6_1_2(_config), 7 ) ; // 6.1.2
00174                         }
00175                         break ;
00176                     }
00177                 case  7 :
00178                     {
00179                         //if( TestFace( test7[_config][0] ) ) _subconfig +=  1 ;
00180                         //if( TestFace( test7[_config][1] ) ) _subconfig +=  2 ;
00181                         //if( TestFace( test7[_config][2] ) ) _subconfig +=  4 ;
00182                         if( TestFace( MCLookUpTable::Test7(_config, 0) ) ) _subconfig +=  1 ;
00183                         if( TestFace( MCLookUpTable::Test7(_config, 1) ) ) _subconfig +=  2 ;
00184                         if( TestFace( MCLookUpTable::Test7(_config, 2) ) ) _subconfig +=  4 ;
00185                         switch( _subconfig )
00186                         {
00187                         case 0 : { AddTriangles( MCLookUpTable::Tiling7_1(_config),             3 ) ; break; }
00188                         case 1 : { AddTriangles( MCLookUpTable::Tiling7_2(_config,0), 5 ) ; break; }
00189                         case 2 : { AddTriangles( MCLookUpTable::Tiling7_2(_config,1), 5 ) ; break; }
00190                         case 3 : { ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling7_3(_config,0), 9, v12 ) ; break ; }
00191                         case 4 : { AddTriangles( MCLookUpTable::Tiling7_2(_config, 2), 5 ) ; break ;}
00192                         case 5 : { ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling7_3(_config,1), 9, v12 ) ; break ; }
00193                         case 6 : { ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling7_3(_config,2), 9, v12 ) ; break ; }
00194                         case 7 :
00195                             {
00196                                 if( TestInterior( MCLookUpTable::Test7(_config, 3) ) )
00197                                     AddTriangles( MCLookUpTable::Tiling7_4_2(_config), 9 ) ;
00198                                 else
00199                                     AddTriangles( MCLookUpTable::Tiling7_4_1(_config), 5 ) ;
00200                                 break ;
00201                             }
00202                         };
00203                         break ;
00204                     } // end of case 7
00205                 case  8 : { AddTriangles( MCLookUpTable::Tiling8(_config), 2 ) ; break ;}
00206                 case  9 : { AddTriangles( MCLookUpTable::Tiling9(_config), 4 ) ; break ;}
00207                 case 10 :
00208                     {
00209                         if( TestFace( MCLookUpTable::Test10(_config, 0)) ) //if( TestFace( test10[_config][0]) )
00210                         {
00211                             if( TestFace( MCLookUpTable::Test10(_config,1) ) )
00212                                 AddTriangles( MCLookUpTable::Tiling10_1_1_(_config), 4 ) ; // 10.1.1
00213                             else
00214                             {
00215                                 ComputeCVertex(v12);
00216                                 AddTriangles( MCLookUpTable::Tiling10_2(_config), 8, v12 ) ; // 10.2
00217                             }
00218                         }
00219                         else
00220                         {
00221                             if( TestFace( MCLookUpTable::Test10(_config, 1) ) )
00222                             {
00223                                 ComputeCVertex(v12) ;
00224                                 AddTriangles( MCLookUpTable::Tiling10_2_(_config), 8, v12 ) ; // 10.2
00225                             }
00226                             else
00227                             {
00228                                 if( TestInterior( MCLookUpTable::Test10(_config, 2) ) )
00229                                     AddTriangles( MCLookUpTable::Tiling10_1_1(_config), 4 ) ; // 10.1.1
00230                                 else
00231                                     AddTriangles( MCLookUpTable::Tiling10_1_2(_config), 8 ) ; // 10.1.2
00232                             }
00233                         }
00234                         break ;
00235                     } // end of case 10
00236                 case 11 : { AddTriangles( MCLookUpTable::Tiling11(_config), 4 ) ; break ; }
00237                 case 12 :
00238                     {
00239                         if( TestFace( MCLookUpTable::Test12(_config, 0) ) ) //if( TestFace( test12[_config][0]) )
00240                         {
00241                             if( TestFace( MCLookUpTable::Test12(_config, 1) ) )
00242                                 AddTriangles( MCLookUpTable::Tiling12_1_1_(_config), 4 ) ; // 12.1.1
00243                             else
00244                             {
00245                                 ComputeCVertex(v12) ;
00246                                 AddTriangles( MCLookUpTable::Tiling12_2(_config), 8, v12 ) ; // 12.2
00247                             }
00248                         }
00249                         else
00250                         {
00251                             if( TestFace( MCLookUpTable::Test12(_config, 1) ) )
00252                             {
00253                                 ComputeCVertex(v12) ;
00254                                 AddTriangles( MCLookUpTable::Tiling12_2_(_config), 8, v12 ) ; // 12.2
00255                             }
00256                             else
00257                             {
00258                                 if( TestInterior( MCLookUpTable::Test12(_config, 2) ) )
00259                                     AddTriangles( MCLookUpTable::Tiling12_1_1(_config), 4 ) ; // 12.1.1
00260                                 else
00261                                     AddTriangles( MCLookUpTable::Tiling12_1_2(_config), 8 ) ; // 12.1.2
00262                             }
00263                         }
00264                         break ;
00265                     } // end of case 12
00266                 case 13 :
00267                     {
00268                         //if( TestFace( test13[_config][0] ) ) _subconfig +=  1 ;
00269                         //if( TestFace( test13[_config][1] ) ) _subconfig +=  2 ;
00270                         //if( TestFace( test13[_config][2] ) ) _subconfig +=  4 ;
00271                         //if( TestFace( test13[_config][3] ) ) _subconfig +=  8 ;
00272                         //if( TestFace( test13[_config][4] ) ) _subconfig += 16 ;
00273                         //if( TestFace( test13[_config][5] ) ) _subconfig += 32 ;
00274                         if( TestFace( MCLookUpTable::Test13(_config, 0) ) ) _subconfig +=  1 ;
00275                         if( TestFace( MCLookUpTable::Test13(_config, 1) ) ) _subconfig +=  2 ;
00276                         if( TestFace( MCLookUpTable::Test13(_config, 2) ) ) _subconfig +=  4 ;
00277                         if( TestFace( MCLookUpTable::Test13(_config, 3) ) ) _subconfig +=  8 ;
00278                         if( TestFace( MCLookUpTable::Test13(_config, 4) ) ) _subconfig += 16 ;
00279                         if( TestFace( MCLookUpTable::Test13(_config, 5) ) ) _subconfig += 32 ;
00280                         switch( MCLookUpTable::Subconfig13(_subconfig) ) //switch( subconfig13[_subconfig] )
00281                         {
00282                         case 0  : { /* 13.1 */ AddTriangles( MCLookUpTable::Tiling13_1(_config)  , 4 ) ; break ; }
00283                         case 1  : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 0), 6 ) ; break ; }
00284                         case 2  : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 1), 6 ) ; break ; }
00285                         case 3  : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 2), 6 ) ; break ; }
00286                         case 4  : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 3), 6 ) ; break ; }
00287                         case 5  : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 4), 6 ) ; break ; }
00288                         case 6  : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2(_config, 5), 6 ) ; break ; }
00289                         case 7  : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 0), 10, v12 ) ; break ; }
00290                         case 8  : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 1), 10, v12 ) ; break ; }
00291                         case 9  : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 2), 10, v12 ) ; break ; }
00292                         case 10 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 3), 10, v12 ) ; break ; }
00293                         case 11 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 4), 10, v12 ) ; break ; }
00294                         case 12 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 5), 10, v12 ) ; break ; }
00295                         case 13 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 6), 10, v12 ) ; break ; }
00296                         case 14 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 7), 10, v12 ) ; break ; }
00297                         case 15 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 8), 10, v12 ) ; break ; }
00298                         case 16 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config, 9), 10, v12 ) ; break ; }
00299                         case 17 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config,10), 10, v12 ) ; break ; }
00300                         case 18 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3(_config,11), 10, v12 ) ; break ; }
00301                         case 19 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 0), 12, v12 ) ; break ; }
00302                         case 20 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 1), 12, v12 ) ; break ; }
00303                         case 21 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 2), 12, v12 ) ; break ; }
00304                         case 22 : { /* 13.4 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_4(_config, 3), 12, v12 ) ; break ; }
00305                         case 23 :
00306                             { /* 13.5 */
00307                                 _subconfig = 0 ;
00308                                 if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
00309                                     AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 0), 6 ) ;
00310                                 else
00311                                     AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 0), 10 ) ;
00312                                 break ;
00313                             }
00314                         case 24 :
00315                             { /* 13.5 */
00316                                 _subconfig = 1 ;
00317                                 if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
00318                                     AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 1), 6 ) ;
00319                                 else
00320                                     AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 1), 10 ) ;
00321                                 break ;
00322                             }
00323                         case 25 :
00324                             {/* 13.5 */
00325                                 _subconfig = 2 ;
00326                                 if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
00327                                     AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 2), 6 ) ;
00328                                 else
00329                                     AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 2), 10 ) ;
00330                                 break ;
00331                             }
00332                         case 26 :
00333                             {/* 13.5 */
00334                                 _subconfig = 3 ;
00335                                 if( TestInterior( MCLookUpTable::Test13(_config, 6) ) )
00336                                     AddTriangles( MCLookUpTable::Tiling13_5_1(_config, 3), 6 ) ;
00337                                 else
00338                                     AddTriangles( MCLookUpTable::Tiling13_5_2(_config, 3), 10 ) ;
00339                                 break ;
00340                             }
00341                         case 27 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 0), 10, v12 ) ; break ; }
00342                         case 28 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 1), 10, v12 ) ; break ; }
00343                         case 29 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 2), 10, v12 ) ; break ; }
00344                         case 30 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 3), 10, v12 ) ; break ; }
00345                         case 31 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 4), 10, v12 ) ; break ; }
00346                         case 32 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 5), 10, v12 ) ; break ; }
00347                         case 33 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 6), 10, v12 ) ; break ; }
00348                         case 34 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 7), 10, v12 ) ; break ; }
00349                         case 35 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 8), 10, v12 ) ; break ; }
00350                         case 36 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config, 9), 10, v12 ) ; break ; }
00351                         case 37 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config,10), 10, v12 ) ; break ; }
00352                         case 38 : { /* 13.3 */ ComputeCVertex(v12); AddTriangles( MCLookUpTable::Tiling13_3_(_config,11), 10, v12 ) ; break ; }
00353                         case 39 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,0), 6 ) ; break ; }
00354                         case 40 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,1), 6 ) ; break ; }
00355                         case 41 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,2), 6 ) ; break ; }
00356                         case 42 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,3), 6 ) ; break ; }
00357                         case 43 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,4), 6 ) ; break ; }
00358                         case 44 : { /* 13.2 */ AddTriangles( MCLookUpTable::Tiling13_2_(_config,5), 6 ) ; break ; }
00359                         case 45 : { /* 13.1 */ AddTriangles( MCLookUpTable::Tiling13_1_(_config)        , 4 ) ; break ; }
00360                         default : { /*Impossible case 13*/  assert(false); }
00361                         }
00362                         break ;
00363                     } // end of case 13
00364 
00365                 case 14 : { AddTriangles( MCLookUpTable::Tiling14(_config), 4 ) ; }
00366                                     break ;
00367                 } //end of switch (_case)
00368 
00369             }; // end of ApplyMC
00370 
00371         private:
00374             WALKER_TYPE         *_walker;
00377             TRIMESH_TYPE        *_mesh;
00378 
00382             ScalarType _field[8];
00383 
00387             vcg::Point3i _corners[8];
00388 
00392             unsigned char _case;
00393 
00397             unsigned char _config;
00398 
00402             unsigned char _subconfig;
00403 
00408             inline bool TestFace(signed char face)
00409             {
00410                 ScalarType A,B,C,D ;
00411 
00412                 switch( face )
00413                 {
00414                 case -1 : case 1 :  A = _field[0] ;  B = _field[4] ;  C = _field[5] ;  D = _field[1] ;  break ;
00415                 case -2 : case 2 :  A = _field[1] ;  B = _field[5] ;  C = _field[6] ;  D = _field[2] ;  break ;
00416                 case -3 : case 3 :  A = _field[2] ;  B = _field[6] ;  C = _field[7] ;  D = _field[3] ;  break ;
00417                 case -4 : case 4 :  A = _field[3] ;  B = _field[7] ;  C = _field[4] ;  D = _field[0] ;  break ;
00418                 case -5 : case 5 :  A = _field[0] ;  B = _field[3] ;  C = _field[2] ;  D = _field[1] ;  break ;
00419                 case -6 : case 6 :  A = _field[4] ;  B = _field[7] ;  C = _field[6] ;  D = _field[5] ;  break ;
00420                 default : assert(false); // Invalid face code
00421                 };
00422 
00423                 return face * A * ( A*C - B*D ) >= 0  ;  // face and A invert signs
00424             } // end of TestFace
00425 
00426 
00431             inline bool TestInterior(signed char s)
00432             {
00433                 ScalarType t, At=0, Bt=0, Ct=0, Dt=0, a, b ;
00434                 char  test =  0 ;
00435                 char  edge = -1 ; // reference edge of the triangulation
00436 
00437                 switch( _case )
00438                 {
00439                 case  4 :
00440                 case 10 :
00441                     {
00442                         a = (_field[4]-_field[0])*(_field[6]-_field[2]) - (_field[7]-_field[3])*(_field[5]-_field[1]);
00443                         b =  _field[2]*(_field[4]-_field[0])+_field[0]*(_field[6]-_field[2])-_field[1]*(_field[7]-_field[3])-_field[3]*(_field[5]-_field[1]);
00444                         t = - b / (2*a) ;
00445                         if( t<0 || t>1 )
00446                             return s>0 ;
00447 
00448                         At = _field[0] + ( _field[4] - _field[0] ) * t ;
00449                         Bt = _field[3] + ( _field[7] - _field[3] ) * t ;
00450                         Ct = _field[2] + ( _field[6] - _field[2] ) * t ;
00451                         Dt = _field[1] + ( _field[5] - _field[1] ) * t ;
00452                         break ;
00453                     }
00454                 case  6 :
00455                 case  7 :
00456                 case 12 :
00457                 case 13 :
00458                     switch( _case )
00459                     {
00460                     case  6 : edge = MCLookUpTable::Test6 (_config, 2) ; break ;
00461                     case  7 : edge = MCLookUpTable::Test7 (_config, 4) ; break ;
00462                     case 12 : edge = MCLookUpTable::Test12(_config, 3) ; break ;
00463                     case 13 : edge = MCLookUpTable::Tiling13_5_1(_config, _subconfig)[0] ; break ;
00464                     }
00465                     switch( edge )
00466                     {
00467                     case  0 :
00468                         t  = _field[0] / ( _field[0] - _field[1] ) ;
00469                         At = 0 ;
00470                         Bt = _field[3] + ( _field[2] - _field[3] ) * t ;
00471                         Ct = _field[7] + ( _field[6] - _field[7] ) * t ;
00472                         Dt = _field[4] + ( _field[5] - _field[4] ) * t ;
00473                         break ;
00474                     case  1 :
00475                         t  = _field[1] / ( _field[1] - _field[2] ) ;
00476                         At = 0 ;
00477                         Bt = _field[0] + ( _field[3] - _field[0] ) * t ;
00478                         Ct = _field[4] + ( _field[7] - _field[4] ) * t ;
00479                         Dt = _field[5] + ( _field[6] - _field[5] ) * t ;
00480                         break ;
00481                     case  2 :
00482                         t  = _field[2] / ( _field[2] - _field[3] ) ;
00483                         At = 0 ;
00484                         Bt = _field[1] + ( _field[0] - _field[1] ) * t ;
00485                         Ct = _field[5] + ( _field[4] - _field[5] ) * t ;
00486                         Dt = _field[6] + ( _field[7] - _field[6] ) * t ;
00487                         break ;
00488                     case  3 :
00489                         t  = _field[3] / ( _field[3] - _field[0] ) ;
00490                         At = 0 ;
00491                         Bt = _field[2] + ( _field[1] - _field[2] ) * t ;
00492                         Ct = _field[6] + ( _field[5] - _field[6] ) * t ;
00493                         Dt = _field[7] + ( _field[4] - _field[7] ) * t ;
00494                         break ;
00495                     case  4 :
00496                         t  = _field[4] / ( _field[4] - _field[5] ) ;
00497                         At = 0 ;
00498                         Bt = _field[7] + ( _field[6] - _field[7] ) * t ;
00499                         Ct = _field[3] + ( _field[2] - _field[3] ) * t ;
00500                         Dt = _field[0] + ( _field[1] - _field[0] ) * t ;
00501                         break ;
00502                     case  5 :
00503                         t  = _field[5] / ( _field[5] - _field[6] ) ;
00504                         At = 0 ;
00505                         Bt = _field[4] + ( _field[7] - _field[4] ) * t ;
00506                         Ct = _field[0] + ( _field[3] - _field[0] ) * t ;
00507                         Dt = _field[1] + ( _field[2] - _field[1] ) * t ;
00508                         break ;
00509                     case  6 :
00510                         t  = _field[6] / ( _field[6] - _field[7] ) ;
00511                         At = 0 ;
00512                         Bt = _field[5] + ( _field[4] - _field[5] ) * t ;
00513                         Ct = _field[1] + ( _field[0] - _field[1] ) * t ;
00514                         Dt = _field[2] + ( _field[3] - _field[2] ) * t ;
00515                         break ;
00516                     case  7 :
00517                         t  = _field[7] / ( _field[7] - _field[4] ) ;
00518                         At = 0 ;
00519                         Bt = _field[6] + ( _field[5] - _field[6] ) * t ;
00520                         Ct = _field[2] + ( _field[1] - _field[2] ) * t ;
00521                         Dt = _field[3] + ( _field[0] - _field[3] ) * t ;
00522                         break ;
00523                     case  8 :
00524                         t  = _field[0] / ( _field[0] - _field[4] ) ;
00525                         At = 0 ;
00526                         Bt = _field[3] + ( _field[7] - _field[3] ) * t ;
00527                         Ct = _field[2] + ( _field[6] - _field[2] ) * t ;
00528                         Dt = _field[1] + ( _field[5] - _field[1] ) * t ;
00529                         break ;
00530                     case  9 :
00531                         t  = _field[1] / ( _field[1] - _field[5] ) ;
00532                         At = 0 ;
00533                         Bt = _field[0] + ( _field[4] - _field[0] ) * t ;
00534                         Ct = _field[3] + ( _field[7] - _field[3] ) * t ;
00535                         Dt = _field[2] + ( _field[6] - _field[2] ) * t ;
00536                         break ;
00537                     case 10 :
00538                         t  = _field[2] / ( _field[2] - _field[6] ) ;
00539                         At = 0 ;
00540                         Bt = _field[1] + ( _field[5] - _field[1] ) * t ;
00541                         Ct = _field[0] + ( _field[4] - _field[0] ) * t ;
00542                         Dt = _field[3] + ( _field[7] - _field[3] ) * t ;
00543                         break ;
00544                     case 11 :
00545                         t  = _field[3] / ( _field[3] - _field[7] ) ;
00546                         At = 0 ;
00547                         Bt = _field[2] + ( _field[6] - _field[2] ) * t ;
00548                         Ct = _field[1] + ( _field[5] - _field[1] ) * t ;
00549                         Dt = _field[0] + ( _field[4] - _field[0] ) * t ;
00550                         break ;
00551                     default: { assert(false); /* Invalid edge */ break ; }
00552                     }
00553                     break ;
00554 
00555                 default : assert(false); /* Invalid ambiguous case */  break;
00556                 }
00557 
00558                 if( At >= 0 ) test ++ ;
00559                 if( Bt >= 0 ) test += 2 ;
00560                 if( Ct >= 0 ) test += 4 ;
00561                 if( Dt >= 0 ) test += 8 ;
00562                 switch( test )
00563                 {
00564                 case  0 : return s>0 ;
00565                 case  1 : return s>0 ;
00566                 case  2 : return s>0 ;
00567                 case  3 : return s>0 ;
00568                 case  4 : return s>0 ;
00569                 case  5 : if( At * Ct <  Bt * Dt ) return s>0 ; break ;
00570                 case  6 : return s>0 ;
00571                 case  7 : return s<0 ;
00572                 case  8 : return s>0 ;
00573                 case  9 : return s>0 ;
00574                 case 10 : if( At * Ct >= Bt * Dt ) return s>0 ; break ;
00575                 case 11 : return s<0 ;
00576                 case 12 : return s>0 ;
00577                 case 13 : return s<0 ;
00578                 case 14 : return s<0 ;
00579                 case 15 : return s<0 ;
00580                 }
00581                 return s<0 ;
00582             } //end of TestInterior
00583 
00588             inline void ComputeCVertex(VertexPointer &v12)
00589             {
00590                 v12 = &*AllocatorType::AddVertices(*_mesh, 1);
00591                 v12->P() = CoordType(0.0, 0.0, 0.0);
00592 
00593                 unsigned int count = 0;
00594                 VertexPointer v = NULL;
00595                 if (_walker->Exist(_corners[0], _corners[1], v) )
00596                 {
00597                     count++;
00598                     v12->P() += v->P();
00599                 }
00600                 if (_walker->Exist(_corners[1], _corners[2], v) )
00601                 {
00602                     count++;
00603                     v12->P() += v->P();
00604                 }
00605                 if (_walker->Exist(_corners[3], _corners[2], v) )
00606                 {
00607                     count++;
00608                     v12->P() += v->P();
00609                 }
00610                 if (_walker->Exist(_corners[0], _corners[3], v) )
00611                 {
00612                     count++;
00613                     v12->P() += v->P();
00614                 }
00615                 if (_walker->Exist(_corners[4], _corners[5], v) )
00616                 {
00617                     count++;
00618                     v12->P() += v->P();
00619                 }
00620                 if (_walker->Exist(_corners[5], _corners[6], v) )
00621                 {
00622                     count++;
00623                     v12->P() += v->P();
00624                 }
00625                 if (_walker->Exist(_corners[7], _corners[6], v) )
00626                 {
00627                     count++;
00628                     v12->P() += v->P();
00629                 }
00630                 if (_walker->Exist(_corners[4], _corners[7], v) )
00631                 {
00632                     count++;
00633                     v12->P() += v->P();
00634                 }
00635                 if (_walker->Exist(_corners[0], _corners[4], v) )
00636                 {
00637                     count++;
00638                     v12->P() += v->P();
00639                 }
00640                 if (_walker->Exist(_corners[1], _corners[5], v) )
00641                 {
00642                     count++;
00643                     v12->P() += v->P();
00644                 }
00645                 if (_walker->Exist(_corners[2], _corners[6], v) )
00646                 {
00647                     count++;
00648                     v12->P() += v->P();
00649                 }
00650                 if (_walker->Exist(_corners[3], _corners[7], v) )
00651                 {
00652                     count++;
00653                     v12->P() += v->P();
00654                 }
00655                 v12->P() /= (float) count;
00656             } // end of AddCVertex
00663             inline void AddTriangles(const char *vertices_list, char n, VertexPointer v12=NULL)
00664             {
00665                 VertexPointer vp                = NULL;
00666                 size_t face_idx                 = _mesh->face.size();
00667                 size_t v12_idx                  = -1;
00668                 size_t vertices_idx[3];
00669                 if (v12 != NULL) v12_idx = v12 - &_mesh->vert[0];
00670                 AllocatorType::AddFaces(*_mesh, (int) n);
00671 
00672                 for (int trig=0; trig<3*n; face_idx++ )
00673                 {
00674                     vp = NULL;
00675                     memset(vertices_idx, -1, 3*sizeof(size_t));
00676                     for (int vert=0; vert<3; vert++, trig++) //ok
00677                     {
00678 
00679                         switch ( vertices_list[trig] )
00680                         {
00681                         case  0: { _walker->GetXIntercept(_corners[0], _corners[1], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00682                         case  1: { _walker->GetYIntercept(_corners[1], _corners[2], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00683                         case  2: { _walker->GetXIntercept(_corners[3], _corners[2], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00684                         case  3: { _walker->GetYIntercept(_corners[0], _corners[3], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00685                         case  4: { _walker->GetXIntercept(_corners[4], _corners[5], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00686                         case  5: { _walker->GetYIntercept(_corners[5], _corners[6], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00687                         case  6: { _walker->GetXIntercept(_corners[7], _corners[6], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00688                         case  7: { _walker->GetYIntercept(_corners[4], _corners[7], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00689                         case  8: { _walker->GetZIntercept(_corners[0], _corners[4], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00690                         case  9: { _walker->GetZIntercept(_corners[1], _corners[5], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00691                         case 10: { _walker->GetZIntercept(_corners[2], _corners[6], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00692                         case 11: { _walker->GetZIntercept(_corners[3], _corners[7], vp); vertices_idx[vert] = tri::Index(*_mesh,vp); break; }
00693                         case 12: { assert(v12 != NULL); vertices_idx[vert] = v12_idx; break; }
00694                         default: { assert(false); /* Invalid edge identifier */ }
00695                         } // end of switch
00696 
00697                         // Note that vp can be zero if we are in case 12 and that vertices_idx is surely >0 so the following assert has to be corrected as below.
00698                         // assert((vp - &_mesh->vert[0])>=0 && vertices_idx[vert]<_mesh->vert.size());
00699                         assert(vertices_idx[vert]<_mesh->vert.size());
00700                     } // end for (int vert=0 ...)
00701 
00702                     _mesh->face[face_idx].V(0) = &_mesh->vert[vertices_idx[0]];
00703                     _mesh->face[face_idx].V(1) = &_mesh->vert[vertices_idx[1]];
00704                     _mesh->face[face_idx].V(2) = &_mesh->vert[vertices_idx[2]];
00705                 } // end for (int trig=0...)
00706             } // end of AddTriangles
00707 
00708 
00709         }; // end of class MarchingCubes
00710 
00712         //end of Doxygen documentation
00713 
00714     } // end of namespace tri
00715 } // end of namespace vcg
00716 
00717 #endif //__VCG_MARCHING_CUBES