AnotherFoxGuy/polyvox
1
0
Fork 0
Code Issues Pull Requests Activity

Merge branch 'feature/extractor-optimizations' into develop

Browse Source
This commit is contained in:
David Williams 2015-05-29 22:39:57 +02:00
commit 5aa631da8f
4 changed files with 321 additions and 780 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

180
include/PolyVox/MarchingCubesSurfaceExtractor.h
View File

@ -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<bool isPrevZAvail>
uint32_t computeBitmaskForSlice(const Array2DUint8& pPreviousBitmask, Array2DUint8& pCurrentBitmask);
//Compute the cell bitmask for a given cell.
template<bool isPrevXAvail, bool isPrevYAvail, bool isPrevZAvail>
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<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py0pz()));
float voxel1px = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py0pz()));
float voxel1ny = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny0pz()));
float voxel1py = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py0pz()));
float voxel1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1nz()));
float voxel1pz = static_cast<float>(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<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1ny1nz()));
const float pVoxel1nx1ny0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1ny0pz()));
const float pVoxel1nx1ny1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1ny1pz()));
const float pVoxel1nx0py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py1nz()));
const float pVoxel1nx0py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py0pz()));
const float pVoxel1nx0py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx0py1pz()));
const float pVoxel1nx1py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1py1nz()));
const float pVoxel1nx1py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1py0pz()));
const float pVoxel1nx1py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1nx1py1pz()));
const float pVoxel0px1ny1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny1nz()));
const float pVoxel0px1ny0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny0pz()));
const float pVoxel0px1ny1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1ny1pz()));
const float pVoxel0px0py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1nz()));
//const float pVoxel0px0py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py0pz()));
const float pVoxel0px0py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px0py1pz()));
const float pVoxel0px1py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py1nz()));
const float pVoxel0px1py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py0pz()));
const float pVoxel0px1py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel0px1py1pz()));
const float pVoxel1px1ny1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1ny1nz()));
const float pVoxel1px1ny0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1ny0pz()));
const float pVoxel1px1ny1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1ny1pz()));
const float pVoxel1px0py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py1nz()));
const float pVoxel1px0py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py0pz()));
const float pVoxel1px0py1pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px0py1pz()));
const float pVoxel1px1py1nz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1py1nz()));
const float pVoxel1px1py0pz = static_cast<float>(m_controller.convertToDensity(volIter.peekVoxel1px1py0pz()));
const float pVoxel1px1py1pz = static_cast<float>(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<typename VolumeType::VoxelType> >
void extractMarchingCubesMeshCustom(VolumeType* volData, Region region, MeshType* result, ControllerType controller = ControllerType())
{
MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType> extractor(volData, region, result, controller);
extractor.execute();
}
void extractMarchingCubesMeshCustom(VolumeType* volData, Region region, MeshType* result, ControllerType controller = ControllerType());
template< typename VolumeType, typename ControllerType = DefaultMarchingCubesController<typename VolumeType::VoxelType> >
Mesh<MarchingCubesVertex<typename VolumeType::VoxelType> > extractMarchingCubesMesh(VolumeType* volData, Region region, ControllerType controller = ControllerType())

896
include/PolyVox/MarchingCubesSurfaceExtractor.inl
View File

@ -25,611 +25,333 @@ freely, subject to the following restrictions:
namespace PolyVox
{
template<typename VolumeType, typename MeshType, typename ControllerType>
MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType>::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<float>(controller.convertToDensity(volIter.peekVoxel1nx0py0pz()));
float voxel1px = static_cast<float>(controller.convertToDensity(volIter.peekVoxel1px0py0pz()));
float voxel1ny = static_cast<float>(controller.convertToDensity(volIter.peekVoxel0px1ny0pz()));
float voxel1py = static_cast<float>(controller.convertToDensity(volIter.peekVoxel0px1py0pz()));
float voxel1nz = static_cast<float>(controller.convertToDensity(volIter.peekVoxel0px0py1nz()));
float voxel1pz = static_cast<float>(controller.convertToDensity(volIter.peekVoxel0px0py1pz()));
return Vector3DFloat
(
voxel1nx - voxel1px,
voxel1ny - voxel1py,
voxel1nz - voxel1pz
);
}
template<typename VolumeType, typename MeshType, typename ControllerType>
void MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType>::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<false>(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<true>(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<float>(tThreshold - v011Density) / static_cast<float>(v111Density - v011Density);
// Compute the position
const Vector3DFloat v3dPosition(static_cast<float>(uXRegSpace - 1) + fInterp, static_cast<float>(uYRegSpace), static_cast<float>(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<typename VolumeType::VoxelType> surfaceVertex;
const Vector3DUint16 v3dScaledPosition(static_cast<uint16_t>(v3dPosition.getX() * 256.0f), static_cast<uint16_t>(v3dPosition.getY() * 256.0f), static_cast<uint16_t>(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<float>(tThreshold - v101Density) / static_cast<float>(v111Density - v101Density);
// Compute the position
const Vector3DFloat v3dPosition(static_cast<float>(uXRegSpace), static_cast<float>(uYRegSpace - 1) + fInterp, static_cast<float>(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<typename VolumeType::VoxelType> surfaceVertex;
const Vector3DUint16 v3dScaledPosition(static_cast<uint16_t>(v3dPosition.getX() * 256.0f), static_cast<uint16_t>(v3dPosition.getY() * 256.0f), static_cast<uint16_t>(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<float>(tThreshold - v110Density) / static_cast<float>(v111Density - v110Density);
// Compute the position
const Vector3DFloat v3dPosition(static_cast<float>(uXRegSpace), static_cast<float>(uYRegSpace), static_cast<float>(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<typename VolumeType::VoxelType> surfaceVertex;
const Vector3DUint16 v3dScaledPosition(static_cast<uint16_t>(v3dPosition.getX() * 256.0f), static_cast<uint16_t>(v3dPosition.getY() * 256.0f), static_cast<uint16_t>(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<typename VolumeType, typename MeshType, typename ControllerType>
template<bool isPrevZAvail>
uint32_t MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType>::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<false, false, isPrevZAvail>(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<false, true, isPrevZAvail>(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<true, false, isPrevZAvail>(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<true, true, isPrevZAvail>(pPreviousBitmask, pCurrentBitmask, uXRegSpace, uYRegSpace);
}
}
return m_uNoOfOccupiedCells;
}
template<typename VolumeType, typename MeshType, typename ControllerType>
template<bool isPrevXAvail, bool isPrevYAvail, bool isPrevZAvail>
void MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType>::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<typename VolumeType, typename MeshType, typename ControllerType>
void MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType>::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<float>(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast<float>(m_controller.convertToDensity(v100) - m_controller.convertToDensity(v000));
const Vector3DFloat v3dPosition(static_cast<float>(iXVolSpace - m_regSizeInVoxels.getLowerX()) + fInterp, static_cast<float>(iYVolSpace - m_regSizeInVoxels.getLowerY()), static_cast<float>(iZVolSpace - m_regSizeInCells.getLowerZ()));
const Vector3DUint16 v3dScaledPosition(static_cast<uint16_t>(v3dPosition.getX() * 256.0f), static_cast<uint16_t>(v3dPosition.getY() * 256.0f), static_cast<uint16_t>(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<typename VolumeType::VoxelType> 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<float>(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast<float>(m_controller.convertToDensity(v010) - m_controller.convertToDensity(v000));
const Vector3DFloat v3dPosition(static_cast<float>(iXVolSpace - m_regSizeInVoxels.getLowerX()), static_cast<float>(iYVolSpace - m_regSizeInVoxels.getLowerY()) + fInterp, static_cast<float>(iZVolSpace - m_regSizeInVoxels.getLowerZ()));
const Vector3DUint16 v3dScaledPosition(static_cast<uint16_t>(v3dPosition.getX() * 256.0f), static_cast<uint16_t>(v3dPosition.getY() * 256.0f), static_cast<uint16_t>(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<typename VolumeType::VoxelType> 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<float>(m_tThreshold - m_controller.convertToDensity(v000)) / static_cast<float>(m_controller.convertToDensity(v001) - m_controller.convertToDensity(v000));
const Vector3DFloat v3dPosition(static_cast<float>(iXVolSpace - m_regSizeInVoxels.getLowerX()), static_cast<float>(iYVolSpace - m_regSizeInVoxels.getLowerY()), static_cast<float>(iZVolSpace - m_regSizeInVoxels.getLowerZ()) + fInterp);
const Vector3DUint16 v3dScaledPosition(static_cast<uint16_t>(v3dPosition.getX() * 256.0f), static_cast<uint16_t>(v3dPosition.getY() * 256.0f), static_cast<uint16_t>(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<typename VolumeType::VoxelType> 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<typename VolumeType, typename MeshType, typename ControllerType>
void MarchingCubesSurfaceExtractor<VolumeType, MeshType, ControllerType>::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(), ")");
}
}

5
include/PolyVox/PolyVoxForwardDeclarations.h
View File

@ -98,11 +98,6 @@ namespace PolyVox
////////////////////////////////////////////////////////////////////////////////
template <typename VoxelType> class FilePager;
////////////////////////////////////////////////////////////////////////////////
// MarchingCubesSurfaceExtractor
////////////////////////////////////////////////////////////////////////////////
template<typename VolumeType, typename MeshType, typename ControllerType> class MarchingCubesSurfaceExtractor;
////////////////////////////////////////////////////////////////////////////////
// MarchingCubesVertex
////////////////////////////////////////////////////////////////////////////////

20
tests/TestSurfaceExtractor.cpp
View File

@ -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<uint8_t> >();
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<float> >();
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<int8_t> >();
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<MaterialDensityPair88> >();
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<float> >(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)