diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h new file mode 100644 index 00000000..3daa6dae --- /dev/null +++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.h @@ -0,0 +1,222 @@ +/******************************************************************************* +Copyright (c) 2005-2009 David Williams + +This software is provided 'as-is', without any express or implied +warranty. In no event will the authors be held liable for any damages +arising from the use of this software. + +Permission is granted to anyone to use this software for any purpose, +including commercial applications, and to alter it and redistribute it +freely, subject to the following restrictions: + + 1. The origin of this software must not be misrepresented; you must not + claim that you wrote the original software. If you use this software + in a product, an acknowledgment in the product documentation would be + appreciated but is not required. + + 2. Altered source versions must be plainly marked as such, and must not be + misrepresented as being the original software. + + 3. This notice may not be removed or altered from any source + distribution. +*******************************************************************************/ + +#ifndef __PolyVox_SurfaceExtractor_H__ +#define __PolyVox_SurfaceExtractor_H__ + +#include "PolyVoxImpl/MarchingCubesTables.h" +#include "PolyVoxImpl/TypeDef.h" + +#include "PolyVoxCore/Array.h" +#include "PolyVoxCore/SurfaceMesh.h" +#include "PolyVoxCore/MarchingCubesController.h" + +namespace PolyVox +{ + template< typename VolumeType, typename Controller = MarchingCubesController > + class MarchingCubesSurfaceExtractor + { + public: + MarchingCubesSurfaceExtractor(VolumeType* volData, Region region, SurfaceMesh* result, Controller controller = Controller()); + + void execute(); + + private: + //Compute the cell bitmask for a particular slice in z. + template + uint32_t computeBitmaskForSlice(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask); + + //Compute the cell bitmask for a given cell. + template + void computeBitmaskForCell(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask); + + //Use the cell bitmasks to generate all the vertices needed for that slice + void generateVerticesForSlice(const Array2DUint8& pCurrentBitmask, + Array2DInt32& m_pCurrentVertexIndicesX, + Array2DInt32& m_pCurrentVertexIndicesY, + Array2DInt32& m_pCurrentVertexIndicesZ); + + //////////////////////////////////////////////////////////////////////////////// + // NOTE: These two functions are in the .h file rather than the .inl due to an apparent bug in VC2010. + //See http://stackoverflow.com/questions/1484885/strange-vc-compile-error-c2244 for details. + //////////////////////////////////////////////////////////////////////////////// + Vector3DFloat computeCentralDifferenceGradient(const typename VolumeType::Sampler& volIter) + { + //FIXME - Should actually use DensityType here, both in principle and because the maths may be + //faster (and to reduce casts). So it would be good to add a way to get DensityType from a voxel. + //But watch out for when the DensityType is unsigned and the difference could be negative. + float voxel1nx = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx0py0pz())); + float voxel1px = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px0py0pz())); + + float voxel1ny = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1ny0pz())); + float voxel1py = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1py0pz())); + + float voxel1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px0py1nz())); + float voxel1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px0py1pz())); + + return Vector3DFloat + ( + voxel1nx - voxel1px, + voxel1ny - voxel1py, + voxel1nz - voxel1pz + ); + } + + Vector3DFloat computeSobelGradient(const typename VolumeType::Sampler& volIter) + { + static const int weights[3][3][3] = { { {2,3,2}, {3,6,3}, {2,3,2} }, { + {3,6,3}, {6,0,6}, {3,6,3} }, { {2,3,2}, {3,6,3}, {2,3,2} } }; + + //FIXME - Should actually use DensityType here, both in principle and because the maths may be + //faster (and to reduce casts). So it would be good to add a way to get DensityType from a voxel. + //But watch out for when the DensityType is unsigned and the difference could be negative. + const float pVoxel1nx1ny1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx1ny1nz())); + const float pVoxel1nx1ny0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx1ny0pz())); + const float pVoxel1nx1ny1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx1ny1pz())); + const float pVoxel1nx0py1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx0py1nz())); + const float pVoxel1nx0py0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx0py0pz())); + const float pVoxel1nx0py1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx0py1pz())); + const float pVoxel1nx1py1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx1py1nz())); + const float pVoxel1nx1py0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx1py0pz())); + const float pVoxel1nx1py1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1nx1py1pz())); + + const float pVoxel0px1ny1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1ny1nz())); + const float pVoxel0px1ny0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1ny0pz())); + const float pVoxel0px1ny1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1ny1pz())); + const float pVoxel0px0py1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px0py1nz())); + //const float pVoxel0px0py0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px0py0pz())); + const float pVoxel0px0py1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px0py1pz())); + const float pVoxel0px1py1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1py1nz())); + const float pVoxel0px1py0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1py0pz())); + const float pVoxel0px1py1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel0px1py1pz())); + + const float pVoxel1px1ny1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px1ny1nz())); + const float pVoxel1px1ny0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px1ny0pz())); + const float pVoxel1px1ny1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px1ny1pz())); + const float pVoxel1px0py1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px0py1nz())); + const float pVoxel1px0py0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px0py0pz())); + const float pVoxel1px0py1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px0py1pz())); + const float pVoxel1px1py1nz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px1py1nz())); + const float pVoxel1px1py0pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px1py0pz())); + const float pVoxel1px1py1pz = static_cast(m_controller.convertToDensity(volIter.peekVoxel1px1py1pz())); + + const float xGrad(- weights[0][0][0] * pVoxel1nx1ny1nz - + weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] * + pVoxel1nx1ny1pz - weights[0][1][0] * pVoxel1nx0py1nz - + weights[1][1][0] * pVoxel1nx0py0pz - weights[2][1][0] * + pVoxel1nx0py1pz - weights[0][2][0] * pVoxel1nx1py1nz - + weights[1][2][0] * pVoxel1nx1py0pz - weights[2][2][0] * + pVoxel1nx1py1pz + weights[0][0][2] * pVoxel1px1ny1nz + + weights[1][0][2] * pVoxel1px1ny0pz + weights[2][0][2] * + pVoxel1px1ny1pz + weights[0][1][2] * pVoxel1px0py1nz + + weights[1][1][2] * pVoxel1px0py0pz + weights[2][1][2] * + pVoxel1px0py1pz + weights[0][2][2] * pVoxel1px1py1nz + + weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] * + pVoxel1px1py1pz); + + const float yGrad(- weights[0][0][0] * pVoxel1nx1ny1nz - + weights[1][0][0] * pVoxel1nx1ny0pz - weights[2][0][0] * + pVoxel1nx1ny1pz + weights[0][2][0] * pVoxel1nx1py1nz + + weights[1][2][0] * pVoxel1nx1py0pz + weights[2][2][0] * + pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz - + weights[1][0][1] * pVoxel0px1ny0pz - weights[2][0][1] * + pVoxel0px1ny1pz + weights[0][2][1] * pVoxel0px1py1nz + + weights[1][2][1] * pVoxel0px1py0pz + weights[2][2][1] * + pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz - + weights[1][0][2] * pVoxel1px1ny0pz - weights[2][0][2] * + pVoxel1px1ny1pz + weights[0][2][2] * pVoxel1px1py1nz + + weights[1][2][2] * pVoxel1px1py0pz + weights[2][2][2] * + pVoxel1px1py1pz); + + const float zGrad(- weights[0][0][0] * pVoxel1nx1ny1nz + + weights[2][0][0] * pVoxel1nx1ny1pz - weights[0][1][0] * + pVoxel1nx0py1nz + weights[2][1][0] * pVoxel1nx0py1pz - + weights[0][2][0] * pVoxel1nx1py1nz + weights[2][2][0] * + pVoxel1nx1py1pz - weights[0][0][1] * pVoxel0px1ny1nz + + weights[2][0][1] * pVoxel0px1ny1pz - weights[0][1][1] * + pVoxel0px0py1nz + weights[2][1][1] * pVoxel0px0py1pz - + weights[0][2][1] * pVoxel0px1py1nz + weights[2][2][1] * + pVoxel0px1py1pz - weights[0][0][2] * pVoxel1px1ny1nz + + weights[2][0][2] * pVoxel1px1ny1pz - weights[0][1][2] * + pVoxel1px0py1nz + weights[2][1][2] * pVoxel1px0py1pz - + weights[0][2][2] * pVoxel1px1py1nz + weights[2][2][2] * + pVoxel1px1py1pz); + + //Note: The above actually give gradients going from low density to high density. + //For our normals we want the the other way around, so we switch the components as we return them. + return Vector3DFloat(-xGrad,-yGrad,-zGrad); + } + //////////////////////////////////////////////////////////////////////////////// + // End of compiler bug workaroumd. + //////////////////////////////////////////////////////////////////////////////// + + //Use the cell bitmasks to generate all the indices needed for that slice + void generateIndicesForSlice(const Array2DUint8& pPreviousBitmask, + const Array2DInt32& m_pPreviousVertexIndicesX, + const Array2DInt32& m_pPreviousVertexIndicesY, + const Array2DInt32& m_pPreviousVertexIndicesZ, + const Array2DInt32& m_pCurrentVertexIndicesX, + const Array2DInt32& m_pCurrentVertexIndicesY); + + //The volume data and a sampler to access it. + VolumeType* m_volData; + typename VolumeType::Sampler m_sampVolume; + + //Holds a position in volume space. + int32_t iXVolSpace; + int32_t iYVolSpace; + int32_t iZVolSpace; + + //Holds a position in region space. + uint32_t uXRegSpace; + uint32_t uYRegSpace; + uint32_t uZRegSpace; + + //Used to return the number of cells in a slice which contain triangles. + uint32_t m_uNoOfOccupiedCells; + + //The surface patch we are currently filling. + SurfaceMesh* m_meshCurrent; + + //Information about the region we are currently processing + Region m_regSizeInVoxels; + Region m_regSizeInCells; + /*Region m_regSizeInVoxelsCropped; + Region m_regSizeInVoxelsUncropped; + Region m_regVolumeCropped;*/ + Region m_regSlicePrevious; + Region m_regSliceCurrent; + + //Our threshold value + //typename VoxelTypeTraits::DensityType m_tThreshold; + typename MarchingCubesController::DensityType m_tThreshold; + + //Used to convert arbitrary voxel types in densities and materials. + //MarchingCubesController m_controller; + Controller m_controller; + }; +} + +#include "PolyVoxCore/MarchingCubesSurfaceExtractor.inl" + +#endif diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl new file mode 100644 index 00000000..e8f3d3ee --- /dev/null +++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl @@ -0,0 +1,624 @@ +/******************************************************************************* +Copyright (c) 2005-2009 David Williams + +This software is provided 'as-is', without any express or implied +warranty. In no event will the authors be held liable for any damages +arising from the use of this software. + +Permission is granted to anyone to use this software for any purpose, +including commercial applications, and to alter it and redistribute it +freely, subject to the following restrictions: + + 1. The origin of this software must not be misrepresented; you must not + claim that you wrote the original software. If you use this software + in a product, an acknowledgment in the product documentation would be + appreciated but is not required. + + 2. Altered source versions must be plainly marked as such, and must not be + misrepresented as being the original software. + + 3. This notice may not be removed or altered from any source + distribution. +*******************************************************************************/ + +namespace PolyVox +{ + template + MarchingCubesSurfaceExtractor::MarchingCubesSurfaceExtractor(VolumeType* volData, Region region, SurfaceMesh* result, Controller controller) + :m_volData(volData) + ,m_sampVolume(volData) + ,m_meshCurrent(result) + ,m_regSizeInVoxels(region) + { + //m_regSizeInVoxels.cropTo(m_volData->getEnclosingRegion()); + m_regSizeInCells = m_regSizeInVoxels; + m_regSizeInCells.setUpperCorner(m_regSizeInCells.getUpperCorner() - Vector3DInt32(1,1,1)); + + m_controller = controller; + m_tThreshold = m_controller.getThreshold(); + } + + template + void MarchingCubesSurfaceExtractor::execute() + { + m_meshCurrent->clear(); + + uint32_t uArrayWidth = m_regSizeInVoxels.getUpperCorner().getX() - m_regSizeInVoxels.getLowerCorner().getX() + 1; + uint32_t uArrayHeight = m_regSizeInVoxels.getUpperCorner().getY() - m_regSizeInVoxels.getLowerCorner().getY() + 1; + uint32_t arraySizes[2]= {uArrayWidth, uArrayHeight}; // Array dimensions + + //For edge indices + Array2DInt32 m_pPreviousVertexIndicesX(arraySizes); + Array2DInt32 m_pPreviousVertexIndicesY(arraySizes); + Array2DInt32 m_pPreviousVertexIndicesZ(arraySizes); + Array2DInt32 m_pCurrentVertexIndicesX(arraySizes); + Array2DInt32 m_pCurrentVertexIndicesY(arraySizes); + Array2DInt32 m_pCurrentVertexIndicesZ(arraySizes); + + Array2DUint8 pPreviousBitmask(arraySizes); + Array2DUint8 pCurrentBitmask(arraySizes); + + //Create a region corresponding to the first slice + m_regSlicePrevious = m_regSizeInVoxels; + Vector3DInt32 v3dUpperCorner = m_regSlicePrevious.getUpperCorner(); + v3dUpperCorner.setZ(m_regSlicePrevious.getLowerCorner().getZ()); //Set the upper z to the lower z to make it one slice thick. + m_regSlicePrevious.setUpperCorner(v3dUpperCorner); + m_regSliceCurrent = m_regSlicePrevious; + + uint32_t uNoOfNonEmptyCellsForSlice0 = 0; + uint32_t uNoOfNonEmptyCellsForSlice1 = 0; + + //Process the first slice (previous slice not available) + computeBitmaskForSlice(pPreviousBitmask, pCurrentBitmask); + uNoOfNonEmptyCellsForSlice1 = m_uNoOfOccupiedCells; + + if(uNoOfNonEmptyCellsForSlice1 != 0) + { + memset(m_pCurrentVertexIndicesX.getRawData(), 0xff, m_pCurrentVertexIndicesX.getNoOfElements() * 4); + memset(m_pCurrentVertexIndicesY.getRawData(), 0xff, m_pCurrentVertexIndicesY.getNoOfElements() * 4); + memset(m_pCurrentVertexIndicesZ.getRawData(), 0xff, m_pCurrentVertexIndicesZ.getNoOfElements() * 4); + generateVerticesForSlice(pCurrentBitmask, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY, m_pCurrentVertexIndicesZ); + } + + std::swap(uNoOfNonEmptyCellsForSlice0, uNoOfNonEmptyCellsForSlice1); + pPreviousBitmask.swap(pCurrentBitmask); + m_pPreviousVertexIndicesX.swap(m_pCurrentVertexIndicesX); + m_pPreviousVertexIndicesY.swap(m_pCurrentVertexIndicesY); + m_pPreviousVertexIndicesZ.swap(m_pCurrentVertexIndicesZ); + + m_regSlicePrevious = m_regSliceCurrent; + m_regSliceCurrent.shift(Vector3DInt32(0,0,1)); + + //Process the other slices (previous slice is available) + for(int32_t uSlice = 1; uSlice <= m_regSizeInVoxels.getUpperCorner().getZ() - m_regSizeInVoxels.getLowerCorner().getZ(); uSlice++) + { + computeBitmaskForSlice(pPreviousBitmask, pCurrentBitmask); + uNoOfNonEmptyCellsForSlice1 = m_uNoOfOccupiedCells; + + if(uNoOfNonEmptyCellsForSlice1 != 0) + { + memset(m_pCurrentVertexIndicesX.getRawData(), 0xff, m_pCurrentVertexIndicesX.getNoOfElements() * 4); + memset(m_pCurrentVertexIndicesY.getRawData(), 0xff, m_pCurrentVertexIndicesY.getNoOfElements() * 4); + memset(m_pCurrentVertexIndicesZ.getRawData(), 0xff, m_pCurrentVertexIndicesZ.getNoOfElements() * 4); + generateVerticesForSlice(pCurrentBitmask, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY, m_pCurrentVertexIndicesZ); + } + + if((uNoOfNonEmptyCellsForSlice0 != 0) || (uNoOfNonEmptyCellsForSlice1 != 0)) + { + generateIndicesForSlice(pPreviousBitmask, m_pPreviousVertexIndicesX, m_pPreviousVertexIndicesY, m_pPreviousVertexIndicesZ, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY); + } + + std::swap(uNoOfNonEmptyCellsForSlice0, uNoOfNonEmptyCellsForSlice1); + pPreviousBitmask.swap(pCurrentBitmask); + m_pPreviousVertexIndicesX.swap(m_pCurrentVertexIndicesX); + m_pPreviousVertexIndicesY.swap(m_pCurrentVertexIndicesY); + m_pPreviousVertexIndicesZ.swap(m_pCurrentVertexIndicesZ); + + m_regSlicePrevious = m_regSliceCurrent; + m_regSliceCurrent.shift(Vector3DInt32(0,0,1)); + } + + m_meshCurrent->m_Region = m_regSizeInVoxels; + + m_meshCurrent->m_vecLodRecords.clear(); + LodRecord lodRecord; + lodRecord.beginIndex = 0; + lodRecord.endIndex = m_meshCurrent->getNoOfIndices(); + m_meshCurrent->m_vecLodRecords.push_back(lodRecord); + } + + template + template + uint32_t MarchingCubesSurfaceExtractor::computeBitmaskForSlice(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask) + { + m_uNoOfOccupiedCells = 0; + + const int32_t iMaxXVolSpace = m_regSliceCurrent.getUpperCorner().getX(); + const int32_t iMaxYVolSpace = m_regSliceCurrent.getUpperCorner().getY(); + + iZVolSpace = m_regSliceCurrent.getLowerCorner().getZ(); + uZRegSpace = iZVolSpace - m_regSizeInVoxels.getLowerCorner().getZ(); + + //Process the lower left corner + iYVolSpace = m_regSliceCurrent.getLowerCorner().getY(); + iXVolSpace = m_regSliceCurrent.getLowerCorner().getX(); + + uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX(); + uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY(); + + m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace); + computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask); + + //Process the edge where x is minimal. + iXVolSpace = m_regSliceCurrent.getLowerCorner().getX(); + m_sampVolume.setPosition(iXVolSpace, m_regSliceCurrent.getLowerCorner().getY(), iZVolSpace); + for(iYVolSpace = m_regSliceCurrent.getLowerCorner().getY() + 1; iYVolSpace <= iMaxYVolSpace; iYVolSpace++) + { + uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX(); + uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY(); + + m_sampVolume.movePositiveY(); + + computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask); + } + + //Process the edge where y is minimal. + iYVolSpace = m_regSliceCurrent.getLowerCorner().getY(); + m_sampVolume.setPosition(m_regSliceCurrent.getLowerCorner().getX(), iYVolSpace, iZVolSpace); + for(iXVolSpace = m_regSliceCurrent.getLowerCorner().getX() + 1; iXVolSpace <= iMaxXVolSpace; iXVolSpace++) + { + uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX(); + uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY(); + + m_sampVolume.movePositiveX(); + + computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask); + } + + //Process all remaining elemnents of the slice. In this case, previous x and y values are always available + for(iYVolSpace = m_regSliceCurrent.getLowerCorner().getY() + 1; iYVolSpace <= iMaxYVolSpace; iYVolSpace++) + { + m_sampVolume.setPosition(m_regSliceCurrent.getLowerCorner().getX(), iYVolSpace, iZVolSpace); + for(iXVolSpace = m_regSliceCurrent.getLowerCorner().getX() + 1; iXVolSpace <= iMaxXVolSpace; iXVolSpace++) + { + uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX(); + uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY(); + + m_sampVolume.movePositiveX(); + + computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask); + } + } + + return m_uNoOfOccupiedCells; + } + + template + template + void MarchingCubesSurfaceExtractor::computeBitmaskForCell(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask) + { + uint8_t iCubeIndex = 0; + + typename VolumeType::VoxelType v000; + typename VolumeType::VoxelType v100; + typename VolumeType::VoxelType v010; + typename VolumeType::VoxelType v110; + typename VolumeType::VoxelType v001; + typename VolumeType::VoxelType v101; + typename VolumeType::VoxelType v011; + typename VolumeType::VoxelType v111; + + if(isPrevZAvail) + { + if(isPrevYAvail) + { + if(isPrevXAvail) + { + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //z + uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace]; + iPreviousCubeIndexZ >>= 4; + + //y + uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1]; + iPreviousCubeIndexY &= 192; //192 = 128 + 64 + iPreviousCubeIndexY >>= 2; + + //x + uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace]; + iPreviousCubeIndexX &= 128; + iPreviousCubeIndexX >>= 1; + + iCubeIndex = iPreviousCubeIndexX | iPreviousCubeIndexY | iPreviousCubeIndexZ; + + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + else //previous X not available + { + v011 = m_sampVolume.peekVoxel0px1py1pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //z + uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace]; + iPreviousCubeIndexZ >>= 4; + + //y + uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1]; + iPreviousCubeIndexY &= 192; //192 = 128 + 64 + iPreviousCubeIndexY >>= 2; + + iCubeIndex = iPreviousCubeIndexY | iPreviousCubeIndexZ; + + if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + } + else //previous Y not available + { + if(isPrevXAvail) + { + v101 = m_sampVolume.peekVoxel1px0py1pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //z + uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace]; + iPreviousCubeIndexZ >>= 4; + + //x + uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace]; + iPreviousCubeIndexX &= 160; //160 = 128+32 + iPreviousCubeIndexX >>= 1; + + iCubeIndex = iPreviousCubeIndexX | iPreviousCubeIndexZ; + + if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + else //previous X not available + { + v001 = m_sampVolume.peekVoxel0px0py1pz(); + v101 = m_sampVolume.peekVoxel1px0py1pz(); + v011 = m_sampVolume.peekVoxel0px1py1pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //z + uint8_t iPreviousCubeIndexZ = pPreviousBitmask[uXRegSpace][uYRegSpace]; + iCubeIndex = iPreviousCubeIndexZ >> 4; + + if (m_controller.convertToDensity(v001) < m_tThreshold) iCubeIndex |= 16; + if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32; + if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + } + } + else //previous Z not available + { + if(isPrevYAvail) + { + if(isPrevXAvail) + { + v110 = m_sampVolume.peekVoxel1px1py0pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //y + uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1]; + iPreviousCubeIndexY &= 204; //204 = 128+64+8+4 + iPreviousCubeIndexY >>= 2; + + //x + uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace]; + iPreviousCubeIndexX &= 170; //170 = 128+32+8+2 + iPreviousCubeIndexX >>= 1; + + iCubeIndex = iPreviousCubeIndexX | iPreviousCubeIndexY; + + if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + else //previous X not available + { + v010 = m_sampVolume.peekVoxel0px1py0pz(); + v110 = m_sampVolume.peekVoxel1px1py0pz(); + + v011 = m_sampVolume.peekVoxel0px1py1pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //y + uint8_t iPreviousCubeIndexY = pCurrentBitmask[uXRegSpace][uYRegSpace-1]; + iPreviousCubeIndexY &= 204; //204 = 128+64+8+4 + iPreviousCubeIndexY >>= 2; + + iCubeIndex = iPreviousCubeIndexY; + + if (m_controller.convertToDensity(v010) < m_tThreshold) iCubeIndex |= 4; + if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8; + if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + } + else //previous Y not available + { + if(isPrevXAvail) + { + v100 = m_sampVolume.peekVoxel1px0py0pz(); + v110 = m_sampVolume.peekVoxel1px1py0pz(); + + v101 = m_sampVolume.peekVoxel1px0py1pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + //x + uint8_t iPreviousCubeIndexX = pCurrentBitmask[uXRegSpace-1][uYRegSpace]; + iPreviousCubeIndexX &= 170; //170 = 128+32+8+2 + iPreviousCubeIndexX >>= 1; + + iCubeIndex = iPreviousCubeIndexX; + + if (m_controller.convertToDensity(v100) < m_tThreshold) iCubeIndex |= 2; + if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8; + if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + else //previous X not available + { + v000 = m_sampVolume.getVoxel(); + v100 = m_sampVolume.peekVoxel1px0py0pz(); + v010 = m_sampVolume.peekVoxel0px1py0pz(); + v110 = m_sampVolume.peekVoxel1px1py0pz(); + + v001 = m_sampVolume.peekVoxel0px0py1pz(); + v101 = m_sampVolume.peekVoxel1px0py1pz(); + v011 = m_sampVolume.peekVoxel0px1py1pz(); + v111 = m_sampVolume.peekVoxel1px1py1pz(); + + if (m_controller.convertToDensity(v000) < m_tThreshold) iCubeIndex |= 1; + if (m_controller.convertToDensity(v100) < m_tThreshold) iCubeIndex |= 2; + if (m_controller.convertToDensity(v010) < m_tThreshold) iCubeIndex |= 4; + if (m_controller.convertToDensity(v110) < m_tThreshold) iCubeIndex |= 8; + if (m_controller.convertToDensity(v001) < m_tThreshold) iCubeIndex |= 16; + if (m_controller.convertToDensity(v101) < m_tThreshold) iCubeIndex |= 32; + if (m_controller.convertToDensity(v011) < m_tThreshold) iCubeIndex |= 64; + if (m_controller.convertToDensity(v111) < m_tThreshold) iCubeIndex |= 128; + } + } + } + + //Save the bitmask + pCurrentBitmask[uXRegSpace][iYVolSpace- m_regSizeInVoxels.getLowerCorner().getY()] = iCubeIndex; + + if(edgeTable[iCubeIndex] != 0) + { + ++m_uNoOfOccupiedCells; + } + } + + template + void MarchingCubesSurfaceExtractor::generateVerticesForSlice(const Array2DUint8& pCurrentBitmask, + Array2DInt32& m_pCurrentVertexIndicesX, + Array2DInt32& m_pCurrentVertexIndicesY, + Array2DInt32& m_pCurrentVertexIndicesZ) + { + int32_t iZVolSpace = m_regSliceCurrent.getLowerCorner().getZ(); + const uint32_t uZRegSpace = iZVolSpace - m_regSizeInVoxels.getLowerCorner().getZ(); + + //Iterate over each cell in the region + for(int32_t iYVolSpace = m_regSliceCurrent.getLowerCorner().getY(); iYVolSpace <= m_regSliceCurrent.getUpperCorner().getY(); iYVolSpace++) + { + const uint32_t uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY(); + + for(int32_t iXVolSpace = m_regSliceCurrent.getLowerCorner().getX(); iXVolSpace <= m_regSliceCurrent.getUpperCorner().getX(); iXVolSpace++) + { + //Current position + const uint32_t uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX(); + + //Determine the index into the edge table which tells us which vertices are inside of the surface + uint8_t iCubeIndex = pCurrentBitmask[uXRegSpace][uYRegSpace]; + + /* Cube is entirely in/out of the surface */ + if (edgeTable[iCubeIndex] == 0) + { + continue; + } + + //Check whether the generated vertex will lie on the edge of the region + + + m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace); + const typename VolumeType::VoxelType v000 = m_sampVolume.getVoxel(); + const Vector3DFloat n000 = computeCentralDifferenceGradient(m_sampVolume); + + /* Find the vertices where the surface intersects the cube */ + if (edgeTable[iCubeIndex] & 1) + { + m_sampVolume.movePositiveX(); + const typename VolumeType::VoxelType v100 = m_sampVolume.getVoxel(); + const Vector3DFloat n100 = computeCentralDifferenceGradient(m_sampVolume); + + float fInterp = static_cast(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast(m_controller.convertToDensity(v100) - m_controller.convertToDensity(v000)); + + const Vector3DFloat v3dPosition(static_cast(iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()) + fInterp, static_cast(iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()), static_cast(iZVolSpace - m_regSizeInCells.getLowerCorner().getZ())); + + Vector3DFloat v3dNormal = (n100*fInterp) + (n000*(1-fInterp)); + v3dNormal.normalise(); + + //Choose one of the two materials to use for the vertex (we don't interpolate as interpolation of + //material IDs does not make sense). We take the largest, so that if we are working on a material-only + //volume we get the one which is non-zero. Both materials can be non-zero if our volume has a density component. + typename MarchingCubesController::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); + typename MarchingCubesController::MaterialType uMaterial100 = m_controller.convertToMaterial(v100); + typename MarchingCubesController::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial100); + + PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, static_cast(uMaterial)); + uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex); + m_pCurrentVertexIndicesX[iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()][iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()] = uLastVertexIndex; + + m_sampVolume.moveNegativeX(); + } + if (edgeTable[iCubeIndex] & 8) + { + m_sampVolume.movePositiveY(); + const typename VolumeType::VoxelType v010 = m_sampVolume.getVoxel(); + const Vector3DFloat n010 = computeCentralDifferenceGradient(m_sampVolume); + + float fInterp = static_cast(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast(m_controller.convertToDensity(v010) - m_controller.convertToDensity(v000)); + + const Vector3DFloat v3dPosition(static_cast(iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()), static_cast(iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()) + fInterp, static_cast(iZVolSpace - m_regSizeInVoxels.getLowerCorner().getZ())); + + Vector3DFloat v3dNormal = (n010*fInterp) + (n000*(1-fInterp)); + v3dNormal.normalise(); + + //Choose one of the two materials to use for the vertex (we don't interpolate as interpolation of + //material IDs does not make sense). We take the largest, so that if we are working on a material-only + //volume we get the one which is non-zero. Both materials can be non-zero if our volume has a density component. + typename MarchingCubesController::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); + typename MarchingCubesController::MaterialType uMaterial010 = m_controller.convertToMaterial(v010); + typename MarchingCubesController::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial010); + + PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, static_cast(uMaterial)); + uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex); + m_pCurrentVertexIndicesY[iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()][iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()] = uLastVertexIndex; + + m_sampVolume.moveNegativeY(); + } + if (edgeTable[iCubeIndex] & 256) + { + m_sampVolume.movePositiveZ(); + const typename VolumeType::VoxelType v001 = m_sampVolume.getVoxel(); + const Vector3DFloat n001 = computeCentralDifferenceGradient(m_sampVolume); + + float fInterp = static_cast(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast(m_controller.convertToDensity(v001) - m_controller.convertToDensity(v000)); + + const Vector3DFloat v3dPosition(static_cast(iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()), static_cast(iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()), static_cast(iZVolSpace - m_regSizeInVoxels.getLowerCorner().getZ()) + fInterp); + + Vector3DFloat v3dNormal = (n001*fInterp) + (n000*(1-fInterp)); + v3dNormal.normalise(); + + //Choose one of the two materials to use for the vertex (we don't interpolate as interpolation of + //material IDs does not make sense). We take the largest, so that if we are working on a material-only + //volume we get the one which is non-zero. Both materials can be non-zero if our volume has a density component. + typename MarchingCubesController::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); + typename MarchingCubesController::MaterialType uMaterial001 = m_controller.convertToMaterial(v001); + typename MarchingCubesController::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial001); + + PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, static_cast(uMaterial)); + uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex); + m_pCurrentVertexIndicesZ[iXVolSpace - m_regSizeInVoxels.getLowerCorner().getX()][iYVolSpace - m_regSizeInVoxels.getLowerCorner().getY()] = uLastVertexIndex; + + m_sampVolume.moveNegativeZ(); + } + }//For each cell + } + } + + template + void MarchingCubesSurfaceExtractor::generateIndicesForSlice(const Array2DUint8& pPreviousBitmask, + const Array2DInt32& m_pPreviousVertexIndicesX, + const Array2DInt32& m_pPreviousVertexIndicesY, + const Array2DInt32& m_pPreviousVertexIndicesZ, + const Array2DInt32& m_pCurrentVertexIndicesX, + const Array2DInt32& m_pCurrentVertexIndicesY) + { + int32_t indlist[12]; + for(int i = 0; i < 12; i++) + { + indlist[i] = -1; + } + + for(int32_t iYVolSpace = m_regSlicePrevious.getLowerCorner().getY(); iYVolSpace <= m_regSizeInCells.getUpperCorner().getY(); iYVolSpace++) + { + for(int32_t iXVolSpace = m_regSlicePrevious.getLowerCorner().getX(); iXVolSpace <= m_regSizeInCells.getUpperCorner().getX(); iXVolSpace++) + { + int32_t iZVolSpace = m_regSlicePrevious.getLowerCorner().getZ(); + m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace); + + //Current position + const uint32_t uXRegSpace = m_sampVolume.getPosition().getX() - m_regSizeInVoxels.getLowerCorner().getX(); + const uint32_t uYRegSpace = m_sampVolume.getPosition().getY() - m_regSizeInVoxels.getLowerCorner().getY(); + + //Determine the index into the edge table which tells us which vertices are inside of the surface + uint8_t iCubeIndex = pPreviousBitmask[uXRegSpace][uYRegSpace]; + + /* Cube is entirely in/out of the surface */ + if (edgeTable[iCubeIndex] == 0) + { + continue; + } + + /* Find the vertices where the surface intersects the cube */ + if (edgeTable[iCubeIndex] & 1) + { + indlist[0] = m_pPreviousVertexIndicesX[uXRegSpace][uYRegSpace]; + //assert(indlist[0] != -1); + } + if (edgeTable[iCubeIndex] & 2) + { + indlist[1] = m_pPreviousVertexIndicesY[uXRegSpace+1][uYRegSpace]; + //assert(indlist[1] != -1); + } + if (edgeTable[iCubeIndex] & 4) + { + indlist[2] = m_pPreviousVertexIndicesX[uXRegSpace][uYRegSpace+1]; + //assert(indlist[2] != -1); + } + if (edgeTable[iCubeIndex] & 8) + { + indlist[3] = m_pPreviousVertexIndicesY[uXRegSpace][uYRegSpace]; + //assert(indlist[3] != -1); + } + if (edgeTable[iCubeIndex] & 16) + { + indlist[4] = m_pCurrentVertexIndicesX[uXRegSpace][uYRegSpace]; + //assert(indlist[4] != -1); + } + if (edgeTable[iCubeIndex] & 32) + { + indlist[5] = m_pCurrentVertexIndicesY[uXRegSpace+1][uYRegSpace]; + //assert(indlist[5] != -1); + } + if (edgeTable[iCubeIndex] & 64) + { + indlist[6] = m_pCurrentVertexIndicesX[uXRegSpace][uYRegSpace+1]; + //assert(indlist[6] != -1); + } + if (edgeTable[iCubeIndex] & 128) + { + indlist[7] = m_pCurrentVertexIndicesY[uXRegSpace][uYRegSpace]; + //assert(indlist[7] != -1); + } + if (edgeTable[iCubeIndex] & 256) + { + indlist[8] = m_pPreviousVertexIndicesZ[uXRegSpace][uYRegSpace]; + //assert(indlist[8] != -1); + } + if (edgeTable[iCubeIndex] & 512) + { + indlist[9] = m_pPreviousVertexIndicesZ[uXRegSpace+1][uYRegSpace]; + //assert(indlist[9] != -1); + } + if (edgeTable[iCubeIndex] & 1024) + { + indlist[10] = m_pPreviousVertexIndicesZ[uXRegSpace+1][uYRegSpace+1]; + //assert(indlist[10] != -1); + } + if (edgeTable[iCubeIndex] & 2048) + { + indlist[11] = m_pPreviousVertexIndicesZ[uXRegSpace][uYRegSpace+1]; + //assert(indlist[11] != -1); + } + + for (int i=0;triTable[iCubeIndex][i]!=-1;i+=3) + { + int32_t ind0 = indlist[triTable[iCubeIndex][i ]]; + int32_t ind1 = indlist[triTable[iCubeIndex][i+1]]; + int32_t ind2 = indlist[triTable[iCubeIndex][i+2]]; + + if((ind0 != -1) && (ind1 != -1) && (ind2 != -1)) + { + m_meshCurrent->addTriangle(ind0, ind1, ind2); + } + }//For each triangle + }//For each cell + } + } +}