OpenWalnut
1.2.5
|
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