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