From eec621859e07fe3577c34d8787b3110d06874d8d Mon Sep 17 00:00:00 2001 From: Revar Desmera Date: Wed, 7 Oct 2026 21:23:18 -0700 Subject: [PATCH] roof(): voronoi method roofs unions of touching pieces The default roof() gave nothing -- or ran away allocating -- for many unions of triangles/rectangles that touch along edges. Three causes: - Clipper returns pieces sharing an edge as outlines a hairline apart, which on the grid overlap or cross. Boost's Voronoi builder needs segments that meet only at their ends: the diagram came out wrong (no roof), or the builder never finished, growing without bound (16 of 6,000 corpus inputs; tens of GB). The outline is now welded and re-unioned first, as method="straight" already does, keeping the original outline (and so the identical roof) when welding changes nothing. - On a clean outline Boost still occasionally builds a wrong diagram where four sites are all but equidistant after rounding to the grid (a cell with no finite edges). The roof is retried on coarser grids (2^28..2^25) when the mesh does not come out a valid, non-empty solid. - An unused vertex -- the far end of a Voronoi edge outside the outline, up to 10^12 away -- stretched the bounding box Manifold sizes its tolerance by; its cleanup then collapsed the roof to nothing, or quietly dropped detail (old volumes off by up to 25% vs OpenSCAD). Such far vertices are now left out of the mesh. The weld's corner-vs-edge scan now only visits corners in each edge's x range, so the extra pass costs 2-10% of roof time at 1000-5000 corners. method="straight" output is byte-identical. Co-Authored-By: Claude Opus 5.5 --- src/builtins/roof.cpp | 101 +++++++++++++++++++++++++++++------- tests/test_extrude_roof.cpp | 81 ++++++++++++++++++++++++++++- 2 files changed, 160 insertions(+), 22 deletions(-) diff --git a/src/builtins/roof.cpp b/src/builtins/roof.cpp index a20c49b..669b425 100644 --- a/src/builtins/roof.cpp +++ b/src/builtins/roof.cpp @@ -114,10 +114,35 @@ class RoofMesh { } } + // A vertex no triangle uses -- the far end of an unbounded Voronoi edge + // outside the outline can lie 10^12 away -- would stretch the bounding + // box Manifold sizes its tolerance by, and its cleanup would then + // collapse the whole roof; such vertices are left out. (Only when one + // lies well beyond the rest: renumbering reorders the output, and one + // nearer changes the tolerance by no more than a few times.) manifold::Manifold build() { mesh_.numProp = 3; + const size_t n = mesh_.vertProperties.size() / 3; + std::vector used(n, false); for (const auto& t : triangles_) { - for (int i : t) mesh_.triVerts.push_back(static_cast(i)); + for (int i : t) used[i] = true; + } + manifold::Box box; + for (size_t i = 0; i < n; ++i) { + if (used[i]) box.Union(manifold::vec3(mesh_.vertProperties[3 * i], mesh_.vertProperties[3 * i + 1], mesh_.vertProperties[3 * i + 2])); + } + std::vector to(n); + std::vector kept; + bool stray = false; + for (size_t i = 0; i < n; ++i) { + const manifold::vec3 p(mesh_.vertProperties[3 * i], mesh_.vertProperties[3 * i + 1], mesh_.vertProperties[3 * i + 2]); + stray = stray || (!used[i] && la::maxelem(la::abs(p)) > 16 * box.Scale()); + to[i] = kept.size() / 3; + if (used[i]) kept.insert(kept.end(), {p.x, p.y, p.z}); + } + if (stray) mesh_.vertProperties = std::move(kept); + for (const auto& t : triangles_) { + for (int i : t) mesh_.triVerts.push_back(stray ? to[i] : static_cast(i)); } manifold::Manifold body(mesh_); if (body.Status() == manifold::Manifold::Error::NoError) return body; @@ -991,16 +1016,17 @@ void dropSpikes(std::vector& c) { } } -// The outline of `cs` snapped to a grid of 2^30 steps across its largest -// coordinate, with repeated and straight corners dropped. Sets `toGrid`. -std::vector> gridContours(const manifold::CrossSection& cs, double& toGrid) { +// The outline of `cs` snapped to a grid of 2^(bits+1) steps across its +// largest coordinate, with repeated and straight corners dropped. Sets +// `toGrid`. +std::vector> gridContours(const manifold::CrossSection& cs, double& toGrid, int bits = 29) { const manifold::Polygons polys = cs.ToPolygons(); double extent = 0; for (const auto& poly : polys) { for (Vec2 p : poly) extent = std::max({extent, std::abs(p.x), std::abs(p.y)}); } if (!(extent > 0) || !std::isfinite(extent)) return {}; - toGrid = std::ldexp(1.0, 29 - std::ilogb(extent)); + toGrid = std::ldexp(1.0, bits - std::ilogb(extent)); std::vector> contours; for (const auto& poly : polys) { @@ -1016,17 +1042,6 @@ std::vector> gridContours(const manifold::CrossSection& cs, do return contours; } -// The roof solid over `cs` (its floor at z = 0), discretizing curved parts -// by `fa`/`fs`. -manifold::Manifold voronoiRoof(const manifold::CrossSection& cs, double fa, double fs) { - double toGrid = 1; - std::vector> contours = gridContours(cs, toGrid); - if (contours.empty()) return manifold::Manifold(); - // $fa/$fs below 0.01 would ask for unbounded detail (computeRoofParams - // has already clamped them, warning). - return RoofBuilder(std::move(contours), 1 / toGrid, std::max(fa, 0.01), std::max(fs, 0.01)).build(); -} - // The outline welded where it touches itself: a union can leave pieces // that share an edge, or a hole against the outline, as outlines a hairline // apart, crossing at nothing Clipper can see. Each corner within kWeld of @@ -1062,8 +1077,11 @@ std::vector> weldContours(std::vector> conto } } } - // ponytail: every corner against every edge, O(n^2); a grid of edges - // when outlines of tens of thousands of corners matter. + // Each edge looks only at the corners in its x range. + // ponytail: still O(n^2) for an outline of long edges all spanning one + // x range; a grid of edges if that ever matters. + const auto byX = [](Vec2 p, Vec2 q) { return p.x != q.x ? p.x < q.x : p.y < q.y; }; + std::sort(corners.begin(), corners.end(), byX); Clipper2Lib::Paths64 paths; for (const auto& c : contours) { Clipper2Lib::Path64 path; @@ -1074,8 +1092,10 @@ std::vector> weldContours(std::vector> conto if (len == 0) continue; const Vec2 lo = la::min(a, b) - kWeld, hi = la::max(a, b) + kWeld; std::vector> on; - for (Vec2 q : corners) { - if (q.x < lo.x || q.y < lo.y || q.x > hi.x || q.y > hi.y || q == a || q == b) continue; + for (auto it = std::lower_bound(corners.begin(), corners.end(), Vec2{lo.x, -std::numeric_limits::infinity()}, byX); + it != corners.end() && it->x <= hi.x; ++it) { + const Vec2 q = *it; + if (q.y < lo.y || q.y > hi.y || q == a || q == b) continue; const double t = la::dot(q - a, d) / len; if (t > 0 && t < len && std::abs(cross2(d, q - a)) / len <= kWeld) on.push_back({t, q}); } @@ -1095,6 +1115,47 @@ std::vector> weldContours(std::vector> conto return welded; } +// Each contour rotated to start at its least corner, and the contours sorted: +// equal for two outlines that differ only in where Clipper starts them. +std::vector> canonical(std::vector> contours) { + const auto less = [](Vec2 a, Vec2 b) { return a.x != b.x ? a.x < b.x : a.y < b.y; }; + for (auto& c : contours) std::rotate(c.begin(), std::min_element(c.begin(), c.end(), less), c.end()); + std::sort(contours.begin(), contours.end(), [&](const auto& a, const auto& b) { return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end(), less); }); + return contours; +} + +// The roof solid over `cs` (its floor at z = 0), discretizing curved parts +// by `fa`/`fs`. +// +// Boost's Voronoi builder needs segments that meet only at their ends. A +// union's pieces sharing an edge come back as outlines a hairline apart, +// which on the grid overlap or cross: the diagram comes out wrong, or the +// builder never finishes, allocating without bound. So the outline is +// welded first, as for the straight skeleton. +// +// Even on valid input the builder occasionally gets a near-degenerate +// configuration wrong -- a cell left with no finite edges, where four sites +// are all but equidistant once the outline is rounded to the grid -- and the +// facets built from it do not close. Rounding to a coarser grid perturbs +// the outline differently, so the roof is retried there. +manifold::Manifold voronoiRoof(const manifold::CrossSection& cs, double fa, double fs) { + manifold::Manifold body; + for (int bits = 29; bits >= 25; --bits) { + double toGrid = 1; + std::vector> contours = gridContours(cs, toGrid, bits); + std::vector> welded = weldContours(contours); + // An outline welding leaves as it was keeps its own order, so its + // roof is the same, vertex for vertex, as without welding. + if (canonical(welded) != canonical(contours)) contours = std::move(welded); + if (contours.empty()) return manifold::Manifold(); + // $fa/$fs below 0.01 would ask for unbounded detail + // (computeRoofParams has already clamped them, warning). + body = RoofBuilder(std::move(contours), 1 / toGrid, std::max(fa, 0.01), std::max(fs, 0.01)).build(); + if (body.Status() == manifold::Manifold::Error::NoError && !body.IsEmpty()) break; + } + return body; +} + manifold::Manifold straightRoof(const manifold::CrossSection& cs) { double toGrid = 1; std::vector> contours = weldContours(gridContours(cs, toGrid)); diff --git a/tests/test_extrude_roof.cpp b/tests/test_extrude_roof.cpp index 555fc04..a5e8f6a 100644 --- a/tests/test_extrude_roof.cpp +++ b/tests/test_extrude_roof.cpp @@ -347,8 +347,8 @@ struct RoofResult { double slopeError; // roof area / (floor area * sqrt 2) - 1: 0 if every face rises at 45 degrees }; -RoofResult straightRoof(const std::string& shape) { - Evaluated e = evalSrc("roof(method=\"straight\") " + shape); +RoofResult straightRoof(const std::string& shape, const std::string& method = "straight") { + Evaluated e = evalSrc("roof(method=\"" + method + "\") " + shape); EXPECT_EQ(e.bodies.size(), 1u); if (e.bodies.empty() || !e.bodies[0].body) return {0, 0, manifold::Manifold::Error::InvalidConstruction, 1}; const manifold::Manifold& m = *e.bodies[0].body; @@ -463,6 +463,83 @@ TEST(Roof, StraightHandlesThousandsOfVertices) { EXPECT_NEAR(g.slopeError, 0, 1e-5); } +// -- method="voronoi" over unions of pieces that touch along edges. These +// all gave no roof at all before: Clipper returns the pieces as outlines a +// hairline apart, which on the grid overlap or cross, and Boost's Voronoi +// builder needs segments that meet only at their ends. + +namespace { + +void expectVoronoiRoof(const std::string& shape, double volume, double height, double tol) { + SCOPED_TRACE(shape); + const RoofResult r = straightRoof(shape, "voronoi"); + EXPECT_EQ(r.status, manifold::Manifold::Error::NoError); + EXPECT_NEAR(r.volume, volume, tol * volume); + EXPECT_NEAR(r.height, height, tol); +} + +} // namespace + +// Four triangles sharing edges. The reference agrees: 23.67785 / 2.39894. +TEST(Roof, VoronoiJoinsTrianglesThatShareEdges) { + expectVoronoiRoof("union() { polygon([[2.764760357, 4.118581106], [5.901400476, 6.372929785], [8.606126287, 1.988496895]]); " + "polygon([[5.901400476, 6.372929785], [8.606126287, 1.988496895], [11.028736411, 6.50946765]]); " + "polygon([[8.606126287, 1.988496895], [2.764760357, 4.118581106], [4.975336337, 1.106200614]]); " + "polygon([[8.606126287, 1.988496895], [2.764760357, 4.118581106], [6.985603416, 6.618989998]]); }", + 23.677855, 2.39894, 1e-5); +} + +// Rotated rectangles sharing (parts of) edges. OpenSCAD 2026.02.01 draws +// nothing for it; the volume is the integral of the distance to the outline +// to within the faceting of the cones (170.68 numerically). +TEST(Roof, VoronoiJoinsRotatedRectanglesThatShareEdges) { + expectVoronoiRoof("union() { polygon([[-3.369828352, 0.413896268], [-7.46233132, 4.010047378], [-11.05848243, -0.08245559], [-6.965979462, -3.6786067]]); " + "polygon([[-4.919727978, -5.476682255], [-9.012230946, -1.880531145], [-14.406457611, -8.019285597], [-10.313954643, -11.615436707]]); " + "polygon([[6.613253139, -4.732154468], [0.474498687, 0.662072197], [-4.919727978, -5.476682255], [1.219026474, -10.87090892]]); " + "polygon([[10.457580178, -4.483978539], [8.411328694, -2.685902984], [6.613253139, -4.732154468], [8.659504623, -6.530230023]]); }", + 170.76468, 4.08601, 1e-5); +} + +// Here the outline as Clipper leaves it crosses itself on the grid, and the +// Voronoi builder never finished, allocating until the process was killed +// (tens of gigabytes). (The integral of the distance to the outline is +// 221.3; the cones over corners meeting head on are faceted coarsely, as by +// the reference.) +TEST(Roof, VoronoiRoofOfACrossingOutlineStaysBounded) { + expectVoronoiRoof("union() { polygon([[8.304902057, -8.407695464], [0.596294574, -9.683488249], [1.872087359, -17.392095732], [9.580694842, -16.116302947]]); " + "polygon([[14.737716755, 0.576704804], [7.029109272, -0.699087981], [8.304902057, -8.407695464], [16.01350954, -7.131902679]]); " + "polygon([[8.304902057, -8.407695464], [-3.258009167, -10.321384641], [-2.620112775, -14.175688383], [8.94279845, -12.261999205]]); " + "polygon([[2.536909139, 2.517319368], [-1.317394603, 1.879422976], [0.596294574, -9.683488249], [4.450598316, -9.045591856]]); }", + 222.74776, 3.90673, 1e-5); +} + +// A staircase of triangles. Outside it, two of its edges are all but +// parallel, so their Voronoi bisectors meet some 10^12 away; that vertex, +// used by no facet, stretched the bounding box Manifold sizes its tolerance +// by, and its cleanup collapsed the whole roof to nothing. (The integral of +// the distance to the outline is 78.02.) +TEST(Roof, VoronoiRoofIgnoresFarVoronoiVertices) { + expectVoronoiRoof("union() { polygon([[-1.015413575, -0.344245093], [-1.006343725, -5.505019653], [4.154430835, -5.495949802]]); " + "polygon([[-1.015413575, -0.344245093], [4.154430835, -5.495949802], [4.145360985, -0.335175242]]); " + "polygon([[-1.006343725, -5.505019653], [4.163500685, -10.656724363], [4.154430835, -5.495949802]]); " + "polygon([[4.154430835, -5.495949802], [4.163500685, -10.656724363], [9.324275246, -10.647654512]]); " + "polygon([[4.163500685, -10.656724363], [9.333345096, -15.808429072], [9.324275246, -10.647654512]]); " + "polygon([[-6.176188135, -0.353314943], [-1.006343725, -5.505019653], [-1.015413575, -0.344245093]]); }", + 78.047622, 3.02312, 1e-5); +} + +// Axis-aligned rectangles: a clean outline, but rounded to the grid four +// sites are all but equidistant at one point, and Boost's builder leaves a +// cell with no finite edges. A coarser grid rounds it differently. +// OpenSCAD 2026.02.01 draws nothing here either; the integral of the +// distance to the outline is 1.348 (the faceted cone sits a little above). +TEST(Roof, VoronoiRetriesWhereTheBuilderFails) { + expectVoronoiRoof("union() { polygon([[0.7, -0.7], [2.1, -0.7], [2.1, 0.0], [0.7, 0.0]]); " + "polygon([[2.1, -0.7], [2.8, -0.7], [2.8, 0.7], [2.1, 0.7]]); " + "polygon([[2.8, -1.4], [4.2, -1.4], [4.2, 0.7], [2.8, 0.7]]); }", + 1.3513906, 0.820101, 1e-5); +} + TEST(Roof, UnknownMethodWarnsAndFallsBackToVoronoi) { std::string lastWarning; Evaluated e = evalSrc("roof(method=\"bogus\") square(4, center=true);", [&](const std::string& msg) { lastWarning = msg; });