Curves: make tangent computation more robust

Previously, when there were multiple curve points at the same or
almost the same position, the computed tangent was unpredictable.

Now, the handling of this case is more explicit, making it more
predictable. In general, when there are duplicate points, it will just use
tangents from neighboring points now.

Also fixes T98209.

Differential Revision: https://developer.blender.org/D15016

2 files changed, 71 insertions, 8 deletions

diff --git a/source/blender/blenkernel/intern/curve_poly.cc b/source/blender/blenkernel/intern/curve_poly.cc
index 2db7cd71ad3..2a546e81825 100644
--- a/source/blender/blenkernel/intern/curve_poly.cc
+++ b/source/blender/blenkernel/intern/curve_poly.cc
@@ -13,15 +13,28 @@
 
 namespace blender::bke::curves::poly {
 
-static float3 direction_bisect(const float3 &prev, const float3 &middle, const float3 &next)
+static float3 direction_bisect(const float3 &prev,
+                               const float3 &middle,
+                               const float3 &next,
+                               bool &r_used_fallback)
 {
+  const float epsilon = 1e-6f;
+  const bool prev_equal = math::almost_equal_relative(prev, middle, epsilon);
+  const bool next_equal = math::almost_equal_relative(middle, next, epsilon);
+  if (prev_equal && next_equal) {
+    r_used_fallback = true;
+    return {0.0f, 0.0f, 0.0f};
+  }
+  if (prev_equal) {
+    return math::normalize(next - middle);
+  }
+  if (next_equal) {
+    return math::normalize(middle - prev);
+  }
+
   const float3 dir_prev = math::normalize(middle - prev);
   const float3 dir_next = math::normalize(next - middle);
-
   const float3 result = math::normalize(dir_prev + dir_next);
-  if (UNLIKELY(math::is_zero(result))) {
-    return float3(0.0f, 0.0f, 1.0f);
-  }
   return result;
 }
 
@@ -36,8 +49,11 @@ void calculate_tangents(const Span<float3> positions,
     return;
   }
 
+  bool used_fallback = false;
+
   for (const int i : IndexRange(1, positions.size() - 2)) {
-    tangents[i] = direction_bisect(positions[i - 1], positions[i], positions[i + 1]);
+    tangents[i] = direction_bisect(
+        positions[i - 1], positions[i], positions[i + 1], used_fallback);
   }
 
   if (is_cyclic) {
@@ -45,13 +61,46 @@ void calculate_tangents(const Span<float3> positions,
     const float3 &last = positions.last();
     const float3 &first = positions.first();
     const float3 &second = positions[1];
-    tangents.first() = direction_bisect(last, first, second);
-    tangents.last() = direction_bisect(second_to_last, last, first);
+    tangents.first() = direction_bisect(last, first, second, used_fallback);
+    tangents.last() = direction_bisect(second_to_last, last, first, used_fallback);
   }
   else {
     tangents.first() = math::normalize(positions[1] - positions.first());
     tangents.last() = math::normalize(positions.last() - positions[positions.size() - 2]);
   }
+
+  if (!used_fallback) {
+    return;
+  }
+
+  /* Find the first tangent that does not use the fallback. */
+  int first_valid_tangent_index = -1;
+  for (const int i : tangents.index_range()) {
+    if (!math::is_zero(tangents[i])) {
+      first_valid_tangent_index = i;
+      break;
+    }
+  }
+  if (first_valid_tangent_index == -1) {
+    /* If all tangents used the fallback, it means that all positions are (almost) the same. Just
+     * use the up-vector as default tangent. */
+    const float3 up_vector{0.0f, 0.0f, 1.0f};
+    tangents.fill(up_vector);
+  }
+  else {
+    const float3 &first_valid_tangent = tangents[first_valid_tangent_index];
+    /* If the first few tangents are invalid, use the tangent from the first point with a valid
+     * tangent. */
+    tangents.take_front(first_valid_tangent_index).fill(first_valid_tangent);
+    /* Use the previous valid tangent for points that had no valid tangent. */
+    for (const int i : tangents.index_range().drop_front(first_valid_tangent_index + 1)) {
+      float3 &tangent = tangents[i];
+      if (math::is_zero(tangent)) {
+        const float3 &prev_tangent = tangents[i - 1];
+        tangent = prev_tangent;
+      }
+    }
+  }
 }
 
 void calculate_normals_z_up(const Span<float3> tangents, MutableSpan<float3> normals)
diff --git a/source/blender/blenlib/BLI_math_vector.hh b/source/blender/blenlib/BLI_math_vector.hh
index b9f0939674e..7983bbccb35 100644
--- a/source/blender/blenlib/BLI_math_vector.hh
+++ b/source/blender/blenlib/BLI_math_vector.hh
@@ -49,6 +49,20 @@ template<typename T, int Size> inline bool is_any_zero(const vec_base<T, Size> &
   return false;
 }
 
+template<typename T, int Size>
+inline bool almost_equal_relative(const vec_base<T, Size> &a,
+                                  const vec_base<T, Size> &b,
+                                  const T &epsilon_factor)
+{
+  for (int i = 0; i < Size; i++) {
+    const float epsilon = epsilon_factor * math::abs(a[i]);
+    if (math::distance(a[i], b[i]) > epsilon) {
+      return false;
+    }
+  }
+  return true;
+}
+
 template<typename T, int Size> inline vec_base<T, Size> abs(const vec_base<T, Size> &a)
 {
   vec_base<T, Size> result;