started implementation

src/manifold/include/manifold.h

@@ -208,6 +208,8 @@ class Manifold {
   Manifold SetProperties(
       int, std::function<void(float*, glm::vec3, const float*)>) const;
   Manifold CalculateCurvature(int gaussianIdx, int meanIdx) const;
+  Manifold CalculateNormals(glm::ivec3 normalIdx,
+                            float minSharpAngle = 180) const;
   Manifold Refine(int) const;
   // Manifold RefineToLength(float);
   // Manifold RefineToPrecision(float);

src/manifold/src/impl.cpp

@@ -817,6 +817,134 @@ void Manifold::Impl::CalculateNormals() {
   for_each(policy, vertNormal_.begin(), vertNormal_.end(), Normalize());
 }
 
+/**
+ * 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(glm::ivec3 normalIdx, float minSharpAngle) {
+  if (IsEmpty()) return;
+  if (glm::any(glm::lessThan(normalIdx, glm::ivec3(0)))) return;
+
+  VecDH<TriRef> triRefOriginal = meshRelation_.triRef;
+  const int oldNumProp = NumProp();
+  // Find flat faces regardless of existing properties
+  CreateFaces(
+      std::vector<float>(oldNumProp, std::numeric_limits<float>::infinity()));
+  const int numTri = NumTri();
+  std::vector<int> faceNumTri(numTri, 0);
+  for (const TriRef triRef : meshRelation_.triRef) {
+    ++faceNumTri[triRef.tri];
+  }
+  std::vector<bool> triIsFlatFace(numTri);
+  for (int tri = 0; tri < numTri; ++tri) {
+    triIsFlatFace[tri] = faceNumTri[meshRelation_.triRef[tri].tri] > 2;
+  }
+
+  const int numProp = glm::max(
+      oldNumProp,
+      glm::max(normalIdx[0], glm::max(normalIdx[1], normalIdx[2])) + 1);
+  const VecDH<float> oldProperties = meshRelation_.properties;
+  VecDH<float> properties(numProp * NumPropVert(), 0);
+  meshRelation_.numProp = numProp;
+  if (meshRelation_.triProperties.size() == 0) {
+    meshRelation_.triProperties.resize(numTri);
+  }
+  VecDH<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 endEdge = current;
+      const Halfedge halfedge = halfedge_[current];
+      const int vert = halfedge.startVert;
+      const glm::vec3 centerPos = vertPos_[vert];
+      int prop = oldNumProp == 0 ? vert : oldTriProp[tri][i];
+
+      // Length degree
+      std::vector<int> props;
+      std::vector<int> edges;
+      std::vector<int> group;
+      // Length number of normals
+      std::vector<glm::vec3> normals = {glm::vec3(0)};
+      std::vector<bool> isFlat = {false};
+
+      int j = 0;
+      const int prevFace =
+          halfedge_[halfedge_[NextHalfedge(NextHalfedge(current))]
+                        .pairedHalfedge]
+              .face;
+      glm::vec3 lastFaceNormal = faceNormal_[prevFace];
+      int lastFace = -1;
+      do {
+        current = halfedge_[current].pairedHalfedge;
+        const glm::vec3 edgeVec =
+            vertPos_[halfedge_[current].startVert] - centerPos;
+
+        current = NextHalfedge(current);
+        edges.push_back(current);
+        const int face = halfedge_[current].face;
+        props.push_back(oldNumProp == 0 ? vert
+                                        : oldTriProp[face][current - 3 * face]);
+
+        const float dihedral = glm::degrees(
+            glm::acos(glm::dot(faceNormal_[face], lastFaceNormal)));
+        if (dihedral > minSharpAngle ||
+            (isFlat[j] && triIsFlatFace[face] &&
+             lastFace != meshRelation_.triRef[tri].tri)) {
+          ++j;
+          normals.push_back(glm::vec3(0));
+          isFlat.push_back(false);
+        }
+        group.push_back(j);
+
+        if (!isFlat[j]) {
+          if (triIsFlatFace[face]) {
+            normals[j] = faceNormal_[face];
+            isFlat[j] = true;
+            lastFace = meshRelation_.triRef[tri].tri;
+          } else {
+            float dot = -glm::dot(
+                edgeVec, vertPos_[halfedge_[current].endVert] - centerPos);
+            const float phi =
+                dot >= 1 ? 0 : (dot <= -1 ? glm::pi<float>() : glm::acos(dot));
+            normals[j] += faceNormal_[face] * phi;
+          }
+        }
+
+        lastFaceNormal = faceNormal_[face];
+      } while (current != endEdge);
+
+      if (normals.size() > 1) {
+        // The first and last are part of the same face
+        normals.front() += normals.back();
+        normals.back() = normals.front();
+      }
+
+      for (auto& normal : normals) {
+        normal = glm::normalize(normal);
+      }
+
+      const int degree = edges.size();
+      for (j = 0; j < degree; ++j) {
+        const int tri = halfedge_[edges[j]].face;
+        // Need to duplicate properties for unique (props[j],group[j]) pairs
+        meshRelation_.triProperties[tri][edges[j] - 3 * tri] = props[j];
+        for (int i = 0; i < oldNumProp; ++i)
+          meshRelation_.properties[props[j] * numProp + i] =
+              oldProperties[props[j] * oldNumProp + i];
+        for (const int i : {0, 1, 2})
+          meshRelation_.properties[props[j] * numProp + normalIdx[i]] =
+              normals[group[j]][i];
+      }
+    }
+  }
+
+  meshRelation_.triRef = triRefOriginal;
+}
+
 /**
  * Remaps all the contained meshIDs to new unique values to represent new
  * instances of these meshes.

src/manifold/src/impl.h

@@ -71,6 +71,7 @@ struct Manifold::Impl {
   void InitializeOriginal();
   void CreateHalfedges(const VecDH<glm::ivec3>& triVerts);
   void CalculateNormals();
+  void SetNormals(glm::ivec3 normalIdx, float minSharpAngle);
   void IncrementMeshIDs();
 
   void Update();

src/manifold/src/manifold.cpp

@@ -629,6 +629,27 @@ 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 channels in which to store the X, Y, and Z
+ * values of the normals. The property set will be automatically expanded to
+ * include the channel indices specified.
+ *
+ * @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(glm::ivec3 normalIdx,
+                                    float minSharpAngle) const {
+  auto pImpl = std::make_shared<Impl>(*GetCsgLeafNode().GetImpl());
+  pImpl->SetNormals(normalIdx, minSharpAngle);
+  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