From 21570d7f1e2eb2030428e0e6d3451ba5d9c8e42c Mon Sep 17 00:00:00 2001 From: Revar Desmera Date: Mon, 10 Aug 2026 18:50:11 -0700 Subject: [PATCH] Triangulate polyhedron faces properly instead of fanning them A fan -- (v0, vi, vi+1) for every i -- is only correct for a polygon that is both convex and planar, and BOSL2's vnf_polyhedron() routinely hands over neither. The end caps of a nurbs_sheet() are 34-gons that are concave and 3.5 units out of plane. Fanning one laid 32 triangles covering 281% of the cap's own area with 15 of them wound inside out, and the finished solid came out 11.8% larger in surface area and 5.5% larger in volume than the reference's from the same input. Ear clipping in the face's own best-fit plane now, with the normal from Newell's method so it stays meaningful when the points are not coplanar. Among the valid ears it takes the best-shaped one rather than the first: any ear gives a correct triangulation, but on a face that is not flat the choice decides how the surface folds, and first-found strung long thin triangles across the caps -- still 29% above the reference. Preferring fat ears tracks the surface, and the cap now comes out at 1102.48 square units against the reference's 1102.48. The whole model matches its surface area exactly and its volume to 0.02%. Degenerate cases fall back to a fan, which cannot make them worse: an all-collinear face has no plane to project onto, and a self-intersecting one has no valid ear to find. Co-Authored-By: Claude Opus 5 (1M context) --- pyproject.toml | 2 +- src/builtins/primitives_3d.cpp | 167 +++++++++++++++++++++++++++++++-- tests/test_booleans.cpp | 90 ++++++++++++++++++ 3 files changed, 250 insertions(+), 9 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 30c08f4..8f29b1e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build" [project] name = "openscad_cpp_evaluator" -version = "0.28.2" +version = "0.29.0" description = "C++ OpenSCAD evaluator with Python bindings" readme = "README.md" requires-python = ">=3.12" diff --git a/src/builtins/primitives_3d.cpp b/src/builtins/primitives_3d.cpp index fce148f..ff52623 100644 --- a/src/builtins/primitives_3d.cpp +++ b/src/builtins/primitives_3d.cpp @@ -239,6 +239,164 @@ std::vector generateCylinder(Evaluator& ev, const CSGParams& params // otherwise produce a NotManifold body). Mirrors // _resolve_polyhedron/_generate_polyhedron. + +namespace { + +// Triangulate one polyhedron face. +// +// A fan -- (v0, vi, vi+1) for every i -- is only correct for a polygon +// that is both convex and planar, and BOSL2's vnf_polyhedron() routinely +// hands over neither. The end caps of a nurbs_sheet() are 34-gons that +// are concave and 3.5 units out of plane; fanning one produced 32 +// triangles covering 281% of the cap's true area with 15 of them wound +// inside out, which inflated the finished solid by 12% of its surface +// area and 5% of its volume. +// +// Ear clipping in the face's own best-fit plane instead. Newell's method +// gives a normal that stays meaningful when the points are not coplanar, +// which is what makes projecting them usable at all here. +// +// ponytail: O(n^2). Faces are a handful of points in nearly every model +// and 34 in the one that prompted this; revisit if a model ever arrives +// with thousand-sided faces. +void triangulateFace(const std::vector>& verts, const std::vector& loop, + std::vector& out) { + const size_t n = loop.size(); + if (n < 3) return; + + // Winding is reversed on the way out throughout: OpenSCAD's faces are + // clockwise seen from outside, Manifold wants counter-clockwise. + auto emit = [&out](size_t a, size_t b, size_t c) { + if (a == b || b == c || a == c) return; + out.push_back(static_cast(a)); + out.push_back(static_cast(c)); + out.push_back(static_cast(b)); + }; + if (n == 3) { + emit(loop[0], loop[1], loop[2]); + return; + } + + std::array nrm = {0.0, 0.0, 0.0}; + for (size_t i = 0; i < n; ++i) { + const std::array& a = verts[loop[i]]; + const std::array& b = verts[loop[(i + 1) % n]]; + nrm[0] += (a[1] - b[1]) * (a[2] + b[2]); + nrm[1] += (a[2] - b[2]) * (a[0] + b[0]); + nrm[2] += (a[0] - b[0]) * (a[1] + b[1]); + } + const double len = std::sqrt(nrm[0] * nrm[0] + nrm[1] * nrm[1] + nrm[2] * nrm[2]); + if (!(len > 1e-12)) { + // Every point collinear, or the loop encloses no area: a fan is as + // good as anything and cannot make it worse. + for (size_t i = 1; i + 1 < n; ++i) emit(loop[0], loop[i], loop[i + 1]); + return; + } + for (double& c : nrm) c /= len; + + // Any two axes spanning the plane will do; take the world axis least + // aligned with the normal so the projection never collapses. + const size_t drop = (std::abs(nrm[0]) > std::abs(nrm[1])) + ? ((std::abs(nrm[0]) > std::abs(nrm[2])) ? 0 : 2) + : ((std::abs(nrm[1]) > std::abs(nrm[2])) ? 1 : 2); + std::array axis = {0.0, 0.0, 0.0}; + axis[(drop + 1) % 3] = 1.0; + std::array u = {axis[1] * nrm[2] - axis[2] * nrm[1], axis[2] * nrm[0] - axis[0] * nrm[2], + axis[0] * nrm[1] - axis[1] * nrm[0]}; + const double ulen = std::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); + if (!(ulen > 1e-12)) { + for (size_t i = 1; i + 1 < n; ++i) emit(loop[0], loop[i], loop[i + 1]); + return; + } + for (double& c : u) c /= ulen; + const std::array w = {nrm[1] * u[2] - nrm[2] * u[1], nrm[2] * u[0] - nrm[0] * u[2], + nrm[0] * u[1] - nrm[1] * u[0]}; + + std::vector> flat(n); + for (size_t i = 0; i < n; ++i) { + const std::array& p = verts[loop[i]]; + flat[i] = {p[0] * u[0] + p[1] * u[1] + p[2] * u[2], p[0] * w[0] + p[1] * w[1] + p[2] * w[2]}; + } + + auto cross2 = [](const std::array& a, const std::array& b, + const std::array& c) { + return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]); + }; + double twiceArea = 0.0; + for (size_t i = 0; i < n; ++i) { + const std::array& a = flat[i]; + const std::array& b = flat[(i + 1) % n]; + twiceArea += a[0] * b[1] - b[0] * a[1]; + } + + std::vector idx(n); + for (size_t i = 0; i < n; ++i) idx[i] = i; + if (twiceArea < 0.0) std::reverse(idx.begin(), idx.end()); // work counter-clockwise + + auto inside = [&](const std::array& a, const std::array& b, + const std::array& c, const std::array& p) { + // Strictly inside, so a vertex sitting exactly on an edge does not + // veto an otherwise good ear. + const double d1 = cross2(a, b, p), d2 = cross2(b, c, p), d3 = cross2(c, a, p); + return d1 > 1e-12 && d2 > 1e-12 && d3 > 1e-12; + }; + + // Take the best-shaped ear available rather than the first one found. + // Any valid ear gives a correct triangulation, but on a face that is + // not flat the choice decides how the surface folds: first-found + // clipping strung long thin triangles across the curved end caps and + // came out 29% larger in area than the reference's tessellation of the + // same polygon. Preferring fat ears tracks the surface instead. + auto squareness = [&](size_t a, size_t b, size_t c) { + // Twice the area over the sum of the squared sides -- highest for + // an equilateral triangle, near zero for a sliver. + const std::array& p = verts[loop[a]]; + const std::array& q = verts[loop[b]]; + const std::array& r = verts[loop[c]]; + const std::array e1 = {q[0] - p[0], q[1] - p[1], q[2] - p[2]}; + const std::array e2 = {r[0] - p[0], r[1] - p[1], r[2] - p[2]}; + const std::array e3 = {r[0] - q[0], r[1] - q[1], r[2] - q[2]}; + const std::array x = {e1[1] * e2[2] - e1[2] * e2[1], e1[2] * e2[0] - e1[0] * e2[2], + e1[0] * e2[1] - e1[1] * e2[0]}; + const double area = std::sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); + const double sides = e1[0] * e1[0] + e1[1] * e1[1] + e1[2] * e1[2] + + e2[0] * e2[0] + e2[1] * e2[1] + e2[2] * e2[2] + + e3[0] * e3[0] + e3[1] * e3[1] + e3[2] * e3[2]; + return sides > 1e-18 ? area / sides : 0.0; + }; + + size_t guard = 0; + while (idx.size() > 3 && guard++ < n * n) { + size_t bestAt = idx.size(); + double bestScore = -1.0; + for (size_t i = 0; i < idx.size(); ++i) { + const size_t pi = idx[(i + idx.size() - 1) % idx.size()]; + const size_t ci = idx[i]; + const size_t ni = idx[(i + 1) % idx.size()]; + if (cross2(flat[pi], flat[ci], flat[ni]) <= 1e-12) continue; // reflex, not an ear + bool empty = true; + for (size_t other : idx) { + if (other == pi || other == ci || other == ni) continue; + if (inside(flat[pi], flat[ci], flat[ni], flat[other])) { empty = false; break; } + } + if (!empty) continue; + const double score = squareness(pi, ci, ni); + if (score > bestScore) { bestScore = score; bestAt = i; } + } + if (bestAt == idx.size()) break; // self-intersecting or otherwise unclippable + const size_t pi = idx[(bestAt + idx.size() - 1) % idx.size()]; + const size_t ci = idx[bestAt]; + const size_t ni = idx[(bestAt + 1) % idx.size()]; + emit(loop[pi], loop[ci], loop[ni]); + idx.erase(idx.begin() + static_cast(bestAt)); + } + // Whatever is left: three points, or a remainder no ear could be found + // in. A fan over the remainder is the best available answer. + for (size_t i = 1; i + 1 < idx.size(); ++i) emit(loop[idx[0]], loop[idx[i]], loop[idx[i + 1]]); +} + +} // namespace + CSGParams resolvePolyhedron(Evaluator& ev, const oscad::ModularCall& node, EvalContext& ctx) { auto [args, effCtx] = resolveCallArgs(ev, node.arguments, ctx); Value pointsArg = getArg(args, 0, "points", Value{}); @@ -292,14 +450,7 @@ CSGParams resolvePolyhedron(Evaluator& ev, const oscad::ModularCall& node, EvalC const size_t idx = static_cast(toDoubleLenient(idxVal)); remapped.push_back(idx < remap.size() ? remap[idx] : 0); } - for (size_t i = 1; i + 1 < remapped.size(); ++i) { - const size_t a = remapped[0], b = remapped[i + 1], c = remapped[i]; - if (a != b && b != c && a != c) { - tris.push_back(static_cast(a)); - tris.push_back(static_cast(b)); - tris.push_back(static_cast(c)); - } - } + triangulateFace(uniqueVerts, remapped, tris); } std::vector vertsValues; diff --git a/tests/test_booleans.cpp b/tests/test_booleans.cpp index 11d21ca..02b0353 100644 --- a/tests/test_booleans.cpp +++ b/tests/test_booleans.cpp @@ -227,3 +227,93 @@ TEST(Modifiers, NoShowOnlyMeansAllRolesPassThrough) { Evaluated e = evaluateSrc("%cube(1); cube(2);"); EXPECT_EQ(e.bodies.size(), 2u); } + +// -- polyhedron face triangulation --------------------------------------- + +// An L-shaped prism: its top and bottom are 6-gons with one reflex corner. +// Fan-triangulating a concave face lays triangles outside the polygon and +// winds some of them inside out, so the solid comes out the wrong size -- +// which is how BOSL2's nurbs_sheet() end caps ended up covering 281% of +// their own area with 15 of 32 triangles reversed. +// A U: vertex 0 cannot see the far arm, so a fan anchored there lays +// triangles across the notch, outside the solid. (An L is not enough on +// its own -- it happens to be star-shaped from its first vertex, so even +// a fan triangulates it correctly, which is why the volume check below +// needs this shape too.) +static const char* kUPrism = R"( +pts = [[0,0],[3,0],[3,3],[2,3],[2,1],[1,1],[1,3],[0,3]]; +polyhedron( + points = concat([for (p=pts) [p[0],p[1],0]], [for (p=pts) [p[0],p[1],1]]), + faces = concat( + [[for (i=[0:len(pts)-1]) i]], + [[for (i=[len(pts)-1:-1:0]) i+len(pts)]], + [for (i=[0:len(pts)-1]) [i, i+len(pts), (i+1)%len(pts)+len(pts), (i+1)%len(pts)]] + ) +); +)"; + +static const char* kLPrism = R"( +polyhedron( + points = [ + [0,0,0], [2,0,0], [2,1,0], [1,1,0], [1,2,0], [0,2,0], + [0,0,1], [2,0,1], [2,1,1], [1,1,1], [1,2,1], [0,2,1] + ], + faces = [ + [0,1,2,3,4,5], + [11,10,9,8,7,6], + [0,6,7,1], [1,7,8,2], [2,8,9,3], + [3,9,10,4], [4,10,11,5], [5,11,6,0] + ] +); +)"; + +TEST(PolyhedronFaces, AConcaveFaceGivesTheRightVolume) { + Evaluated e = evaluateSrc(kLPrism); + ASSERT_EQ(e.bodies.size(), 1u); + ASSERT_TRUE(e.bodies[0].body.has_value()); + // The L is 3 square units, extruded 1 high. + EXPECT_NEAR(e.bodies[0].body->Volume(), 3.0, 1e-9); +} + +TEST(PolyhedronFaces, AConcaveFaceGivesTheRightSurfaceArea) { + Evaluated e = evaluateSrc(kLPrism); + ASSERT_EQ(e.bodies.size(), 1u); + // 2 x 3 for the ends, plus a perimeter of 8 x height 1. + EXPECT_NEAR(e.bodies[0].body->SurfaceArea(), 14.0, 1e-9); +} + +TEST(PolyhedronFaces, EachFaceBecomesExactlyNMinusTwoTriangles) { + Evaluated e = evaluateSrc(kLPrism); + ASSERT_EQ(e.bodies.size(), 1u); + // Two 6-gons (4 each) and six quads (2 each). + const manifold::MeshGL mesh = e.bodies[0].body->GetMeshGL(); + EXPECT_EQ(mesh.triVerts.size() / 3, 4u + 4u + 6u * 2u); +} + +// A convex face has to keep working, and a triangle must not be disturbed. +TEST(PolyhedronFaces, AFaceWithAnUnseeableCornerGivesTheRightVolume) { + Evaluated e = evaluateSrc(kUPrism); + ASSERT_EQ(e.bodies.size(), 1u); + ASSERT_TRUE(e.bodies[0].body.has_value()); + // 3x3 less the 1x2 notch, one unit high. + EXPECT_NEAR(e.bodies[0].body->Volume(), 7.0, 1e-9); + // Two ends of 7, plus a 16-long perimeter one unit high. + EXPECT_NEAR(e.bodies[0].body->SurfaceArea(), 30.0, 1e-9); +} + +TEST(PolyhedronFaces, AConvexFaceIsUnchanged) { + Evaluated e = evaluateSrc("cube(2);"); + ASSERT_EQ(e.bodies.size(), 1u); + EXPECT_NEAR(e.bodies[0].body->Volume(), 8.0, 1e-9); + EXPECT_NEAR(e.bodies[0].body->SurfaceArea(), 24.0, 1e-9); +} + +TEST(PolyhedronFaces, ADegenerateFaceDoesNotThrowOrHang) { + // All points collinear: no plane to project onto. It must fall through + // rather than spin in the ear-clipping loop. + Evaluated e = evaluateSrc(R"( + polyhedron(points=[[0,0,0],[1,0,0],[2,0,0],[0,0,1]], + faces=[[0,1,2],[0,2,3],[0,3,1],[1,3,2]]); + )"); + EXPECT_EQ(e.bodies.size(), 1u); +}