Improved performance of Dual Contouring

Performance is up from 10 times slower than MC to only 3 times.
About a third of that time is spent calculating gradients.
This commit is contained in:
Matt Williams 2014-01-09 20:07:15 +00:00
parent 7877600538
commit 601b2a6d21
3 changed files with 148 additions and 172 deletions

View File

@ -58,6 +58,7 @@ SET(CORE_INC_FILES
include/PolyVoxCore/DefaultIsQuadNeeded.h
include/PolyVoxCore/DefaultMarchingCubesController.h
include/PolyVoxCore/Density.h
include/PolyVoxCore/DualContouringSurfaceExtractor.h
include/PolyVoxCore/FilePager.h
include/PolyVoxCore/GradientEstimators.h
include/PolyVoxCore/GradientEstimators.inl

View File

@ -20,67 +20,42 @@ freely, subject to the following restrictions:
3. This notice may not be removed or altered from any source
distribution.
*******************************************************************************/
#include "SimpleVolume.h"
#include "PolyVoxCore/SurfaceMesh.h"
#include "PolyVoxCore/VertexTypes.h"
#include "Impl/Timer.h"
#include "Impl/QEF.h"
//BUG We will get duplucation of edges if the surface is along region boundaries
namespace PolyVox
{
namespace
{
struct CellData
{
uint32_t vertexIndex;
bool x;
bool y;
bool z;
bool hasVertex;
};
template<typename VoxelType>
struct EdgeData
{
VoxelType intersectionDistance;
VoxelType edgeLength;
bool intersects;
Vector3DFloat normal;
float fraction; ///<fraction (0.0-1.0) along the edge in the positive direction that the intersection happens
bool intersects;
};
template<typename SamplerType>
Vector3DFloat computeCentralDifferenceGradient(const SamplerType& volIter)
template<typename VoxelType>
struct CellData
{
//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>(volIter.peekVoxel1nx0py0pz());
float voxel1px = static_cast<float>(volIter.peekVoxel1px0py0pz());
float voxel1ny = static_cast<float>(volIter.peekVoxel0px1ny0pz());
float voxel1py = static_cast<float>(volIter.peekVoxel0px1py0pz());
float voxel1nz = static_cast<float>(volIter.peekVoxel0px0py1nz());
float voxel1pz = static_cast<float>(volIter.peekVoxel0px0py1pz());
return Vector3DFloat
(
voxel1nx - voxel1px,
voxel1ny - voxel1py,
voxel1nz - voxel1pz
);
}
EdgeData<VoxelType> edges[3];
uint32_t vertexIndex;
};
template<typename VoxelType, typename ThresholdType>
EdgeData<VoxelType> calculateEdge(VoxelType vA, VoxelType vB, Vector3DFloat gA, Vector3DFloat gB, ThresholdType threshold)
EdgeData<VoxelType> calculateEdge(const VoxelType& vA, const VoxelType& vB, const Vector3DFloat& gA, const Vector3DFloat& gB, const ThresholdType& threshold)
{
EdgeData<VoxelType> edge;
edge.intersectionDistance = vA - threshold;
edge.edgeLength = vA - vB;
auto fraction = static_cast<float>(edge.intersectionDistance) / static_cast<float>(edge.edgeLength);
if(fraction > 1.0 || fraction < 0.0)
edge.fraction = static_cast<float>(vA - threshold) / static_cast<float>(vA - vB);
if(edge.fraction > 1.0 || edge.fraction < 0.0)
{
edge.intersects = false;
return edge;
@ -90,7 +65,7 @@ namespace PolyVox
edge.intersects = true;
}
edge.normal = ((gA * static_cast<float>(edge.edgeLength - edge.intersectionDistance) + gB * static_cast<float>(edge.intersectionDistance)) / static_cast<float>(2 * edge.edgeLength));
edge.normal = (gA * edge.fraction + gB * (1.0f-edge.fraction));
if(edge.normal.lengthSquared() > 0.000001f)
{
edge.normal.normalise();
@ -100,29 +75,29 @@ namespace PolyVox
}
template<typename VoxelType>
PositionMaterialNormal computeVertex(std::array<EdgeData<VoxelType>, 12> edges)
PositionMaterialNormal computeVertex(EdgeData<VoxelType>* edges[12])
{
Vector3DFloat massPoint{0,0,0}; //The average of the intersection vertices
std::array<Vector3DFloat, 12> vertices;
Vector3DFloat vertices[12];
vertices[0] = {static_cast<float>(edges[0].intersectionDistance) / static_cast<float>(edges[0].edgeLength), 0, 0};
vertices[1] = {0, static_cast<float>(edges[1].intersectionDistance) / static_cast<float>(edges[1].edgeLength), 0};
vertices[2] = {0, 0, static_cast<float>(edges[2].intersectionDistance) / static_cast<float>(edges[2].edgeLength)};
vertices[3] = {1, static_cast<float>(edges[3].intersectionDistance) / static_cast<float>(edges[3].edgeLength), 0};
vertices[4] = {1, 0, static_cast<float>(edges[4].intersectionDistance) / static_cast<float>(edges[4].edgeLength)};
vertices[5] = {0, 1, static_cast<float>(edges[5].intersectionDistance) / static_cast<float>(edges[5].edgeLength)};
vertices[6] = {static_cast<float>(edges[6].intersectionDistance) / static_cast<float>(edges[6].edgeLength), 1, 0};
vertices[7] = {static_cast<float>(edges[7].intersectionDistance) / static_cast<float>(edges[7].edgeLength), 0, 1};
vertices[8] = {0, static_cast<float>(edges[8].intersectionDistance) / static_cast<float>(edges[8].edgeLength), 1};
vertices[9] = {1, 1, static_cast<float>(edges[9].intersectionDistance) / static_cast<float>(edges[9].edgeLength)};
vertices[10] = {1, static_cast<float>(edges[10].intersectionDistance) / static_cast<float>(edges[10].edgeLength), 1};
vertices[11] = {static_cast<float>(edges[11].intersectionDistance) / static_cast<float>(edges[11].edgeLength), 1, 1};
vertices[0] = {edges[0]->fraction, 0, 0};
vertices[1] = {0, edges[1]->fraction, 0};
vertices[2] = {0, 0, edges[2]->fraction};
vertices[3] = {1, edges[3]->fraction, 0};
vertices[4] = {1, 0, edges[4]->fraction};
vertices[5] = {0, 1, edges[5]->fraction};
vertices[6] = {edges[6]->fraction, 1, 0};
vertices[7] = {edges[7]->fraction, 0, 1};
vertices[8] = {0, edges[8]->fraction, 1};
vertices[9] = {1, 1, edges[9]->fraction};
vertices[10] = {1, edges[10]->fraction, 1};
vertices[11] = {edges[11]->fraction, 1, 1};
int numIntersections = 0;
for(int i = 0; i < 12; ++i)
{
if(!edges[i].intersects)
if(!edges[i]->intersects)
{
continue;
}
@ -141,12 +116,12 @@ namespace PolyVox
for(int i = 0; i < 12; ++i)
{
if(!edges[i].intersects)
if(!edges[i]->intersects)
{
continue;
}
Vector3DFloat normal = edges[i].normal;
Vector3DFloat normal = edges[i]->normal;
matrix[rows][0] = normal.getX();
matrix[rows][1] = normal.getY();
matrix[rows][2] = normal.getZ();
@ -171,6 +146,11 @@ namespace PolyVox
return {vertexPosition, cellVertexNormal, 0};
}
uint32_t convert(uint32_t x, uint32_t y, uint32_t z, uint32_t X, uint32_t Y)
{
return z*Y*X+y*X+x;
}
}
template<typename VolumeType>
@ -178,8 +158,13 @@ namespace PolyVox
{
Timer timer;
Region cellRegion{Vector3DInt32{0,0,0}, region.getDimensionsInVoxels()};
SimpleVolume<CellData> cells{cellRegion}; //Cells should be 1 bigger than 'region' in each dimension
Region cellRegion{Vector3DInt32{0,0,0}, region.getDimensionsInVoxels() + Vector3DInt32{1,1,1}};
auto cellRegionXDimension = cellRegion.getDimensionsInVoxels().getX();
auto cellRegionYDimension = cellRegion.getDimensionsInVoxels().getY();
auto cellRegionZDimension = cellRegion.getDimensionsInVoxels().getZ();
std::vector<CellData<typename VolumeType::VoxelType>> cells;
cells.reserve(cellRegionXDimension * cellRegionYDimension * cellRegionZDimension);
typename VolumeType::Sampler volSampler{volData};
volSampler.setPosition(region.getLowerCorner());
@ -189,139 +174,134 @@ namespace PolyVox
SurfaceMesh<PositionMaterialNormal> mesh;
std::array<EdgeData<typename VolumeType::VoxelType>, 12> edges; //Create this now but it will be overwritten for each cell
EdgeData<typename VolumeType::VoxelType>* edges[12]; //Create this now but it will be overwritten for each cell
for(int32_t cellX = 0; cellX < cellRegion.getDimensionsInVoxels().getX(); cellX++)
const auto lowerCornerX = region.getLowerCorner().getX();
const auto lowerCornerY = region.getLowerCorner().getZ();
const auto lowerCornerZ = region.getLowerCorner().getX();
for(int32_t cellZ = 0; cellZ < cellRegionZDimension; cellZ++)
{
for(int32_t cellY = 0; cellY < cellRegion.getDimensionsInVoxels().getY(); cellY++)
for(int32_t cellY = 0; cellY < cellRegionYDimension; cellY++)
{
for(int32_t cellZ = 0; cellZ < cellRegion.getDimensionsInVoxels().getZ(); cellZ++)
for(int32_t cellX = 0; cellX < cellRegionXDimension; cellX++)
{
//For each cell, calculate the vertex position
volSampler.setPosition(region.getLowerCorner() + Vector3DInt32{cellX, cellY, cellZ} - Vector3DInt32{1,1,1});
typename VolumeType::Sampler gradientSampler{volData};
volSampler.setPosition(lowerCornerX+cellX-1, lowerCornerY+cellY-1, lowerCornerZ+cellZ-1);
typename VolumeType::VoxelType v000{volSampler.getVoxel()};
typename VolumeType::VoxelType v100{volSampler.peekVoxel1px0py0pz()};
typename VolumeType::VoxelType v010{volSampler.peekVoxel0px1py0pz()};
typename VolumeType::VoxelType v110{volSampler.peekVoxel1px1py0pz()};
typename VolumeType::VoxelType v001{volSampler.peekVoxel0px0py1pz()};
typename VolumeType::VoxelType v101{volSampler.peekVoxel1px0py1pz()};
typename VolumeType::VoxelType v011{volSampler.peekVoxel0px1py1pz()};
typename VolumeType::VoxelType v111{volSampler.peekVoxel1px1py1pz()};
gradientSampler.setPosition(volSampler.getPosition());
Vector3DFloat g000 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{1,0,0});
Vector3DFloat g100 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{0,1,0});
Vector3DFloat g010 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{1,1,0});
Vector3DFloat g110 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{0,0,1});
Vector3DFloat g001 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{1,0,1});
Vector3DFloat g101 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{0,1,1});
Vector3DFloat g011 = computeCentralDifferenceGradient(gradientSampler);
gradientSampler.setPosition(volSampler.getPosition() + Vector3DInt32{1,1,1});
Vector3DFloat g111 = computeCentralDifferenceGradient(gradientSampler);
auto voxel1nx = static_cast<float>(volSampler.peekVoxel1nx0py0pz());
//auto voxel1px = static_cast<float>(volSampler.peekVoxel1px0py0pz());
auto voxel1ny = static_cast<float>(volSampler.peekVoxel0px1ny0pz());
//auto voxel1py = static_cast<float>(volSampler.peekVoxel0px1py0pz());
auto voxel1nz = static_cast<float>(volSampler.peekVoxel0px0py1nz());
//auto voxel1pz = static_cast<float>(volSampler.peekVoxel0px0py1pz());
Vector3DFloat g000(voxel1nx - v100, voxel1ny - v010, voxel1nz - v001);
volSampler.movePositiveX();
auto voxel2px = static_cast<float>(volSampler.peekVoxel1px0py0pz());
auto voxel1px1ny = static_cast<float>(volSampler.peekVoxel0px1ny0pz());
auto voxel1px1py = static_cast<float>(volSampler.peekVoxel0px1py0pz());
auto voxel1px1nz = static_cast<float>(volSampler.peekVoxel0px0py1nz());
auto voxel1px1pz = static_cast<float>(volSampler.peekVoxel0px0py1pz());
Vector3DFloat g100(v000 - voxel2px, voxel1px1ny - voxel1px1py, voxel1px1nz - voxel1px1pz);
volSampler.moveNegativeX();
volSampler.movePositiveY();
auto voxel1nx1py = static_cast<float>(volSampler.peekVoxel1nx0py0pz());
auto voxel2py = static_cast<float>(volSampler.peekVoxel0px1py0pz());
auto voxel1py1nz = static_cast<float>(volSampler.peekVoxel0px0py1nz());
auto voxel1py1pz = static_cast<float>(volSampler.peekVoxel0px0py1pz());
Vector3DFloat g010(voxel1nx1py - voxel1px1py, v000 - voxel2py, voxel1py1nz - voxel1py1pz);
volSampler.moveNegativeY();
volSampler.movePositiveZ();
auto voxel1nx1pz = static_cast<float>(volSampler.peekVoxel1nx0py0pz());
auto voxel1ny1pz = static_cast<float>(volSampler.peekVoxel0px1ny0pz());
auto voxel2pz = static_cast<float>(volSampler.peekVoxel0px0py1pz());
Vector3DFloat g001(voxel1nx1pz - voxel1px1pz, voxel1ny1pz - voxel1py1pz, v000 - voxel2pz);
edges[0] = calculateEdge(v000, v100, g000, g100, threshold);
edges[1] = calculateEdge(v000, v010, g000, g010, threshold);
edges[2] = calculateEdge(v000, v001, g000, g001, threshold);
cells.push_back({calculateEdge(v000, v100, g000, g100, threshold), calculateEdge(v000, v010, g000, g010, threshold), calculateEdge(v000, v001, g000, g001, threshold)});
edges[3] = calculateEdge(v100, v110, g100, g110, threshold);
edges[4] = calculateEdge(v100, v101, g100, g101, threshold);
edges[5] = calculateEdge(v010, v011, g010, g011, threshold);
edges[6] = calculateEdge(v010, v110, g010, g110, threshold);
edges[7] = calculateEdge(v001, v101, g001, g101, threshold);
edges[8] = calculateEdge(v001, v011, g001, g011, threshold);
edges[9] = calculateEdge(v110, v111, g110, g111, threshold);
edges[10] = calculateEdge(v101, v111, g101, g111, threshold);
edges[11] = calculateEdge(v011, v111, g011, g111, threshold);
CellData cell;
cell.x = edges[0].intersects;
cell.y = edges[1].intersects;
cell.z = edges[2].intersects;
if(edges[0].intersects || edges[1].intersects || edges[2].intersects || edges[3].intersects || edges[4].intersects || edges[5].intersects || edges[6].intersects || edges[7].intersects || edges[8].intersects || edges[9].intersects || edges[10].intersects || edges[11].intersects)
if(cellZ >= 1 && cellY >= 1 && cellX >= 1)
{
cell.hasVertex = true;
//After the first rows and columns are done, start calculating vertex positions
int32_t cellXVertex = cellX-1;
int32_t cellYVertex = cellY-1;
int32_t cellZVertex = cellZ-1;
auto vertex = computeVertex(edges);
auto& cell = cells[convert(cellXVertex, cellYVertex, cellZVertex, cellRegionXDimension, cellRegionYDimension)];
vertex.setPosition(vertex.getPosition() + Vector3DFloat{static_cast<float>(cellX), static_cast<float>(cellY), static_cast<float>(cellZ)});
edges[0] = &cell.edges[0];
edges[1] = &cell.edges[1];
edges[2] = &cell.edges[2];
mesh.addVertex(vertex);
edges[3] = &cells[convert(cellXVertex+1, cellYVertex, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[1];
edges[4] = &cells[convert(cellXVertex+1, cellYVertex, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[2];
cell.vertexIndex = mesh.getNoOfVertices() - 1; //The index is of the last-added vertex
}
else
{
cell.hasVertex = false;
}
cells.setVoxel(cellX, cellY, cellZ, cell);
}
}
}
edges[5] = &cells[convert(cellXVertex, cellYVertex+1, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[2];
edges[6] = &cells[convert(cellXVertex, cellYVertex+1, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[0];
SimpleVolume<CellData>::Sampler cellSampler{&cells};
cellSampler.setPosition(0,0,0);
cellSampler.setWrapMode(WrapModes::Validate);
edges[7] = &cells[convert(cellXVertex, cellYVertex, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[0];
edges[8] = &cells[convert(cellXVertex, cellYVertex, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[1];
for(int32_t cellX = 1; cellX < cellRegion.getDimensionsInVoxels().getX(); cellX++)
{
for(int32_t cellY = 1; cellY < cellRegion.getDimensionsInVoxels().getY(); cellY++)
{
for(int32_t cellZ = 1; cellZ < cellRegion.getDimensionsInVoxels().getZ(); cellZ++)
{
cellSampler.setPosition(cellX, cellY, cellZ);
edges[9] = &cells[convert(cellXVertex+1, cellYVertex+1, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[2];
auto cell = cellSampler.getVoxel();
edges[10] = &cells[convert(cellXVertex+1, cellYVertex, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[1];
if(cell.x)
{
auto v0 = cell;
auto v1 = cellSampler.peekVoxel0px1ny0pz();
auto v2 = cellSampler.peekVoxel0px0py1nz();
auto v3 = cellSampler.peekVoxel0px1ny1nz();
mesh.addTriangle(v0.vertexIndex, v1.vertexIndex, v2.vertexIndex);
mesh.addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex);
}
edges[11] = &cells[convert(cellXVertex, cellYVertex+1, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[0];
if(cell.y)
{
auto v0 = cell;
auto v1 = cellSampler.peekVoxel1nx0py0pz();
auto v2 = cellSampler.peekVoxel0px0py1nz();
auto v3 = cellSampler.peekVoxel1nx0py1nz();
mesh.addTriangle(v0.vertexIndex, v1.vertexIndex, v2.vertexIndex);
mesh.addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex);
}
if(edges[0]->intersects || edges[1]->intersects || edges[2]->intersects || edges[3]->intersects || edges[4]->intersects || edges[5]->intersects || edges[6]->intersects || edges[7]->intersects || edges[8]->intersects || edges[9]->intersects || edges[10]->intersects || edges[11]->intersects) //'if' Maybe not needed?
{
auto vertex = computeVertex(edges);
if(cell.z)
{
auto v0 = cell;
auto v1 = cellSampler.peekVoxel1nx0py0pz();
auto v2 = cellSampler.peekVoxel0px1ny0pz();
auto v3 = cellSampler.peekVoxel1nx1ny0pz();
mesh.addTriangle(v0.vertexIndex, v1.vertexIndex, v2.vertexIndex);
mesh.addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex);
vertex.setPosition({vertex.getPosition().getX()+cellXVertex, vertex.getPosition().getY()+cellYVertex, vertex.getPosition().getZ()+cellZVertex});
cell.vertexIndex = mesh.addVertex(vertex);
if(cellZVertex >= 1 && cellYVertex >= 1 && cellXVertex >= 1)
{
//Once the second rows and colums are done, start connecting up edges
if(cell.edges[0].intersects)
{
auto v1 = cells[convert(cellXVertex, cellYVertex-1, cellZVertex, cellRegionXDimension, cellRegionYDimension)];
auto v2 = cells[convert(cellXVertex, cellYVertex, cellZVertex-1, cellRegionXDimension, cellRegionYDimension)];
auto v3 = cells[convert(cellXVertex, cellYVertex-1, cellZVertex-1, cellRegionXDimension, cellRegionYDimension)];
mesh.addTriangle(cell.vertexIndex, v1.vertexIndex, v2.vertexIndex);
mesh.addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex);
}
if(cell.edges[1].intersects)
{
auto v1 = cells[convert(cellXVertex-1, cellYVertex, cellZVertex, cellRegionXDimension, cellRegionYDimension)];
auto v2 = cells[convert(cellXVertex, cellYVertex, cellZVertex-1, cellRegionXDimension, cellRegionYDimension)];
auto v3 = cells[convert(cellXVertex-1, cellYVertex, cellZVertex-1, cellRegionXDimension, cellRegionYDimension)];
mesh.addTriangle(cell.vertexIndex, v1.vertexIndex, v2.vertexIndex);
mesh.addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex);
}
if(cell.edges[2].intersects)
{
auto v1 = cells[convert(cellXVertex-1, cellYVertex, cellZVertex, cellRegionXDimension, cellRegionYDimension)];
auto v2 = cells[convert(cellXVertex, cellYVertex-1, cellZVertex, cellRegionXDimension, cellRegionYDimension)];
auto v3 = cells[convert(cellXVertex-1, cellYVertex-1, cellZVertex, cellRegionXDimension, cellRegionYDimension)];
mesh.addTriangle(cell.vertexIndex, v1.vertexIndex, v2.vertexIndex);
mesh.addTriangle(v3.vertexIndex, v2.vertexIndex, v1.vertexIndex);
}
}
}
}
}
}
}
//logTrace() << "Dual contouring surface extraction took " << timer.elapsedTimeInMilliSeconds() << "ms (Region size = " << region.getWidthInVoxels() << "x" << region.getHeightInVoxels() << "x" << region.getDepthInVoxels() << ")";
logTrace() << "Dual contouring surface extraction took " << timer.elapsedTimeInMilliSeconds() << "ms (Region size = " << region.getWidthInVoxels() << "x" << region.getHeightInVoxels() << "x" << region.getDepthInVoxels() << ")";
std::cout << mesh.getNoOfVertices() << std::endl;
return mesh;
}
}

View File

@ -350,11 +350,6 @@ void qrstep(
return;
}
if (cols == 1) {
char *bomb = 0;
*bomb = 0;
}
// handle zeros on the diagonal or at its end
for (i = 0; i < cols - 1; ++i)
if (tau_u[i] == 0.0) {