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()); +}