Skip to content

Commit

Permalink
Fix the spline node algorithm to be continuous across start/end points (
Browse files Browse the repository at this point in the history 
GraphiteEditor#2092)

* Simplify spline node implementation using stroke_bezier_paths

* Improve closed splines

* Code review

---------

Co-authored-by: Keavon Chambers <[email protected]>
  • Loading branch information
@0HyperCube @Keavon
0HyperCube and Keavon authored Nov 7, 2024
1 parent c3b526a commit 320d030
Show file tree
Hide file tree
Showing 2 changed files with 122 additions and 104 deletions.
110 changes: 107 additions & 3 deletions libraries/bezier-rs/src/subpath/core.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -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}"
);
}
}
116 changes: 15 additions & 101 deletions node-graph/gcore/src/vector/vector_nodes.rs
View file
Original file line number Diff line number Diff line change
@@ -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};
Expand All @@ -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
Expand Down Expand Up @@ -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
}
Expand Down

0 comments on commit 320d030

Please sign in to comment.