Fix CalculateNormals (#1724)

elalish · web-flow · commit 8cab4b607c02 · 2026-05-19T23:16:36.000+12:00
* no flat faces for setNormals

* make SmoothOut self-consistent
diff --git a/include/manifold/manifold.h b/include/manifold/manifold.h
@@ -245,7 +245,8 @@ class Manifold {
       int numProp,
       std::function<void(double*, vec3, const double*)> propFunc) const;
   Manifold CalculateCurvature(int gaussianIdx, int meanIdx) const;
-  Manifold CalculateNormals(int normalIdx = 0, double minSharpAngle = 60) const;
+  Manifold CalculateNormals(int normalIdx = 0,
+                            double minSharpAngle = 52.5) const;
   ///@}
 
   /** @name Smoothing
@@ -257,7 +258,8 @@ class Manifold {
   Manifold RefineToLength(double) const;
   Manifold RefineToTolerance(double) const;
   Manifold SmoothByNormals(int normalIdx = 0) const;
-  Manifold SmoothOut(double minSharpAngle = 60, double minSmoothness = 0) const;
+  Manifold SmoothOut(double minSharpAngle = 52.5,
+                     double minSmoothness = 0) const;
   static Manifold Smooth(const MeshGL&,
                          const std::vector<Smoothness>& sharpenedEdges = {});
   static Manifold Smooth(const MeshGL64&,
diff --git a/include/manifold/mesh.h b/include/manifold/mesh.h
@@ -21,6 +21,10 @@
 
 namespace manifold {
 
+/** @addtogroup Core
+ *  @{
+ */
+
 /**
  * @brief Mesh input/output suitable for pushing directly into graphics
  * libraries.
@@ -268,5 +272,5 @@ using MeshGL64 = MeshGLP<double, uint64_t>;
 MeshGL64 ReadOBJ(std::istream& stream);
 bool WriteOBJ(std::ostream& stream, const MeshGL64& mesh);
 #endif
-
+/** @} */
 }  // namespace manifold
