SDF boundary cleanup (#1037)

* tighten test

* fixed SDF eggcrate

* addressing feedback

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<double(vec3)> 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

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);