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