Geometry Nodes: Initial basic curve data support

This patch adds initial curve support to geometry nodes. Currently
there is only one node available, the "Curve to Mesh" node, T87428.

However, the aim of the changes here is larger than just supporting
curve data in nodes-- it also uses the opportunity to add better spline
data structures, intended to replace the existing curve evaluation code.
The curve code in Blender is quite old, and it's generally regarded as
some of the messiest, hardest-to-understand code as well. The classes
in `BKE_spline.hh` aim to be faster, more extensible, and much more
easily understandable. Further explanation can be found in comments in
that file.

Initial builtin spline attributes are supported-- reading and writing
from the `cyclic` and `resolution` attributes works with any of the
attribute nodes. Also, only Z-up normal calculation is implemented
at the moment, and tilts do not apply yet.

**Limitations**
 - For now, you must bring curves into the node tree with an "Object
   Info" node. Changes to the curve modifier stack will come later.
 - Converting to a mesh is necessary to visualize the curve data.

Further progress can be tracked in: T87245
Higher level design document: https://wiki.blender.org/wiki/Modules/Physics_Nodes/Projects/EverythingNodes/CurveNodes

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

11 files changed, 1829 insertions, 5 deletions

diff --git a/source/blender/blenkernel/intern/DerivedMesh.cc b/source/blender/blenkernel/intern/DerivedMesh.cc
index c8416ac897e..fcdfeef5db4 100644
--- a/source/blender/blenkernel/intern/DerivedMesh.cc
+++ b/source/blender/blenkernel/intern/DerivedMesh.cc
@@ -1859,9 +1859,11 @@ static void editbmesh_calc_modifiers(struct Depsgraph *depsgraph,
     BKE_id_free(nullptr, mesh_orco);
   }
 
-  /* Ensure normals calculation below is correct. */
-  BLI_assert((mesh_input->flag & ME_AUTOSMOOTH) == (mesh_final->flag & ME_AUTOSMOOTH));
-  BLI_assert(mesh_input->smoothresh == mesh_final->smoothresh);
+  /* Ensure normals calculation below is correct (normal settings have transferred properly).
+   * However, nodes modifiers might create meshes from scratch or transfer meshes from other
+   * objects with different settings, and in general it doesn't make sense to guarentee that
+   * the settings are the same as the original mesh. If necessary, this could become a modifier
+   * type flag. */
   BLI_assert(mesh_input->smoothresh == mesh_cage->smoothresh);
 
   /* Compute normals. */
