mirrors/Graphite
1
0
Fork 0
mirror of https://github.com/GraphiteEditor/Graphite.git synced 2025-07-08 00:05:00 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions 2

Freehand tool smoothing

Browse source
This commit is contained in:
James Lindsay 2024-07-05 21:42:40 +01:00 • committed by Keavon Chambers
parent 1652c713a6
commit cd49b500d0
5 changed files with 735 additions and 62 deletions
Download patch file Download diff file Expand all files Collapse all files

199
editor/src/messages/tool/tool_messages/freehand_tool.rs
View file

@ -11,7 +11,7 @@ use graph_craft::document::NodeId;
use graphene_core::uuid::generate_uuid;
use graphene_core::vector::VectorModificationType;
use graphene_core::Color;
use graphene_std::vector::{PointId, SegmentId};
use graphene_std::vector::{HandleId, PointId, SegmentId};
use glam::DVec2;
@ -175,13 +175,117 @@ impl ToolTransition for FreehandTool {
#[derive(Clone, Debug, Default)]
struct FreehandToolData {
extend_from_start: bool,
end_point: Option<(DVec2, PointId)>,
positions: Vec<DVec2>,
dragged: bool,
required_tangent: DVec2,
last_tangent: DVec2,
start: Option<(PointId, DVec2)>,
end: Option<(PointId, DVec2)>,
segment: Option<SegmentId>,
last_segment: Option<SegmentId>,
weight: f64,
layer: Option<LayerNodeIdentifier>,
}
impl FreehandToolData {
fn smooth(&mut self, responses: &mut VecDeque<Message>) {
const MAX_POSITIONS: usize = 16;
let Some(layer) = self.layer else { return };
let fit = if self.positions.len() < MAX_POSITIONS {
bezier_rs::Subpath::<PointId>::fit_cubic(&self.positions, 1, self.required_tangent, 0.)
} else {
None
};
let Some(bezier) = fit.and_then(|subpath| subpath.iter().next()) else {
// Start a new segment
if let Some(point) = self.end {
self.positions.clear();
self.positions.push(point.1);
self.start = Some(point);
self.end = None;
self.required_tangent = self.last_tangent;
self.last_segment = self.segment;
self.segment = None;
}
return;
};
let set_point = |point: &mut Option<(PointId, DVec2)>, position: DVec2, layer, responses: &mut VecDeque<Message>| {
if let Some((id, current)) = point {
let delta = position - *current;
*current += delta;
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::ApplyPointDelta { point: *id, delta },
});
*id
} else {
let id = PointId::generate();
*point = Some((id, position));
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::InsertPoint { id, position },
});
id
}
};
let start = set_point(&mut self.start, bezier.start, layer, responses);
let end = set_point(&mut self.end, bezier.end, layer, responses);
let points = [start, end];
let handles = [bezier.handle_start().map(|handle| handle - bezier.start), bezier.handle_end().map(|handle| handle - bezier.end)];
let segment = if let Some(segment) = self.segment {
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::SetHandles { segment, handles },
});
segment
} else {
let id = SegmentId::generate();
self.segment = Some(id);
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::InsertSegment { id, points, handles },
});
id
};
// Mirror handles if appropriate
if let Some(last_segment) = self.last_segment {
if last_segment != segment {
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::SetG1Continuous {
handles: [HandleId::end(last_segment), HandleId::primary(segment)],
enabled: true,
},
});
}
}
self.last_tangent = match bezier.handles {
bezier_rs::BezierHandles::Linear => bezier.end - bezier.start,
bezier_rs::BezierHandles::Quadratic { handle } => bezier.end - handle,
bezier_rs::BezierHandles::Cubic { handle_end, .. } => bezier.end - handle_end,
}
.normalize_or_zero();
}
fn push(&mut self, document: &DocumentMessageHandler, layer: LayerNodeIdentifier, viewport: DVec2) {
let transform = document.metadata().transform_to_viewport(layer);
let position = transform.inverse().transform_point2(viewport);
if self.positions.last() != Some(&position) && position.is_finite() {
self.positions.push(position);
}
}
}
impl Fsm for FreehandToolFsmState {
type ToolData = FreehandToolData;
type ToolOptions = FreehandOptions;
@ -207,48 +311,50 @@ impl Fsm for FreehandToolFsmState {
(FreehandToolFsmState::Ready, FreehandToolMessage::DragStart) => {
responses.add(DocumentMessage::StartTransaction);
tool_data.dragged = false;
tool_data.extend_from_start = false;
tool_data.weight = tool_options.line_weight;
tool_data.positions.clear();
tool_data.required_tangent = DVec2::ZERO;
tool_data.last_tangent = DVec2::ZERO;
tool_data.start = None;
tool_data.end = None;
tool_data.segment = None;
tool_data.last_segment = None;
tool_data.dragged = false;
// Extend an endpoint of the selected path
if let Some((layer, _, position)) = should_extend(document, input.mouse.position, crate::consts::SNAP_POINT_TOLERANCE) {
let layer = if let Some((layer, point, position)) = should_extend(document, input.mouse.position, crate::consts::SNAP_POINT_TOLERANCE) {
tool_data.positions.push(position);
tool_data.start = Some((point, position));
tool_data.layer = Some(layer);
layer
} else {
responses.add(DocumentMessage::DeselectAllLayers);
extend_path_with_next_segment(tool_data, position, responses);
let parent = document.new_layer_parent(true);
return FreehandToolFsmState::Drawing;
}
let nodes = {
let node_type = resolve_document_node_type("Path").expect("Path node does not exist");
let node = node_type.to_document_node_default_inputs([], Default::default());
responses.add(DocumentMessage::DeselectAllLayers);
HashMap::from([(NodeId(0), node)])
};
let parent = document.new_layer_parent(true);
let layer = graph_modification_utils::new_custom(NodeId(generate_uuid()), nodes, parent, responses);
tool_options.fill.apply_fill(layer, responses);
tool_options.stroke.apply_stroke(tool_data.weight, layer, responses);
let nodes = {
let node_type = resolve_document_node_type("Path").expect("Path node does not exist");
let node = node_type.to_document_node_default_inputs([], Default::default());
HashMap::from([(NodeId(0), node)])
tool_data.layer = Some(layer);
parent
};
let layer = graph_modification_utils::new_custom(NodeId(generate_uuid()), nodes, parent, responses);
tool_options.fill.apply_fill(layer, responses);
tool_options.stroke.apply_stroke(tool_data.weight, layer, responses);
tool_data.layer = Some(layer);
let transform = document.metadata().transform_to_viewport(layer);
let position = transform.inverse().transform_point2(input.mouse.position);
extend_path_with_next_segment(tool_data, position, responses);
tool_data.push(&document, layer, input.mouse.position);
tool_data.smooth(responses);
FreehandToolFsmState::Drawing
}
(FreehandToolFsmState::Drawing, FreehandToolMessage::PointerMove) => {
if let Some(layer) = tool_data.layer {
let transform = document.metadata().transform_to_viewport(layer);
let position = transform.inverse().transform_point2(input.mouse.position);
extend_path_with_next_segment(tool_data, position, responses);
tool_data.push(&document, layer, input.mouse.position);
tool_data.smooth(responses);
tool_data.dragged = true;
}
FreehandToolFsmState::Drawing
@ -294,34 +400,3 @@ impl Fsm for FreehandToolFsmState {
responses.add(FrontendMessage::UpdateMouseCursor { cursor: MouseCursorIcon::Default });
}
}
fn extend_path_with_next_segment(tool_data: &mut FreehandToolData, position: DVec2, responses: &mut VecDeque<Message>) {
if !tool_data.end_point.map_or(true, |(last_pos, _)| position != last_pos) || !position.is_finite() {
return;
}
let Some(layer) = tool_data.layer else { return };
let id = PointId::generate();
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::InsertPoint { id, position },
});
if let Some((_, previous_position)) = tool_data.end_point {
let next_id = SegmentId::generate();
let points = [previous_position, id];
responses.add(GraphOperationMessage::Vector {
layer,
modification_type: VectorModificationType::InsertSegment {
id: next_id,
points,
handles: [None, None],
},
});
}
tool_data.dragged = true;
tool_data.end_point = Some((position, id));
}

552
libraries/bezier-rs/src/bezier_fit.rs Normal file
View file

@ -0,0 +1,552 @@
/*
* Modifications and Rust port copyright (C) 2024 by 0Hypercube.
*
* Original version by lib2geom: <https://gitlab.com/inkscape/lib2geom>
*
* The entirety of this file is specially licensed under MPL 1.1 terms:
*
* Original code published in:
* An Algorithm for Automatically Fitting Digitized Curves
* by Philip J. Schneider
* "Graphics Gems", Academic Press, 1990
*
* Authors:
* Philip J. Schneider
* Lauris Kaplinski <lauris@kaplinski.com>
* Peter Moulder <pmoulder@mail.csse.monash.edu.au>
*
* Copyright (C) 1990 Philip J. Schneider
* Copyright (C) 2001 Lauris Kaplinski
* Copyright (C) 2001 Ximian, Inc.
* Copyright (C) 2003,2004 Monash University
* Original authors listed in the version control history of the following files:
* - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp
*
* This file is free software; you can redistribute it and/or modify it
* either under the terms of the Mozilla Public License Version 1.1 (the
* "MPL").
*
* The contents of this file are subject to the Mozilla Public License
* Version 1.1 (the "License"); you may not use this file except in
* compliance with the License. You may obtain a copy of the License at
* https://www.mozilla.org/MPL/1.1/
*
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
* OF ANY KIND, either express or implied. See the MPL for the specific
* language governing rights and limitations.
*/
use glam::DVec2;
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L192
fn bezier_fit_cubic_full(bezier: &mut [DVec2], data: &[DVec2], t_hat_1: DVec2, t_hat_2: DVec2, error: f64, max_beziers: usize) -> Option<usize> {
if data.len() < 2 {
return Some(0);
}
if data.len() == 2 {
// Fit 2 points
bezier[0] = data[0];
bezier[3] = data[data.len() - 1];
let dist = bezier[0].distance(bezier[3]) / 3.;
if dist.is_finite() {
bezier[1] = if t_hat_1 == DVec2::ZERO { (2. * bezier[0] + bezier[3]) / 3. } else { bezier[0] + dist * t_hat_1 };
bezier[2] = if t_hat_2 == DVec2::ZERO { (bezier[0] + 2. * bezier[3]) / 3. } else { bezier[3] + dist * t_hat_2 };
} else {
bezier[1] = bezier[0];
bezier[2] = bezier[3];
}
return Some(1);
}
// Parameterize points, and attempt to fit curve
let mut u = chord_length_parameterize(data);
if u[data.len() - 1] == 0. {
return Some(0);
}
generate_bezier(bezier, data, &u, t_hat_1, t_hat_2, error);
reparameterize(data, &mut u, bezier);
// Find max deviation of points to fitted curve.
let tolerance = (error + 1e-9).sqrt();
let (mut split_point, mut max_error_ratio) = compute_max_error_ratio(data, &u, bezier, tolerance, 0);
if max_error_ratio.abs() <= 1. {
return Some(1);
}
// If error not too large, then try some reparameterization and iteration.
if (0.0..=3.).contains(&max_error_ratio) {
let max_iterations = 4; // Max times to try iterating
for _ in 0..max_iterations {
generate_bezier(bezier, data, &u, t_hat_1, t_hat_2, error);
reparameterize(data, &mut u, bezier);
(split_point, max_error_ratio) = compute_max_error_ratio(data, &u, bezier, tolerance, split_point);
if (max_error_ratio.abs()) <= 1. {
return Some(1);
}
}
}
let is_corner = max_error_ratio < 0.;
if is_corner {
assert!(split_point < data.len());
if split_point == 0 {
if t_hat_1 == DVec2::ZERO {
// Got spike even with unconstrained initial tangent.
split_point += 1;
} else {
return bezier_fit_cubic_full(bezier, data, DVec2::ZERO, t_hat_2, error, max_beziers);
}
} else if split_point == data.len() - 1 {
if t_hat_2 == DVec2::ZERO {
// Got spike even with unconstrained final tangent.
split_point -= 1;
} else {
return bezier_fit_cubic_full(bezier, data, t_hat_1, DVec2::ZERO, error, max_beziers);
}
}
}
if 1 < max_beziers {
// Fitting failed -- split at max error point and fit recursively
let rec_max_beziers1 = max_beziers - 1;
let [rec_t_hat_1, rec_t_hat_2] = if is_corner {
if !(0 < split_point && split_point < data.len() - 1) {
return None;
}
[DVec2::ZERO; 2]
} else {
// Unit tangent vector at splitPoint.
let rec_t_hat_2 = darray_center_tangent(data, split_point, data.len());
[-rec_t_hat_2, rec_t_hat_2]
};
let nsegs1 = bezier_fit_cubic_full(bezier, &data[..split_point + 1], t_hat_1, rec_t_hat_2, error, rec_max_beziers1)?; //
assert!(nsegs1 != 0);
let rec_max_beziers2 = max_beziers - nsegs1;
let nsegs2 = bezier_fit_cubic_full(&mut bezier[nsegs1 * 4..], &data[split_point..], rec_t_hat_1, t_hat_2, error, rec_max_beziers2)?;
Some(nsegs1 + nsegs2)
} else {
None
}
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L538
fn reparameterize(d: &[DVec2], u: &mut [f64], bez_curve: &[DVec2]) {
let len = u.len();
assert!(2 <= len);
let last = len - 1;
assert!(bez_curve[0] == d[0]);
assert!(bez_curve[3] == d[last]);
assert!(u[0] == 0.);
assert!(u[last] == 1.);
// Otherwise, consider including 0 and last in the below loop.
for i in 1..last {
u[i] = newton_raphson_root_find(bez_curve, d[i], u[i]);
}
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L567
fn newton_raphson_root_find(q: &[DVec2], p: DVec2, u: f64) -> f64 {
assert!(0. <= u);
assert!(u <= 1.);
// Generate control vertices for Q'.
let mut q1 = [DVec2::ZERO; 3];
for i in 0..3 {
q1[i] = 3. * (q[i + 1] - q[i]);
}
// Generate control vertices for Q''.
let mut q2 = [DVec2::ZERO; 2];
for i in 0..2 {
q2[i] = 2. * (q1[i + 1] - q1[i]);
}
// Compute Q(u), Q'(u) and Q''(u).
let q_u = crate::Bezier::from_cubic_dvec2(q[0], q[1], q[2], q[3]).evaluate(crate::TValue::Parametric(u));
let q1_u = crate::Bezier::from_quadratic_dvec2(q1[0], q1[1], q1[2]).evaluate(crate::TValue::Parametric(u));
let q2_u = crate::Bezier::from_linear_dvec2(q2[0], q2[1]).evaluate(crate::TValue::Parametric(u));
// Compute f(u)/f'(u), where f is the derivative wrt u of distsq(u) = 0.5 * the square of the
// distance from P to Q(u). Here we're using Newton-Raphson to find a stationary point in the
// distsq(u), hopefully corresponding to a local minimum in distsq (and hence a local minimum
// distance from P to Q(u)).
let diff = q_u - p;
let numerator = diff.dot(q1_u);
let denominator = q1_u.dot(q1_u) + diff.dot(q2_u);
let mut improved_u = if denominator > 0. {
// One iteration of Newton-Raphson:
// improved_u = u - f(u)/f'(u)
u - (numerator / denominator)
} else {
// Using Newton-Raphson would move in the wrong direction (towards a local maximum rather
// than local minimum), so we move an arbitrary amount in the right direction.
if numerator > 0. {
u * 0.98 - 0.001
} else if numerator < 0. {
// Deliberately asymmetrical, to reduce the chance of cycling.
0.031 + u * 0.98
} else {
u
}
};
if !improved_u.is_finite() {
improved_u = u;
} else {
improved_u = improved_u.clamp(0., 1.);
}
// Ensure that improved_u isn't actually worse.
{
let diff_lensq = diff.length_squared();
let mut proportion = 0.125;
loop {
let bezier_pt = crate::Bezier::from_cubic_dvec2(q[0], q[1], q[2], q[3]).evaluate(crate::TValue::Parametric(improved_u));
if (bezier_pt - p).length_squared() > diff_lensq {
if proportion > 1. {
// g_warning("found proportion %g", proportion);
improved_u = u;
break;
}
improved_u = (1. - proportion) * improved_u + proportion * u;
} else {
break;
}
proportion += 0.125;
}
}
improved_u
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L807
fn darray_center_tangent(d: &[DVec2], center: usize, len: usize) -> DVec2 {
assert!(center != 0);
assert!(center < len - 1);
if d[center + 1] == d[center - 1] {
// Rotate 90 degrees in an arbitrary direction.
let diff = d[center] - d[center - 1];
diff.perp()
} else {
d[center - 1] - d[center + 1]
}
.normalize()
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L387
fn estimate_lengths(bezier: &mut [DVec2], data: &[DVec2], u_prime: &[f64], t_hat_1: DVec2, t_hat_2: DVec2) {
let len = data.len();
let mut c = [[0.; 2]; 2]; // Matrix C
let mut x = [0.; 2]; // Matrix X
// First and last control points of the Bezier curve are positioned exactly at the first and last data points.
bezier[0] = data[0];
bezier[3] = data[len - 1];
for i in 0..len {
// Bezier control point coefficients.
let b0 = (1. - u_prime[i]) * (1. - u_prime[i]) * (1. - u_prime[i]);
let b1 = 3. * u_prime[i] * (1. - u_prime[i]) * (1. - u_prime[i]);
let b2 = 3. * u_prime[i] * u_prime[i] * (1. - u_prime[i]);
let b3 = u_prime[i] * u_prime[i] * u_prime[i];
// rhs for eqn
let a1 = b1 * t_hat_1;
let a2 = b2 * t_hat_2;
c[0][0] += a1.dot(a1);
c[0][1] += a1.dot(a2);
c[1][0] = c[0][1];
c[1][1] += a2.dot(a2);
// Additional offset to the data point from the predicted point if we were to set bezier[1] to bezier[0] and bezier[2] to bezier[3].
let shortfall = data[i] - ((b0 + b1) * bezier[0]) - ((b2 + b3) * bezier[3]);
x[0] += a1.dot(shortfall);
x[1] += a2.dot(shortfall);
}
// We've constructed a pair of equations in the form of a matrix product C * alpha = X. Now solve for alpha.
// Compute the determinants of C and X.
let det_c0_c1 = c[0][0] * c[1][1] - c[1][0] * c[0][1];
let [mut alpha_l, mut alpha_r] = if det_c0_c1 != 0. {
// Apparently Kramer's rule.
let det_c0_x = c[0][0] * x[1] - c[0][1] * x[0];
let det_x_c1 = x[0] * c[1][1] - x[1] * c[0][1];
[det_x_c1 / det_c0_c1, det_c0_x / det_c0_c1]
} else {
// The matrix is under-determined. Try requiring alpha_l == alpha_r.
//
// One way of implementing the constraint alpha_l == alpha_r is to treat them as the same
// variable in the equations. We can do this by adding the columns of C to form a single
// column, to be multiplied by alpha to give the column vector X.
//
// We try each row in turn.
let c0 = c[0][0] + c[0][1];
if c0 != 0. {
[x[0] / c0; 2]
} else {
let c1 = c[1][0] + c[1][1];
if c1 != 0. {
[x[1] / c1; 2]
} else {
// Let the below code handle this.
[0.; 2]
}
}
};
// If alpha negative, use the Wu/Barsky heuristic (see text).
// (If alpha is 0, you get coincident control points that lead to divide by zero in any subsequent NewtonRaphsonRootFind() call.)
// TODO: Check whether this special-casing is necessary now that NewtonRaphsonRootFind handles non-positive denominator.
if alpha_l < 1.0e-6 || alpha_r < 1.0e-6 {
alpha_l = data[0].distance(data[len - 1]) / 3.;
alpha_r = alpha_l;
}
// Control points 1 and 2 are positioned an alpha distance out on the tangent vectors, left and right, respectively.
bezier[1] = alpha_l * t_hat_1 + bezier[0];
bezier[2] = alpha_r * t_hat_2 + bezier[3];
}
// ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent:
// Approximate unit tangents at endpoints and "center" of digitized curve
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L706
fn darray_left_tangent_simple(d: &[DVec2]) -> DVec2 {
let len = d.len();
assert!(len >= 2);
assert!(d[0] != d[1]);
(d[1] - d[0]).normalize()
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L724
fn darray_right_tangent_simple(d: &[DVec2]) -> DVec2 {
let len = d.len();
assert!(2 <= len);
let last = len - 1;
let prev = last - 1;
assert!(d[last] != d[prev]);
(d[prev] - d[last]).normalize()
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L745
fn darray_left_tangent(d: &[DVec2], tolerance_sq: f64) -> DVec2 {
let len = d.len();
assert!(2 <= len);
assert!(0. <= tolerance_sq);
let mut i = 1;
loop {
let pi = d[i];
let t = pi - d[0];
let distsq = t.length_squared();
if tolerance_sq < distsq {
return (t).normalize();
}
i += 1;
if i == len {
return if distsq == 0. { darray_left_tangent_simple(d) } else { (t).normalize() };
}
}
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L776
fn darray_right_tangent(d: &[DVec2], tolerance_sq: f64) -> DVec2 {
let len = d.len();
assert!(2 <= len);
assert!(0. <= tolerance_sq);
let last = len - 1;
let mut i = last - 1;
loop {
let pi = d[i];
let t = pi - d[last];
let dist_sq = t.length_squared();
if tolerance_sq < dist_sq {
return t.normalize_or_zero();
}
if i == 0 {
return if dist_sq == 0. { darray_right_tangent_simple(d) } else { (t).normalize_or_zero() };
}
i -= 1;
}
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L360
fn generate_bezier(bezier: &mut [DVec2], data: &[DVec2], u: &[f64], t_hat_1: DVec2, t_hat_2: DVec2, tolerance_sq: f64) {
let est1 = t_hat_1 == DVec2::ZERO;
let est2 = t_hat_2 == DVec2::ZERO;
let mut est_t_hat_1 = if est1 { darray_left_tangent(data, tolerance_sq) } else { t_hat_1 };
let est_t_hat_2 = if est2 { darray_right_tangent(data, tolerance_sq) } else { t_hat_2 };
estimate_lengths(bezier, data, u, est_t_hat_1, est_t_hat_2);
// We find that darray_right_tangent tends to produce better results for our current freehand tool than full estimation.
if est1 {
estimate_bi(bezier, 1, data, u);
if bezier[1] != bezier[0] {
est_t_hat_1 = (bezier[1] - bezier[0]).normalize_or_zero();
}
estimate_lengths(bezier, data, u, est_t_hat_1, est_t_hat_2);
}
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L492
fn estimate_bi(bezier: &mut [DVec2], ei: usize, data: &[DVec2], u: &[f64]) {
if !(1..=2).contains(&ei) {
return;
}
let oi = 3 - ei;
let mut num = [0., 0.];
let mut den = 0.;
for i in 0..data.len() {
let ui = u[i];
let b = [((1. - ui) * (1. - ui) * (1. - ui)), (3. * ui * (1. - ui) * (1. - ui)), (3. * ui * ui * (1. - ui)), (ui * ui * ui)];
for (d, num_d) in num.iter_mut().enumerate() {
*num_d += b[ei] * (b[0] * bezier[0][d] + b[oi] * bezier[oi][d] + b[3] * bezier[3][d] + -data[i][d]);
}
den -= b[ei] * b[ei];
}
if den != 0. {
for (d, num_d) in num.iter().enumerate() {
bezier[ei][d] = num_d / den;
}
} else {
bezier[ei] = (oi as f64 * bezier[0] + ei as f64 * bezier[3]) / 3.;
}
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L898
fn compute_max_error_ratio(d: &[DVec2], u: &[f64], bezier_curve: &[DVec2], tolerance: f64, mut split_point: usize) -> (usize, f64) {
let len = u.len();
assert!(2 <= len);
let last = len - 1;
assert!(bezier_curve[0] == d[0]);
assert!(bezier_curve[3] == d[last]);
assert!(u[0] == 0.);
assert!(u[last] == 1.);
// I.e. assert that the error for the first & last points is zero.
// Otherwise we should include those points in the below loop.
// The assertion is also necessary to ensure 0 < splitPoint < last.
let mut max_distance_sq = 0.; // Maximum error
let mut max_hook_ratio = 0.;
let mut snap_end = 0;
let mut previous_point = bezier_curve[0];
for i in 1..=last {
let current_point = crate::Bezier::from_cubic_dvec2(bezier_curve[0], bezier_curve[1], bezier_curve[2], bezier_curve[3]).evaluate(crate::TValue::Parametric(u[i]));
let distsq = (current_point - d[i]).length_squared();
if distsq > max_distance_sq {
max_distance_sq = distsq;
split_point = i;
}
let hook_ratio = compute_hook(previous_point, current_point, 0.5 * (u[i - 1] + u[i]), bezier_curve, tolerance);
if max_hook_ratio < hook_ratio {
max_hook_ratio = hook_ratio;
snap_end = i;
}
previous_point = current_point;
}
let dist_ratio = (max_distance_sq).sqrt() / tolerance;
let error_ratio = if max_hook_ratio <= dist_ratio {
dist_ratio
} else {
assert!(0 < snap_end);
split_point = snap_end - 1;
-max_hook_ratio
};
assert!(error_ratio == 0. || ((split_point < last) && (split_point != 0 || error_ratio < 0.)));
(split_point, error_ratio)
}
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L969
fn compute_hook(a: DVec2, b: DVec2, u: f64, bezier_curve: &[DVec2], tolerance: f64) -> f64 {
let point = crate::Bezier::from_cubic_dvec2(bezier_curve[0], bezier_curve[1], bezier_curve[2], bezier_curve[3]).evaluate(crate::TValue::Parametric(u));
let distance = ((a + b) / 2.).distance(point);
if distance < tolerance {
return 0.;
}
let allowed = a.distance(b) + tolerance;
distance / allowed
}
/// A value from 0..1 for each point in the path containing the total distance from the start along the path
// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier-utils.cpp#L833
fn chord_length_parameterize(d: &[DVec2]) -> Vec<f64> {
let len = d.len();
if len < 2 {
return Vec::new();
}
// First let u[i] equal the distance travelled along the path from d[0] to d[i]
let mut u = vec![0.; len];
for i in 1..len {
let dist = d[i].distance(d[i - 1]);
u[i] = u[i - 1] + dist;
}
// Scale to 0 ..= 1
let tot_len = u[len - 1];
if tot_len <= 0. {
return Vec::new();
}
if tot_len.is_finite() {
for u in u.iter_mut().skip(1) {
*u /= tot_len;
}
} else {
// fallback to even space
for (i, u) in u.iter_mut().enumerate().skip(1) {
*u = i as f64 / (len - 1) as f64;
}
}
// Ensure last is exactly 1.
u[len - 1] = 1.;
u
}
impl<PointId: crate::Identifier> crate::Subpath<PointId> {
pub fn fit_cubic(points: &[DVec2], max_segs: usize, tangent: DVec2, tolerance_sq: f64) -> Option<Self> {
let mut b = vec![DVec2::ZERO; 4 * max_segs];
let len = bezier_fit_cubic_full(&mut b, points, tangent, DVec2::ZERO, tolerance_sq, max_segs)?;
if len < 1 {
return None;
}
let beziers = (0..len)
.map(|i| crate::Bezier::from_cubic_dvec2(b[i * 4], b[i * 4 + 1], b[i * 4 + 2], b[i * 4 + 3]))
.collect::<Vec<_>>();
Some(Self::from_beziers(&beziers, false))
}
}
#[test]
fn generate_bezier_test() {
let src_b = vec![DVec2::new(5., -3.), DVec2::new(8., 0.), DVec2::new(4., 2.), DVec2::new(3., 3.)];
let t = [
0., 0.001, 0.03, 0.05, 0.09, 0.13, 0.18, 0.25, 0.29, 0.33, 0.39, 0.44, 0.51, 0.57, 0.62, 0.69, 0.75, 0.81, 0.91, 0.93, 0.97, 0.98, 0.999, 1.,
];
let data = t
.iter()
.map(|&t| crate::Bezier::from_cubic_dvec2(src_b[0], src_b[1], src_b[2], src_b[3]).evaluate(crate::TValue::Parametric(t)))
.collect::<Vec<_>>();
let t_hat_1 = (src_b[1] - src_b[0]).normalize();
let t_hat_2 = (src_b[2] - src_b[3]).normalize();
let mut est_b = vec![DVec2::ZERO; 4];
generate_bezier(&mut est_b, &data, &t, t_hat_1, t_hat_2, 1.);
assert_eq!(src_b, est_b);
}

2
libraries/bezier-rs/src/lib.rs
View file

@ -4,6 +4,7 @@
pub(crate) mod compare;
mod bezier;
mod bezier_fit;
mod consts;
mod poisson_disk;
mod polynomial;
@ -12,6 +13,7 @@ mod symmetrical_basis;
mod utils;
pub use bezier::*;
pub use bezier_fit::*;
pub use subpath::*;
pub use symmetrical_basis::*;
pub use utils::{Cap, Join, SubpathTValue, TValue, TValueType};

21
website/other/bezier-rs-demos/src/features/subpath-features.ts
View file

@ -7,6 +7,27 @@ const subpathFeatures = {
name: "Constructor",
callback: (subpath: WasmSubpathInstance): string => subpath.to_svg(),
},
fit: {
name: "Fit Points",
callback: (subpath: WasmSubpathInstance, options: Record<string, number>, _: undefined): string => subpath.fit(options.max_segments, options.error),
inputOptions: [
{
variable: "error",
min: 0,
max: 25,
step: 1,
default: 10,
},
{
variable: "max_segments",
min: 1,
max: 3,
step: 1,
default: 2,
},
],
},
insert: {
name: "Insert",
callback: (subpath: WasmSubpathInstance, options: Record<string, number>, _: undefined): string => subpath.insert(options.t, SUBPATH_T_VALUE_VARIANTS[options.TVariant]),

23
website/other/bezier-rs-demos/wasm/src/subpath.rs
View file

@ -76,6 +76,29 @@ impl WasmSubpath {
subpath_svg
}
pub fn fit(&self, max_segments: f64, tolerance: f64) -> String {
let mut result = String::new();
let subpath = self.0.clone();
let points = subpath.manipulator_groups().iter().map(|group| group.anchor).collect::<Vec<_>>();
let tolerance_sq = 0.02 * tolerance * tolerance * (0.2 * tolerance - 2.).exp();
if let Some(subpath) = Subpath::<EmptyId>::fit_cubic(&points, (max_segments as usize).max(1), DVec2::ZERO, tolerance_sq) {
subpath.to_svg(
&mut result,
CURVE_ATTRIBUTES.to_string().replace(BLACK, RED),
ANCHOR_ATTRIBUTES.to_string().replace(BLACK, RED),
HANDLE_ATTRIBUTES.to_string().replace(GRAY, RED),
HANDLE_LINE_ATTRIBUTES.to_string().replace(GRAY, RED),
);
}
for point in points {
result += &draw_circle(point, 4., BLACK, 1.5, WHITE);
}
wrap_svg_tag(result)
}
fn to_filled_svg(&self) -> String {
let mut subpath_svg = String::new();
self.0.to_svg(