PolyVox  0.2.1
Open source voxel management library
MarchingCubesSurfaceExtractor.inl
Go to the documentation of this file.
1 /*******************************************************************************
2 Copyright (c) 2005-2009 David Williams
3 
4 This software is provided 'as-is', without any express or implied
5 warranty. In no event will the authors be held liable for any damages
6 arising from the use of this software.
7 
8 Permission is granted to anyone to use this software for any purpose,
9 including commercial applications, and to alter it and redistribute it
10 freely, subject to the following restrictions:
11 
12  1. The origin of this software must not be misrepresented; you must not
13  claim that you wrote the original software. If you use this software
14  in a product, an acknowledgment in the product documentation would be
15  appreciated but is not required.
16 
17  2. Altered source versions must be plainly marked as such, and must not be
18  misrepresented as being the original software.
19 
20  3. This notice may not be removed or altered from any source
21  distribution.
22 *******************************************************************************/
23 
24 namespace PolyVox
25 {
26  template<typename VolumeType, typename Controller>
27  MarchingCubesSurfaceExtractor<VolumeType, Controller>::MarchingCubesSurfaceExtractor(VolumeType* volData, Region region, SurfaceMesh<PositionMaterialNormal>* result, Controller controller)
28  :m_volData(volData)
29  ,m_sampVolume(volData)
30  ,m_meshCurrent(result)
31  ,m_regSizeInVoxels(region)
32  {
33  //m_regSizeInVoxels.cropTo(m_volData->getEnclosingRegion());
34  m_regSizeInCells = m_regSizeInVoxels;
35  m_regSizeInCells.setUpperCorner(m_regSizeInCells.getUpperCorner() - Vector3DInt32(1,1,1));
36 
37  m_controller = controller;
38  m_tThreshold = m_controller.getThreshold();
39  }
40 
41  template<typename VolumeType, typename Controller>
42  void MarchingCubesSurfaceExtractor<VolumeType, Controller>::execute()
43  {
44  m_meshCurrent->clear();
45 
46  uint32_t uArrayWidth = m_regSizeInVoxels.getUpperCorner().getX() - m_regSizeInVoxels.getLowerCorner().getX() + 1;
47  uint32_t uArrayHeight = m_regSizeInVoxels.getUpperCorner().getY() - m_regSizeInVoxels.getLowerCorner().getY() + 1;
48  uint32_t arraySizes[2]= {uArrayWidth, uArrayHeight}; // Array dimensions
49 
50  //For edge indices
51  Array2DInt32 m_pPreviousVertexIndicesX(arraySizes);
52  Array2DInt32 m_pPreviousVertexIndicesY(arraySizes);
53  Array2DInt32 m_pPreviousVertexIndicesZ(arraySizes);
54  Array2DInt32 m_pCurrentVertexIndicesX(arraySizes);
55  Array2DInt32 m_pCurrentVertexIndicesY(arraySizes);
56  Array2DInt32 m_pCurrentVertexIndicesZ(arraySizes);
57 
58  Array2DUint8 pPreviousBitmask(arraySizes);
59  Array2DUint8 pCurrentBitmask(arraySizes);
60 
61  //Create a region corresponding to the first slice
62  m_regSlicePrevious = m_regSizeInVoxels;
63  Vector3DInt32 v3dUpperCorner = m_regSlicePrevious.getUpperCorner();
64  v3dUpperCorner.setZ(m_regSlicePrevious.getLowerCorner().getZ()); //Set the upper z to the lower z to make it one slice thick.
65  m_regSlicePrevious.setUpperCorner(v3dUpperCorner);
66  m_regSliceCurrent = m_regSlicePrevious;
67 
68  uint32_t uNoOfNonEmptyCellsForSlice0 = 0;
69  uint32_t uNoOfNonEmptyCellsForSlice1 = 0;
70 
71  //Process the first slice (previous slice not available)
72  computeBitmaskForSlice<false>(pPreviousBitmask, pCurrentBitmask);
73  uNoOfNonEmptyCellsForSlice1 = m_uNoOfOccupiedCells;
74 
75  if(uNoOfNonEmptyCellsForSlice1 != 0)
76  {
77  memset(m_pCurrentVertexIndicesX.getRawData(), 0xff, m_pCurrentVertexIndicesX.getNoOfElements() * 4);
78  memset(m_pCurrentVertexIndicesY.getRawData(), 0xff, m_pCurrentVertexIndicesY.getNoOfElements() * 4);
79  memset(m_pCurrentVertexIndicesZ.getRawData(), 0xff, m_pCurrentVertexIndicesZ.getNoOfElements() * 4);
80  generateVerticesForSlice(pCurrentBitmask, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY, m_pCurrentVertexIndicesZ);
81  }
82 
83  std::swap(uNoOfNonEmptyCellsForSlice0, uNoOfNonEmptyCellsForSlice1);
84  pPreviousBitmask.swap(pCurrentBitmask);
85  m_pPreviousVertexIndicesX.swap(m_pCurrentVertexIndicesX);
86  m_pPreviousVertexIndicesY.swap(m_pCurrentVertexIndicesY);
87  m_pPreviousVertexIndicesZ.swap(m_pCurrentVertexIndicesZ);
88 
89  m_regSlicePrevious = m_regSliceCurrent;
90  m_regSliceCurrent.shift(Vector3DInt32(0,0,1));
91 
92  //Process the other slices (previous slice is available)
93  for(int32_t uSlice = 1; uSlice <= m_regSizeInVoxels.getUpperCorner().getZ() - m_regSizeInVoxels.getLowerCorner().getZ(); uSlice++)
94  {
95  computeBitmaskForSlice<true>(pPreviousBitmask, pCurrentBitmask);
96  uNoOfNonEmptyCellsForSlice1 = m_uNoOfOccupiedCells;
97 
98  if(uNoOfNonEmptyCellsForSlice1 != 0)
99  {
100  memset(m_pCurrentVertexIndicesX.getRawData(), 0xff, m_pCurrentVertexIndicesX.getNoOfElements() * 4);
101  memset(m_pCurrentVertexIndicesY.getRawData(), 0xff, m_pCurrentVertexIndicesY.getNoOfElements() * 4);
102  memset(m_pCurrentVertexIndicesZ.getRawData(), 0xff, m_pCurrentVertexIndicesZ.getNoOfElements() * 4);
103  generateVerticesForSlice(pCurrentBitmask, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY, m_pCurrentVertexIndicesZ);
104  }
105 
106  if((uNoOfNonEmptyCellsForSlice0 != 0) || (uNoOfNonEmptyCellsForSlice1 != 0))
107  {
108  generateIndicesForSlice(pPreviousBitmask, m_pPreviousVertexIndicesX, m_pPreviousVertexIndicesY, m_pPreviousVertexIndicesZ, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY);
109  }
110 
111  std::swap(uNoOfNonEmptyCellsForSlice0, uNoOfNonEmptyCellsForSlice1);
112  pPreviousBitmask.swap(pCurrentBitmask);
113  m_pPreviousVertexIndicesX.swap(m_pCurrentVertexIndicesX);
114  m_pPreviousVertexIndicesY.swap(m_pCurrentVertexIndicesY);
115  m_pPreviousVertexIndicesZ.swap(m_pCurrentVertexIndicesZ);
116 
117  m_regSlicePrevious = m_regSliceCurrent;
118  m_regSliceCurrent.shift(Vector3DInt32(0,0,1));
119  }
120 
121  m_meshCurrent->m_Region = m_regSizeInVoxels;
122 
123  m_meshCurrent->m_vecLodRecords.clear();
124  LodRecord lodRecord;
125  lodRecord.beginIndex = 0;
126  lodRecord.endIndex = m_meshCurrent->getNoOfIndices();
127  m_meshCurrent->m_vecLodRecords.push_back(lodRecord);
128  }
129 
130  template<typename VolumeType, typename Controller>
131  template<bool isPrevZAvail>
132  uint32_t MarchingCubesSurfaceExtractor<VolumeType, Controller>::computeBitmaskForSlice(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask)
133  {
134  m_uNoOfOccupiedCells = 0;
135 
136  const int32_t iMaxXVolSpace = m_regSliceCurrent.getUpperCorner().getX();
137  const int32_t iMaxYVolSpace = m_regSliceCurrent.getUpperCorner().getY();
138 
139  int32_t iZVolSpace = m_regSliceCurrent.getLowerCorner().getZ();
140 
141  //Process the lower left corner
142  int32_t iYVolSpace = m_regSliceCurrent.getLowerCorner().getY();
143  int32_t iXVolSpace = m_regSliceCurrent.getLowerCorner().getX();
144 
145  uint32_t uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX();
146  uint32_t uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY();
147 
148  m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace);
149  computeBitmaskForCell<false, false, isPrevZAvail>(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace);
150 
151  //Process the edge where x is minimal.
152  iXVolSpace = m_regSliceCurrent.getLowerCorner().getX();
153  m_sampVolume.setPosition(iXVolSpace, m_regSliceCurrent.getLowerCorner().getY(), iZVolSpace);
154  for(iYVolSpace = m_regSliceCurrent.getLowerCorner().getY() + 1; iYVolSpace <= iMaxYVolSpace; iYVolSpace++)
155  {
156  uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX();
157  uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY();
158 
159  m_sampVolume.movePositiveY();
160 
161  computeBitmaskForCell<false, true, isPrevZAvail>(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace);
162  }
163 
164  //Process the edge where y is minimal.
165  iYVolSpace = m_regSliceCurrent.getLowerCorner().getY();
166  m_sampVolume.setPosition(m_regSliceCurrent.getLowerCorner().getX(), iYVolSpace, iZVolSpace);
167  for(iXVolSpace = m_regSliceCurrent.getLowerCorner().getX() + 1; iXVolSpace <= iMaxXVolSpace; iXVolSpace++)
168  {
169  uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX();
170  uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY();
171 
172  m_sampVolume.movePositiveX();
173 
174  computeBitmaskForCell<true, false, isPrevZAvail>(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace);
175  }
176 
177  //Process all remaining elemnents of the slice. In this case, previous x and y values are always available
178  for(iYVolSpace = m_regSliceCurrent.getLowerCorner().getY() + 1; iYVolSpace <= iMaxYVolSpace; iYVolSpace++)
179  {
180  m_sampVolume.setPosition(m_regSliceCurrent.getLowerCorner().getX(), iYVolSpace, iZVolSpace);
181  for(iXVolSpace = m_regSliceCurrent.getLowerCorner().getX() + 1; iXVolSpace <= iMaxXVolSpace; iXVolSpace++)
182  {
183  uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX();
184  uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY();
185 
186  m_sampVolume.movePositiveX();
187 
188  computeBitmaskForCell<true, true, isPrevZAvail>(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace);
189  }
190  }
191 
192  return m_uNoOfOccupiedCells;
193  }
194 
195  template<typename VolumeType, typename Controller>
196  template<bool isPrevXAvail, bool isPrevYAvail, bool isPrevZAvail>
197  void MarchingCubesSurfaceExtractor<VolumeType, Controller>::computeBitmaskForCell(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask, uint32_t uXRegSpace, uint32_t uYRegSpace)
198  {
199  uint8_t iCubeIndex = 0;
200 
201  typename VolumeType::VoxelType v000;
202  typename VolumeType::VoxelType v100;
203  typename VolumeType::VoxelType v010;
204  typename VolumeType::VoxelType v110;
205  typename VolumeType::VoxelType v001;
206  typename VolumeType::VoxelType v101;
207  typename VolumeType::VoxelType v011;
208  typename VolumeType::VoxelType v111;
209 
210  if(isPrevZAvail)
211  {
212  if(isPrevYAvail)
213  {
214  if(isPrevXAvail)
215  {
216  v111 = m_sampVolume.peekVoxel1px1py1pz();
217 
218  //z
219  uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace];
220  iPreviousCubeIndexZ >>= 4;
221 
222  //y
223  uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1];
224  iPreviousCubeIndexY &= 192; //192 = 128 + 64
225  iPreviousCubeIndexY >>= 2;
226 
227  //x
228  uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace];
229  iPreviousCubeIndexX &= 128;
230  iPreviousCubeIndexX >>= 1;
231 
232  iCubeIndex = iPreviousCubeIndexX | iPreviousCubeIndexY | iPreviousCubeIndexZ;
233 
234  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
235  }
236  else //previous X not available
237  {
238  v011 = m_sampVolume.peekVoxel0px1py1pz();
239  v111 = m_sampVolume.peekVoxel1px1py1pz();
240 
241  //z
242  uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace];
243  iPreviousCubeIndexZ >>= 4;
244 
245  //y
246  uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1];
247  iPreviousCubeIndexY &= 192; //192 = 128 + 64
248  iPreviousCubeIndexY >>= 2;
249 
250  iCubeIndex = iPreviousCubeIndexY | iPreviousCubeIndexZ;
251 
252  if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64;
253  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
254  }
255  }
256  else //previous Y not available
257  {
258  if(isPrevXAvail)
259  {
260  v101 = m_sampVolume.peekVoxel1px0py1pz();
261  v111 = m_sampVolume.peekVoxel1px1py1pz();
262 
263  //z
264  uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace];
265  iPreviousCubeIndexZ >>= 4;
266 
267  //x
268  uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace];
269  iPreviousCubeIndexX &= 160; //160 = 128+32
270  iPreviousCubeIndexX >>= 1;
271 
272  iCubeIndex = iPreviousCubeIndexX | iPreviousCubeIndexZ;
273 
274  if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32;
275  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
276  }
277  else //previous X not available
278  {
279  v001 = m_sampVolume.peekVoxel0px0py1pz();
280  v101 = m_sampVolume.peekVoxel1px0py1pz();
281  v011 = m_sampVolume.peekVoxel0px1py1pz();
282  v111 = m_sampVolume.peekVoxel1px1py1pz();
283 
284  //z
285  uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace];
286  iCubeIndex = iPreviousCubeIndexZ >> 4;
287 
288  if (m_controller.convertToDensity(v001) < m_tThreshold) iCubeIndex |= 16;
289  if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32;
290  if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64;
291  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
292  }
293  }
294  }
295  else //previous Z not available
296  {
297  if(isPrevYAvail)
298  {
299  if(isPrevXAvail)
300  {
301  v110 = m_sampVolume.peekVoxel1px1py0pz();
302  v111 = m_sampVolume.peekVoxel1px1py1pz();
303 
304  //y
305  uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1];
306  iPreviousCubeIndexY &= 204; //204 = 128+64+8+4
307  iPreviousCubeIndexY >>= 2;
308 
309  //x
310  uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace];
311  iPreviousCubeIndexX &= 170; //170 = 128+32+8+2
312  iPreviousCubeIndexX >>= 1;
313 
314  iCubeIndex = iPreviousCubeIndexX | iPreviousCubeIndexY;
315 
316  if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8;
317  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
318  }
319  else //previous X not available
320  {
321  v010 = m_sampVolume.peekVoxel0px1py0pz();
322  v110 = m_sampVolume.peekVoxel1px1py0pz();
323 
324  v011 = m_sampVolume.peekVoxel0px1py1pz();
325  v111 = m_sampVolume.peekVoxel1px1py1pz();
326 
327  //y
328  uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1];
329  iPreviousCubeIndexY &= 204; //204 = 128+64+8+4
330  iPreviousCubeIndexY >>= 2;
331 
332  iCubeIndex = iPreviousCubeIndexY;
333 
334  if (m_controller.convertToDensity(v010) < m_tThreshold) iCubeIndex |= 4;
335  if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8;
336  if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64;
337  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
338  }
339  }
340  else //previous Y not available
341  {
342  if(isPrevXAvail)
343  {
344  v100 = m_sampVolume.peekVoxel1px0py0pz();
345  v110 = m_sampVolume.peekVoxel1px1py0pz();
346 
347  v101 = m_sampVolume.peekVoxel1px0py1pz();
348  v111 = m_sampVolume.peekVoxel1px1py1pz();
349 
350  //x
351  uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace];
352  iPreviousCubeIndexX &= 170; //170 = 128+32+8+2
353  iPreviousCubeIndexX >>= 1;
354 
355  iCubeIndex = iPreviousCubeIndexX;
356 
357  if (m_controller.convertToDensity(v100) < m_tThreshold) iCubeIndex |= 2;
358  if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8;
359  if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32;
360  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
361  }
362  else //previous X not available
363  {
364  v000 = m_sampVolume.getVoxel();
365  v100 = m_sampVolume.peekVoxel1px0py0pz();
366  v010 = m_sampVolume.peekVoxel0px1py0pz();
367  v110 = m_sampVolume.peekVoxel1px1py0pz();
368 
369  v001 = m_sampVolume.peekVoxel0px0py1pz();
370  v101 = m_sampVolume.peekVoxel1px0py1pz();
371  v011 = m_sampVolume.peekVoxel0px1py1pz();
372  v111 = m_sampVolume.peekVoxel1px1py1pz();
373 
374  if (m_controller.convertToDensity(v000) < m_tThreshold) iCubeIndex |= 1;
375  if (m_controller.convertToDensity(v100) < m_tThreshold) iCubeIndex |= 2;
376  if (m_controller.convertToDensity(v010) < m_tThreshold) iCubeIndex |= 4;
377  if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8;
378  if (m_controller.convertToDensity(v001) < m_tThreshold) iCubeIndex |= 16;
379  if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32;
380  if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64;
381  if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128;
382  }
383  }
384  }
385 
386  //Save the bitmask
387  pCurrentBitmask[uXRegSpace][uYRegSpace] = iCubeIndex;
388 
389  if(edgeTable[iCubeIndex] != 0)
390  {
391  ++m_uNoOfOccupiedCells;
392  }
393  }
394 
395  template<typename VolumeType, typename Controller>
396  void MarchingCubesSurfaceExtractor<VolumeType, Controller>::generateVerticesForSlice(const Array2DUint8& pCurrentBitmask,
397  Array2DInt32& m_pCurrentVertexIndicesX,
398  Array2DInt32& m_pCurrentVertexIndicesY,
399  Array2DInt32& m_pCurrentVertexIndicesZ)
400  {
401  int32_t iZVolSpace = m_regSliceCurrent.getLowerCorner().getZ();
402 
403  //Iterate over each cell in the region
404  for(int32_t iYVolSpace = m_regSliceCurrent.getLowerCorner().getY(); iYVolSpace <= m_regSliceCurrent.getUpperCorner().getY(); iYVolSpace++)
405  {
406  const uint32_t uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY();
407 
408  for(int32_t iXVolSpace = m_regSliceCurrent.getLowerCorner().getX(); iXVolSpace <= m_regSliceCurrent.getUpperCorner().getX(); iXVolSpace++)
409  {
410  //Current position
411  const uint32_t uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX();
412 
413  //Determine the index into the edge table which tells us which vertices are inside of the surface
414  uint8_t iCubeIndex = pCurrentBitmask[uXRegSpace][uYRegSpace];
415 
416  /* Cube is entirely in/out of the surface */
417  if (edgeTable[iCubeIndex] == 0)
418  {
419  continue;
420  }
421 
422  //Check whether the generated vertex will lie on the edge of the region
423 
424 
425  m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace);
426  const typename VolumeType::VoxelType v000 = m_sampVolume.getVoxel();
427  const Vector3DFloat n000 = computeCentralDifferenceGradient(m_sampVolume);
428 
429  /* Find the vertices where the surface intersects the cube */
430  if (edgeTable[iCubeIndex] & 1)
431  {
432  m_sampVolume.movePositiveX();
433  const typename VolumeType::VoxelType v100 = m_sampVolume.getVoxel();
434  const Vector3DFloat n100 = computeCentralDifferenceGradient(m_sampVolume);
435 
436  float fInterp = static_cast<float>(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast<float>(m_controller.convertToDensity(v100) - m_controller.convertToDensity(v000));
437 
438  const Vector3DFloat v3dPosition(static_cast<float>(iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()) + fInterp, static_cast<float>(iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()), static_cast<float>(iZVolSpace - m_regSizeInCells.getLowerCorner().getZ()));
439 
440  Vector3DFloat v3dNormal = (n100*fInterp) + (n000*(1-fInterp));
441  v3dNormal.normalise();
442 
443  //Choose one of the two materials to use for the vertex (we don't interpolate as interpolation of
444  //material IDs does not make sense). We take the largest, so that if we are working on a material-only
445  //volume we get the one which is non-zero. Both materials can be non-zero if our volume has a density component.
446  typename Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000);
447  typename Controller::MaterialType uMaterial100 = m_controller.convertToMaterial(v100);
448  typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial100);
449 
450  PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, static_cast<float>(uMaterial));
451  uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex);
452  m_pCurrentVertexIndicesX[iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()][iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()] = uLastVertexIndex;
453 
454  m_sampVolume.moveNegativeX();
455  }
456  if (edgeTable[iCubeIndex] & 8)
457  {
458  m_sampVolume.movePositiveY();
459  const typename VolumeType::VoxelType v010 = m_sampVolume.getVoxel();
460  const Vector3DFloat n010 = computeCentralDifferenceGradient(m_sampVolume);
461 
462  float fInterp = static_cast<float>(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast<float>(m_controller.convertToDensity(v010) - m_controller.convertToDensity(v000));
463 
464  const Vector3DFloat v3dPosition(static_cast<float>(iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()), static_cast<float>(iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()) + fInterp, static_cast<float>(iZVolSpace - m_regSizeInVoxels.getLowerCorner().getZ()));
465 
466  Vector3DFloat v3dNormal = (n010*fInterp) + (n000*(1-fInterp));
467  v3dNormal.normalise();
468 
469  //Choose one of the two materials to use for the vertex (we don't interpolate as interpolation of
470  //material IDs does not make sense). We take the largest, so that if we are working on a material-only
471  //volume we get the one which is non-zero. Both materials can be non-zero if our volume has a density component.
472  typename Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000);
473  typename Controller::MaterialType uMaterial010 = m_controller.convertToMaterial(v010);
474  typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial010);
475 
476  PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, static_cast<float>(uMaterial));
477  uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex);
478  m_pCurrentVertexIndicesY[iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()][iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()] = uLastVertexIndex;
479 
480  m_sampVolume.moveNegativeY();
481  }
482  if (edgeTable[iCubeIndex] & 256)
483  {
484  m_sampVolume.movePositiveZ();
485  const typename VolumeType::VoxelType v001 = m_sampVolume.getVoxel();
486  const Vector3DFloat n001 = computeCentralDifferenceGradient(m_sampVolume);
487 
488  float fInterp = static_cast<float>(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast<float>(m_controller.convertToDensity(v001) - m_controller.convertToDensity(v000));
489 
490  const Vector3DFloat v3dPosition(static_cast<float>(iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()), static_cast<float>(iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()), static_cast<float>(iZVolSpace - m_regSizeInVoxels.getLowerCorner().getZ()) + fInterp);
491 
492  Vector3DFloat v3dNormal = (n001*fInterp) + (n000*(1-fInterp));
493  v3dNormal.normalise();
494 
495  //Choose one of the two materials to use for the vertex (we don't interpolate as interpolation of
496  //material IDs does not make sense). We take the largest, so that if we are working on a material-only
497  //volume we get the one which is non-zero. Both materials can be non-zero if our volume has a density component.
498  typename Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000);
499  typename Controller::MaterialType uMaterial001 = m_controller.convertToMaterial(v001);
500  typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial001);
501 
502  PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, static_cast<float>(uMaterial));
503  uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex);
504  m_pCurrentVertexIndicesZ[iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()][iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()] = uLastVertexIndex;
505 
506  m_sampVolume.moveNegativeZ();
507  }
508  }//For each cell
509  }
510  }
511 
512  template<typename VolumeType, typename Controller>
513  void MarchingCubesSurfaceExtractor<VolumeType, Controller>::generateIndicesForSlice(const Array2DUint8& pPreviousBitmask,
514  const Array2DInt32& m_pPreviousVertexIndicesX,
515  const Array2DInt32& m_pPreviousVertexIndicesY,
516  const Array2DInt32& m_pPreviousVertexIndicesZ,
517  const Array2DInt32& m_pCurrentVertexIndicesX,
518  const Array2DInt32& m_pCurrentVertexIndicesY)
519  {
520  int32_t indlist[12];
521  for(int i = 0; i < 12; i++)
522  {
523  indlist[i] = -1;
524  }
525 
526  for(int32_t iYVolSpace = m_regSlicePrevious.getLowerCorner().getY(); iYVolSpace <= m_regSizeInCells.getUpperCorner().getY(); iYVolSpace++)
527  {
528  for(int32_t iXVolSpace = m_regSlicePrevious.getLowerCorner().getX(); iXVolSpace <= m_regSizeInCells.getUpperCorner().getX(); iXVolSpace++)
529  {
530  int32_t iZVolSpace = m_regSlicePrevious.getLowerCorner().getZ();
531  m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace);
532 
533  //Current position
534  const uint32_t uXRegSpace = m_sampVolume.getPosition().getX() - m_regSizeInVoxels.getLowerCorner().getX();
535  const uint32_t uYRegSpace = m_sampVolume.getPosition().getY() - m_regSizeInVoxels.getLowerCorner().getY();
536 
537  //Determine the index into the edge table which tells us which vertices are inside of the surface
538  uint8_t iCubeIndex = pPreviousBitmask[uXRegSpace][uYRegSpace];
539 
540  /* Cube is entirely in/out of the surface */
541  if (edgeTable[iCubeIndex] == 0)
542  {
543  continue;
544  }
545 
546  /* Find the vertices where the surface intersects the cube */
547  if (edgeTable[iCubeIndex] & 1)
548  {
549  indlist[0] = m_pPreviousVertexIndicesX[uXRegSpace][uYRegSpace];
550  //assert(indlist[0] != -1);
551  }
552  if (edgeTable[iCubeIndex] & 2)
553  {
554  indlist[1] = m_pPreviousVertexIndicesY[uXRegSpace+1][uYRegSpace];
555  //assert(indlist[1] != -1);
556  }
557  if (edgeTable[iCubeIndex] & 4)
558  {
559  indlist[2] = m_pPreviousVertexIndicesX[uXRegSpace][uYRegSpace+1];
560  //assert(indlist[2] != -1);
561  }
562  if (edgeTable[iCubeIndex] & 8)
563  {
564  indlist[3] = m_pPreviousVertexIndicesY[uXRegSpace][uYRegSpace];
565  //assert(indlist[3] != -1);
566  }
567  if (edgeTable[iCubeIndex] & 16)
568  {
569  indlist[4] = m_pCurrentVertexIndicesX[uXRegSpace][uYRegSpace];
570  //assert(indlist[4] != -1);
571  }
572  if (edgeTable[iCubeIndex] & 32)
573  {
574  indlist[5] = m_pCurrentVertexIndicesY[uXRegSpace+1][uYRegSpace];
575  //assert(indlist[5] != -1);
576  }
577  if (edgeTable[iCubeIndex] & 64)
578  {
579  indlist[6] = m_pCurrentVertexIndicesX[uXRegSpace][uYRegSpace+1];
580  //assert(indlist[6] != -1);
581  }
582  if (edgeTable[iCubeIndex] & 128)
583  {
584  indlist[7] = m_pCurrentVertexIndicesY[uXRegSpace][uYRegSpace];
585  //assert(indlist[7] != -1);
586  }
587  if (edgeTable[iCubeIndex] & 256)
588  {
589  indlist[8] = m_pPreviousVertexIndicesZ[uXRegSpace][uYRegSpace];
590  //assert(indlist[8] != -1);
591  }
592  if (edgeTable[iCubeIndex] & 512)
593  {
594  indlist[9] = m_pPreviousVertexIndicesZ[uXRegSpace+1][uYRegSpace];
595  //assert(indlist[9] != -1);
596  }
597  if (edgeTable[iCubeIndex] & 1024)
598  {
599  indlist[10] = m_pPreviousVertexIndicesZ[uXRegSpace+1][uYRegSpace+1];
600  //assert(indlist[10] != -1);
601  }
602  if (edgeTable[iCubeIndex] & 2048)
603  {
604  indlist[11] = m_pPreviousVertexIndicesZ[uXRegSpace][uYRegSpace+1];
605  //assert(indlist[11] != -1);
606  }
607 
608  for (int i=0;triTable[iCubeIndex][i]!=-1;i+=3)
609  {
610  int32_t ind0 = indlist[triTable[iCubeIndex][i ]];
611  int32_t ind1 = indlist[triTable[iCubeIndex][i+1]];
612  int32_t ind2 = indlist[triTable[iCubeIndex][i+2]];
613 
614  if((ind0 != -1) && (ind1 != -1) && (ind2 != -1))
615  {
616  m_meshCurrent->addTriangle(ind0, ind1, ind2);
617  }
618  }//For each triangle
619  }//For each cell
620  }
621  }
622 }
Generated by   doxygen 1.8.1