Refactor the 'Spline' node to use Kurbo instead of Bezier-rs (#2701)

* copy and refactor spline algorithm from Bezier-rs

* add comments

* refactor spline node

* fix spline when path is closed

* use refactored version of spline algorithm

* revert changes to 'Spline' node funtion

* use manipulator groups instead of bezpath

* Fix comment

---------

Co-authored-by: Keavon Chambers <keavon@keavon.com>

node-graph/gcore/src/vector/algorithms/mod.rs

@@ -3,3 +3,4 @@ mod instance;
 mod merge_by_distance;
 pub mod offset_subpath;
 mod poisson_disk;
+pub mod spline;

node-graph/gcore/src/vector/algorithms/spline.rs

@@ -0,0 +1,178 @@
+use glam::DVec2;
+
+/// Solve for the first handle of an open spline. (The opposite handle can be found by mirroring the result about the anchor.)
+pub fn solve_spline_first_handle_open(points: &[DVec2]) -> Vec<DVec2> {
+	let len_points = points.len();
+	if len_points == 0 {
+		return Vec::new();
+	}
+	if len_points == 1 {
+		return vec![points[0]];
+	}
+
+	// Matrix coefficients a, b and c (see https://mathworld.wolfram.com/CubicSpline.html).
+	// Because the `a` coefficients are all 1, they need not be stored.
+	// This algorithm does a variation of the above algorithm.
+	// Instead of using the traditional cubic (a + bt + ct^2 + dt^3), we use the bezier cubic.
+
+	let mut b = vec![DVec2::new(4., 4.); len_points];
+	b[0] = DVec2::new(2., 2.);
+	b[len_points - 1] = DVec2::new(2., 2.);
+
+	let mut c = vec![DVec2::new(1., 1.); len_points];
+
+	// 'd' is the the second point in a cubic bezier, which is what we solve for
+	let mut d = vec![DVec2::ZERO; len_points];
+
+	d[0] = DVec2::new(2. * points[1].x + points[0].x, 2. * points[1].y + points[0].y);
+	d[len_points - 1] = DVec2::new(3. * points[len_points - 1].x, 3. * points[len_points - 1].y);
+	for idx in 1..(len_points - 1) {
+		d[idx] = DVec2::new(4. * points[idx].x + 2. * points[idx + 1].x, 4. * points[idx].y + 2. * points[idx + 1].y);
+	}
+
+	// Solve with Thomas algorithm (see https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)
+	// Now we do row operations to eliminate `a` coefficients.
+	c[0] /= -b[0];
+	d[0] /= -b[0];
+	#[allow(clippy::assign_op_pattern)]
+	for i in 1..len_points {
+		b[i] += c[i - 1];
+		// For some reason this `+=` version makes the borrow checker mad:
+		// d[i] += d[i-1]
+		d[i] = d[i] + d[i - 1];
+		c[i] /= -b[i];
+		d[i] /= -b[i];
+	}
+
+	// At this point b[i] == -a[i + 1] and a[i] == 0.
+	// Now we do row operations to eliminate 'c' coefficients and solve.
+	d[len_points - 1] *= -1.;
+	#[allow(clippy::assign_op_pattern)]
+	for i in (0..len_points - 1).rev() {
+		d[i] = d[i] - (c[i] * d[i + 1]);
+		d[i] *= -1.; // d[i] /= b[i]
+	}
+
+	d
+}
+
+/// Solve for the first handle of a closed spline. (The opposite handle can be found by mirroring the result about the anchor.)
+/// If called with fewer than 3 points, this function will return an empty result.
+pub fn solve_spline_first_handle_closed(points: &[DVec2]) -> Vec<DVec2> {
+	let len_points = points.len();
+	if len_points < 3 {
+		return Vec::new();
+	}
+
+	// Matrix coefficients `a`, `b` and `c` (see https://mathworld.wolfram.com/CubicSpline.html).
+	// We don't really need to allocate them but it keeps the maths understandable.
+	let a = vec![DVec2::ONE; len_points];
+	let b = vec![DVec2::splat(4.); len_points];
+	let c = vec![DVec2::ONE; len_points];
+
+	let mut cmod = vec![DVec2::ZERO; len_points];
+	let mut u = vec![DVec2::ZERO; len_points];
+
+	// `x` is initially the output of the matrix multiplication, but is converted to the second value.
+	let mut x = vec![DVec2::ZERO; len_points];
+
+	for (i, point) in x.iter_mut().enumerate() {
+		let previous_i = i.checked_sub(1).unwrap_or(len_points - 1);
+		let next_i = (i + 1) % len_points;
+		*point = 3. * (points[next_i] - points[previous_i]);
+	}
+
+	// Solve using https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm#Variants (the variant using periodic boundary conditions).
+	// This code below is based on the reference C language implementation provided in that section of the article.
+	let alpha = a[0];
+	let beta = c[len_points - 1];
+
+	// Arbitrary, but chosen such that division by zero is avoided.
+	let gamma = -b[0];
+
+	cmod[0] = alpha / (b[0] - gamma);
+	u[0] = gamma / (b[0] - gamma);
+	x[0] /= b[0] - gamma;
+
+	// Handle from from `1` to `len_points - 2` (inclusive).
+	for ix in 1..=(len_points - 2) {
+		let m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
+		cmod[ix] = c[ix] * m;
+		u[ix] = (0.0 - a[ix] * u[ix - 1]) * m;
+		x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
+	}
+
+	// Handle `len_points - 1`.
+	let m = 1.0 / (b[len_points - 1] - alpha * beta / gamma - beta * cmod[len_points - 2]);
+	u[len_points - 1] = (alpha - a[len_points - 1] * u[len_points - 2]) * m;
+	x[len_points - 1] = (x[len_points - 1] - a[len_points - 1] * x[len_points - 2]) * m;
+
+	// Loop from `len_points - 2` to `0` (inclusive).
+	for ix in (0..=(len_points - 2)).rev() {
+		u[ix] = u[ix] - cmod[ix] * u[ix + 1];
+		x[ix] = x[ix] - cmod[ix] * x[ix + 1];
+	}
+
+	let fact = (x[0] + x[len_points - 1] * beta / gamma) / (1.0 + u[0] + u[len_points - 1] * beta / gamma);
+
+	for ix in 0..(len_points) {
+		x[ix] -= fact * u[ix];
+	}
+
+	let mut real = vec![DVec2::ZERO; len_points];
+	for i in 0..len_points {
+		let previous = i.checked_sub(1).unwrap_or(len_points - 1);
+		let next = (i + 1) % len_points;
+		real[i] = x[previous] * a[next] + x[i] * b[i] + x[next] * c[i];
+	}
+
+	// The matrix is now solved.
+
+	// Since we have computed the derivative, work back to find the start handle.
+	for i in 0..len_points {
+		x[i] = (x[i] / 3.) + points[i];
+	}
+
+	x
+}
+
+#[test]
+fn closed_spline() {
+	use crate::vector::misc::{dvec2_to_point, point_to_dvec2};
+	use kurbo::{BezPath, ParamCurve, ParamCurveDeriv};
+
+	// These points are just chosen arbitrary
+	let points = [DVec2::new(0., 0.), DVec2::new(0., 0.), DVec2::new(6., 5.), DVec2::new(7., 9.), DVec2::new(2., 3.)];
+
+	// List of first handle or second point in a cubic bezier curve.
+	let first_handles = solve_spline_first_handle_closed(&points);
+
+	// Construct the Subpath
+	let mut bezpath = BezPath::new();
+	bezpath.move_to(dvec2_to_point(points[0]));
+
+	for i in 0..first_handles.len() {
+		let next_i = i + 1;
+		let next_i = if next_i == first_handles.len() { 0 } else { next_i };
+
+		// First handle or second point of a cubic Bezier curve.
+		let p1 = dvec2_to_point(first_handles[i]);
+		// Second handle or third point of a cubic Bezier curve.
+		let p2 = dvec2_to_point(2. * points[next_i] - first_handles[next_i]);
+		// Endpoint or fourth point of a cubic Bezier curve.
+		let p3 = dvec2_to_point(points[next_i]);
+
+		bezpath.curve_to(p1, p2, p3);
+	}
+
+	// For each pair of bézier curves, ensure that the second derivative is continuous
+	for (bézier_a, bézier_b) in bezpath.segments().zip(bezpath.segments().skip(1).chain(bezpath.segments().take(1))) {
+		let derivative2_end_a = point_to_dvec2(bézier_a.to_cubic().deriv().eval(1.));
+		let derivative2_start_b = point_to_dvec2(bézier_b.to_cubic().deriv().eval(0.));
+
+		assert!(
+			derivative2_end_a.abs_diff_eq(derivative2_start_b, 1e-10),
+			"second derivative at the end of a {derivative2_end_a} is equal to the second derivative at the start of b {derivative2_start_b}"
+		);
+	}
+}

node-graph/gcore/src/vector/vector_data/attributes.rs

@@ -1,6 +1,6 @@
 use crate::vector::misc::dvec2_to_point;
 use crate::vector::vector_data::{HandleId, VectorData};
-use bezier_rs::BezierHandles;
+use bezier_rs::{BezierHandles, ManipulatorGroup};
 use core::iter::zip;
 use dyn_any::DynAny;
 use glam::{DAffine2, DVec2};
@@ -849,7 +849,7 @@ impl VectorData {
 			})
 	}
 
