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