Curves: Port fillet node to the new data-block

This commit ports the fillet curves node to the new curves data-block,
and moves the fillet node implementation to the geometry module to help
separate the implementation from the node.

The changes are similar to the subdivide node or resample node. I've
resused common utilities where it makes sense, though some things like
the iteration over attributes can be generalized further. The node
is now multi-threaded per-curve and inside each curve, and some buffers
are reused per curve to avoid many allocations.

The code is more explicit now, and though there is more boilerplate to
pass around many spans, the more complex logic should be more readable.

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

11 files changed, 722 insertions, 607 deletions

diff --git a/source/blender/blenkernel/BKE_attribute.hh b/source/blender/blenkernel/BKE_attribute.hh
index 05ab4f1f1f1..108993d91c0 100644
--- a/source/blender/blenkernel/BKE_attribute.hh
+++ b/source/blender/blenkernel/BKE_attribute.hh
@@ -666,6 +666,22 @@ class MutableAttributeAccessor : public AttributeAccessor {
   void remove_anonymous();
 };
 
+struct AttributeTransferData {
+  /* Expect that if an attribute exists, it is stored as a contiguous array internally anyway. */
+  GVArraySpan src;
+  AttributeMetaData meta_data;
+  bke::GSpanAttributeWriter dst;
+};
+/**
+ * Retrieve attribute arrays and writers for attributes that should be transferred between
+ * data-blocks of the same type.
+ */
+Vector<AttributeTransferData> retrieve_attributes_for_transfer(
+    const bke::AttributeAccessor &src_attributes,
+    bke::MutableAttributeAccessor &dst_attributes,
+    eAttrDomainMask domain_mask,
+    const Set<std::string> &skip = {});
+
 bool allow_procedural_attribute_access(StringRef attribute_name);
 extern const char *no_procedural_access_message;
 
diff --git a/source/blender/blenkernel/intern/attribute_access.cc b/source/blender/blenkernel/intern/attribute_access.cc
index ac1ee19927c..a834b77d65e 100644
--- a/source/blender/blenkernel/intern/attribute_access.cc
+++ b/source/blender/blenkernel/intern/attribute_access.cc
@@ -1011,6 +1011,37 @@ GSpanAttributeWriter MutableAttributeAccessor::lookup_or_add_for_write_only_span
   return {};
 }
 
+Vector<AttributeTransferData> retrieve_attributes_for_transfer(
+    const bke::AttributeAccessor &src_attributes,
+    bke::MutableAttributeAccessor &dst_attributes,
+    const eAttrDomainMask domain_mask,
+    const Set<std::string> &skip)
+{
+  Vector<AttributeTransferData> attributes;
+  src_attributes.for_all(
+      [&](const bke::AttributeIDRef &id, const bke::AttributeMetaData meta_data) {
+        if (!(ATTR_DOMAIN_AS_MASK(meta_data.domain) & domain_mask)) {
+          return true;
+        }
+        if (id.is_named() && skip.contains(id.name())) {
+          return true;
+        }
+        if (!id.should_be_kept()) {
+          return true;
+        }
+
+        GVArray src = src_attributes.lookup(id, meta_data.domain);
+        BLI_assert(src);
+        bke::GSpanAttributeWriter dst = dst_attributes.lookup_or_add_for_write_only_span(
+            id, meta_data.domain, meta_data.data_type);
+        BLI_assert(dst);
+        attributes.append({std::move(src), meta_data, std::move(dst)});
+
+        return true;
+      });
+  return attributes;
+}
+
 }  // namespace blender::bke
 
 /** \} */
diff --git a/source/blender/blenlib/BLI_index_range.hh b/source/blender/blenlib/BLI_index_range.hh
index bd0a7e5bb7a..6fcc560d856 100644
--- a/source/blender/blenlib/BLI_index_range.hh
+++ b/source/blender/blenlib/BLI_index_range.hh
@@ -198,6 +198,16 @@ class IndexRange {
   }
 
   /**
+   * Get the element one before the beginning. The returned value is undefined when the range is
+   * empty, and the range must start after zero already.
+   */
+  constexpr int64_t one_before_start() const
+  {
+    BLI_assert(start_ > 0);
+    return start_ - 1;
+  }
+
+  /**
    * Get the element one after the end. The returned value is undefined when the range is empty.
    */
   constexpr int64_t one_after_last() const
diff --git a/source/blender/blenlib/BLI_math_rotation.hh b/source/blender/blenlib/BLI_math_rotation.hh
index e8b746b34df..50a062162ad 100644
--- a/source/blender/blenlib/BLI_math_rotation.hh
+++ b/source/blender/blenlib/BLI_math_rotation.hh
@@ -15,4 +15,13 @@ namespace blender::math {
  */
 float3 rotate_direction_around_axis(const float3 &direction, const float3 &axis, float angle);
 
+/**
+ * Rotate any arbitrary \a vector around the \a center position, with a unit-length \a axis
+ * and the specified \a angle.
+ */
+float3 rotate_around_axis(const float3 &vector,
+                          const float3 &center,
+                          const float3 &axis,
+                          float angle);
+
 }  // namespace blender::math
diff --git a/source/blender/blenlib/intern/math_rotation.cc b/source/blender/blenlib/intern/math_rotation.cc
index 74300d55954..091e8af85d9 100644
--- a/source/blender/blenlib/intern/math_rotation.cc
+++ b/source/blender/blenlib/intern/math_rotation.cc
@@ -23,4 +23,17 @@ float3 rotate_direction_around_axis(const float3 &direction, const float3 &axis,
   return axis_scaled + diff * std::cos(angle) + cross * std::sin(angle);
 }
 
+float3 rotate_around_axis(const float3 &vector,
+                          const float3 &center,
+                          const float3 &axis,
+                          const float angle)
+
+{
+  float3 result = vector - center;
+  float mat[3][3];
+  axis_angle_normalized_to_mat3(mat, axis, angle);
+  mul_m3_v3(mat, result);
+  return result + center;
+}
+
 }  // namespace blender::math
diff --git a/source/blender/blenlib/tests/BLI_index_range_test.cc b/source/blender/blenlib/tests/BLI_index_range_test.cc
index 10f6784cd44..f5b994d409a 100644
--- a/source/blender/blenlib/tests/BLI_index_range_test.cc
+++ b/source/blender/blenlib/tests/BLI_index_range_test.cc
@@ -105,6 +105,12 @@ TEST(index_range, OneAfterEnd)
   EXPECT_EQ(range.one_after_last(), 8);
 }
 
