Maintain manifoldness when converting between PolySet and ManifoldGeo…

…metry (#5020)

* Add constResult()/mutableResult() for explicit result creation
* Test case had polygon winding wrong
* Remove obsolete trustedPolySetToManifold()

src/geometry/GeometryEvaluator.cc

@@ -108,7 +108,7 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren(const Abstrac
   for (const auto& item : this->visitedchildren[node.index()]) {
     if (!isValidDim(item, dim)) break;
   }
-  if (dim == 2) return {std::shared_ptr<Geometry>(applyToChildren2D(node, op))};
+  if (dim == 2) return ResultObject::mutableResult(std::shared_ptr<Geometry>(applyToChildren2D(node, op)));
   else if (dim == 3) return applyToChildren3D(node, op);
   return {};
 }
@@ -124,15 +124,15 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren3D(const Abstr
   if (children.empty()) return {};
 
   if (op == OpenSCADOperator::HULL) {
-    return {std::shared_ptr<Geometry>(applyHull(children))};
+    return ResultObject::mutableResult(std::shared_ptr<Geometry>(applyHull(children)));
   } else if (op == OpenSCADOperator::FILL) {
     for (const auto& item : children) {
       LOG(message_group::Warning, item.first->modinst->location(), this->tree.getDocumentPath(), "fill() not yet implemented for 3D");
     }
   }
 
   // Only one child -> this is a noop
-  if (children.size() == 1) return {children.front().second};
+  if (children.size() == 1) return ResultObject::constResult(children.front().second);
 
   switch (op) {
   case OpenSCADOperator::MINKOWSKI:
@@ -142,8 +142,8 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren3D(const Abstr
       if (item.second && !item.second->isEmpty()) actualchildren.push_back(item);
     }
     if (actualchildren.empty()) return {};
-    if (actualchildren.size() == 1) return {actualchildren.front().second};
-    return {applyMinkowski(actualchildren)};
+    if (actualchildren.size() == 1) return ResultObject::constResult(actualchildren.front().second);
+    return ResultObject::constResult(applyMinkowski(actualchildren));
     break;
   }
   case OpenSCADOperator::UNION:
@@ -153,17 +153,17 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren3D(const Abstr
       if (item.second && !item.second->isEmpty()) actualchildren.push_back(item);
     }
     if (actualchildren.empty()) return {};
-    if (actualchildren.size() == 1) return {actualchildren.front().second};
+    if (actualchildren.size() == 1) return ResultObject::constResult(actualchildren.front().second);
 #ifdef ENABLE_MANIFOLD
     if (Feature::ExperimentalManifold.is_enabled()) {
-      return {ManifoldUtils::applyOperator3DManifold(actualchildren, op)};
+      return ResultObject::mutableResult(ManifoldUtils::applyOperator3DManifold(actualchildren, op));
     }
 #endif
 #ifdef ENABLE_CGAL
     else if (Feature::ExperimentalFastCsg.is_enabled()) {
-      return {std::shared_ptr<Geometry>(CGALUtils::applyUnion3DHybrid(actualchildren.begin(), actualchildren.end()))};
+      return ResultObject::mutableResult(std::shared_ptr<Geometry>(CGALUtils::applyUnion3DHybrid(actualchildren.begin(), actualchildren.end())));
     }
-    return {CGALUtils::applyUnion3D(actualchildren.begin(), actualchildren.end())};
+    return ResultObject::constResult(std::shared_ptr<const Geometry>(CGALUtils::applyUnion3D(actualchildren.begin(), actualchildren.end())));
 #else
     assert(false && "No boolean backend available");
 #endif
@@ -173,15 +173,15 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren3D(const Abstr
   {
 #ifdef ENABLE_MANIFOLD
     if (Feature::ExperimentalManifold.is_enabled()) {
-      return {ManifoldUtils::applyOperator3DManifold(children, op)};
+      return ResultObject::mutableResult(ManifoldUtils::applyOperator3DManifold(children, op));
     }
 #endif
 #ifdef ENABLE_CGAL
     if (Feature::ExperimentalFastCsg.is_enabled()) {
       // FIXME: It's annoying to have to disambiguate here:
-      return {std::shared_ptr<Geometry>(CGALUtils::applyOperator3DHybrid(children, op))};
+      return ResultObject::mutableResult(std::shared_ptr<Geometry>(CGALUtils::applyOperator3DHybrid(children, op)));
     }
-    return {CGALUtils::applyOperator3D(children, op)};
+    return ResultObject::constResult(CGALUtils::applyOperator3D(children, op));
 #else
     assert(false && "No boolean backend available");
     #endif
@@ -697,30 +697,21 @@ Response GeometryEvaluator::visit(State& state, const TransformNode& node)
         ResultObject res = applyToChildren(node, OpenSCADOperator::UNION);
         if ((geom = res.constptr())) {
           if (geom->getDimension() == 2) {
-            std::shared_ptr<const Polygon2d> polygons = std::dynamic_pointer_cast<const Polygon2d>(geom);
+            auto polygons =  std::dynamic_pointer_cast<Polygon2d>(res.asMutableGeometry());
             assert(polygons);
 
-            // If we got a const object, make a copy
-            std::shared_ptr<Polygon2d> newpoly;
-            if (res.isConst()) {
-              newpoly = std::make_shared<Polygon2d>(*polygons);
-            }
-            else {
-              newpoly = std::dynamic_pointer_cast<Polygon2d>(res.ptr());
-            }
-
             Transform2d mat2;
             mat2.matrix() <<
               node.matrix(0, 0), node.matrix(0, 1), node.matrix(0, 3),
               node.matrix(1, 0), node.matrix(1, 1), node.matrix(1, 3),
               node.matrix(3, 0), node.matrix(3, 1), node.matrix(3, 3);
-            newpoly->transform(mat2);
+            polygons->transform(mat2);
             // FIXME: We lose the transform if we copied a const geometry above. Probably similar issue in multiple places
             // A 2D transformation may flip the winding order of a polygon.
             // If that happens with a sanitized polygon, we need to reverse
             // the winding order for it to be correct.
-            if (newpoly->isSanitized() && mat2.matrix().determinant() <= 0) {
-              geom = ClipperUtils::sanitize(*newpoly);
+            if (polygons->isSanitized() && mat2.matrix().determinant() <= 0) {
+              geom = ClipperUtils::sanitize(*polygons);
             }
           } else if (geom->getDimension() == 3) {
             auto mutableGeom = res.asMutableGeometry();
@@ -1033,6 +1024,7 @@ static Outline2d splitOutlineByFn(
    Input to extrude should be sanitized. This means non-intersecting, correct winding order
    etc., the input coming from a library like Clipper.
  */
+ // FIXME: What happens if the input Polygon isn't manifold, or has coincident vertices?
 static std::unique_ptr<Geometry> extrudePolygon(const LinearExtrudeNode& node, const Polygon2d& poly)
 {
   bool non_linear = node.twist != 0 || node.scale_x != node.scale_y;
@@ -1478,12 +1470,9 @@ Response GeometryEvaluator::visit(State& state, const CgalAdvNode& node)
         geom = res.constptr();
         // If we added convexity, we need to pass it on
         if (geom && geom->getConvexity() != node.convexity) {
-          std::shared_ptr<Geometry> editablegeom;
-          // If we got a const object, make a copy
-          if (res.isConst()) editablegeom = geom->copy();
-          else editablegeom = res.ptr();
-          geom = editablegeom;
+          auto editablegeom = res.asMutableGeometry();
           editablegeom->setConvexity(node.convexity);
+          geom = editablegeom;
         }
         break;
       }
@@ -1537,6 +1526,7 @@ Response GeometryEvaluator::visit(State& state, const AbstractIntersectionNode&
 }
 
 #if defined(ENABLE_EXPERIMENTAL) && defined(ENABLE_CGAL)
+// FIXME: What is the convex/manifold situation of the resulting PolySet?
 static std::unique_ptr<Geometry> roofOverPolygon(const RoofNode& node, const Polygon2d& poly)
 {
   std::unique_ptr<PolySet> roof;

src/geometry/GeometryEvaluator.h

@@ -50,21 +50,30 @@ class GeometryEvaluator : public NodeVisitor
   class ResultObject
   {
 public:
+    // This makes it explicit if we want a const vs. non-const result.
+    // This is important to avoid inadvertently tagging a geometry as const when
+    // the underlying geometry is actually mutable. 
+    // The template trick, combined with private constructors, makes it possible
+    // to create a ResultObject containing a const, _only_ from const objects
+    // (i.e. no implicit conversion from non-const to const).
+    template<class T> static ResultObject constResult(std::shared_ptr<const T> geom) {return {geom};}
+    template<class T> static ResultObject mutableResult(std::shared_ptr<T> geom) {return {geom};}
+
     // Default constructor with nullptr can be used to represent empty geometry,
     // for example union() with no children, etc.
     ResultObject() : is_const(true) {}
-    ResultObject(std::shared_ptr<const Geometry> g) : is_const(true), const_pointer(std::move(g)) {}
-    ResultObject(std::shared_ptr<Geometry> g) : is_const(false), pointer(std::move(g)) {}
-    [[nodiscard]] bool isConst() const { return is_const; }
     std::shared_ptr<Geometry> ptr() { assert(!is_const); return pointer; }
     [[nodiscard]] std::shared_ptr<const Geometry> constptr() const {
       return is_const ? const_pointer : std::static_pointer_cast<const Geometry>(pointer);
     }
     std::shared_ptr<Geometry> asMutableGeometry() {
-      if (isConst()) return std::shared_ptr<Geometry>(constptr() ? constptr()->copy() : nullptr);
+      if (is_const) return {constptr() ? constptr()->copy() : nullptr};
       else return ptr();
     }
 private:
+    template<class T> ResultObject(std::shared_ptr<const T> g) : is_const(true), const_pointer(std::move(g)) {}
+    template<class T> ResultObject(std::shared_ptr<T> g) : is_const(false), pointer(std::move(g)) {}
+
     bool is_const;
     std::shared_ptr<Geometry> pointer;
     std::shared_ptr<const Geometry> const_pointer;

src/geometry/boolean_utils.cc

@@ -11,11 +11,11 @@
 #include <CGAL/version.h>
 #include <CGAL/convex_hull_3.h>
 #include "cgalutils.h"
-#endif
+#endif  // ENABLE_CGAL
 #ifdef ENABLE_MANIFOLD
 #include "ManifoldGeometry.h"
 #include "manifoldutils.h"
-#endif
+#endif  // ENABLE_MANIFOLD
 
 #include "Feature.h"
 #include "PolySet.h"
@@ -65,7 +65,7 @@ std::unique_ptr<PolySet> applyHull(const Geometry::Geometries& children)
           addPoint(CGALUtils::vector_convert<K::Point_3>(p));
           return false;
         });
-#endif
+#endif  // ENABLE_MANIFOLD
     } else if (const auto *ps = dynamic_cast<const PolySet*>(chgeom.get())) {
       addCapacity(ps->indices.size() * 3);
       for (const auto& p : ps->indices) {
@@ -88,6 +88,8 @@ std::unique_ptr<PolySet> applyHull(const Geometry::Geometries& children)
       PRINTDB("After hull facets: %d", r.size_of_facets());
       PRINTDB("After hull closed: %d", r.is_closed());
       PRINTDB("After hull valid: %d", r.is_valid());
+      // FIXME: Make sure PolySet is set to convex.
+      // FIXME: Can we guarantee a manifold PolySet here?
       return CGALUtils::createPolySetFromPolyhedron(r);
     } catch (const CGAL::Failure_exception& e) {
       LOG(message_group::Error, "CGAL error in applyHull(): %1$s", e.what());
@@ -98,14 +100,16 @@ std::unique_ptr<PolySet> applyHull(const Geometry::Geometries& children)
 
 /*!
    children cannot contain nullptr objects
+
+  FIXME: This shouldn't return const, but it does due to internal implementation details
  */
 std::shared_ptr<const Geometry> applyMinkowski(const Geometry::Geometries& children)
 {
 #if ENABLE_MANIFOLD
   if (Feature::ExperimentalManifold.is_enabled()) {
     return ManifoldUtils::applyMinkowskiManifold(children);
   }
-#endif
+#endif  // ENABLE_MANIFOLD
   if (Feature::ExperimentalFastCsg.is_enabled()) {
     return CGALUtils::applyMinkowskiHybrid(children);
   }
@@ -309,7 +313,7 @@ std::shared_ptr<const Geometry> applyMinkowski(const Geometry::Geometries& child
     return N;
   }
 }
-#else
+#else  // ENABLE_CGAL
 bool applyHull(const Geometry::Geometries& children, PolySet& result)
 {
   return false;
@@ -319,4 +323,4 @@ std::shared_ptr<const Geometry> applyMinkowski(const Geometry::Geometries& child
 {
   return std::make_shared<PolySet>(3);
 }
-#endif
+#endif  // ENABLE_CGAL

src/geometry/cgal/cgalutils-mesh.cc

@@ -41,9 +41,11 @@ bool createMeshFromPolySet(const PolySet& ps, TriangleMesh& mesh)
 template bool createMeshFromPolySet(const PolySet& ps, CGAL_HybridMesh& mesh);
 template bool createMeshFromPolySet(const PolySet& ps, CGAL_DoubleMesh& mesh);
 
+
 template <class TriangleMesh>
 std::unique_ptr<PolySet> createPolySetFromMesh(const TriangleMesh& mesh)
 {
+  //  FIXME: We may want to convert directly, without PolySetBuilder here, to maintain manifoldness, if possible.
   PolySetBuilder builder(0, mesh.number_of_faces()+ mesh.number_of_faces());
   for (const auto& f : mesh.faces()) {
     builder.appendPoly(mesh.degree(f));

src/geometry/manifold/ManifoldGeometry.cc

@@ -165,31 +165,7 @@ std::shared_ptr<manifold::Manifold> binOp(const manifold::Manifold& lhs, const m
   return std::make_shared<manifold::Manifold>(lhs.Boolean(rhs, opType));
 }
 
-ManifoldGeometry ManifoldGeometry::operator+(const ManifoldGeometry& other) const {
-  return {binOp(*this->manifold_, *other.manifold_, manifold::OpType::Add)};
-}
-
-void ManifoldGeometry::operator+=(ManifoldGeometry& other) {
-  manifold_ = binOp(*this->manifold_, *other.manifold_, manifold::OpType::Add);
-}
-
-ManifoldGeometry ManifoldGeometry::operator*(const ManifoldGeometry& other) const {
-  return {binOp(*this->manifold_, *other.manifold_, manifold::OpType::Intersect)};
-}
-
-void ManifoldGeometry::operator*=(ManifoldGeometry& other) {
-  manifold_ = binOp(*this->manifold_, *other.manifold_, manifold::OpType::Intersect);
-}
-
-ManifoldGeometry ManifoldGeometry::operator-(const ManifoldGeometry& other) const {
-  return {binOp(*this->manifold_, *other.manifold_, manifold::OpType::Subtract)};
-}
-
-void ManifoldGeometry::operator-=(ManifoldGeometry& other) {
-  manifold_ = binOp(*this->manifold_, *other.manifold_, manifold::OpType::Subtract);
-}
-
-std::shared_ptr<manifold::Manifold> minkowskiOp(const ManifoldGeometry& lhs, const ManifoldGeometry& rhs) {
+std::shared_ptr<ManifoldGeometry> minkowskiOp(const ManifoldGeometry& lhs, const ManifoldGeometry& rhs) {
 // FIXME: How to deal with operation not supported?
 #ifdef ENABLE_CGAL
   auto lhs_nef = std::shared_ptr<CGAL_Nef_polyhedron>(CGALUtils::createNefPolyhedronFromPolySet(*lhs.toPolySet()));
@@ -202,17 +178,25 @@ std::shared_ptr<manifold::Manifold> minkowskiOp(const ManifoldGeometry& lhs, con
   auto ps = PolySetUtils::getGeometryAsPolySet(lhs_nef);
   if (!ps) return {};
   else {
-    return ManifoldUtils::trustedPolySetToManifold(*ps);
+    return ManifoldUtils::createManifoldFromPolySet(*ps);
   }
 #endif
 }
 
-void ManifoldGeometry::minkowski(ManifoldGeometry& other) {
-  manifold_ = minkowskiOp(*this, other);
+ManifoldGeometry ManifoldGeometry::operator+(const ManifoldGeometry& other) const {
+  return {binOp(*this->manifold_, *other.manifold_, manifold::OpType::Add)};
+}
+
+ManifoldGeometry ManifoldGeometry::operator*(const ManifoldGeometry& other) const {
+  return {binOp(*this->manifold_, *other.manifold_, manifold::OpType::Intersect)};
+}
+
+ManifoldGeometry ManifoldGeometry::operator-(const ManifoldGeometry& other) const {
+  return {binOp(*this->manifold_, *other.manifold_, manifold::OpType::Subtract)};
 }
 
 ManifoldGeometry ManifoldGeometry::minkowski(const ManifoldGeometry& other) const {
-  return {minkowskiOp(*this, other)};
+  return {*minkowskiOp(*this, other)};
 }
 
 void ManifoldGeometry::transform(const Transform3d& mat) {

src/geometry/manifold/ManifoldGeometry.h

@@ -39,15 +39,6 @@ class ManifoldGeometry : public Geometry
   template <class Polyhedron>
   [[nodiscard]] std::shared_ptr<Polyhedron> toPolyhedron() const;
 
-  /*! In-place union. */
-  void operator+=(ManifoldGeometry& other);
-  /*! In-place intersection. */
-  void operator*=(ManifoldGeometry& other);
-  /*! In-place difference. */
-  void operator-=(ManifoldGeometry& other);
-  /*! In-place minkowksi operation. */
-  void minkowski(ManifoldGeometry& other);
-
   /*! union. */
   ManifoldGeometry operator+(const ManifoldGeometry& other) const;
   /*! intersection. */

src/geometry/manifold/manifold-applyops-minkowski.cc

@@ -119,7 +119,7 @@ std::shared_ptr<const Geometry> applyMinkowskiManifold(const Geometry::Geometrie
 
         if (minkowski_points.size() <= 3) {
           t.stop();
-          return std::make_shared<const ManifoldGeometry>();
+          return std::make_shared<ManifoldGeometry>();
         }
 
         t.stop();

src/geometry/manifold/manifold-applyops.cc

@@ -20,7 +20,7 @@ Location getLocation(const std::shared_ptr<const AbstractNode>& node)
    Applies op to all children and returns the result.
    The child list should be guaranteed to contain non-NULL 3D or empty Geometry objects
  */
-std::shared_ptr<const ManifoldGeometry> applyOperator3DManifold(const Geometry::Geometries& children, OpenSCADOperator op)
+std::shared_ptr<ManifoldGeometry> applyOperator3DManifold(const Geometry::Geometries& children, OpenSCADOperator op)
 {
   std::shared_ptr<ManifoldGeometry> geom;