Fix #142485: Shading artifacts with free custom normals and scale transform

When transforming a geometry, we often apply the transposed inverse
to normals / custom normal data. However, that matrix can still contain
scale from the original matrix. That scale has to be removed so we can
avoid also scaling the normals.

I used the opportunity to remove the duplication between mesh and curves
processing of the custom normals, and to formalize an optimization to
skip the final normalization of each vector if the transform is such
that it isn't necessary. The new functions don't fit beautifully into
their public headers, but I don't know of a better place for them.

Pull Request: https://projects.blender.org/blender/blender/pulls/142896

source/blender/blenkernel/BKE_attribute.hh

@@ -14,6 +14,7 @@
 #include "BLI_function_ref.hh"
 #include "BLI_generic_pointer.hh"
 #include "BLI_generic_virtual_array.hh"
+#include "BLI_math_matrix_types.hh"
 #include "BLI_offset_indices.hh"
 #include "BLI_set.hh"
 #include "BLI_struct_equality_utils.hh"
@@ -987,4 +988,10 @@ void fill_attribute_range_default(MutableAttributeAccessor dst_attributes,
                                   const AttributeFilter &attribute_filter,
                                   IndexRange range);
 
+/**
+ * Apply a transform to the "custom_normal" attribute.
+ */
+void transform_custom_normal_attribute(const float4x4 &transform,
+                                       MutableAttributeAccessor &attributes);
+
 }  // namespace blender::bke

source/blender/blenkernel/intern/attribute_access.cc

@@ -1169,4 +1169,32 @@ void fill_attribute_range_default(MutableAttributeAccessor attributes,
   });
 }
 
+void transform_custom_normal_attribute(const float4x4 &transform,
+                                       MutableAttributeAccessor &attributes)
+{
+  const GAttributeReader normals = attributes.lookup("custom_normal");
+  if (!normals) {
+    return;
+  }
+  if (!normals.varray.type().is<float3>()) {
+    return;
+  }
+  if (normals.sharing_info->is_mutable()) {
+    SpanAttributeWriter<float3> normals = attributes.lookup_for_write_span<float3>(
+        "custom_normal");
+    math::transform_normals(float3x3(transform), normals.span);
+    normals.finish();
+  }
+  else {
+    /* It's a bit faster to combine transforming and copying the attribute if it's shared. */
+    float3 *new_data = MEM_malloc_arrayN<float3>(size_t(normals.varray.size()), __func__);
+    math::transform_normals(VArraySpan(normals.varray.typed<float3>()),
+                            float3x3(transform),
+                            {new_data, normals.varray.size()});
+    const AttrDomain domain = normals.domain;
+    attributes.remove("custom_normal");
+    attributes.add<float3>("custom_normal", domain, AttributeInitMoveArray(new_data));
+  }
+}
+
 }  // namespace blender::bke

source/blender/blenkernel/intern/curves_geometry.cc

@@ -1246,16 +1246,6 @@ static void transform_positions(MutableSpan<float3> positions, const float4x4 &m
   });
 }
 
-static void transform_normals(MutableSpan<float3> normals, const float4x4 &matrix)
-{
-  const float3x3 normal_transform = math::transpose(math::invert(float3x3(matrix)));
-  threading::parallel_for(normals.index_range(), 1024, [&](const IndexRange range) {
-    for (float3 &normal : normals.slice(range)) {
-      normal = normal_transform * normal;
-    }
-  });
-}
-
 void CurvesGeometry::calculate_bezier_auto_handles()
 {
   if (!this->has_curve_with_type(CURVE_TYPE_BEZIER)) {
@@ -1325,10 +1315,7 @@ void CurvesGeometry::transform(const float4x4 &matrix)
     transform_positions(this->handle_positions_right_for_write(), matrix);
   }
   MutableAttributeAccessor attributes = this->attributes_for_write();
-  if (SpanAttributeWriter normals = attributes.lookup_for_write_span<float3>("custom_normal")) {
-    transform_normals(normals.span, matrix);
-    normals.finish();
-  }
+  transform_custom_normal_attribute(matrix, attributes);
   this->tag_positions_changed();
 }
 

source/blender/blenkernel/intern/mesh.cc

@@ -1849,16 +1849,6 @@ static void translate_positions(MutableSpan<float3> positions, const float3 &tra
   });
 }
 
-static void transform_normals(MutableSpan<float3> normals, const float4x4 &matrix)
-{
-  const float3x3 normal_transform = math::transpose(math::invert(float3x3(matrix)));
-  threading::parallel_for(normals.index_range(), 1024, [&](const IndexRange range) {
-    for (float3 &normal : normals.slice(range)) {
-      normal = normal_transform * normal;
-    }
-  });
-}
-
 void mesh_translate(Mesh &mesh, const float3 &translation, const bool do_shape_keys)
 {
   if (math::is_zero(translation)) {
@@ -1897,15 +1887,7 @@ void mesh_transform(Mesh &mesh, const float4x4 &transform, bool do_shape_keys)
     }
   }
   MutableAttributeAccessor attributes = mesh.attributes_for_write();
