OpenWalnut  1.4.0
WMarchingLegoAlgorithm.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 WMARCHINGLEGOALGORITHM_H
00026 #define WMARCHINGLEGOALGORITHM_H
00027 
00028 #include <vector>
00029 #include <map>
00030 
00031 #include "../math/WMatrix.h"
00032 #include "../WProgressCombiner.h"
00033 
00034 #include "core/graphicsEngine/WTriangleMesh.h"
00035 
00036 /**
00037  * A point consisting of its coordinates and ID
00038  */
00039 struct WMLPointXYZId
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, WMLPointXYZId > ID2WMLPointXYZId;
00048 
00049 /**
00050  * Encapsulated ids representing a triangle.
00051  */
00052 struct WMLTriangle
00053 {
00054     unsigned int pointID[3]; //!< The IDs of the vertices of the triangle.
00055 };
00056 
00057 typedef std::vector<WMLTriangle> WMLTriangleVECTOR;
00058 
00059 
00060 /**
00061  * Creates a non interpolated triangulation of an isosurface
00062  */
00063 class WMarchingLegoAlgorithm
00064 {
00065 public:
00066     /**
00067      * standard constructor
00068      */
00069     WMarchingLegoAlgorithm();
00070 
00071     /**
00072      * destructor
00073      */
00074     ~WMarchingLegoAlgorithm();
00075 
00076     /**
00077      * Generate the triangles for the surface on the given dataSet (inGrid, vals). The texture coordinates in the resulting mesh are relative to
00078      * the grid. This means they are NOT transformed. This ensure faster grid matrix updates in texture space.
00079      * This might be useful where texture transformation matrices are used.
00080      *
00081      * \param nbCoordsX number of vertices in X direction
00082      * \param nbCoordsY number of vertices in Y direction
00083      * \param nbCoordsZ number of vertices in Z direction
00084      * \param mat the matrix transforming the vertices from canonical space
00085      * \param vals the values at the vertices
00086      * \param isoValue The surface will run through all positions with this value.
00087      * \param mainProgress Pointer to the parent's progress reporter. Leave empty if no progress should be shown
00088      *
00089      * \return the created triangle mesh
00090      */
00091     template< typename T >
00092     boost::shared_ptr< WTriangleMesh > generateSurface( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
00093                                                         const WMatrix< double >& mat,
00094                                                         const std::vector< T >* vals,
00095                                                         double isoValue,
00096                                                         boost::shared_ptr<WProgressCombiner> mainProgress
00097                                                             = boost::shared_ptr < WProgressCombiner >() );
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 progress parent's WProgressCombiner. May be empty if no status report is requested
00111      *
00112      * \return the created triangle mesh
00113      */
00114     boost::shared_ptr< WTriangleMesh > genSurfaceOneValue( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
00115                                                            const WMatrix< double >& mat,
00116                                                            const std::vector< size_t >* vals,
00117                                                            size_t isoValue,
00118                                                            boost::shared_ptr<WProgressCombiner> progress
00119                                                                 = boost::shared_ptr < WProgressCombiner >() );
00120 
00121 protected:
00122 private:
00123     /**
00124      * adds 2 triangles for a given face of the voxel
00125      * \param x position of the voxel
00126      * \param y position of the voxel
00127      * \param z position of the voxel
00128      * \param surface which side of the voxel to paint
00129      */
00130     void addSurface( size_t x, size_t y, size_t z, size_t surface );
00131 
00132     /**
00133      * returns a vertex id for a given grid point
00134      * \param nX x position in space
00135      * \param nY y position in space
00136      * \param nZ z position in space
00137      * \return the id
00138      */
00139     size_t getVertexID( size_t nX, size_t nY, size_t nZ );
00140 
00141     unsigned int m_nCellsX;  //!< No. of cells in x direction.
00142     unsigned int m_nCellsY;  //!< No. of cells in y direction.
00143     unsigned int m_nCellsZ;  //!< No. of cells in z direction.
00144 
00145     double m_tIsoLevel;  //!< The isovalue.
00146 
00147     WMatrix< double > m_matrix; //!< The 4x4 transformation matrix for the triangle vertices.
00148 
00149     ID2WMLPointXYZId m_idToVertices;  //!< List of WPointXYZIds which form the isosurface.
00150     WMLTriangleVECTOR m_trivecTriangles;  //!< List of WMCTriangleS which form the triangulation of the isosurface.
00151 };
00152 
00153 template<typename T> boost::shared_ptr<WTriangleMesh>
00154 WMarchingLegoAlgorithm::generateSurface( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
00155                                          const WMatrix< double >& mat,
00156                                          const std::vector< T >* vals,
00157                                          double isoValue,
00158                                          boost::shared_ptr<WProgressCombiner> mainProgress )
00159 {
00160     WAssert( vals, "No value set provided." );
00161 
00162     m_idToVertices.clear();
00163     m_trivecTriangles.clear();
00164 
00165     m_nCellsX = nbCoordsX - 1;
00166     m_nCellsY = nbCoordsY - 1;
00167     m_nCellsZ = nbCoordsZ - 1;
00168 
00169     m_matrix = mat;
00170 
00171     m_tIsoLevel = isoValue;
00172 
00173     size_t nX = nbCoordsX;
00174     size_t nY = nbCoordsY;
00175 
00176     size_t nPointsInSlice = nX * nY;
00177 
00178     boost::shared_ptr< WProgress > progress;
00179     if( mainProgress )
00180     {
00181         progress = boost::shared_ptr< WProgress >( new WProgress( "Marching Cubes", m_nCellsZ ) );
00182         mainProgress->addSubProgress( progress );
00183     }
00184 
00185     // Generate isosurface.
00186     for( size_t z = 0; z < m_nCellsZ; z++ )
00187     {
00188         if( progress )
00189         {
00190             ++*progress;
00191         }
00192         for( size_t y = 0; y < m_nCellsY; y++ )
00193         {
00194             for( size_t x = 0; x < m_nCellsX; x++ )
00195             {
00196                 if( ( *vals )[ z * nPointsInSlice + y * nX + x ] < m_tIsoLevel )
00197                 {
00198                     continue;
00199                 }
00200 
00201                 if( x > 0 && ( ( *vals )[ z * nPointsInSlice + y * nX + x - 1 ] < m_tIsoLevel ) )
00202                 {
00203                     addSurface( x, y, z, 1 );
00204                 }
00205                 if( x < m_nCellsX - 1 && ( ( *vals )[ z * nPointsInSlice + y * nX + x + 1 ] < m_tIsoLevel ) )
00206                 {
00207                     addSurface( x, y, z, 2 );
00208                 }
00209 
00210                 if( y > 0 && ( ( *vals )[ z * nPointsInSlice + ( y - 1 ) * nX + x ] < m_tIsoLevel ) )
00211                 {
00212                     addSurface( x, y, z, 3 );
00213                 }
00214 
00215                 if( y < m_nCellsY - 1 && ( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + x ] < m_tIsoLevel ) )
00216                 {
00217                     addSurface( x, y, z, 4 );
00218                 }
00219 
00220                 if( z > 0 && ( ( *vals )[ ( z - 1 ) * nPointsInSlice + y * nX + x ] < m_tIsoLevel ) )
00221                 {
00222                     addSurface( x, y, z, 5 );
00223                 }
00224 
00225                 if( z < m_nCellsZ - 1 && ( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + x ] < m_tIsoLevel ) )
00226                 {
00227                     addSurface( x, y, z, 6 );
00228                 }
00229 
00230                 if( x == 0 )
00231                 {
00232                     addSurface( x, y, z, 1 );
00233                 }
00234                 if( x == m_nCellsX - 1 )
00235                 {
00236                     addSurface( x, y, z, 2 );
00237                 }
00238 
00239                 if( y == 0 )
00240                 {
00241                     addSurface( x, y, z, 3 );
00242                 }
00243 
00244                 if( y == m_nCellsY - 1 )
00245                 {
00246                     addSurface( x, y, z, 4 );
00247                 }
00248 
00249                 if( z == 0 )
00250                 {
00251                     addSurface( x, y, z, 5 );
00252                 }
00253 
00254                 if( z == m_nCellsZ - 1 )
00255                 {
00256                     addSurface( x, y, z, 6 );
00257                 }
00258             }
00259         }
00260     }
00261     unsigned int nextID = 0;
00262     boost::shared_ptr< WTriangleMesh > triMesh( new WTriangleMesh( m_idToVertices.size(), m_trivecTriangles.size() ) );
00263 
00264     // Rename vertices.
00265     ID2WMLPointXYZId::iterator mapIterator = m_idToVertices.begin();
00266     while( mapIterator != m_idToVertices.end() )
00267     {
00268         WPosition texCoord = WPosition( mapIterator->second.x / nbCoordsX,
00269                                                       mapIterator->second.y / nbCoordsY,
00270                                                       mapIterator->second.z / nbCoordsZ );
00271 
00272         // transform from grid coordinate system to world coordinates
00273         WPosition pos = WPosition( mapIterator->second.x, mapIterator->second.y, mapIterator->second.z );
00274 
00275         std::vector< double > resultPos4D( 4 );
00276         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;
00277         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;
00278         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;
00279         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;
00280 
00281         ( *mapIterator ).second.newID = nextID;
00282         triMesh->addVertex( resultPos4D[0] / resultPos4D[3],
00283                             resultPos4D[1] / resultPos4D[3],
00284                             resultPos4D[2] / resultPos4D[3] );
00285         triMesh->addTextureCoordinate( texCoord );
00286         nextID++;
00287         mapIterator++;
00288     }
00289 
00290     // Now rename triangles.
00291     WMLTriangleVECTOR::iterator vecIterator = m_trivecTriangles.begin();
00292     while( vecIterator != m_trivecTriangles.end() )
00293     {
00294         for( unsigned int i = 0; i < 3; i++ )
00295         {
00296             unsigned int newID = m_idToVertices[( *vecIterator ).pointID[i]].newID;
00297             ( *vecIterator ).pointID[i] = newID;
00298         }
00299         triMesh->addTriangle( ( *vecIterator ).pointID[0], ( *vecIterator ).pointID[1], ( *vecIterator ).pointID[2] );
00300         vecIterator++;
00301     }
00302     if( progress )
00303     {
00304         progress->finish();
00305     }
00306     return triMesh;
00307 }
00308 #endif  // WMARCHINGLEGOALGORITHM_H