levelset(): a solid from a grid or a function(x,y,z)

CLAUDE.md

@@ -460,6 +460,43 @@ grep for `ponytail:`.
   unknown arguments -- `children(separate=true)` renders the wrong shape there with no warning --
   and both read as `undef` in OpenSCAD, so a guard written against them works in both.
 
+- **`levelset(field, bounds, isovalue)`** (`src/builtins/topology.cpp`, feature name `levelset`) —
+  a solid from a grid of sampled values: `field[i][j][k]` with i,j,k mapping to x,y,z. Wraps
+  `Manifold::LevelSet` with a C++ trilinear-sampling lambda; no new Manifold API.
+  **The field is a grid OR a `function(x,y,z)`, and the two trade off against each other** —
+  neither is simply better:
+
+  | | sampling | parallel | accuracy | ceiling |
+  |---|---|---|---|---|
+  | grid `field[i][j][k]` | 0.41 µs | **yes** | capped by the grid | ~550 MB at 200³ |
+  | `function(x,y,z)` | 0.96 µs | **no** | Manifold snaps toward the true surface | none |
+
+  Measured at matched resolution against the analytic sphere: **50³ grid 0.23 s (−0.244%) vs
+  function 0.32 s (−0.103%); 100³ grid 1.06 s (−0.059%) vs function 1.91 s (−0.025%)**. So the
+  function form is ~1.5× slower and ~2.4× more accurate. Both still beat BOSL2.
+  `canParallel` is the crux and is **not a tuning knob**: true for a grid because the lambda is
+  pure C++, false for a function because Manifold warns that parallel policies "will crash
+  language runtimes with runtime locks that expect to not be called back by unregistered threads".
+  The function form also needs `edge=` explicitly — there is no grid to infer spacing from and the
+  cost is cubic in it, so guessing would silently produce either a useless mesh or a ten-minute
+  one.
+  Measured against BOSL2's `isosurface.scad` (217 KB of marching cubes in script), same field:
+  **50³ 0.26 s vs 2.07 s (8.0×), 100³ 1.15 s vs 14.59 s (12.7×)**.
+  The trade is **memory**: ~67 bytes/cell, so 200³ ≈ 550 MB is the practical ceiling. A callback
+  samples lazily and has none.
+  **Accuracy is bounded by the grid, not by Manifold.** It samples on a body-centred cubic lattice
+  (`sdf.cpp:440-446`) so a plain cubic grid does not line up and the lambda interpolates — same
+  accuracy as grid-based marching cubes in script, worse than `LevelSet` given a true SDF. Hence
+  `tolerance = -1` is forced: a positive value makes Manifold do extra evaluations per vertex that
+  would only re-interpolate data already used. `edgeLength` defaults to the grid spacing for the
+  same reason.
+  Calling the field function from the generate pass has no live `EvalContext`, so a root is built
+  from the closure's own scope — `$`-variables sit at their defaults inside it.
+  Sign convention: Manifold takes positive as inside; a distance field is the other way round, so
+  the default flips it and `invert=true` flips back.
+  Two behaviours that look alike and are opposite: an isovalue **nothing reaches** gives empty
+  geometry; one **everything satisfies** gives the whole bounding box. Both are tested.
+
 - **`linear_solve(A, b)`** (`src/builtins/function_builtins.cpp`, feature name `linear-solve`) —
   returns `.x` (solution), `.det`, `.singular`. Dispatches on shape:
 

pyproject.toml

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

src/builtins/builtins.hpp

