levelset: cut cleanly where a 3D surface meets the box

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "openscad_cpp_evaluator"
-version = "0.51.0"
+version = "0.51.1"
 description = "C++ OpenSCAD evaluator with Python bindings"
 readme = "README.md"
 requires-python = ">=3.12"

src/builtins/topology.cpp

@@ -970,6 +970,19 @@ std::vector<ColoredBody> generateLevelSet(Evaluator& ev, const CSGParams& params
             if (i0[a] > n[a] - 2) i0[a] = n[a] - 2;
             t[a] = g[a] - static_cast<double>(i0[a]);
         }
+        // Outside the sampled block the answer is "outside", falling away
+        // with distance. Clamping instead (which is what the interpolation
+        // above does on its own) would smear the boundary values outward and
+        // extend the surface rather than close it.
+        double outside = 0.0;
+        for (int a = 0; a < 3; ++a) {
+            const double maxg = static_cast<double>(n[a] - 1);
+            const double raw = (pos[a] - origin[a]) / spacing[a];
+            if (raw < 0.0) outside = std::max(outside, -raw);
+            else if (raw > maxg) outside = std::max(outside, raw - maxg);
+        }
+        if (outside > 0.0) return invert ? (1.0 + outside) : -(1.0 + outside);
+
         double acc = 0.0;
         for (int c = 0; c < 8; ++c) {
             const size_t di = static_cast<size_t>(c & 1), dj = static_cast<size_t>((c >> 1) & 1),
@@ -980,8 +993,18 @@ std::vector<ColoredBody> generateLevelSet(Evaluator& ev, const CSGParams& params
         return bandDistance(acc, isoLo, isoHi, invert);
     };
 
-    manifold::Box bounds(manifold::vec3(origin[0], origin[1], origin[2]),
-                          manifold::vec3((*hi)[0], (*hi)[1], (*hi)[2]));
+    // Manifold closes the mesh where the surface meets the edge of the box it
+    // is given, and that closure follows the sample lattice -- a visible
+    // staircase on anything that reaches the boundary, which a gyroid does
+    // everywhere. So sample a PADDED box, let the ragged closure happen out
+    // there, and cut back to the box actually asked for. The cut is then a
+    // clean plane.
+    //
+    // Exactly the fix already applied to the 2D path (clipToBounds); it was
+    // not applied here at the time because the 2D error was the one measured.
+    const double pad = 2.0 * edge;
+    manifold::Box bounds(manifold::vec3(origin[0] - pad, origin[1] - pad, origin[2] - pad),
+                          manifold::vec3((*hi)[0] + pad, (*hi)[1] + pad, (*hi)[2] + pad));
     // tolerance -1: a positive value makes Manifold do EXTRA evaluations per
     // output vertex to snap nearer the true surface. Against a fixed grid
     // those only re-interpolate data already used -- cost, no information.
@@ -996,6 +1019,10 @@ std::vector<ColoredBody> generateLevelSet(Evaluator& ev, const CSGParams& params
                  : manifold::Manifold::LevelSet(sampleGrid, bounds, edge, 0.0, -1.0, /*canParallel=*/true);
 
     if (solid.IsEmpty()) return {};
+    solid = solid ^ manifold::Manifold::Cube(
+                        manifold::vec3((*hi)[0] - origin[0], (*hi)[1] - origin[1], (*hi)[2] - origin[2]))
+                        .Translate(manifold::vec3(origin[0], origin[1], origin[2]));
+    if (solid.IsEmpty()) return {};
     ColoredBody b;
     b.body = std::move(solid);
     std::vector<ColoredBody> out;

tests/test_booleans.cpp

@@ -1143,3 +1143,49 @@ TEST(LevelSet2d, A3dArrayWithA2dBoundsWarns) {
     EXPECT_FALSE(levelsetWarnings("levelset([[[1,2],[3,4]],[[5,6],[7,8]]], bounds=[[0,0],[1,1]]);")
                      .empty());
 }
+
+TEST(LevelSet, ASurfaceMeetingTheBoxIsCutCleanly) {
+    // Manifold closes the mesh where the surface reaches the edge of the box
+    // it is given, and that closure follows the SAMPLE LATTICE -- a visible
+    // staircase on anything that touches the boundary. A gyroid touches it
+    // everywhere, which is how this was found: the shape looked ragged next
+    // to BOSL2's isosurface() of the same field.
+    //
+    // Sampling a padded box and cutting back makes the cut a plane. The
+    // numeric form of that: a half-space is EXACTLY half the box, at every
+    // resolution. Before the fix it was resolution-dependent.
+    //
+    // The 2D path had this fixed already (clipToBounds); 3D did not, because
+    // 2D was the case that got measured at the time.
+    for (const char* edge : {"2", "1", "0.5"}) {
+        Evaluated e = evalSrc(std::string("levelset(function(x,y,z) x, "
+                                           "bounds=[[-20,-20,-20],[20,20,20]], "
+                                           "isovalue=[-1e18,0], edge=") +
+                               edge + ");");
+        ASSERT_EQ(e.bodies.size(), 1u) << "edge=" << edge;
+        EXPECT_NEAR(soleBody(e).Volume(), 32000.0, 1e-6) << "edge=" << edge;
+    }
+}
+
+TEST(LevelSet, PaddingDoesNotDisturbAShapeInsideTheBox) {
+    // The padding is sampled outside the requested bounds, so a shape that
+    // never reaches the boundary must come out exactly as before.
+    Evaluated e = evalSrc(sphereFieldSrc(50, 30) +
+                           "levelset(f, bounds=[[-30,-30,-30],[30,30,30]], isovalue=20);");
+    ASSERT_EQ(e.bodies.size(), 1u);
+    const double analytic = 4.0 / 3.0 * 3.14159265358979 * 8000;
+    EXPECT_NEAR(soleBody(e).Volume(), analytic, 0.01 * analytic);
+    EXPECT_EQ(soleBody(e).Genus(), 0);
+}
+
+TEST(LevelSet, AGridFieldIsAlsoCutCleanlyAtTheBox) {
+    // The grid intake cannot sample beyond its data, so the sampler reads
+    // "outside" past the block rather than clamping and smearing the edge
+    // values outward. Same clean cut as the function form.
+    Evaluated e = evalSrc(
+        "N = 41;\nfunction co(t) = -20 + 40*t/(N-1);\n"
+        "f = [for (i=[0:N-1]) [for (j=[0:N-1]) [for (k=[0:N-1]) co(i) ]]];\n"
+        "levelset(f, bounds=[[-20,-20,-20],[20,20,20]], isovalue=[-1e18,0]);");
+    ASSERT_EQ(e.bodies.size(), 1u);
+    EXPECT_NEAR(soleBody(e).Volume(), 32000.0, 40.0);   // one sample layer of slack
+}