diff --git a/include/PolyVox/MarchingCubesSurfaceExtractor.h b/include/PolyVox/MarchingCubesSurfaceExtractor.h index efc6eebf..0bd2df59 100644 --- a/include/PolyVox/MarchingCubesSurfaceExtractor.h +++ b/include/PolyVox/MarchingCubesSurfaceExtractor.h @@ -148,184 +148,12 @@ namespace PolyVox return result; } - /// Do not use this class directly. Use the 'extractMarchingCubesSurface' function instead (see examples). - template< typename VolumeType, typename MeshType, typename ControllerType> - class MarchingCubesSurfaceExtractor - { - public: - MarchingCubesSurfaceExtractor(VolumeType* volData, Region region, MeshType* result, ControllerType 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, uint32_t uXRegSpace, uint32_t uYRegSpace); - - //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; - - //Used to return the number of cells in a slice which contain triangles. - uint32_t m_uNoOfOccupiedCells; - - //The surface patch we are currently filling. - MeshType* 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; - - //Used to convert arbitrary voxel types in densities and materials. - ControllerType m_controller; - - //Our threshold value - typename ControllerType::DensityType m_tThreshold; - }; - // This version of the function performs the extraction into a user-provided mesh rather than allocating a mesh automatically. // There are a few reasons why this might be useful to more advanced users: // // 1. It leaves the user in control of memory allocation and would allow them to implement e.g. a mesh pooling system. // 2. The user-provided mesh could have a different index type (e.g. 16-bit indices) to reduce memory usage. - // 3. The user could provide a custom mesh class, e.g a thin wrapper around an openGL VBO to allow direct writing into this structure. + // 3. The user could provide a custom mesh class, e.g a thin wrapper around an OpenGL VBO to allow direct writing into this structure. // // We don't provide a default MeshType here. If the user doesn't want to provide a MeshType then it probably makes // more sense to use the other variant of this function where the mesh is a return value rather than a parameter. @@ -334,11 +162,7 @@ namespace PolyVox // are provided (would the third parameter be a controller or a mesh?). It seems this can be fixed by using enable_if/static_assert to emulate concepts, // but this is relatively complex and I haven't done it yet. Could always add it later as another overload. template< typename VolumeType, typename MeshType, typename ControllerType = DefaultMarchingCubesController > - void extractMarchingCubesMeshCustom(VolumeType* volData, Region region, MeshType* result, ControllerType controller = ControllerType()) - { - MarchingCubesSurfaceExtractor extractor(volData, region, result, controller); - extractor.execute(); - } + void extractMarchingCubesMeshCustom(VolumeType* volData, Region region, MeshType* result, ControllerType controller = ControllerType()); template< typename VolumeType, typename ControllerType = DefaultMarchingCubesController > Mesh > extractMarchingCubesMesh(VolumeType* volData, Region region, ControllerType controller = ControllerType()) diff --git a/include/PolyVox/MarchingCubesSurfaceExtractor.inl b/include/PolyVox/MarchingCubesSurfaceExtractor.inl index cef3ea13..e349e726 100644 --- a/include/PolyVox/MarchingCubesSurfaceExtractor.inl +++ b/include/PolyVox/MarchingCubesSurfaceExtractor.inl @@ -25,611 +25,333 @@ freely, subject to the following restrictions: namespace PolyVox { - template - MarchingCubesSurfaceExtractor::MarchingCubesSurfaceExtractor(VolumeType* volData, Region region, MeshType* result, ControllerType controller) - :m_volData(volData) - ,m_sampVolume(volData) - ,m_meshCurrent(result) - ,m_regSizeInVoxels(region) - ,m_controller(controller) - ,m_tThreshold(m_controller.getThreshold()) + template< typename Sampler, typename ControllerType> + Vector3DFloat computeCentralDifferenceGradient(const Sampler& volIter, ControllerType& controller) { - POLYVOX_THROW_IF(m_meshCurrent == nullptr, std::invalid_argument, "Provided mesh cannot be null"); - //m_regSizeInVoxels.cropTo(m_volData->getEnclosingRegion()); - m_regSizeInCells = m_regSizeInVoxels; - m_regSizeInCells.setUpperCorner(m_regSizeInCells.getUpperCorner() - Vector3DInt32(1,1,1)); + //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(controller.convertToDensity(volIter.peekVoxel1nx0py0pz())); + float voxel1px = static_cast(controller.convertToDensity(volIter.peekVoxel1px0py0pz())); + + float voxel1ny = static_cast(controller.convertToDensity(volIter.peekVoxel0px1ny0pz())); + float voxel1py = static_cast(controller.convertToDensity(volIter.peekVoxel0px1py0pz())); + + float voxel1nz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py1nz())); + float voxel1pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py1pz())); + + return Vector3DFloat + ( + voxel1nx - voxel1px, + voxel1ny - voxel1py, + voxel1nz - voxel1pz + ); } - template - void MarchingCubesSurfaceExtractor::execute() + template< typename VolumeType, typename MeshType, typename ControllerType > + void extractMarchingCubesMeshCustom(VolumeType* volData, Region region, MeshType* result, ControllerType controller) { + // Validate parameters + POLYVOX_THROW_IF(volData == nullptr, std::invalid_argument, "Provided volume cannot be null"); + POLYVOX_THROW_IF(result == nullptr, std::invalid_argument, "Provided mesh cannot be null"); + + // For profiling this function Timer timer; - m_meshCurrent->clear(); + result->clear(); - const uint32_t uArrayWidth = m_regSizeInVoxels.getUpperX() - m_regSizeInVoxels.getLowerX() + 1; - const uint32_t uArrayHeight = m_regSizeInVoxels.getUpperY() - m_regSizeInVoxels.getLowerY() + 1; + // Store some commonly used values for performance and convienience + const uint32_t uRegionWidthInVoxels = region.getWidthInVoxels(); + const uint32_t uRegionHeightInVoxels = region.getHeightInVoxels(); + const uint32_t uRegionDepthInVoxels = region.getDepthInVoxels(); - //For edge indices - Array2DInt32 m_pPreviousVertexIndicesX(uArrayWidth, uArrayHeight); - Array2DInt32 m_pPreviousVertexIndicesY(uArrayWidth, uArrayHeight); - Array2DInt32 m_pPreviousVertexIndicesZ(uArrayWidth, uArrayHeight); - Array2DInt32 m_pCurrentVertexIndicesX(uArrayWidth, uArrayHeight); - Array2DInt32 m_pCurrentVertexIndicesY(uArrayWidth, uArrayHeight); - Array2DInt32 m_pCurrentVertexIndicesZ(uArrayWidth, uArrayHeight); + typename ControllerType::DensityType tThreshold = controller.getThreshold(); - Array2DUint8 pPreviousBitmask(uArrayWidth, uArrayHeight); - Array2DUint8 pCurrentBitmask(uArrayWidth, uArrayHeight); + // A naive implemetation of Marching Cubes might sample the eight corner voxels of every cell to determine the cell index. + // However, when processing the cells sequentially we cn observe that many of the voxels are shared with previous adjacent + // cells, and so we can obtain these by careful bit-shifting. These variables keep track of previous cells for this purpose. + // We don't clear the arrays because the algorithm ensures that we only read from elements we have previously written to. + uint8_t uPreviousCellIndex = 0; + Array1DUint8 pPreviousRowCellIndices(uRegionWidthInVoxels); + Array2DUint8 pPreviousSliceCellIndices(uRegionWidthInVoxels, uRegionHeightInVoxels); - //Create a region corresponding to the first slice - m_regSlicePrevious = m_regSizeInVoxels; - Vector3DInt32 v3dUpperCorner = m_regSlicePrevious.getUpperCorner(); - v3dUpperCorner.setZ(m_regSlicePrevious.getLowerZ()); //Set the upper z to the lower z to make it one slice thick. - m_regSlicePrevious.setUpperCorner(v3dUpperCorner); - m_regSliceCurrent = m_regSlicePrevious; + // A given vertex may be shared by multiple triangles, so we need to keep track of the indices into the vertex array. + // We don't clear the arrays because the algorithm ensures that we only read from elements we have previously written to. + Array<2, Vector3DInt32> pIndices(uRegionWidthInVoxels, uRegionHeightInVoxels); + Array<2, Vector3DInt32> pPreviousIndices(uRegionWidthInVoxels, uRegionHeightInVoxels); - uint32_t uNoOfNonEmptyCellsForSlice0 = 0; - uint32_t uNoOfNonEmptyCellsForSlice1 = 0; + // A sampler pointing at the beginning of the region, which gets incremented to always point at the beginning of a slice. + typename VolumeType::Sampler startOfSlice(volData); + startOfSlice.setPosition(region.getLowerX(), region.getLowerY(), region.getLowerZ()); - //Process the first slice (previous slice not available) - computeBitmaskForSlice(pPreviousBitmask, pCurrentBitmask); - uNoOfNonEmptyCellsForSlice1 = m_uNoOfOccupiedCells; - - if(uNoOfNonEmptyCellsForSlice1 != 0) + for (uint32_t uZRegSpace = 0; uZRegSpace < uRegionDepthInVoxels; uZRegSpace++) { - 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); - } + // A sampler pointing at the beginning of the slice, which gets incremented to always point at the beginning of a row. + typename VolumeType::Sampler startOfRow = startOfSlice; - 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.getUpperZ() - m_regSizeInVoxels.getLowerZ(); uSlice++) - { - computeBitmaskForSlice(pPreviousBitmask, pCurrentBitmask); - uNoOfNonEmptyCellsForSlice1 = m_uNoOfOccupiedCells; - - if(uNoOfNonEmptyCellsForSlice1 != 0) + for (uint32_t uYRegSpace = 0; uYRegSpace < uRegionHeightInVoxels; uYRegSpace++) { - 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); - } + // Copying a sampler which is already pointing at the correct location seems (slightly) faster than + // calling setPosition(). Therefore we make use of 'startOfRow' and 'startOfSlice' to reset the sampler. + typename VolumeType::Sampler sampler = startOfRow; - if((uNoOfNonEmptyCellsForSlice0 != 0) || (uNoOfNonEmptyCellsForSlice1 != 0)) - { - generateIndicesForSlice(pPreviousBitmask, m_pPreviousVertexIndicesX, m_pPreviousVertexIndicesY, m_pPreviousVertexIndicesZ, m_pCurrentVertexIndicesX, m_pCurrentVertexIndicesY); - } + for (uint32_t uXRegSpace = 0; uXRegSpace < uRegionWidthInVoxels; uXRegSpace++) + { + // Note: In many cases the provided region will be (mostly) empty which means mesh vertices/indices + // are not generated and the only thing that is done for each cell is the computation of uCellIndex. + // It appears that retriving the voxel value is not so expensive and that it is the bitwise combining + // which actually carries the cost. + // + // If we really need to speed this up more then it may be possible to pack 4 8-bit cell indices into + // a single 32-bit value and then perform the bitwise logic on all four of them at the same time. + // However, this complicates the code and there would still be the cost of packing/unpacking so it's + // not clear if there is really a benefit. It's something to consider in the future. - std::swap(uNoOfNonEmptyCellsForSlice0, uNoOfNonEmptyCellsForSlice1); - pPreviousBitmask.swap(pCurrentBitmask); - m_pPreviousVertexIndicesX.swap(m_pCurrentVertexIndicesX); - m_pPreviousVertexIndicesY.swap(m_pCurrentVertexIndicesY); - m_pPreviousVertexIndicesZ.swap(m_pCurrentVertexIndicesZ); + // Each bit of the cell index specifies whether a given corner of the cell is above or below the threshold. + uint8_t uCellIndex = 0; - m_regSlicePrevious = m_regSliceCurrent; - m_regSliceCurrent.shift(Vector3DInt32(0,0,1)); - } + // Four bits of our cube index are obtained by looking at the cube index for + // the previous slice and copying four of those bits into their new positions. + uint8_t uPreviousCellIndexZ = pPreviousSliceCellIndices(uXRegSpace, uYRegSpace); + uPreviousCellIndexZ >>= 4; + uCellIndex |= uPreviousCellIndexZ; - m_meshCurrent->setOffset(m_regSizeInVoxels.getLowerCorner()); + // Two bits of our cube index are obtained by looking at the cube index for + // the previous row and copying two of those bits into their new positions. + uint8_t uPreviousCellIndexY = pPreviousRowCellIndices(uXRegSpace); + uPreviousCellIndexY &= 204; //204 = 128+64+8+4 + uPreviousCellIndexY >>= 2; + uCellIndex |= uPreviousCellIndexY; + + // One bit of our cube index are obtained by looking at the cube index for + // the previous cell and copying one of those bits into it's new position. + uint8_t UPreviousCellIndexX = uPreviousCellIndex; + UPreviousCellIndexX &= 170; //170 = 128+32+8+2 + UPreviousCellIndexX >>= 1; + uCellIndex |= UPreviousCellIndexX; + + // The last bit of our cube index is obtained by looking + // at the relevant voxel and comparing it to the threshold + typename VolumeType::VoxelType v111 = sampler.getVoxel(); + if (controller.convertToDensity(v111) < tThreshold) uCellIndex |= 128; + + // The current value becomes the previous value, ready for the next iteration. + uPreviousCellIndex = uCellIndex; + pPreviousRowCellIndices(uXRegSpace) = uCellIndex; + pPreviousSliceCellIndices(uXRegSpace, uYRegSpace) = uCellIndex; + + // 12 bits of uEdge determine whether a vertex is placed on each of the 12 edges of the cell. + uint16_t uEdge = edgeTable[uCellIndex]; + + // Test whether any vertices and indices should be generated for the current cell (i.e. it is occupied). + // Performance note: This condition is usually false because most cells in a volume are completely above + // or below the threshold and hence unoccupied. However, even when it is always false (testing on an empty + // volume) it still incurs significant overhead, probably because the code is large and bloats the for loop + // which contains it. On my empty volume test case the code as given runs in 34ms, but if I replace the + // condition with 'false' it runs in 24ms and gives the same output (i.e. none). + // + // An improvement is to move the code into a seperate function which does speed things up (30ms), but this + // is messy as the function needs to be passed about 10 differnt parameters, probably adding some overhead + // in its self. This does indeed seem to slow down the case when cells are occupied, by about 10-20%. + // + // Overall I don't know the right solution, but I'm leaving the code as-is to avoid making it messy. If we + // can reduce the number of parameters which need to be passed then it might be worth moving it into a + // function, or otherwise it may simply be worth trying to shorten the code (e.g. adding other function + // calls). For now we will leave it as-is, until we have more information from real-world profiling. + if (uEdge != 0) + { + auto v111Density = controller.convertToDensity(v111); + const Vector3DFloat n000 = computeCentralDifferenceGradient(sampler, controller); + + /* Find the vertices where the surface intersects the cube */ + if ((uEdge & 64) && (uXRegSpace > 0)) + { + sampler.moveNegativeX(); + typename VolumeType::VoxelType v011 = sampler.getVoxel(); + auto v011Density = controller.convertToDensity(v011); + const float fInterp = static_cast(tThreshold - v011Density) / static_cast(v111Density - v011Density); + + // Compute the position + const Vector3DFloat v3dPosition(static_cast(uXRegSpace - 1) + fInterp, static_cast(uYRegSpace), static_cast(uZRegSpace)); + + // Compute the normal + const Vector3DFloat n100 = computeCentralDifferenceGradient(sampler, controller); + Vector3DFloat v3dNormal = (n100*fInterp) + (n000*(1 - fInterp)); + + // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so + // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). + if (v3dNormal.lengthSquared() > 0.000001f) + { + v3dNormal.normalise(); + } + + // Allow the controller to decide how the material should be derived from the voxels. + const typename VolumeType::VoxelType uMaterial = controller.blendMaterials(v011, v111, fInterp); + + MarchingCubesVertex surfaceVertex; + const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); + surfaceVertex.encodedPosition = v3dScaledPosition; + surfaceVertex.encodedNormal = encodeNormal(v3dNormal); + surfaceVertex.data = uMaterial; + + const uint32_t uLastVertexIndex = result->addVertex(surfaceVertex); + pIndices(uXRegSpace, uYRegSpace).setX(uLastVertexIndex); + + sampler.movePositiveX(); + } + if ((uEdge & 32) && (uYRegSpace > 0)) + { + sampler.moveNegativeY(); + typename VolumeType::VoxelType v101 = sampler.getVoxel(); + auto v101Density = controller.convertToDensity(v101); + const float fInterp = static_cast(tThreshold - v101Density) / static_cast(v111Density - v101Density); + + // Compute the position + const Vector3DFloat v3dPosition(static_cast(uXRegSpace), static_cast(uYRegSpace - 1) + fInterp, static_cast(uZRegSpace)); + + // Compute the normal + const Vector3DFloat n010 = computeCentralDifferenceGradient(sampler, controller); + Vector3DFloat v3dNormal = (n010*fInterp) + (n000*(1 - fInterp)); + + // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so + // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). + if (v3dNormal.lengthSquared() > 0.000001f) + { + v3dNormal.normalise(); + } + + // Allow the controller to decide how the material should be derived from the voxels. + const typename VolumeType::VoxelType uMaterial = controller.blendMaterials(v101, v111, fInterp); + + MarchingCubesVertex surfaceVertex; + const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); + surfaceVertex.encodedPosition = v3dScaledPosition; + surfaceVertex.encodedNormal = encodeNormal(v3dNormal); + surfaceVertex.data = uMaterial; + + uint32_t uLastVertexIndex = result->addVertex(surfaceVertex); + pIndices(uXRegSpace, uYRegSpace).setY(uLastVertexIndex); + + sampler.movePositiveY(); + } + if ((uEdge & 1024) && (uZRegSpace > 0)) + { + sampler.moveNegativeZ(); + typename VolumeType::VoxelType v110 = sampler.getVoxel(); + auto v110Density = controller.convertToDensity(v110); + const float fInterp = static_cast(tThreshold - v110Density) / static_cast(v111Density - v110Density); + + // Compute the position + const Vector3DFloat v3dPosition(static_cast(uXRegSpace), static_cast(uYRegSpace), static_cast(uZRegSpace - 1) + fInterp); + + // Compute the normal + const Vector3DFloat n001 = computeCentralDifferenceGradient(sampler, controller); + Vector3DFloat v3dNormal = (n001*fInterp) + (n000*(1 - fInterp)); + + // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so + // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). + if (v3dNormal.lengthSquared() > 0.000001f) + { + v3dNormal.normalise(); + } + + // Allow the controller to decide how the material should be derived from the voxels. + const typename VolumeType::VoxelType uMaterial = controller.blendMaterials(v110, v111, fInterp); + + MarchingCubesVertex surfaceVertex; + const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); + surfaceVertex.encodedPosition = v3dScaledPosition; + surfaceVertex.encodedNormal = encodeNormal(v3dNormal); + surfaceVertex.data = uMaterial; + + const uint32_t uLastVertexIndex = result->addVertex(surfaceVertex); + pIndices(uXRegSpace, uYRegSpace).setZ(uLastVertexIndex); + + sampler.movePositiveZ(); + } + + // Now output the indices. For the first row, column or slice there aren't + // any (the region size in cells is one less than the region size in voxels) + if ((uXRegSpace != 0) && (uYRegSpace != 0) && (uZRegSpace != 0)) + { + + int32_t indlist[12]; + + /* Find the vertices where the surface intersects the cube */ + if (uEdge & 1) + { + indlist[0] = pPreviousIndices(uXRegSpace, uYRegSpace - 1).getX(); + } + if (uEdge & 2) + { + indlist[1] = pPreviousIndices(uXRegSpace, uYRegSpace).getY(); + } + if (uEdge & 4) + { + indlist[2] = pPreviousIndices(uXRegSpace, uYRegSpace).getX(); + } + if (uEdge & 8) + { + indlist[3] = pPreviousIndices(uXRegSpace - 1, uYRegSpace).getY(); + } + if (uEdge & 16) + { + indlist[4] = pIndices(uXRegSpace, uYRegSpace - 1).getX(); + } + if (uEdge & 32) + { + indlist[5] = pIndices(uXRegSpace, uYRegSpace).getY(); + } + if (uEdge & 64) + { + indlist[6] = pIndices(uXRegSpace, uYRegSpace).getX(); + } + if (uEdge & 128) + { + indlist[7] = pIndices(uXRegSpace - 1, uYRegSpace).getY(); + } + if (uEdge & 256) + { + indlist[8] = pIndices(uXRegSpace - 1, uYRegSpace - 1).getZ(); + } + if (uEdge & 512) + { + indlist[9] = pIndices(uXRegSpace, uYRegSpace - 1).getZ(); + } + if (uEdge & 1024) + { + indlist[10] = pIndices(uXRegSpace, uYRegSpace).getZ(); + } + if (uEdge & 2048) + { + indlist[11] = pIndices(uXRegSpace - 1, uYRegSpace).getZ(); + } + + for (int i = 0; triTable[uCellIndex][i] != -1; i += 3) + { + const int32_t ind0 = indlist[triTable[uCellIndex][i]]; + const int32_t ind1 = indlist[triTable[uCellIndex][i + 1]]; + const int32_t ind2 = indlist[triTable[uCellIndex][i + 2]]; + + if ((ind0 != -1) && (ind1 != -1) && (ind2 != -1)) + { + result->addTriangle(ind0, ind1, ind2); + } + } // For each triangle + } + } // For each cell + sampler.movePositiveX(); + } // For X + startOfRow.movePositiveY(); + } // For Y + startOfSlice.movePositiveZ(); + + pIndices.swap(pPreviousIndices); + } // For Z + + result->setOffset(region.getLowerCorner()); POLYVOX_LOG_TRACE("Marching cubes surface extraction took ", timer.elapsedTimeInMilliSeconds(), - "ms (Region size = ", m_regSizeInVoxels.getWidthInVoxels(), "x", m_regSizeInVoxels.getHeightInVoxels(), - "x", m_regSizeInVoxels.getDepthInVoxels(), ")"); - } - - template - template - uint32_t MarchingCubesSurfaceExtractor::computeBitmaskForSlice(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask) - { - m_uNoOfOccupiedCells = 0; - - const int32_t iMaxXVolSpace = m_regSliceCurrent.getUpperX(); - const int32_t iMaxYVolSpace = m_regSliceCurrent.getUpperY(); - - const int32_t iZVolSpace = m_regSliceCurrent.getLowerZ(); - - //Process the lower left corner - int32_t iYVolSpace = m_regSliceCurrent.getLowerY(); - int32_t iXVolSpace = m_regSliceCurrent.getLowerX(); - - uint32_t uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerX(); - uint32_t uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerY(); - - - m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace); - computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace); - - //Process the edge where x is minimal. - iXVolSpace = m_regSliceCurrent.getLowerX(); - m_sampVolume.setPosition(iXVolSpace, m_regSliceCurrent.getLowerY(), iZVolSpace); - for(iYVolSpace = m_regSliceCurrent.getLowerY() + 1; iYVolSpace <= iMaxYVolSpace; iYVolSpace++) - { - uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerX(); - uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerY(); - - m_sampVolume.movePositiveY(); - - computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace); - } - - //Process the edge where y is minimal. - iYVolSpace = m_regSliceCurrent.getLowerY(); - m_sampVolume.setPosition(m_regSliceCurrent.getLowerX(), iYVolSpace, iZVolSpace); - for(iXVolSpace = m_regSliceCurrent.getLowerX() + 1; iXVolSpace <= iMaxXVolSpace; iXVolSpace++) - { - uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerX(); - uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerY(); - - m_sampVolume.movePositiveX(); - - computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace); - } - - //Process all remaining elemnents of the slice. In this case, previous x and y values are always available - for(iYVolSpace = m_regSliceCurrent.getLowerY() + 1; iYVolSpace <= iMaxYVolSpace; iYVolSpace++) - { - m_sampVolume.setPosition(m_regSliceCurrent.getLowerX(), iYVolSpace, iZVolSpace); - for(iXVolSpace = m_regSliceCurrent.getLowerX() + 1; iXVolSpace <= iMaxXVolSpace; iXVolSpace++) - { - uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerX(); - uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerY(); - - m_sampVolume.movePositiveX(); - - computeBitmaskForCell(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace); - } - } - - return m_uNoOfOccupiedCells; - } - - template - template - void MarchingCubesSurfaceExtractor::computeBitmaskForCell(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask, uint32_t uXRegSpace, uint32_t uYRegSpace) - { - 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, uYRegSpace) = iCubeIndex; - - if(edgeTable[iCubeIndex] != 0) - { - ++m_uNoOfOccupiedCells; - } - } - - template - void MarchingCubesSurfaceExtractor::generateVerticesForSlice(const Array2DUint8& pCurrentBitmask, - Array2DInt32& m_pCurrentVertexIndicesX, - Array2DInt32& m_pCurrentVertexIndicesY, - Array2DInt32& m_pCurrentVertexIndicesZ) - { - const int32_t iZVolSpace = m_regSliceCurrent.getLowerZ(); - - //Iterate over each cell in the region - for(int32_t iYVolSpace = m_regSliceCurrent.getLowerY(); iYVolSpace <= m_regSliceCurrent.getUpperY(); iYVolSpace++) - { - const uint32_t uYRegSpace = iYVolSpace - m_regSizeInVoxels.getLowerY(); - - for(int32_t iXVolSpace = m_regSliceCurrent.getLowerX(); iXVolSpace <= m_regSliceCurrent.getUpperX(); iXVolSpace++) - { - //Current position - const uint32_t uXRegSpace = iXVolSpace - m_regSizeInVoxels.getLowerX(); - - //Determine the index into the edge table which tells us which vertices are inside of the surface - const 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(); - POLYVOX_ASSERT(v000 != v100, "Attempting to insert vertex between two voxels with the same value"); - const Vector3DFloat n100 = computeCentralDifferenceGradient(m_sampVolume); - - const 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.getLowerX()) + fInterp, static_cast(iYVolSpace - m_regSizeInVoxels.getLowerY()), static_cast(iZVolSpace - m_regSizeInCells.getLowerZ())); - const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); - - Vector3DFloat v3dNormal = (n100*fInterp) + (n000*(1-fInterp)); - - // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so - // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). - if(v3dNormal.lengthSquared() > 0.000001f) - { - v3dNormal.normalise(); - } - - // Allow the controller to decide how the material should be derived from the voxels. - const typename VolumeType::VoxelType uMaterial = m_controller.blendMaterials(v000, v100, fInterp); - - MarchingCubesVertex surfaceVertex; - surfaceVertex.encodedPosition = v3dScaledPosition; - surfaceVertex.encodedNormal = encodeNormal(v3dNormal); - surfaceVertex.data = uMaterial; - - const uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex); - m_pCurrentVertexIndicesX(iXVolSpace - m_regSizeInVoxels.getLowerX(), iYVolSpace - m_regSizeInVoxels.getLowerY()) = uLastVertexIndex; - - m_sampVolume.moveNegativeX(); - } - if (edgeTable[iCubeIndex] & 8) - { - m_sampVolume.movePositiveY(); - const typename VolumeType::VoxelType v010 = m_sampVolume.getVoxel(); - POLYVOX_ASSERT(v000 != v010, "Attempting to insert vertex between two voxels with the same value"); - const Vector3DFloat n010 = computeCentralDifferenceGradient(m_sampVolume); - - const 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.getLowerX()), static_cast(iYVolSpace - m_regSizeInVoxels.getLowerY()) + fInterp, static_cast(iZVolSpace - m_regSizeInVoxels.getLowerZ())); - const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); - - Vector3DFloat v3dNormal = (n010*fInterp) + (n000*(1-fInterp)); - - // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so - // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). - if(v3dNormal.lengthSquared() > 0.000001f) - { - v3dNormal.normalise(); - } - - // Allow the controller to decide how the material should be derived from the voxels. - const typename VolumeType::VoxelType uMaterial = m_controller.blendMaterials(v000, v010, fInterp); - - MarchingCubesVertex surfaceVertex; - surfaceVertex.encodedPosition = v3dScaledPosition; - surfaceVertex.encodedNormal = encodeNormal(v3dNormal); - surfaceVertex.data = uMaterial; - - uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex); - m_pCurrentVertexIndicesY(iXVolSpace - m_regSizeInVoxels.getLowerX(), iYVolSpace - m_regSizeInVoxels.getLowerY()) = uLastVertexIndex; - - m_sampVolume.moveNegativeY(); - } - if (edgeTable[iCubeIndex] & 256) - { - m_sampVolume.movePositiveZ(); - const typename VolumeType::VoxelType v001 = m_sampVolume.getVoxel(); - POLYVOX_ASSERT(v000 != v001, "Attempting to insert vertex between two voxels with the same value"); - const Vector3DFloat n001 = computeCentralDifferenceGradient(m_sampVolume); - - const 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.getLowerX()), static_cast(iYVolSpace - m_regSizeInVoxels.getLowerY()), static_cast(iZVolSpace - m_regSizeInVoxels.getLowerZ()) + fInterp); - const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); - - Vector3DFloat v3dNormal = (n001*fInterp) + (n000*(1-fInterp)); - // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so - // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). - if(v3dNormal.lengthSquared() > 0.000001f) - { - v3dNormal.normalise(); - } - - // Allow the controller to decide how the material should be derived from the voxels. - const typename VolumeType::VoxelType uMaterial = m_controller.blendMaterials(v000, v001, fInterp); - - MarchingCubesVertex surfaceVertex; - surfaceVertex.encodedPosition = v3dScaledPosition; - surfaceVertex.encodedNormal = encodeNormal(v3dNormal); - surfaceVertex.data = uMaterial; - - const uint32_t uLastVertexIndex = m_meshCurrent->addVertex(surfaceVertex); - m_pCurrentVertexIndicesZ(iXVolSpace - m_regSizeInVoxels.getLowerX(), iYVolSpace - m_regSizeInVoxels.getLowerY()) = 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; - } - - const int32_t iZVolSpace = m_regSlicePrevious.getLowerZ(); - - for(int32_t iYVolSpace = m_regSlicePrevious.getLowerY(); iYVolSpace <= m_regSizeInCells.getUpperY(); iYVolSpace++) - { - for(int32_t iXVolSpace = m_regSlicePrevious.getLowerX(); iXVolSpace <= m_regSizeInCells.getUpperX(); iXVolSpace++) - { - m_sampVolume.setPosition(iXVolSpace,iYVolSpace,iZVolSpace); - - //Current position - const uint32_t uXRegSpace = m_sampVolume.getPosition().getX() - m_regSizeInVoxels.getLowerX(); - const uint32_t uYRegSpace = m_sampVolume.getPosition().getY() - m_regSizeInVoxels.getLowerY(); - - //Determine the index into the edge table which tells us which vertices are inside of the surface - const 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); - } - if (edgeTable[iCubeIndex] & 2) - { - indlist[1] = m_pPreviousVertexIndicesY(uXRegSpace + 1, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 4) - { - indlist[2] = m_pPreviousVertexIndicesX(uXRegSpace, uYRegSpace + 1); - } - if (edgeTable[iCubeIndex] & 8) - { - indlist[3] = m_pPreviousVertexIndicesY(uXRegSpace, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 16) - { - indlist[4] = m_pCurrentVertexIndicesX(uXRegSpace, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 32) - { - indlist[5] = m_pCurrentVertexIndicesY(uXRegSpace + 1, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 64) - { - indlist[6] = m_pCurrentVertexIndicesX(uXRegSpace, uYRegSpace + 1); - } - if (edgeTable[iCubeIndex] & 128) - { - indlist[7] = m_pCurrentVertexIndicesY(uXRegSpace, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 256) - { - indlist[8] = m_pPreviousVertexIndicesZ(uXRegSpace, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 512) - { - indlist[9] = m_pPreviousVertexIndicesZ(uXRegSpace + 1, uYRegSpace); - } - if (edgeTable[iCubeIndex] & 1024) - { - indlist[10] = m_pPreviousVertexIndicesZ(uXRegSpace + 1, uYRegSpace + 1); - } - if (edgeTable[iCubeIndex] & 2048) - { - indlist[11] = m_pPreviousVertexIndicesZ(uXRegSpace, uYRegSpace + 1); - } - - for (int i=0;triTable[iCubeIndex][i]!=-1;i+=3) - { - const int32_t ind0 = indlist[triTable[iCubeIndex][i ]]; - const int32_t ind1 = indlist[triTable[iCubeIndex][i+1]]; - const 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 - } - } + "ms (Region size = ", region.getWidthInVoxels(), "x", region.getHeightInVoxels(), + "x", region.getDepthInVoxels(), ")"); + } } diff --git a/include/PolyVox/PolyVoxForwardDeclarations.h b/include/PolyVox/PolyVoxForwardDeclarations.h index 9dbc44ca..c0c1f872 100644 --- a/include/PolyVox/PolyVoxForwardDeclarations.h +++ b/include/PolyVox/PolyVoxForwardDeclarations.h @@ -98,11 +98,6 @@ namespace PolyVox //////////////////////////////////////////////////////////////////////////////// template class FilePager; - //////////////////////////////////////////////////////////////////////////////// - // MarchingCubesSurfaceExtractor - //////////////////////////////////////////////////////////////////////////////// - template class MarchingCubesSurfaceExtractor; - //////////////////////////////////////////////////////////////////////////////// // MarchingCubesVertex //////////////////////////////////////////////////////////////////////////////// diff --git a/tests/TestSurfaceExtractor.cpp b/tests/TestSurfaceExtractor.cpp index 3850b9cb..dd800d94 100644 --- a/tests/TestSurfaceExtractor.cpp +++ b/tests/TestSurfaceExtractor.cpp @@ -168,18 +168,18 @@ void TestSurfaceExtractor::testBehaviour() // This basic test just uses the default controller and automatically generates a mesh of the appropriate type. auto uintVol = createAndFillVolume< RawVolume >(); auto uintMesh = extractMarchingCubesMesh(uintVol, uintVol->getEnclosingRegion()); - QCOMPARE(uintMesh.getNoOfVertices(), uint32_t(12096)); // Verifies size of mesh and that we have 32-bit indices + QCOMPARE(uintMesh.getNoOfVertices(), uint32_t(6048)); // Verifies size of mesh and that we have 32-bit indices QCOMPARE(uintMesh.getNoOfIndices(), uint32_t(35157)); // Verifies size of mesh - QCOMPARE(uintMesh.getIndex(100), uint32_t(44)); // Verifies that we have 32-bit indices + QCOMPARE(uintMesh.getIndex(100), uint32_t(24)); // Verifies that we have 32-bit indices QCOMPARE(uintMesh.getVertex(100).data, uint8_t(1)); // Not really meaningful for a primative type // This test makes use of a custom controller auto floatVol = createAndFillVolume< RawVolume >(); CustomMarchingCubesController floatCustomController; auto floatMesh = extractMarchingCubesMesh(floatVol, floatVol->getEnclosingRegion(), floatCustomController); - QCOMPARE(floatMesh.getNoOfVertices(), uint32_t(16113)); // Verifies size of mesh and that we have 32-bit indices + QCOMPARE(floatMesh.getNoOfVertices(), uint32_t(3825)); // Verifies size of mesh and that we have 32-bit indices QCOMPARE(floatMesh.getNoOfIndices(), uint32_t(22053)); // Verifies size of mesh - QCOMPARE(floatMesh.getIndex(100), uint32_t(26)); // Verifies that we have 32-bit indices + QCOMPARE(floatMesh.getIndex(100), uint32_t(119)); // Verifies that we have 32-bit indices QCOMPARE(floatMesh.getVertex(100).data, float(1.0f)); // Not really meaningful for a primative type // This test makes use of a user provided mesh. It uses the default controller, but we have to explicitly provide this because C++ won't let us @@ -187,7 +187,7 @@ void TestSurfaceExtractor::testBehaviour() auto intVol = createAndFillVolume< RawVolume >(); Mesh< MarchingCubesVertex< int8_t >, uint16_t > intMesh; extractMarchingCubesMeshCustom(intVol, intVol->getEnclosingRegion(), &intMesh); - QCOMPARE(intMesh.getNoOfVertices(), uint16_t(11718)); // Verifies size of mesh and that we have 16-bit indices + QCOMPARE(intMesh.getNoOfVertices(), uint16_t(5859)); // Verifies size of mesh and that we have 16-bit indices QCOMPARE(intMesh.getNoOfIndices(), uint32_t(34041)); // Verifies size of mesh QCOMPARE(intMesh.getIndex(100), uint16_t(29)); // Verifies that we have 16-bit indices QCOMPARE(intMesh.getVertex(100).data, int8_t(1)); // Not really meaningful for a primative type @@ -197,17 +197,17 @@ void TestSurfaceExtractor::testBehaviour() CustomMarchingCubesController doubleCustomController; Mesh< MarchingCubesVertex< double >, uint16_t > doubleMesh; extractMarchingCubesMeshCustom(doubleVol, doubleVol->getEnclosingRegion(), &doubleMesh, doubleCustomController); - QCOMPARE(doubleMesh.getNoOfVertices(), uint16_t(16113)); // Verifies size of mesh and that we have 32-bit indices + QCOMPARE(doubleMesh.getNoOfVertices(), uint16_t(3825)); // Verifies size of mesh and that we have 32-bit indices QCOMPARE(doubleMesh.getNoOfIndices(), uint32_t(22053)); // Verifies size of mesh - QCOMPARE(doubleMesh.getIndex(100), uint16_t(26)); // Verifies that we have 32-bit indices + QCOMPARE(doubleMesh.getIndex(100), uint16_t(119)); // Verifies that we have 32-bit indices QCOMPARE(doubleMesh.getVertex(100).data, double(1.0f)); // Not really meaningful for a primative type // This test ensures the extractor works on a non-primitive voxel type. auto materialVol = createAndFillVolume< RawVolume >(); auto materialMesh = extractMarchingCubesMesh(materialVol, materialVol->getEnclosingRegion()); - QCOMPARE(materialMesh.getNoOfVertices(), uint32_t(12096)); // Verifies size of mesh and that we have 32-bit indices + QCOMPARE(materialMesh.getNoOfVertices(), uint32_t(6048)); // Verifies size of mesh and that we have 32-bit indices QCOMPARE(materialMesh.getNoOfIndices(), uint32_t(35157)); // Verifies size of mesh - QCOMPARE(materialMesh.getIndex(100), uint32_t(44)); // Verifies that we have 32-bit indices + QCOMPARE(materialMesh.getIndex(100), uint32_t(24)); // Verifies that we have 32-bit indices QCOMPARE(materialMesh.getVertex(100).data.getMaterial(), uint16_t(79)); // Verify the data attached to the vertex } @@ -224,7 +224,7 @@ void TestSurfaceExtractor::testNoiseVolumePerformance() auto noiseVol = createAndFillVolumeWithNoise< PagedVolume >(128, 128, -1.0f, 1.0f); Mesh< MarchingCubesVertex< float >, uint16_t > noiseMesh; QBENCHMARK{ extractMarchingCubesMeshCustom(noiseVol, Region(32, 32, 32, 63, 63, 63), &noiseMesh); } - QCOMPARE(noiseMesh.getNoOfVertices(), uint16_t(36755)); + QCOMPARE(noiseMesh.getNoOfVertices(), uint16_t(35672)); } QTEST_MAIN(TestSurfaceExtractor)