OpenWalnut  1.4.0
WMarchingLegoAlgorithm.cpp
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 #include <vector>
26 
27 #include "WMarchingLegoAlgorithm.h"
28 
30  : m_matrix( 4, 4 )
31 {
32 }
33 
35 {
36 }
37 
38 void WMarchingLegoAlgorithm::addSurface( size_t x, size_t y, size_t z, size_t surface )
39 {
40  WMLPointXYZId pt1;
41  WMLPointXYZId pt2;
42  WMLPointXYZId pt3;
43  WMLPointXYZId pt4;
44 
45  pt1.newID = 0;
46  pt2.newID = 0;
47  pt3.newID = 0;
48  pt4.newID = 0;
49 
50  switch( surface )
51  {
52  case 1:
53  {
54  pt1.x = x;
55  pt1.y = y;
56  pt1.z = z;
57  unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
58  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
59 
60  pt2.x = x;
61  pt2.y = y + 1;
62  pt2.z = z;
63  unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
64  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
65 
66  pt3.x = x;
67  pt3.y = y + 1;
68  pt3.z = z + 1;
69  unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
70  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
71 
72  pt4.x = x;
73  pt4.y = y;
74  pt4.z = z + 1;
75  unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
76  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
77 
78  WMLTriangle triangle1;
79  triangle1.pointID[0] = id1;
80  triangle1.pointID[1] = id2;
81  triangle1.pointID[2] = id3;
82  WMLTriangle triangle2;
83  triangle2.pointID[0] = id3;
84  triangle2.pointID[1] = id4;
85  triangle2.pointID[2] = id1;
86  m_trivecTriangles.push_back( triangle1 );
87  m_trivecTriangles.push_back( triangle2 );
88  break;
89  }
90  case 2:
91  {
92  pt1.x = x + 1;
93  pt1.y = y;
94  pt1.z = z;
95  unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
96  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
97 
98  pt2.x = x + 1;
99  pt2.y = y;
100  pt2.z = z + 1;
101  unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
102  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
103 
104  pt3.x = x + 1;
105  pt3.y = y + 1;
106  pt3.z = z + 1;
107  unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
108  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
109 
110  pt4.x = x + 1;
111  pt4.y = y + 1;
112  pt4.z = z;
113  unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
114  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
115 
116  WMLTriangle triangle1;
117  triangle1.pointID[0] = id1;
118  triangle1.pointID[1] = id2;
119  triangle1.pointID[2] = id3;
120  WMLTriangle triangle2;
121  triangle2.pointID[0] = id3;
122  triangle2.pointID[1] = id4;
123  triangle2.pointID[2] = id1;
124  m_trivecTriangles.push_back( triangle1 );
125  m_trivecTriangles.push_back( triangle2 );
126  break;
127  }
128  case 3:
129  {
130  pt1.x = x;
131  pt1.y = y;
132  pt1.z = z;
133  unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
134  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
135 
136  pt2.x = x;
137  pt2.y = y;
138  pt2.z = z + 1;
139  unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
140  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
141 
142  pt3.x = x + 1;
143  pt3.y = y;
144  pt3.z = z + 1;
145  unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
146  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
147 
148  pt4.x = x + 1;
149  pt4.y = y;
150  pt4.z = z;
151  unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
152  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
153 
154  WMLTriangle triangle1;
155  triangle1.pointID[0] = id1;
156  triangle1.pointID[1] = id2;
157  triangle1.pointID[2] = id3;
158  WMLTriangle triangle2;
159  triangle2.pointID[0] = id3;
160  triangle2.pointID[1] = id4;
161  triangle2.pointID[2] = id1;
162  m_trivecTriangles.push_back( triangle1 );
163  m_trivecTriangles.push_back( triangle2 );
164  break;
165  }
166  case 4:
167  {
168  pt1.x = x;
169  pt1.y = y + 1;
170  pt1.z = z;
171  unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
172  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
173 
174  pt2.x = x + 1;
175  pt2.y = y + 1;
176  pt2.z = z;
177  unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
178  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
179 
180  pt3.x = x + 1;
181  pt3.y = y + 1;
182  pt3.z = z + 1;
183  unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
184  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
185 
186  pt4.x = x;
187  pt4.y = y + 1;
188  pt4.z = z + 1;
189  unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
190  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
191 
192  WMLTriangle triangle1;
193  triangle1.pointID[0] = id1;
194  triangle1.pointID[1] = id2;
195  triangle1.pointID[2] = id3;
196  WMLTriangle triangle2;
197  triangle2.pointID[0] = id3;
198  triangle2.pointID[1] = id4;
199  triangle2.pointID[2] = id1;
200  m_trivecTriangles.push_back( triangle1 );
201  m_trivecTriangles.push_back( triangle2 );
202  break;
203  }
204  case 5:
205  {
206  pt1.x = x;
207  pt1.y = y;
208  pt1.z = z;
209  unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
210  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
211 
212  pt2.x = x + 1;
213  pt2.y = y;
214  pt2.z = z;
215  unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
216  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
217 
218  pt3.x = x + 1;
219  pt3.y = y + 1;
220  pt3.z = z;
221  unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
222  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
223 
224  pt4.x = x;
225  pt4.y = y + 1;
226  pt4.z = z;
227  unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
228  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
229 
230  WMLTriangle triangle1;
231  triangle1.pointID[0] = id1;
232  triangle1.pointID[1] = id2;
233  triangle1.pointID[2] = id3;
234  WMLTriangle triangle2;
235  triangle2.pointID[0] = id3;
236  triangle2.pointID[1] = id4;
237  triangle2.pointID[2] = id1;
238  m_trivecTriangles.push_back( triangle1 );
239  m_trivecTriangles.push_back( triangle2 );
240  break;
241  }
242  case 6:
243  {
244  pt1.x = x;
245  pt1.y = y;
246  pt1.z = z + 1;
247  unsigned int id1 = getVertexID( pt1.x , pt1.y, pt1.z );
248  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id1, pt1 ) );
249 
250  pt2.x = x;
251  pt2.y = y + 1;
252  pt2.z = z + 1;
253  unsigned int id2 = getVertexID( pt2.x , pt2.y, pt2.z );
254  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id2, pt2 ) );
255 
256  pt3.x = x + 1;
257  pt3.y = y + 1;
258  pt3.z = z + 1;
259  unsigned int id3 = getVertexID( pt3.x , pt3.y, pt3.z );
260  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id3, pt3 ) );
261 
262  pt4.x = x + 1;
263  pt4.y = y;
264  pt4.z = z + 1;
265  unsigned int id4 = getVertexID( pt4.x , pt4.y, pt4.z );
266  m_idToVertices.insert( ID2WMLPointXYZId::value_type( id4, pt4 ) );
267 
268  WMLTriangle triangle1;
269  triangle1.pointID[0] = id1;
270  triangle1.pointID[1] = id2;
271  triangle1.pointID[2] = id3;
272  WMLTriangle triangle2;
273  triangle2.pointID[0] = id3;
274  triangle2.pointID[1] = id4;
275  triangle2.pointID[2] = id1;
276  m_trivecTriangles.push_back( triangle1 );
277  m_trivecTriangles.push_back( triangle2 );
278  break;
279  }
280  default:
281  break;
282  }
283 }
284 
285 size_t WMarchingLegoAlgorithm::getVertexID( size_t nX, size_t nY, size_t nZ )
286 {
287  return nZ * ( m_nCellsY + 1 ) * ( m_nCellsX + 1) + nY * ( m_nCellsX + 1 ) + nX;
288 }
289 
290 boost::shared_ptr<WTriangleMesh> WMarchingLegoAlgorithm::genSurfaceOneValue( size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ,
291  const WMatrix< double >& mat,
292  const std::vector< size_t >* vals,
293  size_t isoValue,
294  boost::shared_ptr< WProgressCombiner > mainProgress )
295 {
296  WAssert( vals, "No value set provided." );
297 
298  m_idToVertices.clear();
299  m_trivecTriangles.clear();
300 
301  m_nCellsX = nbCoordsX - 1;
302  m_nCellsY = nbCoordsY - 1;
303  m_nCellsZ = nbCoordsZ - 1;
304 
305  m_matrix = mat;
306 
307  size_t nX = nbCoordsX;
308  size_t nY = nbCoordsY;
309 
310  size_t nPointsInSlice = nX * nY;
311 
312  boost::shared_ptr< WProgress > progress;
313  if( mainProgress )
314  {
315  progress = boost::shared_ptr< WProgress >( new WProgress( "Marching Legos", m_nCellsZ ) );
316  mainProgress->addSubProgress( progress );
317  }
318 
319  // Generate isosurface.
320  for( size_t z = 0; z < m_nCellsZ; z++ )
321  {
322  if( progress )
323  {
324  ++*progress;
325  }
326  for( size_t y = 0; y < m_nCellsY; y++ )
327  {
328  for( size_t x = 0; x < m_nCellsX; x++ )
329  {
330  if( ( *vals )[ z * nPointsInSlice + y * nX + x ] != isoValue )
331  {
332  continue;
333  }
334 
335  if( x > 0 && ( ( *vals )[ z * nPointsInSlice + y * nX + x - 1 ] != isoValue ) )
336  {
337  addSurface( x, y, z, 1 );
338  }
339  if( x < m_nCellsX - 1 && ( ( *vals )[ z * nPointsInSlice + y * nX + x + 1 ] != isoValue ) )
340  {
341  addSurface( x, y, z, 2 );
342  }
343 
344  if( y > 0 && ( ( *vals )[ z * nPointsInSlice + ( y - 1 ) * nX + x ] != isoValue ) )
345  {
346  addSurface( x, y, z, 3 );
347  }
348 
349  if( y < m_nCellsY - 1 && ( ( *vals )[ z * nPointsInSlice + ( y + 1 ) * nX + x ] != isoValue ) )
350  {
351  addSurface( x, y, z, 4 );
352  }
353 
354  if( z > 0 && ( ( *vals )[ ( z - 1 ) * nPointsInSlice + y * nX + x ] != isoValue ) )
355  {
356  addSurface( x, y, z, 5 );
357  }
358 
359  if( z < m_nCellsZ - 1 && ( ( *vals )[ ( z + 1 ) * nPointsInSlice + y * nX + x ] != isoValue ) )
360  {
361  addSurface( x, y, z, 6 );
362  }
363 
364  if( x == 0 )
365  {
366  addSurface( x, y, z, 1 );
367  }
368  if( x == m_nCellsX - 1 )
369  {
370  addSurface( x, y, z, 2 );
371  }
372 
373  if( y == 0 )
374  {
375  addSurface( x, y, z, 3 );
376  }
377 
378  if( y == m_nCellsY - 1 )
379  {
380  addSurface( x, y, z, 4 );
381  }
382 
383  if( z == 0 )
384  {
385  addSurface( x, y, z, 5 );
386  }
387 
388  if( z == m_nCellsZ - 1 )
389  {
390  addSurface( x, y, z, 6 );
391  }
392  }
393  }
394  }
395  unsigned int nextID = 0;
396  boost::shared_ptr< WTriangleMesh > triMesh( new WTriangleMesh( m_idToVertices.size(), m_trivecTriangles.size() ) );
397 
398  // Rename vertices.
399  ID2WMLPointXYZId::iterator mapIterator = m_idToVertices.begin();
400  while( mapIterator != m_idToVertices.end() )
401  {
402  WPosition texCoord = WPosition( mapIterator->second.x / nbCoordsX,
403  mapIterator->second.y / nbCoordsY,
404  mapIterator->second.z / nbCoordsZ );
405 
406  // transform from grid coordinate system to world coordinates
407  WPosition pos = WPosition( mapIterator->second.x, mapIterator->second.y, mapIterator->second.z );
408 
409  std::vector< double > resultPos4D( 4 );
410  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;
411  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;
412  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;
413  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;
414 
415  ( *mapIterator ).second.newID = nextID;
416  triMesh->addVertex( resultPos4D[0] / resultPos4D[3],
417  resultPos4D[1] / resultPos4D[3],
418  resultPos4D[2] / resultPos4D[3] );
419  triMesh->addTextureCoordinate( texCoord );
420  nextID++;
421  mapIterator++;
422  }
423 
424  // Now rename triangles.
425  WMLTriangleVECTOR::iterator vecIterator = m_trivecTriangles.begin();
426  while( vecIterator != m_trivecTriangles.end() )
427  {
428  for( unsigned int i = 0; i < 3; i++ )
429  {
430  unsigned int newID = m_idToVertices[( *vecIterator ).pointID[i]].newID;
431  ( *vecIterator ).pointID[i] = newID;
432  }
433  triMesh->addTriangle( ( *vecIterator ).pointID[0], ( *vecIterator ).pointID[1], ( *vecIterator ).pointID[2] );
434  vecIterator++;
435  }
436 
437  if( progress )
438  {
439  progress->finish();
440  }
441  return triMesh;
442 }
unsigned int newID
ID of the point.
unsigned int pointID[3]
The IDs of the vertices of the triangle.
ID2WMLPointXYZId m_idToVertices
List of WPointXYZIds which form the isosurface.
A point consisting of its coordinates and ID.
Class managing progress inside of modules.
Definition: WProgress.h:43
Encapsulated ids representing a triangle.
boost::shared_ptr< WTriangleMesh > genSurfaceOneValue(size_t nbCoordsX, size_t nbCoordsY, size_t nbCoordsZ, const WMatrix< double > &mat, const std::vector< size_t > *vals, size_t isoValue, boost::shared_ptr< WProgressCombiner > progress=boost::shared_ptr< WProgressCombiner >())
Generate the triangles for the surface on the given dataSet (inGrid, vals).
void addSurface(size_t x, size_t y, size_t z, size_t surface)
adds 2 triangles for a given face of the voxel
double y
y coordinates of the point.
unsigned int m_nCellsY
No. of cells in y direction.
This only is a 3d double vector.
unsigned int m_nCellsZ
No. of cells in z direction.
double z
z coordinates of the point.
WMLTriangleVECTOR m_trivecTriangles
List of WMCTriangleS which form the triangulation of the isosurface.
Triangle mesh data structure allowing for convenient access of the elements.
Definition: WTriangleMesh.h:44
WMarchingLegoAlgorithm()
standard constructor
WMatrix< double > m_matrix
The 4x4 transformation matrix for the triangle vertices.
unsigned int m_nCellsX
No. of cells in x direction.
size_t getVertexID(size_t nX, size_t nY, size_t nZ)
returns a vertex id for a given grid point
double x
x coordinates of the point.