/******************************************************************************* * The MIT License (MIT) * * Copyright (c) 2015 David Williams and Matthew Williams * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in all * copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. *******************************************************************************/ #include "Impl/Timer.h" namespace PolyVox { //////////////////////////////////////////////////////////////////////////////// // Vertex encoding/decoding //////////////////////////////////////////////////////////////////////////////// inline Vector3DFloat decodePosition(const Vector3DUint16& encodedPosition) { Vector3DFloat result(encodedPosition.getX(), encodedPosition.getY(), encodedPosition.getZ()); result *= (1.0f / 256.0f); // Division is compile-time constant return result; } inline uint16_t encodeNormal(const Vector3DFloat& normal) { // The first part of this function is based off the code in Listing 1 of http://jcgt.org/published/0003/02/01/ // It was rewritten in C++ and is restructued for the CPU rather than the GPU. // Get the input components float vx = normal.getX(); float vy = normal.getY(); float vz = normal.getZ(); // Project the sphere onto the octahedron, and then onto the xy plane float px = vx * (1.0f / (std::abs(vx) + std::abs(vy) + std::abs(vz))); float py = vy * (1.0f / (std::abs(vx) + std::abs(vy) + std::abs(vz))); // Reflect the folds of the lower hemisphere over the diagonals. if (vz <= 0.0f) { float refx = ((1.0f - std::abs(py)) * (px >= 0.0f ? +1.0f : -1.0f)); float refy = ((1.0f - std::abs(px)) * (py >= 0.0f ? +1.0f : -1.0f)); px = refx; py = refy; } // The next part was not given in the paper. We map our two // floats into two bytes and store them in a single uint16_t // Move from range [-1.0f, 1.0f] to [0.0f, 255.0f] px = (px + 1.0f) * 127.5f; py = (py + 1.0f) * 127.5f; // Convert to uints uint16_t resultX = static_cast(px + 0.5f); uint16_t resultY = static_cast(py + 0.5f); // Make sure only the lower bits are set. Probably // not necessary but we're just being careful really. resultX &= 0xFF; resultY &= 0xFF; // Contatenate the bytes and return the result. return (resultX << 8) | resultY; } inline Vector3DFloat decodeNormal(const uint16_t& encodedNormal) { // Extract the two bytes from the uint16_t. uint16_t ux = (encodedNormal >> 8) & 0xFF; uint16_t uy = (encodedNormal)& 0xFF; // Convert to floats in the range [-1.0f, +1.0f]. float ex = ux / 127.5f - 1.0f; float ey = uy / 127.5f - 1.0f; // Reconstruct the origninal vector. This is a C++ implementation // of Listing 2 of http://jcgt.org/published/0003/02/01/ float vx = ex; float vy = ey; float vz = 1.0f - std::abs(ex) - std::abs(ey); if (vz < 0.0f) { float refX = ((1.0f - std::abs(vy)) * (vx >= 0.0f ? +1.0f : -1.0f)); float refY = ((1.0f - std::abs(vx)) * (vy >= 0.0f ? +1.0f : -1.0f)); vx = refX; vy = refY; } // Normalise and return the result. Vector3DFloat v(vx, vy, vz); v.normalise(); return v; } template Vertex decodeVertex(const MarchingCubesVertex& marchingCubesVertex) { Vertex result; result.position = decodePosition(marchingCubesVertex.encodedPosition); result.normal = decodeNormal(marchingCubesVertex.encodedNormal); result.data = marchingCubesVertex.data; // Data is not encoded return result; } //////////////////////////////////////////////////////////////////////////////// // Gradient estimation //////////////////////////////////////////////////////////////////////////////// template< typename Sampler, typename ControllerType> Vector3DFloat computeCentralDifferenceGradient(const Sampler& volIter, ControllerType& controller) { //FIXME - Should actually use DensityType here, both in principle and because the maths may be //faster (and to reduce casts). So it would be good to add a way to get DensityType from a voxel. //But watch out for when the DensityType is unsigned and the difference could be negative. float voxel1nx = static_cast(controller.convertToDensity(volIter.peekVoxel1nx0py0pz())); float voxel1px = static_cast(controller.convertToDensity(volIter.peekVoxel1px0py0pz())); float voxel1ny = static_cast(controller.convertToDensity(volIter.peekVoxel0px1ny0pz())); float voxel1py = static_cast(controller.convertToDensity(volIter.peekVoxel0px1py0pz())); float voxel1nz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py1nz())); float voxel1pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py1pz())); return Vector3DFloat ( voxel1nx - voxel1px, voxel1ny - voxel1py, voxel1nz - voxel1pz ); } // This 'sobel' version of gradient estimation provides better (smoother) normals than the central difference version. // Even with the 16-bit normal encoding it does seem to make a difference, so is probably worth keeping. However, there // is no way to call it at the moment beyond modifying the main Marching Cubes function below to call this function // instead of the central difference one. We should provide a way to control the normal generation method, perhaps // including *no* normals incase the user wants to generate them afterwards (e.g. from the mesh). template< typename Sampler, typename ControllerType> Vector3DFloat computeSobelGradient(const Sampler& volIter, ControllerType& controller) { 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(controller.convertToDensity(volIter.peekVoxel1nx1ny1nz())); const float pVoxel1nx1ny0pz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx1ny0pz())); const float pVoxel1nx1ny1pz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx1ny1pz())); const float pVoxel1nx0py1nz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx0py1nz())); const float pVoxel1nx0py0pz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx0py0pz())); const float pVoxel1nx0py1pz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx0py1pz())); const float pVoxel1nx1py1nz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx1py1nz())); const float pVoxel1nx1py0pz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx1py0pz())); const float pVoxel1nx1py1pz = static_cast(controller.convertToDensity(volIter.peekVoxel1nx1py1pz())); const float pVoxel0px1ny1nz = static_cast(controller.convertToDensity(volIter.peekVoxel0px1ny1nz())); const float pVoxel0px1ny0pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px1ny0pz())); const float pVoxel0px1ny1pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px1ny1pz())); const float pVoxel0px0py1nz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py1nz())); //const float pVoxel0px0py0pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py0pz())); const float pVoxel0px0py1pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px0py1pz())); const float pVoxel0px1py1nz = static_cast(controller.convertToDensity(volIter.peekVoxel0px1py1nz())); const float pVoxel0px1py0pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px1py0pz())); const float pVoxel0px1py1pz = static_cast(controller.convertToDensity(volIter.peekVoxel0px1py1pz())); const float pVoxel1px1ny1nz = static_cast(controller.convertToDensity(volIter.peekVoxel1px1ny1nz())); const float pVoxel1px1ny0pz = static_cast(controller.convertToDensity(volIter.peekVoxel1px1ny0pz())); const float pVoxel1px1ny1pz = static_cast(controller.convertToDensity(volIter.peekVoxel1px1ny1pz())); const float pVoxel1px0py1nz = static_cast(controller.convertToDensity(volIter.peekVoxel1px0py1nz())); const float pVoxel1px0py0pz = static_cast(controller.convertToDensity(volIter.peekVoxel1px0py0pz())); const float pVoxel1px0py1pz = static_cast(controller.convertToDensity(volIter.peekVoxel1px0py1pz())); const float pVoxel1px1py1nz = static_cast(controller.convertToDensity(volIter.peekVoxel1px1py1nz())); const float pVoxel1px1py0pz = static_cast(controller.convertToDensity(volIter.peekVoxel1px1py0pz())); const float pVoxel1px1py1pz = static_cast(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); } //////////////////////////////////////////////////////////////////////////////// // Surface extraction //////////////////////////////////////////////////////////////////////////////// /// This is probably the version of Marching Cubes extraction which you will want to use initially, at least /// until you determine you have a need for the extra functionality provied by extractMarchingCubesMeshCustom(). template< typename VolumeType, typename ControllerType > Mesh > extractMarchingCubesMesh(VolumeType* volData, Region region, ControllerType controller) { Mesh > result; extractMarchingCubesMeshCustom, DefaultIndexType > >(volData, region, &result, controller); return result; } /// 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. /// /// 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. /// /// Note: This function is called 'extractMarchingCubesMeshCustom' rather than 'extractMarchingCubesMesh' to avoid ambiguity when only three parameters /// 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 > 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; // Performance note: Profiling indicates that simply adding vertices and indices to the std::vector is one // of the bottlenecks when generating the mesh. Reserving space in advance helps here but is wasteful in the // common case that no/few vertices are generated. Maybe it's worth reserving a couple of thousand or so? // Alternatively, maybe the docs should suggest the user reserves some space in the mesh they pass in? result->clear(); // 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(); typename ControllerType::DensityType tThreshold = controller.getThreshold(); // 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); // 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); // 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()); for (uint32_t uZRegSpace = 0; uZRegSpace < uRegionDepthInVoxels; uZRegSpace++) { // 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; for (uint32_t uYRegSpace = 0; uYRegSpace < uRegionHeightInVoxels; uYRegSpace++) { // 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; 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. // Each bit of the cell index specifies whether a given corner of the cell is above or below the threshold. uint8_t uCellIndex = 0; // 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; // 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); // Performance note: Computing normals is one of the bottlencks in the mesh generation process. The // central difference approach actually samples the same voxel more than once as we call it on two // adjacent voxels. Perhaps we could expand this and eliminate dupicates in the future. Alternatively, // we could compute vertex normals from adjacent face normals instead of via central differencing, // but not for vertices on the edge of the region (as this causes visual discontinities). const Vector3DFloat n111 = computeCentralDifferenceGradient(sampler, controller); /* Find the vertices where the surface intersects the cube */ if ((uEdge & 64) && (uXRegSpace > 0)) { sampler.moveNegativeX(); typename VolumeType::VoxelType v011 = sampler.getVoxel(); auto v011Density = controller.convertToDensity(v011); const float fInterp = static_cast(tThreshold - v011Density) / static_cast(v111Density - v011Density); // Compute the position const Vector3DFloat v3dPosition(static_cast(uXRegSpace - 1) + fInterp, static_cast(uYRegSpace), static_cast(uZRegSpace)); // Compute the normal const Vector3DFloat n011 = computeCentralDifferenceGradient(sampler, controller); Vector3DFloat v3dNormal = (n111*fInterp) + (n011*(1 - fInterp)); // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). if (v3dNormal.lengthSquared() > 0.000001f) { v3dNormal.normalise(); } // Allow the controller to decide how the material should be derived from the voxels. const typename VolumeType::VoxelType uMaterial = controller.blendMaterials(v011, v111, fInterp); MarchingCubesVertex surfaceVertex; const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); surfaceVertex.encodedPosition = v3dScaledPosition; surfaceVertex.encodedNormal = encodeNormal(v3dNormal); surfaceVertex.data = uMaterial; const uint32_t uLastVertexIndex = result->addVertex(surfaceVertex); pIndices(uXRegSpace, uYRegSpace).setX(uLastVertexIndex); sampler.movePositiveX(); } if ((uEdge & 32) && (uYRegSpace > 0)) { sampler.moveNegativeY(); typename VolumeType::VoxelType v101 = sampler.getVoxel(); auto v101Density = controller.convertToDensity(v101); const float fInterp = static_cast(tThreshold - v101Density) / static_cast(v111Density - v101Density); // Compute the position const Vector3DFloat v3dPosition(static_cast(uXRegSpace), static_cast(uYRegSpace - 1) + fInterp, static_cast(uZRegSpace)); // Compute the normal const Vector3DFloat n101 = computeCentralDifferenceGradient(sampler, controller); Vector3DFloat v3dNormal = (n111*fInterp) + (n101*(1 - fInterp)); // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). if (v3dNormal.lengthSquared() > 0.000001f) { v3dNormal.normalise(); } // Allow the controller to decide how the material should be derived from the voxels. const typename VolumeType::VoxelType uMaterial = controller.blendMaterials(v101, v111, fInterp); MarchingCubesVertex surfaceVertex; const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); surfaceVertex.encodedPosition = v3dScaledPosition; surfaceVertex.encodedNormal = encodeNormal(v3dNormal); surfaceVertex.data = uMaterial; uint32_t uLastVertexIndex = result->addVertex(surfaceVertex); pIndices(uXRegSpace, uYRegSpace).setY(uLastVertexIndex); sampler.movePositiveY(); } if ((uEdge & 1024) && (uZRegSpace > 0)) { sampler.moveNegativeZ(); typename VolumeType::VoxelType v110 = sampler.getVoxel(); auto v110Density = controller.convertToDensity(v110); const float fInterp = static_cast(tThreshold - v110Density) / static_cast(v111Density - v110Density); // Compute the position const Vector3DFloat v3dPosition(static_cast(uXRegSpace), static_cast(uYRegSpace), static_cast(uZRegSpace - 1) + fInterp); // Compute the normal const Vector3DFloat n110 = computeCentralDifferenceGradient(sampler, controller); Vector3DFloat v3dNormal = (n111*fInterp) + (n110*(1 - fInterp)); // The gradient for a voxel can be zero (e.g. solid voxel surrounded by empty ones) and so // the interpolated normal can also be zero (e.g. a grid of alternating solid and empty voxels). if (v3dNormal.lengthSquared() > 0.000001f) { v3dNormal.normalise(); } // Allow the controller to decide how the material should be derived from the voxels. const typename VolumeType::VoxelType uMaterial = controller.blendMaterials(v110, v111, fInterp); MarchingCubesVertex surfaceVertex; const Vector3DUint16 v3dScaledPosition(static_cast(v3dPosition.getX() * 256.0f), static_cast(v3dPosition.getY() * 256.0f), static_cast(v3dPosition.getZ() * 256.0f)); surfaceVertex.encodedPosition = v3dScaledPosition; surfaceVertex.encodedNormal = encodeNormal(v3dNormal); surfaceVertex.data = uMaterial; const uint32_t uLastVertexIndex = result->addVertex(surfaceVertex); pIndices(uXRegSpace, uYRegSpace).setZ(uLastVertexIndex); sampler.movePositiveZ(); } // Now output the indices. For the first row, column or slice there aren't // any (the region size in cells is one less than the region size in voxels) if ((uXRegSpace != 0) && (uYRegSpace != 0) && (uZRegSpace != 0)) { int32_t indlist[12]; /* Find the vertices where the surface intersects the cube */ if (uEdge & 1) { indlist[0] = pPreviousIndices(uXRegSpace, uYRegSpace - 1).getX(); } if (uEdge & 2) { indlist[1] = pPreviousIndices(uXRegSpace, uYRegSpace).getY(); } if (uEdge & 4) { indlist[2] = pPreviousIndices(uXRegSpace, uYRegSpace).getX(); } if (uEdge & 8) { indlist[3] = pPreviousIndices(uXRegSpace - 1, uYRegSpace).getY(); } if (uEdge & 16) { indlist[4] = pIndices(uXRegSpace, uYRegSpace - 1).getX(); } if (uEdge & 32) { indlist[5] = pIndices(uXRegSpace, uYRegSpace).getY(); } if (uEdge & 64) { indlist[6] = pIndices(uXRegSpace, uYRegSpace).getX(); } if (uEdge & 128) { indlist[7] = pIndices(uXRegSpace - 1, uYRegSpace).getY(); } if (uEdge & 256) { indlist[8] = pIndices(uXRegSpace - 1, uYRegSpace - 1).getZ(); } if (uEdge & 512) { indlist[9] = pIndices(uXRegSpace, uYRegSpace - 1).getZ(); } if (uEdge & 1024) { indlist[10] = pIndices(uXRegSpace, uYRegSpace).getZ(); } if (uEdge & 2048) { indlist[11] = pIndices(uXRegSpace - 1, uYRegSpace).getZ(); } for (int i = 0; triTable[uCellIndex][i] != -1; i += 3) { const int32_t ind0 = indlist[triTable[uCellIndex][i]]; const int32_t ind1 = indlist[triTable[uCellIndex][i + 1]]; const int32_t ind2 = indlist[triTable[uCellIndex][i + 2]]; if ((ind0 != -1) && (ind1 != -1) && (ind2 != -1)) { result->addTriangle(ind0, ind1, ind2); } } // For each triangle } } // For each cell sampler.movePositiveX(); } // For X startOfRow.movePositiveY(); } // For Y startOfSlice.movePositiveZ(); pIndices.swap(pPreviousIndices); } // For Z result->setOffset(region.getLowerCorner()); POLYVOX_LOG_TRACE("Marching cubes surface extraction took ", timer.elapsedTimeInMilliSeconds(), "ms (Region size = ", region.getWidthInVoxels(), "x", region.getHeightInVoxels(), "x", region.getDepthInVoxels(), ")"); } }