diff --git a/source/blender/blenkernel/intern/attribute_access.cc b/source/blender/blenkernel/intern/attribute_access.cc
index 5b2f588959f..c24630c5666 100644
--- a/source/blender/blenkernel/intern/attribute_access.cc
+++ b/source/blender/blenkernel/intern/attribute_access.cc
@@ -143,10 +143,8 @@ CustomDataType attribute_data_type_highest_complexity(Span<CustomDataType> data_
 static int attribute_domain_priority(const AttributeDomain domain)
 {
   switch (domain) {
-#if 0
     case ATTR_DOMAIN_CURVE:
       return 0;
-#endif
     case ATTR_DOMAIN_FACE:
       return 1;
     case ATTR_DOMAIN_EDGE:
diff --git a/source/blender/blenkernel/intern/curve_eval.cc b/source/blender/blenkernel/intern/curve_eval.cc
new file mode 100644
index 00000000000..4c63c7f05ee
--- /dev/null
+++ b/source/blender/blenkernel/intern/curve_eval.cc
@@ -0,0 +1,184 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "BLI_array.hh"
+#include "BLI_listbase.h"
+#include "BLI_span.hh"
+
+#include "DNA_curve_types.h"
+
+#include "BKE_curve.h"
+#include "BKE_spline.hh"
+
+using blender::float3;
+using blender::float4x4;
+using blender::Span;
+
+CurveEval *CurveEval::copy()
+{
+  CurveEval *new_curve = new CurveEval();
+
+  for (SplinePtr &spline : this->splines) {
+    new_curve->splines.append(spline->copy());
+  }
+
+  return new_curve;
+}
+
+void CurveEval::translate(const float3 &translation)
+{
+  for (SplinePtr &spline : this->splines) {
+    spline->translate(translation);
+    spline->mark_cache_invalid();
+  }
+}
+
+void CurveEval::transform(const float4x4 &matrix)
+{
+  for (SplinePtr &spline : this->splines) {
+    spline->transform(matrix);
+  }
+}
+
+void CurveEval::bounds_min_max(float3 &min, float3 &max, const bool use_evaluated) const
+{
+  for (const SplinePtr &spline : this->splines) {
+    spline->bounds_min_max(min, max, use_evaluated);
+  }
+}
+
+static BezierSpline::HandleType handle_type_from_dna_bezt(const eBezTriple_Handle dna_handle_type)
+{
+  switch (dna_handle_type) {
+    case HD_FREE:
+      return BezierSpline::HandleType::Free;
+    case HD_AUTO:
+      return BezierSpline::HandleType::Auto;
+    case HD_VECT:
+      return BezierSpline::HandleType::Vector;
+    case HD_ALIGN:
+      return BezierSpline::HandleType::Align;
+    case HD_AUTO_ANIM:
+      return BezierSpline::HandleType::Auto;
+    case HD_ALIGN_DOUBLESIDE:
+      return BezierSpline::HandleType::Align;
+  }
+  BLI_assert_unreachable();
+  return BezierSpline::HandleType::Auto;
+}
+
+static Spline::NormalCalculationMode normal_mode_from_dna_curve(const int twist_mode)
+{
+  switch (twist_mode) {
+    case CU_TWIST_Z_UP:
+      return Spline::NormalCalculationMode::ZUp;
+    case CU_TWIST_MINIMUM:
+      return Spline::NormalCalculationMode::Minimum;
+    case CU_TWIST_TANGENT:
+      return Spline::NormalCalculationMode::Tangent;
+  }
+  BLI_assert_unreachable();
+  return Spline::NormalCalculationMode::Minimum;
+}
+
+static NURBSpline::KnotsMode knots_mode_from_dna_nurb(const short flag)
+{
+  switch (flag & (CU_NURB_ENDPOINT | CU_NURB_BEZIER)) {
+    case CU_NURB_ENDPOINT:
+      return NURBSpline::KnotsMode::EndPoint;
+    case CU_NURB_BEZIER:
+      return NURBSpline::KnotsMode::Bezier;
+    default:
+      return NURBSpline::KnotsMode::Normal;
+  }
+
+  BLI_assert_unreachable();
+  return NURBSpline::KnotsMode::Normal;
+}
+
+std::unique_ptr<CurveEval> curve_eval_from_dna_curve(const Curve &dna_curve)
+{
+  std::unique_ptr<CurveEval> curve = std::make_unique<CurveEval>();
+
+  const ListBase *nurbs = BKE_curve_nurbs_get(&const_cast<Curve &>(dna_curve));
+
+  curve->splines.reserve(BLI_listbase_count(nurbs));
+
+  /* TODO: Optimize by reserving the correct points size. */
+  LISTBASE_FOREACH (const Nurb *, nurb, nurbs) {
+    switch (nurb->type) {
+      case CU_BEZIER: {
+        std::unique_ptr<BezierSpline> spline = std::make_unique<BezierSpline>();
+        spline->set_resolution(nurb->resolu);
+        spline->set_cyclic(nurb->flagu & CU_NURB_CYCLIC);
+
+        for (const BezTriple &bezt : Span(nurb->bezt, nurb->pntsu)) {
+          spline->add_point(bezt.vec[1],
+                            handle_type_from_dna_bezt((eBezTriple_Handle)bezt.h1),
+                            bezt.vec[0],
+                            handle_type_from_dna_bezt((eBezTriple_Handle)bezt.h2),
+                            bezt.vec[2],
+                            bezt.radius,
+                            bezt.tilt);
+        }
+
+        curve->splines.append(std::move(spline));
+        break;
+      }
+      case CU_NURBS: {
+        std::unique_ptr<NURBSpline> spline = std::make_unique<NURBSpline>();
+        spline->set_resolution(nurb->resolu);
+        spline->set_cyclic(nurb->flagu & CU_NURB_CYCLIC);
+        spline->set_order(nurb->orderu);
+        spline->knots_mode = knots_mode_from_dna_nurb(nurb->flagu);
+
+        for (const BPoint &bp : Span(nurb->bp, nurb->pntsu)) {
+          spline->add_point(bp.vec, bp.radius, bp.tilt, bp.vec[3]);
+        }
+
+        curve->splines.append(std::move(spline));
+        break;
+      }
+      case CU_POLY: {
+        std::unique_ptr<PolySpline> spline = std::make_unique<PolySpline>();
+        spline->set_cyclic(nurb->flagu & CU_NURB_CYCLIC);
+
+        for (const BPoint &bp : Span(nurb->bp, nurb->pntsu)) {
+          spline->add_point(bp.vec, bp.radius, bp.tilt);
+        }
+
+        curve->splines.append(std::move(spline));
+        break;
+      }
+      default: {
+        BLI_assert_unreachable();
+        break;
+      }
+    }
+  }
+
+  /* Note: Normal mode is stored separately in each spline to facilitate combining splines
+   * from multiple curve objects, where the value may be different. */
+  const Spline::NormalCalculationMode normal_mode = normal_mode_from_dna_curve(
+      dna_curve.twist_mode);
+  for (SplinePtr &spline : curve->splines) {
+    spline->normal_mode = normal_mode;
+  }
+
+  return curve;
+}
+
+/** \} */
diff --git a/source/blender/blenkernel/intern/geometry_component_curve.cc b/source/blender/blenkernel/intern/geometry_component_curve.cc
new file mode 100644
index 00000000000..0490d577b88
--- /dev/null
+++ b/source/blender/blenkernel/intern/geometry_component_curve.cc
@@ -0,0 +1,299 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "BKE_spline.hh"
+
+#include "BKE_attribute_access.hh"
+#include "BKE_attribute_math.hh"
+#include "BKE_geometry_set.hh"
+
+#include "attribute_access_intern.hh"
+
+/* -------------------------------------------------------------------- */
+/** \name Geometry Component Implementation
+ * \{ */
+
+CurveComponent::CurveComponent() : GeometryComponent(GEO_COMPONENT_TYPE_CURVE)
+{
+}
+
+CurveComponent::~CurveComponent()
+{
+  this->clear();
+}
+
+GeometryComponent *CurveComponent::copy() const
+{
+  CurveComponent *new_component = new CurveComponent();
+  if (curve_ != nullptr) {
+    new_component->curve_ = curve_->copy();
+    new_component->ownership_ = GeometryOwnershipType::Owned;
+  }
+  return new_component;
+}
+
+void CurveComponent::clear()
+{
+  BLI_assert(this->is_mutable());
+  if (curve_ != nullptr) {
+    if (ownership_ == GeometryOwnershipType::Owned) {
+      delete curve_;
+    }
+    curve_ = nullptr;
+  }
+}
+
+bool CurveComponent::has_curve() const
+{
+  return curve_ != nullptr;
+}
+
+/* Clear the component and replace it with the new curve. */
+void CurveComponent::replace(CurveEval *curve, GeometryOwnershipType ownership)
+{
+  BLI_assert(this->is_mutable());
+  this->clear();
+  curve_ = curve;
+  ownership_ = ownership;
+}
+
+CurveEval *CurveComponent::release()
+{
+  BLI_assert(this->is_mutable());
+  CurveEval *curve = curve_;
+  curve_ = nullptr;
+  return curve;
+}
+
+const CurveEval *CurveComponent::get_for_read() const
+{
+  return curve_;
+}
+
+CurveEval *CurveComponent::get_for_write()
+{
+  BLI_assert(this->is_mutable());
+  if (ownership_ == GeometryOwnershipType::ReadOnly) {
+    curve_ = curve_->copy();
+    ownership_ = GeometryOwnershipType::Owned;
+  }
+  return curve_;
+}
+
+bool CurveComponent::is_empty() const
+{
+  return curve_ == nullptr;
+}
+
+bool CurveComponent::owns_direct_data() const
+{
+  return ownership_ == GeometryOwnershipType::Owned;
+}
+
+void CurveComponent::ensure_owns_direct_data()
+{
+  BLI_assert(this->is_mutable());
+  if (ownership_ != GeometryOwnershipType::Owned) {
+    curve_ = curve_->copy();
+    ownership_ = GeometryOwnershipType::Owned;
+  }
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Attribute Access
+ * \{ */
+
+int CurveComponent::attribute_domain_size(const AttributeDomain domain) const
+{
+  if (curve_ == nullptr) {
+    return 0;
+  }
+  if (domain == ATTR_DOMAIN_POINT) {
+    int total = 0;
+    for (const SplinePtr &spline : curve_->splines) {
+      total += spline->size();
+    }
+    return total;
+  }
+  if (domain == ATTR_DOMAIN_CURVE) {
+    return curve_->splines.size();
+  }
+  return 0;
+}
+
+namespace blender::bke {
+
+class BuiltinSplineAttributeProvider final : public BuiltinAttributeProvider {
+  using AsReadAttribute = GVArrayPtr (*)(const CurveEval &data);
+  using AsWriteAttribute = GVMutableArrayPtr (*)(CurveEval &data);
+  using UpdateOnWrite = void (*)(Spline &spline);
+  const AsReadAttribute as_read_attribute_;
+  const AsWriteAttribute as_write_attribute_;
+
+ public:
+  BuiltinSplineAttributeProvider(std::string attribute_name,
+                                 const CustomDataType attribute_type,
+                                 const WritableEnum writable,
+                                 const AsReadAttribute as_read_attribute,
+                                 const AsWriteAttribute as_write_attribute)
+      : BuiltinAttributeProvider(std::move(attribute_name),
+                                 ATTR_DOMAIN_CURVE,
+                                 attribute_type,
+                                 BuiltinAttributeProvider::NonCreatable,
+                                 writable,
+                                 BuiltinAttributeProvider::NonDeletable),
+        as_read_attribute_(as_read_attribute),
+        as_write_attribute_(as_write_attribute)
+  {
+  }
+
+  GVArrayPtr try_get_for_read(const GeometryComponent &component) const final
+  {
+    const CurveComponent &curve_component = static_cast<const CurveComponent &>(component);
+    const CurveEval *curve = curve_component.get_for_read();
+    if (curve == nullptr) {
+      return {};
+    }
+
+    return as_read_attribute_(*curve);
+  }
+
+  GVMutableArrayPtr try_get_for_write(GeometryComponent &component) const final
+  {
+    if (writable_ != Writable) {
+      return {};
+    }
+    CurveComponent &curve_component = static_cast<CurveComponent &>(component);
+    CurveEval *curve = curve_component.get_for_write();
+    if (curve == nullptr) {
+      return {};
+    }
+
+    return as_write_attribute_(*curve);
+  }
+
+  bool try_delete(GeometryComponent &UNUSED(component)) const final
+  {
+    return false;
+  }
+
+  bool try_create(GeometryComponent &UNUSED(component),
+                  const AttributeInit &UNUSED(initializer)) const final
+  {
+    return false;
+  }
+
+  bool exists(const GeometryComponent &component) const final
+  {
+    return component.attribute_domain_size(ATTR_DOMAIN_CURVE) != 0;
+  }
+};
+
+static int get_spline_resolution(const SplinePtr &spline)
+{
+  if (const BezierSpline *bezier_spline = dynamic_cast<const BezierSpline *>(spline.get())) {
+    return bezier_spline->resolution();
+  }
+  if (const NURBSpline *nurb_spline = dynamic_cast<const NURBSpline *>(spline.get())) {
+    return nurb_spline->resolution();
+  }
+  return 1;
+}
+
+static void set_spline_resolution(SplinePtr &spline, const int resolution)
+{
+  if (BezierSpline *bezier_spline = dynamic_cast<BezierSpline *>(spline.get())) {
+    bezier_spline->set_resolution(std::max(resolution, 1));
+  }
+  if (NURBSpline *nurb_spline = dynamic_cast<NURBSpline *>(spline.get())) {
+    nurb_spline->set_resolution(std::max(resolution, 1));
+  }
+}
+
+static GVArrayPtr make_resolution_read_attribute(const CurveEval &curve)
+{
+  return std::make_unique<fn::GVArray_For_DerivedSpan<SplinePtr, int, get_spline_resolution>>(
+      curve.splines.as_span());
+}
+
+static GVMutableArrayPtr make_resolution_write_attribute(CurveEval &curve)
+{
+  return std::make_unique<fn::GVMutableArray_For_DerivedSpan<SplinePtr,
+                                                             int,
+                                                             get_spline_resolution,
+                                                             set_spline_resolution>>(
+      curve.splines.as_mutable_span());
+}
+
+static bool get_cyclic_value(const SplinePtr &spline)
+{
+  return spline->is_cyclic();
+}
+
+static void set_cyclic_value(SplinePtr &spline, const bool value)
+{
+  if (spline->is_cyclic() != value) {
+    spline->set_cyclic(value);
+    spline->mark_cache_invalid();
+  }
+}
+
+static GVArrayPtr make_cyclic_read_attribute(const CurveEval &curve)
+{
+  return std::make_unique<fn::GVArray_For_DerivedSpan<SplinePtr, bool, get_cyclic_value>>(
+      curve.splines.as_span());
+}
+
+static GVMutableArrayPtr make_cyclic_write_attribute(CurveEval &curve)
+{
+  return std::make_unique<
+      fn::GVMutableArray_For_DerivedSpan<SplinePtr, bool, get_cyclic_value, set_cyclic_value>>(
+      curve.splines.as_mutable_span());
+}
+
+/**
+ * In this function all the attribute providers for a curve component are created. Most data
+ * in this function is statically allocated, because it does not change over time.
+ */
+static ComponentAttributeProviders create_attribute_providers_for_curve()
+{
+  static BuiltinSplineAttributeProvider resolution("resolution",
+                                                   CD_PROP_INT32,
+                                                   BuiltinAttributeProvider::Writable,
+                                                   make_resolution_read_attribute,
+                                                   make_resolution_write_attribute);
+
+  static BuiltinSplineAttributeProvider cyclic("cyclic",
+                                               CD_PROP_BOOL,
+                                               BuiltinAttributeProvider::Writable,
+                                               make_cyclic_read_attribute,
+                                               make_cyclic_write_attribute);
+
+  return ComponentAttributeProviders({&resolution, &cyclic}, {});
+}
+
+}  // namespace blender::bke
+
+const blender::bke::ComponentAttributeProviders *CurveComponent::get_attribute_providers() const
+{
+  static blender::bke::ComponentAttributeProviders providers =
+      blender::bke::create_attribute_providers_for_curve();
+  return &providers;
+}
+
+/** \} */
diff --git a/source/blender/blenkernel/intern/geometry_set.cc b/source/blender/blenkernel/intern/geometry_set.cc
index 6ae78343fea..3d85118deee 100644
--- a/source/blender/blenkernel/intern/geometry_set.cc
+++ b/source/blender/blenkernel/intern/geometry_set.cc
@@ -24,6 +24,7 @@
 #include "BKE_mesh_wrapper.h"
 #include "BKE_modifier.h"
 #include "BKE_pointcloud.h"
+#include "BKE_spline.hh"
 #include "BKE_volume.h"
 
 #include "DNA_collection_types.h"
@@ -60,6 +61,8 @@ GeometryComponent *GeometryComponent::create(GeometryComponentType component_typ
       return new InstancesComponent();
     case GEO_COMPONENT_TYPE_VOLUME:
       return new VolumeComponent();
+    case GEO_COMPONENT_TYPE_CURVE:
+      return new CurveComponent();
   }
   BLI_assert_unreachable();
   return nullptr;
@@ -182,6 +185,11 @@ void GeometrySet::compute_boundbox_without_instances(float3 *r_min, float3 *r_ma
   if (volume != nullptr) {
     BKE_volume_min_max(volume, *r_min, *r_max);
   }
+  const CurveEval *curve = this->get_curve_for_read();
+  if (curve != nullptr) {
+    /* Using the evaluated positions is somewhat arbitrary, but it is probably expected. */
+    curve->bounds_min_max(*r_min, *r_max, true);
+  }
 }
 
 std::ostream &operator<<(std::ostream &stream, const GeometrySet &geometry_set)
@@ -252,6 +260,13 @@ const Volume *GeometrySet::get_volume_for_read() const
   return (component == nullptr) ? nullptr : component->get_for_read();
 }
 
+/* Returns a read-only curve or null. */
+const CurveEval *GeometrySet::get_curve_for_read() const
+{
+  const CurveComponent *component = this->get_component_for_read<CurveComponent>();
+  return (component == nullptr) ? nullptr : component->get_for_read();
+}
+
 /* Returns true when the geometry set has a point cloud component that has a point cloud. */
 bool GeometrySet::has_pointcloud() const
 {
@@ -273,6 +288,13 @@ bool GeometrySet::has_volume() const
   return component != nullptr && component->has_volume();
 }
 
+/* Returns true when the geometry set has a curve component that has a curve. */
+bool GeometrySet::has_curve() const
+{
+  const CurveComponent *component = this->get_component_for_read<CurveComponent>();
+  return component != nullptr && component->has_curve();
+}
+
 /* Create a new geometry set that only contains the given mesh. */
 GeometrySet GeometrySet::create_with_mesh(Mesh *mesh, GeometryOwnershipType ownership)
 {
@@ -292,6 +314,15 @@ GeometrySet GeometrySet::create_with_pointcloud(PointCloud *pointcloud,
   return geometry_set;
 }
 
+/* Create a new geometry set that only contains the given curve. */
+GeometrySet GeometrySet::create_with_curve(CurveEval *curve, GeometryOwnershipType ownership)
+{
+  GeometrySet geometry_set;
+  CurveComponent &component = geometry_set.get_component_for_write<CurveComponent>();
+  component.replace(curve, ownership);
+  return geometry_set;
+}
+
 /* Clear the existing mesh and replace it with the given one. */
 void GeometrySet::replace_mesh(Mesh *mesh, GeometryOwnershipType ownership)
 {
@@ -299,6 +330,13 @@ void GeometrySet::replace_mesh(Mesh *mesh, GeometryOwnershipType ownership)
   component.replace(mesh, ownership);
 }
 
+/* Clear the existing curve and replace it with the given one. */
+void GeometrySet::replace_curve(CurveEval *curve, GeometryOwnershipType ownership)
+{
+  CurveComponent &component = this->get_component_for_write<CurveComponent>();
+  component.replace(curve, ownership);
+}
+
 /* Clear the existing point cloud and replace with the given one. */
 void GeometrySet::replace_pointcloud(PointCloud *pointcloud, GeometryOwnershipType ownership)
 {
@@ -334,6 +372,13 @@ Volume *GeometrySet::get_volume_for_write()
   return component.get_for_write();
 }
 
+/* Returns a mutable curve or null. No ownership is transferred. */
+CurveEval *GeometrySet::get_curve_for_write()
+{
+  CurveComponent &component = this->get_component_for_write<CurveComponent>();
+  return component.get_for_write();
+}
+
 /** \} */
 
 /* -------------------------------------------------------------------- */
diff --git a/source/blender/blenkernel/intern/geometry_set_instances.cc b/source/blender/blenkernel/intern/geometry_set_instances.cc
index 47db5d1f901..7e994b32ded 100644
--- a/source/blender/blenkernel/intern/geometry_set_instances.cc
+++ b/source/blender/blenkernel/intern/geometry_set_instances.cc
@@ -19,6 +19,7 @@
 #include "BKE_mesh_wrapper.h"
 #include "BKE_modifier.h"
 #include "BKE_pointcloud.h"
+#include "BKE_spline.hh"
 
 #include "DNA_collection_types.h"
 #include "DNA_mesh_types.h"
@@ -50,6 +51,16 @@ static void add_final_mesh_as_geometry_component(const Object &object, GeometryS
   }
 }
 
+static void add_curve_data_as_geometry_component(const Object &object, GeometrySet &geometry_set)
+{
+  BLI_assert(object.type == OB_CURVE);
+  if (object.data != nullptr) {
+    std::unique_ptr<CurveEval> curve = curve_eval_from_dna_curve(*(Curve *)object.data);
+    CurveComponent &curve_component = geometry_set.get_component_for_write<CurveComponent>();
+    curve_component.replace(curve.release(), GeometryOwnershipType::Owned);
+  }
+}
+
 /**
  * \note This doesn't extract instances from the "dupli" system for non-geometry-nodes instances.
  */
@@ -73,6 +84,9 @@ static GeometrySet object_get_geometry_set_for_read(const Object &object)
   if (object.type == OB_MESH) {
     add_final_mesh_as_geometry_component(object, geometry_set);
   }
+  else if (object.type == OB_CURVE) {
+    add_curve_data_as_geometry_component(object, geometry_set);
+  }
 
   /* TODO: Cover the case of point-clouds without modifiers-- they may not be covered by the
    * #geometry_set_eval case above. */
@@ -492,6 +506,28 @@ static void join_attributes(Span<GeometryInstanceGroup> set_groups,
   }
 }
 
+static void join_curve_splines(Span<GeometryInstanceGroup> set_groups, CurveComponent &result)
+{
+  CurveEval *new_curve = new CurveEval();
+  for (const GeometryInstanceGroup &set_group : set_groups) {
+    const GeometrySet &set = set_group.geometry_set;
+    if (!set.has_curve()) {
+      continue;
+    }
+
+    const CurveEval &source_curve = *set.get_curve_for_read();
+    for (const SplinePtr &source_spline : source_curve.splines) {
+      for (const float4x4 &transform : set_group.transforms) {
+        SplinePtr new_spline = source_spline->copy();
+        new_spline->transform(transform);
+        new_curve->splines.append(std::move(new_spline));
+      }
+    }
+  }
+
+  result.replace(new_curve);
+}
+
 static void join_instance_groups_mesh(Span<GeometryInstanceGroup> set_groups,
                                       bool convert_points_to_vertices,
                                       GeometrySet &result)
@@ -558,6 +594,12 @@ static void join_instance_groups_volume(Span<GeometryInstanceGroup> set_groups,
   UNUSED_VARS(set_groups, dst_component);
 }
 
+static void join_instance_groups_curve(Span<GeometryInstanceGroup> set_groups, GeometrySet &result)
+{
+  CurveComponent &dst_component = result.get_component_for_write<CurveComponent>();
+  join_curve_splines(set_groups, dst_component);
+}
+
 GeometrySet geometry_set_realize_mesh_for_modifier(const GeometrySet &geometry_set)
 {
   if (!geometry_set.has_instances() && !geometry_set.has_pointcloud()) {
@@ -589,6 +631,7 @@ GeometrySet geometry_set_realize_instances(const GeometrySet &geometry_set)
   join_instance_groups_mesh(set_groups, false, new_geometry_set);
   join_instance_groups_pointcloud(set_groups, new_geometry_set);
   join_instance_groups_volume(set_groups, new_geometry_set);
+  join_instance_groups_curve(set_groups, new_geometry_set);
 
   return new_geometry_set;
 }
diff --git a/source/blender/blenkernel/intern/node.cc b/source/blender/blenkernel/intern/node.cc
index d940326cc2c..6634c42a274 100644
--- a/source/blender/blenkernel/intern/node.cc
+++ b/source/blender/blenkernel/intern/node.cc
@@ -4947,6 +4947,7 @@ static void registerGeometryNodes()
   register_node_type_geo_boolean();
   register_node_type_geo_bounding_box();
   register_node_type_geo_collection_info();
+  register_node_type_geo_curve_to_mesh();
   register_node_type_geo_edge_split();
   register_node_type_geo_is_viewport();
   register_node_type_geo_join_geometry();
diff --git a/source/blender/blenkernel/intern/spline_base.cc b/source/blender/blenkernel/intern/spline_base.cc
new file mode 100644
index 00000000000..acfd3eba912
--- /dev/null
+++ b/source/blender/blenkernel/intern/spline_base.cc
@@ -0,0 +1,252 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "BLI_array.hh"
+#include "BLI_span.hh"
+
+#include "BKE_spline.hh"
+
+using blender::float3;
+using blender::IndexRange;
+using blender::MutableSpan;
+using blender::Span;
+
+Spline::Type Spline::type() const
+{
+  return type_;
+}
+
+void Spline::translate(const blender::float3 &translation)
+{
+  for (float3 &position : this->positions()) {
+    position += translation;
+  }
+  this->mark_cache_invalid();
+}
+
+void Spline::transform(const blender::float4x4 &matrix)
+{
+  for (float3 &position : this->positions()) {
+    position = matrix * position;
+  }
+  this->mark_cache_invalid();
+}
+
+int Spline::evaluated_edges_size() const
+{
+  const int eval_size = this->evaluated_points_size();
+  if (eval_size == 1) {
+    return 0;
+  }
+
+  return this->is_cyclic_ ? eval_size : eval_size - 1;
+}
+
+float Spline::length() const
+{
+  return this->evaluated_lengths().last();
+}
+
+int Spline::segments_size() const
+{
+  const int points_len = this->size();
+
+  return is_cyclic_ ? points_len : points_len - 1;
+}
+
+bool Spline::is_cyclic() const
+{
+  return is_cyclic_;
+}
+
+void Spline::set_cyclic(const bool value)
+{
+  is_cyclic_ = value;
+}
+
+static void accumulate_lengths(Span<float3> positions,
+                               const bool is_cyclic,
+                               MutableSpan<float> lengths)
+{
+  float length = 0.0f;
+  for (const int i : IndexRange(positions.size() - 1)) {
+    length += float3::distance(positions[i], positions[i + 1]);
+    lengths[i] = length;
+  }
+  if (is_cyclic) {
+    lengths.last() = length + float3::distance(positions.last(), positions.first());
+  }
+}
+
+/**
+ * Return non-owning access to the cache of accumulated lengths along the spline. Each item is the
+ * length of the subsequent segment, i.e. 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.
+ * However, the results are used often, so it makes sense to cache it.
+ */
+Span<float> Spline::evaluated_lengths() const
+{
+  if (!length_cache_dirty_) {
+    return evaluated_lengths_cache_;
+  }
+
+  std::lock_guard lock{length_cache_mutex_};
+  if (!length_cache_dirty_) {
+    return evaluated_lengths_cache_;
+  }
+
+  const int total = evaluated_edges_size();
+  evaluated_lengths_cache_.resize(total);
+
+  Span<float3> positions = this->evaluated_positions();
+  accumulate_lengths(positions, is_cyclic_, evaluated_lengths_cache_);
+
+  length_cache_dirty_ = false;
+  return evaluated_lengths_cache_;
+}
+
+static float3 direction_bisect(const float3 &prev, const float3 &middle, const float3 &next)
+{
+  const float3 dir_prev = (middle - prev).normalized();
+  const float3 dir_next = (next - middle).normalized();
+
+  return (dir_prev + dir_next).normalized();
+}
+
+static void calculate_tangents(Span<float3> positions,
+                               const bool is_cyclic,
+                               MutableSpan<float3> tangents)
+{
+  if (positions.size() == 1) {
+    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() = (positions[1] - positions[0]).normalized();
+    tangents.last() = (positions.last() - positions[positions.size() - 2]).normalized();
+  }
+}
+
+/**
+ * Return non-owning access to the direction of the curve at each evaluated point.
+ */
+Span<float3> Spline::evaluated_tangents() const
+{
+  if (!tangent_cache_dirty_) {
+    return evaluated_tangents_cache_;
+  }
+
+  std::lock_guard lock{tangent_cache_mutex_};
+  if (!tangent_cache_dirty_) {
+    return evaluated_tangents_cache_;
+  }
+
+  const int eval_size = this->evaluated_points_size();
+  evaluated_tangents_cache_.resize(eval_size);
+
+  Span<float3> positions = this->evaluated_positions();
+
+  if (eval_size == 1) {
+    evaluated_tangents_cache_.first() = float3(1.0f, 0.0f, 0.0f);
+  }
+  else {
+    calculate_tangents(positions, is_cyclic_, evaluated_tangents_cache_);
+    this->correct_end_tangents();
+  }
+
+  tangent_cache_dirty_ = false;
+  return evaluated_tangents_cache_;
+}
+
+static void calculate_normals_z_up(Span<float3> tangents, MutableSpan<float3> normals)
+{
+  for (const int i : normals.index_range()) {
+    normals[i] = float3::cross(tangents[i], float3(0.0f, 0.0f, 1.0f)).normalized();
+  }
+}
+
+/**
+ * Return non-owning access to the direction vectors perpendicular to the tangents at every
+ * evaluated point. The method used to generate the normal vectors depends on Spline.normal_mode.
+ *
+ * TODO: Support changing the normal based on tilts.
+ */
+Span<float3> Spline::evaluated_normals() const
+{
+  if (!normal_cache_dirty_) {
+    return evaluated_normals_cache_;
+  }
+
+  std::lock_guard lock{normal_cache_mutex_};
+  if (!normal_cache_dirty_) {
+    return evaluated_normals_cache_;
+  }
+
+  const int eval_size = this->evaluated_points_size();
+  evaluated_normals_cache_.resize(eval_size);
+
+  Span<float3> tangents = evaluated_tangents();
+
+  /* Only Z up normals are supported at the moment. */
+  calculate_normals_z_up(tangents, evaluated_normals_cache_);
+
+  normal_cache_dirty_ = false;
+  return evaluated_normals_cache_;
+}
+
+Spline::LookupResult Spline::lookup_evaluated_factor(const float factor) const
+{
+  return this->lookup_evaluated_length(this->length() * factor);
+}
+
+/**
+ * \note This does not support extrapolation currently.
+ */
+Spline::LookupResult Spline::lookup_evaluated_length(const float length) const
+{
+  BLI_assert(length >= 0.0f && length <= this->length());
+
+  Span<float> lengths = this->evaluated_lengths();
+
+  const float *offset = std::lower_bound(lengths.begin(), lengths.end(), length);
+  const int index = offset - lengths.begin();
+  const int next_index = (index == this->size() - 1) ? 0 : index + 1;
+
+  const float previous_length = (index == 0) ? 0.0f : lengths[index - 1];
+  const float factor = (length - previous_length) / (lengths[index] - previous_length);
+
+  return LookupResult{index, next_index, factor};
+}
+
+void Spline::bounds_min_max(float3 &min, float3 &max, const bool use_evaluated) const
+{
+  Span<float3> positions = use_evaluated ? this->evaluated_positions() : this->positions();
+  for (const float3 &position : positions) {
+    minmax_v3v3_v3(min, max, position);
+  }
+}
diff --git a/source/blender/blenkernel/intern/spline_bezier.cc b/source/blender/blenkernel/intern/spline_bezier.cc
new file mode 100644
index 00000000000..4981f441190
--- /dev/null
+++ b/source/blender/blenkernel/intern/spline_bezier.cc
@@ -0,0 +1,478 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "BLI_array.hh"
+#include "BLI_span.hh"
+#include "BLI_task.hh"
+
+#include "BKE_spline.hh"
+
+using blender::Array;
+using blender::float3;
+using blender::IndexRange;
+using blender::MutableSpan;
+using blender::Span;
+
+SplinePtr BezierSpline::copy() const
+{
+  return std::make_unique<BezierSpline>(*this);
+}
+
+int BezierSpline::size() const
+{
+  const int size = positions_.size();
+  BLI_assert(size == handle_types_left_.size());
+  BLI_assert(size == handle_positions_left_.size());
+  BLI_assert(size == handle_types_right_.size());
+  BLI_assert(size == handle_positions_right_.size());
+  BLI_assert(size == radii_.size());
+  BLI_assert(size == tilts_.size());
+  return size;
+}
+
+int BezierSpline::resolution() const
+{
+  return resolution_;
+}
+
+void BezierSpline::set_resolution(const int value)
+{
+  BLI_assert(value > 0);
+  resolution_ = value;
+  this->mark_cache_invalid();
+}
+
+void BezierSpline::add_point(const float3 position,
+                             const HandleType handle_type_start,
+                             const float3 handle_position_start,
+                             const HandleType handle_type_end,
+                             const float3 handle_position_end,
+                             const float radius,
+                             const float tilt)
+{
+  handle_types_left_.append(handle_type_start);
+  handle_positions_left_.append(handle_position_start);
+  positions_.append(position);
+  handle_types_right_.append(handle_type_end);
+  handle_positions_right_.append(handle_position_end);
+  radii_.append(radius);
+  tilts_.append(tilt);
+  this->mark_cache_invalid();
+}
+
+MutableSpan<float3> BezierSpline::positions()
+{
+  return positions_;
+}
+Span<float3> BezierSpline::positions() const
+{
+  return positions_;
+}
+MutableSpan<float> BezierSpline::radii()
+{
+  return radii_;
+}
+Span<float> BezierSpline::radii() const
+{
+  return radii_;
+}
+MutableSpan<float> BezierSpline::tilts()
+{
+  return tilts_;
+}
+Span<float> BezierSpline::tilts() const
+{
+  return tilts_;
+}
+Span<BezierSpline::HandleType> BezierSpline::handle_types_left() const
+{
+  return handle_types_left_;
+}
+MutableSpan<BezierSpline::HandleType> BezierSpline::handle_types_left()
+{
+  return handle_types_left_;
+}
+Span<float3> BezierSpline::handle_positions_left() const
+{
+  return handle_positions_left_;
+}
+MutableSpan<float3> BezierSpline::handle_positions_left()
+{
+  return handle_positions_left_;
+}
+Span<BezierSpline::HandleType> BezierSpline::handle_types_right() const
+{
+  return handle_types_right_;
+}
+MutableSpan<BezierSpline::HandleType> BezierSpline::handle_types_right()
+{
+  return handle_types_right_;
+}
+Span<float3> BezierSpline::handle_positions_right() const
+{
+  return handle_positions_right_;
+}
+MutableSpan<float3> BezierSpline::handle_positions_right()
+{
+  return handle_positions_right_;
+}
+
+void BezierSpline::translate(const blender::float3 &translation)
+{
+  for (float3 &position : this->positions()) {
+    position += translation;
+  }
+  for (float3 &handle_position : this->handle_positions_left()) {
+    handle_position += translation;
+  }
+  for (float3 &handle_position : this->handle_positions_right()) {
+    handle_position += translation;
+  }
+  this->mark_cache_invalid();
+}
+
+void BezierSpline::transform(const blender::float4x4 &matrix)
+{
+  for (float3 &position : this->positions()) {
+    position = matrix * position;
+  }
+  for (float3 &handle_position : this->handle_positions_left()) {
+    handle_position = matrix * handle_position;
+  }
+  for (float3 &handle_position : this->handle_positions_right()) {
+    handle_position = matrix * handle_position;
+  }
+  this->mark_cache_invalid();
+}
+
+bool BezierSpline::point_is_sharp(const int index) const
+{
+  return ELEM(handle_types_left_[index], HandleType::Vector, HandleType::Free) ||
+         ELEM(handle_types_right_[index], HandleType::Vector, HandleType::Free);
+}
+
+bool BezierSpline::segment_is_vector(const int index) const
+{
+  if (index == this->size() - 1) {
+    BLI_assert(is_cyclic_);
+    return handle_types_right_.last() == HandleType::Vector &&
+           handle_types_left_.first() == HandleType::Vector;
+  }
+  return handle_types_right_[index] == HandleType::Vector &&
+         handle_types_left_[index + 1] == HandleType::Vector;
+}
+
+void BezierSpline::mark_cache_invalid()
+{
+  offset_cache_dirty_ = true;
+  position_cache_dirty_ = true;
+  mapping_cache_dirty_ = true;
+  tangent_cache_dirty_ = true;
+  normal_cache_dirty_ = true;
+  length_cache_dirty_ = true;
+}
+
+int BezierSpline::evaluated_points_size() const
+{
+  const int points_len = this->size();
+  BLI_assert(points_len > 0);
+
+  const int last_offset = this->control_point_offsets().last();
+  if (is_cyclic_ && points_len > 1) {
+    return last_offset + (this->segment_is_vector(points_len - 1) ? 0 : resolution_);
+  }
+
+  return last_offset + 1;
+}
+
+/**
+ * If the spline is not cyclic, the direction for the first and last points is just the
+ * direction formed by the corresponding handles and control points. In the unlikely situation
+ * that the handles define a zero direction, fallback to using the direction defined by the
+ * first and last evaluated segments already calculated in #Spline::evaluated_tangents().
+ */
+void BezierSpline::correct_end_tangents() const
+{
+  if (is_cyclic_) {
+    return;
+  }
+
+  MutableSpan<float3> tangents(evaluated_tangents_cache_);
+
+  if (handle_positions_left_.first() != positions_.first()) {
+    tangents.first() = (positions_.first() - handle_positions_left_.first()).normalized();
+  }
+  if (handle_positions_right_.last() != positions_.last()) {
+    tangents.last() = (handle_positions_right_.last() - positions_.last()).normalized();
+  }
+}
+
+static void bezier_forward_difference_3d(const float3 &point_0,
+                                         const float3 &point_1,
+                                         const float3 &point_2,
+                                         const float3 &point_3,
+                                         MutableSpan<float3> result)
+{
+  BLI_assert(result.size() > 0);
+  const float inv_len = 1.0f / static_cast<float>(result.size());
+  const float inv_len_squared = inv_len * inv_len;
+  const float inv_len_cubed = inv_len_squared * inv_len;
+
+  const float3 rt1 = 3.0f * (point_1 - point_0) * inv_len;
+  const float3 rt2 = 3.0f * (point_0 - 2.0f * point_1 + point_2) * inv_len_squared;
+  const float3 rt3 = (point_3 - point_0 + 3.0f * (point_1 - point_2)) * inv_len_cubed;
+
+  float3 q0 = point_0;
+  float3 q1 = rt1 + rt2 + rt3;
+  float3 q2 = 2.0f * rt2 + 6.0f * rt3;
+  float3 q3 = 6.0f * rt3;
+  for (const int i : result.index_range()) {
+    result[i] = q0;
+    q0 += q1;
+    q1 += q2;
+    q2 += q3;
+  }
+}
+
+void BezierSpline::evaluate_bezier_segment(const int index,
+                                           const int next_index,
+                                           MutableSpan<float3> positions) const
+{
+  if (this->segment_is_vector(index)) {
+    BLI_assert(positions.size() == 1);
+    positions.first() = positions_[index];
+  }
+  else {
+    bezier_forward_difference_3d(positions_[index],
+                                 handle_positions_right_[index],
+                                 handle_positions_left_[next_index],
+                                 positions_[next_index],
+                                 positions);
+  }
+}
+
+/**
+ * Returns access to a cache of offsets into the evaluated point array for each control point.
+ * This is important because while most control point edges generate the number of edges specified
+ * by the resolution, vector segments only generate one edge.
+ */
+Span<int> BezierSpline::control_point_offsets() const
+{
+  if (!offset_cache_dirty_) {
+    return offset_cache_;
+  }
+
+  std::lock_guard lock{offset_cache_mutex_};
+  if (!offset_cache_dirty_) {
+    return offset_cache_;
+  }
+
+  const int points_len = this->size();
+  offset_cache_.resize(points_len);
+
+  MutableSpan<int> offsets = offset_cache_;
+
+  int offset = 0;
+  for (const int i : IndexRange(points_len - 1)) {
+    offsets[i] = offset;
+    offset += this->segment_is_vector(i) ? 1 : resolution_;
+  }
+  offsets.last() = offset;
+
+  offset_cache_dirty_ = false;
+  return offsets;
+}
+
+/**
+ * Returns non-owning access to an array of values containing the information necessary to
+ * interpolate values from the original control points to evaluated points. The control point
+ * index is the integer part of each value, and the factor used for interpolating to the next
+ * control point is the remaining factional part.
+ */
+Span<float> BezierSpline::evaluated_mappings() const
+{
+  if (!mapping_cache_dirty_) {
+    return evaluated_mapping_cache_;
+  }
+
+  std::lock_guard lock{mapping_cache_mutex_};
+  if (!mapping_cache_dirty_) {
+    return evaluated_mapping_cache_;
+  }
+
+  const int size = this->size();
+  const int eval_size = this->evaluated_points_size();
+  evaluated_mapping_cache_.resize(eval_size);
+  MutableSpan<float> mappings = evaluated_mapping_cache_;
+
+  if (eval_size == 1) {
+    mappings.first() = 0.0f;
+    mapping_cache_dirty_ = false;
+    return mappings;
+  }
+
+  Span<int> offsets = this->control_point_offsets();
+  Span<float> lengths = this->evaluated_lengths();
+
+  /* Subtract one from the index into the lengths array to get the length
+   * at the start point rather than the length at the end of the edge. */
+
+  const float first_segment_len = lengths[offsets[1] - 1];
+  for (const int eval_index : IndexRange(0, offsets[1])) {
+    const float point_len = eval_index == 0 ? 0.0f : lengths[eval_index - 1];
+    const float length_factor = (first_segment_len == 0.0f) ? 0.0f : 1.0f / first_segment_len;
+
+    mappings[eval_index] = point_len * length_factor;
+  }
+
+  const int grain_size = std::max(512 / resolution_, 1);
+  blender::parallel_for(IndexRange(1, size - 2), grain_size, [&](IndexRange range) {
+    for (const int i : range) {
+      const float segment_start_len = lengths[offsets[i] - 1];
+      const float segment_end_len = lengths[offsets[i + 1] - 1];
+      const float segment_len = segment_end_len - segment_start_len;
+      const float length_factor = (segment_len == 0.0f) ? 0.0f : 1.0f / segment_len;
+
+      for (const int eval_index : IndexRange(offsets[i], offsets[i + 1] - offsets[i])) {
+        const float factor = (lengths[eval_index - 1] - segment_start_len) * length_factor;
+        mappings[eval_index] = i + factor;
+      }
+    }
+  });
+
+  if (is_cyclic_) {
+    const float segment_start_len = lengths[offsets.last() - 1];
+    const float segment_end_len = this->length();
+    const float segment_len = segment_end_len - segment_start_len;
+    const float length_factor = (segment_len == 0.0f) ? 0.0f : 1.0f / segment_len;
+
+    for (const int eval_index : IndexRange(offsets.last(), eval_size - offsets.last())) {
+      const float factor = (lengths[eval_index - 1] - segment_start_len) * length_factor;
+      mappings[eval_index] = size - 1 + factor;
+    }
+    mappings.last() = 0.0f;
+  }
+  else {
+    mappings.last() = size - 1;
+  }
+
+  mapping_cache_dirty_ = false;
+  return mappings;
+}
+
+Span<float3> BezierSpline::evaluated_positions() const
+{
+  if (!position_cache_dirty_) {
+    return evaluated_position_cache_;
+  }
+
+  std::lock_guard lock{position_cache_mutex_};
+  if (!position_cache_dirty_) {
+    return evaluated_position_cache_;
+  }
+
+  const int eval_size = this->evaluated_points_size();
+  evaluated_position_cache_.resize(eval_size);
+
+  MutableSpan<float3> positions = evaluated_position_cache_;
+
+  Span<int> offsets = this->control_point_offsets();
+  BLI_assert(offsets.last() <= eval_size);
+
+  const int grain_size = std::max(512 / resolution_, 1);
+  blender::parallel_for(IndexRange(this->size() - 1), grain_size, [&](IndexRange range) {
+    for (const int i : range) {
+      this->evaluate_bezier_segment(
+          i, i + 1, positions.slice(offsets[i], offsets[i + 1] - offsets[i]));
+    }
+  });
+
+  const int i_last = this->size() - 1;
+  if (is_cyclic_) {
+    this->evaluate_bezier_segment(i_last, 0, positions.slice(offsets.last(), resolution_));
+  }
+  else {
+    /* Since evualating the bezier segment doesn't add the final point,
+     * it must be added manually in the non-cyclic case. */
+    positions.last() = positions_.last();
+  }
+
+  position_cache_dirty_ = false;
+  return positions;
+}
+
+/**
+ * Convert the data encoded in #evaulated_mappings into its parts-- the information necessary
+ * to interpolate data from control points to evaluated points between them. The next control
+ * point index result will not overflow the size of the vector.
+ */
+BezierSpline::InterpolationData BezierSpline::interpolation_data_from_index_factor(
+    const float index_factor) const
+{
+  const int points_len = this->size();
+  const int index = std::floor(index_factor);
+  if (index == points_len) {
+    BLI_assert(is_cyclic_);
+    return InterpolationData{points_len - 1, 0, 1.0f};
+  }
+  if (index == points_len - 1) {
+    return InterpolationData{points_len - 2, points_len - 1, 1.0f};
+  }
+  return InterpolationData{index, index + 1, index_factor - index};
+}
+
+/* Use a spline argument to avoid adding this to the header. */
+template<typename T>
+static void interpolate_to_evaluated_points_impl(const BezierSpline &spline,
+                                                 const blender::VArray<T> &source_data,
+                                                 MutableSpan<T> result_data)
+{
+  Span<float> mappings = spline.evaluated_mappings();
+
+  for (const int i : result_data.index_range()) {
+    BezierSpline::InterpolationData interp = spline.interpolation_data_from_index_factor(
+        mappings[i]);
+
+    const T &value = source_data[interp.control_point_index];
+    const T &next_value = source_data[interp.next_control_point_index];
+
+    result_data[i] = blender::attribute_math::mix2(interp.factor, value, next_value);
+  }
+}
+
+blender::fn::GVArrayPtr BezierSpline::interpolate_to_evaluated_points(
+    const blender::fn::GVArray &source_data) const
+{
+  BLI_assert(source_data.size() == this->size());
+
+  const int eval_size = this->evaluated_points_size();
+  if (eval_size == 1) {
+    return source_data.shallow_copy();
+  }
+
+  blender::fn::GVArrayPtr new_varray;
+  blender::attribute_math::convert_to_static_type(source_data.type(), [&](auto dummy) {
+    using T = decltype(dummy);
+    if constexpr (!std::is_void_v<blender::attribute_math::DefaultMixer<T>>) {
+      Array<T> values(eval_size);
+      interpolate_to_evaluated_points_impl<T>(*this, source_data.typed<T>(), values);
+      new_varray = std::make_unique<blender::fn::GVArray_For_ArrayContainer<Array<T>>>(
+          std::move(values));
+    }
+  });
+
+  return new_varray;
+}
diff --git a/source/blender/blenkernel/intern/spline_nurbs.cc b/source/blender/blenkernel/intern/spline_nurbs.cc
new file mode 100644
index 00000000000..37d1232cfeb
--- /dev/null
+++ b/source/blender/blenkernel/intern/spline_nurbs.cc
@@ -0,0 +1,417 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "BLI_array.hh"
+#include "BLI_span.hh"
+#include "BLI_virtual_array.hh"
+
+#include "BKE_attribute_math.hh"
+#include "BKE_spline.hh"
+
+using blender::Array;
+using blender::float3;
+using blender::IndexRange;
+using blender::MutableSpan;
+using blender::Span;
+
+SplinePtr NURBSpline::copy() const
+{
+  return std::make_unique<NURBSpline>(*this);
+}
+
+int NURBSpline::size() const
+{
+  const int size = positions_.size();
+  BLI_assert(size == radii_.size());
+  BLI_assert(size == tilts_.size());
+  BLI_assert(size == weights_.size());
+  return size;
+}
+
+int NURBSpline::resolution() const
+{
+  return resolution_;
+}
+
+void NURBSpline::set_resolution(const int value)
+{
+  BLI_assert(value > 0);
+  resolution_ = value;
+  this->mark_cache_invalid();
+}
+
+uint8_t NURBSpline::order() const
+{
+  return order_;
+}
+
+void NURBSpline::set_order(const uint8_t value)
+{
+  BLI_assert(value >= 2 && value <= 6);
+  order_ = value;
+  this->mark_cache_invalid();
+}
+
+void NURBSpline::add_point(const float3 position,
+                           const float radius,
+                           const float tilt,
+                           const float weight)
+{
+  positions_.append(position);
+  radii_.append(radius);
+  tilts_.append(tilt);
+  weights_.append(weight);
+  knots_dirty_ = true;
+  this->mark_cache_invalid();
+}
+
+MutableSpan<float3> NURBSpline::positions()
+{
+  return positions_;
+}
+Span<float3> NURBSpline::positions() const
+{
+  return positions_;
+}
+MutableSpan<float> NURBSpline::radii()
+{
+  return radii_;
+}
+Span<float> NURBSpline::radii() const
+{
+  return radii_;
+}
+MutableSpan<float> NURBSpline::tilts()
+{
+  return tilts_;
+}
+Span<float> NURBSpline::tilts() const
+{
+  return tilts_;
+}
+MutableSpan<float> NURBSpline::weights()
+{
+  return weights_;
+}
+Span<float> NURBSpline::weights() const
+{
+  return weights_;
+}
+
+void NURBSpline::mark_cache_invalid()
+{
+  basis_cache_dirty_ = true;
+  position_cache_dirty_ = true;
+  tangent_cache_dirty_ = true;
+  normal_cache_dirty_ = true;
+  length_cache_dirty_ = true;
+}
+
+int NURBSpline::evaluated_points_size() const
+{
+  if (!this->check_valid_size_and_order()) {
+    return 0;
+  }
+  return resolution_ * this->segments_size();
+}
+
+void NURBSpline::correct_end_tangents() const
+{
+}
+
+bool NURBSpline::check_valid_size_and_order() const
+{
+  if (this->size() < order_) {
+    return false;
+  }
+
+  if (!is_cyclic_ && this->knots_mode == KnotsMode::Bezier) {
+    if (order_ == 4) {
+      if (this->size() < 5) {
+        return false;
+      }
+    }
+    else if (order_ != 3) {
+      return false;
+    }
+  }
+
+  return true;
+}
+
+int NURBSpline::knots_size() const
+{
+  const int size = this->size() + order_;
+  return is_cyclic_ ? size + order_ - 1 : size;
+}
+
+void NURBSpline::calculate_knots() const
+{
+  const KnotsMode mode = this->knots_mode;
+  const int length = this->size();
+  const int order = order_;
+
+  knots_.resize(this->knots_size());
+
+  MutableSpan<float> knots = knots_;
+
+  if (mode == NURBSpline::KnotsMode::Normal || is_cyclic_) {
+    for (const int i : knots.index_range()) {
+      knots[i] = static_cast<float>(i);
+    }
+  }
+  else if (mode == NURBSpline::KnotsMode::EndPoint) {
+    float k = 0.0f;
+    for (const int i : IndexRange(1, knots.size())) {
+      knots[i - 1] = k;
+      if (i >= order && i <= length) {
+        k += 1.0f;
+      }
+    }
+  }
+  else if (mode == NURBSpline::KnotsMode::Bezier) {
+    BLI_assert(ELEM(order, 3, 4));
+    if (order == 3) {
+      float k = 0.6f;
+      for (const int i : knots.index_range()) {
+        if (i >= order && i <= length) {
+          k += 0.5f;
+        }
+        knots[i] = std::floor(k);
+      }
+    }
+    else {
+      float k = 0.34f;
+      for (const int i : knots.index_range()) {
+        knots[i] = std::floor(k);
+        k += 1.0f / 3.0f;
+      }
+    }
+  }
+
+  if (is_cyclic_) {
+    const int b = length + order - 1;
+    if (order > 2) {
+      for (const int i : IndexRange(1, order - 2)) {
+        if (knots[b] != knots[b - i]) {
+          if (i == order - 1) {
+            knots[length + order - 2] += 1.0f;
+            break;
+          }
+        }
+      }
+    }
+
+    int c = order;
+    for (int i = b; i < this->knots_size(); i++) {
+      knots[i] = knots[i - 1] + (knots[c] - knots[c - 1]);
+      c--;
+    }
+  }
+}
+
+Span<float> NURBSpline::knots() const
+{
+  if (!knots_dirty_) {
+    BLI_assert(knots_.size() == this->size() + order_);
+    return knots_;
+  }
+
+  std::lock_guard lock{knots_mutex_};
+  if (!knots_dirty_) {
+    BLI_assert(knots_.size() == this->size() + order_);
+    return knots_;
+  }
+
+  this->calculate_knots();
+
+  knots_dirty_ = false;
+
+  return knots_;
+}
+
+static void calculate_basis_for_point(const float parameter,
+                                      const int points_len,
+                                      const int order,
+                                      Span<float> knots,
+                                      MutableSpan<float> basis_buffer,
+                                      NURBSpline::BasisCache &basis_cache)
+{
+  /* Clamp parameter due to floating point inaccuracy. TODO: Look into using doubles. */
+  const float t = std::clamp(parameter, knots[0], knots[points_len + order - 1]);
+
+  int start = 0;
+  int end = 0;
+  for (const int i : IndexRange(points_len + order - 1)) {
+    const bool knots_equal = knots[i] == knots[i + 1];
+    if (knots_equal || t < knots[i] || t > knots[i + 1]) {
+      basis_buffer[i] = 0.0f;
+      continue;
+    }
+
+    basis_buffer[i] = 1.0f;
+    start = std::max(i - order - 1, 0);
+    end = i;
+    basis_buffer.slice(i + 1, points_len + order - 1 - i).fill(0.0f);
+    break;
+  }
+  basis_buffer[points_len + order - 1] = 0.0f;
+
+  for (const int i_order : IndexRange(2, order - 1)) {
+    if (end + i_order >= points_len + order) {
+      end = points_len + order - 1 - i_order;
+    }
+    for (const int i : IndexRange(start, end - start + 1)) {
+      float new_basis = 0.0f;
+      if (basis_buffer[i] != 0.0f) {
+        new_basis += ((t - knots[i]) * basis_buffer[i]) / (knots[i + i_order - 1] - knots[i]);
+      }
+
+      if (basis_buffer[i + 1] != 0.0f) {
+        new_basis += ((knots[i + i_order] - t) * basis_buffer[i + 1]) /
+                     (knots[i + i_order] - knots[i + 1]);
+      }
+
+      basis_buffer[i] = new_basis;
+    }
+  }
+
+  /* Shrink the range of calculated values to avoid storing unnecessary zeros. */
+  while (basis_buffer[start] == 0.0f && start < end) {
+    start++;
+  }
+  while (basis_buffer[end] == 0.0f && end > start) {
+    end--;
+  }
+
+  basis_cache.weights.clear();
+  basis_cache.weights.extend(basis_buffer.slice(start, end - start + 1));
+  basis_cache.start_index = start;
+}
+
+void NURBSpline::calculate_basis_cache() const
+{
+  if (!basis_cache_dirty_) {
+    return;
+  }
+
+  std::lock_guard lock{basis_cache_mutex_};
+  if (!basis_cache_dirty_) {
+    return;
+  }
+
+  const int points_len = this->size();
+  const int eval_size = this->evaluated_points_size();
+  BLI_assert(this->evaluated_edges_size() > 0);
+  basis_cache_.resize(eval_size);
+
+  const int order = this->order();
+  Span<float> control_weights = this->weights();
+  Span<float> knots = this->knots();
+
+  MutableSpan<BasisCache> basis_cache(basis_cache_);
+
+  /* This buffer is reused by each basis calculation to store temporary values.
+   * Theoretically it could be optimized away in the future. */
+  Array<float> basis_buffer(this->knots_size());
+
+  const float start = knots[order - 1];
+  const float end = is_cyclic_ ? knots[points_len + order - 1] : knots[points_len];
+  const float step = (end - start) / this->evaluated_edges_size();
+  float parameter = start;
+  for (const int i : IndexRange(eval_size)) {
+    BasisCache &basis = basis_cache[i];
+    calculate_basis_for_point(
+        parameter, points_len + (is_cyclic_ ? order - 1 : 0), order, knots, basis_buffer, basis);
+    BLI_assert(basis.weights.size() <= order);
+
+    for (const int j : basis.weights.index_range()) {
+      const int point_index = (basis.start_index + j) % points_len;
+      basis.weights[j] *= control_weights[point_index];
+    }
+
+    parameter += step;
+  }
+
+  basis_cache_dirty_ = false;
+}
+
+template<typename T>
+void interpolate_to_evaluated_points_impl(Span<NURBSpline::BasisCache> weights,
+                                          const blender::VArray<T> &source_data,
+                                          MutableSpan<T> result_data)
+{
+  const int points_len = source_data.size();
+  BLI_assert(result_data.size() == weights.size());
+  blender::attribute_math::DefaultMixer<T> mixer(result_data);
+
+  for (const int i : result_data.index_range()) {
+    Span<float> point_weights = weights[i].weights;
+    const int start_index = weights[i].start_index;
+
+    for (const int j : point_weights.index_range()) {
+      const int point_index = (start_index + j) % points_len;
+      mixer.mix_in(i, source_data[point_index], point_weights[j]);
+    }
+  }
+
+  mixer.finalize();
+}
+
+blender::fn::GVArrayPtr NURBSpline::interpolate_to_evaluated_points(
+    const blender::fn::GVArray &source_data) const
+{
+  BLI_assert(source_data.size() == this->size());
+
+  this->calculate_basis_cache();
+  Span<BasisCache> weights(basis_cache_);
+
+  blender::fn::GVArrayPtr new_varray;
+  blender::attribute_math::convert_to_static_type(source_data.type(), [&](auto dummy) {
+    using T = decltype(dummy);
+    if constexpr (!std::is_void_v<blender::attribute_math::DefaultMixer<T>>) {
+      Array<T> values(this->evaluated_points_size());
+      interpolate_to_evaluated_points_impl<T>(weights, source_data.typed<T>(), values);
+      new_varray = std::make_unique<blender::fn::GVArray_For_ArrayContainer<Array<T>>>(
+          std::move(values));
+    }
+  });
+
+  return new_varray;
+}
+
+Span<float3> NURBSpline::evaluated_positions() const
+{
+  if (!position_cache_dirty_) {
+    return evaluated_position_cache_;
+  }
+
+  std::lock_guard lock{position_cache_mutex_};
+  if (!position_cache_dirty_) {
+    return evaluated_position_cache_;
+  }
+
+  const int eval_size = this->evaluated_points_size();
+  evaluated_position_cache_.resize(eval_size);
+
+  blender::fn::GVArray_Typed<float3> evaluated_positions{
+      this->interpolate_to_evaluated_points(blender::fn::GVArray_For_Span<float3>(positions_))};
+
+  evaluated_positions->materialize(evaluated_position_cache_);
+
+  position_cache_dirty_ = false;
+  return evaluated_position_cache_;
+}
diff --git a/source/blender/blenkernel/intern/spline_poly.cc b/source/blender/blenkernel/intern/spline_poly.cc
new file mode 100644
index 00000000000..09c578c0503
--- /dev/null
+++ b/source/blender/blenkernel/intern/spline_poly.cc
@@ -0,0 +1,105 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include "BLI_span.hh"
+#include "BLI_virtual_array.hh"
+
+#include "BKE_spline.hh"
+
+using blender::float3;
+using blender::MutableSpan;
+using blender::Span;
+
+SplinePtr PolySpline::copy() const
+{
+  return std::make_unique<PolySpline>(*this);
+}
+
+int PolySpline::size() const
+{
+  const int size = positions_.size();
+  BLI_assert(size == radii_.size());
+  BLI_assert(size == tilts_.size());
+  return size;
+}
+
+void PolySpline::add_point(const float3 position, const float radius, const float tilt)
+{
+  positions_.append(position);
+  radii_.append(radius);
+  tilts_.append(tilt);
+  this->mark_cache_invalid();
+}
+
+MutableSpan<float3> PolySpline::positions()
+{
+  return positions_;
+}
+Span<float3> PolySpline::positions() const
+{
+  return positions_;
+}
+MutableSpan<float> PolySpline::radii()
+{
+  return radii_;
+}
+Span<float> PolySpline::radii() const
+{
+  return radii_;
+}
+MutableSpan<float> PolySpline::tilts()
+{
+  return tilts_;
+}
+Span<float> PolySpline::tilts() const
+{
+  return tilts_;
+}
+
+void PolySpline::mark_cache_invalid()
+{
+  tangent_cache_dirty_ = true;
+  normal_cache_dirty_ = true;
+  length_cache_dirty_ = true;
+}
+
+int PolySpline::evaluated_points_size() const
+{
+  return this->size();
+}
+
+void PolySpline::correct_end_tangents() const
+{
+}
+
+Span<float3> PolySpline::evaluated_positions() const
+{
+  return this->positions();
+}
+
+/**
+ * Poly spline interpolation from control points to evaluated points is a special case, since
+ * the result data is the same as the input data. This function returns a GVArray that points to
+ * the original data. Therefore the lifetime of the returned virtual array must not be longer than
+ * the source data.
+ */
+blender::fn::GVArrayPtr PolySpline::interpolate_to_evaluated_points(
+    const blender::fn::GVArray &source_data) const
+{
+  BLI_assert(source_data.size() == this->size());
+
+  return source_data.shallow_copy();
+}