Skip to content

roof(): voronoi method roofs unions of touching pieces #268

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
revarbat merged 1 commit into from
Oct 8, 2026
Merged

roof(): voronoi method roofs unions of touching pieces #268

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 81 additions & 20 deletions src/builtins/roof.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool> used(n, false);
for (const auto& t : triangles_) {
for (int i : t) mesh_.triVerts.push_back(static_cast<uint64_t>(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<uint64_t> to(n);
std::vector<double> 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<uint64_t>(i));
}
manifold::Manifold body(mesh_);
if (body.Status() == manifold::Manifold::Error::NoError) return body;
Expand Down Expand Up @@ -991,16 +1016,17 @@ void dropSpikes(std::vector<Vec2>& 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<std::vector<Vec2>> 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<std::vector<Vec2>> 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<std::vector<Vec2>> contours;
for (const auto& poly : polys) {
Expand All @@ -1016,17 +1042,6 @@ std::vector<std::vector<Vec2>> 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<std::vector<Vec2>> 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
Expand Down Expand Up @@ -1062,8 +1077,11 @@ std::vector<std::vector<Vec2>> weldContours(std::vector<std::vector<Vec2>> 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;
Expand All @@ -1074,8 +1092,10 @@ std::vector<std::vector<Vec2>> weldContours(std::vector<std::vector<Vec2>> conto
if (len == 0) continue;
const Vec2 lo = la::min(a, b) - kWeld, hi = la::max(a, b) + kWeld;
std::vector<std::pair<double, Vec2>> 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<double>::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});
}
Expand All @@ -1095,6 +1115,47 @@ std::vector<std::vector<Vec2>> weldContours(std::vector<std::vector<Vec2>> 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<std::vector<Vec2>> canonical(std::vector<std::vector<Vec2>> 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<std::vector<Vec2>> contours = gridContours(cs, toGrid, bits);
std::vector<std::vector<Vec2>> 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<std::vector<Vec2>> contours = weldContours(gridContours(cs, toGrid));
Expand Down
81 changes: 79 additions & 2 deletions tests/test_extrude_roof.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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; });
Expand Down
Loading