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