diff --git a/bindings/c/manifoldc.cpp b/bindings/c/manifoldc.cpp
index 20930cc20..2887f5683 100644
--- a/bindings/c/manifoldc.cpp
+++ b/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); }
diff --git a/bindings/python/examples/all_apis.py b/bindings/python/examples/all_apis.py
index 027b4acf2..c9fcff831 100644
--- a/bindings/python/examples/all_apis.py
+++ b/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)
diff --git a/bindings/python/manifold3d.cpp b/bindings/python/manifold3d.cpp
index c7884efeb..cafa89a72 100644
--- a/bindings/python/manifold3d.cpp
+++ b/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)
diff --git a/bindings/wasm/bindings.cpp b/bindings/wasm/bindings.cpp
index c75edd20f..d8d4b1815 100644
--- a/bindings/wasm/bindings.cpp
+++ b/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)
diff --git a/bindings/wasm/bindings.js b/bindings/wasm/bindings.js
index 7d32deb4a..fd2b23df6 100644
--- a/bindings/wasm/bindings.js
+++ b/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) {
diff --git a/include/manifold/manifold.h b/include/manifold/manifold.h
index ffe8bef26..10597a618 100644
--- a/include/manifold/manifold.h
+++ b/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;
diff --git a/src/boolean_result.cpp b/src/boolean_result.cpp
index 4e68f1b24..63876b300 100644
--- a/src/boolean_result.cpp
+++ b/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].
diff --git a/src/constructors.cpp b/src/constructors.cpp
index 34b89eddd..f90ab8834 100644
--- a/src/constructors.cpp
+++ b/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 =
diff --git a/src/csg_tree.cpp b/src/csg_tree.cpp
index 3e6cbe777..a01804a59 100644
--- a/src/csg_tree.cpp
+++ b/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);
diff --git a/src/edge_op.cpp b/src/edge_op.cpp
index 9398deade..197585e58 100644
--- a/src/edge_op.cpp
+++ b/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
diff --git a/src/face_op.cpp b/src/face_op.cpp
index eb65d278e..15cab42a8 100644
--- a/src/face_op.cpp
+++ b/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;
diff --git a/src/impl.cpp b/src/impl.cpp
index 12eb7bdb9..69fee504d 100644
--- a/src/impl.cpp
+++ b/src/impl.cpp
@@ -51,7 +51,6 @@ struct AssignNormals {
   VecView<vec3> vertNormal;
   VecView<const vec3> vertPos;
   VecView<const Halfedge> halfedges;
-  const double precision;
 
   void operator()(const int face) {
     vec3& triNormal = faceNormal[face];
@@ -99,7 +98,8 @@ struct CoplanarEdge {
   VecView<const TriRef> triRef;
   VecView<const ivec3> triProp;
   const int numProp;
-  const double precision;
+  const double epsilon;
+  const double tolerance;
 
   // FIXME: race condition
   void operator()(const int edgeIdx) {
@@ -139,11 +139,11 @@ struct CoplanarEdge {
     triArea[edgeFace] = area;
     triArea[pairFace] = areaPair;
     // Don't link degenerate triangles
-    if (area < length * precision || areaPair < lengthPair * precision) return;
+    if (area < length * epsilon || areaPair < lengthPair * epsilon) return;
 
     const double volume = std::abs(la::dot(normal, pairVec));
     // Only operate on coplanar triangles
-    if (volume > std::max(area, areaPair) * precision) return;
+    if (volume > std::max(area, areaPair) * tolerance) return;
 
     face2face[edgeIdx] = std::make_pair(edgeFace, pairFace);
   }
@@ -154,7 +154,7 @@ struct CheckCoplanarity {
   VecView<const Halfedge> halfedge;
   VecView<const vec3> vertPos;
   std::vector<int>* components;
-  const double precision;
+  const double tolerance;
 
   void operator()(int tri) {
     const int component = (*components)[tri];
@@ -171,7 +171,7 @@ struct CheckCoplanarity {
       // If any component vertex is not coplanar with the component's reference
       // triangle, unmark the entire component so that none of its triangles are
       // marked coplanar.
-      if (std::abs(la::dot(normal, vert - origin)) > precision) {
+      if (std::abs(la::dot(normal, vert - origin)) > tolerance) {
         reinterpret_cast<std::atomic<int>*>(&comp2tri[component])
             ->store(-1, std::memory_order_relaxed);
         break;
@@ -359,7 +359,7 @@ void Manifold::Impl::CreateFaces() {
   for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), halfedge_.size(),
              CoplanarEdge({face2face, triArea, halfedge_, vertPos_,
                            meshRelation_.triRef, meshRelation_.triProperties,
-                           meshRelation_.numProp, precision_}));
+                           meshRelation_.numProp, epsilon_, tolerance_}));
 
   std::vector<int> components;
   const int numComponent = GetLabels(components, face2face, NumTri());
@@ -376,7 +376,7 @@ void Manifold::Impl::CreateFaces() {
 
   for_each_n(autoPolicy(halfedge_.size(), 1e4), countAt(0), NumTri(),
              CheckCoplanarity(
-                 {comp2tri, halfedge_, vertPos_, &components, precision_}));
+                 {comp2tri, halfedge_, vertPos_, &components, tolerance_}));
 
   Vec<TriRef>& triRef = meshRelation_.triRef;
   for (size_t tri = 0; tri < NumTri(); ++tri) {
@@ -498,7 +498,7 @@ void Manifold::Impl::WarpBatch(std::function<void(VecView<vec3>)> warpFunc) {
   Update();
   faceNormal_.resize(0);  // force recalculation of triNormal
   CalculateNormals();
-  SetPrecision();
+  SetEpsilon();
   Finish();
   CreateFaces();
   meshRelation_.originalID = -1;
@@ -515,7 +515,8 @@ Manifold::Impl Manifold::Impl::Transform(const mat3x4& transform_) const {
   }
   result.collider_ = collider_;
   result.meshRelation_ = meshRelation_;
-  result.precision_ = precision_;
+  result.epsilon_ = epsilon_;
+  result.tolerance_ = tolerance_;
   result.bBox_ = bBox_;
   result.halfedge_ = halfedge_;
   result.halfedgeTangent_.resize(halfedgeTangent_.size());
@@ -556,9 +557,9 @@ Manifold::Impl Manifold::Impl::Transform(const mat3x4& transform_) const {
 
   result.CalculateBBox();
   // Scale the precision by the norm of the 3x3 portion of the transform.
-  result.precision_ *= SpectralNorm(mat3(transform_));
+  result.epsilon_ *= SpectralNorm(mat3(transform_));
   // Maximum of inherited precision loss and translational precision loss.
-  result.SetPrecision(result.precision_);
+  result.SetEpsilon(result.epsilon_);
   return result;
 }
 
@@ -566,8 +567,9 @@ Manifold::Impl Manifold::Impl::Transform(const mat3x4& transform_) const {
  * Sets the precision based on the bounding box, and limits its minimum value
  * by the optional input.
  */
-void Manifold::Impl::SetPrecision(double minPrecision) {
-  precision_ = MaxPrecision(minPrecision, bBox_);
+void Manifold::Impl::SetEpsilon(double minEpsilon) {
+  epsilon_ = MaxEpsilon(minEpsilon, bBox_);
+  tolerance_ = std::max(tolerance_, epsilon_);
 }
 
 /**
@@ -593,13 +595,13 @@ void Manifold::Impl::CalculateNormals() {
     calculateTriNormal = true;
   }
   if (calculateTriNormal)
-    for_each_n(policy, countAt(0), NumTri(),
-               AssignNormals<true>({faceNormal_, vertNormal_, vertPos_,
-                                    halfedge_, precision_}));
+    for_each_n(
+        policy, countAt(0), NumTri(),
+        AssignNormals<true>({faceNormal_, vertNormal_, vertPos_, halfedge_}));
   else
-    for_each_n(policy, countAt(0), NumTri(),
-               AssignNormals<false>({faceNormal_, vertNormal_, vertPos_,
-                                     halfedge_, precision_}));
+    for_each_n(
+        policy, countAt(0), NumTri(),
+        AssignNormals<false>({faceNormal_, vertNormal_, vertPos_, halfedge_}));
   for_each(policy, vertNormal_.begin(), vertNormal_.end(),
            [](vec3& v) { v = SafeNormalize(v); });
 }
diff --git a/src/impl.h b/src/impl.h
index a27d8f90f..f8671cef4 100644
--- a/src/impl.h
+++ b/src/impl.h
@@ -45,7 +45,8 @@ struct Manifold::Impl {
   };
 
   Box bBox_;
-  double precision_ = -1;
+  double epsilon_ = -1;
+  double tolerance_ = -1;
   Error status_ = Error::NoError;
   Vec<vec3> vertPos_;
   Vec<Halfedge> halfedge_;
@@ -110,6 +111,7 @@ struct Manifold::Impl {
     const auto numProp = meshGL.numProp - 3;
     meshRelation_.numProp = numProp;
     meshRelation_.properties.resize(meshGL.NumVert() * numProp);
+    tolerance_ = meshGL.tolerance;
     // This will have unreferenced duplicate positions that will be removed by
     // Impl::RemoveUnreferencedVerts().
     vertPos_.resize(meshGL.NumVert());
@@ -200,7 +202,7 @@ struct Manifold::Impl {
       MarkFailure(Error::NonFiniteVertex);
       return;
     }
-    SetPrecision(meshGL.precision);
+    SetEpsilon();
 
     SplitPinchedVerts();
 
@@ -275,7 +277,7 @@ struct Manifold::Impl {
   void CalculateBBox();
   bool IsFinite() const;
   bool IsIndexInBounds(VecView<const ivec3> triVerts) const;
-  void SetPrecision(double minPrecision = -1);
+  void SetEpsilon(double minEpsilon = -1);
   bool IsManifold() const;
   bool Is2Manifold() const;
   bool MatchesTriNormals() const;
diff --git a/src/manifold.cpp b/src/manifold.cpp
index 02ddf9188..565437436 100644
--- a/src/manifold.cpp
+++ b/src/manifold.cpp
@@ -70,10 +70,9 @@ MeshGLP<Precision, I> GetMeshGLImpl(const manifold::Manifold::Impl& impl,
       !isOriginal && la::all(la::greater(normalIdx, ivec3(2)));
 
   MeshGLP<Precision, I> out;
-  out.precision =
-      std::max(impl.precision_,
-               std::numeric_limits<Precision>::epsilon() * impl.bBox_.Scale());
   out.numProp = 3 + numProp;
+  out.tolerance = std::max(std::numeric_limits<Precision>::epsilon(),
+                           static_cast<Precision>(impl.tolerance_));
   out.triVerts.resize(3 * numTri);
 
   const int numHalfedge = impl.halfedgeTangent_.size();
@@ -374,14 +373,46 @@ size_t Manifold::NumPropVert() const {
 Box Manifold::BoundingBox() const { return GetCsgLeafNode().GetImpl()->bBox_; }
 
 /**
- * Returns the precision of this Manifold's vertices, which tracks the
+ * Returns the epsilon value of this Manifold's vertices, which tracks the
  * approximate rounding error over all the transforms and operations that have
- * led to this state. Any triangles that are colinear within this precision are
- * considered degenerate and removed. This is the value of &epsilon; defining
+ * led to this state. This is the value of &epsilon; defining
  * [&epsilon;-valid](https://github.com/elalish/manifold/wiki/Manifold-Library#definition-of-%CE%B5-valid).
  */
-double Manifold::Precision() const {
-  return GetCsgLeafNode().GetImpl()->precision_;
+double Manifold::GetEpsilon() const {
+  return GetCsgLeafNode().GetImpl()->epsilon_;
+}
+
+Manifold Manifold::SetEpsilon(double epsilon) const {
+  auto impl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  auto oldTolerance = impl->tolerance_;
+  impl->SetEpsilon(epsilon);
+  if (impl->tolerance_ > oldTolerance) impl->SimplifyTopology();
+  return Manifold(impl);
+}
+
+/**
+ * Returns the tolerance of this Manifold's vertices.
+ * Edges shorter than this tolerance value will be collapsed.
+ */
+double Manifold::GetTolerance() const {
+  return GetCsgLeafNode().GetImpl()->tolerance_;
+}
+
+/**
+ * Return a copy of the manifold with the set tolerance value.
+ * This performs mesh simplification when the tolerance value is increased.
+ */
+Manifold Manifold::SetTolerance(double tolerance) const {
+  auto impl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  if (tolerance > impl->tolerance_) {
+    impl->tolerance_ = tolerance;
+    impl->SimplifyTopology();
+  } else {
+    // for reducing tolerance, we need to make sure it is still at least
+    // equal to epsilon.
+    impl->tolerance_ = std::max(impl->epsilon_, tolerance);
+  }
+  return Manifold(impl);
 }
 
 /**
diff --git a/src/polygon.cpp b/src/polygon.cpp
index 8f8fcdaa9..ae0ef5b3a 100644
--- a/src/polygon.cpp
+++ b/src/polygon.cpp
@@ -945,39 +945,39 @@ namespace manifold {
  * @param polys The set of polygons, wound CCW and representing multiple
  * polygons and/or holes. These have 2D-projected positions as well as
  * references back to the original vertices.
- * @param precision The value of &epsilon;, bounding the uncertainty of the
+ * @param epsilon The value of &epsilon;, bounding the uncertainty of the
  * input.
  * @return std::vector<ivec3> The triangles, referencing the original
  * vertex indicies.
  */
-std::vector<ivec3> TriangulateIdx(const PolygonsIdx &polys, double precision) {
+std::vector<ivec3> TriangulateIdx(const PolygonsIdx &polys, double epsilon) {
   std::vector<ivec3> triangles;
-  double updatedPrecision = precision;
+  double updatedEpsilon = epsilon;
 #ifdef MANIFOLD_EXCEPTIONS
   try {
 #endif
-    if (IsConvex(polys, precision)) {  // fast path
+    if (IsConvex(polys, epsilon)) {  // fast path
       triangles = TriangulateConvex(polys);
     } else {
-      EarClip triangulator(polys, precision);
+      EarClip triangulator(polys, epsilon);
       triangles = triangulator.Triangulate();
-      updatedPrecision = triangulator.GetPrecision();
+      updatedEpsilon = triangulator.GetPrecision();
     }
 #ifdef MANIFOLD_EXCEPTIONS
 #ifdef MANIFOLD_DEBUG
     if (params.intermediateChecks) {
       CheckTopology(triangles, polys);
       if (!params.processOverlaps) {
-        CheckGeometry(triangles, polys, 2 * updatedPrecision);
+        CheckGeometry(triangles, polys, 2 * updatedEpsilon);
       }
     }
   } catch (const geometryErr &e) {
     if (!params.suppressErrors) {
-      PrintFailure(e, polys, triangles, updatedPrecision);
+      PrintFailure(e, polys, triangles, updatedEpsilon);
     }
     throw;
   } catch (const std::exception &e) {
-    PrintFailure(e, polys, triangles, updatedPrecision);
+    PrintFailure(e, polys, triangles, updatedEpsilon);
     throw;
 #else
   } catch (const std::exception &e) {
@@ -994,12 +994,12 @@ std::vector<ivec3> TriangulateIdx(const PolygonsIdx &polys, double precision) {
  *
  * @param polygons The set of polygons, wound CCW and representing multiple
  * polygons and/or holes.
- * @param precision The value of &epsilon;, bounding the uncertainty of the
+ * @param epsilon The value of &epsilon;, bounding the uncertainty of the
  * input.
  * @return std::vector<ivec3> The triangles, referencing the original
  * polygon points in order.
  */
-std::vector<ivec3> Triangulate(const Polygons &polygons, double precision) {
+std::vector<ivec3> Triangulate(const Polygons &polygons, double epsilon) {
   int idx = 0;
   PolygonsIdx polygonsIndexed;
   for (const auto &poly : polygons) {
@@ -1009,7 +1009,7 @@ std::vector<ivec3> Triangulate(const Polygons &polygons, double precision) {
     }
     polygonsIndexed.push_back(simpleIndexed);
   }
-  return TriangulateIdx(polygonsIndexed, precision);
+  return TriangulateIdx(polygonsIndexed, epsilon);
 }
 
 ExecutionParams &PolygonParams() { return params; }
diff --git a/src/properties.cpp b/src/properties.cpp
index ab42c597b..17f3658d1 100644
--- a/src/properties.cpp
+++ b/src/properties.cpp
@@ -229,7 +229,7 @@ bool Manifold::Impl::Is2Manifold() const {
 bool Manifold::Impl::MatchesTriNormals() const {
   if (halfedge_.size() == 0 || faceNormal_.size() != NumTri()) return true;
   return all_of(countAt(0_uz), countAt(NumTri()),
-                CheckCCW({halfedge_, vertPos_, faceNormal_, 2 * precision_}));
+                CheckCCW({halfedge_, vertPos_, faceNormal_, 2 * epsilon_}));
 }
 
 /**
@@ -239,7 +239,7 @@ int Manifold::Impl::NumDegenerateTris() const {
   if (halfedge_.size() == 0 || faceNormal_.size() != NumTri()) return true;
   return count_if(
       countAt(0_uz), countAt(NumTri()),
-      CheckCCW({halfedge_, vertPos_, faceNormal_, -1 * precision_ / 2}));
+      CheckCCW({halfedge_, vertPos_, faceNormal_, -1 * epsilon_ / 2}));
 }
 
 Properties Manifold::Impl::GetProperties() const {
@@ -251,8 +251,7 @@ Properties Manifold::Impl::GetProperties() const {
   double areaCompensation = 0;
   double volumeCompensation = 0;
   for (size_t i = 0; i < NumTri(); ++i) {
-    auto [area1, volume1] =
-        FaceAreaVolume({halfedge_, vertPos_, precision_})(i);
+    auto [area1, volume1] = FaceAreaVolume({halfedge_, vertPos_, epsilon_})(i);
     const double t1 = area + area1;
     const double t2 = volume + volume1;
     areaCompensation += (area - t1) + area1;
diff --git a/src/quickhull.cpp b/src/quickhull.cpp
index 2a5f27e95..8c4dc6a0d 100644
--- a/src/quickhull.cpp
+++ b/src/quickhull.cpp
@@ -847,7 +847,7 @@ void Manifold::Impl::Hull(VecView<vec3> vertPos) {
   QuickHull qh(vertPos);
   std::tie(halfedge_, vertPos_) = qh.buildMesh();
   CalculateBBox();
-  SetPrecision(bBox_.Scale() * kTolerance);
+  SetEpsilon(bBox_.Scale() * kTolerance);
   CalculateNormals();
   InitializeOriginal();
   Finish();
diff --git a/src/shared.h b/src/shared.h
index ebf42348c..4bd72e40c 100644
--- a/src/shared.h
+++ b/src/shared.h
@@ -29,9 +29,9 @@ inline vec3 SafeNormalize(vec3 v) {
   return std::isfinite(v.x) ? v : vec3(0.0);
 }
 
-inline double MaxPrecision(double minPrecision, const Box& bBox) {
-  double precision = std::max(minPrecision, kTolerance * bBox.Scale());
-  return std::isfinite(precision) ? precision : -1;
+inline double MaxEpsilon(double minEpsilon, const Box& bBox) {
+  double epsilon = std::max(minEpsilon, kTolerance * bBox.Scale());
+  return std::isfinite(epsilon) ? epsilon : -1;
 }
 
 inline int NextHalfedge(int current) {
@@ -70,7 +70,7 @@ inline mat2x3 GetAxisAlignedProjection(vec3 normal) {
 }
 
 inline vec3 GetBarycentric(const vec3& v, const mat3& triPos,
-                           double precision) {
+                           double tolerance) {
   const mat3 edges(triPos[2] - triPos[1], triPos[0] - triPos[2],
                    triPos[1] - triPos[0]);
   const vec3 d2(la::dot(edges[0], edges[0]), la::dot(edges[1], edges[1]),
@@ -80,7 +80,7 @@ inline vec3 GetBarycentric(const vec3& v, const mat3& triPos,
                                                       : 2;
   const vec3 crossP = la::cross(edges[0], edges[1]);
   const double area2 = la::dot(crossP, crossP);
-  const double tol2 = precision * precision;
+  const double tol2 = tolerance * tolerance;
 
   vec3 uvw(0.0);
   for (const int i : {0, 1, 2}) {
diff --git a/src/sort.cpp b/src/sort.cpp
index ae32c147b..1ef2f6c25 100644
--- a/src/sort.cpp
+++ b/src/sort.cpp
@@ -143,25 +143,25 @@ bool MergeMeshGLP(MeshGLP<Precision, I>& mesh) {
     bBox.min[i] = minMax.first;
     bBox.max[i] = minMax.second;
   }
-  mesh.precision = MaxPrecision(mesh.precision, bBox);
-  if (mesh.precision < 0) return false;
+  auto epsilon = MaxEpsilon(0, bBox);
   auto policy = autoPolicy(numOpenVert, 1e5);
   Vec<Box> vertBox(numOpenVert);
   Vec<uint32_t> vertMorton(numOpenVert);
 
-  for_each_n(policy, countAt(0), numOpenVert,
-             [&vertMorton, &vertBox, &openVerts, &bBox, &mesh](const int i) {
-               int vert = openVerts[i];
+  for_each_n(
+      policy, countAt(0), numOpenVert,
+      [&vertMorton, &vertBox, &openVerts, &bBox, &mesh, epsilon](const int i) {
+        int vert = openVerts[i];
 
-               const vec3 center(mesh.vertProperties[mesh.numProp * vert],
-                                 mesh.vertProperties[mesh.numProp * vert + 1],
-                                 mesh.vertProperties[mesh.numProp * vert + 2]);
+        const vec3 center(mesh.vertProperties[mesh.numProp * vert],
+                          mesh.vertProperties[mesh.numProp * vert + 1],
+                          mesh.vertProperties[mesh.numProp * vert + 2]);
 
-               vertBox[i].min = center - mesh.precision / 2.0;
-               vertBox[i].max = center + mesh.precision / 2.0;
+        vertBox[i].min = center - epsilon / 2.0;
+        vertBox[i].max = center + epsilon / 2.0;
 
-               vertMorton[i] = MortonCode(center, bBox);
-             });
+        vertMorton[i] = MortonCode(center, bBox);
+      });
 
   Vec<int> vertNew2Old(numOpenVert);
   sequence(vertNew2Old.begin(), vertNew2Old.end());
@@ -212,7 +212,7 @@ void Manifold::Impl::Finish() {
   if (halfedge_.size() == 0) return;
 
   CalculateBBox();
-  SetPrecision(precision_);
+  SetEpsilon(epsilon_);
   if (!bBox_.IsFinite()) {
     // Decimated out of existence - early out.
     MarkFailure(Error::NoError);
diff --git a/test/boolean_complex_test.cpp b/test/boolean_complex_test.cpp
index 61aea86f4..8aaaf00c5 100644
--- a/test/boolean_complex_test.cpp
+++ b/test/boolean_complex_test.cpp
@@ -236,7 +236,7 @@ TEST(BooleanComplex, Close) {
   Manifold result = a;
   for (int i = 0; i < 10; i++) {
     // std::cout << i << std::endl;
-    result ^= a.Translate({a.Precision() / 10 * i, 0.0, 0.0});
+    result ^= a.Translate({a.GetEpsilon() / 10 * i, 0.0, 0.0});
   }
   auto prop = result.GetProperties();
   const double tol = 0.004;
diff --git a/test/properties_test.cpp b/test/properties_test.cpp
index ada1a2efe..0ce338122 100644
--- a/test/properties_test.cpp
+++ b/test/properties_test.cpp
@@ -37,17 +37,17 @@ TEST(Properties, GetProperties) {
 
 TEST(Properties, Precision) {
   Manifold cube = Manifold::Cube();
-  EXPECT_FLOAT_EQ(cube.Precision(), kTolerance);
+  EXPECT_FLOAT_EQ(cube.GetEpsilon(), kTolerance);
   cube = cube.Scale({0.1, 1, 10});
-  EXPECT_FLOAT_EQ(cube.Precision(), 10 * kTolerance);
+  EXPECT_FLOAT_EQ(cube.GetEpsilon(), 10 * kTolerance);
   cube = cube.Translate({-100, -10, -1});
-  EXPECT_FLOAT_EQ(cube.Precision(), 100 * kTolerance);
+  EXPECT_FLOAT_EQ(cube.GetEpsilon(), 100 * kTolerance);
 }
 
 TEST(Properties, Precision2) {
   Manifold cube = Manifold::Cube();
   cube = cube.Translate({-0.5, 0, 0}).Scale({2, 1, 1});
-  EXPECT_FLOAT_EQ(cube.Precision(), 2 * kTolerance);
+  EXPECT_FLOAT_EQ(cube.GetEpsilon(), 2 * kTolerance);
 }
 
 TEST(Properties, Precision3) {
@@ -55,9 +55,8 @@ TEST(Properties, Precision3) {
   const auto prop = cylinder.GetProperties();
 
   MeshGL mesh = cylinder.GetMeshGL();
-  mesh.precision = 0.001;
   mesh.faceID.clear();
-  Manifold cylinder2(mesh);
+  Manifold cylinder2 = Manifold(mesh).SetEpsilon(0.001);
 
   const auto prop2 = cylinder2.GetProperties();
   EXPECT_NEAR(prop.volume, prop2.volume, 0.001);
diff --git a/test/samples_test.cpp b/test/samples_test.cpp
index dc1dd6548..f2c2d2ec8 100644
--- a/test/samples_test.cpp
+++ b/test/samples_test.cpp
@@ -221,7 +221,7 @@ TEST(Samples, GyroidModule) {
   CheckGL(gyroid);
 
   const Box bounds = gyroid.BoundingBox();
-  const double precision = gyroid.Precision();
+  const double precision = gyroid.GetEpsilon();
   EXPECT_NEAR(bounds.min.z, 0, precision);
   EXPECT_NEAR(bounds.max.z, size * std::sqrt(2.0), precision);
 
@@ -267,7 +267,7 @@ TEST(Samples, Sponge4) {
   EXPECT_EQ(cutSponge.second.Genus(), 13394);
 
   CrossSection projection(cutSponge.first.Project());
-  projection = projection.Simplify(cutSponge.first.Precision());
+  projection = projection.Simplify(cutSponge.first.GetEpsilon());
   Rect rect = projection.Bounds();
   Box box = cutSponge.first.BoundingBox();
   EXPECT_EQ(rect.min.x, box.min.x);
diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp
index f5fe4a1b7..21ebcfae7 100644
--- a/test/sdf_test.cpp
+++ b/test/sdf_test.cpp
@@ -70,7 +70,7 @@ TEST(SDF, Bounds) {
   Manifold cubeVoid = Manifold::LevelSet(
       CubeVoid(), {vec3(-size / 2), vec3(size / 2)}, edgeLength);
   Box bounds = cubeVoid.BoundingBox();
-  const double precision = cubeVoid.Precision();
+  const double precision = cubeVoid.GetEpsilon();
 #ifdef MANIFOLD_EXPORT
   if (options.exportModels)
     ExportMesh("cubeVoid.glb", cubeVoid.GetMeshGL(), {});
@@ -94,7 +94,7 @@ TEST(SDF, Bounds2) {
   Manifold cubeVoid = Manifold::LevelSet(
       CubeVoid(), {vec3(-size / 2), vec3(size / 2)}, edgeLength);
   Box bounds = cubeVoid.BoundingBox();
-  const double precision = cubeVoid.Precision();
+  const double precision = cubeVoid.GetEpsilon();
 #ifdef MANIFOLD_EXPORT
   if (options.exportModels)
     ExportMesh("cubeVoid2.glb", cubeVoid.GetMeshGL(), {});
@@ -121,7 +121,7 @@ TEST(SDF, Surface) {
   Manifold cube = Manifold::Cube(vec3(size), true);
   cube -= cubeVoid;
   Box bounds = cube.BoundingBox();
-  const double precision = cube.Precision();
+  const double precision = cube.GetEpsilon();
 #ifdef MANIFOLD_EXPORT
   if (options.exportModels) ExportMesh("cube.gltf", cube.GetMeshGL(), {});
 #endif
diff --git a/test/test_main.cpp b/test/test_main.cpp
index 13c45d013..cda4be654 100644
--- a/test/test_main.cpp
+++ b/test/test_main.cpp
@@ -299,6 +299,11 @@ void RelatedGL(const Manifold& out, const std::vector<MeshGL>& originals,
   ASSERT_FALSE(out.IsEmpty());
   const ivec3 normalIdx = updateNormals ? ivec3(3, 4, 5) : ivec3(0);
   MeshGL output = out.GetMeshGL(normalIdx);
+
+  float epsilon = std::max(static_cast<float>(out.GetEpsilon()),
+                           std::numeric_limits<float>::epsilon() *
+                               static_cast<float>(out.BoundingBox().Scale()));
+
   for (size_t run = 0; run < output.runOriginalID.size(); ++run) {
     const float* m = output.runTransform.data() + 12 * run;
     const mat3x4 transform =
@@ -350,7 +355,7 @@ void RelatedGL(const Manifold& out, const std::vector<MeshGL>& originals,
         vec3 edges[3];
         for (int k : {0, 1, 2}) edges[k] = inTriPos[k] - outTriPos[j];
         const double volume = la::dot(edges[0], la::cross(edges[1], edges[2]));
-        ASSERT_LE(volume, area * output.precision);
+        ASSERT_LE(volume, area * epsilon);
 
         if (checkNormals) {
           vec3 normal;
@@ -373,7 +378,7 @@ void RelatedGL(const Manifold& out, const std::vector<MeshGL>& originals,
             const double volumeP =
                 la::dot(edgesP[0], la::cross(edgesP[1], edgesP[2]));
 
-            ASSERT_LE(volumeP, area * output.precision);
+            ASSERT_LE(volumeP, area * epsilon);
           }
         }
       }