consistent transformations of normals (#1602)

* consistent transformations of normals

* added backside API to MeshGL

* fix GCC compile

* addressing feedback

* stray struct

.vscode/launch.json

@@ -12,7 +12,7 @@
       "program": "${workspaceFolder}/build/test/manifold_test",
       "args": [
         "--gtest_catch_exceptions=0",
-        "--gtest_filter=SDF.SineSurface"
+        "--gtest_filter=Manifold.MeshRelation"
       ],
       "stopAtEntry": false,
       "cwd": "${workspaceFolder}/build/test",

CMakeLists.txt

@@ -159,6 +159,7 @@ else()
     -Wno-unknown-warning-option
     -Wno-unused
     -Wno-shorten-64-to-32
+    -Wno-sign-compare
     -Wno-alloc-size-larger-than
   )
   if(CMAKE_SYSTEM_NAME STREQUAL "Windows")

include/manifold/manifold.h

@@ -115,6 +115,8 @@ struct MeshGLP {
   I NumVert() const { return vertProperties.size() / numProp; };
   /// Number of triangles
   I NumTri() const { return triVerts.size() / 3; };
+  /// Number of triangle runs
+  I NumRun() const { return runOriginalID.size(); };
   /// Number of properties per vertex, always >= 3.
   I numProp = 3;
   /// Flat, GL-style interleaved list of all vertex properties: propVal =
@@ -150,6 +152,10 @@ struct MeshGLP {
   /// This matrix is stored in column-major order and the length of the overall
   /// vector is 12 * runOriginalID.size().
   std::vector<Precision> runTransform;
+  /// Optional: For each run, defines a set of flags giving extra information
+  /// about the run. See the corresponding getter functions for details on the
+  /// specific flags. These are primarily used on output.
+  std::vector<uint8_t> runFlags;
   /// Optional: Length NumTri, contains the source face ID this triangle comes
   /// from. Simplification will maintain all edges between triangles with
   /// different faceIDs. Input faceIDs will be maintained to the outputs, but if
@@ -170,18 +176,30 @@ struct MeshGLP {
   MeshGLP() = default;
 
   /**
-   * Updates the mergeFromVert and mergeToVert vectors in order to create a
-   * manifold solid. If the MeshGL is already manifold, no change will occur and
-   * the function will return false. Otherwise, this will merge verts along open
-   * edges within tolerance (the maximum of the MeshGL tolerance and the
-   * baseline bounding-box tolerance), keeping any from the existing merge
-   * vectors, and return true.
+   * Updates the normals of the mesh based on the runTransform matrices and
+   * backside flags, which are then cleared to avoid double-applying them when
+   * round-tripping.
    *
-   * There is no guarantee the result will be manifold - this is a best-effort
-   * helper function designed primarily to aid in the case where a manifold
-   * multi-material MeshGL was produced, but its merge vectors were lost due to
-   * a round-trip through a file format. Constructing a Manifold from the result
-   * will report an error status if it is not manifold.
+   * @param normalIdx Specifies the first of the three consecutive property
+   * channels forming the (x, y, z) normals to update. NumProp must be at least
+   * normalIdx + 3 and normalIdx must be >= 3.
+   */
+  void UpdateNormals(int normalIdx);
+
+  /**
+   * Updates the mergeFromVert and mergeToVert vectors in order to create a
+   * manifold solid. If the MeshGL is already manifold, no change will occur
+   * and the function will return false. Otherwise, this will merge verts
+   * along open edges within tolerance (the maximum of the MeshGL tolerance
+   * and the baseline bounding-box tolerance), keeping any from the existing
+   * merge vectors, and return true.
+   *
+   * There is no guarantee the result will be manifold - this is a
+   * best-effort helper function designed primarily to aid in the case where
+   * a manifold multi-material MeshGL was produced, but its merge vectors
+   * were lost due to a round-trip through a file format. Constructing a
+   * Manifold from the result will report an error status if it is not
+   * manifold.
    */
   bool Merge();
 
@@ -219,6 +237,30 @@ struct MeshGLP {
         halfedgeTangent[offset], halfedgeTangent[offset + 1],
         halfedgeTangent[offset + 2], halfedgeTangent[offset + 3]);
   }
+
+  /**
+   * Returns the transformation matrix for the specified run.
+   *
+   * @param run The index of the triangle run (0 <= run < runOriginalID.size()).
+   */
+  mat3x4 GetRunTransform(size_t run) const {
+    size_t offset = 12 * run;
+    if (offset + 12 > runTransform.size()) {
+      return la::identity;
+    }
+    return mat3x4(la::mat<Precision, 3, 4>(&runTransform[offset]));
+  }
+
+  /**
+   * Returns true if this triangle run is on the backside compared to the
+   * original mesh, e.g. from a subtraction. In this case vertex normals will
+   * need to be flipped. UpdateNormals() will take care of this.
+   *
+   * @param run The index of the triangle run (0 <= run < runFlags.size()).
+   */
+  bool Backside(size_t run) const {
+    return run < runFlags.size() && runFlags[run] == 1;
+  }
 };
 
 /**

src/impl.h

@@ -29,6 +29,14 @@ struct Manifold::Impl {
     int originalID = -1;
     mat3x4 transform = la::identity;
     bool backSide = false;
+
+    mat3 GetNormalTransform() const {
+      return NormalTransform(transform) * (backSide ? -1.0 : 1.0);
+    }
+
+    mat3 GetInverseNormalTransform() const {
+      return InverseNormalTransform(transform) * (backSide ? -1.0 : 1.0);
+    }
   };
   struct MeshRelationD {
     /// The originalID of this Manifold if it is an original; -1 otherwise.
@@ -184,6 +192,7 @@ struct Manifold::Impl {
     for (size_t i = 0; i < runOriginalID.size(); ++i) {
       const int meshID = startID + i;
       const int originalID = runOriginalID[i];
+      const bool backside = meshGL.Backside(i);
       for (size_t tri = runIndex[i] / 3; tri < runIndex[i + 1] / 3; ++tri) {
         TriRef& ref = triRef[tri];
         ref.meshID = meshID;
@@ -200,7 +209,8 @@ struct Manifold::Impl {
                                                  {{m[0], m[1], m[2]},
                                                   {m[3], m[4], m[5]},
                                                   {m[6], m[7], m[8]},
-                                                  {m[9], m[10], m[11]}}};
+                                                  {m[9], m[10], m[11]}},
+                                                 backside};
       }
     }
 
@@ -453,11 +463,12 @@ inline MeshGLP<Precision, I> GetMeshGLImpl(const manifold::Manifold::Impl& impl,
                     const manifold::Manifold::Impl::Relation& rel) {
     out.runIndex.push_back(3 * tri);
     out.runOriginalID.push_back(rel.originalID);
+    out.runFlags.push_back(rel.backSide ? 1 : 0);
     if (updateNormals) {
-      runNormalTransform.push_back(NormalTransform(rel.transform) *
-                                   (rel.backSide ? -1.0 : 1.0));
+      runNormalTransform.push_back(rel.GetNormalTransform());
     }
-    if (!isOriginal) {
+    // clear transforms if applying them to normals
+    if (!isOriginal && !updateNormals) {
       for (const int col : {0, 1, 2, 3}) {
         for (const int row : {0, 1, 2}) {
           out.runTransform.push_back(rel.transform[col][row]);

src/manifold.cpp

@@ -182,18 +182,19 @@ Manifold::Manifold(const MeshGL64& meshGL64)
     : pNode_(std::make_shared<CsgLeafNode>(std::make_shared<Impl>(meshGL64))) {}
 
 /**
- * The most complete output of this library, returning a MeshGL that is designed
+ * Returns a MeshGL that is designed
  * to easily push into a renderer, including all interleaved vertex properties
  * that may have been input. It also includes relations to all the input meshes
  * that form a part of this result and the transforms applied to each.
  *
- * @param normalIdx If the original MeshGL inputs that formed this manifold had
- * properties corresponding to normal vectors, you can specify the first of the
- * three consecutive property channels forming the (x, y, z) normals, which will
- * cause this output MeshGL to automatically update these normals according to
- * the applied transforms and front/back side. normalIdx + 3 must be <=
- * numProp, and all original MeshGLs must use the same channels for their
- * normals.
+ * @param normalIdx If this manifold has properties corresponding to normal
+ * vectors, you can specify the first of the three consecutive property channels
+ * forming the (x, y, z) normals, which will cause this output MeshGL to
+ * automatically update these normals according to the applied transforms and
+ * front/back side. normalIdx + 3 must be <= numProp, and all original meshes
+ * must use the same channels for their normals. Default is -1, meaning no
+ * normals. If normals are selected, the runTransform matrices will be removed
+ * from the output, to avoid them being double-applied when round-tripped.
  */
 MeshGL Manifold::GetMeshGL(int normalIdx) const {
   const Impl& impl = *GetCsgLeafNode().GetImpl();
@@ -201,18 +202,20 @@ MeshGL Manifold::GetMeshGL(int normalIdx) const {
 }
 
 /**
- * The most complete output of this library, returning a MeshGL that is designed
+ * Returns a MeshGL that is designed
  * to easily push into a renderer, including all interleaved vertex properties
  * that may have been input. It also includes relations to all the input meshes
  * that form a part of this result and the transforms applied to each.
  *
- * @param normalIdx If the original MeshGL inputs that formed this manifold had
- * properties corresponding to normal vectors, you can specify the first of the
- * three consecutive property channels forming the (x, y, z) normals, which will
- * cause this output MeshGL to automatically update these normals according to
- * the applied transforms and front/back side. normalIdx + 3 must be <=
- * numProp, and all original MeshGLs must use the same channels for their
- * normals.
+ * @param normalIdx If this manifold has properties corresponding to normal
+ * vectors, you can specify the first of the three consecutive property channels
+ * forming the (x, y, z) normals, which will cause this output MeshGL to
+ * automatically update these normals according to the applied transforms and
+ * front/back side. normalIdx + 3 must be <= numProp, and all original meshes
+ * must use the same channels for their normals. Default is -1, meaning no
+ * normals. If normals are selected, the runTransform matrices will be removed
+ * from the output, to avoid them being double-applied when round-tripped.
+
  */
 MeshGL64 Manifold::GetMeshGL64(int normalIdx) const {
   const Impl& impl = *GetCsgLeafNode().GetImpl();

src/shared.h

@@ -36,10 +36,23 @@ inline int NextHalfedge(int current) {
   return current;
 }
 
+/**
+When a transform is applied to the verts of an object, and different transform
+is needed for the normals. Note that applying this transform will stretch their
+length, so they should be renormalized after.
+*/
 inline mat3 NormalTransform(const mat3x4& transform) {
   return la::inverse(la::transpose(mat3(transform)));
 }
 
+/**
+This is the corresponding normal transform for the inverse of the input
+transform.
+*/
+inline mat3 InverseNormalTransform(const mat3x4& transform) {
+  return la::inverse(la::transpose(la::inverse(mat3(transform))));
+}
+
 /**
  * By using the closest axis-aligned projection to the normal instead of a
  * projection along the normal, we avoid introducing any rounding error.

src/smoothing.cpp

@@ -256,12 +256,17 @@ namespace manifold {
  * normalIdx shows the beginning of where normals are stored in the properties.
  */
 vec3 Manifold::Impl::GetNormal(int halfedge, int normalIdx) const {
+  const int meshID = meshRelation_.triRef[halfedge / 3].meshID;
+
+  const mat3 transform =
+      meshRelation_.meshIDtransform.find(meshID)->second.GetNormalTransform();
+
   const int prop = halfedge_[halfedge].propVert;
   vec3 normal;
   for (const int i : {0, 1, 2}) {
     normal[i] = properties_[prop * numProp_ + normalIdx + i];
   }
-  return normal;
+  return transform * normal;
 }
 
 /**
@@ -479,15 +484,24 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
                halfedge_[i].propVert = -1;
              });
 
+  std::map<int, mat3> meshIDtoNormalTransform;
+
   const int numEdge = halfedge_.size();
   for (int startEdge = 0; startEdge < numEdge; ++startEdge) {
     if (halfedge_[startEdge].propVert >= 0) continue;
     const int vert = halfedge_[startEdge].startVert;
 
+    const int meshID = meshRelation_.triRef[startEdge / 3].meshID;
+    if (meshIDtoNormalTransform.find(meshID) == meshIDtoNormalTransform.end()) {
+      meshIDtoNormalTransform[meshID] =
+          meshRelation_.meshIDtransform[meshID].GetInverseNormalTransform();
+    }
+    const mat3 transform = meshIDtoNormalTransform[meshID];
+
     if (vertNumSharp[vert] < 2) {  // vertex has single normal
-      const vec3 normal = vertFlatFace[vert] >= 0
-                              ? faceNormal_[vertFlatFace[vert]]
-                              : vertNormal_[vert];
+      const vec3 normal =
+          transform * (vertFlatFace[vert] >= 0 ? faceNormal_[vertFlatFace[vert]]
+                                               : vertNormal_[vert]);
       int lastProp = -1;
       ForVert(startEdge, [&](int current) {
         const int prop = oldHalfedgeProp[current];
@@ -579,7 +593,7 @@ void Manifold::Impl::SetNormals(int normalIdx, double minSharpAngle) {
           });
 
       for (auto& normal : normals) {
-        normal = SafeNormalize(normal);
+        normal = transform * SafeNormalize(normal);
       }
 
       int lastGroup = 0;

src/sort.cpp

@@ -178,6 +178,36 @@ bool MergeMeshGLP(MeshGLP<Precision, I>& mesh) {
 
   return true;
 }
+
+template <typename Precision, typename I>
+void UpdateNormalsMeshGLP(MeshGLP<Precision, I>& mesh, int normalIdx) {
+  ZoneScoped;
+  if (normalIdx < 3 || normalIdx + 3 > mesh.numProp) {
+    return;
+  }
+  std::vector<bool> vertUpdated(mesh.NumVert(), false);
+  for (size_t run = 0; run < mesh.NumRun(); ++run) {
+    const mat3 transform = NormalTransform(mesh.GetRunTransform(run)) *
+                           (mesh.Backside(run) ? -1.0 : 1.0);
+    for (I* itr = &mesh.triVerts[mesh.runIndex[run]];
+         itr < &mesh.triVerts[mesh.runIndex[run + 1]]; ++itr) {
+      const I vert = *itr;
+      if (vertUpdated[vert]) continue;
+      vertUpdated[vert] = true;
+      vec3 normal;
+      const int start = vert * mesh.numProp + normalIdx;
+      for (int i : {0, 1, 2}) {
+        normal[i] = mesh.vertProperties[start + i];
+      }
+      normal = SafeNormalize(transform * normal);
+      for (int i : {0, 1, 2}) {
+        mesh.vertProperties[start + i] = normal[i];
+      }
+    }
+  }
+  mesh.runTransform.clear();
+  mesh.runFlags.clear();
+}
 }  // namespace
 
 namespace manifold {
@@ -529,4 +559,14 @@ template <>
 bool MeshGL64::Merge() {
   return MergeMeshGLP(*this);
 }
+
+template <>
+void MeshGL::UpdateNormals(int normalIdx) {
+  UpdateNormalsMeshGLP(*this, normalIdx);
+}
+
+template <>
+void MeshGL64::UpdateNormals(int normalIdx) {
+  UpdateNormalsMeshGLP(*this, normalIdx);
+}
 }  // namespace manifold

test/boolean_complex_test.cpp

@@ -850,7 +850,7 @@ TEST(BooleanComplex, InterpolatedNormals) {
 
   auto aMinusB = aManifold - bManifold;
 
-  RelatedGL(aMinusB, {a, b}, false, false);
+  RelatedGL(aMinusB, {a, b}, true, true);
 }
 
 MeshGL mgl_0() {

test/boolean_test.cpp

@@ -72,7 +72,7 @@ TEST(Boolean, Normals) {
 
   RelatedGL(result, {cubeGL, sphereGL}, true, true);
 
-  MeshGL output = result.GetMeshGL(0);
+  MeshGL output = result.GetMeshGL();
 
   if (options.exportModels) WriteTestOBJ("normals.obj", result);
 
@@ -81,7 +81,7 @@ TEST(Boolean, Normals) {
   output.Merge();
   Manifold roundTrip(output);
 
-  RelatedGL(roundTrip, {cubeGL, sphereGL}, true, false);
+  RelatedGL(roundTrip, {cubeGL, sphereGL}, true, true);
 }
 
 TEST(Boolean, MissingNormals) {

test/smooth_test.cpp

@@ -169,6 +169,19 @@ TEST(Smooth, Normals) {
   if (options.exportModels) WriteTestOBJ("smoothNormals.obj", byNormals);
 }
 
+TEST(Smooth, NormalTransform) {
+  Manifold cube1 = Manifold::Cube().Rotate(30).CalculateNormals(0);
+  Manifold cube2 =
+      Manifold::Cube().CalculateNormals(0).Rotate(30).Translate({3, 0, 0});
+  Manifold combo = cube1 + cube2;
+  Manifold out1 = combo.SmoothByNormals(0).Refine(10);
+  EXPECT_FLOAT_EQ(out1.Volume(), 2);
+  EXPECT_FLOAT_EQ(out1.SurfaceArea(), 12);
+  Manifold out2 = Manifold(combo.GetMeshGL(0)).SmoothByNormals(0).Refine(10);
+  EXPECT_FLOAT_EQ(out2.Volume(), 2);
+  EXPECT_FLOAT_EQ(out2.SurfaceArea(), 12);
+}
+
 TEST(Smooth, Manual) {
   // Unit Octahedron
   const auto oct = Manifold::Sphere(1, 4).GetMeshGL();

test/test_main.cpp

@@ -323,19 +323,19 @@ void Identical(const MeshGL& mesh1, const MeshGL& mesh2) {
 void RelatedGL(const Manifold& out, const std::vector<MeshGL>& originals,
                bool checkNormals, bool updateNormals) {
   ASSERT_FALSE(out.IsEmpty());
-  const int normalIdx = updateNormals ? 0 : -1;
-  MeshGL output = out.GetMeshGL(normalIdx);
+  MeshGL output = out.GetMeshGL();
+  std::vector<mat3x4> runTransforms(output.NumRun());
+  for (size_t run = 0; run < output.NumRun(); ++run) {
+    runTransforms[run] = output.GetRunTransform(run);
+  }
+  if (updateNormals) {
+    output.UpdateNormals(3);
+  }
 
-  for (size_t run = 0; run < output.runOriginalID.size(); ++run) {
-    const float* m = output.runTransform.data() + 12 * run;
-    const mat3x4 transform =
-        output.runTransform.empty()
-            ? la::identity
-            : mat3x4({m[0], m[1], m[2]}, {m[3], m[4], m[5]}, {m[6], m[7], m[8]},
-                     {m[9], m[10], m[11]});
+  for (size_t run = 0; run < output.NumRun(); ++run) {
     size_t i = 0;
     for (; i < originals.size(); ++i) {
-      ASSERT_EQ(originals[i].runOriginalID.size(), 1);
+      ASSERT_EQ(originals[i].NumRun(), 1);
       if (originals[i].runOriginalID[0] == output.runOriginalID[run]) break;
     }
     ASSERT_LT(i, originals.size());
@@ -365,7 +365,7 @@ void RelatedGL(const Manifold& out, const std::vector<MeshGL>& originals,
           outTriPos[j][k] = output.vertProperties[vert * output.numProp + k];
         }
         pos[3] = 1;
-        inTriPos[j] = transform * pos;
+        inTriPos[j] = runTransforms[run] * pos;
       }
       vec3 outNormal =
           la::cross(outTriPos[1] - outTriPos[0], outTriPos[2] - outTriPos[0]);