Repository: GraphiteEditor/Graphite (language: rust)

Issue:
Fix the spline node algorithm to be continuous across start/end points
Currently, the Spline from Points node doesn't handle continuity across the start/end point portion of the edge flow. This is what a square looks like when the node is applied, which fails in the top left corner where the start/end points connect:

![capture](https://github.com/user-attachments/assets/81d94ba5-9e56-4fcc-b907-75ea40ac012b)

Instead, this example should look like a symmetrical blob, with the same roundness around all corners.

This task involves updating the spline calculation algorithm, as pointed out by [this todo comment](https://github.com/GraphiteEditor/Graphite/commit/3f17e83e80d345a373f0df07ed642fc95de871e2#diff-e7a79f7e09ecb50a47584450e9b37dda9b437bbf0dbda93360c8ad1cdb8445c0R955).


Discussion / hints:
Hi @Keavon I was interested into contributing this code , I have already setup the project locally by following all the guidelines listed on the contribution page.

This will be the first time I'll be contributing to graphite, so I'm not totally aware about the codebase, still I wanted  to solve this issue I need to add a calculate_spline_control_points function that will calculate the spline control points instead of solve_spline_first_handle

and Update the spline node code to handle closed paths.

But I'm little confused to be honest any feedback on how should I approach this problem.

Looking forward for your response , 

Divv
Welcome @DivvSaxena! We can continue this discussion through Discord.
Just tried my first contribution, and have a working fix for this now on my computer. My method was to "overshoot" by using the two next nodes in the loop, and then copy the Bezier handle from the first "extra" node to the first node.
This was also one of my first "real" steps in Rust.
Thanks for working on it, @kalle95data! Could you please open a PR for it? (Let me know on Discord if you need help with the Git aspects.) I look forward to testing out your solution to see whether this produces a sufficiently close result to the mathematically "correct" approach. It'll certainly be a good improvement, even if we end up revisiting this later to implement the mathematically correct approach.
@0HyperCube wrote in Discord:

> Currently we implement equation 18 from https://mathworld.wolfram.com/CubicSpline.html#eqn19, if we use equation 19 with the periodic boundary condition (which requires an adapted version of the Thomas algorithm) we would solve it perfectly.

Fix patch:
diff --git a/libraries/bezier-rs/src/subpath/core.rs b/libraries/bezier-rs/src/subpath/core.rs
index de3575fe49..b48db37a37 100644
--- a/libraries/bezier-rs/src/subpath/core.rs
+++ b/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}"
+		);
+	}
+}
diff --git a/node-graph/gcore/src/vector/vector_nodes.rs b/node-graph/gcore/src/vector/vector_nodes.rs
index c82613eb08..a502ddf680 100644
--- a/node-graph/gcore/src/vector/vector_nodes.rs
+++ 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
 }


Test patch: