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