added global coplanarity check (elalish#408)

* added global coplanarity check

* fix early-out

* added test, fixed precision

* one more precision

Author: elalish

.vscode/settings.json

@@ -120,6 +120,7 @@
   "editor.detectIndentation": false,
   "cSpell.words": [
     "Bary",
+    "Coplanarity",
     "csaszar",
     "Gyroid",
     "halfedge",

src/manifold/src/impl.cpp

@@ -274,6 +274,36 @@ struct CoplanarEdge {
   }
 };
 
+struct CheckCoplanarity {
+  int* comp2tri;
+  const Halfedge* halfedge;
+  const glm::vec3* vertPos;
+  const int* components;
+  const float precision;
+
+  __host__ __device__ void operator()(int tri) {
+    const int component = components[tri];
+    const int referenceTri = comp2tri[component];
+    if (referenceTri < 0 || referenceTri == tri) return;
+
+    const glm::vec3 origin = vertPos[halfedge[3 * referenceTri].startVert];
+    const glm::vec3 normal = glm::normalize(
+        glm::cross(vertPos[halfedge[3 * referenceTri + 1].startVert] - origin,
+                   vertPos[halfedge[3 * referenceTri + 2].startVert] - origin));
+
+    for (const int i : {0, 1, 2}) {
+      const glm::vec3 vert = vertPos[halfedge[3 * tri + i].startVert];
+      // 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 (glm::abs(glm::dot(normal, vert - origin)) > precision) {
+        comp2tri[component] = -1;
+        break;
+      }
+    }
+  }
+};
+
 struct EdgeBox {
   const glm::vec3* vertPos;
 
@@ -326,6 +356,7 @@ uint32_t Manifold::Impl::ReserveIDs(uint32_t n) {
 Manifold::Impl::Impl(const MeshGL& meshGL,
                      std::vector<float> propertyTolerance) {
   Mesh mesh;
+  mesh.precision = meshGL.precision;
   const int numVert = meshGL.NumVert();
   const int numTri = meshGL.NumTri();
 
@@ -472,7 +503,7 @@ Manifold::Impl::Impl(const Mesh& mesh, const MeshRelationD& relation,
     MarkFailure(Error::NonFiniteVertex);
     return;
   }
-  SetPrecision();
+  SetPrecision(mesh.precision);
 
   CreateHalfedges(triVerts);
   if (!IsManifold()) {
@@ -609,9 +640,20 @@ void Manifold::Impl::CreateFaces(const std::vector<float>& propertyTolerance) {
     }
   }
 
+  VecDH<int> componentsD(components);
+  VecDH<int> comp2triD(comp2tri);
+  for_each_n(
+      autoPolicy(halfedge_.size()), countAt(0), NumTri(),
+      CheckCoplanarity({comp2triD.ptrD(), halfedge_.cptrD(), vertPos_.cptrD(),
+                        componentsD.cptrD(), precision_}));
+
   VecDH<TriRef>& triRef = meshRelation_.triRef;
-  for (int tri = 0; tri < NumTri(); ++tri)
-    triRef[tri].tri = comp2tri[components[tri]];
+  for (int tri = 0; tri < NumTri(); ++tri) {
+    const int referenceTri = comp2triD[components[tri]];
+    if (referenceTri >= 0) {
+      triRef[tri].tri = referenceTri;
+    }
+  }
 }
 
 /**

src/manifold/src/manifold.cpp

@@ -137,6 +137,7 @@ Mesh Manifold::GetMesh() const {
   const Impl& impl = *GetCsgLeafNode().GetImpl();
 
   Mesh result;
+  result.precision = Precision();
   result.vertPos.insert(result.vertPos.end(), impl.vertPos_.begin(),
                         impl.vertPos_.end());
   result.vertNormal.insert(result.vertNormal.end(), impl.vertNormal_.begin(),
@@ -178,6 +179,7 @@ MeshGL Manifold::GetMeshGL(glm::ivec3 normalIdx) const {
       !isOriginal && glm::all(glm::greaterThan(normalIdx, glm::ivec3(2)));
 
   MeshGL out;
+  out.precision = Precision();
   out.numProp = 3 + numProp;
   out.triVerts.resize(3 * numTri);
 
@@ -379,10 +381,7 @@ int Manifold::Genus() const {
 }
 
 /**
- * Returns the surface area and volume of the manifold. These properties are
- * clamped to zero for a given face if they are within the Precision(). This
- * means degenerate manifolds can by identified by testing these properties as
- * == 0.
+ * Returns the surface area and volume of the manifold.
  */
 Properties Manifold::GetProperties() const {
   return GetCsgLeafNode().GetImpl()->GetProperties();

src/manifold/src/properties.cpp

@@ -43,9 +43,7 @@ struct FaceAreaVolume {
     float area = glm::length(crossP);
     float volume = glm::dot(crossP, vertPos[halfedges[3 * face].startVert]);
 
-    return area > perimeter * precision
-               ? thrust::make_pair(area / 2.0f, volume / 6.0f)
-               : thrust::make_pair(0.0f, 0.0f);
+    return thrust::make_pair(area / 2.0f, volume / 6.0f);
   }
 };
 

src/manifold/src/sort.cpp

@@ -480,6 +480,7 @@ void Manifold::Impl::GatherFaces(const Impl& old,
 /// Constructs a position-only MeshGL from the input Mesh.
 MeshGL::MeshGL(const Mesh& mesh) {
   numProp = 3;
+  precision = mesh.precision;
   vertProperties.resize(numProp * mesh.vertPos.size());
   for (int i = 0; i < mesh.vertPos.size(); ++i) {
     for (int j : {0, 1, 2}) vertProperties[3 * i + j] = mesh.vertPos[i][j];

src/utilities/include/public.h

@@ -145,6 +145,11 @@ struct Mesh {
   /// as 3 * tri + i, representing the tangent from Mesh.triVerts[tri][i] along
   /// the CCW edge. If empty, mesh is faceted.
   std::vector<glm::vec4> halfedgeTangent;
+  /// The absolute precision of the vertex positions, based on accrued rounding
+  /// errors. 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.
+  float precision = 0;
 };
 
 /**

test/cross_section_test.cpp

@@ -61,7 +61,7 @@ TEST(CrossSection, RoundOffset) {
 #endif
 
   EXPECT_EQ(result.Genus(), 0);
-  EXPECT_NEAR(result.GetProperties().volume, 4373, 1);
+  EXPECT_NEAR(result.GetProperties().volume, 4393, 1);
 }
 
 TEST(CrossSection, Empty) {

test/manifold_test.cpp

@@ -362,6 +362,20 @@ TEST(Manifold, Precision2) {
   EXPECT_FLOAT_EQ(cube.Precision(), 2 * kTolerance);
 }
 
+TEST(Manifold, Precision3) {
+  Manifold cylinder = Manifold::Cylinder(1, 1, 1, 1000);
+  const auto prop = cylinder.GetProperties();
+
+  MeshGL mesh = cylinder.GetMeshGL();
+  mesh.precision = 0.001;
+  mesh.faceID.clear();
+  Manifold cylinder2(mesh);
+
+  const auto prop2 = cylinder2.GetProperties();
+  EXPECT_NEAR(prop.volume, prop2.volume, 0.001);
+  EXPECT_NEAR(prop.surfaceArea, prop2.surfaceArea, 0.001);
+}
+
 /**
  * Curvature is the inverse of the radius of curvature, and signed such that
  * positive is convex and negative is concave. There are two orthogonal