Curves: Add length cache, length paramerterize utility

This commit adds calculation of lengths along the curve for each
evaluated point. This is used for sampling, resampling, the "curve
parameter" node, and potentially more places in the future.

This commit also includes a utility for calculation of uniform samples
in blenlib. It can find evenlyspaced samples along a sequence of points
and use linear interpolation to move data from those points to the
samples. Making the utility more general aligns better with the more
functional approach of the new curves code and makes the behavior
available elsewhere.

A "color math" header is added to allow very basic interpolation
between two colors in the `blender::math` namespace.

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

7 files changed, 480 insertions, 0 deletions

diff --git a/source/blender/blenkernel/BKE_curves.hh b/source/blender/blenkernel/BKE_curves.hh
index 82f77d83bec..05a20917552 100644
--- a/source/blender/blenkernel/BKE_curves.hh
+++ b/source/blender/blenkernel/BKE_curves.hh
@@ -77,6 +77,15 @@ class CurvesGeometryRuntime {
   mutable std::mutex position_cache_mutex;
   mutable bool position_cache_dirty = true;
 
+  /**
+   * Cache of lengths along each evaluated curve for for each evaluated point. If a curve is
+   * cyclic, it needs one more length value to correspond to the last segment, so in order to
+   * make slicing this array for a curve fast, an extra float is stored for every curve.
+   */
+  mutable Vector<float> evaluated_length_cache;
+  mutable std::mutex length_cache_mutex;
+  mutable bool length_cache_dirty = true;
+
   /** Direction of the spline at each evaluated point. */
   mutable Vector<float3> evaluated_tangents_cache;
   mutable std::mutex tangent_cache_mutex;
@@ -267,6 +276,20 @@ class CurvesGeometry : public ::CurvesGeometry {
   Span<float3> evaluated_positions() const;
 
   /**
+   * Return a cache of accumulated lengths along the curve. Each item is the length of the
+   * subsequent segment (the first value is the length of the first segment rather than 0).
+   * This calculation is rather trivial, and only depends on the evaluated positions, but
+   * the results are used often, and it is necessarily single threaded per curve, so it is cached.
+   *
+   * \param cyclic: This argument is redundant with the data stored for the curve,
+   * but is passed for performance reasons to avoid looking up the attribute.
+   */
+  Span<float> evaluated_lengths_for_curve(int curve_index, bool cyclic) const;
+
+  /** Calculates the data described by #evaluated_lengths_for_curve if necessary. */
+  void ensure_evaluated_lengths() const;
+
+  /**
    * Evaluate a generic data to the standard evaluated points of a specific curve,
    * defined by the resolution attribute or other factors, depending on the curve type.
    *
@@ -281,6 +304,9 @@ class CurvesGeometry : public ::CurvesGeometry {
    */
   void ensure_nurbs_basis_cache() const;
 
+  /** Return the slice of #evaluated_length_cache that corresponds to this curve index. */
+  IndexRange lengths_range_for_curve(int curve_index, bool cyclic) const;
+
   /* --------------------------------------------------------------------
    * Operations.
    */
diff --git a/source/blender/blenkernel/intern/curves_geometry.cc b/source/blender/blenkernel/intern/curves_geometry.cc
index 7ceaa8f0f37..c4e9a06aad0 100644
--- a/source/blender/blenkernel/intern/curves_geometry.cc
+++ b/source/blender/blenkernel/intern/curves_geometry.cc
@@ -11,6 +11,7 @@
 
 #include "BLI_bounds.hh"
 #include "BLI_index_mask_ops.hh"
+#include "BLI_length_parameterize.hh"
 
 #include "DNA_curves_types.h"
 
@@ -721,6 +722,63 @@ void CurvesGeometry::interpolate_to_evaluated(const int curve_index,
   BLI_assert_unreachable();
 }
 
+IndexRange CurvesGeometry::lengths_range_for_curve(const int curve_index, const bool cyclic) const
+{
+  BLI_assert(cyclic == this->cyclic()[curve_index]);
+  const IndexRange points = this->evaluated_points_for_curve(curve_index);
+  const int start = points.start() + curve_index;
+  const int size = curves::curve_segment_size(points.size(), cyclic);
+  return {start, size};
+}
+
+void CurvesGeometry::ensure_evaluated_lengths() const
+{
+  if (!this->runtime->length_cache_dirty) {
+    return;
+  }
+
+  /* A double checked lock. */
+  std::scoped_lock lock{this->runtime->length_cache_mutex};
+  if (!this->runtime->length_cache_dirty) {
+    return;
+  }
+
+  threading::isolate_task([&]() {
+    /* Use an extra length value for the final cyclic segment for a consistent size
+     * (see comment on #evaluated_length_cache). */
+    const int total_size = this->evaluated_points_num() + this->curves_num();
+    this->runtime->evaluated_length_cache.resize(total_size);
+    MutableSpan<float> evaluated_lengths = this->runtime->evaluated_length_cache;
+
+    Span<float3> evaluated_positions = this->evaluated_positions();
+    VArray<bool> curves_cyclic = this->cyclic();
+
+    threading::parallel_for(this->curves_range(), 128, [&](IndexRange curves_range) {
+      for (const int curve_index : curves_range) {
+        const bool cyclic = curves_cyclic[curve_index];
+        const IndexRange evaluated_points = this->evaluated_points_for_curve(curve_index);
+        if (UNLIKELY(evaluated_points.is_empty())) {
+          continue;
+        }
+        const IndexRange lengths_range = this->lengths_range_for_curve(curve_index, cyclic);
+        length_parameterize::accumulate_lengths(evaluated_positions.slice(evaluated_points),
+                                                cyclic,
+                                                evaluated_lengths.slice(lengths_range));
+      }
+    });
+  });
+
+  this->runtime->length_cache_dirty = false;
+}
+
+Span<float> CurvesGeometry::evaluated_lengths_for_curve(const int curve_index,
+                                                        const bool cyclic) const
+{
+  BLI_assert(!this->runtime->length_cache_dirty);
+  const IndexRange range = this->lengths_range_for_curve(curve_index, cyclic);
+  return this->runtime->evaluated_length_cache.as_span().slice(range);
+}
+
 /** \} */
 
 /* -------------------------------------------------------------------- */
@@ -747,6 +805,7 @@ void CurvesGeometry::tag_positions_changed()
   this->runtime->position_cache_dirty = true;
   this->runtime->tangent_cache_dirty = true;
   this->runtime->normal_cache_dirty = true;
+  this->runtime->length_cache_dirty = true;
 }
 void CurvesGeometry::tag_topology_changed()
 {
@@ -755,6 +814,7 @@ void CurvesGeometry::tag_topology_changed()
   this->runtime->normal_cache_dirty = true;
   this->runtime->offsets_cache_dirty = true;
   this->runtime->nurbs_basis_cache_dirty = true;
+  this->runtime->length_cache_dirty = true;
 }
 void CurvesGeometry::tag_normals_changed()
 {
diff --git a/source/blender/blenlib/BLI_length_parameterize.hh b/source/blender/blenlib/BLI_length_parameterize.hh
new file mode 100644
index 00000000000..e4a4e9cbb9c
--- /dev/null
+++ b/source/blender/blenlib/BLI_length_parameterize.hh
@@ -0,0 +1,80 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math_base.hh"
+#include "BLI_math_color.hh"
+#include "BLI_math_vector.hh"
+#include "BLI_vector.hh"
+
+namespace blender::length_parameterize {
+
+/**
+ * Return the size of the necessary lengths array for a group of points, taking into account the
+ * possible last cyclic segment.
+ *
+ * \note This is the same as #bke::curves::curve_segment_size.
+ */
+inline int lengths_num(const int points_num, const bool cyclic)
+{
+  return cyclic ? points_num : points_num - 1;
+}
+
+/**
+ * Accumulate the length of the next segment into each point.
+ */
+template<typename T>
+void accumulate_lengths(const Span<T> values, const bool cyclic, MutableSpan<float> lengths)
+{
+  BLI_assert(lengths.size() == lengths_num(values.size(), cyclic));
+  float length = 0.0f;
+  for (const int i : IndexRange(values.size() - 1)) {
+    length += math::distance(values[i], values[i + 1]);
+    lengths[i] = length;
+  }
+  if (cyclic) {
+    lengths.last() = length + math::distance(values.last(), values.first());
+  }
+}
+
+template<typename T>
+void linear_interpolation(const Span<T> src,
+                          const Span<int> indices,
+                          const Span<float> factors,
+                          MutableSpan<T> dst)
+{
+  BLI_assert(indices.size() == factors.size());
+  BLI_assert(indices.size() == dst.size());
+  const int last_src_index = src.index_range().last();
+
+  int cyclic_sample_count = 0;
+  for (int i = indices.index_range().last(); i > 0; i--) {
+    if (indices[i] != last_src_index) {
+      break;
+    }
+    dst[i] = math::interpolate(src.last(), src.first(), factors[i]);
+    cyclic_sample_count++;
+  }
+
+  for (const int i : dst.index_range().drop_back(cyclic_sample_count)) {
+    dst[i] = math::interpolate(src[indices[i]], src[indices[i] + 1], factors[i]);
+  }
+}
+
+/**
+ * Find the given number of points, evenly spaced along the provided length. For non-cyclic
+ * sequences, the first point will always be included, and last point will always be included if
+ * the #count is greater than zero. For cyclic sequences, the first point will always be included.
+ *
+ * \warning The #count argument must be greater than zero.
+ */
+void create_uniform_samples(Span<float> lengths,
+                            bool cyclic,
+                            MutableSpan<int> indices,
+                            MutableSpan<float> factors);
+
+}  // namespace blender::length_parameterize
diff --git a/source/blender/blenlib/BLI_math_color.hh b/source/blender/blenlib/BLI_math_color.hh
new file mode 100644
index 00000000000..5195cfb6238
--- /dev/null
+++ b/source/blender/blenlib/BLI_math_color.hh
@@ -0,0 +1,28 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2022 Blender Foundation. */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include <cmath>
+#include <type_traits>
+
+#include "BLI_color.hh"
+#include "BLI_math_base.hh"
+
+namespace blender::math {
+
+inline ColorGeometry4f interpolate(const ColorGeometry4f &a,
+                                   const ColorGeometry4f &b,
+                                   const float t)
+{
+  return {math::interpolate(a.r, b.r, t),
+          math::interpolate(a.g, b.g, t),
+          math::interpolate(a.b, b.b, t),
+          math::interpolate(a.a, b.a, t)};
+}
+
+}  // namespace blender::math
diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt
index 647726722b1..99e07264276 100644
--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -80,6 +80,7 @@ set(SRC
   intern/kdtree_4d.c
   intern/lasso_2d.c
   intern/listbase.c
+  intern/length_parameterize.cc
   intern/math_base.c
   intern/math_base_inline.c
   intern/math_base_safe_inline.c
@@ -225,6 +226,7 @@ set(SRC
   BLI_kdtree.h
   BLI_kdtree_impl.h
   BLI_lasso_2d.h
+  BLI_length_parameterize.hh
   BLI_linear_allocator.hh
   BLI_link_utils.h
   BLI_linklist.h
@@ -241,6 +243,7 @@ set(SRC
   BLI_math_bits.h
   BLI_math_boolean.hh
   BLI_math_color.h
+  BLI_math_color.hh
   BLI_math_color_blend.h
   BLI_math_geom.h
   BLI_math_inline.h
@@ -435,6 +438,7 @@ if(WITH_GTESTS)
     tests/BLI_index_range_test.cc
     tests/BLI_inplace_priority_queue_test.cc
     tests/BLI_kdopbvh_test.cc
+    tests/BLI_length_parameterize_test.cc
     tests/BLI_linear_allocator_test.cc
     tests/BLI_linklist_lockfree_test.cc
     tests/BLI_listbase_test.cc
diff --git a/source/blender/blenlib/intern/length_parameterize.cc b/source/blender/blenlib/intern/length_parameterize.cc
new file mode 100644
index 00000000000..4947740ba8c
--- /dev/null
+++ b/source/blender/blenlib/intern/length_parameterize.cc
@@ -0,0 +1,80 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BLI_length_parameterize.hh"
+
+namespace blender::length_parameterize {
+
+void create_uniform_samples(const Span<float> lengths,
+                            const bool cyclic,
+                            MutableSpan<int> indices,
+                            MutableSpan<float> factors)
+{
+  const int count = indices.size();
+  BLI_assert(count > 0);
+  BLI_assert(lengths.size() >= 1);
+  BLI_assert(std::is_sorted(lengths.begin(), lengths.end()));
+  const int segments_num = lengths.size();
+  const int points_num = cyclic ? segments_num : segments_num + 1;
+
+  indices.first() = 0;
+  factors.first() = 0.0f;
+  if (count == 1) {
+    return;
+  }
+
+  const float total_length = lengths.last();
+  if (total_length == 0.0f) {
+    indices.fill(0);
+    factors.fill(0.0f);
+    return;
+  }
+
+  const float step_length = total_length / (count - (cyclic ? 0 : 1));
+  const float step_length_inv = 1.0f / step_length;
+
+  int i_dst = 1;
+  /* Store the length at the previous point in a variable so it can start out at zero
+   * (the lengths array doesn't contain 0 for the first point). */
+  float prev_length = 0.0f;
+  for (const int i_src : IndexRange(points_num - 1)) {
+    const float next_length = lengths[i_src];
+    const float segment_length = next_length - prev_length;
+    if (segment_length == 0.0f) {
+      continue;
+    }
+    /* Add every sample that fits in this segment. */
+    const float segment_length_inv = 1.0f / segment_length;
+    const int segment_samples_num = std::ceil(next_length * step_length_inv - i_dst);
+    indices.slice(i_dst, segment_samples_num).fill(i_src);
+
+    for (const int i : factors.index_range().slice(i_dst, segment_samples_num)) {
+      const float length_in_segment = step_length * i - prev_length;
+      factors[i] = length_in_segment * segment_length_inv;
+    }
+
+    i_dst += segment_samples_num;
+
+    prev_length = next_length;
+  }
+
+  /* Add the samples on the last cyclic segment if necessary, and also the samples
+   * that weren't created in the previous loop due to floating point inacuracy. */
+  if (cyclic && lengths.size() > 1) {
+    indices.drop_front(i_dst).fill(points_num - 1);
+    const float segment_length = lengths.last() - lengths.last(1);
+    if (segment_length == 0.0f) {
+      return;
+    }
+    const float segment_length_inv = 1.0f / segment_length;
+    for (const int i : indices.index_range().drop_front(i_dst)) {
+      const float length_in_segment = step_length * i - prev_length;
+      factors[i] = length_in_segment * segment_length_inv;
+    }
+  }
+  else {
+    indices.drop_front(i_dst).fill(points_num - 2);
+    factors.drop_front(i_dst).fill(1.0f);
+  }
+}
+
+}  // namespace blender::length_parameterize
diff --git a/source/blender/blenlib/tests/BLI_length_parameterize_test.cc b/source/blender/blenlib/tests/BLI_length_parameterize_test.cc
new file mode 100644
index 00000000000..4a8b7095888
--- /dev/null
+++ b/source/blender/blenlib/tests/BLI_length_parameterize_test.cc
@@ -0,0 +1,202 @@
+/* SPDX-License-Identifier: Apache-2.0 */
+
+#include "BLI_array.hh"
+#include "BLI_length_parameterize.hh"
+#include "BLI_vector.hh"
+
+#include "testing/testing.h"
+
+namespace blender::length_parameterize::tests {
+
+template<typename T> Array<float> calculate_lengths(const Span<T> values, const bool cyclic)
+{
+  Array<float> lengths(lengths_num(values.size(), cyclic));
+  accumulate_lengths<T>(values, cyclic, lengths);
+  return lengths;
+}
+
+template<typename T> void test_uniform_lengths(const Span<T> values)
+{
+  const float segment_length = math::distance(values.first(), values.last()) / (values.size() - 1);
+  for (const int i : values.index_range().drop_back(1)) {
+    EXPECT_NEAR(math::distance(values[i], values[i + 1]), segment_length, 1e-5);
+  }
+}
+
+TEST(length_parameterize, FloatSimple)
+{
+  Array<float> values{{0, 1, 4}};
+  Array<float> lengths = calculate_lengths(values.as_span(), false);
+
+  Array<int> indices(4);
+  Array<float> factors(4);
+  create_uniform_samples(lengths, false, indices, factors);
+  Array<float> results(4);
+  linear_interpolation<float>(values, indices, factors, results);
+  Array<float> expected({
+      0.0f,
+      1.33333f,
+      2.66667f,
+      4.0f,
+  });
+  for (const int i : results.index_range()) {
+    EXPECT_NEAR(results[i], expected[i], 1e-5);
+  }
+  test_uniform_lengths(results.as_span());
+}
+
+TEST(length_parameterize, Float)
+{
+  Array<float> values{{1, 2, 3, 5, 10}};
+  Array<float> lengths = calculate_lengths(values.as_span(), false);
+
+  Array<int> indices(20);
+  Array<float> factors(20);
+  create_uniform_samples(lengths, false, indices, factors);
+  Array<float> results(20);
+  linear_interpolation<float>(values, indices, factors, results);
+  Array<float> expected({
+      1.0f,     1.47368f, 1.94737f, 2.42105f, 2.89474f, 3.36842f, 3.84211f,
+      4.31579f, 4.78947f, 5.26316f, 5.73684f, 6.21053f, 6.68421f, 7.1579f,
+      7.63158f, 8.10526f, 8.57895f, 9.05263f, 9.52632f, 10.0f,
+  });
+  for (const int i : results.index_range()) {
+    EXPECT_NEAR(results[i], expected[i], 1e-5);
+  }
+  test_uniform_lengths(results.as_span());
+}
+
+TEST(length_parameterize, Float2)
+{
+  Array<float2> values{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}};
+  Array<float> lengths = calculate_lengths(values.as_span(), false);
+
+  Array<int> indices(12);
+  Array<float> factors(12);
+  create_uniform_samples(lengths, false, indices, factors);
+  Array<float2> results(12);
+  linear_interpolation<float2>(values, indices, factors, results);
+  Array<float2> expected({
+      {0.0f, 0.0f},
+      {0.272727f, 0.0f},
+      {0.545455f, 0.0f},
+      {0.818182f, 0.0f},
+      {1.0f, 0.0909091f},
+      {1.0f, 0.363636f},
+      {1.0f, 0.636364f},
+      {1.0f, 0.909091f},
+      {0.818182f, 1.0f},
+      {0.545455f, 1.0f},
+      {0.272727f, 1.0f},
+      {0.0f, 1.0f},
+  });
+  for (const int i : results.index_range()) {
+    EXPECT_NEAR(results[i].x, expected[i].x, 1e-5);
+    EXPECT_NEAR(results[i].y, expected[i].y, 1e-5);
+  }
+}
+
+TEST(length_parameterize, Float2Cyclic)
+{
+  Array<float2> values{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}};
+  Array<float> lengths = calculate_lengths(values.as_span(), true);
+
+  Array<int> indices(12);
+  Array<float> factors(12);
+  create_uniform_samples(lengths, true, indices, factors);
+  Array<float2> results(12);
+  linear_interpolation<float2>(values, indices, factors, results);
+  Array<float2> expected({
+      {0.0f, 0.0f},
+      {0.333333f, 0.0f},
+      {0.666667f, 0.0f},
+      {1.0f, 0.0f},
+      {1.0f, 0.333333f},
+      {1.0f, 0.666667f},
+      {1.0f, 1.0f},
+      {0.666667f, 1.0f},
+      {0.333333f, 1.0f},
+      {0.0f, 1.0f},
+      {0.0f, 0.666667f},
+      {0.0f, 0.333333f},
+  });
+  for (const int i : results.index_range()) {
+    EXPECT_NEAR(results[i].x, expected[i].x, 1e-5);
+    EXPECT_NEAR(results[i].y, expected[i].y, 1e-5);
+  }
+}
+
+TEST(length_parameterize, LineMany)
+{
+  Array<float> values{{1, 2}};
+  Array<float> lengths = calculate_lengths(values.as_span(), false);
+
+  Array<int> indices(5007);
+  Array<float> factors(5007);
+  create_uniform_samples(lengths, false, indices, factors);
+  Array<float> results(5007);
+  linear_interpolation<float>(values, indices, factors, results);
+  Array<float> expected({
+      1.9962f, 1.9964f, 1.9966f, 1.9968f, 1.997f, 1.9972f, 1.9974f, 1.9976f, 1.9978f, 1.998f,
+      1.9982f, 1.9984f, 1.9986f, 1.9988f, 1.999f, 1.9992f, 1.9994f, 1.9996f, 1.9998f, 2.0f,
+  });
+  for (const int i : expected.index_range()) {
+    EXPECT_NEAR(results.as_span().take_back(20)[i], expected[i], 1e-5);
+  }
+}
+
+TEST(length_parameterize, CyclicMany)
+{
+  Array<float2> values{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}};
+  Array<float> lengths = calculate_lengths(values.as_span(), true);
+
+  Array<int> indices(5007);
+  Array<float> factors(5007);
+  create_uniform_samples(lengths, true, indices, factors);
+  Array<float2> results(5007);
+  linear_interpolation<float2>(values, indices, factors, results);
+  Array<float2> expected({
+      {0, 0.0159776},  {0, 0.0151787},  {0, 0.0143797},  {0, 0.013581},   {0, 0.0127821},
+      {0, 0.0119832},  {0, 0.0111842},  {0, 0.0103855},  {0, 0.00958657}, {0, 0.00878763},
+      {0, 0.00798869}, {0, 0.00718999}, {0, 0.00639105}, {0, 0.00559211}, {0, 0.00479317},
+      {0, 0.00399446}, {0, 0.00319552}, {0, 0.00239658}, {0, 0.00159764}, {0, 0.000798941},
+  });
+  for (const int i : expected.index_range()) {
+    EXPECT_NEAR(results.as_span().take_back(20)[i].x, expected[i].x, 1e-5);
+    EXPECT_NEAR(results.as_span().take_back(20)[i].y, expected[i].y, 1e-5);
+  }
+}
+
+TEST(length_parameterize, InterpolateColor)
+{
+  Array<float2> values{{{0, 0}, {1, 0}, {1, 1}, {0, 1}}};
+  Array<float> lengths = calculate_lengths(values.as_span(), true);
+
+  Array<ColorGeometry4f> colors{{{0, 0, 0, 1}, {1, 0, 0, 1}, {1, 1, 0, 1}, {0, 1, 0, 1}}};
+
+  Array<int> indices(10);
+  Array<float> factors(10);
+  create_uniform_samples(lengths, true, indices, factors);
+  Array<ColorGeometry4f> results(10);
+  linear_interpolation<ColorGeometry4f>(colors, indices, factors, results);
+  Array<ColorGeometry4f> expected({
+      {0, 0, 0, 1},
+      {0.4, 0, 0, 1},
+      {0.8, 0, 0, 1},
+      {1, 0.2, 0, 1},
+      {1, 0.6, 0, 1},
+      {1, 1, 0, 1},
+      {0.6, 1, 0, 1},
+      {0.2, 1, 0, 1},
+      {0, 0.8, 0, 1},
+      {0, 0.4, 0, 1},
+  });
+  for (const int i : results.index_range()) {
+    EXPECT_NEAR(results[i].r, expected[i].r, 1e-6);
+    EXPECT_NEAR(results[i].g, expected[i].g, 1e-6);
+    EXPECT_NEAR(results[i].b, expected[i].b, 1e-6);
+    EXPECT_NEAR(results[i].a, expected[i].a, 1e-6);
+  }
+}
+
+}  // namespace blender::length_parameterize::tests