OpenWalnut 1.2.5
WMarchingCubesAlgorithm.h
00001 //---------------------------------------------------------------------------
00002 //
00003 // Project: OpenWalnut ( http://www.openwalnut.org )
00004 //
00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
00006 // For more information see http://www.openwalnut.org/copying
00007 //
00008 // This file is part of OpenWalnut.
00009 //
00010 // OpenWalnut is free software: you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as published by
00012 // the Free Software Foundation, either version 3 of the License, or
00013 // (at your option) any later version.
00014 //
00015 // OpenWalnut is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 // GNU Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public License
00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
00022 //
00023 //---------------------------------------------------------------------------
00024 
00025 #ifndef WMARCHINGCUBESALGORITHM_H
00026 #define WMARCHINGCUBESALGORITHM_H
00027 
00028 #include <vector>
00029 #include <map>
00030 #include "../../common/math/WMatrix.h"
00031 #include "../../common/WProgressCombiner.h"
00032 #include "../WTriangleMesh.h"
00033 #include "marchingCubesCaseTables.h"
00034 #include "../WExportWGE.h"
00035 
00036 /**
00037  * A point consisting of its coordinates and ID
00038  */
00039 struct WPointXYZId
00040 {
00041     unsigned int newID; //!< ID of the point
00042     double x; //!< x coordinates of the point.
00043     double y; //!< y coordinates of the point.
00044     double z; //!< z coordinates of the point.
00045 };
00046 
00047 typedef std::map< unsigned int, WPointXYZId > ID2WPointXYZId;
00048 
00049 /**
00050  * Encapsulated ids representing a triangle.
00051  */
00052 struct WMCTriangle
00053 {
00054     unsigned int pointID[3]; //!< The IDs of the vertices of the triangle.
00055 };
00056 
00057 typedef std::vector<WMCTriangle> WMCTriangleVECTOR;
00058 
00059 // -------------------------------------------------------
00060 //
00061 // Numbering of edges (0..B) and vertices (0..7) per cube.
00062 //
00063 //      5--5--6
00064 //     /|    /|
00065 //    4 |   6 |    A=10
00066 //   /  9  /  A
00067 //  4--7--7   |
00068 //  |   | |   |
00069 //  |   1-|1--2
00070 //  8  /  B  /
00071 //  | 0   | 2      B=11
00072 //  |/    |/
00073 //  0--3--3
00074 //
00075 //  |  /
00076 //  z y
00077 //  |/
00078 //  0--x--
00079 //
00080 // -------------------------------------------------------
00081 
00082 /**
00083  * This class does the actual computation of marching cubes.
00084  */
00085 class WGE_EXPORT WMarchingCubesAlgorithm
00086 {
00087 /**
00088  * Only UnitTests may be friends.
00089  */
00090 friend class WMarchingCubesAlgorithmTest;
00091 
00092 public:
00093     /**
00094      * Constructor needed for matrix initalization.
00095      */
00096      WMarchingCubesAlgorithm();
00097 
00098     /**
00099      * Generate the triangles for the surface on the given dataSet (inGrid, vals). The texture coordinates in the resulting mesh are relative to
00100      * the grid. This means they are NOT transformed. This ensure faster grid matrix updates in texture space.
00101      * This might be useful where texture transformation matrices are used.
00102      *
00103      * \param nbCoordsX number of vertices in X direction
00104      * \param nbCoordsY number of vertices in Y direction
00105      * \param nbCoordsZ number of vertices in Z direction
00106      * \param mat the matrix transforming the vertices from canonical space
00107      * \param vals the values at the vertices
00108      * \param isoValue The surface will run through all positions with this value.
00109      * \param mainProgress progress combiner used to report our progress to
00110      *
00111      * \return the genereated surface
00112      */
00113     template< typename T >
00114     boost::shared_ptr< WTriangleMesh > generateSurface(  size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
00115                                                           const WMatrix< double >& mat,
00116                                                           const std::vector< T >* vals,
00117                                                           double isoValue,
00118                                                           boost::shared_ptr< WProgressCombiner > mainProgress );
00119 
00120 protected:
00121 private:
00122     /**
00123      * Calculates the intersection point id of the isosurface with an
00124      * edge.
00125      *
00126      * \param vals the values at the vertices
00127      * \param nX id of cell in x direction
00128      * \param nY id of cell in y direction
00129      * \param nZ id of cell in z direction
00130      * \param nEdgeNo id of the edge the point that will be interpolates lies on
00131      *
00132      * \return intersection point id
00133      */
00134     template< typename T > WPointXYZId calculateIntersection( const std::vector< T >* vals,
00135                                                               unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo );
00136 
00137     /**
00138      * Interpolates between two grid points to produce the point at which
00139      * the isosurface intersects an edge.
00140      *
00141      * \param fX1 x coordinate of first position
00142      * \param fY1 y coordinate of first position
00143      * \param fZ1 z coordinate of first position
00144      * \param fX2 x coordinate of second position
00145      * \param fY2 y coordinate of first position
00146      * \param fZ2 z coordinate of first position
00147      * \param tVal1 scalar value at first position
00148      * \param tVal2 scalar value at second position
00149      *
00150      * \return interpolated point
00151      */
00152     WPointXYZId interpolate( double fX1, double fY1, double fZ1, double fX2, double fY2, double fZ2, double tVal1, double tVal2 );
00153 
00154     /**
00155      * Returns the edge ID.
00156      *
00157      * \param nX ID of desired cell along x axis
00158      * \param nY ID of desired cell along y axis
00159      * \param nZ ID of desired cell along z axis
00160      * \param nEdgeNo id of edge inside cell
00161      *
00162      * \return The id of the edge in the large array.
00163      */
00164     int getEdgeID( unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo );
00165 
00166     /**
00167      * Returns the ID of the vertex given by by the IDs along the axis
00168      * \param nX ID of desired vertex along x axis
00169      * \param nY ID of desired vertex along y axis
00170      * \param nZ ID of desired vertex along z axis
00171      *
00172      * \return ID of vertex with the given coordinates
00173      */
00174     unsigned int getVertexID( unsigned int nX, unsigned int nY, unsigned int nZ );
00175 
00176     unsigned int m_nCellsX;  //!< No. of cells in x direction.
00177     unsigned int m_nCellsY;  //!< No. of cells in y direction.
00178     unsigned int m_nCellsZ;  //!< No. of cells in z direction.
00179 
00180     double m_tIsoLevel;  //!< The isovalue.
00181 
00182     WMatrix< double > m_matrix; //!< The 4x4 transformation matrix for the triangle vertices.
00183 
00184     ID2WPointXYZId m_idToVertices;  //!< List of WPointXYZIds which form the isosurface.
00185     WMCTriangleVECTOR m_trivecTriangles;  //!< List of WMCTriangleS which form the triangulation of the isosurface.
00186 };
00187 
00188 
00189 template<typename T> boost::shared_ptr<WTriangleMesh> WMarchingCubesAlgorithm::generateSurface( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
00190                                                                                                  const WMatrix< double >& mat,
00191                                                                                                  const std::vector< T >* vals,
00192                                                                                                  double isoValue,
00193                                                                                                  boost::shared_ptr< WProgressCombiner > mainProgress )
00194 {
00195     WAssert( vals, "No value set provided." );
00196 
00197     m_idToVertices.clear();
00198     m_trivecTriangles.clear();
00199 
00200     m_nCellsX = nbCoordsX - 1;
00201     m_nCellsY = nbCoordsY - 1;
00202     m_nCellsZ = nbCoordsZ - 1;
00203 
00204     m_matrix = mat;
00205 
00206     m_tIsoLevel = isoValue;
00207 
00208     unsigned int nX = m_nCellsX + 1;
00209     unsigned int nY = m_nCellsY + 1;
00210 
00211 
00212     unsigned int nPointsInSlice = nX * nY;
00213 
00214     boost::shared_ptr< WProgress > progress = boost::shared_ptr< WProgress >( new WProgress( "Marching Cubes", m_nCellsZ ) );
00215     mainProgress->addSubProgress( progress );
00216     // Generate isosurface.
00217     for( unsigned int z = 0; z < m_nCellsZ; z++ )
00218     {
00219         ++*progress;
00220         for( unsigned int y = 0; y < m_nCellsY; y++ )
00221         {
00222             for( unsigned int x = 0; x < m_nCellsX; x++ )
00223             {
00224                 // Calculate table lookup index from those
00225                 // vertices which are below the isolevel.
00226                 unsigned int tableIndex = 0;
00227                 if( ( *vals )[ z * nPointsInSlice + y * nX + x ] < m_tIsoLevel )
00228                     tableIndex |= 1;
00229                 if( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + x ] < m_tIsoLevel )
00230                     tableIndex |= 2;
00231                 if( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + ( x + 1 ) ] < m_tIsoLevel )
00232                     tableIndex |= 4;
00233                 if( ( *vals )[ z * nPointsInSlice + y * nX + ( x + 1 ) ] < m_tIsoLevel )
00234                     tableIndex |= 8;
00235                 if( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + x ] < m_tIsoLevel )
00236                     tableIndex |= 16;
00237                 if( ( *vals )[ ( z + 1 ) * nPointsInSlice + ( y + 1 ) * nX + x ] < m_tIsoLevel )
00238                     tableIndex |= 32;
00239                 if( ( *vals )[ ( z + 1 ) * nPointsInSlice + ( y + 1 ) * nX + ( x + 1 ) ] < m_tIsoLevel )
00240                     tableIndex |= 64;
00241                 if( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + ( x + 1 ) ] < m_tIsoLevel )
00242                     tableIndex |= 128;
00243 
00244                 // Now create a triangulation of the isosurface in this
00245                 // cell.
00246                 if( wMarchingCubesCaseTables::edgeTable[tableIndex] != 0 )
00247                 {
00248                     if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 8 )
00249                     {
00250                         WPointXYZId pt = calculateIntersection( vals, x, y, z, 3 );
00251                         unsigned int id = getEdgeID( x, y, z, 3 );
00252                         m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00253                     }
00254                     if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 1 )
00255                     {
00256                         WPointXYZId pt = calculateIntersection( vals, x, y, z, 0 );
00257                         unsigned int id = getEdgeID( x, y, z, 0 );
00258                         m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00259                     }
00260                     if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 256 )
00261                     {
00262                         WPointXYZId pt = calculateIntersection( vals, x, y, z, 8 );
00263                         unsigned int id = getEdgeID( x, y, z, 8 );
00264                         m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00265                     }
00266 
00267                     if( x == m_nCellsX - 1 )
00268                     {
00269                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 4 )
00270                         {
00271                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 2 );
00272                             unsigned int id = getEdgeID( x, y, z, 2 );
00273                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00274                         }
00275                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 2048 )
00276                         {
00277                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 11 );
00278                             unsigned int id = getEdgeID( x, y, z, 11 );
00279                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00280                         }
00281                     }
00282                     if( y == m_nCellsY - 1 )
00283                     {
00284                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 2 )
00285                         {
00286                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 1 );
00287                             unsigned int id = getEdgeID( x, y, z, 1 );
00288                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00289                         }
00290                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 512 )
00291                         {
00292                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 9 );
00293                             unsigned int id = getEdgeID( x, y, z, 9 );
00294                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00295                         }
00296                     }
00297                     if( z == m_nCellsZ - 1 )
00298                     {
00299                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 16 )
00300                         {
00301                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 4 );
00302                             unsigned int id = getEdgeID( x, y, z, 4 );
00303                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00304                         }
00305                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 128 )
00306                         {
00307                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 7 );
00308                             unsigned int id = getEdgeID( x, y, z, 7 );
00309                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00310                         }
00311                     }
00312                     if( ( x == m_nCellsX - 1 ) && ( y == m_nCellsY - 1 ) )
00313                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 1024 )
00314                         {
00315                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 10 );
00316                             unsigned int id = getEdgeID( x, y, z, 10 );
00317                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00318                         }
00319                     if( ( x == m_nCellsX - 1 ) && ( z == m_nCellsZ - 1 ) )
00320                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 64 )
00321                         {
00322                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 6 );
00323                             unsigned int id = getEdgeID( x, y, z, 6 );
00324                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00325                         }
00326                     if( ( y == m_nCellsY - 1 ) && ( z == m_nCellsZ - 1 ) )
00327                         if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 32 )
00328                         {
00329                             WPointXYZId pt = calculateIntersection( vals, x, y, z, 5 );
00330                             unsigned int id = getEdgeID( x, y, z, 5 );
00331                             m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
00332                         }
00333 
00334                     for( int i = 0; wMarchingCubesCaseTables::triTable[tableIndex][i] != -1; i += 3 )
00335                     {
00336                         WMCTriangle triangle;
00337                         unsigned int pointID0, pointID1, pointID2;
00338                         pointID0 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i] );
00339                         pointID1 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i + 1] );
00340                         pointID2 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i + 2] );
00341                         triangle.pointID[0] = pointID0;
00342                         triangle.pointID[1] = pointID1;
00343                         triangle.pointID[2] = pointID2;
00344                         m_trivecTriangles.push_back( triangle );
00345                     }
00346                 }
00347             }
00348         }
00349     }
00350 
00351     unsigned int nextID = 0;
00352     boost::shared_ptr< WTriangleMesh > triMesh( new WTriangleMesh( m_idToVertices.size(), m_trivecTriangles.size() ) );
00353 
00354     // Rename vertices.
00355     ID2WPointXYZId::iterator mapIterator = m_idToVertices.begin();
00356     while( mapIterator != m_idToVertices.end() )
00357     {
00358         WPosition texCoord = WPosition( mapIterator->second.x / nbCoordsX,
00359                                                       mapIterator->second.y / nbCoordsY,
00360                                                       mapIterator->second.z / nbCoordsZ );
00361 
00362         // transform from grid coordinate system to world coordinates
00363         WPosition pos = WPosition( mapIterator->second.x, mapIterator->second.y, mapIterator->second.z );
00364 
00365         std::vector< double > resultPos4D( 4 );
00366         resultPos4D[0] = m_matrix( 0, 0 ) * pos[0] + m_matrix( 0, 1 ) * pos[1] + m_matrix( 0, 2 ) * pos[2] + m_matrix( 0, 3 ) * 1;
00367         resultPos4D[1] = m_matrix( 1, 0 ) * pos[0] + m_matrix( 1, 1 ) * pos[1] + m_matrix( 1, 2 ) * pos[2] + m_matrix( 1, 3 ) * 1;
00368         resultPos4D[2] = m_matrix( 2, 0 ) * pos[0] + m_matrix( 2, 1 ) * pos[1] + m_matrix( 2, 2 ) * pos[2] + m_matrix( 2, 3 ) * 1;
00369         resultPos4D[3] = m_matrix( 3, 0 ) * pos[0] + m_matrix( 3, 1 ) * pos[1] + m_matrix( 3, 2 ) * pos[2] + m_matrix( 3, 3 ) * 1;
00370 
00371         ( *mapIterator ).second.newID = nextID;
00372         triMesh->addVertex( resultPos4D[0] / resultPos4D[3], resultPos4D[1] / resultPos4D[3], resultPos4D[2] / resultPos4D[3] );
00373         triMesh->addTextureCoordinate( texCoord );
00374         nextID++;
00375         mapIterator++;
00376     }
00377 
00378     // Now rename triangles.
00379     WMCTriangleVECTOR::iterator vecIterator = m_trivecTriangles.begin();
00380     while( vecIterator != m_trivecTriangles.end() )
00381     {
00382         for( unsigned int i = 0; i < 3; i++ )
00383         {
00384             unsigned int newID = m_idToVertices[( *vecIterator ).pointID[i]].newID;
00385             ( *vecIterator ).pointID[i] = newID;
00386         }
00387         triMesh->addTriangle( ( *vecIterator ).pointID[0], ( *vecIterator ).pointID[1], ( *vecIterator ).pointID[2] );
00388         vecIterator++;
00389     }
00390 
00391     progress->finish();
00392     return triMesh;
00393 }
00394 
00395 template< typename T > WPointXYZId WMarchingCubesAlgorithm::calculateIntersection( const std::vector< T >* vals,
00396                                                                                    unsigned int nX, unsigned int nY, unsigned int nZ,
00397                                                                                    unsigned int nEdgeNo )
00398 {
00399     double x1;
00400     double y1;
00401     double z1;
00402 
00403     double x2;
00404     double y2;
00405     double z2;
00406 
00407     unsigned int v1x = nX;
00408     unsigned int v1y = nY;
00409     unsigned int v1z = nZ;
00410 
00411     unsigned int v2x = nX;
00412     unsigned int v2y = nY;
00413     unsigned int v2z = nZ;
00414 
00415     switch ( nEdgeNo )
00416     {
00417         case 0:
00418             v2y += 1;
00419             break;
00420         case 1:
00421             v1y += 1;
00422             v2x += 1;
00423             v2y += 1;
00424             break;
00425         case 2:
00426             v1x += 1;
00427             v1y += 1;
00428             v2x += 1;
00429             break;
00430         case 3:
00431             v1x += 1;
00432             break;
00433         case 4:
00434             v1z += 1;
00435             v2y += 1;
00436             v2z += 1;
00437             break;
00438         case 5:
00439             v1y += 1;
00440             v1z += 1;
00441             v2x += 1;
00442             v2y += 1;
00443             v2z += 1;
00444             break;
00445         case 6:
00446             v1x += 1;
00447             v1y += 1;
00448             v1z += 1;
00449             v2x += 1;
00450             v2z += 1;
00451             break;
00452         case 7:
00453             v1x += 1;
00454             v1z += 1;
00455             v2z += 1;
00456             break;
00457         case 8:
00458             v2z += 1;
00459             break;
00460         case 9:
00461             v1y += 1;
00462             v2y += 1;
00463             v2z += 1;
00464             break;
00465         case 10:
00466             v1x += 1;
00467             v1y += 1;
00468             v2x += 1;
00469             v2y += 1;
00470             v2z += 1;
00471             break;
00472         case 11:
00473             v1x += 1;
00474             v2x += 1;
00475             v2z += 1;
00476             break;
00477     }
00478 
00479     x1 = v1x;
00480     y1 = v1y;
00481     z1 = v1z;
00482     x2 = v2x;
00483     y2 = v2y;
00484     z2 = v2z;
00485 
00486     unsigned int nPointsInSlice = ( m_nCellsX + 1 ) * ( m_nCellsY + 1 );
00487     double val1 = ( *vals )[ v1z * nPointsInSlice + v1y * ( m_nCellsX + 1 ) + v1x ];
00488     double val2 = ( *vals )[ v2z * nPointsInSlice + v2y * ( m_nCellsX + 1 ) + v2x ];
00489 
00490     WPointXYZId intersection = interpolate( x1, y1, z1, x2, y2, z2, val1, val2 );
00491     intersection.newID = 0;
00492     return intersection;
00493 }
00494 
00495 #endif  // WMARCHINGCUBESALGORITHM_H
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends