[WIP] Add tolerance value (#997)

* renaming precision to uncertainty, and hiding

* add tolerance

* rename uncertainty to epsilon

* use epsilon instead of tolerance

* guard simplify topology for SetEpsilon

* format

* remove conflict markers

* fix builds

* use epsilon for barycentric

* prefer using epsilon

* update python bindings

* address comments

* format

* remove SetTolerance

Tolerance value is set before the mesh is created, when SimplfiyTopology
cannot be called.

* read/write MeshGL tolerance field

bindings/c/manifoldc.cpp

@@ -525,7 +525,7 @@ ManifoldBox *manifold_bounding_box(void *mem, ManifoldManifold *m) {
 }
 
 double manifold_precision(ManifoldManifold *m) {
-  return from_c(m)->Precision();
+  return from_c(m)->GetEpsilon();
 }
 
 uint32_t manifold_reserve_ids(uint32_t n) { return Manifold::ReserveIDs(n); }

bindings/python/examples/all_apis.py

@@ -90,7 +90,8 @@ def all_manifold():
     n = m.num_tri()
     n = m.num_vert()
     i = m.original_id()
-    p = m.precision()
+    p = m.get_tolerance()
+    pp = m.set_tolerance(0.0001)
     c = m.project()
     m = m.refine(2)
     m = m.refine_to_length(0.1)

bindings/python/manifold3d.cpp

@@ -216,7 +216,7 @@ NB_MODULE(manifold3d, m) {
         nb::arg("radius"), get_circular_segments__radius);
 
   m.def("triangulate", &Triangulate, nb::arg("polygons"),
-        nb::arg("precision") = -1, triangulate__polygons__precision);
+        nb::arg("epsilon") = -1, triangulate__polygons__epsilon);
 
   nb::class_<Manifold>(m, "Manifold")
       .def(nb::init<>(), manifold__manifold)
@@ -338,19 +338,21 @@ NB_MODULE(manifold3d, m) {
       .def("num_tri", &Manifold::NumTri, manifold__num_tri)
       .def("num_prop", &Manifold::NumProp, manifold__num_prop)
       .def("num_prop_vert", &Manifold::NumPropVert, manifold__num_prop_vert)
-      .def("precision", &Manifold::Precision, manifold__precision)
       .def("genus", &Manifold::Genus, manifold__genus)
       .def(
           "volume",
           [](const Manifold &self) { return self.GetProperties().volume; },
           "Get the volume of the manifold\n This is clamped to zero for a "
-          "given face if they are within the Precision().")
+          "given face if they are within the Epsilon().")
       .def(
           "surface_area",
           [](const Manifold &self) { return self.GetProperties().surfaceArea; },
           "Get the surface area of the manifold\n This is clamped to zero for "
-          "a given face if they are within the Precision().")
+          "a given face if they are within the Epsilon().")
       .def("original_id", &Manifold::OriginalID, manifold__original_id)
+      .def("get_tolerance", &Manifold::GetTolerance, manifold__get_tolerance)
+      .def("set_tolerance", &Manifold::SetTolerance,
+           manifold__set_tolerance__tolerance)
       .def("as_original", &Manifold::AsOriginal, manifold__as_original)
       .def("is_empty", &Manifold::IsEmpty, manifold__is_empty)
       .def("decompose", &Manifold::Decompose, manifold__decompose)