@@ -183,6 +183,11 @@ std::vector<ColoredBody> generateHull(Evaluator& ev, const CSGParams& params,
 Value builtinDxfDim(Evaluator& ev, const CallArgs& args, const oscad::ASTNode& node);
 Value builtinDxfCross(Evaluator& ev, const CallArgs& args, const oscad::ASTNode& node);
 
+CSGParams resolveLevelSet(Evaluator& ev, const oscad::ModularCall& node, EvalContext& ctx);
+std::vector<ColoredBody> generateLevelSet(Evaluator& ev, const CSGParams& params,
+                                           const std::vector<std::unique_ptr<CSGNode>>& children,
+                                           const oscad::ASTNode& node);
+
 CSGParams resolveSimplify(Evaluator& ev, const oscad::ModularCall& node, EvalContext& ctx);
 std::vector<ColoredBody> generateSimplify(Evaluator& ev, const CSGParams& params,
                                            const std::vector<std::unique_ptr<CSGNode>>& children,

src/builtins/function_builtins.cpp

@@ -913,6 +913,7 @@ const std::unordered_map<std::string, double>& featureLevels() {
     static const std::unordered_map<std::string, double> levels = {
         {"render-expr", 1.0},        // render() in expression position
         {"linear-solve", 1.0},       // linear_solve(A, b) -> {x, det, singular}
+        {"levelset", 1.0},           // levelset(field, bounds, isovalue)
         {"polyhedron-vnf", 1.0},     // polyhedron(vnf) / polyhedron(object)
         {"separate-children", 1.0},  // children(..., separate=true)
         {"minkowski-diff", 1.0},     // minkowski_difference()

src/builtins/registry.cpp

@@ -37,6 +37,7 @@ const std::unordered_map<std::string_view, ResolveFn>& resolveDispatch() {
         {"fill", &resolveFill},
         {"minkowski_difference", &resolveMinkowskiDifference},
         {"simplify", &resolveSimplify},
+        {"levelset", &resolveLevelSet},
         {"offset", &resolveOffset},
         {"surface", &resolveSurface},
         {"text", &resolveText},
@@ -78,6 +79,7 @@ const std::unordered_map<std::string_view, GenerateFn>& generateDispatch() {
         {"fill", &generateFill},
         {"minkowski_difference", &generateMinkowskiDifference},
         {"simplify", &generateSimplify},
+        {"levelset", &generateLevelSet},
         {"offset", &generateOffset},
         {"surface", &generateSurface},
         {"text", &generateText},
@@ -125,6 +127,7 @@ const std::vector<std::string>* builtinParamNames(const std::string& name) {
         {"fill", {}},
         {"minkowski_difference", {}},
         {"simplify", {"tolerance"}},
+        {"levelset", {"field", "bounds", "isovalue", "invert", "edge"}},
         {"offset", {"r", "delta", "chamfer"}},
         {"surface", {"file", "center", "convexity", "invert"}},
         {"text", {"text", "size", "font", "direction", "language", "script", "halign", "valign", "spacing"}},

src/builtins/topology.cpp

@@ -451,4 +451,253 @@ std::vector<ColoredBody> generateSimplify(Evaluator& ev, const CSGParams& params
     return out;
 }
 
+// -- levelset(field, bounds, isovalue) -------------------------------------
+//
+// A solid from a grid of sampled values. Wraps Manifold::LevelSet with a C++
+// trilinear-sampling lambda over the array the script handed us.
+//
+// Why a grid rather than Manifold's own SDF-callback shape, which would be
+// the obvious mapping: measured, building the field in OpenSCAD costs
+// 0.41 us/sample against 0.96 us for a closure call per sample, because the
+// arithmetic inlines into the list comprehension. More importantly a grid
+// means the C++ side never re-enters the evaluator, so canParallel can be
+// TRUE -- Manifold's docs warn that parallel policies "will crash language
+// runtimes with runtime locks that expect to not be called back by
+// unregistered threads", which is exactly what a callback design forfeits.
+//
+// The trade is memory: ~67 bytes per cell, so 200^3 is about 550 MB and that
+// is the practical ceiling. A callback samples lazily and has no ceiling.
+//
+// Accuracy is bounded by the GRID, not by Manifold. Manifold samples on a
+// body-centred cubic lattice (two interleaved cubic grids, sdf.cpp:440-446)
+// and a plain cubic grid does not line up with those points, so the lambda
+// interpolates. Same accuracy as a grid-based marching-cubes implementation
+// in script; worse than LevelSet given a true SDF. Do not imply otherwise.
+
+namespace {
+
+struct ScalarField {
+    std::vector<double> v;      // flat, x fastest
+    size_t nx = 0, ny = 0, nz = 0;
+    double at(size_t i, size_t j, size_t k) const { return v[(k * ny + j) * nx + i]; }
+};
+
+// field[i][j][k] -> flat, checking it is a full rectangular block.
+std::optional<ScalarField> readField(const Value& val) {
+    const ListPtr* xs = std::get_if<ListPtr>(&val);
+    if (!xs || !*xs || (*xs)->items.empty()) return std::nullopt;
+    ScalarField f;
+    f.nx = (*xs)->items.size();
+    for (size_t i = 0; i < f.nx; ++i) {
+        const ListPtr* ys = std::get_if<ListPtr>(&(*xs)->items[i]);
+        if (!ys || !*ys || (*ys)->items.empty()) return std::nullopt;
+        if (i == 0) f.ny = (*ys)->items.size();
+        else if ((*ys)->items.size() != f.ny) return std::nullopt;
+        for (size_t j = 0; j < f.ny; ++j) {
+            const ListPtr* zs = std::get_if<ListPtr>(&(*ys)->items[j]);
+            if (!zs || !*zs || (*zs)->items.empty()) return std::nullopt;
+            if (i == 0 && j == 0) {
+                f.nz = (*zs)->items.size();
+                f.v.assign(f.nx * f.ny * f.nz, 0.0);
+            } else if ((*zs)->items.size() != f.nz) {
+                return std::nullopt;
+            }
+            for (size_t k = 0; k < f.nz; ++k) {
+                const double* d = std::get_if<double>(&(*zs)->items[k]);
+                if (!d || !std::isfinite(*d)) return std::nullopt;
+                f.v[(k * f.ny + j) * f.nx + i] = *d;
+            }
+        }
+    }
+    return f;
+}
+
+std::optional<std::array<double, 3>> readVec3(const Value& v) {
+    const ListPtr* l = std::get_if<ListPtr>(&v);
+    if (!l || !*l || (*l)->items.size() != 3) return std::nullopt;
+    std::array<double, 3> out{};
+    for (size_t i = 0; i < 3; ++i) {
+        const double* d = std::get_if<double>(&(*l)->items[i]);
+        if (!d || !std::isfinite(*d)) return std::nullopt;
+        out[i] = *d;
+    }
+    return out;
+}
+
+} // namespace
+
+CSGParams resolveLevelSet(Evaluator& ev, const oscad::ModularCall& node, EvalContext& ctx) {
+    auto [args, effCtx] = resolveCallArgs(ev, node.arguments, ctx);
+    CSGParams params;
+    params["field"] = getArg(args, 0, "field", Value{});
+    params["bounds"] = getArg(args, 1, "bounds", Value{});
+    params["isovalue"] = getArg(args, 2, "isovalue", Value{0.0});
+    params["invert"] = getArg(args, 3, "invert", Value{false});
+    params["edge"] = getArg(args, 4, "edge", Value{});
+    ev.evalChildren(node.children, effCtx);
+    return params;
+}
+
+std::vector<ColoredBody> generateLevelSet(Evaluator& ev, const CSGParams& params,
+                                           const std::vector<std::unique_ptr<CSGNode>>&,
+                                           const oscad::ASTNode& node) {
+    // Two shapes of field, and the choice matters beyond taste:
+    //
+    //   a GRID    -- faster to sample (0.41 vs 0.96 us) and meshes in
+    //                PARALLEL, because nothing re-enters the evaluator.
+    //                Accuracy is capped by the grid; memory is the ceiling.
+    //   a FUNCTION -- Manifold picks its own sample points and can snap to
+    //                the true surface, so it is MORE accurate, with no
+    //                memory ceiling. But every sample is a closure call, and
+    //                canParallel must be false or the evaluator gets called
+    //                back from threads it knows nothing about.
+    const Value& fieldArg = params.at("field");
+    const ClosurePtr* fieldFn = std::get_if<ClosurePtr>(&fieldArg);
+    std::optional<ScalarField> field;
+    if (!fieldFn) {
+        field = readField(fieldArg);
+        if (!field) {
+            ev.warn("levelset(): field must be a function(x,y,z) or a rectangular field[i][j][k] of numbers",
+                    &node.position());
+            return {};
+        }
+        if (field->nx < 2 || field->ny < 2 || field->nz < 2) {
+            ev.warn("levelset(): the field needs at least 2 samples along each axis", &node.position());
+            return {};
+        }
+    } else if (!*fieldFn || (*fieldFn)->node == nullptr || (*fieldFn)->node->parameters.size() < 3) {
+        ev.warn("levelset(): the field function needs three parameters, as in function(x,y,z) ...",
+                &node.position());
+        return {};
+    }
+
+    const ListPtr* bb = std::get_if<ListPtr>(&params.at("bounds"));
+    std::optional<std::array<double, 3>> lo, hi;
+    if (bb && *bb && (*bb)->items.size() == 2) {
+        lo = readVec3((*bb)->items[0]);
+        hi = readVec3((*bb)->items[1]);
+    }
+    if (!lo || !hi) {
+        ev.warn("levelset(): bounds must be [[x0,y0,z0],[x1,y1,z1]]", &node.position());
+        return {};
+    }
+    for (int a = 0; a < 3; ++a) {
+        if (!((*hi)[a] > (*lo)[a])) {
+            ev.warn("levelset(): bounds must be increasing along every axis", &node.position());
+            return {};
+        }
+    }
+
+    const double* isoArg = std::get_if<double>(&params.at("isovalue"));
+    if (!isoArg) {
+        ev.warn("levelset(): isovalue must be a number", &node.position());
+        return {};
+    }
+    const double isovalue = *isoArg;
+    const bool invert = truthy(params.at("invert"));
+
+    const std::array<double, 3> origin = *lo;
+    std::array<double, 3> spacing{1.0, 1.0, 1.0};
+    double edge = 0.0;
+    if (field) {
+        // Grid spacing is corner-sample to corner-sample, so n-1 intervals.
+        spacing = {((*hi)[0] - (*lo)[0]) / static_cast<double>(field->nx - 1),
+                    ((*hi)[1] - (*lo)[1]) / static_cast<double>(field->ny - 1),
+                    ((*hi)[2] - (*lo)[2]) / static_cast<double>(field->nz - 1)};
+        // Finer than the grid buys nothing -- there is no information between
+        // samples -- and coarser throws away what the script paid to compute.
+        edge = std::min({spacing[0], spacing[1], spacing[2]});
+    }
+    if (const double* e = std::get_if<double>(&params.at("edge"))) {
+        if (*e > 0.0) edge = *e;
+        else ev.warn("levelset(): edge must be positive", &node.position());
+    }
+    if (edge <= 0.0) {
+        // No grid to infer it from. Guessing would silently pick either a
+        // useless mesh or a ten-minute one -- the cost is cubic in this
+        // number, so it is the caller's call to make.
+        ev.warn("levelset(): a function field needs edge= (the sample spacing)", &node.position());
+        return {};
+    }
+
+    // The function form: one closure call per sample. No live EvalContext at
+    // generate time, so a root is built from the closure's own scope -- the
+    // same approach the warp experiment used. $-variables are therefore at
+    // their defaults inside a field function.
+    std::optional<EvalContext> fnCtx;
+    std::array<std::string, 3> fnParams;
+    if (fieldFn) {
+        fnCtx = EvalContext::makeRoot((*fieldFn)->node->scope());
+        for (int a = 0; a < 3; ++a) fnParams[static_cast<size_t>(a)] = (*fieldFn)->node->parameters[static_cast<size_t>(a)]->name->name;
+    }
+
+    const auto sampleFn = [&](manifold::vec3 p) -> double {
+        BoundArgs bound;
+        const double xyz[3] = {p.x, p.y, p.z};
+        for (int a = 0; a < 3; ++a) bound.set(fnParams[static_cast<size_t>(a)], Value{xyz[a]});
+        const Value out = ev.evalFunctionLiteralFromBound(**fieldFn, std::move(bound), *fnCtx, &node.position());
+        const double* d = std::get_if<double>(&out);
+        // A field that returns junk at one point must not abort the whole
+        // build; treat it as far outside instead.
+        const double v = (d && std::isfinite(*d)) ? *d : std::numeric_limits<double>::max();
+        const double diff = v - isovalue;
+        return invert ? diff : -diff;
+    };
+
+    // Never dereferenced on the function path, but it must still be a valid
+    // reference -- binding one to a null pointer is undefined behaviour even
+    // where the lambda is never called.
+    const ScalarField emptyField;
+    const ScalarField& f = field ? *field : emptyField;
+    const auto sampleGrid = [&f, origin, spacing, isovalue, invert](manifold::vec3 p) -> double {
+        double g[3];
+        size_t i0[3];
+        double t[3];
+        const double pos[3] = {p.x, p.y, p.z};
+        const size_t n[3] = {f.nx, f.ny, f.nz};
+        for (int a = 0; a < 3; ++a) {
+            g[a] = (pos[a] - origin[a]) / spacing[a];
+            if (g[a] < 0.0) g[a] = 0.0;
+            const double maxg = static_cast<double>(n[a] - 1);
+            if (g[a] > maxg) g[a] = maxg;
+            i0[a] = static_cast<size_t>(g[a]);
+            if (i0[a] > n[a] - 2) i0[a] = n[a] - 2;
+            t[a] = g[a] - static_cast<double>(i0[a]);
+        }
+        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),
+                          dk = static_cast<size_t>((c >> 2) & 1);
+            const double w = (di ? t[0] : 1.0 - t[0]) * (dj ? t[1] : 1.0 - t[1]) * (dk ? t[2] : 1.0 - t[2]);
+            acc += w * f.at(i0[0] + di, i0[1] + dj, i0[2] + dk);
+        }
+        // Manifold takes POSITIVE as inside. A distance field is the other
+        // way round, which is the common case, so the default flips it.
+        const double d = acc - isovalue;
+        return invert ? d : -d;
+    };
+
+    manifold::Box bounds(manifold::vec3(origin[0], origin[1], origin[2]),
+                          manifold::vec3((*hi)[0], (*hi)[1], (*hi)[2]));
+    // 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.
+    // canParallel true: the lambda is pure C++ and never re-enters the
+    // evaluator, which is the whole point of taking a grid.
+    // canParallel is the crux: safe for a grid because the lambda is pure
+    // C++, and NOT safe for a function, which re-enters the evaluator.
+    // Manifold: parallel policies "will crash language runtimes with runtime
+    // locks that expect to not be called back by unregistered threads".
+    manifold::Manifold solid =
+        fieldFn ? manifold::Manifold::LevelSet(sampleFn, bounds, edge, 0.0, -1.0, /*canParallel=*/false)
+                 : manifold::Manifold::LevelSet(sampleGrid, bounds, edge, 0.0, -1.0, /*canParallel=*/true);
+
+    if (solid.IsEmpty()) return {};
+    ColoredBody b;
+    b.body = std::move(solid);
+    std::vector<ColoredBody> out;
+    out.push_back(std::move(b));
+    return out;
+}
+
 } // namespace oscadeval