diff --git a/src/manifold.cpp b/src/manifold.cpp
@@ -680,7 +680,7 @@ Manifold Manifold::CalculateCurvature(int gaussianIdx, int meanIdx) const {
 
 /**
  * Fills in vertex properties for normal vectors, calculated from the mesh
- * geometry. Flat faces composed of three or more triangles will remain flat.
+ * geometry.
  *
  * @param normalIdx The property channel in which to store the X values of the
  * normals. The X, Y, and Z channels will be sequential. The property set will
@@ -718,11 +718,9 @@ Manifold Manifold::CalculateNormals(int normalIdx, double minSharpAngle) const {
  * Smooths out the Manifold by filling in the halfedgeTangent vectors. The
  * geometry will remain unchanged until Refine, RefineToLength, or
  * RefineToTolerance is called to interpolate the surface. This version uses the
- * supplied vertex normal properties to define the tangent vectors. Faces of two
- * coplanar triangles will be marked as quads, while faces with three or more
- * will be flat. Zero-length normals are considered missing and will defer to
- * their neighboring normals instead. If all normals are missing, the vertex
- * pseudonormal will be used.
+ * supplied vertex normal properties to define the tangent vectors. Zero-length
+ * normals are considered missing and will defer to their neighboring normals
+ * instead. If all normals are missing, the vertex pseudonormal will be used.
  *
  * @param normalIdx The first property channel of the normals. NumProp must be
  * at least normalIdx + 3. Any vertex where multiple normals exist and don't
@@ -746,13 +744,13 @@ Manifold Manifold::SmoothByNormals(int normalIdx) const {
  * geometry will remain unchanged until Refine, RefineToLength, or
  * RefineToTolerance is called to interpolate the surface. This version uses the
  * geometry of the triangles and pseudo-normals to define the tangent vectors.
- * Faces of two coplanar triangles will be marked as quads.
+ * Faces of two coplanar triangles will be marked as quads, while faces with
+ * three or more will be flat.
  *
- * @param minSharpAngle degrees, default 60. Any edges with angles greater than
- * this value will remain sharp. The rest will be smoothed to G1 continuity,
- * with the caveat that flat faces of three or more triangles will always remain
- * flat. With a value of zero, the model is faceted, but in this case there is
- * no point in smoothing.
+ * @param minSharpAngle degrees, default 52.5. Any edges with angles greater
+ * than this value will remain sharp. The rest will be smoothed to G1
+ * continuity. With a value of zero, the model is faceted, but in this case
+ * there is no point in smoothing.
  *
  * @param minSmoothness range: 0 - 1, default 0. The smoothness applied to sharp
  * angles. The default gives a hard edge, while values > 0 will give a small
@@ -765,19 +763,7 @@ Manifold Manifold::SmoothOut(double minSharpAngle, double minSmoothness) const {
     return PropagateStatus(leafImpl->status_);
   auto pImpl = std::make_shared<Impl>(*leafImpl);
   if (!IsEmpty()) {
-    if (minSmoothness == 0) {
-      const int numProp = pImpl->numProp_;
-      Vec<double> properties = pImpl->properties_;
-      Halfedges halfedge(pImpl->halfedge_.ToData());
-      pImpl->SetNormals(0, minSharpAngle);
-      pImpl->CreateTangents(0);
-      // Reset the properties to the original values, removing temporary normals
-      pImpl->numProp_ = numProp;
-      pImpl->properties_ = std::move(properties);
-      pImpl->halfedge_ = std::move(halfedge);
-    } else {
-      pImpl->CreateTangents(pImpl->SharpenEdges(minSharpAngle, minSmoothness));
-    }
+    pImpl->CreateTangents(pImpl->SharpenEdges(minSharpAngle, minSmoothness));
   }
   return Manifold(std::make_shared<CsgLeafNode>(pImpl));
 }
diff --git a/src/smoothing.cpp b/src/smoothing.cpp
@@ -466,8 +466,6 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
 
   const int oldNumProp = NumProp();
 
-  Vec<bool> triIsFlatFace = FlatFaces();
-  Vec<int> vertFlatFace = VertFlatFace(triIsFlatFace);
   Vec<int> vertNumSharp(NumVert(), 0);
   for (size_t e = 0; e < halfedge_.size(); ++e) {
     if (!halfedge_.IsForward(e)) continue;
@@ -479,17 +477,6 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
     if (dihedral > minSharpAngle) {
       ++vertNumSharp[halfedge_.Start(e)];
       ++vertNumSharp[halfedge_.End(e)];
-    } else {
-      const bool faceSplit =
-          triIsFlatFace[tri1] != triIsFlatFace[tri2] ||
-          (triIsFlatFace[tri1] && triIsFlatFace[tri2] &&
-           !meshRelation_.triRef[tri1].SameFace(meshRelation_.triRef[tri2]));
-      if (vertFlatFace[halfedge_.Start(e)] == -2 && faceSplit) {
-        ++vertNumSharp[halfedge_.Start(e)];
-      }
-      if (vertFlatFace[halfedge_.End(e)] == -2 && faceSplit) {
-        ++vertNumSharp[halfedge_.End(e)];
-      }
     }
   }
 
@@ -525,9 +512,7 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
     const int vert = halfedge_.Start(startEdge);
 
     if (vertNumSharp[vert] < 2) {  // vertex has single normal
-      const vec3 worldNormal = vertFlatFace[vert] >= 0
-                                   ? faceNormal_[vertFlatFace[vert]]
-                                   : vertNormal_[vert];
+      const vec3 worldNormal = vertNormal_[vert];
       // Non-zero normalIdx is the legacy deferred-transform path: store in
       // per-mesh frame so GetMeshGL's runTransform application on export
       // recovers world frame even after later transforms. Standard slot 0
@@ -574,11 +559,7 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
 
       const double dihedral =
           degrees(AngleBetween(faceNormal_[face], faceNormal_[prevFace]));
-      if (dihedral > minSharpAngle ||
-          triIsFlatFace[face] != triIsFlatFace[prevFace] ||
-          (triIsFlatFace[face] && triIsFlatFace[prevFace] &&
-           !meshRelation_.triRef[face].SameFace(
-               meshRelation_.triRef[prevFace]))) {
+      if (dihedral > minSharpAngle) {
         break;
       }
       current = next;
@@ -596,44 +577,21 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
     ForVert<FaceEdge>(
         endEdge,
         [&](int current) {
-          if (IsInsideQuad(current)) {
-            return FaceEdge({current / 3, vec3(NAN)});
-          }
           const int vert = halfedge_.End(current);
-          vec3 pos = vertPos_[vert];
-          if (vertNumSharp[vert] < 2) {
-            // opposite vert has fixed normal
-            const vec3 normal = vertFlatFace[vert] >= 0
-                                    ? faceNormal_[vertFlatFace[vert]]
-                                    : vertNormal_[vert];
-            // Flair out the normal we're calculating to give the edge a
-            // more constant curvature to meet the opposite normal. Achieve
-            // this by pointing the tangent toward the opposite bezier
-            // control point instead of the vert itself.
-            pos += vec3(TangentFromNormal(normal, halfedge_.Pair(current)));
-          }
-          return FaceEdge({current / 3, SafeNormalize(pos - centerPos)});
+          return FaceEdge(
+              {current / 3, SafeNormalize(vertPos_[vert] - centerPos)});
         },
         [&](int, const FaceEdge& here, FaceEdge& next) {
           const double dihedral = degrees(
               AngleBetween(faceNormal_[here.face], faceNormal_[next.face]));
-          if (dihedral > minSharpAngle ||
-              triIsFlatFace[here.face] != triIsFlatFace[next.face] ||
-              (triIsFlatFace[here.face] && triIsFlatFace[next.face] &&
-               !meshRelation_.triRef[here.face].SameFace(
-                   meshRelation_.triRef[next.face]))) {
+          if (dihedral > minSharpAngle) {
             normals.push_back(vec3(0.0));
             meshIds.push_back(meshRelation_.triRef[next.face].meshID);
           }
           groups.push_back(normals.size() - 1);
           if (std::isfinite(next.normalizedEdge.x)) {
-            // Boolean subtractee triangles keep their original winding but
-            // have faceNormal_ negated, so the edge-cross points opposite to
-            // faceNormal_ there - flip to keep the accumulation outward-from-
-            // result.
             vec3 dir = SafeNormalize(
                 la::cross(next.normalizedEdge, here.normalizedEdge));
-            if (la::dot(dir, faceNormal_[next.face]) < 0) dir = -dir;
             normals.back() +=
                 dir * AngleBetween(here.normalizedEdge, next.normalizedEdge);
           } else {
@@ -642,10 +600,10 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
         });
 
     for (int i = 0; i < normals.size(); ++i) {
-      vec3 n = SafeNormalize(normals[i]);
+      vec3 n = normals[i];
       // Same frame-storage rule as the single-normal path above.
       if (normalIdx != 0) n = getTransform(meshIds[i]) * n;
-      normals[i] = n;
+      normals[i] = SafeNormalize(n);
     }
 
     int lastGroup = 0;
diff --git a/test/properties_test.cpp b/test/properties_test.cpp
@@ -121,13 +121,14 @@ TEST(Properties, CalculateCurvature) {
 
 TEST(Properties, CalculateNormals) {
   Manifold sphere = Manifold::Sphere(10);
-  Manifold cut = (sphere - sphere.Translate(vec3(10)));
+  vec3 center(10);
+  Manifold cut = (sphere - sphere.Translate(center));
   cut.Status();
   Manifold cut2(cut.GetMeshGL64());
   EXPECT_TRUE(cut.MatchesTriNormals());
   EXPECT_TRUE(cut2.MatchesTriNormals());
-  MeshGL64 out = cut.CalculateNormals(0).GetMeshGL64(0);
-  MeshGL64 out2 = cut2.CalculateNormals(0).GetMeshGL64(0);
+  MeshGL64 out = cut.CalculateNormals().GetMeshGL64();
+  MeshGL64 out2 = cut2.CalculateNormals().GetMeshGL64();
   ASSERT_EQ(out.NumTri(), out2.NumTri());
   ASSERT_EQ(out.NumVert(), out2.NumVert());
   ASSERT_EQ(out.numProp, out2.numProp);
@@ -143,12 +144,26 @@ TEST(Properties, CalculateNormals) {
       norm[j] = out.vertProperties[np * v + 3 + j];
       norm2[j] = out2.vertProperties[np * v + 3 + j];
       ASSERT_FLOAT_EQ(pos[j], pos2[j]);
+      ASSERT_NEAR(norm[j], norm2[j], 1e-14);
     }
     if (dot(pos, norm) <= 0) ++numBad;
     if (dot(pos2, norm2) <= 0) ++numBad2;
+    EXPECT_FLOAT_EQ(la::length(norm), 1.0);
+    EXPECT_TRUE(dot(la::normalize(pos), norm) > 0.99 ||
+                dot(la::normalize(center - pos), norm) > 0.99);
   }
   EXPECT_EQ(numBad, 0);
   EXPECT_EQ(numBad2, 0);
+  for (int tv = out2.runIndex[0]; tv < out2.runIndex[1]; ++tv) {
+    const int v = out2.triVerts[tv];
+    auto pos = out2.GetVertPos(v);
+    auto norm = pos;
+    for (int j : {0, 1, 2}) {
+      norm[j] = out2.vertProperties[np * v + 3 + j];
+    }
+    EXPECT_FLOAT_EQ(la::length(norm), 1.0);
+    EXPECT_GT(dot(la::normalize(pos), norm), 0.99);
+  }
 }
 
 TEST(Properties, Coplanar) {
diff --git a/test/sdf_test.cpp b/test/sdf_test.cpp
@@ -164,22 +164,19 @@ TEST(SDF, Resize) {
 }
 
 TEST(SDF, SineSurface) {
-  Manifold surface =
-      Manifold::LevelSet(
-          [](vec3 p) {
-            double mid = la::sin(p.x) + la::sin(p.y);
-            return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0f : -1.0f;
-          },
-          {vec3(-1.75 * kPi), vec3(1.75 * kPi)}, 1)
-          .Simplify();
-  Manifold smoothed = surface.SmoothOut(180).RefineToLength(0.05);
-
-  EXPECT_EQ(smoothed.Status(), Manifold::Error::NoError);
-  EXPECT_EQ(smoothed.Genus(), 38);
-  EXPECT_NEAR(smoothed.Volume(), 107.1, 0.1);
-  EXPECT_NEAR(smoothed.SurfaceArea(), 394.0, 0.1);
-
-  if (options.exportModels) WriteTestOBJ("sinesurface.obj", smoothed);
+  Manifold surface = Manifold::LevelSet(
+      [](vec3 p) {
+        double mid = la::sin(p.x) + la::sin(p.y);
+        return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0f : -1.0f;
+      },
+      {vec3(-1.75 * kPi), vec3(1.75 * kPi)}, 1);
+
+  EXPECT_EQ(surface.Status(), Manifold::Error::NoError);
+  EXPECT_EQ(surface.Genus(), 38);
+  EXPECT_NEAR(surface.Volume(), 102.4, 0.1);
+  EXPECT_NEAR(surface.SurfaceArea(), 392.4, 0.1);
+
+  if (options.exportModels) WriteTestOBJ("sinesurface.obj", surface);
 }
 
 TEST(SDF, Blobs) {
diff --git a/test/smooth_test.cpp b/test/smooth_test.cpp
@@ -69,8 +69,8 @@ TEST(Smooth, RefineQuads) {
 TEST(Smooth, TruncatedCone) {
   Manifold cone = Manifold::Cylinder(5, 10, 5, 12);
   Manifold smooth = cone.SmoothOut().RefineToLength(0.5).CalculateNormals(0);
-  EXPECT_NEAR(smooth.Volume(), 1158.61, 0.01);
-  EXPECT_NEAR(smooth.SurfaceArea(), 768.12, 0.01);
+  EXPECT_NEAR(smooth.Volume(), 1163.53, 0.01);
+  EXPECT_NEAR(smooth.SurfaceArea(), 769.33, 0.01);
   CheckGL(smooth, false);
 
   Manifold smooth1 = cone.SmoothOut(180, 1).RefineToLength(0.5);
@@ -90,16 +90,16 @@ TEST(Smooth, ToLength) {
   Manifold smooth =
       cone.AsOriginal().Simplify().SmoothOut(180).RefineToLength(0.1);
   ExpectMeshes(smooth, {{85250, 170496}});
-  EXPECT_NEAR(smooth.Volume(), 4505, 1);
-  EXPECT_NEAR(smooth.SurfaceArea(), 1337, 1);
+  EXPECT_NEAR(smooth.Volume(), 4570, 1);
+  EXPECT_NEAR(smooth.SurfaceArea(), 1348, 1);
 
   MeshGL out = smooth.CalculateCurvature(-1, 0).GetMeshGL();
   float maxMeanCurvature = 0;
   for (size_t i = 3; i < out.vertProperties.size(); i += 4) {
     maxMeanCurvature =
         std::max(maxMeanCurvature, std::abs(out.vertProperties[i]));
   }
-  EXPECT_NEAR(maxMeanCurvature, 2.32, 0.01);
+  EXPECT_NEAR(maxMeanCurvature, 1.63, 0.01);
 
   if (options.exportModels) WriteTestOBJ("smoothToLength.obj", smooth);
 }
@@ -213,11 +213,11 @@ TEST(Smooth, MissingNormals) {
 }
 
 TEST(Smooth, MissingNormalsCone) {
-  Manifold cone = Manifold::Cylinder(10, 10, 0, 5).CalculateNormals(0);
+  Manifold cone = Manifold::Cylinder(10, 10, 0, 5).CalculateNormals(0, 60);
   Manifold diff = cone - Manifold::Cube(vec3(10), true).Translate({0, 0, 10});
   Manifold out = diff.SmoothByNormals(0).Refine(20);
-  EXPECT_NEAR(out.Volume(), 1092, 1);
-  EXPECT_NEAR(out.SurfaceArea(), 748, 1);
+  EXPECT_NEAR(out.Volume(), 1009, 1);
+  EXPECT_NEAR(out.SurfaceArea(), 736, 1);
   if (options.exportModels) WriteTestOBJ("missingNormalsCone.obj", out);
 }
 
@@ -383,48 +383,6 @@ TEST(Smooth, Torus) {
 }
 #endif
 
-TEST(Smooth, SineSurface) {
-  Manifold surface =
-      Manifold::LevelSet(
-          [](vec3 p) {
-            double mid = la::sin(p.x) + la::sin(p.y);
-            return (p.z > mid - 0.5 && p.z < mid + 0.5) ? 1.0 : -1.0;
-          },
-          {vec3(-2 * kPi + 0.2), vec3(0 * kPi - 0.2)}, 1)
-          .Simplify();
-
-  Manifold smoothed =
-      surface.CalculateNormals(0, 50).SmoothByNormals(0).Refine(8);
-  EXPECT_NEAR(smoothed.Volume(), 8.07, 0.01);
-  EXPECT_NEAR(smoothed.SurfaceArea(), 30.87, 0.01);
-  EXPECT_EQ(smoothed.Genus(), 0);
-  EXPECT_NEAR(smoothed.TrimByPlane({0, 1, 1}, -3.19487).Volume(),
-              smoothed.Volume(), 1e-5);
-
-  Manifold smoothed1 = surface.SmoothOut(50).Refine(8);
-  EXPECT_FLOAT_EQ(smoothed1.Volume(), smoothed.Volume());
-  EXPECT_FLOAT_EQ(smoothed1.SurfaceArea(), smoothed.SurfaceArea());
-  EXPECT_EQ(smoothed1.Genus(), 0);
-  EXPECT_NEAR(smoothed1.TrimByPlane({0, 1, 1}, -3.19487).Volume(),
-              smoothed1.Volume(), 1e-5);
-
-  Manifold smoothed2 = surface.SmoothOut(180, 1).Refine(8);
-  EXPECT_NEAR(smoothed2.Volume(), 8.95, 0.01);
-  EXPECT_NEAR(smoothed2.SurfaceArea(), 33.35, 0.01);
-  EXPECT_EQ(smoothed2.Genus(), 0);
-  EXPECT_NEAR(smoothed2.TrimByPlane({0, 1, 1}, -3.19487).Volume(),
-              smoothed2.Volume(), 1e-3);
-
-  Manifold smoothed3 = surface.SmoothOut(50, 0.5).Refine(8);
-  EXPECT_NEAR(smoothed3.Volume(), 8.38, 0.01);
-  EXPECT_NEAR(smoothed3.SurfaceArea(), 31.55, 0.02);
-  EXPECT_EQ(smoothed3.Genus(), 0);
-  EXPECT_NEAR(smoothed3.TrimByPlane({0, 1, 1}, -3.19487).Volume(),
-              smoothed3.Volume(), 1e-5);
-
-  if (options.exportModels) WriteTestOBJ("smoothSineSurface.obj", smoothed);
-}
-
 TEST(Smooth, SDF) {
   const double r = 10;
   const double extra = 2;