@@ -371,7 +373,7 @@ NB_MODULE(manifold3d, m) {
       .def(
           "project",
           [](const Manifold &self) {
-            return CrossSection(self.Project()).Simplify(self.Precision());
+            return CrossSection(self.Project()).Simplify(self.GetEpsilon());
           },
           manifold__project)
       .def("status", &Manifold::Status, manifold__status)

bindings/wasm/bindings.cpp

@@ -161,7 +161,7 @@ EMSCRIPTEN_BINDINGS(whatever) {
       .function("numProp", &Manifold::NumProp)
       .function("numPropVert", &Manifold::NumPropVert)
       .function("_boundingBox", &Manifold::BoundingBox)
-      .function("precision", &Manifold::Precision)
+      .function("tolerance", &Manifold::GetTolerance)
       .function("genus", &Manifold::Genus)
       .function("getProperties", &Manifold::GetProperties)
       .function("minGap", &Manifold::MinGap)

bindings/wasm/bindings.js

@@ -294,7 +294,7 @@ Module.setup = function() {
     const polygonsVec = this._Project();
     const result = new CrossSectionCtor(polygonsVec, fillRuleToInt('Positive'));
     disposePolygons(polygonsVec);
-    return result.simplify(this.precision);
+    return result.simplify(this.tolerance);
   };
 
   Module.Manifold.prototype.split = function(manifold) {

include/manifold/manifold.h

@@ -88,11 +88,12 @@ struct MeshGLP {
   /// as 4 * (3 * tri + i) + j, i < 3, j < 4, representing the tangent value
   /// Mesh.triVerts[tri][i] along the CCW edge. If empty, mesh is faceted.
   std::vector<Precision> halfedgeTangent;
-  /// The absolute precision of the vertex positions, based on accrued rounding
-  /// errors. When creating a Manifold, the precision used will be the maximum
+  /// Tolerance for mesh simplification.
+  /// When creating a Manifold, the precision used will be the maximum
   /// of this and a baseline precision from the size of the bounding box. Any
-  /// edge shorter than precision may be collapsed.
-  Precision precision = 0;
+  /// edge shorter than tolerance may be collapsed.
+  /// Tolerance may be enlarged when floating point error accumulates.
+  Precision tolerance = 0;
 
   MeshGLP() = default;
 
@@ -230,10 +231,10 @@ class Manifold {
   size_t NumProp() const;
   size_t NumPropVert() const;
   Box BoundingBox() const;
-  double Precision() const;
   int Genus() const;
   Properties GetProperties() const;
   double MinGap(const Manifold& other, double searchLength) const;
+  double GetTolerance() const;
   ///@}
 
   /** @name Mesh ID
@@ -266,6 +267,7 @@ class Manifold {
   Manifold Refine(int) const;
   Manifold RefineToLength(double) const;
   Manifold RefineToPrecision(double) const;
+  Manifold SetTolerance(double) const;
   ///@}
 
   /** @name Boolean
@@ -310,6 +312,8 @@ class Manifold {
   bool MatchesTriNormals() const;
   size_t NumDegenerateTris() const;
   size_t NumOverlaps(const Manifold& second) const;
+  double GetEpsilon() const;
+  Manifold SetEpsilon(double epsilon) const;
   ///@}
 
   struct Impl;

src/boolean_result.cpp

@@ -526,7 +526,7 @@ struct Barycentric {
   VecView<const Halfedge> halfedgeP;
   VecView<const Halfedge> halfedgeQ;
   VecView<const Halfedge> halfedgeR;
-  const double precision;
+  const double epsilon;
 
   void operator()(const int tri) {
     const TriRef refPQ = ref[tri];
@@ -543,7 +543,7 @@ struct Barycentric {
 
     for (const int i : {0, 1, 2}) {
       const int vert = halfedgeR[3 * tri + i].startVert;
-      uvw[3 * tri + i] = GetBarycentric(vertPosR[vert], triPos, precision);
+      uvw[3 * tri + i] = GetBarycentric(vertPosR[vert], triPos, epsilon);
     }
   }
 };
@@ -564,7 +564,7 @@ void CreateProperties(Manifold::Impl &outR, const Manifold::Impl &inP,
   for_each_n(autoPolicy(numTri, 1e4), countAt(0), numTri,
              Barycentric({bary, outR.meshRelation_.triRef, inP.vertPos_,
                           inQ.vertPos_, outR.vertPos_, inP.halfedge_,
-                          inQ.halfedge_, outR.halfedge_, outR.precision_}));
+                          inQ.halfedge_, outR.halfedge_, outR.epsilon_}));
 
   using Entry = std::pair<ivec3, int>;
   int idMissProp = outR.NumVert();
@@ -738,7 +738,8 @@ Manifold::Impl Boolean3::Result(OpType op) const {
 
   if (numVertR == 0) return outR;
 
-  outR.precision_ = std::max(inP_.precision_, inQ_.precision_);
+  outR.epsilon_ = std::max(inP_.epsilon_, inQ_.epsilon_);
+  outR.tolerance_ = std::max(inP_.tolerance_, inQ_.tolerance_);
 
   outR.vertPos_.resize(numVertR);
   // Add vertices, duplicating for inclusion numbers not in [-1, 1].

src/constructors.cpp

@@ -407,8 +407,7 @@ Manifold Manifold::Revolve(const Polygons& crossSection, int circularSegments,
 
   // Add front and back triangles if not a full revolution.
   if (!isFullRevolution) {
-    std::vector<ivec3> frontTriangles =
-        Triangulate(polygons, pImpl_->precision_);
+    std::vector<ivec3> frontTriangles = Triangulate(polygons, pImpl_->epsilon_);
     for (auto& t : frontTriangles) {
       triVerts.push_back({startPoses[t.x], startPoses[t.y], startPoses[t.z]});
     }
@@ -468,7 +467,8 @@ std::vector<Manifold> Manifold::Decompose() const {
   for (int i = 0; i < numComponents; ++i) {
     auto impl = std::make_shared<Impl>();
     // inherit original object's precision
-    impl->precision_ = pImpl_->precision_;
+    impl->epsilon_ = pImpl_->epsilon_;
+    impl->tolerance_ = pImpl_->tolerance_;
 
     Vec<int> vertNew2Old(numVert);
     const int nVert =

src/csg_tree.cpp

@@ -172,7 +172,8 @@ CsgNodeType CsgLeafNode::GetNodeType() const { return CsgNodeType::Leaf; }
 Manifold::Impl CsgLeafNode::Compose(
     const std::vector<std::shared_ptr<CsgLeafNode>> &nodes) {
   ZoneScoped;
-  double precision = -1;
+  double epsilon = -1;
+  double tolerance = -1;
   int numVert = 0;
   int numEdge = 0;
   int numTri = 0;
@@ -191,11 +192,12 @@ Manifold::Impl CsgLeafNode::Compose(
     double nodeOldScale = node->pImpl_->bBox_.Scale();
     double nodeNewScale =
         node->pImpl_->bBox_.Transform(node->transform_).Scale();
-    double nodePrecision = node->pImpl_->precision_;
-    nodePrecision *= std::max(1.0, nodeNewScale / nodeOldScale);
-    nodePrecision = std::max(nodePrecision, kTolerance * nodeNewScale);
-    if (!std::isfinite(nodePrecision)) nodePrecision = -1;
-    precision = std::max(precision, nodePrecision);
+    double nodeEpsilon = node->pImpl_->epsilon_;
+    nodeEpsilon *= std::max(1.0, nodeNewScale / nodeOldScale);
+    nodeEpsilon = std::max(nodeEpsilon, kTolerance * nodeNewScale);
+    if (!std::isfinite(nodeEpsilon)) nodeEpsilon = -1;
+    epsilon = std::max(epsilon, nodeEpsilon);
+    tolerance = std::max(tolerance, node->pImpl_->tolerance_);
 
     vertIndices.push_back(numVert);
     edgeIndices.push_back(numEdge * 2);
@@ -212,7 +214,8 @@ Manifold::Impl CsgLeafNode::Compose(
   }
 
   Manifold::Impl combined;
-  combined.precision_ = precision;
+  combined.epsilon_ = epsilon;
+  combined.tolerance_ = tolerance;
   combined.vertPos_.resize(numVert);
   combined.halfedge_.resize(2 * numEdge);
   combined.faceNormal_.resize(numTri);

src/edge_op.cpp

@@ -206,7 +206,7 @@ void Manifold::Impl::SimplifyTopology() {
   {
     ZoneScopedN("CollapseShortEdge");
     numFlagged = 0;
-    ShortEdge se{halfedge_, vertPos_, precision_};
+    ShortEdge se{halfedge_, vertPos_, tolerance_};
     for_each_n(policy, countAt(0_uz), nbEdges,
                [&](size_t i) { bFlags[i] = se(i); });
     for (size_t i = 0; i < nbEdges; ++i) {
@@ -250,7 +250,7 @@ void Manifold::Impl::SimplifyTopology() {
   {
     ZoneScopedN("RecursiveEdgeSwap");
     numFlagged = 0;
-    SwappableEdge se{halfedge_, vertPos_, faceNormal_, precision_};
+    SwappableEdge se{halfedge_, vertPos_, faceNormal_, tolerance_};
     for_each_n(policy, countAt(0_uz), nbEdges,
                [&](size_t i) { bFlags[i] = se(i); });
     std::vector<int> edgeSwapStack;
@@ -470,7 +470,7 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
   const vec3 pNew = vertPos_[endVert];
   const vec3 pOld = vertPos_[toRemove.startVert];
   const vec3 delta = pNew - pOld;
-  const bool shortEdge = la::dot(delta, delta) < precision_ * precision_;
+  const bool shortEdge = la::dot(delta, delta) < tolerance_ * tolerance_;
 
   // Orbit endVert
   int current = halfedge_[tri0edge[1]].pairedHalfedge;
@@ -502,14 +502,14 @@ void Manifold::Impl::CollapseEdge(const int edge, std::vector<int>& edges) {
           // Don't collapse if the edges separating the faces are not colinear
           // (can happen when the two faces are coplanar).
           if (CCW(projection * pOld, projection * pLast, projection * pNew,
-                  precision_) != 0)
+                  epsilon_) != 0)
             return;
         }
       }
 
       // Don't collapse edge if it would cause a triangle to invert.
       if (CCW(projection * pNext, projection * pLast, projection * pNew,
-              precision_) < 0)
+              epsilon_) < 0)
         return;
 
       pLast = pNext;
@@ -582,7 +582,7 @@ void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag,
   for (int i : {0, 1, 2})
     v[i] = projection * vertPos_[halfedge_[tri0edge[i]].startVert];
   // Only operate on the long edge of a degenerate triangle.
-  if (CCW(v[0], v[1], v[2], precision_) > 0 || !Is01Longest(v[0], v[1], v[2]))
+  if (CCW(v[0], v[1], v[2], tolerance_) > 0 || !Is01Longest(v[0], v[1], v[2]))
     return;
 
   // Switch to neighbor's projection.
@@ -644,12 +644,12 @@ void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag,
   };
 
   // Only operate if the other triangles are not degenerate.
-  if (CCW(v[1], v[0], v[3], precision_) <= 0) {
+  if (CCW(v[1], v[0], v[3], tolerance_) <= 0) {
     if (!Is01Longest(v[1], v[0], v[3])) return;
     // Two facing, long-edge degenerates can swap.
     SwapEdge();
     const vec2 e23 = v[3] - v[2];
-    if (la::dot(e23, e23) < precision_ * precision_) {
+    if (la::dot(e23, e23) < tolerance_ * tolerance_) {
       tag++;
       CollapseEdge(tri0edge[2], edges);
       edges.resize(0);
@@ -660,8 +660,8 @@ void Manifold::Impl::RecursiveEdgeSwap(const int edge, int& tag,
                                                  tri0edge[1], tri0edge[0]});
     }
     return;
-  } else if (CCW(v[0], v[3], v[2], precision_) <= 0 ||
-             CCW(v[1], v[2], v[3], precision_) <= 0) {
+  } else if (CCW(v[0], v[3], v[2], tolerance_) <= 0 ||
+             CCW(v[1], v[2], v[3], tolerance_) <= 0) {
     return;
   }
   // Normal path

src/face_op.cpp

@@ -79,7 +79,7 @@ void Manifold::Impl::Face2Tri(const Vec<int>& faceEdge,
       auto triCCW = [&projection, this](const ivec3 tri) {
         return CCW(projection * this->vertPos_[tri[0]],
                    projection * this->vertPos_[tri[1]],
-                   projection * this->vertPos_[tri[2]], precision_) >= 0;
+                   projection * this->vertPos_[tri[2]], epsilon_) >= 0;
       };
 
       ivec3 tri0(halfedge_[firstEdge].startVert, halfedge_[firstEdge].endVert,
@@ -130,7 +130,7 @@ void Manifold::Impl::Face2Tri(const Vec<int>& faceEdge,
     const PolygonsIdx polys =
         Face2Polygons(halfedge_.cbegin() + faceEdge[face],
                       halfedge_.cbegin() + faceEdge[face + 1], projection);
-    return TriangulateIdx(polys, precision_);
+    return TriangulateIdx(polys, epsilon_);
   };
 #if (MANIFOLD_PAR == 1) && __has_include(<tbb/tbb.h>)
   tbb::task_group group;