-	fn build_stroke_path_iter(&self) -> StrokePathIter {
+	pub fn build_stroke_path_iter(&self) -> StrokePathIter {
 		let mut points = vec![StrokePathIterPointMetadata::default(); self.point_domain.ids().len()];
 		for (segment_index, (&start, &end)) in self.segment_domain.start_point.iter().zip(&self.segment_domain.end_point).enumerate() {
 			points[start].set(StrokePathIterPointSegmentMetadata::new(segment_index, false));
@@ -869,6 +869,12 @@ impl VectorData {
 		self.build_stroke_path_iter().map(|(group, closed)| bezier_rs::Subpath::new(group, closed))
 	}
 
+	/// Construct and return an iterator of Vec of `(bezier_rs::ManipulatorGroup<PointId>], bool)` for stroke.
+	/// The boolean in the tuple indicates if the path is closed.
+	pub fn stroke_manipulator_groups(&self) -> impl Iterator<Item = (Vec<ManipulatorGroup<PointId>>, bool)> {
+		self.build_stroke_path_iter()
+	}
+
 	/// Construct a [`kurbo::BezPath`] curve for stroke.
 	pub fn stroke_bezpath_iter(&self) -> impl Iterator<Item = kurbo::BezPath> {
 		self.build_stroke_path_iter().map(|(group, closed)| {

node-graph/gcore/src/vector/vector_nodes.rs

@@ -1,5 +1,6 @@
 use super::algorithms::bezpath_algorithms::{self, position_on_bezpath, sample_points_on_bezpath, tangent_on_bezpath};
 use super::algorithms::offset_subpath::offset_subpath;
+use super::algorithms::spline::{solve_spline_first_handle_closed, solve_spline_first_handle_open};
 use super::misc::{CentroidType, point_to_dvec2};
 use super::style::{Fill, Gradient, GradientStops, Stroke};
 use super::{PointId, SegmentDomain, SegmentId, StrokeId, VectorData, VectorDataTable};
@@ -1436,15 +1437,15 @@ async fn spline(_: impl Ctx, vector_data: VectorDataTable) -> VectorDataTable {
 		}
 
 		let mut segment_domain = SegmentDomain::default();
-		for subpath in vector_data_instance.instance.stroke_bezier_paths() {
-			let positions = subpath.manipulator_groups().iter().map(|group| group.anchor).collect::<Vec<_>>();
-			let closed = subpath.closed() && positions.len() > 2;
+		for (manipulator_groups, closed) in vector_data_instance.instance.stroke_manipulator_groups() {
+			let positions = manipulator_groups.iter().map(|group| group.anchor).collect::<Vec<_>>();
+			let closed = closed && positions.len() > 2;
 
 			// Compute control point handles for Bezier spline.
 			let first_handles = if closed {
-				bezier_rs::solve_spline_first_handle_closed(&positions)
+				solve_spline_first_handle_closed(&positions)
 			} else {
-				bezier_rs::solve_spline_first_handle_open(&positions)
+				solve_spline_first_handle_open(&positions)
 			};
 
 			let stroke_id = StrokeId::ZERO;
@@ -1453,8 +1454,8 @@ async fn spline(_: impl Ctx, vector_data: VectorDataTable) -> VectorDataTable {
 			for i in 0..(positions.len() - if closed { 0 } else { 1 }) {
 				let next_index = (i + 1) % positions.len();
 
-				let start_index = vector_data_instance.instance.point_domain.resolve_id(subpath.manipulator_groups()[i].id).unwrap();
-				let end_index = vector_data_instance.instance.point_domain.resolve_id(subpath.manipulator_groups()[next_index].id).unwrap();
+				let start_index = vector_data_instance.instance.point_domain.resolve_id(manipulator_groups[i].id).unwrap();
+				let end_index = vector_data_instance.instance.point_domain.resolve_id(manipulator_groups[next_index].id).unwrap();
 
 				let handle_start = first_handles[i];
 				let handle_end = positions[next_index] * 2. - first_handles[next_index];