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