00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
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
00038
00039 struct WPointXYZId
00040 {
00041 unsigned int newID;
00042 double x;
00043 double y;
00044 double z;
00045 };
00046
00047 typedef std::map< unsigned int, WPointXYZId > ID2WPointXYZId;
00048
00049
00050
00051
00052 struct WMCTriangle
00053 {
00054 unsigned int pointID[3];
00055 };
00056
00057 typedef std::vector<WMCTriangle> WMCTriangleVECTOR;
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085 class WGE_EXPORT WMarchingCubesAlgorithm
00086 {
00087
00088
00089
00090 friend class WMarchingCubesAlgorithmTest;
00091
00092 public:
00093
00094
00095
00096 WMarchingCubesAlgorithm();
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
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
00124
00125
00126
00127
00128
00129
00130
00131
00132
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
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 WPointXYZId interpolate( double fX1, double fY1, double fZ1, double fX2, double fY2, double fZ2, double tVal1, double tVal2 );
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164 int getEdgeID( unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo );
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 unsigned int getVertexID( unsigned int nX, unsigned int nY, unsigned int nZ );
00175
00176 unsigned int m_nCellsX;
00177 unsigned int m_nCellsY;
00178 unsigned int m_nCellsZ;
00179
00180 double m_tIsoLevel;
00181
00182 WMatrix< double > m_matrix;
00183
00184 ID2WPointXYZId m_idToVertices;
00185 WMCTriangleVECTOR m_trivecTriangles;
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
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
00225
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
00245
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
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
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
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