diff --git a/pyproject.toml b/pyproject.toml index d6692b7..59c9231 100644 --- a/pyproject.toml +++ b/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" diff --git a/src/builtins/topology.cpp b/src/builtins/topology.cpp index d620da7..f78f5b6 100644 --- a/src/builtins/topology.cpp +++ b/src/builtins/topology.cpp @@ -970,6 +970,19 @@ std::vector generateLevelSet(Evaluator& ev, const CSGParams& params if (i0[a] > n[a] - 2) i0[a] = n[a] - 2; t[a] = g[a] - static_cast(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(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(c & 1), dj = static_cast((c >> 1) & 1), @@ -980,8 +993,18 @@ std::vector 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 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 out; diff --git a/tests/test_booleans.cpp b/tests/test_booleans.cpp index c4bc9e6..86cd861 100644 --- a/tests/test_booleans.cpp +++ b/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 +}