Curves: Port tangent and normal calculation to the new data-block

Port the "Normal" and "Curve Tangent" nodes to the new curves data-block
to avoid the conversion to `CurveEval`. This should make them faster by
avoiding all that copying, but otherwise nothing else has changed.

This also includes a fix to move the normal mode as a built-in curve
attribute when converting to and from `CurveEval`. The attribute is
needed because the option is used implicitly in many nodes currently.

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

1 files changed, 154 insertions, 0 deletions

diff --git a/source/blender/blenkernel/intern/curve_poly.cc b/source/blender/blenkernel/intern/curve_poly.cc
new file mode 100644
index 00000000000..b0ed62d38dd
--- /dev/null
+++ b/source/blender/blenkernel/intern/curve_poly.cc
@@ -0,0 +1,154 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+/** \file
+ * \ingroup bke
+ */
+
+#include <algorithm>
+
+#include "BLI_math_vector.h"
+#include "BLI_math_vector.hh"
+
+#include "BKE_curves.hh"
+
+namespace blender::bke::curves::poly {
+
+static float3 direction_bisect(const float3 &prev, const float3 &middle, const float3 &next)
+{
+  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;
+}
+
+void calculate_tangents(const Span<float3> positions,
+                        const bool is_cyclic,
+                        MutableSpan<float3> tangents)
+{
+  BLI_assert(positions.size() == tangents.size());
+
+  if (positions.size() == 1) {
+    tangents.first() = float3(0.0f, 0.0f, 1.0f);
+    return;
+  }
+
+  for (const int i : IndexRange(1, positions.size() - 2)) {
+    tangents[i] = direction_bisect(positions[i - 1], positions[i], positions[i + 1]);
+  }
+
+  if (is_cyclic) {
+    const float3 &second_to_last = positions[positions.size() - 2];
+    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);
+  }
+  else {
+    tangents.first() = math::normalize(positions[1] - positions.first());
+    tangents.last() = math::normalize(positions.last() - positions[positions.size() - 2]);
+  }
+}
+
+static float3 rotate_direction_around_axis(const float3 &direction,
+                                           const float3 &axis,
+                                           const float angle)
+{
+  BLI_ASSERT_UNIT_V3(direction);
+  BLI_ASSERT_UNIT_V3(axis);
+
+  const float3 axis_scaled = axis * math::dot(direction, axis);
+  const float3 diff = direction - axis_scaled;
+  const float3 cross = math::cross(axis, diff);
+
+  return axis_scaled + diff * std::cos(angle) + cross * std::sin(angle);
+}
+
+void calculate_normals_z_up(const Span<float3> tangents, MutableSpan<float3> normals)
+{
+  BLI_assert(normals.size() == tangents.size());
+
+  /* Same as in `vec_to_quat`. */
+  const float epsilon = 1e-4f;
+  for (const int i : normals.index_range()) {
+    const float3 &tangent = tangents[i];
+    if (std::abs(tangent.x) + std::abs(tangent.y) < epsilon) {
+      normals[i] = {1.0f, 0.0f, 0.0f};
+    }
+    else {
+      normals[i] = math::normalize(float3(tangent.y, -tangent.x, 0.0f));
+    }
+  }
+}
+
+/**
+ * Rotate the last normal in the same way the tangent has been rotated.
+ */
+static float3 calculate_next_normal(const float3 &last_normal,
+                                    const float3 &last_tangent,
+                                    const float3 &current_tangent)
+{
+  if (math::is_zero(last_tangent) || math::is_zero(current_tangent)) {
+    return last_normal;
+  }
+  const float angle = angle_normalized_v3v3(last_tangent, current_tangent);
+  if (angle != 0.0) {
+    const float3 axis = math::normalize(math::cross(last_tangent, current_tangent));
+    return rotate_direction_around_axis(last_normal, axis, angle);
+  }
+  return last_normal;
+}
+
+void calculate_normals_minimum(const Span<float3> tangents,
+                               const bool cyclic,
+                               MutableSpan<float3> normals)
+{
+  BLI_assert(normals.size() == tangents.size());
+
+  if (normals.is_empty()) {
+    return;
+  }
+
+  const float epsilon = 1e-4f;
+
+  /* Set initial normal. */
+  const float3 &first_tangent = tangents.first();
+  if (fabs(first_tangent.x) + fabs(first_tangent.y) < epsilon) {
+    normals.first() = {1.0f, 0.0f, 0.0f};
+  }
+  else {
+    normals.first() = math::normalize(float3(first_tangent.y, -first_tangent.x, 0.0f));
+  }
+
+  /* Forward normal with minimum twist along the entire spline. */
+  for (const int i : IndexRange(1, normals.size() - 1)) {
+    normals[i] = calculate_next_normal(normals[i - 1], tangents[i - 1], tangents[i]);
+  }
+
+  if (!cyclic) {
+    return;
+  }
+
+  /* Compute how much the first normal deviates from the normal that has been forwarded along the
+   * entire cyclic spline. */
+  const float3 uncorrected_last_normal = calculate_next_normal(
+      normals.last(), tangents.last(), tangents.first());
+  float correction_angle = angle_signed_on_axis_v3v3_v3(
+      normals.first(), uncorrected_last_normal, tangents.first());
+  if (correction_angle > M_PI) {
+    correction_angle = correction_angle - 2 * M_PI;
+  }
+
+  /* Gradually apply correction by rotating all normals slightly. */
+  const float angle_step = correction_angle / normals.size();
+  for (const int i : normals.index_range()) {
+    const float angle = angle_step * i;
+    normals[i] = rotate_direction_around_axis(normals[i], tangents[i], angle);
+  }
+}
+
+}  // namespace blender::bke::curves::poly