Smoothing fixes (#796)

* fixed nan normals

* re-identify flat faces after refine

* fix color interpolation

* fix build

Author: elalish

src/manifold/src/impl.h

@@ -178,6 +178,7 @@ struct Manifold::Impl {
   std::vector<Smoothness> SharpenEdges(float minSharpAngle,
                                        float minSmoothness) const;
   void SetNormals(int normalIdx, float minSharpAngle);
+  void LinearizeFlatTangents();
   void CreateTangents(int normalIdx);
   void CreateTangents(std::vector<Smoothness>);
   Vec<Barycentric> Subdivide(std::function<int(glm::vec3)>);

src/manifold/src/smoothing.cpp

@@ -84,9 +84,8 @@ struct SmoothBezier {
     const glm::vec3 edgeVec =
         impl->vertPos_[edge.endVert] - impl->vertPos_[edge.startVert];
     const glm::vec3 edgeNormal =
-        (impl->faceNormal_[edge.face] +
-         impl->faceNormal_[impl->halfedge_[edge.pairedHalfedge].face]) /
-        2.0f;
+        impl->faceNormal_[edge.face] +
+        impl->faceNormal_[impl->halfedge_[edge.pairedHalfedge].face];
     glm::vec3 dir =
         glm::cross(glm::cross(edgeNormal, edgeVec), vertNormal[edge.startVert]);
     tangent = CircularTangent(dir, edgeVec);
@@ -966,7 +965,7 @@ void Manifold::Impl::SetNormals(int normalIdx, float minSharpAngle) {
               group.push_back(normals.size() - 1);
               float dot = glm::dot(here.edgeVec, next.edgeVec);
               const float phi =
-                  dot >= 1 ? 0
+                  dot >= 1 ? kTolerance
                            : (dot <= -1 ? glm::pi<float>() : glm::acos(dot));
               normals.back() += faceNormal_[next.face] * phi;
             });
@@ -1014,6 +1013,35 @@ void Manifold::Impl::SetNormals(int normalIdx, float minSharpAngle) {
   }
 }
 
+/**
+ * Tangents get flattened to create sharp edges by setting their weight to zero.
+ * This is the natural limit of reducing the weight to increase the sharpness
+ * smoothly. This limit gives a decent shape, but it causes the parameterization
+ * to be stretched and compresses it near the edges, which is good for resolving
+ * tight curvature, but bad for property interpolation. This function fixes the
+ * parameter stretch at the limit for sharp edges, since there is no curvature
+ * to resolve. Note this also changes the overall shape - making it more evenly
+ * curved.
+ */
+void Manifold::Impl::LinearizeFlatTangents() {
+  const int n = halfedgeTangent_.size();
+  for_each_n(autoPolicy(n), zip(halfedgeTangent_.begin(), countAt(0)), n,
+             [this](thrust::tuple<glm::vec4&, int> inOut) {
+               glm::vec4& tangent = thrust::get<0>(inOut);
+               const int halfedge = thrust::get<1>(inOut);
+               if (tangent.w != 0) {
+                 return;
+               }
+
+               const glm::vec4 otherTangent =
+                   halfedgeTangent_[halfedge_[halfedge].pairedHalfedge];
+               glm::vec3 edge = vertPos_[halfedge_[halfedge].endVert] +
+                                glm::vec3(otherTangent) -
+                                vertPos_[halfedge_[halfedge].startVert];
+               tangent = glm::vec4(edge / 3.0f, 1);
+             });
+}
+
 /**
  * Calculates halfedgeTangent_, allowing the manifold to be refined and
  * smoothed. The tangents form weighted cubic Beziers along each edge. This
@@ -1094,6 +1122,7 @@ void Manifold::Impl::CreateTangents(int normalIdx) {
       });
     }
   }
+  LinearizeFlatTangents();
 }
 
 /**
@@ -1207,6 +1236,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
       });
     }
   }
+  LinearizeFlatTangents();
 }
 
 /**
@@ -1517,6 +1547,7 @@ void Manifold::Impl::Refine(std::function<int(glm::vec3)> edgeDivisions) {
     // being non-coplanar, and hence not being related to the original faces.
     meshRelation_.originalID = ReserveIDs(1);
     InitializeOriginal();
+    CreateFaces();
   }
 
   halfedgeTangent_.resize(0);

test/manifold_test.cpp

@@ -357,30 +357,44 @@ TEST(Manifold, RefineQuads) {
   EXPECT_EQ(cylinder.NumTri(), 16892);
   auto prop = cylinder.GetProperties();
   EXPECT_NEAR(prop.volume, 2 * glm::pi<float>(), 0.003);
-  EXPECT_NEAR(prop.surfaceArea, 6 * glm::pi<float>(), 0.0003);
+  EXPECT_NEAR(prop.surfaceArea, 6 * glm::pi<float>(), 0.004);
+  const MeshGL out = cylinder.GetMeshGL();
+  CheckGL(out);
+
+  const MeshGL baseline = WithPositionColors(cylinder);
+  float maxDiff = 0;
+  for (int i = 0; i < out.vertProperties.size(); ++i) {
+    maxDiff = glm::max(
+        maxDiff, glm::abs(out.vertProperties[i] - baseline.vertProperties[i]));
+  }
+  // This has a wide tolerance because the triangle colors on the ends are still
+  // being stretched out into circular arcs, which introduces unavoidable error.
+  EXPECT_LE(maxDiff, 0.07);
 
 #ifdef MANIFOLD_EXPORT
   ExportOptions options2;
+  options2.mat.metalness = 0;
+  options2.mat.roughness = 0.5;
   options2.mat.colorChannels = {3, 4, 5, -1};
-  if (options.exportModels)
-    ExportMesh("refinedCylinder.glb", cylinder.GetMeshGL(), options2);
+  if (options.exportModels) ExportMesh("refinedCylinder.glb", out, options2);
 #endif
 }
 
 TEST(Manifold, SmoothFlat) {
   Manifold cone = Manifold::Cylinder(5, 10, 5, 12).SmoothOut();
   Manifold smooth = cone.RefineToLength(0.1).CalculateNormals(0);
   auto prop = smooth.GetProperties();
-  EXPECT_NEAR(prop.volume, 1159.02, 0.01);
-  EXPECT_NEAR(prop.surfaceArea, 771.45, 0.01);
+  EXPECT_NEAR(prop.volume, 1142.9, 0.01);
+  EXPECT_NEAR(prop.surfaceArea, 764.28, 0.01);
+  MeshGL out = smooth.GetMeshGL();
+  CheckGL(out);
 
 #ifdef MANIFOLD_EXPORT
   ExportOptions options2;
   options2.faceted = false;
   options2.mat.normalChannels = {3, 4, 5};
   options2.mat.roughness = 0;
-  if (options.exportModels)
-    ExportMesh("smoothCone.glb", smooth.GetMeshGL(), options2);
+  if (options.exportModels) ExportMesh("smoothCone.glb", out, options2);
 #endif
 }
 

test/test.h

@@ -57,6 +57,7 @@ MeshGL WithPositionColors(const Manifold& in);
 MeshGL WithNormals(const Manifold& in);
 float GetMaxProperty(const MeshGL& mesh, int channel);
 float GetMinProperty(const MeshGL& mesh, int channel);
+void CheckFinite(const MeshGL& mesh);
 void Identical(const Mesh& mesh1, const Mesh& mesh2);
 void RelatedGL(const Manifold& out, const std::vector<MeshGL>& originals,
                bool checkNormals = false, bool updateNormals = false);

test/test_main.cpp

@@ -305,6 +305,18 @@ float GetMinProperty(const MeshGL& mesh, int channel) {
   return min;
 }
 
+void CheckFinite(const MeshGL& mesh) {
+  for (float v : mesh.vertProperties) {
+    ASSERT_TRUE(glm::isfinite(v));
+  }
+  for (float v : mesh.runTransform) {
+    ASSERT_TRUE(glm::isfinite(v));
+  }
+  for (float v : mesh.halfedgeTangent) {
+    ASSERT_TRUE(glm::isfinite(v));
+  }
+}
+
 void Identical(const Mesh& mesh1, const Mesh& mesh2) {
   ASSERT_EQ(mesh1.vertPos.size(), mesh2.vertPos.size());
   for (int i = 0; i < mesh1.vertPos.size(); ++i)
@@ -455,6 +467,7 @@ void CheckGL(const Manifold& manifold) {
     EXPECT_EQ(meshGL.runTransform.size(), 12 * meshGL.runOriginalID.size());
   }
   EXPECT_EQ(meshGL.faceID.size(), meshGL.NumTri());
+  CheckFinite(meshGL);
 }
 
 #ifdef MANIFOLD_EXPORT