CalculateNormals and Smooth (#771)

* started implementation

* update to Vec

* refactor

* snap normals to faces

* added flat face sharpening

* added Smooth member

* Smooth working

* refactor SetNormals

* calculateNormals runs

* CalculateNormals works

* update test

* add bindings

* fix python bindings

* fix format

* Smooth -> SmoothOut

* a few fixes

* weird format bug

* addressing feedback

* update docs

.github/workflows/check_format.yml

@@ -12,7 +12,7 @@ jobs:
 
     steps:
     - uses: actions/checkout@v4
-    - uses: DoozyX/clang-format-lint-action@v0.14
+    - uses: DoozyX/clang-format-lint-action@v0.17
       with:
         source: '.'
         exclude: '*/third_party'

.vscode/settings.json

@@ -120,6 +120,7 @@
   "cmake.configureArgs": [
     "-DBUILD_SHARED_LIBS=ON",
     "-DMANIFOLD_CBIND=ON",
+    "-DMANIFOLD_PYBIND=ON",
     "-DMANIFOLD_EXPORT=ON",
     "-DMANIFOLD_DEBUG=ON",
     "-DMANIFOLD_PAR=TBB",

bindings/c/include/manifoldc.h

@@ -139,6 +139,8 @@ ManifoldManifold *manifold_mirror(void *mem, ManifoldManifold *m, float nx,
                                   float ny, float nz);
 ManifoldManifold *manifold_warp(void *mem, ManifoldManifold *m,
                                 ManifoldVec3 (*fun)(float, float, float));
+ManifoldManifold *manifold_smooth_out(void *mem, ManifoldManifold *m,
+                                      float minSharpAngle, float minSmoothness);
 ManifoldManifold *manifold_refine(void *mem, ManifoldManifold *m, int refine);
 ManifoldManifold *manifold_refine_to_length(void *mem, ManifoldManifold *m,
                                             float length);
@@ -188,6 +190,9 @@ ManifoldManifold *manifold_set_properties(
     void (*fun)(float *new_prop, ManifoldVec3 position, const float *old_prop));
 ManifoldManifold *manifold_calculate_curvature(void *mem, ManifoldManifold *m,
                                                int gaussian_idx, int mean_idx);
+ManifoldManifold *manifold_calculate_normals(void *mem, ManifoldManifold *m,
+                                             int normal_idx,
+                                             int min_sharp_angle);
 
 // CrossSection Shapes/Constructors
 

bindings/c/manifoldc.cpp

@@ -300,6 +300,13 @@ ManifoldMeshGL *manifold_level_set_seq_context(
   return level_set_context(mem, sdf, bounds, edge_length, level, true, ctx);
 }
 
+ManifoldManifold *manifold_smooth_out(void *mem, ManifoldManifold *m,
+                                      float minSharpAngle,
+                                      float minSmoothness) {
+  auto smoothed = from_c(m)->SmoothOut(minSharpAngle, minSmoothness);
+  return to_c(new (mem) Manifold(smoothed));
+}
+
 ManifoldManifold *manifold_refine(void *mem, ManifoldManifold *m, int refine) {
   auto refined = from_c(m)->Refine(refine);
   return to_c(new (mem) Manifold(refined));
@@ -543,6 +550,13 @@ ManifoldManifold *manifold_calculate_curvature(void *mem, ManifoldManifold *m,
   return to_c(new (mem) Manifold(man));
 }
 
+ManifoldManifold *manifold_calculate_normals(void *mem, ManifoldManifold *m,
+                                             int normal_idx,
+                                             int min_sharp_angle) {
+  auto man = from_c(m)->CalculateNormals(normal_idx, min_sharp_angle);
+  return to_c(new (mem) Manifold(man));
+}
+
 // Static Quality Globals
 
 void manifold_set_min_circular_angle(float degrees) {

bindings/python/examples/all_apis.py

@@ -69,6 +69,7 @@ def all_manifold():
     m = Manifold.batch_hull([m, m.translate((0, 0, 1))])
     b = m.bounding_box()
     m = m.calculate_curvature(4, 5)
+    m = m.calculate_normals(0)
     m = Manifold.compose([m, m.translate((5, 0, 0))])
     m = Manifold.cube((1, 1, 1))
     m = Manifold.cylinder(1, 1)
@@ -92,6 +93,7 @@ def all_manifold():
     c = m.project()
     m = m.refine(2)
     m = m.refine_to_length(0.1)
+    m = m.smooth_out()
     i = Manifold.reserve_ids(1)
     m = m.scale((1, 2, 3))
     m = m.set_properties(3, lambda pos, prop: pos)

bindings/python/manifold3d.cpp

@@ -315,6 +315,12 @@ NB_MODULE(manifold3d, m) {
       .def("calculate_curvature", &Manifold::CalculateCurvature,
            nb::arg("gaussian_idx"), nb::arg("mean_idx"),
            manifold__calculate_curvature__gaussian_idx__mean_idx)
+      .def("calculate_normals", &Manifold::CalculateNormals,
+           nb::arg("normal_idx"), nb::arg("min_sharp_angle") = 60,
+           manifold__calculate_normals__normal_idx__min_sharp_angle)
+      .def("smooth_out", &Manifold::SmoothOut, nb::arg("min_sharp_angle") = 60,
+           nb::arg("min_smoothness") = 0,
+           manifold__smooth_out__min_sharp_angle__min_smoothness)
       .def("refine", &Manifold::Refine, nb::arg("n"), manifold__refine__n)
       .def("refine_to_length", &Manifold::RefineToLength, nb::arg("length"),
            manifold__refine_to_length__length)

bindings/wasm/bindings.cpp

@@ -144,6 +144,7 @@ EMSCRIPTEN_BINDINGS(whatever) {
       .function("_GetMeshJS", &js::GetMeshJS)
       .function("refine", &Manifold::Refine)
       .function("refineToLength", &Manifold::RefineToLength)
+      .function("smoothOut", &Manifold::SmoothOut)
       .function("_Warp", &man_js::Warp)
       .function("_SetProperties", &man_js::SetProperties)
       .function("transform", &man_js::Transform)
@@ -165,6 +166,7 @@ EMSCRIPTEN_BINDINGS(whatever) {
       .function("genus", &Manifold::Genus)
       .function("getProperties", &Manifold::GetProperties)
       .function("calculateCurvature", &Manifold::CalculateCurvature)
+      .function("calculateNormals", &Manifold::CalculateNormals)
       .function("originalID", &Manifold::OriginalID)
       .function("asOriginal", &Manifold::AsOriginal);
 

bindings/wasm/examples/worker.test.js

@@ -116,8 +116,8 @@ suite('Examples', () => {
   test('Scallop', async () => {
     const result = await runExample('Scallop');
     expect(result.genus).to.equal(0, 'Genus');
-    expect(result.volume).to.be.closeTo(41284, 1, 'Volume');
-    expect(result.surfaceArea).to.be.closeTo(7810, 1, 'Surface Area');
+    expect(result.volume).to.be.closeTo(41100, 100, 'Volume');
+    expect(result.surfaceArea).to.be.closeTo(7790, 10, 'Surface Area');
   });
 
   test('Torus Knot', async () => {

bindings/wasm/examples/worker.ts

@@ -59,10 +59,29 @@ const manifoldStaticFunctions = [
 ];
 // manifold member functions (that return a new manifold)
 const manifoldMemberFunctions = [
-  'add',       'subtract',       'intersect',     'decompose',  'warp',
-  'transform', 'translate',      'rotate',        'scale',      'mirror',
-  'refine',    'refineToLength', 'setProperties', 'asOriginal', 'trimByPlane',
-  'split',     'splitByPlane',   'slice',         'project',    'hull'
+  'add',
+  'subtract',
+  'intersect',
+  'decompose',
+  'warp',
+  'transform',
+  'translate',
+  'rotate',
+  'scale',
+  'mirror',
+  'calculateCurvature',
+  'calculateNormals',
+  'smoothOut',
+  'refine',
+  'refineToLength',
+  'setProperties',
+  'asOriginal',
+  'trimByPlane',
+  'split',
+  'splitByPlane',
+  'slice',
+  'project',
+  'hull'
 ];
 // CrossSection static methods (that return a new cross-section)
 const crossSectionStaticFunctions = [

bindings/wasm/manifold-encapsulated-types.d.ts

@@ -559,6 +559,24 @@ export class Manifold {
    */
   warp(warpFunc: (vert: Vec3) => void): Manifold;
 
+  /**
+   * Smooths out the Manifold by filling in the halfedgeTangent vectors. The
+   * geometry will remain unchanged until Refine or RefineToLength is called to
+   * interpolate the surface.
+   *
+   * @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 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 fillet on these sharp edges. A value of 1 is equivalent to a
+   * minSharpAngle of 180 - all edges will be smooth.
+   */
+  smoothOut(minSharpAngle?: number, minSmoothness?: number): Manifold;
+
   /**
    * Increase the density of the mesh by splitting every edge into n pieces. For
    * instance, with n = 2, each triangle will be split into 4 triangles. These
@@ -618,6 +636,24 @@ export class Manifold {
    */
   calculateCurvature(gaussianIdx: number, meanIdx: number): Manifold;
 
+  /**
+   * Fills in vertex properties for normal vectors, calculated from the mesh
+   * geometry. Flat faces composed of three or more triangles will remain flat.
+   *
+   * @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 be automatically expanded to include up through normalIdx
+   * + 2.
+   *
+   * @param minSharpAngle Any edges with angles greater than this value will
+   * remain sharp, getting different normal vector properties on each side of
+   * the edge. By default, no edges are sharp and all normals are shared. With a
+   * value of zero, the model is faceted and all normals match their triangle
+   * normals, but in this case it would be better not to calculate normals at
+   * all.
+   */
+  calculateNormals(normalIdx: number, minSharpAngle: number): Manifold;
+
   // Boolean Operations
 
   /**

src/manifold/include/manifold.h

@@ -128,8 +128,8 @@ struct MeshGL {
  * potentially Mesh if only basic geometry is required.
  *
  * In addition to storing geometric data, a Manifold can also store an arbitrary
- * number of vertex properties. These could be anything, e.g. UV coordinates,
- * colors, bone weights, etc, but this library is completely agnostic. All
+ * number of vertex properties. These could be anything, e.g. normals, UV
+ * coordinates, colors, etc, but this library is completely agnostic. All
  * properties are merely float values indexed by channel number. It is up to the
  * user to associate channel numbers with meaning.
  *
@@ -244,6 +244,8 @@ class Manifold {
   Manifold SetProperties(
       int, std::function<void(float*, glm::vec3, const float*)>) const;
   Manifold CalculateCurvature(int gaussianIdx, int meanIdx) const;
+  Manifold CalculateNormals(int normalIdx, float minSharpAngle = 60) const;
+  Manifold SmoothOut(float minSharpAngle = 60, float minSmoothness = 0) const;
   Manifold Refine(int) const;
   Manifold RefineToLength(float) const;
   // Manifold RefineToPrecision(float);

src/manifold/src/constructors.cpp

@@ -88,7 +88,7 @@ Manifold Manifold::Smooth(const MeshGL& meshGL,
   std::vector<float> propertyTolerance(meshGL.numProp - 3, -1);
   std::shared_ptr<Impl> impl =
       std::make_shared<Impl>(meshGL, propertyTolerance);
-  impl->CreateTangents(sharpenedEdges);
+  impl->CreateTangents(impl->UpdateSharpenedEdges(sharpenedEdges));
   return Manifold(impl);
 }
 
@@ -128,7 +128,7 @@ Manifold Manifold::Smooth(const Mesh& mesh,
 
   Impl::MeshRelationD relation = {(int)ReserveIDs(1)};
   std::shared_ptr<Impl> impl = std::make_shared<Impl>(mesh, relation);
-  impl->CreateTangents(sharpenedEdges);
+  impl->CreateTangents(impl->UpdateSharpenedEdges(sharpenedEdges));
   return Manifold(impl);
 }
 

src/manifold/src/impl.h

@@ -137,7 +137,14 @@ struct Manifold::Impl {
   void SplitPinchedVerts();
 
   // smoothing.cu
-  void CreateTangents(const std::vector<Smoothness>&);
+  std::vector<Smoothness> UpdateSharpenedEdges(
+      const std::vector<Smoothness>&) const;
+  Vec<bool> FlatFaces() const;
+  Vec<int> VertFlatFace(const Vec<bool>&) const;
+  std::vector<Smoothness> SharpenEdges(float minSharpAngle,
+                                       float minSmoothness) const;
+  void SetNormals(int normalIdx, float minSharpAngle);
+  void CreateTangents(std::vector<Smoothness>);
   Vec<Barycentric> Subdivide(std::function<int(glm::vec3)>);
   void Refine(std::function<int(glm::vec3)>);
 };

src/manifold/src/manifold.cpp

@@ -655,6 +655,51 @@ Manifold Manifold::CalculateCurvature(int gaussianIdx, int meanIdx) const {
   return Manifold(std::make_shared<CsgLeafNode>(pImpl));
 }
 
+/**
+ * Fills in vertex properties for normal vectors, calculated from the mesh
+ * geometry. Flat faces composed of three or more triangles will remain flat.
+ *
+ * @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 be automatically expanded to include up through normalIdx
+ * + 2.
+ *
+ * @param minSharpAngle Any edges with angles greater than this value will
+ * remain sharp, getting different normal vector properties on each side of the
+ * edge. By default, no edges are sharp and all normals are shared. With a value
+ * of zero, the model is faceted and all normals match their triangle normals,
+ * but in this case it would be better not to calculate normals at all.
+ */
+Manifold Manifold::CalculateNormals(int normalIdx, float minSharpAngle) const {
+  auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  pImpl->SetNormals(normalIdx, minSharpAngle);
+  return Manifold(std::make_shared<CsgLeafNode>(pImpl));
+}
+
+/**
+ * Smooths out the Manifold by filling in the halfedgeTangent vectors. The
+ * geometry will remain unchanged until Refine or RefineToLength is called to
+ * interpolate the surface.
+ *
+ * @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 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
+ * fillet on these sharp edges. A value of 1 is equivalent to a minSharpAngle of
+ * 180 - all edges will be smooth.
+ */
+Manifold Manifold::SmoothOut(float minSharpAngle, float minSmoothness) const {
+  auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  if (!IsEmpty()) {
+    pImpl->CreateTangents(pImpl->SharpenEdges(minSharpAngle, minSmoothness));
+  }
+  return Manifold(std::make_shared<CsgLeafNode>(pImpl));
+}
+
 /**
  * Increase the density of the mesh by splitting every edge into n pieces. For
  * instance, with n = 2, each triangle will be split into 4 triangles. These

src/manifold/src/shared.h

@@ -157,6 +157,10 @@ struct TriRef {
   /// The triangle index of the original triangle this was part of:
   /// Mesh.triVerts[tri].
   int tri;
+
+  bool SameFace(const TriRef& other) const {
+    return meshID == other.meshID && tri == other.tri;
+  }
 };
 
 /**

src/manifold/src/smoothing.cpp

@@ -501,6 +501,271 @@ class Partition {
 
 namespace manifold {
 
+// sharpenedEdges are referenced to the input Mesh, but the triangles have
+// been sorted in creating the Manifold, so the indices are converted using
+// meshRelation_.
+std::vector<Smoothness> Manifold::Impl::UpdateSharpenedEdges(
+    const std::vector<Smoothness>& sharpenedEdges) const {
+  std::unordered_map<int, int> oldHalfedge2New;
+  for (int tri = 0; tri < NumTri(); ++tri) {
+    int oldTri = meshRelation_.triRef[tri].tri;
+    for (int i : {0, 1, 2}) oldHalfedge2New[3 * oldTri + i] = 3 * tri + i;
+  }
+  std::vector<Smoothness> newSharp = sharpenedEdges;
+  for (Smoothness& edge : newSharp) {
+    edge.halfedge = oldHalfedge2New[edge.halfedge];
+  }
+  return newSharp;
+}
+
+// Find faces containing at least 3 triangles - these will not have
+// interpolated normals - all their vert normals must match their face normal.
+Vec<bool> Manifold::Impl::FlatFaces() const {
+  const int numTri = NumTri();
+  Vec<bool> triIsFlatFace(numTri, false);
+  for_each_n(
+      autoPolicy(numTri), countAt(0), numTri,
+      [this, &triIsFlatFace](const int tri) {
+        const TriRef& ref = this->meshRelation_.triRef[tri];
+        int faceNeighbors = 0;
+        glm::ivec3 faceTris = {-1, -1, -1};
+        for (const int j : {0, 1, 2}) {
+          const int neighborTri =
+              this->halfedge_[this->halfedge_[3 * tri + j].pairedHalfedge].face;
+          const TriRef& jRef = this->meshRelation_.triRef[neighborTri];
+          if (jRef.SameFace(ref)) {
+            ++faceNeighbors;
+            faceTris[j] = neighborTri;
+          }
+        }
+        if (faceNeighbors > 1) {
+          triIsFlatFace[tri] = true;
+          for (const int j : {0, 1, 2}) {
+            if (faceTris[j] >= 0) {
+              triIsFlatFace[faceTris[j]] = true;
+            }
+          }
+        }
+      });
+  return triIsFlatFace;
+}
+
+// Returns a vector of length numVert that has a tri that is part of a
+// neighboring flat face if there is only one flat face. If there are none it
+// gets -1, and if there are more than one it gets -2.
+Vec<int> Manifold::Impl::VertFlatFace(const Vec<bool>& flatFaces) const {
+  Vec<int> vertFlatFace(NumVert(), -1);
+  Vec<TriRef> vertRef(NumVert(), {-1, -1, -1});
+  for (int tri = 0; tri < NumTri(); ++tri) {
+    if (flatFaces[tri]) {
+      for (const int j : {0, 1, 2}) {
+        const int vert = halfedge_[3 * tri + j].startVert;
+        if (vertRef[vert].SameFace(meshRelation_.triRef[tri])) continue;
+        vertRef[vert] = meshRelation_.triRef[tri];
+        vertFlatFace[vert] = vertFlatFace[vert] == -1 ? tri : -2;
+      }
+    }
+  }
+  return vertFlatFace;
+}
+
+std::vector<Smoothness> Manifold::Impl::SharpenEdges(
+    float minSharpAngle, float minSmoothness) const {
+  std::vector<Smoothness> sharpenedEdges;
+  const float minRadians = glm::radians(minSharpAngle);
+  for (int e = 0; e < halfedge_.size(); ++e) {
+    if (!halfedge_[e].IsForward()) continue;
+    const int pair = halfedge_[e].pairedHalfedge;
+    const float dihedral =
+        glm::acos(glm::dot(faceNormal_[e / 3], faceNormal_[pair / 3]));
+    if (dihedral > minRadians) {
+      sharpenedEdges.push_back({e, minSmoothness});
+      sharpenedEdges.push_back({pair, minSmoothness});
+    }
+  }
+  return sharpenedEdges;
+}
+
+/**
+ * Instead of calculating the internal shared normals like CalculateNormals
+ * does, this method fills in vertex properties, unshared across edges that
+ * are bent more than minSharpAngle.
+ */
+void Manifold::Impl::SetNormals(int normalIdx, float minSharpAngle) {
+  if (IsEmpty()) return;
+  if (normalIdx < 0) return;
+
+  const int oldNumProp = NumProp();
+  const int numTri = NumTri();
+
+  Vec<bool> triIsFlatFace = FlatFaces();
+  Vec<int> vertFlatFace = VertFlatFace(triIsFlatFace);
+  Vec<int> vertNumSharp(NumVert(), 0);
+  for (int e = 0; e < halfedge_.size(); ++e) {
+    if (!halfedge_[e].IsForward()) continue;
+    const int pair = halfedge_[e].pairedHalfedge;
+    const int tri1 = e / 3;
+    const int tri2 = pair / 3;
+    const float dihedral =
+        glm::degrees(glm::acos(glm::dot(faceNormal_[tri1], faceNormal_[tri2])));
+    if (dihedral > minSharpAngle) {
+      ++vertNumSharp[halfedge_[e].startVert];
+      ++vertNumSharp[halfedge_[e].endVert];
+    } else {
+      const bool faceSplit =
+          triIsFlatFace[tri1] != triIsFlatFace[tri2] ||
+          (triIsFlatFace[tri1] && triIsFlatFace[tri2] &&
+           !meshRelation_.triRef[tri1].SameFace(meshRelation_.triRef[tri2]));
+      if (vertFlatFace[halfedge_[e].startVert] == -2 && faceSplit) {
+        ++vertNumSharp[halfedge_[e].startVert];
+      }
+      if (vertFlatFace[halfedge_[e].endVert] == -2 && faceSplit) {
+        ++vertNumSharp[halfedge_[e].endVert];
+      }
+    }
+  }
+
+  const int numProp = glm::max(oldNumProp, normalIdx + 3);
+  Vec<float> oldProperties(numProp * NumPropVert(), 0);
+  meshRelation_.properties.swap(oldProperties);
+  meshRelation_.numProp = numProp;
+  if (meshRelation_.triProperties.size() == 0) {
+    meshRelation_.triProperties.resize(numTri);
+    for_each_n(autoPolicy(numTri), countAt(0), numTri, [this](int tri) {
+      for (const int j : {0, 1, 2})
+        this->meshRelation_.triProperties[tri][j] =
+            this->halfedge_[3 * tri + j].startVert;
+    });
+  }
+  Vec<glm::ivec3> oldTriProp(numTri, {-1, -1, -1});
+  meshRelation_.triProperties.swap(oldTriProp);
+
+  for (int tri = 0; tri < numTri; ++tri) {
+    for (const int i : {0, 1, 2}) {
+      if (meshRelation_.triProperties[tri][i] >= 0) continue;
+      int current = 3 * tri + i;
+      const int startEdge = current;
+      const int vert = halfedge_[current].startVert;
+
+      if (vertNumSharp[vert] < 2) {
+        const glm::vec3 normal = vertFlatFace[vert] >= 0
+                                     ? faceNormal_[vertFlatFace[vert]]
+                                     : vertNormal_[vert];
+        int lastProp = -1;
+        do {
+          current = NextHalfedge(halfedge_[current].pairedHalfedge);
+          const int thisTri = current / 3;
+          const int j = current - 3 * thisTri;
+          const int prop = oldTriProp[thisTri][j];
+          meshRelation_.triProperties[thisTri][j] = prop;
+          if (prop == lastProp) continue;
+          lastProp = prop;
+          auto start = oldProperties.begin() + prop * oldNumProp;
+          std::copy(start, start + oldNumProp,
+                    meshRelation_.properties.begin() + prop * numProp);
+          for (const int i : {0, 1, 2})
+            meshRelation_.properties[prop * numProp + normalIdx + i] =
+                normal[i];
+        } while (current != startEdge);
+      } else {
+        const glm::vec3 centerPos = vertPos_[vert];
+        // Length degree
+        std::vector<int> group;
+        // Length number of normals
+        std::vector<glm::vec3> normals;
+        int prevFace = halfedge_[current].face;
+
+        do {
+          int next = NextHalfedge(halfedge_[current].pairedHalfedge);
+          const int face = halfedge_[next].face;
+
+          const float dihedral = glm::degrees(
+              glm::acos(glm::dot(faceNormal_[face], faceNormal_[prevFace])));
+          if (dihedral > minSharpAngle ||
+              triIsFlatFace[face] != triIsFlatFace[prevFace] ||
+              (triIsFlatFace[face] && triIsFlatFace[prevFace] &&
+               !meshRelation_.triRef[face].SameFace(
+                   meshRelation_.triRef[prevFace]))) {
+            break;
+          }
+          current = next;
+          prevFace = face;
+        } while (current != startEdge);
+
+        const int endEdge = current;
+        glm::vec3 prevEdgeVec =
+            glm::normalize(vertPos_[halfedge_[current].endVert] - centerPos);
+
+        do {
+          current = NextHalfedge(halfedge_[current].pairedHalfedge);
+          const int face = halfedge_[current].face;
+
+          const float dihedral = glm::degrees(
+              glm::acos(glm::dot(faceNormal_[face], faceNormal_[prevFace])));
+          if (dihedral > minSharpAngle ||
+              triIsFlatFace[face] != triIsFlatFace[prevFace] ||
+              (triIsFlatFace[face] && triIsFlatFace[prevFace] &&
+               !meshRelation_.triRef[face].SameFace(
+                   meshRelation_.triRef[prevFace]))) {
+            normals.push_back(glm::vec3(0));
+          }
+          group.push_back(normals.size() - 1);
+
+          const glm::vec3 edgeVec =
+              glm::normalize(vertPos_[halfedge_[current].endVert] - centerPos);
+          float dot = glm::dot(prevEdgeVec, edgeVec);
+          const float phi =
+              dot >= 1 ? 0 : (dot <= -1 ? glm::pi<float>() : glm::acos(dot));
+          normals.back() += faceNormal_[face] * phi;
+
+          prevFace = face;
+          prevEdgeVec = edgeVec;
+        } while (current != endEdge);
+
+        for (auto& normal : normals) {
+          normal = glm::normalize(normal);
+        }
+
+        int lastGroup = 0;
+        int lastProp = -1;
+        int newProp = -1;
+        int idx = 0;
+        do {
+          current = NextHalfedge(halfedge_[current].pairedHalfedge);
+          const int thisTri = current / 3;
+          const int j = current - 3 * thisTri;
+          const int prop = oldTriProp[thisTri][j];
+          auto start = oldProperties.begin() + prop * oldNumProp;
+
+          if (group[idx] != lastGroup && group[idx] != 0 && prop == lastProp) {
+            lastGroup = group[idx];
+            newProp = NumPropVert();
+            meshRelation_.properties.resize(meshRelation_.properties.size() +
+                                            numProp);
+            std::copy(start, start + oldNumProp,
+                      meshRelation_.properties.begin() + newProp * numProp);
+            for (const int i : {0, 1, 2}) {
+              meshRelation_.properties[newProp * numProp + normalIdx + i] =
+                  normals[group[idx]][i];
+            }
+          } else if (prop != lastProp) {
+            lastProp = prop;
+            newProp = prop;
+            std::copy(start, start + oldNumProp,
+                      meshRelation_.properties.begin() + prop * numProp);
+            for (const int i : {0, 1, 2})
+              meshRelation_.properties[prop * numProp + normalIdx + i] =
+                  normals[group[idx]][i];
+          }
+
+          meshRelation_.triProperties[thisTri][j] = newProp;
+          ++idx;
+        } while (current != endEdge);
+      }
+    }
+  }
+}
+
 /**
  * Calculates halfedgeTangent_, allowing the manifold to be refined and
  * smoothed. The tangents form weighted cubic Beziers along each edge. This
@@ -510,101 +775,109 @@ namespace manifold {
  * curvature there, while the tangents of the sharp edges themselves are aligned
  * for continuity.
  */
-void Manifold::Impl::CreateTangents(
-    const std::vector<Smoothness>& sharpenedEdges) {
+void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
   ZoneScoped;
   const int numHalfedge = halfedge_.size();
   halfedgeTangent_.resize(numHalfedge);
 
+  Vec<bool> triIsFlatFace = FlatFaces();
+  Vec<int> vertFlatFace = VertFlatFace(triIsFlatFace);
+  Vec<glm::vec3> vertNormal = vertNormal_;
+  for (int v = 0; v < NumVert(); ++v) {
+    if (vertFlatFace[v] >= 0) {
+      vertNormal[v] = faceNormal_[vertFlatFace[v]];
+    }
+  }
+
   for_each_n(autoPolicy(numHalfedge),
              zip(halfedgeTangent_.begin(), halfedge_.cbegin()), numHalfedge,
-             SmoothBezier({vertPos_, faceNormal_, vertNormal_, halfedge_}));
+             SmoothBezier({vertPos_, faceNormal_, vertNormal, halfedge_}));
 
-  if (!sharpenedEdges.empty()) {
-    const Vec<TriRef>& triRef = meshRelation_.triRef;
-
-    // sharpenedEdges are referenced to the input Mesh, but the triangles have
-    // been sorted in creating the Manifold, so the indices are converted using
-    // meshRelation_.
-    std::vector<int> oldHalfedge2New(halfedge_.size());
-    for (int tri = 0; tri < NumTri(); ++tri) {
-      int oldTri = triRef[tri].tri;
-      for (int i : {0, 1, 2}) oldHalfedge2New[3 * oldTri + i] = 3 * tri + i;
-    }
-
-    using Pair = std::pair<Smoothness, Smoothness>;
-    // Fill in missing pairs with default smoothness = 1.
-    std::map<int, Pair> edges;
-    for (Smoothness edge : sharpenedEdges) {
-      if (edge.smoothness == 1) continue;
-      edge.halfedge = oldHalfedge2New[edge.halfedge];
-      int pair = halfedge_[edge.halfedge].pairedHalfedge;
-      if (edges.find(pair) == edges.end()) {
-        edges[edge.halfedge] = {edge, {pair, 1}};
-      } else {
-        edges[pair].second = edge;
+  // Add sharpened edges around faces, just on the face side.
+  for (int tri = 0; tri < NumTri(); ++tri) {
+    if (!triIsFlatFace[tri]) continue;
+    for (const int j : {0, 1, 2}) {
+      const int tri2 = halfedge_[3 * tri + j].pairedHalfedge / 3;
+      if (!triIsFlatFace[tri2] ||
+          !meshRelation_.triRef[tri].SameFace(meshRelation_.triRef[tri2])) {
+        sharpenedEdges.push_back({3 * tri + j, 0});
       }
     }
+  }
 
-    std::map<int, std::vector<Pair>> vertTangents;
-    for (const auto& value : edges) {
-      const Pair edge = value.second;
-      vertTangents[halfedge_[edge.first.halfedge].startVert].push_back(edge);
-      vertTangents[halfedge_[edge.second.halfedge].startVert].push_back(
-          {edge.second, edge.first});
+  if (sharpenedEdges.empty()) return;
+
+  using Pair = std::pair<Smoothness, Smoothness>;
+  // Fill in missing pairs with default smoothness = 1.
+  std::map<int, Pair> edges;
+  for (Smoothness edge : sharpenedEdges) {
+    if (edge.smoothness >= 1) continue;
+    const bool forward = halfedge_[edge.halfedge].IsForward();
+    const int pair = halfedge_[edge.halfedge].pairedHalfedge;
+    const int idx = forward ? edge.halfedge : pair;
+    if (edges.find(idx) == edges.end()) {
+      edges[idx] = {edge, {pair, 1}};
+      if (!forward) std::swap(edges[idx].first, edges[idx].second);
+    } else {
+      Smoothness& e = forward ? edges[idx].first : edges[idx].second;
+      e.smoothness = glm::min(edge.smoothness, e.smoothness);
     }
+  }
 
-    Vec<glm::vec4>& tangent = halfedgeTangent_;
-    for (const auto& value : vertTangents) {
-      const std::vector<Pair>& vert = value.second;
-      // Sharp edges that end are smooth at their terminal vert.
-      if (vert.size() == 1) continue;
-      if (vert.size() == 2) {  // Make continuous edge
-        const int first = vert[0].first.halfedge;
-        const int second = vert[1].first.halfedge;
-        const glm::vec3 newTangent = glm::normalize(glm::vec3(tangent[first]) -
-                                                    glm::vec3(tangent[second]));
-        tangent[first] =
-            glm::vec4(glm::length(glm::vec3(tangent[first])) * newTangent,
-                      tangent[first].w);
-        tangent[second] =
-            glm::vec4(-glm::length(glm::vec3(tangent[second])) * newTangent,
-                      tangent[second].w);
+  std::map<int, std::vector<Pair>> vertTangents;
+  for (const auto& value : edges) {
+    const Pair edge = value.second;
+    vertTangents[halfedge_[edge.first.halfedge].startVert].push_back(edge);
+    vertTangents[halfedge_[edge.second.halfedge].startVert].push_back(
+        {edge.second, edge.first});
+  }
 
-        auto SmoothHalf = [&](int first, int last, float smoothness) {
-          int current = NextHalfedge(halfedge_[first].pairedHalfedge);
-          while (current != last) {
-            const float cosBeta = glm::dot(
-                newTangent, glm::normalize(glm::vec3(tangent[current])));
-            const float factor =
-                (1 - smoothness) * cosBeta * cosBeta + smoothness;
-            tangent[current] = glm::vec4(factor * glm::vec3(tangent[current]),
-                                         tangent[current].w);
-            current = NextHalfedge(halfedge_[current].pairedHalfedge);
-          }
-        };
+  Vec<glm::vec4>& tangent = halfedgeTangent_;
+  for (const auto& value : vertTangents) {
+    const std::vector<Pair>& vert = value.second;
+    // Sharp edges that end are smooth at their terminal vert.
+    if (vert.size() == 1) continue;
+    if (vert.size() == 2) {  // Make continuous edge
+      const int first = vert[0].first.halfedge;
+      const int second = vert[1].first.halfedge;
+      const glm::vec3 newTangent = glm::normalize(glm::vec3(tangent[first]) -
+                                                  glm::vec3(tangent[second]));
+      tangent[first] =
+          glm::vec4(glm::length(glm::vec3(tangent[first])) * newTangent,
+                    tangent[first].w);
+      tangent[second] =
+          glm::vec4(-glm::length(glm::vec3(tangent[second])) * newTangent,
+                    tangent[second].w);
 
-        SmoothHalf(first, second,
-                   (vert[0].second.smoothness + vert[1].first.smoothness) / 2);
-        SmoothHalf(second, first,
-                   (vert[1].second.smoothness + vert[0].first.smoothness) / 2);
-
-      } else {  // Sharpen vertex uniformly
-        float smoothness = 0;
-        for (const Pair& pair : vert) {
-          smoothness += pair.first.smoothness;
-          smoothness += pair.second.smoothness;
-        }
-        smoothness /= 2 * vert.size();
-
-        const int start = vert[0].first.halfedge;
-        int current = start;
-        do {
+      auto SmoothHalf = [&](int first, int last, float smoothness) {
+        int current = NextHalfedge(halfedge_[first].pairedHalfedge);
+        while (current != last) {
           tangent[current] = glm::vec4(smoothness * glm::vec3(tangent[current]),
                                        tangent[current].w);
           current = NextHalfedge(halfedge_[current].pairedHalfedge);
-        } while (current != start);
+        }
+      };
+
+      SmoothHalf(first, second,
+                 (vert[0].second.smoothness + vert[1].first.smoothness) / 2);
+      SmoothHalf(second, first,
+                 (vert[1].second.smoothness + vert[0].first.smoothness) / 2);
+
+    } else {  // Sharpen vertex uniformly
+      float smoothness = 0;
+      for (const Pair& pair : vert) {
+        smoothness += pair.first.smoothness;
+        smoothness += pair.second.smoothness;
       }
+      smoothness /= 2 * vert.size();
+
+      const int start = vert[0].first.halfedge;
+      int current = start;
+      do {
+        tangent[current] = glm::vec4(smoothness * glm::vec3(tangent[current]),
+                                     tangent[current].w);
+        current = NextHalfedge(halfedge_[current].pairedHalfedge);
+      } while (current != start);
     }
   }
 }
@@ -612,8 +885,8 @@ void Manifold::Impl::CreateTangents(
 /**
  * Split each edge into n pieces as defined by calling the edgeDivisions
  * function, and sub-triangulate each triangle accordingly. This function
- * doesn't run Finish(), as that is expensive and it'll need to be run after the
- * new vertices have moved, which is a likely scenario after refinement
+ * doesn't run Finish(), as that is expensive and it'll need to be run after
+ * the new vertices have moved, which is a likely scenario after refinement
  * (smoothing).
  */
 Vec<Barycentric> Manifold::Impl::Subdivide(
@@ -846,6 +1119,7 @@ Vec<Barycentric> Manifold::Impl::Subdivide(
 }
 
 void Manifold::Impl::Refine(std::function<int(glm::vec3)> edgeDivisions) {
+  if (IsEmpty()) return;
   Manifold::Impl old = *this;
   Vec<Barycentric> vertBary = Subdivide(edgeDivisions);
   if (vertBary.size() == 0) return;

test/manifold_test.cpp

@@ -343,6 +343,22 @@ TEST(Manifold, Smooth) {
 #endif
 }
 
+TEST(Manifold, SmoothFlat) {
+  Manifold cone = Manifold::Cylinder(5, 10, 5).SmoothOut().CalculateNormals(0);
+  Manifold smooth = cone.RefineToLength(0.1);
+  auto prop = smooth.GetProperties();
+  EXPECT_NEAR(prop.volume, 1105, 1);
+  EXPECT_NEAR(prop.surfaceArea, 759, 1);
+
+#ifdef MANIFOLD_EXPORT
+  ExportOptions options2;
+  options2.faceted = false;
+  options2.mat.normalChannels = {3, 4, 5};
+  if (options.exportModels)
+    ExportMesh("smoothCone.glb", smooth.GetMeshGL(), options2);
+#endif
+}
+
 TEST(Manifold, Smooth2Length) {
   Manifold cone = Manifold::Extrude(
       CrossSection::Circle(10, 10).Translate({10, 0}), 2, 0, 0, {0, 0});

test/samples_test.cpp

@@ -120,8 +120,8 @@ TEST(Samples, Scallop) {
       3, colorCurvature);
   CheckNormals(scallop);
   auto prop = scallop.GetProperties();
-  EXPECT_NEAR(prop.volume, 41.3, 0.1);
-  EXPECT_NEAR(prop.surfaceArea, 78.1, 0.1);
+  EXPECT_NEAR(prop.volume, 41.1, 0.1);
+  EXPECT_NEAR(prop.surfaceArea, 77.9, 0.1);
   CheckGL(scallop);
 
 #ifdef MANIFOLD_EXPORT