OpenWalnut  1.5.0dev
WMarchingCubesAlgorithm.h
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #ifndef WMARCHINGCUBESALGORITHM_H
26 #define WMARCHINGCUBESALGORITHM_H
27 
28 #include <map>
29 #include <memory>
30 #include <vector>
31 
32 #include "../WProgressCombiner.h"
33 #include "../math/WMatrix.h"
34 #include "WMarchingCubesCaseTables.h"
35 #include "core/graphicsEngine/WTriangleMesh.h"
36 
37 /**
38  * A point consisting of its coordinates and ID
39  */
41 {
42  unsigned int newID; //!< ID of the point
43  double x; //!< x coordinates of the point.
44  double y; //!< y coordinates of the point.
45  double z; //!< z coordinates of the point.
46 };
47 
48 typedef std::map< unsigned int, WPointXYZId > ID2WPointXYZId;
49 
50 /**
51  * Encapsulated ids representing a triangle.
52  */
54 {
55  unsigned int pointID[3]; //!< The IDs of the vertices of the triangle.
56 };
57 
58 typedef std::vector<WMCTriangle> WMCTriangleVECTOR;
59 
60 // -------------------------------------------------------
61 //
62 // Numbering of edges (0..B) and vertices (0..7) per cube.
63 //
64 // 5--5--6
65 // /| /|
66 // 4 | 6 | A=10
67 // / 9 / A
68 // 4--7--7 |
69 // | | | |
70 // | 1-|1--2
71 // 8 / B /
72 // | 0 | 2 B=11
73 // |/ |/
74 // 0--3--3
75 //
76 // | /
77 // z y
78 // |/
79 // 0--x--
80 //
81 // -------------------------------------------------------
82 
83 /**
84  * This class does the actual computation of marching cubes.
85  */
87 {
88 /**
89  * Only UnitTests may be friends.
90  */
92 
93 public:
94  /**
95  * Constructor needed for matrix initalization.
96  */
98 
99  /**
100  * Generate the triangles for the surface on the given dataSet (inGrid, vals). The texture coordinates in the resulting mesh are relative to
101  * the grid. This means they are NOT transformed. This ensure faster grid matrix updates in texture space.
102  * This might be useful where texture transformation matrices are used.
103  *
104  * \param nbCoordsX number of vertices in X direction
105  * \param nbCoordsY number of vertices in Y direction
106  * \param nbCoordsZ number of vertices in Z direction
107  * \param mat the matrix transforming the vertices from canonical space
108  * \param vals the values at the vertices
109  * \param isoValue The surface will run through all positions with this value.
110  * \param mainProgress progress combiner used to report our progress to
111  *
112  * \return the genereated surface
113  */
114  template< typename T >
115  std::shared_ptr< WTriangleMesh > generateSurface( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
116  const WMatrix< double >& mat,
117  const std::vector< T >* vals,
118  double isoValue,
119  std::shared_ptr< WProgressCombiner > mainProgress );
120 
121 protected:
122 private:
123  /**
124  * Calculates the intersection point id of the isosurface with an
125  * edge.
126  *
127  * \param vals the values at the vertices
128  * \param nX id of cell in x direction
129  * \param nY id of cell in y direction
130  * \param nZ id of cell in z direction
131  * \param nEdgeNo id of the edge the point that will be interpolates lies on
132  *
133  * \return intersection point id
134  */
135  template< typename T > WPointXYZId calculateIntersection( const std::vector< T >* vals,
136  unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo );
137 
138  /**
139  * Interpolates between two grid points to produce the point at which
140  * the isosurface intersects an edge.
141  *
142  * \param fX1 x coordinate of first position
143  * \param fY1 y coordinate of first position
144  * \param fZ1 z coordinate of first position
145  * \param fX2 x coordinate of second position
146  * \param fY2 y coordinate of first position
147  * \param fZ2 z coordinate of first position
148  * \param tVal1 scalar value at first position
149  * \param tVal2 scalar value at second position
150  *
151  * \return interpolated point
152  */
153  WPointXYZId interpolate( double fX1, double fY1, double fZ1, double fX2, double fY2, double fZ2, double tVal1, double tVal2 );
154 
155  /**
156  * Returns the edge ID.
157  *
158  * \param nX ID of desired cell along x axis
159  * \param nY ID of desired cell along y axis
160  * \param nZ ID of desired cell along z axis
161  * \param nEdgeNo id of edge inside cell
162  *
163  * \return The id of the edge in the large array.
164  */
165  int getEdgeID( unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo );
166 
167  /**
168  * Returns the ID of the vertex given by by the IDs along the axis
169  * \param nX ID of desired vertex along x axis
170  * \param nY ID of desired vertex along y axis
171  * \param nZ ID of desired vertex along z axis
172  *
173  * \return ID of vertex with the given coordinates
174  */
175  unsigned int getVertexID( unsigned int nX, unsigned int nY, unsigned int nZ );
176 
177  unsigned int m_nCellsX; //!< No. of cells in x direction.
178  unsigned int m_nCellsY; //!< No. of cells in y direction.
179  unsigned int m_nCellsZ; //!< No. of cells in z direction.
180 
181  double m_tIsoLevel; //!< The isovalue.
182 
183  WMatrix< double > m_matrix; //!< The 4x4 transformation matrix for the triangle vertices.
184 
185  ID2WPointXYZId m_idToVertices; //!< List of WPointXYZIds which form the isosurface.
186  WMCTriangleVECTOR m_trivecTriangles; //!< List of WMCTriangleS which form the triangulation of the isosurface.
187 };
188 
189 
190 template<typename T> std::shared_ptr<WTriangleMesh> WMarchingCubesAlgorithm::generateSurface( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
191  const WMatrix< double >& mat,
192  const std::vector< T >* vals,
193  double isoValue,
194  std::shared_ptr< WProgressCombiner > mainProgress )
195 {
196  WAssert( vals, "No value set provided." );
197 
198  m_idToVertices.clear();
199  m_trivecTriangles.clear();
200 
201  m_nCellsX = nbCoordsX - 1;
202  m_nCellsY = nbCoordsY - 1;
203  m_nCellsZ = nbCoordsZ - 1;
204 
205  m_matrix = mat;
206 
207  m_tIsoLevel = isoValue;
208 
209  unsigned int nX = m_nCellsX + 1;
210  unsigned int nY = m_nCellsY + 1;
211 
212 
213  unsigned int nPointsInSlice = nX * nY;
214 
215  std::shared_ptr< WProgress > progress( new WProgress( "Marching Cubes", m_nCellsZ ) );
216  mainProgress->addSubProgress( progress );
217  // Generate isosurface.
218  for( unsigned int z = 0; z < m_nCellsZ; z++ )
219  {
220  ++*progress;
221  for( unsigned int y = 0; y < m_nCellsY; y++ )
222  {
223  for( unsigned int x = 0; x < m_nCellsX; x++ )
224  {
225  // Calculate table lookup index from those
226  // vertices which are below the isolevel.
227  unsigned int tableIndex = 0;
228  if( ( *vals )[ z * nPointsInSlice + y * nX + x ] < m_tIsoLevel )
229  tableIndex |= 1;
230  if( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + x ] < m_tIsoLevel )
231  tableIndex |= 2;
232  if( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + ( x + 1 ) ] < m_tIsoLevel )
233  tableIndex |= 4;
234  if( ( *vals )[ z * nPointsInSlice + y * nX + ( x + 1 ) ] < m_tIsoLevel )
235  tableIndex |= 8;
236  if( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + x ] < m_tIsoLevel )
237  tableIndex |= 16;
238  if( ( *vals )[ ( z + 1 ) * nPointsInSlice + ( y + 1 ) * nX + x ] < m_tIsoLevel )
239  tableIndex |= 32;
240  if( ( *vals )[ ( z + 1 ) * nPointsInSlice + ( y + 1 ) * nX + ( x + 1 ) ] < m_tIsoLevel )
241  tableIndex |= 64;
242  if( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + ( x + 1 ) ] < m_tIsoLevel )
243  tableIndex |= 128;
244 
245  // Now create a triangulation of the isosurface in this
246  // cell.
247  if( wMarchingCubesCaseTables::edgeTable[tableIndex] != 0 )
248  {
249  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 8 )
250  {
251  WPointXYZId pt = calculateIntersection( vals, x, y, z, 3 );
252  unsigned int id = getEdgeID( x, y, z, 3 );
253  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
254  }
255  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 1 )
256  {
257  WPointXYZId pt = calculateIntersection( vals, x, y, z, 0 );
258  unsigned int id = getEdgeID( x, y, z, 0 );
259  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
260  }
261  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 256 )
262  {
263  WPointXYZId pt = calculateIntersection( vals, x, y, z, 8 );
264  unsigned int id = getEdgeID( x, y, z, 8 );
265  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
266  }
267 
268  if( x == m_nCellsX - 1 )
269  {
270  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 4 )
271  {
272  WPointXYZId pt = calculateIntersection( vals, x, y, z, 2 );
273  unsigned int id = getEdgeID( x, y, z, 2 );
274  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
275  }
276  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 2048 )
277  {
278  WPointXYZId pt = calculateIntersection( vals, x, y, z, 11 );
279  unsigned int id = getEdgeID( x, y, z, 11 );
280  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
281  }
282  }
283  if( y == m_nCellsY - 1 )
284  {
285  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 2 )
286  {
287  WPointXYZId pt = calculateIntersection( vals, x, y, z, 1 );
288  unsigned int id = getEdgeID( x, y, z, 1 );
289  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
290  }
291  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 512 )
292  {
293  WPointXYZId pt = calculateIntersection( vals, x, y, z, 9 );
294  unsigned int id = getEdgeID( x, y, z, 9 );
295  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
296  }
297  }
298  if( z == m_nCellsZ - 1 )
299  {
300  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 16 )
301  {
302  WPointXYZId pt = calculateIntersection( vals, x, y, z, 4 );
303  unsigned int id = getEdgeID( x, y, z, 4 );
304  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
305  }
306  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 128 )
307  {
308  WPointXYZId pt = calculateIntersection( vals, x, y, z, 7 );
309  unsigned int id = getEdgeID( x, y, z, 7 );
310  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
311  }
312  }
313  if( ( x == m_nCellsX - 1 ) && ( y == m_nCellsY - 1 ) )
314  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 1024 )
315  {
316  WPointXYZId pt = calculateIntersection( vals, x, y, z, 10 );
317  unsigned int id = getEdgeID( x, y, z, 10 );
318  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
319  }
320  if( ( x == m_nCellsX - 1 ) && ( z == m_nCellsZ - 1 ) )
321  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 64 )
322  {
323  WPointXYZId pt = calculateIntersection( vals, x, y, z, 6 );
324  unsigned int id = getEdgeID( x, y, z, 6 );
325  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
326  }
327  if( ( y == m_nCellsY - 1 ) && ( z == m_nCellsZ - 1 ) )
328  if( wMarchingCubesCaseTables::edgeTable[tableIndex] & 32 )
329  {
330  WPointXYZId pt = calculateIntersection( vals, x, y, z, 5 );
331  unsigned int id = getEdgeID( x, y, z, 5 );
332  m_idToVertices.insert( ID2WPointXYZId::value_type( id, pt ) );
333  }
334 
335  for( int i = 0; wMarchingCubesCaseTables::triTable[tableIndex][i] != -1; i += 3 )
336  {
337  WMCTriangle triangle;
338  unsigned int pointID0, pointID1, pointID2;
339  pointID0 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i] );
340  pointID1 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i + 1] );
341  pointID2 = getEdgeID( x, y, z, wMarchingCubesCaseTables::triTable[tableIndex][i + 2] );
342  triangle.pointID[0] = pointID0;
343  triangle.pointID[1] = pointID1;
344  triangle.pointID[2] = pointID2;
345  m_trivecTriangles.push_back( triangle );
346  }
347  }
348  }
349  }
350  }
351 
352  unsigned int nextID = 0;
353  std::shared_ptr< WTriangleMesh > triMesh( new WTriangleMesh( m_idToVertices.size(), m_trivecTriangles.size() ) );
354 
355  // Rename vertices.
356  ID2WPointXYZId::iterator mapIterator = m_idToVertices.begin();
357  while( mapIterator != m_idToVertices.end() )
358  {
359  WPosition texCoord = WPosition( mapIterator->second.x / nbCoordsX,
360  mapIterator->second.y / nbCoordsY,
361  mapIterator->second.z / nbCoordsZ );
362 
363  // transform from grid coordinate system to world coordinates
364  WPosition pos = WPosition( mapIterator->second.x, mapIterator->second.y, mapIterator->second.z );
365 
366  std::vector< double > resultPos4D( 4 );
367  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;
368  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;
369  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;
370  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;
371 
372  ( *mapIterator ).second.newID = nextID;
373  triMesh->addVertex( resultPos4D[0] / resultPos4D[3], resultPos4D[1] / resultPos4D[3], resultPos4D[2] / resultPos4D[3] );
374  triMesh->addTextureCoordinate( texCoord );
375  nextID++;
376  mapIterator++;
377  }
378 
379  // Now rename triangles.
380  WMCTriangleVECTOR::iterator vecIterator = m_trivecTriangles.begin();
381  while( vecIterator != m_trivecTriangles.end() )
382  {
383  for( unsigned int i = 0; i < 3; i++ )
384  {
385  unsigned int newID = m_idToVertices[( *vecIterator ).pointID[i]].newID;
386  ( *vecIterator ).pointID[i] = newID;
387  }
388  triMesh->addTriangle( ( *vecIterator ).pointID[0], ( *vecIterator ).pointID[1], ( *vecIterator ).pointID[2] );
389  vecIterator++;
390  }
391 
392  progress->finish();
393  return triMesh;
394 }
395 
396 template< typename T > WPointXYZId WMarchingCubesAlgorithm::calculateIntersection( const std::vector< T >* vals,
397  unsigned int nX, unsigned int nY, unsigned int nZ,
398  unsigned int nEdgeNo )
399 {
400  double x1;
401  double y1;
402  double z1;
403 
404  double x2;
405  double y2;
406  double z2;
407 
408  unsigned int v1x = nX;
409  unsigned int v1y = nY;
410  unsigned int v1z = nZ;
411 
412  unsigned int v2x = nX;
413  unsigned int v2y = nY;
414  unsigned int v2z = nZ;
415 
416  switch( nEdgeNo )
417  {
418  case 0:
419  v2y += 1;
420  break;
421  case 1:
422  v1y += 1;
423  v2x += 1;
424  v2y += 1;
425  break;
426  case 2:
427  v1x += 1;
428  v1y += 1;
429  v2x += 1;
430  break;
431  case 3:
432  v1x += 1;
433  break;
434  case 4:
435  v1z += 1;
436  v2y += 1;
437  v2z += 1;
438  break;
439  case 5:
440  v1y += 1;
441  v1z += 1;
442  v2x += 1;
443  v2y += 1;
444  v2z += 1;
445  break;
446  case 6:
447  v1x += 1;
448  v1y += 1;
449  v1z += 1;
450  v2x += 1;
451  v2z += 1;
452  break;
453  case 7:
454  v1x += 1;
455  v1z += 1;
456  v2z += 1;
457  break;
458  case 8:
459  v2z += 1;
460  break;
461  case 9:
462  v1y += 1;
463  v2y += 1;
464  v2z += 1;
465  break;
466  case 10:
467  v1x += 1;
468  v1y += 1;
469  v2x += 1;
470  v2y += 1;
471  v2z += 1;
472  break;
473  case 11:
474  v1x += 1;
475  v2x += 1;
476  v2z += 1;
477  break;
478  }
479 
480  x1 = v1x;
481  y1 = v1y;
482  z1 = v1z;
483  x2 = v2x;
484  y2 = v2y;
485  z2 = v2z;
486 
487  unsigned int nPointsInSlice = ( m_nCellsX + 1 ) * ( m_nCellsY + 1 );
488  double val1 = ( *vals )[ v1z * nPointsInSlice + v1y * ( m_nCellsX + 1 ) + v1x ];
489  double val2 = ( *vals )[ v2z * nPointsInSlice + v2y * ( m_nCellsX + 1 ) + v2x ];
490 
491  WPointXYZId intersection = interpolate( x1, y1, z1, x2, y2, z2, val1, val2 );
492  intersection.newID = 0;
493  return intersection;
494 }
495 
496 #endif // WMARCHINGCUBESALGORITHM_H
Tests for the class computing the actual marching cubes.
This class does the actual computation of marching cubes.
unsigned int m_nCellsZ
No. of cells in z direction.
WPointXYZId interpolate(double fX1, double fY1, double fZ1, double fX2, double fY2, double fZ2, double tVal1, double tVal2)
Interpolates between two grid points to produce the point at which the isosurface intersects an edge.
WMCTriangleVECTOR m_trivecTriangles
List of WMCTriangleS which form the triangulation of the isosurface.
unsigned int getVertexID(unsigned int nX, unsigned int nY, unsigned int nZ)
Returns the ID of the vertex given by by the IDs along the axis.
WMatrix< double > m_matrix
The 4x4 transformation matrix for the triangle vertices.
int getEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
Returns the edge ID.
WMarchingCubesAlgorithm()
Constructor needed for matrix initalization.
ID2WPointXYZId m_idToVertices
List of WPointXYZIds which form the isosurface.
unsigned int m_nCellsY
No. of cells in y direction.
WPointXYZId calculateIntersection(const std::vector< T > *vals, unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo)
Calculates the intersection point id of the isosurface with an edge.
std::shared_ptr< WTriangleMesh > generateSurface(size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ, const WMatrix< double > &mat, const std::vector< T > *vals, double isoValue, std::shared_ptr< WProgressCombiner > mainProgress)
Generate the triangles for the surface on the given dataSet (inGrid, vals).
unsigned int m_nCellsX
No. of cells in x direction.
This only is a 3d double vector.
Class managing progress inside of modules.
Definition: WProgress.h:42
Triangle mesh data structure allowing for convenient access of the elements.
Definition: WTriangleMesh.h:46
Encapsulated ids representing a triangle.
unsigned int pointID[3]
The IDs of the vertices of the triangle.
A point consisting of its coordinates and ID.
double z
z coordinates of the point.
double x
x coordinates of the point.
double y
y coordinates of the point.
unsigned int newID
ID of the point.