diff --git a/library/PolyVoxCore/include/PolyVoxCore/DefaultMarchingCubesController.h b/library/PolyVoxCore/include/PolyVoxCore/DefaultMarchingCubesController.h index 43aa0956..fd0d00ef 100644 --- a/library/PolyVoxCore/include/PolyVoxCore/DefaultMarchingCubesController.h +++ b/library/PolyVoxCore/include/PolyVoxCore/DefaultMarchingCubesController.h @@ -30,34 +30,34 @@ freely, subject to the following restrictions: namespace PolyVox { - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// This class provides a default implementation of a controller for the MarchingCubesSurfaceExtractor. It controls the behaviour of the - /// MarchingCubesSurfaceExtractor and provides the required properties from the underlying voxel type. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// PolyVox does not enforce any requirements regarding what data must be present in a voxel, and instead allows any primitive or user-defined - /// type to be used. However, the Marching Cubes algorithm does have some requirents about the underlying data in that conceptually it operates - /// on a density field. In addition, the PolyVox implementation of the Marching Cubes algorithm also understands the idea of each voxel - /// having a material which is copied into the vertex data. - /// - /// Because we want the MarchingCubesSurfaceExtractor to work on any voxel type, we use a Marching Cubes controller (passed as - /// a parameter of the MarchingCubesSurfaceExtractor) to expose the required properties. This parameter defaults to the DefaultMarchingCubesController. - /// The main implementation of this class is designed to work with primitives data types, and the class is also specialised for the Material, - /// Density and MaterialdensityPair classes. - /// - /// If you create a custom class for your voxel data then you probably want to include a specialisation of DefaultMarchingCubesController, - /// though you don't have to if you don't want to use the Marching Cubes algorithm or if you prefer to define a seperate Marching Cubes controller - /// and pass it as an explicit parameter (rather than relying on the default). - /// - /// For primitive types, the DefaultMarchingCubesController considers the value of the voxel to represent it's density and just returns a constant - /// for the material. So you can, for example, run the MarchingCubesSurfaceExtractor on a volume of floats or ints. - /// - /// It is possible to customise the behaviour of the controller by providing a threshold value through the constructor. The extracted surface - /// will pass through the density value specified by the threshold, and so you should make sure that the threshold value you choose is between - /// the minimum and maximum values found in your volume data. By default it is in the middle of the representable range of the underlying type. - /// - /// \sa MarchingCubesSurfaceExtractor - /// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** + * This class provides a default implementation of a controller for the MarchingCubesSurfaceExtractor. It controls the behaviour of the + * MarchingCubesSurfaceExtractor and provides the required properties from the underlying voxel type. + * + * PolyVox does not enforce any requirements regarding what data must be present in a voxel, and instead allows any primitive or user-defined + * type to be used. However, the Marching Cubes algorithm does have some requirents about the underlying data in that conceptually it operates + * on a density field. In addition, the PolyVox implementation of the Marching Cubes algorithm also understands the idea of each voxel + * having a material which is copied into the vertex data. + * + * Because we want the MarchingCubesSurfaceExtractor to work on any voxel type, we use a Marching Cubes controller (passed as + * a parameter of the MarchingCubesSurfaceExtractor) to expose the required properties. This parameter defaults to the DefaultMarchingCubesController. + * The main implementation of this class is designed to work with primitives data types, and the class is also specialised for the Material, + * Density and MaterialdensityPair classes. + * + * If you create a custom class for your voxel data then you probably want to include a specialisation of DefaultMarchingCubesController, + * though you don't have to if you don't want to use the Marching Cubes algorithm or if you prefer to define a seperate Marching Cubes controller + * and pass it as an explicit parameter (rather than relying on the default). + * + * For primitive types, the DefaultMarchingCubesController considers the value of the voxel to represent it's density and just returns a constant + * for the material. So you can, for example, run the MarchingCubesSurfaceExtractor on a volume of floats or ints. + * + * It is possible to customise the behaviour of the controller by providing a threshold value through the constructor. The extracted surface + * will pass through the density value specified by the threshold, and so you should make sure that the threshold value you choose is between + * the minimum and maximum values found in your volume data. By default it is in the middle of the representable range of the underlying type. + * + * \sa MarchingCubesSurfaceExtractor + * + */ template class DefaultMarchingCubesController { @@ -69,65 +69,79 @@ namespace PolyVox /// but this is not really desirable on modern hardware. We'll probably come back to material representation in the future. typedef float MaterialType; - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// Constructor - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// This version of the constructor takes no parameters and sets the threshold to the middle of the representable range of the underlying type. - /// For example, if the voxel type is 'uint8_t' then the representable range is 0-255, and the threshold will be set to 127. On the other hand, - /// if the voxel type is 'float' then the representable range is -FLT_MAX to FLT_MAX and the threshold will be set to zero. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** + * Constructor + * + * This version of the constructor takes no parameters and sets the threshold to the middle of the representable range of the underlying type. + * For example, if the voxel type is 'uint8_t' then the representable range is 0-255, and the threshold will be set to 127. On the other hand, + * if the voxel type is 'float' then the representable range is -FLT_MAX to FLT_MAX and the threshold will be set to zero. + */ DefaultMarchingCubesController(void) - { - m_tThreshold = ((std::numeric_limits::min)() + (std::numeric_limits::max)()) / 2; + :m_tThreshold(((std::numeric_limits::min)() + (std::numeric_limits::max)()) / 2) + ,m_eWrapMode(WrapModes::Border) + ,m_tBorder(VoxelType(0)) + { } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// Constructor - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// This version of the constructor allows you to set a custom threshold. - /// \param tThreshold The threshold to use. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - DefaultMarchingCubesController(DensityType tThreshold) - { - m_tThreshold = tThreshold; - } - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// Converts the underlying voxel type into a density value. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// The default implementation of this function just returns the voxel type directly and is suitable for primitives types. Specialisations of - /// this class can modify this behaviour. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** + * Converts the underlying voxel type into a density value. + * + * The default implementation of this function just returns the voxel type directly and is suitable for primitives types. Specialisations of + * this class can modify this behaviour. + */ DensityType convertToDensity(VoxelType voxel) { return voxel; } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// Converts the underlying voxel type into a material value. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// The default implementation of this function just returns the constant '1'. There's not much else it can do, as it needs to work with primitive - /// types and the actual value of the type is already being considered to be the density. Specialisations of this class can modify this behaviour. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - MaterialType convertToMaterial(VoxelType voxel) + /** + * Converts the underlying voxel type into a material value. + * + * The default implementation of this function just returns the constant '1'. There's not much else it can do, as it needs to work with primitive + * types and the actual value of the type is already being considered to be the density. Specialisations of this class can modify this behaviour. + */ + MaterialType convertToMaterial(VoxelType /*voxel*/) { return 1; } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// Returns the density value which was passed to the constructor. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// As mentioned in the class description, the extracted surface will pass through the density value specified by the threshold, and so you - /// should make sure that the threshold value you choose is between the minimum and maximum values found in your volume data. By default it - ///is in the middle of the representable range of the underlying type. - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + VoxelType getBorderValue(void) + { + return m_tBorder; + } + + /** + * Returns the density value which was passed to the constructor. + * + * As mentioned in the class description, the extracted surface will pass through the density value specified by the threshold, and so you + * should make sure that the threshold value you choose is between the minimum and maximum values found in your volume data. By default it + * is in the middle of the representable range of the underlying type. + */ DensityType getThreshold(void) { return m_tThreshold; } + WrapMode getWrapMode(void) + { + return m_eWrapMode; + } + + void setThreshold(DensityType tThreshold) + { + m_tThreshold = tThreshold; + } + + void setWrapMode(WrapMode eWrapMode, VoxelType tBorder = VoxelType(0)) + { + m_eWrapMode = eWrapMode; + m_tBorder = tBorder; + } + private: DensityType m_tThreshold; + WrapMode m_eWrapMode; + VoxelType m_tBorder; }; } diff --git a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl index e9761d1e..be2a6674 100644 --- a/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl +++ b/library/PolyVoxCore/include/PolyVoxCore/MarchingCubesSurfaceExtractor.inl @@ -1,626 +1,628 @@ -/******************************************************************************* -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(); - - //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 = computeSobelGradient(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 = computeSobelGradient(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 Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); - typename Controller::MaterialType uMaterial100 = m_controller.convertToMaterial(v100); - //typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial100); - typename Controller::MaterialType uMaterial = m_controller.blendMaterials(uMaterial000, uMaterial100, fInterp); - - PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, 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 = computeSobelGradient(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 Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); - typename Controller::MaterialType uMaterial010 = m_controller.convertToMaterial(v010); - //typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial010); - typename Controller::MaterialType uMaterial = m_controller.blendMaterials(uMaterial000, uMaterial010, fInterp); - - PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, 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 = computeSobelGradient(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 Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); - typename Controller::MaterialType uMaterial001 = m_controller.convertToMaterial(v001); - //typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial001); - typename Controller::MaterialType uMaterial = m_controller.blendMaterials(uMaterial000, uMaterial001, fInterp); - - PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, 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 - } - } -} +/******************************************************************************* +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(); + + m_sampVolume.setWrapMode(m_controller.getWrapMode(), m_controller.getBorderValue()); + } + + 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(); + + //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 = computeSobelGradient(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 = computeSobelGradient(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 Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); + typename Controller::MaterialType uMaterial100 = m_controller.convertToMaterial(v100); + //typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial100); + typename Controller::MaterialType uMaterial = m_controller.blendMaterials(uMaterial000, uMaterial100, fInterp); + + PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, 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 = computeSobelGradient(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 Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); + typename Controller::MaterialType uMaterial010 = m_controller.convertToMaterial(v010); + //typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial010); + typename Controller::MaterialType uMaterial = m_controller.blendMaterials(uMaterial000, uMaterial010, fInterp); + + PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, 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 = computeSobelGradient(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 Controller::MaterialType uMaterial000 = m_controller.convertToMaterial(v000); + typename Controller::MaterialType uMaterial001 = m_controller.convertToMaterial(v001); + //typename Controller::MaterialType uMaterial = (std::max)(uMaterial000, uMaterial001); + typename Controller::MaterialType uMaterial = m_controller.blendMaterials(uMaterial000, uMaterial001, fInterp); + + PositionMaterialNormal surfaceVertex(v3dPosition, v3dNormal, 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 + } + } +}