-  if (const std::optional<AttributeMetaData> meta_data = attributes.lookup_meta_data(
-          "custom_normal"))
-  {
-    if (meta_data->data_type == bke::AttrType::Float3) {
-      bke::SpanAttributeWriter normals = attributes.lookup_for_write_span<float3>("custom_normal");
-      transform_normals(normals.span, transform);
-      normals.finish();
-    }
-  }
+  transform_custom_normal_attribute(transform, attributes);
 
   mesh.tag_positions_changed();
 }

source/blender/blenlib/BLI_math_matrix.hh

@@ -1806,4 +1806,11 @@ extern template float4x4 perspective(
 
 /** \} */
 
+/**
+ * Transform normal vectors, maintaining their unit length status, but implementing some
+ * optimizations for identity matrics and uniform scaling.
+ */
+void transform_normals(const float3x3 &transform, MutableSpan<float3> normals);
+void transform_normals(Span<float3> src, const float3x3 &transform, MutableSpan<float3> dst);
+
 }  // namespace blender::math

source/blender/blenlib/intern/math_matrix.cc

@@ -10,6 +10,7 @@
 
 #include "BLI_math_rotation.hh"
 #include "BLI_simd.hh"
+#include "BLI_task.hh"
 
 #include <Eigen/Core>
 #include <Eigen/Dense>
@@ -517,4 +518,82 @@ template float4x4 perspective_infinite(
 
 /** \} */
 
+/**
+ * Check that each column is orthogonal to the others, and that each column is the same length.
+ * In other words, there is no shear, and any scaling is uniform.
+ */
+template<typename T>
+bool is_similarity_transform(const MatBase<T, 3, 3> &matrix, const T &epsilon = 1e-6)
+{
+  if (math::abs(math::dot(matrix[0], matrix[1])) > epsilon) {
+    return false;
+  }
+  if (math::abs(math::dot(matrix[0], matrix[2])) > epsilon) {
+    return false;
+  }
+  if (math::abs(math::dot(matrix[1], matrix[2])) > epsilon) {
+    return false;
+  }
+  const float length_0 = math::length_squared(matrix[0]);
+  const float length_1 = math::length_squared(matrix[1]);
+  const float length_2 = math::length_squared(matrix[2]);
+  if (math::abs(length_0 - length_1) > epsilon) {
+    return false;
+  }
+  if (math::abs(length_0 - length_2) > epsilon) {
+    return false;
+  }
+  if (math::abs(length_1 - length_2) > epsilon) {
+    return false;
+  }
+  return true;
+}
+
+void transform_normals(const float3x3 &transform, MutableSpan<float3> normals)
+{
+  if (math::is_equal(transform, float3x3::identity(), 1e-6f)) {
+    return;
+  }
+  const float3x3 normal_transform = math::transpose(math::invert(transform));
+  if (is_similarity_transform(normal_transform)) {
+    const float3x3 normalized_transform = math::normalize(normal_transform);
+    threading::parallel_for(normals.index_range(), 1024, [&](const IndexRange range) {
+      for (float3 &normal : normals.slice(range)) {
+        normal = normalized_transform * normal;
+      }
+    });
+  }
+  else {
+    threading::parallel_for(normals.index_range(), 1024, [&](const IndexRange range) {
+      for (float3 &normal : normals.slice(range)) {
+        normal = math::normalize(normal_transform * normal);
+      }
+    });
+  }
+}
+
+void transform_normals(Span<float3> src, const float3x3 &transform, MutableSpan<float3> dst)
+{
+  if (math::is_equal(transform, float3x3::identity(), 1e-6f)) {
+    dst.copy_from(src);
+    return;
+  }
+  const float3x3 normal_transform = math::transpose(math::invert(transform));
+  if (is_similarity_transform(normal_transform)) {
+    const float3x3 normalized_transform = math::normalize(normal_transform);
+    threading::parallel_for(src.index_range(), 1024, [&](const IndexRange range) {
+      for (const int i : range) {
+        dst[i] = normalized_transform * src[i];
+      }
+    });
+  }
+  else {
+    threading::parallel_for(src.index_range(), 1024, [&](const IndexRange range) {
+      for (const int i : range) {
+        dst[i] = math::normalize(normal_transform * src[i]);
+      }
+    });
+  }
+}
+
 }  // namespace blender::math

source/blender/blenlib/tests/BLI_math_matrix_test.cc

@@ -4,11 +4,13 @@
 
 #include "testing/testing.h"
 