+TEST(index_range, OneBeforeStart)
+{
+  IndexRange range = IndexRange(5, 3);
+  EXPECT_EQ(range.one_before_start(), 4);
+}
+
 TEST(index_range, Start)
 {
   IndexRange range = IndexRange(6, 2);
diff --git a/source/blender/geometry/CMakeLists.txt b/source/blender/geometry/CMakeLists.txt
index df66a806c16..da83d9e8957 100644
--- a/source/blender/geometry/CMakeLists.txt
+++ b/source/blender/geometry/CMakeLists.txt
@@ -16,6 +16,7 @@ set(INC
 
 set(SRC
   intern/add_curves_on_mesh.cc
+  intern/fillet_curves.cc
   intern/mesh_merge_by_distance.cc
   intern/mesh_primitive_cuboid.cc
   intern/mesh_to_curve_convert.cc
@@ -29,6 +30,7 @@ set(SRC
   intern/uv_parametrizer.c
 
   GEO_add_curves_on_mesh.hh
+  GEO_fillet_curves.hh
   GEO_mesh_merge_by_distance.hh
   GEO_mesh_primitive_cuboid.hh
   GEO_mesh_to_curve.hh
diff --git a/source/blender/geometry/GEO_fillet_curves.hh b/source/blender/geometry/GEO_fillet_curves.hh
new file mode 100644
index 00000000000..1f832f8b6cc
--- /dev/null
+++ b/source/blender/geometry/GEO_fillet_curves.hh
@@ -0,0 +1,23 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "BLI_function_ref.hh"
+#include "BLI_index_mask.hh"
+
+#include "BKE_curves.hh"
+
+namespace blender::geometry {
+
+bke::CurvesGeometry fillet_curves_poly(const bke::CurvesGeometry &src_curves,
+                                       IndexMask curve_selection,
+                                       const VArray<float> &radius,
+                                       const VArray<int> &counts,
+                                       bool limit_radius);
+
+bke::CurvesGeometry fillet_curves_bezier(const bke::CurvesGeometry &src_curves,
+                                         IndexMask curve_selection,
+                                         const VArray<float> &radius,
+                                         bool limit_radius);
+
+}  // namespace blender::geometry
diff --git a/source/blender/geometry/intern/fillet_curves.cc b/source/blender/geometry/intern/fillet_curves.cc
new file mode 100644
index 00000000000..2cca91f40ae
--- /dev/null
+++ b/source/blender/geometry/intern/fillet_curves.cc
@@ -0,0 +1,561 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BKE_attribute_math.hh"
+#include "BKE_curves.hh"
+#include "BKE_curves_utils.hh"
+#include "BKE_geometry_set.hh"
+
+#include "BLI_devirtualize_parameters.hh"
+#include "BLI_math_geom.h"
+#include "BLI_math_rotation.hh"
+#include "BLI_task.hh"
+
+#include "GEO_fillet_curves.hh"
+
+namespace blender::geometry {
+
+/**
+ * Return a range used to retrieve values from an array of values stored per point, but with an
+ * extra element at the end of each curve. This is useful for offsets within curves, where it is
+ * convenient to store the first 0 and have the last offset be the total result curve size.
+ */
+static IndexRange curve_dst_offsets(const IndexRange points, const int curve_index)
+{
+  return {curve_index + points.start(), points.size() + 1};
+}
+
+template<typename T>
+static void threaded_slice_fill(const Span<T> src, const Span<int> offsets, MutableSpan<T> dst)
+{
+  threading::parallel_for(src.index_range(), 512, [&](IndexRange range) {
+    for (const int i : range) {
+      dst.slice(bke::offsets_to_range(offsets, i)).fill(src[i]);
+    }
+  });
+}
+
+template<typename T>
+static void duplicate_fillet_point_data(const bke::CurvesGeometry &src_curves,
+                                        const bke::CurvesGeometry &dst_curves,
+                                        const IndexMask curve_selection,
+                                        const Span<int> point_offsets,
+                                        const Span<T> src,
+                                        MutableSpan<T> dst)
+{
+  threading::parallel_for(curve_selection.index_range(), 512, [&](IndexRange range) {
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange src_points = src_curves.points_for_curve(curve_i);
+      const IndexRange dst_points = dst_curves.points_for_curve(curve_i);
+      const Span<int> offsets = point_offsets.slice(curve_dst_offsets(src_points, curve_i));
+      threaded_slice_fill(src.slice(src_points), offsets, dst.slice(dst_points));
+    }
+  });
+}
+
+static void duplicate_fillet_point_data(const bke::CurvesGeometry &src_curves,
+                                        const bke::CurvesGeometry &dst_curves,
+                                        const IndexMask selection,
+                                        const Span<int> point_offsets,
+                                        const GSpan src,
+                                        GMutableSpan dst)
+{
+  attribute_math::convert_to_static_type(dst.type(), [&](auto dummy) {
+    using T = decltype(dummy);
+    duplicate_fillet_point_data(
+        src_curves, dst_curves, selection, point_offsets, src.typed<T>(), dst.typed<T>());
+  });
+}
+
+static void calculate_result_offsets(const bke::CurvesGeometry &src_curves,
+                                     const IndexMask selection,
+                                     const Span<IndexRange> unselected_ranges,
+                                     const VArray<float> &radii,
+                                     const VArray<int> &counts,
+                                     const Span<bool> cyclic,
+                                     MutableSpan<int> dst_curve_offsets,
+                                     MutableSpan<int> dst_point_offsets)
+{
+  /* Fill the offsets array with the curve point counts, then accumulate them to form offsets. */
+  bke::curves::fill_curve_counts(src_curves, unselected_ranges, dst_curve_offsets);
+  threading::parallel_for(selection.index_range(), 512, [&](IndexRange range) {
+    for (const int curve_i : selection.slice(range)) {
+      const IndexRange src_points = src_curves.points_for_curve(curve_i);
+      const IndexRange offsets_range = curve_dst_offsets(src_points, curve_i);
+
+      MutableSpan<int> point_offsets = dst_point_offsets.slice(offsets_range);
+      MutableSpan<int> point_counts = point_offsets.drop_back(1);
+
+      counts.materialize_compressed(src_points, point_counts);
+      for (int &count : point_counts) {
+        /* Make sure the number of cuts is greater than zero and add one for the existing point. */
+        count = std::max(count, 0) + 1;
+      }
+      if (!cyclic[curve_i]) {
+        /* Endpoints on non-cyclic curves cannot be filleted. */
+        point_counts.first() = 1;
+        point_counts.last() = 1;
+      }
+      /* Implicitly "deselect" points with zero radius. */
+      devirtualize_varray(radii, [&](const auto radii) {
+        for (const int i : IndexRange(src_points.size())) {
+          if (radii[i] == 0.0f) {
+            point_counts[i] = 1;
+          }
+        }
+      });
+
+      bke::curves::accumulate_counts_to_offsets(point_offsets);
+
+      dst_curve_offsets[curve_i] = point_offsets.last();
+    }
+  });
+  bke::curves::accumulate_counts_to_offsets(dst_curve_offsets);
+}
+
+static void calculate_directions(const Span<float3> positions, MutableSpan<float3> directions)
+{
+  for (const int i : positions.index_range().drop_back(1)) {
+    directions[i] = math::normalize(positions[i + 1] - positions[i]);
+  }
+  directions.last() = math::normalize(positions.first() - positions.last());
+}
+
+static void calculate_angles(const Span<float3> directions, MutableSpan<float> angles)
+{
+  angles.first() = M_PI - angle_v3v3(-directions.last(), directions.first());
+  for (const int i : directions.index_range().drop_front(1)) {
+    angles[i] = M_PI - angle_v3v3(-directions[i - 1], directions[i]);
+  }
+}
+
+/**
+ * Find the portion of the previous and next segments used by the current and next point fillets.
+ * If more than the total length of the segment would be used, scale the current point's radius
+ * just enough to make the two points meet in the middle.
+ */
+static float limit_radius(const float3 &position_prev,
+                          const float3 &position,
+                          const float3 &position_next,
+                          const float angle_prev,
+                          const float angle,
+                          const float angle_next,
+                          const float radius_prev,
+                          const float radius,
+                          const float radius_next)
+{
+  const float displacement = radius * std::tan(angle / 2.0f);
+
+  const float displacement_prev = radius_prev * std::tan(angle_prev / 2.0f);
+  const float segment_length_prev = math::distance(position, position_prev);
+  const float total_displacement_prev = displacement_prev + displacement;
+  const float factor_prev = std::clamp(segment_length_prev / total_displacement_prev, 0.0f, 1.0f);
+
+  const float displacement_next = radius_next * std::tan(angle_next / 2.0f);
+  const float segment_length_next = math::distance(position, position_next);
+  const float total_displacement_next = displacement_next + displacement;
+  const float factor_next = std::clamp(segment_length_next / total_displacement_next, 0.0f, 1.0f);
+
+  return radius * std::min(factor_prev, factor_next);
+}
+
+static void limit_radii(const Span<float3> positions,
+                        const Span<float> angles,
+                        const Span<float> radii,
+                        const bool cyclic,
+                        MutableSpan<float> radii_clamped)
+{
+  if (cyclic) {
+    /* First point. */
+    radii_clamped.first() = limit_radius(positions.last(),
+                                         positions.first(),
+                                         positions[1],
+                                         angles.last(),
+                                         angles.first(),
+                                         angles[1],
+                                         radii.last(),
+                                         radii.first(),
+                                         radii[1]);
+    /* All middle points. */
+    for (const int i : positions.index_range().drop_back(1).drop_front(1)) {
+      const int i_prev = i - 1;
+      const int i_next = i + 1;
+      radii_clamped[i] = limit_radius(positions[i_prev],
+                                      positions[i],
+                                      positions[i_next],
+                                      angles[i_prev],
+                                      angles[i],
+                                      angles[i_next],
+                                      radii[i_prev],
+                                      radii[i],
+                                      radii[i_next]);
+    }
+    /* Last point. */
+    radii_clamped.last() = limit_radius(positions.last(1),
+                                        positions.last(),
+                                        positions.first(),
+                                        angles.last(1),
+                                        angles.last(),
+                                        angles.first(),
+                                        radii.last(1),
+                                        radii.last(),
+                                        radii.first());
+  }
+  else {
+    const int i_last = positions.index_range().last();
+    /* First point. */
+    radii_clamped.first() = 0.0f;
+    /* All middle points. */
+    for (const int i : positions.index_range().drop_back(1).drop_front(1)) {
+      const int i_prev = i - 1;
+      const int i_next = i + 1;
+      /* Use a zero radius for the first and last points, because they don't have fillets.
+       * This logic could potentially be unrolled, but it doesn't seem worth it. */
+      const float radius_prev = i_prev == 0 ? 0.0f : radii[i_prev];
+      const float radius_next = i_next == i_last ? 0.0f : radii[i_next];
+      radii_clamped[i] = limit_radius(positions[i_prev],
+                                      positions[i],
+                                      positions[i_next],
+                                      angles[i_prev],
+                                      angles[i],
+                                      angles[i_next],
+                                      radius_prev,
+                                      radii[i],
+                                      radius_next);
+    }
+    /* Last point. */
+    radii_clamped.last() = 0.0f;
+  }
+}
+
+static void calculate_fillet_positions(const Span<float3> src_positions,
+                                       const Span<float> angles,
+                                       const Span<float> radii,
+                                       const Span<float3> directions,
+                                       const Span<int> dst_offsets,
+                                       MutableSpan<float3> dst)
+{
+  const int i_src_last = src_positions.index_range().last();
+  threading::parallel_for(src_positions.index_range(), 512, [&](IndexRange range) {
+    for (const int i_src : range) {
+      const IndexRange arc = bke::offsets_to_range(dst_offsets, i_src);
+      const float3 &src = src_positions[i_src];
+      if (arc.size() == 1) {
+        dst[arc.first()] = src;
+        continue;
+      }
+
+      const int i_src_prev = i_src == 0 ? i_src_last : i_src - 1;
+      const float angle = angles[i_src];
+      const float radius = radii[i_src];
+      const float displacement = radius * std::tan(angle / 2.0f);
+      const float3 prev_dir = -directions[i_src_prev];
+      const float3 &next_dir = directions[i_src];
+      const float3 arc_start = src + prev_dir * displacement;
+      const float3 arc_end = src + next_dir * displacement;
+
+      dst[arc.first()] = arc_start;
+      dst[arc.last()] = arc_end;
+
+      const IndexRange middle = arc.drop_front(1).drop_back(1);
+      if (middle.is_empty()) {
+        continue;
+      }
+
+      const float3 axis = -math::normalize(math::cross(prev_dir, next_dir));
+      const float3 center_direction = math::normalize(math::midpoint(next_dir, prev_dir));
+      const float distance_to_center = std::sqrt(pow2f(radius) + pow2f(displacement));
+      const float3 center = src + center_direction * distance_to_center;
+
+      /* Rotate each middle fillet point around the center. */
+      const float segment_angle = angle / (middle.size() + 1);
+      for (const int i : IndexRange(middle.size())) {
+        const int point_i = middle[i];
+        dst[point_i] = math::rotate_around_axis(arc_start, center, axis, segment_angle * (i + 1));
+      }
+    }
+  });
+}
+
+/**
+ * Set handles for the "Bezier" mode where we rely on setting the inner handles to approximate a
+ * circular arc. The outer (previous and next) handles outside the result fillet segment are set
+ * to vector handles.
+ */
+static void calculate_bezier_handles_bezier_mode(const Span<float3> src_handles_l,
+                                                 const Span<float3> src_handles_r,
+                                                 const Span<int8_t> src_types_l,
+                                                 const Span<int8_t> src_types_r,
+                                                 const Span<float> angles,
+                                                 const Span<float> radii,
+                                                 const Span<float3> directions,
+                                                 const Span<int> dst_offsets,
+                                                 const Span<float3> dst_positions,
+                                                 MutableSpan<float3> dst_handles_l,
+                                                 MutableSpan<float3> dst_handles_r,
+                                                 MutableSpan<int8_t> dst_types_l,
+                                                 MutableSpan<int8_t> dst_types_r)
+{
+  const int i_src_last = src_handles_l.index_range().last();
+  const int i_dst_last = dst_positions.index_range().last();
+  threading::parallel_for(src_handles_l.index_range(), 512, [&](IndexRange range) {
+    for (const int i_src : range) {
+      const IndexRange arc = bke::offsets_to_range(dst_offsets, i_src);
+      if (arc.size() == 1) {
+        dst_handles_l[arc.first()] = src_handles_l[i_src];
+        dst_handles_r[arc.first()] = src_handles_r[i_src];
+        dst_types_l[arc.first()] = src_types_l[i_src];
+        dst_types_r[arc.first()] = src_types_r[i_src];
+        continue;
+      }
+      BLI_assert(arc.size() == 2);
+      const int i_dst_a = arc.first();
+      const int i_dst_b = arc.last();
+
+      const int i_src_prev = i_src == 0 ? i_src_last : i_src - 1;
+      const float angle = angles[i_src];
+      const float radius = radii[i_src];
+      const float3 prev_dir = -directions[i_src_prev];
+      const float3 &next_dir = directions[i_src];
+
+      const float3 &arc_start = dst_positions[arc.first()];
+      const float3 &arc_end = dst_positions[arc.last()];
+
+      /* Calculate the point's handles on the outside of the fillet segment,
+       * connecting to the next or previous result points. */
+      const int i_dst_prev = i_dst_a == 0 ? i_dst_last : i_dst_a - 1;
+      const int i_dst_next = i_dst_b == i_dst_last ? 0 : i_dst_b + 1;
+      dst_handles_l[i_dst_a] = bke::curves::bezier::calculate_vector_handle(
+          dst_positions[i_dst_a], dst_positions[i_dst_prev]);
+      dst_handles_r[i_dst_b] = bke::curves::bezier::calculate_vector_handle(
+          dst_positions[i_dst_b], dst_positions[i_dst_next]);
+      dst_types_l[i_dst_a] = BEZIER_HANDLE_VECTOR;
+      dst_types_r[i_dst_b] = BEZIER_HANDLE_VECTOR;
+
+      /* The inner handles are aligned with the aligned with the outer vector
+       * handles, but have a specific length to best approximate a circle. */
+      const float handle_length = (4.0f / 3.0f) * radius * std::tan(angle / 4.0f);
+      dst_handles_r[i_dst_a] = arc_start - prev_dir * handle_length;
+      dst_handles_l[i_dst_b] = arc_end - next_dir * handle_length;
+      dst_types_r[i_dst_a] = BEZIER_HANDLE_ALIGN;
+      dst_types_l[i_dst_b] = BEZIER_HANDLE_ALIGN;
+    }
+  });
+}
+
+/**
+ * In the poly fillet mode, all the inner handles are set to vector handles, along with the "outer"
+ * (previous and next) handles at each fillet.
+ */
+static void calculate_bezier_handles_poly_mode(const Span<float3> src_handles_l,
+                                               const Span<float3> src_handles_r,
+                                               const Span<int8_t> src_types_l,
+                                               const Span<int8_t> src_types_r,
+                                               const Span<int> dst_offsets,
+                                               const Span<float3> dst_positions,
+                                               MutableSpan<float3> dst_handles_l,
+                                               MutableSpan<float3> dst_handles_r,
+                                               MutableSpan<int8_t> dst_types_l,
+                                               MutableSpan<int8_t> dst_types_r)
+{
+  const int i_dst_last = dst_positions.index_range().last();
+  threading::parallel_for(src_handles_l.index_range(), 512, [&](IndexRange range) {
+    for (const int i_src : range) {
+      const IndexRange arc = bke::offsets_to_range(dst_offsets, i_src);
+      if (arc.size() == 1) {
+        dst_handles_l[arc.first()] = src_handles_l[i_src];
+        dst_handles_r[arc.first()] = src_handles_r[i_src];
+        dst_types_l[arc.first()] = src_types_l[i_src];
+        dst_types_r[arc.first()] = src_types_r[i_src];
+        continue;
+      }
+
+      /* The fillet's next and previous handles are vector handles, as are the inner handles. */
+      dst_types_l.slice(arc).fill(BEZIER_HANDLE_VECTOR);
+      dst_types_r.slice(arc).fill(BEZIER_HANDLE_VECTOR);
+
+      /* Calculate the point's handles on the outside of the fillet segment. This point
+       * won't be selected for a fillet if it is the first or last in a non-cyclic curve. */
+
+      const int i_dst_prev = arc.first() == 0 ? i_dst_last : arc.one_before_start();
+      const int i_dst_next = arc.last() == i_dst_last ? 0 : arc.one_after_last();
+      dst_handles_l[arc.first()] = bke::curves::bezier::calculate_vector_handle(
+          dst_positions[arc.first()], dst_positions[i_dst_prev]);
+      dst_handles_r[arc.last()] = bke::curves::bezier::calculate_vector_handle(
+          dst_positions[arc.last()], dst_positions[i_dst_next]);
+
+      /* Set the values for the inner handles. */
+      const IndexRange middle = arc.drop_front(1).drop_back(1);
+      for (const int i : middle) {
+        dst_handles_r[i] = bke::curves::bezier::calculate_vector_handle(dst_positions[i],
+                                                                        dst_positions[i - 1]);
+        dst_handles_l[i] = bke::curves::bezier::calculate_vector_handle(dst_positions[i],
+                                                                        dst_positions[i + 1]);
+      }
+    }
+  });
+}
+
+static bke::CurvesGeometry fillet_curves(const bke::CurvesGeometry &src_curves,
+                                         const IndexMask curve_selection,
+                                         const VArray<float> &radius_input,
+                                         const VArray<int> &counts,
+                                         const bool limit_radius,
+                                         const bool use_bezier_mode)
+{
+  const Vector<IndexRange> unselected_ranges = curve_selection.extract_ranges_invert(
+      src_curves.curves_range());
+
+  const Span<float3> positions = src_curves.positions();
+  const VArraySpan<bool> cyclic{src_curves.cyclic()};
+  const bke::AttributeAccessor src_attributes = src_curves.attributes();
+
+  bke::CurvesGeometry dst_curves = bke::curves::copy_only_curve_domain(src_curves);
+  /* Stores the offset of every result point for every original point.
+   * The extra length is used in order to store an extra zero for every curve. */
+  Array<int> dst_point_offsets(src_curves.points_num() + src_curves.curves_num());
+  calculate_result_offsets(src_curves,
+                           curve_selection,
+                           unselected_ranges,
+                           radius_input,
+                           counts,
+                           cyclic,
+                           dst_curves.offsets_for_write(),
+                           dst_point_offsets);
+  const Span<int> point_offsets = dst_point_offsets.as_span();
+
+  dst_curves.resize(dst_curves.offsets().last(), dst_curves.curves_num());
+  bke::MutableAttributeAccessor dst_attributes = dst_curves.attributes_for_write();
+  MutableSpan<float3> dst_positions = dst_curves.positions_for_write();
+
+  VArraySpan<int8_t> src_types_l;
+  VArraySpan<int8_t> src_types_r;
+  Span<float3> src_handles_l;
+  Span<float3> src_handles_r;
+  MutableSpan<int8_t> dst_types_l;
+  MutableSpan<int8_t> dst_types_r;
+  MutableSpan<float3> dst_handles_l;
+  MutableSpan<float3> dst_handles_r;
+  if (src_curves.has_curve_with_type(CURVE_TYPE_BEZIER)) {
+    src_types_l = src_curves.handle_types_left();
+    src_types_r = src_curves.handle_types_right();
+    src_handles_l = src_curves.handle_positions_left();
+    src_handles_r = src_curves.handle_positions_right();
+
+    dst_types_l = dst_curves.handle_types_left_for_write();
+    dst_types_r = dst_curves.handle_types_right_for_write();
+    dst_handles_l = dst_curves.handle_positions_left_for_write();
+    dst_handles_r = dst_curves.handle_positions_right_for_write();
+  }
+
+  threading::parallel_for(curve_selection.index_range(), 512, [&](IndexRange range) {
+    Array<float3> directions;
+    Array<float> angles;
+    Array<float> radii;
+    Array<float> input_radii_buffer;
+
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange src_points = src_curves.points_for_curve(curve_i);
+      const Span<int> offsets = point_offsets.slice(curve_dst_offsets(src_points, curve_i));
+      const IndexRange dst_points = dst_curves.points_for_curve(curve_i);
+      const Span<float3> src_positions = positions.slice(src_points);
+
+      directions.reinitialize(src_points.size());
+      calculate_directions(src_positions, directions);
+
+      angles.reinitialize(src_points.size());
+      calculate_angles(directions, angles);
+
+      radii.reinitialize(src_points.size());
+      if (limit_radius) {
+        input_radii_buffer.reinitialize(src_points.size());
+        radius_input.materialize_compressed(src_points, input_radii_buffer);
+        limit_radii(src_positions, angles, input_radii_buffer, cyclic[curve_i], radii);
+      }
+      else {
+        radius_input.materialize_compressed(src_points, radii);
+      }
+
+      calculate_fillet_positions(positions.slice(src_points),
+                                 angles,
+                                 radii,
+                                 directions,
+                                 offsets,
+                                 dst_positions.slice(dst_points));
+
+      if (src_curves.has_curve_with_type(CURVE_TYPE_BEZIER)) {
+        if (use_bezier_mode) {
+          calculate_bezier_handles_bezier_mode(src_handles_l.slice(src_points),
+                                               src_handles_r.slice(src_points),
+                                               src_types_l.slice(src_points),
+                                               src_types_r.slice(src_points),
+                                               angles,
+                                               radii,
+                                               directions,
+                                               offsets,
+                                               dst_positions.slice(dst_points),
+                                               dst_handles_l.slice(dst_points),
+                                               dst_handles_r.slice(dst_points),
+                                               dst_types_l.slice(dst_points),
+                                               dst_types_r.slice(dst_points));
+        }
+        else {
+          calculate_bezier_handles_poly_mode(src_handles_l.slice(src_points),
+                                             src_handles_r.slice(src_points),
+                                             src_types_l.slice(src_points),
+                                             src_types_r.slice(src_points),
+                                             offsets,
+                                             dst_positions.slice(dst_points),
+                                             dst_handles_l.slice(dst_points),
+                                             dst_handles_r.slice(dst_points),
+                                             dst_types_l.slice(dst_points),
+                                             dst_types_r.slice(dst_points));
+        }
+      }
+    }
+  });
+
+  for (auto &attribute : bke::retrieve_attributes_for_transfer(
+           src_attributes,
+           dst_attributes,
+           ATTR_DOMAIN_MASK_POINT,
+           {"position", "handle_type_left", "handle_type_right", "handle_right", "handle_left"})) {
+    duplicate_fillet_point_data(
+        src_curves, dst_curves, curve_selection, point_offsets, attribute.src, attribute.dst.span);
+    attribute.dst.finish();
+  }
+
+  if (!unselected_ranges.is_empty()) {
+    for (auto &attribute : bke::retrieve_attributes_for_transfer(
+             src_attributes, dst_attributes, ATTR_DOMAIN_MASK_POINT)) {
+      bke::curves::copy_point_data(
+          src_curves, dst_curves, unselected_ranges, attribute.src, attribute.dst.span);
+      attribute.dst.finish();
+    }
+  }
+
+  return dst_curves;
+}
+
+bke::CurvesGeometry fillet_curves_poly(const bke::CurvesGeometry &src_curves,
+                                       const IndexMask curve_selection,
+                                       const VArray<float> &radius,
+                                       const VArray<int> &count,
+                                       const bool limit_radius)
+{
+  return fillet_curves(src_curves, curve_selection, radius, count, limit_radius, false);
+}
+
+bke::CurvesGeometry fillet_curves_bezier(const bke::CurvesGeometry &src_curves,
+                                         const IndexMask curve_selection,
+                                         const VArray<float> &radius,
+                                         const bool limit_radius)
+{
+  return fillet_curves(src_curves,
+                       curve_selection,
+                       radius,
+                       VArray<int>::ForSingle(1, src_curves.points_num()),
+                       limit_radius,
+                       true);
+}
+
+}  // namespace blender::geometry
diff --git a/source/blender/geometry/intern/subdivide_curves.cc b/source/blender/geometry/intern/subdivide_curves.cc
index b6476d19818..8a6f3cbd5e3 100644
--- a/source/blender/geometry/intern/subdivide_curves.cc
+++ b/source/blender/geometry/intern/subdivide_curves.cc
@@ -56,40 +56,6 @@ static void calculate_result_offsets(const bke::CurvesGeometry &src_curves,
   bke::curves::accumulate_counts_to_offsets(dst_curve_offsets);
 }
 
-struct AttributeTransferData {
-  /* Expect that if an attribute exists, it is stored as a contiguous array internally anyway. */
-  GVArraySpan src;
-  bke::GSpanAttributeWriter dst;
-};
-
-static Vector<AttributeTransferData> retrieve_point_attributes(
-    const bke::AttributeAccessor &src_attributes,
-    bke::MutableAttributeAccessor &dst_attributes,
-    const Set<std::string> &skip = {})
-{
-  Vector<AttributeTransferData> attributes;
-  src_attributes.for_all(
-      [&](const bke::AttributeIDRef &id, const bke::AttributeMetaData meta_data) {
-        if (meta_data.domain != ATTR_DOMAIN_POINT) {
-          /* Curve domain attributes are all copied directly to the result in one step. */
-          return true;
-        }
-        if (id.is_named() && skip.contains(id.name())) {
-          return true;
-        }
-
-        GVArray src = src_attributes.lookup(id, ATTR_DOMAIN_POINT);
-        BLI_assert(src);
-        bke::GSpanAttributeWriter dst = dst_attributes.lookup_or_add_for_write_only_span(
-            id, ATTR_DOMAIN_POINT, meta_data.data_type);
-        BLI_assert(dst);
-        attributes.append({std::move(src), std::move(dst)});
-
-        return true;
-      });
-  return attributes;
-}
-
 template<typename T>
 static inline void linear_interpolation(const T &a, const T &b, MutableSpan<T> dst)
 {
@@ -365,7 +331,8 @@ bke::CurvesGeometry subdivide_curves(const bke::CurvesGeometry &src_curves,
   bke::MutableAttributeAccessor dst_attributes = dst_curves.attributes_for_write();
 
   auto subdivide_catmull_rom = [&](IndexMask selection) {
-    for (auto &attribute : retrieve_point_attributes(src_attributes, dst_attributes)) {
+    for (auto &attribute : bke::retrieve_attributes_for_transfer(
+             src_attributes, dst_attributes, ATTR_DOMAIN_MASK_POINT)) {
       subdivide_attribute_catmull_rom(src_curves,
                                       dst_curves,
                                       selection,
@@ -378,7 +345,8 @@ bke::CurvesGeometry subdivide_curves(const bke::CurvesGeometry &src_curves,
   };
 
   auto subdivide_poly = [&](IndexMask selection) {
-    for (auto &attribute : retrieve_point_attributes(src_attributes, dst_attributes)) {
+    for (auto &attribute : bke::retrieve_attributes_for_transfer(
+             src_attributes, dst_attributes, ATTR_DOMAIN_MASK_POINT)) {
       subdivide_attribute_linear(
           src_curves, dst_curves, selection, point_offsets, attribute.src, attribute.dst.span);
       attribute.dst.finish();
@@ -419,13 +387,14 @@ bke::CurvesGeometry subdivide_curves(const bke::CurvesGeometry &src_curves,
       }
     });
 
-    for (auto &attribute : retrieve_point_attributes(src_attributes,
-                                                     dst_attributes,
-                                                     {"position",
-                                                      "handle_type_left",
-                                                      "handle_type_right",
-                                                      "handle_right",
-                                                      "handle_left"})) {
+    for (auto &attribute : bke::retrieve_attributes_for_transfer(src_attributes,
+                                                                 dst_attributes,
+                                                                 ATTR_DOMAIN_MASK_POINT,
+                                                                 {"position",
+                                                                  "handle_type_left",
+                                                                  "handle_type_right",
+                                                                  "handle_right",
+                                                                  "handle_left"})) {
       subdivide_attribute_linear(
           src_curves, dst_curves, selection, point_offsets, attribute.src, attribute.dst.span);
       attribute.dst.finish();
@@ -445,7 +414,8 @@ bke::CurvesGeometry subdivide_curves(const bke::CurvesGeometry &src_curves,
                                      subdivide_nurbs);
 
   if (!unselected_ranges.is_empty()) {
-    for (auto &attribute : retrieve_point_attributes(src_attributes, dst_attributes)) {
+    for (auto &attribute : bke::retrieve_attributes_for_transfer(
+             src_attributes, dst_attributes, ATTR_DOMAIN_MASK_POINT)) {
       bke::curves::copy_point_data(
           src_curves, dst_curves, unselected_ranges, attribute.src, attribute.dst.span);
       attribute.dst.finish();
diff --git a/source/blender/nodes/geometry/nodes/node_geo_curve_fillet.cc b/source/blender/nodes/geometry/nodes/node_geo_curve_fillet.cc
index dd8471d2dac..ab1f8269c39 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_curve_fillet.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_curve_fillet.cc
@@ -1,17 +1,12 @@
 /* SPDX-License-Identifier: GPL-2.0-or-later */
 
-#include "BLI_task.hh"
-
 #include "UI_interface.h"
 #include "UI_resources.h"
 
-#include "DNA_node_types.h"
+#include "GEO_fillet_curves.hh"
 
 #include "node_geometry_util.hh"
 
-#include "BKE_curves.hh"
-#include "BKE_spline.hh"
-
 namespace blender::nodes::node_geo_curve_fillet_cc {
 
 NODE_STORAGE_FUNCS(NodeGeometryCurveFillet)
@@ -45,574 +40,18 @@ static void node_layout(uiLayout *layout, bContext *UNUSED(C), PointerRNA *ptr)
 static void node_init(bNodeTree *UNUSED(tree), bNode *node)
 {
   NodeGeometryCurveFillet *data = MEM_cnew<NodeGeometryCurveFillet>(__func__);
-
   data->mode = GEO_NODE_CURVE_FILLET_BEZIER;
   node->storage = data;
 }
 
-struct FilletParam {
-  GeometryNodeCurveFilletMode mode;
-
-  /* Number of points to be added. */
-  VArray<int> counts;
-
-  /* Radii for fillet arc at all vertices. */
-  VArray<float> radii;
-
-  /* Whether or not fillets are allowed to overlap. */
-  bool limit_radius;
-};
-
-/* A data structure used to store fillet data about all vertices to be filleted. */
-struct FilletData {
-  Span<float3> positions;
-  Array<float3> directions, axes;
-  Array<float> radii, angles;
-  Array<int> counts;
-};
-
 static void node_update(bNodeTree *ntree, bNode *node)
 {
   const NodeGeometryCurveFillet &storage = node_storage(*node);
   const GeometryNodeCurveFilletMode mode = (GeometryNodeCurveFilletMode)storage.mode;
-
   bNodeSocket *poly_socket = ((bNodeSocket *)node->inputs.first)->next;
-
   nodeSetSocketAvailability(ntree, poly_socket, mode == GEO_NODE_CURVE_FILLET_POLY);
 }
 
-/* Function to get the center of a fillet. */
-static float3 get_center(const float3 vec_pos2prev,
-                         const float3 pos,
-                         const float3 axis,
-                         const float angle)
-{
-  float3 vec_pos2center;
-  rotate_normalized_v3_v3v3fl(vec_pos2center, vec_pos2prev, axis, M_PI_2 - angle / 2.0f);
-  vec_pos2center *= 1.0f / sinf(angle / 2.0f);
-
-  return vec_pos2center + pos;
-}
-
-/* Function to get the center of the fillet using fillet data */
-static float3 get_center(const float3 vec_pos2prev, const FilletData &fd, const int index)
-{
-  const float angle = fd.angles[index];
-  const float3 axis = fd.axes[index];
-  const float3 pos = fd.positions[index];
-
-  return get_center(vec_pos2prev, pos, axis, angle);
-}
-
-/* Calculate the direction vectors from each vertex to their previous vertex. */
-static Array<float3> calculate_directions(const Span<float3> positions)
-{
-  const int num = positions.size();
-  Array<float3> directions(num);
-
-  for (const int i : IndexRange(num - 1)) {
-    directions[i] = math::normalize(positions[i + 1] - positions[i]);
-  }
-  directions[num - 1] = math::normalize(positions[0] - positions[num - 1]);
-
-  return directions;
-}
-
-/* Calculate the axes around which the fillet is built. */
-static Array<float3> calculate_axes(const Span<float3> directions)
-{
-  const int num = directions.size();
-  Array<float3> axes(num);
-
-  axes[0] = math::normalize(math::cross(-directions[num - 1], directions[0]));
-  for (const int i : IndexRange(1, num - 1)) {
-    axes[i] = math::normalize(math::cross(-directions[i - 1], directions[i]));
-  }
-
-  return axes;
-}
-
-/* Calculate the angle of the arc formed by the fillet. */
-static Array<float> calculate_angles(const Span<float3> directions)
-{
-  const int num = directions.size();
-  Array<float> angles(num);
-
-  angles[0] = M_PI - angle_v3v3(-directions[num - 1], directions[0]);
-  for (const int i : IndexRange(1, num - 1)) {
-    angles[i] = M_PI - angle_v3v3(-directions[i - 1], directions[i]);
-  }
-
-  return angles;
-}
-
-/* Calculate the segment count in each filleted arc. */
-static Array<int> calculate_counts(const FilletParam &fillet_param,
-                                   const int num,
-                                   const int spline_offset,
-                                   const bool cyclic)
-{
-  Array<int> counts(num, 1);
-  if (fillet_param.mode == GEO_NODE_CURVE_FILLET_POLY) {
-    for (const int i : IndexRange(num)) {
-      counts[i] = fillet_param.counts[spline_offset + i];
-    }
-  }
-  if (!cyclic) {
-    counts[0] = counts[num - 1] = 0;
-  }
-
-  return counts;
-}
-
-/* Calculate the radii for the vertices to be filleted. */
-static Array<float> calculate_radii(const FilletParam &fillet_param,
-                                    const int num,
-                                    const int spline_offset)
-{
-  Array<float> radii(num, 0.0f);
-  if (fillet_param.limit_radius) {
-    for (const int i : IndexRange(num)) {
-      radii[i] = std::max(fillet_param.radii[spline_offset + i], 0.0f);
-    }
-  }
-  else {
-    for (const int i : IndexRange(num)) {
-      radii[i] = fillet_param.radii[spline_offset + i];
-    }
-  }
-
-  return radii;
-}
-
-/* Calculate the number of vertices added per vertex on the source spline. */
-static int calculate_point_counts(MutableSpan<int> point_counts,
-                                  const Span<float> radii,
-                                  const Span<int> counts)
-{
-  int added_count = 0;
-  for (const int i : IndexRange(point_counts.size())) {
-    /* Calculate number of points to be added for the vertex. */
-    if (radii[i] != 0.0f) {
-      added_count += counts[i];
-      point_counts[i] = counts[i] + 1;
-    }
-  }
-
-  return added_count;
-}
-
-static FilletData calculate_fillet_data(const Spline &spline,
-                                        const FilletParam &fillet_param,
-                                        int &added_count,
-                                        MutableSpan<int> point_counts,
-                                        const int spline_offset)
-{
-  const int num = spline.size();
-
-  FilletData fd;
-  fd.directions = calculate_directions(spline.positions());
-  fd.positions = spline.positions();
-  fd.axes = calculate_axes(fd.directions);
-  fd.angles = calculate_angles(fd.directions);
-  fd.counts = calculate_counts(fillet_param, num, spline_offset, spline.is_cyclic());
-  fd.radii = calculate_radii(fillet_param, num, spline_offset);
-
-  added_count = calculate_point_counts(point_counts, fd.radii, fd.counts);
-
-  return fd;
-}
-
-/* Limit the radius based on angle and radii to prevent overlapping. */
-static void limit_radii(FilletData &fd, const bool cyclic)
-{
-  MutableSpan<float> radii(fd.radii);
-  Span<float> angles(fd.angles);
-  Span<float3> positions(fd.positions);
-
-  const int num = radii.size();
-  const int fillet_count = cyclic ? num : num - 2;
-  const int start = cyclic ? 0 : 1;
-  Array<float> max_radii(num, FLT_MAX);
-
-  if (cyclic) {
-    /* Calculate lengths between adjacent control points. */
-    const float len_prev = math::distance(positions[0], positions[num - 1]);
-    const float len_next = math::distance(positions[0], positions[1]);
-
-    /* Calculate tangent lengths of fillets in control points. */
-    const float tan_len = radii[0] * tan(angles[0] / 2.0f);
-    const float tan_len_prev = radii[num - 1] * tan(angles[num - 1] / 2.0f);
-    const float tan_len_next = radii[1] * tan(angles[1] / 2.0f);
-
-    float factor_prev = 1.0f, factor_next = 1.0f;
-    if (tan_len + tan_len_prev > len_prev) {
-      factor_prev = len_prev / (tan_len + tan_len_prev);
-    }
-    if (tan_len + tan_len_next > len_next) {
-      factor_next = len_next / (tan_len + tan_len_next);
-    }
-
-    /* Scale max radii by calculated factors. */
-    max_radii[0] = radii[0] * std::min(factor_next, factor_prev);
-    max_radii[1] = radii[1] * factor_next;
-    max_radii[num - 1] = radii[num - 1] * factor_prev;
-  }
-
-  /* Initialize max_radii to largest possible radii. */
-  float prev_dist = math::distance(positions[1], positions[0]);
-  for (const int i : IndexRange(1, num - 2)) {
-    const float temp_dist = math::distance(positions[i], positions[i + 1]);
-    max_radii[i] = std::min(prev_dist, temp_dist) / tan(angles[i] / 2.0f);
-    prev_dist = temp_dist;
-  }
-
-  /* Max radii calculations for each index. */
-  for (const int i : IndexRange(start, fillet_count - 1)) {
-    const float len_next = math::distance(positions[i], positions[i + 1]);
-    const float tan_len = radii[i] * tan(angles[i] / 2.0f);
-    const float tan_len_next = radii[i + 1] * tan(angles[i + 1] / 2.0f);
-
-    /* Scale down radii if too large for segment. */
-    float factor = 1.0f;
-    if (tan_len + tan_len_next > len_next) {
-      factor = len_next / (tan_len + tan_len_next);
-    }
-    max_radii[i] = std::min(max_radii[i], radii[i] * factor);
-    max_radii[i + 1] = std::min(max_radii[i + 1], radii[i + 1] * factor);
-  }
-
-  /* Assign the max_radii to the fillet data's radii. */
-  for (const int i : IndexRange(num)) {
-    radii[i] = std::min(radii[i], max_radii[i]);
-  }
-}
-
-/*
- * Create a mapping from each vertex in the destination spline to that of the source spline.
- * Used for copying the data from the source spline.
- */
-static Array<int> create_dst_to_src_map(const Span<int> point_counts, const int total_points)
-{
-  Array<int> map(total_points);
-  MutableSpan<int> map_span{map};
-  int index = 0;
-
-  for (const int i : point_counts.index_range()) {
-    map_span.slice(index, point_counts[i]).fill(i);
-    index += point_counts[i];
-  }
-
-  BLI_assert(index == total_points);
-
-  return map;
-}
-
-template<typename T>
-static void copy_attribute_by_mapping(const Span<T> src,
-                                      MutableSpan<T> dst,
-                                      const Span<int> mapping)
-{
-  for (const int i : dst.index_range()) {
-    dst[i] = src[mapping[i]];
-  }
-}
-
-/* Copy radii and tilts from source spline to destination. Positions are handled later in update
- * positions methods. */
-static void copy_common_attributes_by_mapping(const Spline &src,
-                                              Spline &dst,
-                                              const Span<int> mapping)
-{
-  copy_attribute_by_mapping(src.radii(), dst.radii(), mapping);
-  copy_attribute_by_mapping(src.tilts(), dst.tilts(), mapping);
-
-  src.attributes.foreach_attribute(
-      [&](const AttributeIDRef &attribute_id, const AttributeMetaData &meta_data) {
-        std::optional<GSpan> src_attribute = src.attributes.get_for_read(attribute_id);
-        if (dst.attributes.create(attribute_id, meta_data.data_type)) {
-          std::optional<GMutableSpan> dst_attribute = dst.attributes.get_for_write(attribute_id);
-          if (dst_attribute) {
-            attribute_math::convert_to_static_type(dst_attribute->type(), [&](auto dummy) {
-              using T = decltype(dummy);
-              copy_attribute_by_mapping(
-                  src_attribute->typed<T>(), dst_attribute->typed<T>(), mapping);
-            });
-            return true;
-          }
-        }
-        BLI_assert_unreachable();
-        return false;
-      },
-      ATTR_DOMAIN_POINT);
-}
-
-/* Update the vertex positions and handle positions of a Bezier spline based on fillet data. */
-static void update_bezier_positions(const FilletData &fd,
-                                    BezierSpline &dst_spline,
-                                    const BezierSpline &src_spline,
-                                    const Span<int> point_counts)
-{
-  Span<float> radii(fd.radii);
-  Span<float> angles(fd.angles);
-  Span<float3> axes(fd.axes);
-  Span<float3> positions(fd.positions);
-  Span<float3> directions(fd.directions);
-
-  const int num = radii.size();
-
-  int i_dst = 0;
-  for (const int i_src : IndexRange(num)) {
-    const int count = point_counts[i_src];
-
-    /* Skip if the point count for the vertex is 1. */
-    if (count == 1) {
-      dst_spline.positions()[i_dst] = src_spline.positions()[i_src];
-      dst_spline.handle_types_left()[i_dst] = src_spline.handle_types_left()[i_src];
-      dst_spline.handle_types_right()[i_dst] = src_spline.handle_types_right()[i_src];
-      dst_spline.handle_positions_left()[i_dst] = src_spline.handle_positions_left()[i_src];
-      dst_spline.handle_positions_right()[i_dst] = src_spline.handle_positions_right()[i_src];
-      i_dst++;
-      continue;
-    }
-
-    /* Calculate the angle to be formed between any 2 adjacent vertices within the fillet. */
-    const float segment_angle = angles[i_src] / (count - 1);
-    /* Calculate the handle length for each added vertex. Equation: L = 4R/3 * tan(A/4) */
-    const float handle_length = 4.0f * radii[i_src] / 3.0f * tan(segment_angle / 4.0f);
-    /* Calculate the distance by which each vertex should be displaced from their initial position.
-     */
-    const float displacement = radii[i_src] * tan(angles[i_src] / 2.0f);
-
-    /* Position the end points of the arc and their handles. */
-    const int end_i = i_dst + count - 1;
-    const float3 prev_dir = i_src == 0 ? -directions[num - 1] : -directions[i_src - 1];
-    const float3 next_dir = directions[i_src];
-    dst_spline.positions()[i_dst] = positions[i_src] + displacement * prev_dir;
-    dst_spline.positions()[end_i] = positions[i_src] + displacement * next_dir;
-    dst_spline.handle_positions_right()[i_dst] = dst_spline.positions()[i_dst] -
-                                                 handle_length * prev_dir;
-    dst_spline.handle_positions_left()[end_i] = dst_spline.positions()[end_i] -
-                                                handle_length * next_dir;
-    dst_spline.handle_types_right()[i_dst] = dst_spline.handle_types_left()[end_i] =
-        BEZIER_HANDLE_ALIGN;
-    dst_spline.handle_types_left()[i_dst] = dst_spline.handle_types_right()[end_i] =
-        BEZIER_HANDLE_VECTOR;
-    dst_spline.mark_cache_invalid();
-
-    /* Calculate the center of the radius to be formed. */
-    const float3 center = get_center(dst_spline.positions()[i_dst] - positions[i_src], fd, i_src);
-    /* Calculate the vector of the radius formed by the first vertex. */
-    float3 radius_vec = dst_spline.positions()[i_dst] - center;
-    float radius;
-    radius_vec = math::normalize_and_get_length(radius_vec, radius);
-
-    dst_spline.handle_types_right().slice(1, count - 2).fill(BEZIER_HANDLE_ALIGN);
-    dst_spline.handle_types_left().slice(1, count - 2).fill(BEZIER_HANDLE_ALIGN);
-
-    /* For each of the vertices in between the end points. */
-    for (const int j : IndexRange(1, count - 2)) {
-      int index = i_dst + j;
-      /* Rotate the radius by the segment angle and determine its tangent (used for getting handle
-       * directions). */
-      float3 new_radius_vec, tangent_vec;
-      rotate_normalized_v3_v3v3fl(new_radius_vec, radius_vec, -axes[i_src], segment_angle);
-      rotate_normalized_v3_v3v3fl(tangent_vec, new_radius_vec, axes[i_src], M_PI_2);
-      radius_vec = new_radius_vec;
-      tangent_vec *= handle_length;
-
-      /* Adjust the positions of the respective vertex and its handles. */
-      dst_spline.positions()[index] = center + new_radius_vec * radius;
-      dst_spline.handle_positions_left()[index] = dst_spline.positions()[index] + tangent_vec;
-      dst_spline.handle_positions_right()[index] = dst_spline.positions()[index] - tangent_vec;
-    }
-
-    i_dst += count;
-  }
-}
-
-/* Update the vertex positions of a Poly spline based on fillet data. */
-static void update_poly_positions(const FilletData &fd,
-                                  Spline &dst_spline,
-                                  const Spline &src_spline,
-                                  const Span<int> point_counts)
-{
-  Span<float> radii(fd.radii);
-  Span<float> angles(fd.angles);
-  Span<float3> axes(fd.axes);
-  Span<float3> positions(fd.positions);
-  Span<float3> directions(fd.directions);
-
-  const int num = radii.size();
-
-  int i_dst = 0;
-  for (const int i_src : IndexRange(num)) {
-    const int count = point_counts[i_src];
-
-    /* Skip if the point count for the vertex is 1. */
-    if (count == 1) {
-      dst_spline.positions()[i_dst] = src_spline.positions()[i_src];
-      i_dst++;
-      continue;
-    }
-
-    const float segment_angle = angles[i_src] / (count - 1);
-    const float displacement = radii[i_src] * tan(angles[i_src] / 2.0f);
-
-    /* Position the end points of the arc. */
-    const int end_i = i_dst + count - 1;
-    const float3 prev_dir = i_src == 0 ? -directions[num - 1] : -directions[i_src - 1];
-    const float3 next_dir = directions[i_src];
-    dst_spline.positions()[i_dst] = positions[i_src] + displacement * prev_dir;
-    dst_spline.positions()[end_i] = positions[i_src] + displacement * next_dir;
-
-    /* Calculate the center of the radius to be formed. */
-    const float3 center = get_center(dst_spline.positions()[i_dst] - positions[i_src], fd, i_src);
-    /* Calculate the vector of the radius formed by the first vertex. */
-    float3 radius_vec = dst_spline.positions()[i_dst] - center;
-
-    for (const int j : IndexRange(1, count - 2)) {
-      /* Rotate the radius by the segment angle */
-      float3 new_radius_vec;
-      rotate_normalized_v3_v3v3fl(new_radius_vec, radius_vec, -axes[i_src], segment_angle);
-      radius_vec = new_radius_vec;
-
-      dst_spline.positions()[i_dst + j] = center + new_radius_vec;
-    }
-
-    i_dst += count;
-  }
-}
-
-static SplinePtr fillet_spline(const Spline &spline,
-                               const FilletParam &fillet_param,
-                               const int spline_offset)
-{
-  const int num = spline.size();
-  const bool cyclic = spline.is_cyclic();
-
-  if (num < 3) {
-    return spline.copy();
-  }
-
-  /* Initialize the point_counts with 1s (at least one vertex on dst for each vertex on src). */
-  Array<int> point_counts(num, 1);
-
-  int added_count = 0;
-  /* Update point_counts array and added_count. */
-  FilletData fd = calculate_fillet_data(
-      spline, fillet_param, added_count, point_counts, spline_offset);
-  if (fillet_param.limit_radius) {
-    limit_radii(fd, cyclic);
-  }
-
-  const int total_points = added_count + num;
-  const Array<int> dst_to_src = create_dst_to_src_map(point_counts, total_points);
-  SplinePtr dst_spline_ptr = spline.copy_only_settings();
-  (*dst_spline_ptr).resize(total_points);
-  copy_common_attributes_by_mapping(spline, *dst_spline_ptr, dst_to_src);
-
-  switch (spline.type()) {
-    case CURVE_TYPE_BEZIER: {
-      const BezierSpline &src_spline = static_cast<const BezierSpline &>(spline);
-      BezierSpline &dst_spline = static_cast<BezierSpline &>(*dst_spline_ptr);
-      if (fillet_param.mode == GEO_NODE_CURVE_FILLET_POLY) {
-        dst_spline.handle_types_left().fill(BEZIER_HANDLE_VECTOR);
-        dst_spline.handle_types_right().fill(BEZIER_HANDLE_VECTOR);
-        update_poly_positions(fd, dst_spline, src_spline, point_counts);
-      }
-      else {
-        update_bezier_positions(fd, dst_spline, src_spline, point_counts);
-      }
-      break;
-    }
-    case CURVE_TYPE_POLY: {
-      update_poly_positions(fd, *dst_spline_ptr, spline, point_counts);
-      break;
-    }
-    case CURVE_TYPE_NURBS: {
-      const NURBSpline &src_spline = static_cast<const NURBSpline &>(spline);
-      NURBSpline &dst_spline = static_cast<NURBSpline &>(*dst_spline_ptr);
-      copy_attribute_by_mapping(src_spline.weights(), dst_spline.weights(), dst_to_src);
-      update_poly_positions(fd, dst_spline, src_spline, point_counts);
-      break;
-    }
-    case CURVE_TYPE_CATMULL_ROM: {
-      BLI_assert_unreachable();
-      break;
-    }
-  }
-
-  return dst_spline_ptr;
-}
-
-static std::unique_ptr<CurveEval> fillet_curve(const CurveEval &input_curve,
-                                               const FilletParam &fillet_param)
-{
-  Span<SplinePtr> input_splines = input_curve.splines();
-
-  std::unique_ptr<CurveEval> output_curve = std::make_unique<CurveEval>();
-  const int splines_num = input_splines.size();
-  output_curve->resize(splines_num);
-  MutableSpan<SplinePtr> output_splines = output_curve->splines();
-  Array<int> spline_offsets = input_curve.control_point_offsets();
-
-  threading::parallel_for(input_splines.index_range(), 128, [&](IndexRange range) {
-    for (const int i : range) {
-      output_splines[i] = fillet_spline(*input_splines[i], fillet_param, spline_offsets[i]);
-    }
-  });
-  output_curve->attributes = input_curve.attributes;
-
-  return output_curve;
-}
-
-static void calculate_curve_fillet(GeometrySet &geometry_set,
-                                   const GeometryNodeCurveFilletMode mode,
-                                   const Field<float> &radius_field,
-                                   const std::optional<Field<int>> &count_field,
-                                   const bool limit_radius)
-{
-  if (!geometry_set.has_curves()) {
-    return;
-  }
-
-  FilletParam fillet_param;
-  fillet_param.mode = mode;
-
-  CurveComponent &component = geometry_set.get_component_for_write<CurveComponent>();
-  GeometryComponentFieldContext field_context{component, ATTR_DOMAIN_POINT};
-  const int domain_size = component.attribute_domain_size(ATTR_DOMAIN_POINT);
-  fn::FieldEvaluator field_evaluator{field_context, domain_size};
-
-  field_evaluator.add(radius_field);
-
-  if (mode == GEO_NODE_CURVE_FILLET_POLY) {
-    field_evaluator.add(*count_field);
-  }
-
-  field_evaluator.evaluate();
-
-  fillet_param.radii = field_evaluator.get_evaluated<float>(0);
-  if (fillet_param.radii.is_single() && fillet_param.radii.get_internal_single() < 0.0f) {
-    return;
-  }
-
-  if (mode == GEO_NODE_CURVE_FILLET_POLY) {
-    fillet_param.counts = field_evaluator.get_evaluated<int>(1);
-  }
-
-  fillet_param.limit_radius = limit_radius;
-
-  const Curves &src_curves_id = *geometry_set.get_curves_for_read();
-  const std::unique_ptr<CurveEval> input_curve = curves_to_curve_eval(*component.get_for_read());
-  std::unique_ptr<CurveEval> output_curve = fillet_curve(*input_curve, fillet_param);
-
-  Curves *dst_curves_id = curve_eval_to_curves(*output_curve);
-  bke::curves_copy_parameters(src_curves_id, *dst_curves_id);
-  geometry_set.replace_curves(dst_curves_id);
-}
-
 static void node_geo_exec(GeoNodeExecParams params)
 {
   GeometrySet geometry_set = params.extract_input<GeometrySet>("Curve");
@@ -629,7 +68,42 @@ static void node_geo_exec(GeoNodeExecParams params)
   }
 
   geometry_set.modify_geometry_sets([&](GeometrySet &geometry_set) {
-    calculate_curve_fillet(geometry_set, mode, radius_field, count_field, limit_radius);
+    if (!geometry_set.has_curves()) {
+      return;
+    }
+
+    const CurveComponent &component = *geometry_set.get_component_for_read<CurveComponent>();
+    const Curves &curves_id = *component.get_for_read();
+    const bke::CurvesGeometry &curves = bke::CurvesGeometry::wrap(curves_id.geometry);
+    GeometryComponentFieldContext context{component, ATTR_DOMAIN_POINT};
+    fn::FieldEvaluator evaluator{context, curves.points_num()};
+    evaluator.add(radius_field);
+
+    switch (mode) {
+      case GEO_NODE_CURVE_FILLET_BEZIER: {
+        evaluator.evaluate();
+        bke::CurvesGeometry dst_curves = geometry::fillet_curves_bezier(
+            curves, curves.curves_range(), evaluator.get_evaluated<float>(0), limit_radius);
+        Curves *dst_curves_id = bke::curves_new_nomain(std::move(dst_curves));
+        bke::curves_copy_parameters(curves_id, *dst_curves_id);
+        geometry_set.replace_curves(dst_curves_id);
+        break;
+      }
+      case GEO_NODE_CURVE_FILLET_POLY: {
+        evaluator.add(*count_field);
+        evaluator.evaluate();
+        bke::CurvesGeometry dst_curves = geometry::fillet_curves_poly(
+            curves,
+            curves.curves_range(),
+            evaluator.get_evaluated<float>(0),
+            evaluator.get_evaluated<int>(1),
+            limit_radius);
+        Curves *dst_curves_id = bke::curves_new_nomain(std::move(dst_curves));
+        bke::curves_copy_parameters(curves_id, *dst_curves_id);
+        geometry_set.replace_curves(dst_curves_id);
+        break;
+      }
+    }
   });
 
   params.set_output("Curve", std::move(geometry_set));