OpenWalnut  1.4.0
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 <vector>
29 #include <map>
30 
31 #include "../math/WMatrix.h"
32 #include "../WProgressCombiner.h"
33 #include "core/graphicsEngine/WTriangleMesh.h"
34 
35 #include "WMarchingCubesCaseTables.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  boost::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  boost::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> boost::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  boost::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  boost::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  boost::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