BRL-CAD/GeometricTools
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed problems in SurfaceExtractorMC.h

Browse Source 
When the single-voxel Extract returns false, the full-image extraction was supposed to ignore the voxel. Instead it deleted the entire mesh, which is a disastrous bug. The unit tests verified only that the Marching Cubes 256 triangle-lookup table entries are correct, and none of these provided code coverage that would have found the bug. Also, floating-point rounding errors in mesh vertex construction led to small sets of vertices, each set containing vertices that are nearly the same but should have been the same without the rounding errors. MakeUnique is unable to identify these as duplicates. I modified the single-voxel Extract function to generate the same vertex on an edge shared by multiple voxels, using a consistent ordering of the voxel corners to ensure that the all constructions of the vertex have the same rounding errors. Finally, the interfaces were modified to have a new parameter 'perturb'
that allows you to avoid the ignoring of voxels for which at least one corner of the voxel is exactly the input 'level' value. Thanks to Scott Johnson for reporting these problems.
This commit is contained in:
deber
2023-09-23 14:09:14 -07:00
parent 2dfbd18579
commit 3d97d20043
3 changed files with 7658 additions and 7486 deletions
Download Patch File Download Diff File Expand all files Collapse all files
GTE/Gte6UpdateHistory.pdf
BIN
View File
Binary file not shown.
GTE/GteFullUpdateHistory.pdf
+7485 -7373
View File
File diff suppressed because it is too large Load Diff
GTE/Mathematics/SurfaceExtractorMC.h
+173 -113
View File
@@ -3,7 +3,7 @@
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 6.0.2023.08.08
// Version: 6.0.2023.09.22
#pragma once
@@ -21,7 +21,7 @@
namespace gte
{
template <typename Real>
template <typename T>
class SurfaceExtractorMC : public MarchingCubes
{
public:
@@ -30,7 +30,7 @@ namespace gte
{
}
SurfaceExtractorMC(Image3<Real> const& image)
SurfaceExtractorMC(Image3<T> const& image)
:
mImage(image)
{
@@ -46,34 +46,66 @@ namespace gte
// All members are set to zeros.
Mesh()
{
std::array<Real, 3> zero = { (Real)0, (Real)0, (Real)0 };
std::array<T, 3> zero = { (T)0, (T)0, (T)0 };
std::fill(vertices.begin(), vertices.end(), zero);
}
Topology topology;
std::array<Vector3<Real>, MAX_VERTICES> vertices;
std::array<Vector3<T>, MAX_VERTICES> vertices;
};
// Extract the triangle mesh approximating F = 0 for a single voxel.
// The input function values must be stored as
// F[0] = function(0,0,0), F[4] = function(0,0,1),
// F[1] = function(1,0,0), F[5] = function(1,0,1),
// F[2] = function(0,1,0), F[6] = function(0,1,1),
// F[3] = function(1,1,0), F[7] = function(1,1,1).
// Thus, F[k] = function(k & 1, (k & 2) >> 1, (k & 4) >> 2).
// The return value is 'true' iff the F[] values are all nonzero.
// If they are not, the returned 'mesh' has no vertices and no
// triangles--as if F[] had all positive or all negative values.
bool Extract(std::array<Real, 8> const& F, Mesh& mesh) const
// Extract the triangle mesh approximating F = level for a single
// voxel whose origin corner is (x,y,z). The image has dimensions
// d0, d1, and d2, and the origin corner satisfies the conditions
// 0 <= x < d0-1, 0 <= y < d1-1, and 0 <= z < d2-1. The input
// function values must be stored as
// F[0] = image(x ,y ,z), F[4] = image(x ,y ,z+1),
// F[1] = image(x+1,y ,z), F[5] = image(x+1,y ,z+1),
// F[2] = image(x ,y+1,z), F[6] = image(x ,y+1,z+1),
// F[3] = image(x+1,y+1,z), F[7] = image(x+1,y+1,z+1)
// In local coordinates where the corners are (0,0,0), (1,0,0),
// (0,1,0), (1,1,0), (0,0,1), (1,0,1), (0,1,1), and (1,1,1),
// Thus, F[k] = imageLocal(k & 1, (k & 2) >> 1, (k & 4) >> 2).
// The caller of this function must add in the (x,y,z) origin corner
// to the mesh.vertices[] to obtain the global coordinates of the
// vertices.
//
// The return value is 'true' iff the F[] values are all not equal
// to 'level'. If at least one of F[] is 'level', the returned
// 'mesh' has no vertices and no triangles. If you want this behavior,
// set 'perturb' to zero.
//
// If you want to avoid the case when F[i] = level for some i, set
// 'perturb' to a small nonzero number whose absolute value is smaller
// than the minimum absolute value of the differences between voxel
// values and 'level'.
bool Extract(T const& level, T const& perturb, std::array<T, 8> const& F,
Mesh& mesh) const
{
T const zero = static_cast<T>(0);
T const one = static_cast<T>(1);
std::array<T, 8> localF{};
int32_t entry = 0;
for (int32_t i = 0, mask = 1; i < 8; ++i, mask <<= 1)
{
if (F[i] < (Real)0)
localF[i] = F[i] - level;
if (localF[i] < zero)
{
entry |= mask;
}
else if (F[i] == (Real)0)
else if (localF[i] > zero)
{
continue;
}
// If 'perturb' is zero, the function will report that no
// geometry is generated for this voxel. If 'perturb' is
// not zero, the comparison to zero still needs to be made
// after the perturbation in case floating-point rounding
// errors cause localF[i] still to be zero.
localF[i] += perturb;
if (localF[i] == zero)
{
return false;
}
@@ -86,63 +118,97 @@ namespace gte
int32_t j0 = mesh.topology.vpair[i][0];
int32_t j1 = mesh.topology.vpair[i][1];
std::array<Real, 3> corner0{};
corner0[0] = static_cast<Real>(j0 & 1);
corner0[1] = static_cast<Real>((j0 & 2) >> 1);
corner0[2] = static_cast<Real>((j0 & 4) >> 2);
std::array<Real, 3> corner1{};
corner1[0] = static_cast<Real>(j1 & 1);
corner1[1] = static_cast<Real>((j1 & 2) >> 1);
corner1[2] = static_cast<Real>((j1 & 4) >> 2);
Real invDenom = ((Real)1) / (F[j0] - F[j1]);
for (int32_t k = 0; k < 3; ++k)
// The vertex can be computed with 3D-only computations as
// V = (F[j0] * k1 - F[j1] * k0) / (F[j1] - F[j0])
// but floating-point rounding errors can cause at least one
// of the integer-valued components of V not to be 0 or 1 as
// the case may be. Such errors in turn can lead to 2 nearly
// identical vertices in the mesh, and MakeUnique will not be
// able to characterize them as the same vertex. The
// componentwise computations avoid these floating-point
// rounding errors. It is guaranteed that j0 < j1, so multiple
// voxels sharing the same edge will generate the same vertex.
auto& vertex = mesh.vertices[i];
std::array<int32_t, 3> k0 = { j0 & 1, (j0 & 2) >> 1, (j0 & 4) >> 2 };
std::array<int32_t, 3> k1 = { j1 & 1, (j1 & 2) >> 1, (j1 & 4) >> 2 };
for (int32_t index = 0; index < 3; ++index)
{
Real numer = F[j0] * corner1[k] - F[j1] * corner0[k];
mesh.vertices[i][k] = numer * invDenom;
if (k0[index] == 0)
{
if (k1[index] == 0)
{
vertex[index] = zero;
}
else // k1[index] = 1
{
vertex[index] = F[j0] / (F[j0] - F[j1]);
}
}
else // k0[index] = 1
{
if (k1[index] == 0)
{
vertex[index] = F[j1] / (F[j1] - F[j0]);
}
else // k1[index] = 1
{
vertex[index] = one;
}
}
}
}
return true;
}
// Extract the triangle mesh approximating F = 0 for all the voxels in
// a 3D image. The input image must be stored in a 1-dimensional
// array with lexicographical order; that is, image[i] corresponds to
// voxel location (x,y,z) where i = x + bound0 * (y + bound1 * z).
// The output 'indices' consists indices.size()/3 triangles, each a
// triple of indices into 'vertices'
bool Extract(Real level, std::vector<Vector3<Real>>& vertices, std::vector<int32_t>& indices) const
// Extract the triangle mesh approximating F = level for all the
// voxels in a 3D image. The input image must be stored in a
// 1-dimensional array with lexicographical order. If the image
// has dimensions d0, d1, and d2, then voxel location (x, y, z)
// contains value image[x + d0 * (y + d1 * z)]. The output 'indices'
// consists indices.size()/3 triangles, each a triple of indices into
// 'vertices'
//
// The triangle table lookups depend on voxel values never being
// exactly equal to 'level'. Set 'perturb' to zero so that any
// voxel cube with at least one corner value equal to 'level' is
// ignored in the final mesh; that is, such a voxel does not generate
// any triangles in the final mesh. If you want triangles from such
// voxels, set 'perturb' to a small nonzero number whose absolute
// value is smaller than the minimum absolute value of the differences
// between voxel values and 'level'.
void Extract(T const& level, T const& perturb,
std::vector<Vector3<T>>& vertices, std::vector<int32_t>& indices) const
{
vertices.clear();
indices.clear();
for (int32_t z = 0; z + 1 < mImage.GetDimension(2); ++z)
for (int32_t z0 = 0, z1 = 1; z1 < mImage.GetDimension(2); z0 = z1++)
{
for (int32_t y = 0; y + 1 < mImage.GetDimension(1); ++y)
for (int32_t y0 = 0, y1 = 1; y1 < mImage.GetDimension(1); y0 = y1++)
{
for (int32_t x = 0; x + 1 < mImage.GetDimension(0); ++x)
for (int32_t x0 = 0, x1 = 1; x1 < mImage.GetDimension(0); x0 = x1++)
{
std::array<size_t, 8> corners{};
mImage.GetCorners(x, y, z, corners);
std::array<Real, 8> F{};
for (int32_t k = 0; k < 8; ++k)
{
F[k] = mImage[corners[k]] - level;
}
std::array<T, 8> F{};
F[0] = mImage(x0, y0, z0);
F[1] = mImage(x1, y0, z0);
F[2] = mImage(x0, y1, z0);
F[3] = mImage(x1, y1, z0);
F[4] = mImage(x0, y0, z1);
F[5] = mImage(x1, y0, z1);
F[6] = mImage(x0, y1, z1);
F[7] = mImage(x1, y1, z1);
Mesh mesh{};
if (Extract(F, mesh))
if (Extract(level, perturb, F, mesh))
{
int32_t vbase = static_cast<int32_t>(vertices.size());
for (int32_t i = 0; i < mesh.topology.numVertices; ++i)
{
Vector3<Real> position = mesh.vertices[i];
position[0] += static_cast<Real>(x);
position[1] += static_cast<Real>(y);
position[2] += static_cast<Real>(z);
Vector3<T> position = mesh.vertices[i];
position[0] += static_cast<T>(x0);
position[1] += static_cast<T>(y0);
position[2] += static_cast<T>(z0);
vertices.push_back(position);
}
@@ -154,26 +220,18 @@ namespace gte
}
}
}
else
{
vertices.clear();
indices.clear();
return false;
}
}
}
}
return true;
}
// The extraction has duplicate vertices on edges shared by voxels.
// This function will eliminate the duplication.
void MakeUnique(std::vector<Vector3<Real>>& vertices, std::vector<int32_t>& indices) const
void MakeUnique(std::vector<Vector3<T>>& vertices, std::vector<int32_t>& indices) const
{
std::vector<Vector3<Real>> outVertices;
std::vector<int32_t> outIndices;
UniqueVerticesSimplices<Vector3<Real>, int32_t, 3> uvt;
std::vector<Vector3<T>> outVertices{};
std::vector<int32_t> outIndices{};
UniqueVerticesSimplices<Vector3<T>, int32_t, 3> uvt{};
uvt.RemoveDuplicateVertices(vertices, indices, outVertices, outIndices);
vertices = std::move(outVertices);
indices = std::move(outIndices);
@@ -190,37 +248,38 @@ namespace gte
// build a table of vertex, edge, and face adjacencies, but the
// resulting data structure is somewhat expensive to process to
// reorient triangles.
void OrientTriangles(std::vector<Vector3<Real>> const& vertices, std::vector<int32_t>& indices, bool sameDir) const
void OrientTriangles(std::vector<Vector3<T>> const& vertices, std::vector<int32_t>& indices, bool sameDir) const
{
T const zero = static_cast<T>(0);
int32_t const numTriangles = static_cast<int32_t>(indices.size() / 3);
int32_t* triangle = indices.data();
for (int32_t t = 0; t < numTriangles; ++t, triangle += 3)
{
// Get triangle vertices.
Vector3<Real> v0 = vertices[triangle[0]];
Vector3<Real> v1 = vertices[triangle[1]];
Vector3<Real> v2 = vertices[triangle[2]];
Vector3<T> v0 = vertices[triangle[0]];
Vector3<T> v1 = vertices[triangle[1]];
Vector3<T> v2 = vertices[triangle[2]];
// Construct triangle normal based on current orientation.
Vector3<Real> edge1 = v1 - v0;
Vector3<Real> edge2 = v2 - v0;
Vector3<Real> normal = Cross(edge1, edge2);
Vector3<T> edge1 = v1 - v0;
Vector3<T> edge2 = v2 - v0;
Vector3<T> normal = Cross(edge1, edge2);
// Get the image gradient at the vertices.
Vector3<Real> gradient0 = GetGradient(v0);
Vector3<Real> gradient1 = GetGradient(v1);
Vector3<Real> gradient2 = GetGradient(v2);
Vector3<T> gradient0 = GetGradient(v0);
Vector3<T> gradient1 = GetGradient(v1);
Vector3<T> gradient2 = GetGradient(v2);
// Compute the average gradient.
Vector3<Real> gradientAvr = (gradient0 + gradient1 + gradient2) / (Real)3;
Vector3<T> gradientAvr = (gradient0 + gradient1 + gradient2) / (T)3;
// Compute the dot product of normal and average gradient.
Real dot = Dot(gradientAvr, normal);
T dot = Dot(gradientAvr, normal);
// Choose triangle orientation based on gradient direction.
if (sameDir)
{
if (dot < (Real)0)
if (dot < zero)
{
// Wrong orientation, reorder it.
std::swap(triangle[1], triangle[2]);
@@ -228,7 +287,7 @@ namespace gte
}
else
{
if (dot > (Real)0)
if (dot > zero)
{
// Wrong orientation, reorder it.
std::swap(triangle[1], triangle[2]);
@@ -238,13 +297,13 @@ namespace gte
}
// Compute vertex normals for the mesh.
void ComputeNormals(std::vector<Vector3<Real>> const& vertices, std::vector<int32_t> const& indices,
std::vector<Vector3<Real>>& normals) const
void ComputeNormals(std::vector<Vector3<T>> const& vertices, std::vector<int32_t> const& indices,
std::vector<Vector3<T>>& normals) const
{
// Maintain a running sum of triangle normals at each vertex.
int32_t const numVertices = static_cast<int32_t>(vertices.size());
normals.resize(numVertices);
Vector3<Real> zero = Vector3<Real>::Zero();
Vector3<T> zero = Vector3<T>::Zero();
std::fill(normals.begin(), normals.end(), zero);
int32_t const numTriangles = static_cast<int32_t>(indices.size() / 3);
@@ -254,14 +313,14 @@ namespace gte
int32_t i0 = *current++;
int32_t i1 = *current++;
int32_t i2 = *current++;
Vector3<Real> v0 = vertices[i0];
Vector3<Real> v1 = vertices[i1];
Vector3<Real> v2 = vertices[i2];
Vector3<T> v0 = vertices[i0];
Vector3<T> v1 = vertices[i1];
Vector3<T> v2 = vertices[i2];
// Construct triangle normal.
Vector3<Real> edge1 = v1 - v0;
Vector3<Real> edge2 = v2 - v0;
Vector3<Real> normal = Cross(edge1, edge2);
Vector3<T> edge1 = v1 - v0;
Vector3<T> edge2 = v2 - v0;
Vector3<T> normal = Cross(edge1, edge2);
// Maintain the sum of normals at each vertex.
normals[i0] += normal;
@@ -270,7 +329,7 @@ namespace gte
}
// The normal vector storage was used to accumulate the sum of
// triangle normals. Now these vectors must be rescaled to be
// triangle normals. Now these vectors must be rescaled to be
// unit length.
for (auto& normal : normals)
{
@@ -279,49 +338,50 @@ namespace gte
}
protected:
Vector3<Real> GetGradient(Vector3<Real> position) const
Vector3<T> GetGradient(Vector3<T> position) const
{
int32_t x = static_cast<int32_t>(std::floor(position[0]));
if (x < 0 || x >= mImage.GetDimension(0) - 1)
{
return Vector3<Real>::Zero();
return Vector3<T>::Zero();
}
int32_t y = static_cast<int32_t>(std::floor(position[1]));
if (y < 0 || y >= mImage.GetDimension(1) - 1)
{
return Vector3<Real>::Zero();
return Vector3<T>::Zero();
}
int32_t z = static_cast<int32_t>(std::floor(position[2]));
if (z < 0 || z >= mImage.GetDimension(2) - 1)
{
return Vector3<Real>::Zero();
return Vector3<T>::Zero();
}
position[0] -= static_cast<Real>(x);
position[1] -= static_cast<Real>(y);
position[2] -= static_cast<Real>(z);
Real oneMX = (Real)1 - position[0];
Real oneMY = (Real)1 - position[1];
Real oneMZ = (Real)1 - position[2];
T const one = static_cast<T>(1);
position[0] -= static_cast<T>(x);
position[1] -= static_cast<T>(y);
position[2] -= static_cast<T>(z);
T oneMX = one - position[0];
T oneMY = one - position[1];
T oneMZ = one - position[2];
// Get image values at corners of voxel.
std::array<size_t, 8> corners{};
mImage.GetCorners(x, y, z, corners);
Real f000 = mImage[corners[0]];
Real f100 = mImage[corners[1]];
Real f010 = mImage[corners[2]];
Real f110 = mImage[corners[3]];
Real f001 = mImage[corners[4]];
Real f101 = mImage[corners[5]];
Real f011 = mImage[corners[6]];
Real f111 = mImage[corners[7]];
T f000 = mImage[corners[0]];
T f100 = mImage[corners[1]];
T f010 = mImage[corners[2]];
T f110 = mImage[corners[3]];
T f001 = mImage[corners[4]];
T f101 = mImage[corners[5]];
T f011 = mImage[corners[6]];
T f111 = mImage[corners[7]];
Vector3<Real> gradient;
Vector3<T> gradient{};
Real tmp0 = oneMY * (f100 - f000) + position[1] * (f110 - f010);
Real tmp1 = oneMY * (f101 - f001) + position[1] * (f111 - f011);
T tmp0 = oneMY * (f100 - f000) + position[1] * (f110 - f010);
T tmp1 = oneMY * (f101 - f001) + position[1] * (f111 - f011);
gradient[0] = oneMZ * tmp0 + position[2] * tmp1;
tmp0 = oneMX * (f010 - f000) + position[0] * (f110 - f100);
@@ -335,6 +395,6 @@ namespace gte
return gradient;
}
Image3<Real> const& mImage;
Image3<T> const& mImage;
};
}