PolyVox  0.2.1
Open source voxel management library
MarchingCubesSurfaceExtractor.h
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 #ifndef __PolyVox_SurfaceExtractor_H__
25 #define __PolyVox_SurfaceExtractor_H__
26 
28 #include "Impl/TypeDef.h"
29 
30 #include "PolyVoxCore/Array.h"
33 
34 namespace PolyVox
35 {
36  template< typename VolumeType, typename Controller = DefaultMarchingCubesController<typename VolumeType::VoxelType> >
38  {
39  public:
40  MarchingCubesSurfaceExtractor(VolumeType* volData, Region region, SurfaceMesh<PositionMaterialNormal>* result, Controller controller = Controller());
41 
42  void execute();
43 
44  private:
45  //Compute the cell bitmask for a particular slice in z.
46  template<bool isPrevZAvail>
47  uint32_t computeBitmaskForSlice(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask);
48 
49  //Compute the cell bitmask for a given cell.
50  template<bool isPrevXAvail, bool isPrevYAvail, bool isPrevZAvail>
51  void computeBitmaskForCell(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask, uint32_t uXRegSpace, uint32_t uYRegSpace);
52 
53  //Use the cell bitmasks to generate all the vertices needed for that slice
54  void generateVerticesForSlice(const Array2DUint8& pCurrentBitmask,
55  Array2DInt32& m_pCurrentVertexIndicesX,
56  Array2DInt32& m_pCurrentVertexIndicesY,
57  Array2DInt32& m_pCurrentVertexIndicesZ);
58 
60  // NOTE: These two functions are in the .h file rather than the .inl due to an apparent bug in VC2010.
61  //See http://stackoverflow.com/questions/1484885/strange-vc-compile-error-c2244 for details.
63  Vector3DFloat computeCentralDifferenceGradient(const typename VolumeType::Sampler& volIter)
64  {
65  //FIXME - Should actually use DensityType here, both in principle and because the maths may be
66  //faster (and to reduce casts). So it would be good to add a way to get DensityType from a voxel.
67  //But watch out for when the DensityType is unsigned and the difference could be negative.
68  float voxel1nx = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py0pz()));
69  float voxel1px = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py0pz()));
70 
71  float voxel1ny = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny0pz()));
72  float voxel1py = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py0pz()));
73 
74  float voxel1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1nz()));
75  float voxel1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1pz()));
76 
77  return Vector3DFloat
78  (
79  voxel1nx - voxel1px,
80  voxel1ny - voxel1py,
81  voxel1nz - voxel1pz
82  );
83  }
84 
85  Vector3DFloat computeSobelGradient(const typename VolumeType::Sampler& volIter)
86  {
87  static const int weights[3][3][3] = { { {2,3,2}, {3,6,3}, {2,3,2} }, {
88  {3,6,3}, {6,0,6}, {3,6,3} }, { {2,3,2}, {3,6,3}, {2,3,2} } };
89 
90  //FIXME - Should actually use DensityType here, both in principle and because the maths may be
91  //faster (and to reduce casts). So it would be good to add a way to get DensityType from a voxel.
92  //But watch out for when the DensityType is unsigned and the difference could be negative.
93  const float pVoxel1nx1ny1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1ny1nz()));
94  const float pVoxel1nx1ny0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1ny0pz()));
95  const float pVoxel1nx1ny1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1ny1pz()));
96  const float pVoxel1nx0py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py1nz()));
97  const float pVoxel1nx0py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py0pz()));
98  const float pVoxel1nx0py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py1pz()));
99  const float pVoxel1nx1py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1py1nz()));
100  const float pVoxel1nx1py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1py0pz()));
101  const float pVoxel1nx1py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1py1pz()));
102 
103  const float pVoxel0px1ny1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny1nz()));
104  const float pVoxel0px1ny0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny0pz()));
105  const float pVoxel0px1ny1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny1pz()));
106  const float pVoxel0px0py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1nz()));
107  //const float pVoxel0px0py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py0pz()));
108  const float pVoxel0px0py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1pz()));
109  const float pVoxel0px1py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py1nz()));
110  const float pVoxel0px1py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py0pz()));
111  const float pVoxel0px1py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py1pz()));
112 
113  const float pVoxel1px1ny1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1ny1nz()));
114  const float pVoxel1px1ny0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1ny0pz()));
115  const float pVoxel1px1ny1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1ny1pz()));
116  const float pVoxel1px0py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py1nz()));
117  const float pVoxel1px0py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py0pz()));
118  const float pVoxel1px0py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py1pz()));
119  const float pVoxel1px1py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1py1nz()));
120  const float pVoxel1px1py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1py0pz()));
121  const float pVoxel1px1py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1py1pz()));
122 
123  const float xGrad(- weights[0][0][0] * pVoxel1nx1ny1nz -
124  weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] *
125  pVoxel1nx1ny1pz - weights[0][1][0] * pVoxel1nx0py1nz -
126  weights[1][1][0] * pVoxel1nx0py0pz - weights[2][1][0] *
127  pVoxel1nx0py1pz - weights[0][2][0] * pVoxel1nx1py1nz -
128  weights[1][2][0] * pVoxel1nx1py0pz - weights[2][2][0] *
129  pVoxel1nx1py1pz + weights[0][0][2] * pVoxel1px1ny1nz +
130  weights[1][0][2] * pVoxel1px1ny0pz + weights[2][0][2] *
131  pVoxel1px1ny1pz + weights[0][1][2] * pVoxel1px0py1nz +
132  weights[1][1][2] * pVoxel1px0py0pz + weights[2][1][2] *
133  pVoxel1px0py1pz + weights[0][2][2] * pVoxel1px1py1nz +
134  weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] *
135  pVoxel1px1py1pz);
136 
137  const float yGrad(- weights[0][0][0] * pVoxel1nx1ny1nz -
138  weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] *
139  pVoxel1nx1ny1pz + weights[0][2][0] * pVoxel1nx1py1nz +
140  weights[1][2][0] * pVoxel1nx1py0pz + weights[2][2][0] *
141  pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz -
142  weights[1][0][1] * pVoxel0px1ny0pz - weights[2][0][1] *
143  pVoxel0px1ny1pz + weights[0][2][1] * pVoxel0px1py1nz +
144  weights[1][2][1] * pVoxel0px1py0pz + weights[2][2][1] *
145  pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz -
146  weights[1][0][2] * pVoxel1px1ny0pz - weights[2][0][2] *
147  pVoxel1px1ny1pz + weights[0][2][2] * pVoxel1px1py1nz +
148  weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] *
149  pVoxel1px1py1pz);
150 
151  const float zGrad(- weights[0][0][0] * pVoxel1nx1ny1nz +
152  weights[2][0][0] * pVoxel1nx1ny1pz - weights[0][1][0] *
153  pVoxel1nx0py1nz + weights[2][1][0] * pVoxel1nx0py1pz -
154  weights[0][2][0] * pVoxel1nx1py1nz + weights[2][2][0] *
155  pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz +
156  weights[2][0][1] * pVoxel0px1ny1pz - weights[0][1][1] *
157  pVoxel0px0py1nz + weights[2][1][1] * pVoxel0px0py1pz -
158  weights[0][2][1] * pVoxel0px1py1nz + weights[2][2][1] *
159  pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz +
160  weights[2][0][2] * pVoxel1px1ny1pz - weights[0][1][2] *
161  pVoxel1px0py1nz + weights[2][1][2] * pVoxel1px0py1pz -
162  weights[0][2][2] * pVoxel1px1py1nz + weights[2][2][2] *
163  pVoxel1px1py1pz);
164 
165  //Note: The above actually give gradients going from low density to high density.
166  //For our normals we want the the other way around, so we switch the components as we return them.
167  return Vector3DFloat(-xGrad,-yGrad,-zGrad);
168  }
170  // End of compiler bug workaroumd.
172 
173  //Use the cell bitmasks to generate all the indices needed for that slice
174  void generateIndicesForSlice(const Array2DUint8& pPreviousBitmask,
175  const Array2DInt32& m_pPreviousVertexIndicesX,
176  const Array2DInt32& m_pPreviousVertexIndicesY,
177  const Array2DInt32& m_pPreviousVertexIndicesZ,
178  const Array2DInt32& m_pCurrentVertexIndicesX,
179  const Array2DInt32& m_pCurrentVertexIndicesY);
180 
181  //The volume data and a sampler to access it.
182  VolumeType* m_volData;
183  typename VolumeType::Sampler m_sampVolume;
184 
185  //Used to return the number of cells in a slice which contain triangles.
186  uint32_t m_uNoOfOccupiedCells;
187 
188  //The surface patch we are currently filling.
190 
191  //Information about the region we are currently processing
192  Region m_regSizeInVoxels;
193  Region m_regSizeInCells;
194  /*Region m_regSizeInVoxelsCropped;
195  Region m_regSizeInVoxelsUncropped;
196  Region m_regVolumeCropped;*/
197  Region m_regSlicePrevious;
198  Region m_regSliceCurrent;
199 
200  //Our threshold value
201  typename Controller::DensityType m_tThreshold;
202 
203  //Used to convert arbitrary voxel types in densities and materials.
204  Controller m_controller;
205  };
206 }
207 
209 
210 #endif