diff --git a/library/PolyVoxCore/CMakeLists.txt b/library/PolyVoxCore/CMakeLists.txt index 4d4fbd43..580e6427 100644 --- a/library/PolyVoxCore/CMakeLists.txt +++ b/library/PolyVoxCore/CMakeLists.txt @@ -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 diff --git a/library/PolyVoxCore/include/PolyVoxCore/DualContouringSurfaceExtractor.inl b/library/PolyVoxCore/include/PolyVoxCore/DualContouringSurfaceExtractor.inl index b180c1a3..29b8180a 100644 --- a/library/PolyVoxCore/include/PolyVoxCore/DualContouringSurfaceExtractor.inl +++ b/library/PolyVoxCore/include/PolyVoxCore/DualContouringSurfaceExtractor.inl @@ -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 struct EdgeData { - VoxelType intersectionDistance; - VoxelType edgeLength; - bool intersects; Vector3DFloat normal; + float fraction; /// - Vector3DFloat computeCentralDifferenceGradient(const SamplerType& volIter) + template + 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(volIter.peekVoxel1nx0py0pz()); - float voxel1px = static_cast(volIter.peekVoxel1px0py0pz()); - - float voxel1ny = static_cast(volIter.peekVoxel0px1ny0pz()); - float voxel1py = static_cast(volIter.peekVoxel0px1py0pz()); - - float voxel1nz = static_cast(volIter.peekVoxel0px0py1nz()); - float voxel1pz = static_cast(volIter.peekVoxel0px0py1pz()); - - return Vector3DFloat - ( - voxel1nx - voxel1px, - voxel1ny - voxel1py, - voxel1nz - voxel1pz - ); - } + EdgeData edges[3]; + uint32_t vertexIndex; + }; template - EdgeData calculateEdge(VoxelType vA, VoxelType vB, Vector3DFloat gA, Vector3DFloat gB, ThresholdType threshold) + EdgeData calculateEdge(const VoxelType& vA, const VoxelType& vB, const Vector3DFloat& gA, const Vector3DFloat& gB, const ThresholdType& threshold) { EdgeData edge; - edge.intersectionDistance = vA - threshold; - edge.edgeLength = vA - vB; - auto fraction = static_cast(edge.intersectionDistance) / static_cast(edge.edgeLength); - if(fraction > 1.0 || fraction < 0.0) + edge.fraction = static_cast(vA - threshold) / static_cast(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(edge.edgeLength - edge.intersectionDistance) + gB * static_cast(edge.intersectionDistance)) / static_cast(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 - PositionMaterialNormal computeVertex(std::array, 12> edges) + PositionMaterialNormal computeVertex(EdgeData* edges[12]) { Vector3DFloat massPoint{0,0,0}; //The average of the intersection vertices - std::array vertices; + Vector3DFloat vertices[12]; - vertices[0] = {static_cast(edges[0].intersectionDistance) / static_cast(edges[0].edgeLength), 0, 0}; - vertices[1] = {0, static_cast(edges[1].intersectionDistance) / static_cast(edges[1].edgeLength), 0}; - vertices[2] = {0, 0, static_cast(edges[2].intersectionDistance) / static_cast(edges[2].edgeLength)}; - vertices[3] = {1, static_cast(edges[3].intersectionDistance) / static_cast(edges[3].edgeLength), 0}; - vertices[4] = {1, 0, static_cast(edges[4].intersectionDistance) / static_cast(edges[4].edgeLength)}; - vertices[5] = {0, 1, static_cast(edges[5].intersectionDistance) / static_cast(edges[5].edgeLength)}; - vertices[6] = {static_cast(edges[6].intersectionDistance) / static_cast(edges[6].edgeLength), 1, 0}; - vertices[7] = {static_cast(edges[7].intersectionDistance) / static_cast(edges[7].edgeLength), 0, 1}; - vertices[8] = {0, static_cast(edges[8].intersectionDistance) / static_cast(edges[8].edgeLength), 1}; - vertices[9] = {1, 1, static_cast(edges[9].intersectionDistance) / static_cast(edges[9].edgeLength)}; - vertices[10] = {1, static_cast(edges[10].intersectionDistance) / static_cast(edges[10].edgeLength), 1}; - vertices[11] = {static_cast(edges[11].intersectionDistance) / static_cast(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 @@ -178,8 +158,13 @@ namespace PolyVox { Timer timer; - Region cellRegion{Vector3DInt32{0,0,0}, region.getDimensionsInVoxels()}; - SimpleVolume 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> cells; + cells.reserve(cellRegionXDimension * cellRegionYDimension * cellRegionZDimension); typename VolumeType::Sampler volSampler{volData}; volSampler.setPosition(region.getLowerCorner()); @@ -189,139 +174,134 @@ namespace PolyVox SurfaceMesh mesh; - std::array, 12> edges; //Create this now but it will be overwritten for each cell + EdgeData* 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(volSampler.peekVoxel1nx0py0pz()); + //auto voxel1px = static_cast(volSampler.peekVoxel1px0py0pz()); + auto voxel1ny = static_cast(volSampler.peekVoxel0px1ny0pz()); + //auto voxel1py = static_cast(volSampler.peekVoxel0px1py0pz()); + auto voxel1nz = static_cast(volSampler.peekVoxel0px0py1nz()); + //auto voxel1pz = static_cast(volSampler.peekVoxel0px0py1pz()); + Vector3DFloat g000(voxel1nx - v100, voxel1ny - v010, voxel1nz - v001); + volSampler.movePositiveX(); + auto voxel2px = static_cast(volSampler.peekVoxel1px0py0pz()); + auto voxel1px1ny = static_cast(volSampler.peekVoxel0px1ny0pz()); + auto voxel1px1py = static_cast(volSampler.peekVoxel0px1py0pz()); + auto voxel1px1nz = static_cast(volSampler.peekVoxel0px0py1nz()); + auto voxel1px1pz = static_cast(volSampler.peekVoxel0px0py1pz()); + Vector3DFloat g100(v000 - voxel2px, voxel1px1ny - voxel1px1py, voxel1px1nz - voxel1px1pz); + volSampler.moveNegativeX(); + volSampler.movePositiveY(); + auto voxel1nx1py = static_cast(volSampler.peekVoxel1nx0py0pz()); + auto voxel2py = static_cast(volSampler.peekVoxel0px1py0pz()); + auto voxel1py1nz = static_cast(volSampler.peekVoxel0px0py1nz()); + auto voxel1py1pz = static_cast(volSampler.peekVoxel0px0py1pz()); + Vector3DFloat g010(voxel1nx1py - voxel1px1py, v000 - voxel2py, voxel1py1nz - voxel1py1pz); + volSampler.moveNegativeY(); + volSampler.movePositiveZ(); + auto voxel1nx1pz = static_cast(volSampler.peekVoxel1nx0py0pz()); + auto voxel1ny1pz = static_cast(volSampler.peekVoxel0px1ny0pz()); + auto voxel2pz = static_cast(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(cellX), static_cast(cellY), static_cast(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 + edges[5] = &cells[convert(cellXVertex, cellYVertex+1, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[2]; + edges[6] = &cells[convert(cellXVertex, cellYVertex+1, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[0]; + + edges[7] = &cells[convert(cellXVertex, cellYVertex, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[0]; + edges[8] = &cells[convert(cellXVertex, cellYVertex, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[1]; + + edges[9] = &cells[convert(cellXVertex+1, cellYVertex+1, cellZVertex, cellRegionXDimension, cellRegionYDimension)].edges[2]; + + edges[10] = &cells[convert(cellXVertex+1, cellYVertex, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[1]; + + edges[11] = &cells[convert(cellXVertex, cellYVertex+1, cellZVertex+1, cellRegionXDimension, cellRegionYDimension)].edges[0]; + + 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); + + 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); + } + } + } } - else - { - cell.hasVertex = false; - } - cells.setVoxel(cellX, cellY, cellZ, cell); } } } - SimpleVolume::Sampler cellSampler{&cells}; - cellSampler.setPosition(0,0,0); - cellSampler.setWrapMode(WrapModes::Validate); - - 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); - - auto cell = cellSampler.getVoxel(); - - 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); - } - - 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(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); - } - } - } - } + //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; } } diff --git a/library/PolyVoxCore/include/PolyVoxCore/Impl/QEF.inl b/library/PolyVoxCore/include/PolyVoxCore/Impl/QEF.inl index 9cbb65db..ad4baed3 100644 --- a/library/PolyVoxCore/include/PolyVoxCore/Impl/QEF.inl +++ b/library/PolyVoxCore/include/PolyVoxCore/Impl/QEF.inl @@ -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) {