Fix the spline node algorithm to be continuous across start/end points (#2092)

* Simplify spline node implementation using stroke_bezier_paths

* Improve closed splines

* Code review

---------

Co-authored-by: Keavon Chambers <[email protected]>

Author: 0HyperCube

libraries/bezier-rs/src/subpath/core.rs

@@ -325,7 +325,7 @@ impl<PointId: crate::Identifier> Subpath<PointId> {
 		// Number of points = number of points to find handles for
 		let len_points = points.len();
 
-		let out_handles = solve_spline_first_handle(&points);
+		let out_handles = solve_spline_first_handle_open(&points);
 
 		let mut subpath = Subpath::new(Vec::new(), false);
 
@@ -366,14 +366,14 @@ impl<PointId: crate::Identifier> Subpath<PointId> {
 	}
 }
 
-pub fn solve_spline_first_handle(points: &[DVec2]) -> Vec<DVec2> {
+pub fn solve_spline_first_handle_open(points: &[DVec2]) -> Vec<DVec2> {
 	let len_points = points.len();
 	if len_points == 0 {
 		return Vec::new();
 	}
 
 	// 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.
+	// 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.
 
@@ -417,3 +417,107 @@ pub fn solve_spline_first_handle(points: &[DVec2]) -> Vec<DVec2> {
 
 	d
 }
+
+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::splat(1.); len_points];
+	let b = vec![DVec2::splat(4.); len_points];
+	let c = vec![DVec2::splat(1.); 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() {
+	// 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.)];
+
+	let out_handles = solve_spline_first_handle_closed(&points);
+
+	// Construct the Subpath
+	let mut manipulator_groups = Vec::new();
+	for i in 0..out_handles.len() {
+		manipulator_groups.push(ManipulatorGroup::<EmptyId>::new(points[i], Some(2. * points[i] - out_handles[i]), Some(out_handles[i])));
+	}
+	let subpath = Subpath::new(manipulator_groups, true);
+
+	// For each pair of bézier curves, ensure that the second derivative is continuous
+	for (bézier_a, bézier_b) in subpath.iter().zip(subpath.iter().skip(1).chain(subpath.iter().take(1))) {
+		let derivative2_end_a = bézier_a.derivative().unwrap().derivative().unwrap().evaluate(crate::TValue::Parametric(1.));
+		let derivative2_start_b = bézier_b.derivative().unwrap().derivative().unwrap().evaluate(crate::TValue::Parametric(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_nodes.rs

@@ -1,6 +1,6 @@
 use super::misc::CentroidType;
 use super::style::{Fill, Gradient, GradientStops, Stroke};
-use super::{PointId, SegmentId, StrokeId, VectorData};
+use super::{PointId, SegmentDomain, SegmentId, StrokeId, VectorData};
 use crate::registry::types::{Angle, Fraction, IntegerCount, Length, SeedValue};
 use crate::renderer::GraphicElementRendered;
 use crate::transform::{Footprint, Transform, TransformMut};
@@ -11,7 +11,6 @@ use crate::{Color, GraphicElement, GraphicGroup};
 use bezier_rs::{Cap, Join, Subpath, SubpathTValue, TValue};
 use glam::{DAffine2, DVec2};
 use rand::{Rng, SeedableRng};
-use std::collections::{BTreeMap, BTreeSet, VecDeque};
 
 /// Implemented for types that can be converted to an iterator of vector data.
 /// Used for the fill and stroke node so they can be used on VectorData or GraphicGroup
@@ -857,120 +856,35 @@ async fn splines_from_points<F: 'n + Send>(
 		return vector_data;
 	}
 
-	// Extract points and take ownership of the segment domain for processing.
-	let points = &vector_data.point_domain;
-	let segments = std::mem::take(&mut vector_data.segment_domain);
-
-	// Map segment IDs to their indices using BTreeMap for deterministic ordering.
-	let segment_id_to_index = segments.ids().iter().copied().enumerate().map(|(i, id)| (id, i)).collect::<BTreeMap<_, _>>();
-
-	// Iterate over all segments to generate splines.
-	let mut visited_segments = BTreeSet::new();
-	for (segment_index, &segment_id) in segments.ids().iter().enumerate() {
-		// Skip segments that have already been visited.
-		if visited_segments.contains(&segment_id) {
-			continue;
-		}
-
-		let mut current_subpath_segments = Vec::new();
-		let mut queue = VecDeque::new();
-		queue.push_back(segment_index);
-
-		// Traverse the connected segments to form a subpath.
-		while let Some(segment_index) = queue.pop_front() {
-			// Skip segments that have already been visited, otherwise add them to the visited set and the current subpath.
-			let seg_id = segments.ids()[segment_index];
-			if visited_segments.contains(&seg_id) {
-				continue;
-			}
-			visited_segments.insert(seg_id);
-			current_subpath_segments.push(segment_index);
-
-			// Get the start and end points of the segment.
-			let start_point_index = segments.start_point()[segment_index];
-			let end_point_index = segments.end_point()[segment_index];
-
-			// For both start and end points, find and enqueue connected segments.
-			for point_index in [start_point_index, end_point_index] {
-				let mut connected_seg_ids = segments.start_connected(point_index).chain(segments.end_connected(point_index)).collect::<Vec<_>>();
-				connected_seg_ids.sort_unstable(); // Ensure deterministic order
-				for connected_seg_id in connected_seg_ids {
-					let connected_seg_index = *segment_id_to_index.get(&connected_seg_id).unwrap_or(&usize::MAX);
-					if connected_seg_index != usize::MAX && !visited_segments.contains(&connected_seg_id) {
-						queue.push_back(connected_seg_index);
-					}
-				}
-			}
-		}
-
-		// Build a mapping from each point to its connected points using BTreeMap for deterministic ordering.
-		let mut point_connections: BTreeMap<usize, Vec<usize>> = BTreeMap::new();
-		for &seg_index in &current_subpath_segments {
-			let start = segments.start_point()[seg_index];
-			let end = segments.end_point()[seg_index];
-			point_connections.entry(start).or_default().push(end);
-			point_connections.entry(end).or_default().push(start);
-		}
-
-		// Sort connected points for deterministic traversal.
-		for neighbors in point_connections.values_mut() {
-			neighbors.sort_unstable();
-		}
-
-		// Identify endpoints.
-		let endpoints = point_connections
-			.iter()
-			.filter(|(_, neighbors)| neighbors.len() == 1)
-			.map(|(&point_index, _)| point_index)
-			.collect::<Vec<_>>();
-
-		let mut ordered_point_indices = Vec::new();
-
-		// Start with the first endpoint or the first point if there are no endpoints because it's a closed subpath.
-		let start_point_index = endpoints.first().copied().unwrap_or_else(|| *point_connections.keys().next().unwrap());
-
-		// Traverse points to order them into a path.
-		let mut visited_points = BTreeSet::new();
-		let mut current_point = start_point_index;
-		loop {
-			ordered_point_indices.push(current_point);
-			visited_points.insert(current_point);
-
-			let Some(neighbors) = point_connections.get(&current_point) else { break };
-			let next_point = neighbors.iter().find(|&pt| !visited_points.contains(pt));
-			let Some(&next_point) = next_point else { break };
-			current_point = next_point;
-		}
-
-		// If it's a closed subpath, close the spline loop by adding the start point at the end.
-		let closed = endpoints.is_empty();
-		if closed {
-			ordered_point_indices.push(start_point_index);
-		}
-
-		// Collect the positions of the ordered points.
-		let positions = ordered_point_indices.iter().map(|&index| points.positions()[index]).collect::<Vec<_>>();
+	let mut segment_domain = SegmentDomain::default();
+	for subpath in vector_data.stroke_bezier_paths() {
+		let positions = subpath.manipulator_groups().iter().map(|group| group.anchor).collect::<Vec<_>>();
+		let closed = subpath.closed();
 
 		// Compute control point handles for Bezier spline.
-		// TODO: Make this support wrapping around between start and end points for closed subpaths.
-		let first_handles = bezier_rs::solve_spline_first_handle(&positions);
+		let first_handles = if closed {
+			bezier_rs::solve_spline_first_handle_closed(&positions)
+		} else {
+			bezier_rs::solve_spline_first_handle_open(&positions)
+		};
 
 		let stroke_id = StrokeId::ZERO;
 
 		// Create segments with computed Bezier handles and add them to vector data.
-		for i in 0..(positions.len() - 1) {
+		for i in 0..(positions.len() - if closed { 0 } else { 1 }) {
 			let next_index = (i + 1) % positions.len();
 
-			let start_index = ordered_point_indices[i];
-			let end_index = ordered_point_indices[next_index];
+			let start_index = vector_data.point_domain.resolve_id(subpath.manipulator_groups()[i].id).unwrap();
+			let end_index = vector_data.point_domain.resolve_id(subpath.manipulator_groups()[next_index].id).unwrap();
 
 			let handle_start = first_handles[i];
 			let handle_end = positions[next_index] * 2. - first_handles[next_index];
 			let handles = bezier_rs::BezierHandles::Cubic { handle_start, handle_end };
 
-			vector_data.segment_domain.push(SegmentId::generate(), start_index, end_index, handles, stroke_id);
+			segment_domain.push(SegmentId::generate(), start_index, end_index, handles, stroke_id);
 		}
 	}
+	vector_data.segment_domain = segment_domain;
 
 	vector_data
 }