+#include "BLI_array.hh"
 #include "BLI_math_base.h"
 #include "BLI_math_matrix.h"
 #include "BLI_math_matrix.hh"
 #include "BLI_math_rotation.h"
 #include "BLI_math_rotation.hh"
+#include "BLI_rand.hh"
 
 TEST(math_matrix, interp_m4_m4m4_regular)
 {
@@ -627,4 +629,43 @@ TEST(math_matrix, ToQuaternionSafe)
   EXPECT_M3_NEAR(from_rotation<float3x3>(rotation), expect, 1e-5);
 }
 
+static void transform_normals_test(const float3x3 &transform, const Span<float3> normals)
+{
+  Array<float3, 10> transformed_normals(normals.size());
+  transform_normals(normals, transform, transformed_normals);
+  /* Input unit vectors should still be unit length. */
+  for (const int64_t i : transformed_normals.index_range()) {
+    EXPECT_NEAR(math::length(transformed_normals[i]), 1.0f, 1e-6f);
+  }
+  /* Make sure any optimization inside #transform_normals that skips the final vector normalization
+   * doesn't change the result. */
+  const float3x3 normal_transform = math::transpose(math::invert(transform));
+  for (const int64_t i : normals.index_range()) {
+    const float3 transformed_then_normalized = math::normalize(normal_transform * normals[i]);
+    EXPECT_V3_NEAR(transformed_normals[i], transformed_then_normalized, 1e-6f);
+  }
+}
+
+TEST(math_matrix, TransformNormals)
+{
+  RandomNumberGenerator rng;
+  Array<float3, 10> normals(10);
+  for (const int64_t i : normals.index_range()) {
+    normals[i] = rng.get_unit_float3();
+  }
+  for (const int64_t i : normals.index_range()) {
+    EXPECT_NEAR(math::length(normals[i]), 1.0f, 1e-6f);
+  }
+
+  transform_normals_test(float3x3::identity(), normals);
+  transform_normals_test(from_rotation<float3x3>(EulerXYZ(4.2f, 0.2f, -3.2f)), normals);
+  transform_normals_test(float3x3::diagonal(0.5f), normals);
+  transform_normals_test(float3x3::diagonal(10.5f), normals);
+  transform_normals_test(from_scale<float3x3>(float3(0.5f, 1.5f, 5.7f)), normals);
+  transform_normals_test(from_rot_scale<float3x3>(EulerXYZ(0.0f, 7.1f, 3.14f), float3(2.15f)),
+                         normals);
+  transform_normals_test(
+      from_rot_scale<float3x3>(EulerXYZ(0.0f, 7.1f, 3.14f), float3(0.5f, 1.5f, 5.7f)), normals);
+}
+
 }  // namespace blender::tests

source/blender/geometry/intern/realize_instances.cc

@@ -468,23 +468,6 @@ static void transform_positions(const float4x4 &transform, MutableSpan<float3> p
   });
 }
 
-static void copy_transformed_normals(const Span<float3> src,
-                                     const float4x4 &transform,
-                                     MutableSpan<float3> dst)
-{
-  const float3x3 normal_transform = math::transpose(math::invert(float3x3(transform)));
-  if (math::is_equal(normal_transform, float3x3::identity(), 1e-6f)) {
-    dst.copy_from(src);
-  }
-  else {
-    threading::parallel_for(src.index_range(), 1024, [&](const IndexRange range) {
-      for (const int i : range) {
-        dst[i] = normal_transform * src[i];
-      }
-    });
-  }
-}
-
 static void threaded_copy(const GSpan src, GMutableSpan dst)
 {
   BLI_assert(src.size() == dst.size());
@@ -1666,9 +1649,9 @@ static void execute_realize_mesh_task(const RealizeInstancesOptions &options,
     }
     else {
       const IndexRange dst_range = domain_to_range(all_dst_custom_normals.domain);
-      copy_transformed_normals(mesh_info.custom_normal.typed<float3>(),
-                               task.transform,
-                               all_dst_custom_normals.span.typed<float3>().slice(dst_range));
+      math::transform_normals(mesh_info.custom_normal.typed<float3>(),
+                              float3x3(task.transform),
+                              all_dst_custom_normals.span.typed<float3>().slice(dst_range));
     }
   }
 
@@ -2033,8 +2016,9 @@ static void execute_realize_curve_task(const RealizeInstancesOptions &options,
       all_custom_normals.slice(dst_point_range).fill(float3(0, 0, 1));
     }
     else {
-      copy_transformed_normals(
-          curves_info.custom_normal, task.transform, all_custom_normals.slice(dst_point_range));
+      math::transform_normals(curves_info.custom_normal,
+                              float3x3(task.transform),
+                              all_custom_normals.slice(dst_point_range));
     }
   }