diff --git a/src/sdf.cpp b/src/sdf.cpp index 0e644fb19..51d880580 100644 --- a/src/sdf.cpp +++ b/src/sdf.cpp @@ -117,13 +117,15 @@ vec3 Position(ivec4 gridIndex, vec3 origin, vec3 spacing) { return origin + spacing * (vec3(gridIndex) + (gridIndex.w == 1 ? 0.0 : -0.5)); } +vec3 Bound(vec3 pos, vec3 origin, vec3 spacing, ivec3 gridSize) { + return min(max(pos, origin), origin + spacing * (vec3(gridSize) - 1)); +} + double BoundedSDF(ivec4 gridIndex, vec3 origin, vec3 spacing, ivec3 gridSize, double level, std::function sdf) { - auto Min = [](ivec3 p) { return std::min(p.x, std::min(p.y, p.z)); }; - const ivec3 xyz(gridIndex); - const int lowerBoundDist = Min(xyz); - const int upperBoundDist = Min(gridSize - xyz); + const int lowerBoundDist = minelem(xyz); + const int upperBoundDist = minelem(gridSize - xyz); const int boundDist = std::min(lowerBoundDist, upperBoundDist - gridIndex.w); if (boundDist < 0) { @@ -264,7 +266,7 @@ struct NearSurface { // Bound the delta of each vert to ensure the tetrahedron cannot invert. if (la::all(la::less(la::abs(pos - gridPos), kS * spacing))) { const int idx = AtomicAdd(vertIndex[0], 1); - vertPos[idx] = pos; + vertPos[idx] = Bound(pos, origin, spacing, gridSize); gridVert.movedVert = idx; for (int j = 0; j < 7; ++j) { if (gridVert.edgeVerts[j] == kCrossing) gridVert.edgeVerts[j] = idx; @@ -323,9 +325,10 @@ struct ComputeVerts { } const int idx = AtomicAdd(vertIndex[0], 1); - vertPos[idx] = FindSurface(position, gridVert.distance, - Position(neighborIndex, origin, spacing), val, - tol, level, sdf); + const vec3 pos = FindSurface(position, gridVert.distance, + Position(neighborIndex, origin, spacing), + val, tol, level, sdf); + vertPos[idx] = Bound(pos, origin, spacing, gridSize); gridVert.edgeVerts[i] = idx; } } @@ -423,23 +426,23 @@ namespace manifold { /** * Constructs a level-set manifold from the input Signed-Distance Function - * (SDF). This uses a form of Marching Tetrahedra (akin to Marching Cubes, but - * better for manifoldness). Instead of using a cubic grid, it uses a - * body-centered cubic grid (two shifted cubic grids). This means if your - * function's interior exceeds the given bounds, you will see a kind of - * egg-crate shape closing off the manifold, which is due to the underlying - * grid. + * (SDF). This uses a form of Marching Tetrahedra (akin to Marching + * Cubes, but better for manifoldness). Instead of using a cubic grid, it uses a + * body-centered cubic grid (two shifted cubic grids). These grid points are + * snapped to the surface where possible to keep short edges from forming. * * @param sdf The signed-distance functor, containing this function signature: * `double operator()(vec3 point)`, which returns the * signed distance of a given point in R^3. Positive values are inside, - * negative outside. + * negative outside. There is no requirement that the function be a true + * distance, or even continuous. * @param bounds An axis-aligned box that defines the extent of the grid. * @param edgeLength Approximate maximum edge length of the triangles in the * final result. This affects grid spacing, and hence has a strong effect on * performance. - * @param level You can inset your Mesh by using a positive value, or outset - * it with a negative value. + * @param level Extract the surface at this value of your sdf; defaults to + * zero. You can inset your mesh by using a positive value, or outset it with a + * negative value. * @param tolerance Ensure each vertex is within this distance of the true * surface. Defaults to -1, which will return the interpolated * crossing-point based on the two nearest grid points. Small positive values diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp index b6a27ccda..bac87ab45 100644 --- a/test/sdf_test.cpp +++ b/test/sdf_test.cpp @@ -78,7 +78,7 @@ TEST(SDF, Bounds) { EXPECT_EQ(cubeVoid.Status(), Manifold::Error::NoError); EXPECT_EQ(cubeVoid.Genus(), -1); - const double outerBound = size / 2 + edgeLength / 2; + const double outerBound = size / 2; EXPECT_NEAR(bounds.min.x, -outerBound, epsilon); EXPECT_NEAR(bounds.min.y, -outerBound, epsilon); EXPECT_NEAR(bounds.min.z, -outerBound, epsilon); @@ -102,7 +102,7 @@ TEST(SDF, Bounds2) { EXPECT_EQ(cubeVoid.Status(), Manifold::Error::NoError); EXPECT_EQ(cubeVoid.Genus(), -1); - const double outerBound = size / 2 + edgeLength / 2; + const double outerBound = size / 2; EXPECT_NEAR(bounds.min.x, -outerBound, epsilon); EXPECT_NEAR(bounds.min.y, -outerBound, epsilon); EXPECT_NEAR(bounds.min.z, -outerBound, epsilon); @@ -111,7 +111,28 @@ TEST(SDF, Bounds2) { EXPECT_NEAR(bounds.max.z, outerBound, epsilon); } -TEST(SDF, Surface) { +TEST(SDF, Bounds3) { + const double radius = 1.2; + Manifold sphere = + Manifold::LevelSet([radius](vec3 pos) { return radius - length(pos); }, + {vec3(-1), vec3(1)}, 0.1); +#ifdef MANIFOLD_EXPORT + if (options.exportModels) ExportMesh("sphere.glb", sphere.GetMeshGL(), {}); +#endif + + EXPECT_EQ(sphere.Status(), Manifold::Error::NoError); + EXPECT_EQ(sphere.Genus(), 0); + const double epsilon = sphere.GetEpsilon(); + Box bounds = sphere.BoundingBox(); + EXPECT_NEAR(bounds.min.x, -1, epsilon); + EXPECT_NEAR(bounds.min.y, -1, epsilon); + EXPECT_NEAR(bounds.min.z, -1, epsilon); + EXPECT_NEAR(bounds.max.x, 1, epsilon); + EXPECT_NEAR(bounds.max.y, 1, epsilon); + EXPECT_NEAR(bounds.max.z, 1, epsilon); +} + +TEST(SDF, Void) { const double size = 4; const double edgeLength = 0.5; @@ -142,20 +163,31 @@ TEST(SDF, Resize) { const double size = 20; Manifold layers = Manifold::LevelSet(Layers(), {vec3(0.0), vec3(size)}, 1); #ifdef MANIFOLD_EXPORT - if (options.exportModels) ExportMesh("layers.gltf", layers.GetMeshGL(), {}); + if (options.exportModels) ExportMesh("layers.glb", layers.GetMeshGL(), {}); #endif EXPECT_EQ(layers.Status(), Manifold::Error::NoError); EXPECT_EQ(layers.Genus(), -8); + const double outerBound = size / 2; + const double epsilon = layers.GetEpsilon(); + Box bounds = layers.BoundingBox(); + EXPECT_NEAR(bounds.min.x, 0, epsilon); + EXPECT_NEAR(bounds.min.y, 0, epsilon); + EXPECT_NEAR(bounds.min.z, 1.5, epsilon); + EXPECT_NEAR(bounds.max.x, size, epsilon); + EXPECT_NEAR(bounds.max.y, size, epsilon); + EXPECT_NEAR(bounds.max.z, size - 1.5, epsilon); } TEST(SDF, SineSurface) { - Manifold surface = Manifold::LevelSet( - [](vec3 p) { - double mid = la::sin(p.x) + la::sin(p.y); - return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0f : -1.0f; - }, - {vec3(-1.75 * kPi), vec3(1.75 * kPi)}, 1); + Manifold surface = + Manifold::LevelSet( + [](vec3 p) { + double mid = la::sin(p.x) + la::sin(p.y); + return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0f : -1.0f; + }, + {vec3(-1.75 * kPi), vec3(1.75 * kPi)}, 1) + .AsOriginal(); Manifold smoothed = surface.SmoothOut(180).RefineToLength(0.05); EXPECT_EQ(smoothed.Status(), Manifold::Error::NoError);