From de6d4a7c4bb4f371d7698403e62c7206aa898708 Mon Sep 17 00:00:00 2001 From: Joe Neeman Date: Wed, 6 Aug 2025 22:46:12 -0500 Subject: [PATCH 1/4] Convert tangents_to_point --- Cargo.lock | 10 ++++++++++ Cargo.toml | 1 + libraries/bezier-rs/Cargo.toml | 4 ++++ libraries/bezier-rs/src/bezier/solvers.rs | 19 +++++++++++++++---- 4 files changed, 30 insertions(+), 4 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index b9a131eba5..34334b397e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -468,6 +468,7 @@ dependencies = [ "dyn-any", "glam", "kurbo", + "poly-cool", "serde", ] @@ -4040,6 +4041,15 @@ version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2f3a9f18d041e6d0e102a0a46750538147e5e8992d3b4873aaafee2520b00ce3" +[[package]] +name = "poly-cool" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6113ca52ade3ae52044cf0519e2c78ffa82120f6fa82f5099c8a4fd3ec8de43" +dependencies = [ + "arrayvec", +] + [[package]] name = "portable-atomic" version = "1.11.1" diff --git a/Cargo.toml b/Cargo.toml index f6dc8df0fc..99ead82fb6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -165,6 +165,7 @@ tracing-subscriber = { version = "0.3.19", features = ["env-filter"] } tracing = "0.1.41" rfd = "0.15.4" open = "5.3.2" +poly-cool = "0.2.0" [profile.dev] opt-level = 1 diff --git a/libraries/bezier-rs/Cargo.toml b/libraries/bezier-rs/Cargo.toml index ec9d6499f7..e030c5d8a7 100644 --- a/libraries/bezier-rs/Cargo.toml +++ b/libraries/bezier-rs/Cargo.toml @@ -13,9 +13,13 @@ homepage = "https://github.com/GraphiteEditor/Graphite/tree/master/libraries/bez repository = "https://github.com/GraphiteEditor/Graphite/tree/master/libraries/bezier-rs" documentation = "https://graphite.rs/libraries/bezier-rs/" +[features] +std = ["glam/std"] + [dependencies] # Required dependencies glam = { workspace = true } +poly-cool = { workspace = true } # Optional local dependencies dyn-any = { version = "0.3.0", path = "../dyn-any", optional = true } diff --git a/libraries/bezier-rs/src/bezier/solvers.rs b/libraries/bezier-rs/src/bezier/solvers.rs index 320e772d9f..b764bf6bd6 100644 --- a/libraries/bezier-rs/src/bezier/solvers.rs +++ b/libraries/bezier-rs/src/bezier/solvers.rs @@ -105,10 +105,21 @@ impl Bezier { /// #[must_use] pub fn tangents_to_point(self, point: DVec2) -> Vec { - let sbasis: crate::SymmetricalBasisPair = to_symmetrical_basis_pair(self); - let derivative = sbasis.derivative(); - let cross = (sbasis - point).cross(&derivative); - SymmetricalBasis::roots(&cross) + // We solve deriv(t) \times (self(t) - point) = 0. In principle, this is a quintic. + // In fact, the highest-order term cancels out so it's at most a quartic. + let (mut x, mut y) = self.parametric_polynomial(); + let x = x.coefficients_mut(); + let y = y.coefficients_mut(); + x[0] -= point.x; + y[0] -= point.y; + let poly = poly_cool::Poly::new([ + x[0] * y[1] - y[0] * x[1], + 2.0 * (x[0] * y[2] - y[0] * x[2]), + x[2] * y[1] - y[2] * x[1] + 2.0 * (x[1] * y[2] - y[1] * x[2]) + 3.0 * (x[0] * y[3] - y[0] * x[3]), + x[3] * y[1] - y[3] * x[1] + 3.0 * (x[1] * y[3] - y[1] * x[3]), + 2.0 * (x[3] * y[2] - y[3] * x[2]) + 3.0 * (x[2] * y[3] - y[2] * x[3]), + ]); + poly.roots_between(0.0, 1.0, 1e-8) } /// Returns a normalized unit vector representing the direction of the normal at the point `t` along the curve. From e4ca4f7280467938971c465a39a1ac87d34a4757 Mon Sep 17 00:00:00 2001 From: Joe Neeman Date: Wed, 6 Aug 2025 22:55:31 -0500 Subject: [PATCH 2/4] Convert normals_to_point --- libraries/bezier-rs/src/bezier/solvers.rs | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/libraries/bezier-rs/src/bezier/solvers.rs b/libraries/bezier-rs/src/bezier/solvers.rs index b764bf6bd6..11f36fdb61 100644 --- a/libraries/bezier-rs/src/bezier/solvers.rs +++ b/libraries/bezier-rs/src/bezier/solvers.rs @@ -132,10 +132,21 @@ impl Bezier { /// #[must_use] pub fn normals_to_point(self, point: DVec2) -> Vec { - let sbasis = to_symmetrical_basis_pair(self); - let derivative = sbasis.derivative(); - let cross = (sbasis - point).dot(&derivative); - SymmetricalBasis::roots(&cross) + // We solve deriv(t) \cdot (self(t) - point) = 0. + let (mut x, mut y) = self.parametric_polynomial(); + let x = x.coefficients_mut(); + let y = y.coefficients_mut(); + x[0] -= point.x; + y[0] -= point.y; + let poly = poly_cool::Poly::new([ + x[0] * x[1] + y[0] * y[1], + x[1] * x[1] + y[1] * y[1] + 2.0 * (x[0] * x[2] + y[0] * y[2]), + 3.0 * (x[2] * x[1] + y[2] * y[1]) + 3.0 * (x[0] * x[3] + y[0] * y[3]), + 4.0 * (x[3] * x[1] + y[3] * y[1]) + 2.0 * (x[2] * x[2] + y[2] * y[2]), + 5.0 * (x[3] * x[2] + y[3] * y[2]), + 3.0 * (x[3] * x[3] + y[3] * y[3]), + ]); + poly.roots_between(0.0, 1.0, 1e-8) } /// Returns the curvature, a scalar value for the derivative at the point `t` along the curve. From f080c2a6161ae6699cadc2693ccee698439f9f09 Mon Sep 17 00:00:00 2001 From: Joe Neeman Date: Thu, 7 Aug 2025 22:22:40 -0500 Subject: [PATCH 3/4] Port over roots and project --- libraries/bezier-rs/src/bezier/lookup.rs | 15 ++++++++------- libraries/bezier-rs/src/bezier/solvers.rs | 7 ++++--- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/libraries/bezier-rs/src/bezier/lookup.rs b/libraries/bezier-rs/src/bezier/lookup.rs index 499c933f03..1ddbe42fc5 100644 --- a/libraries/bezier-rs/src/bezier/lookup.rs +++ b/libraries/bezier-rs/src/bezier/lookup.rs @@ -250,23 +250,23 @@ impl Bezier { /// Returns the parametric `t`-value that corresponds to the closest point on the curve to the provided point. /// pub fn project(&self, point: DVec2) -> f64 { - let sbasis = crate::symmetrical_basis::to_symmetrical_basis_pair(*self); - let derivative = sbasis.derivative(); - let dd = (sbasis - point).dot(&derivative); - let roots = dd.roots(); + // The points at which the line from us to `point` is perpendicular + // to our curve are the critical points of the distance function. + let critical = self.normals_to_point(point); let mut closest = 0.; let mut min_dist_squared = self.evaluate(TValue::Parametric(0.)).distance_squared(point); - for time in roots { + for time in critical { let distance = self.evaluate(TValue::Parametric(time)).distance_squared(point); + dbg!(distance, time); if distance < min_dist_squared { closest = time; min_dist_squared = distance; } } - if self.evaluate(TValue::Parametric(1.)).distance_squared(point) < min_dist_squared { + if dbg!(self.evaluate(TValue::Parametric(1.)).distance_squared(point)) < min_dist_squared { closest = 1.; } closest @@ -354,7 +354,8 @@ mod tests { assert_eq!(bezier1.project(DVec2::new(100., 100.)), 1.); let bezier2 = Bezier::from_quadratic_coordinates(0., 0., 0., 100., 100., 100.); - assert_eq!(bezier2.project(DVec2::new(100., 0.)), 0.); + assert_eq!(bezier2.project(DVec2::new(99.99, 0.)), 0.); + assert!((bezier2.project(DVec2::new(-50., 150.)) - 0.5).abs() <= 1e-8); let bezier3 = Bezier::from_cubic_coordinates(-50., -50., -50., -50., 50., -50., 50., -50.); assert_eq!(DVec2::new(0., -50.), bezier3.evaluate(TValue::Parametric(bezier3.project(DVec2::new(0., -50.))))); diff --git a/libraries/bezier-rs/src/bezier/solvers.rs b/libraries/bezier-rs/src/bezier/solvers.rs index 11f36fdb61..a45191ea3c 100644 --- a/libraries/bezier-rs/src/bezier/solvers.rs +++ b/libraries/bezier-rs/src/bezier/solvers.rs @@ -1,7 +1,6 @@ use super::*; use crate::polynomial::Polynomial; use crate::utils::{TValue, solve_cubic, solve_quadratic}; -use crate::{SymmetricalBasis, to_symmetrical_basis_pair}; use glam::DMat2; use std::ops::Range; @@ -10,8 +9,10 @@ impl Bezier { /// Get roots as [[x], [y]] #[must_use] pub fn roots(self) -> [Vec; 2] { - let s_basis = to_symmetrical_basis_pair(self); - [s_basis.x.roots(), s_basis.y.roots()] + let (x, y) = self.parametric_polynomial(); + let x = poly_cool::Poly::new(x.coefficients().iter().copied()); + let y = poly_cool::Poly::new(y.coefficients().iter().copied()); + [x.roots_between(0.0, 1.0, 1e-8), y.roots_between(0.0, 1.0, 1e-8)] } /// Returns a list of lists of points representing the De Casteljau points for all iterations at the point `t` along the curve using De Casteljau's algorithm. From 1e9fd29d223edb7c182df5e36b1f3cc4792654ff Mon Sep 17 00:00:00 2001 From: Joe Neeman Date: Thu, 7 Aug 2025 22:25:28 -0500 Subject: [PATCH 4/4] Remove symmetrical_basis --- libraries/bezier-rs/src/lib.rs | 2 - libraries/bezier-rs/src/symmetrical_basis.rs | 622 ------------------- 2 files changed, 624 deletions(-) delete mode 100644 libraries/bezier-rs/src/symmetrical_basis.rs diff --git a/libraries/bezier-rs/src/lib.rs b/libraries/bezier-rs/src/lib.rs index 06faf15d84..d9442ba688 100644 --- a/libraries/bezier-rs/src/lib.rs +++ b/libraries/bezier-rs/src/lib.rs @@ -8,10 +8,8 @@ mod consts; mod poisson_disk; mod polynomial; mod subpath; -mod symmetrical_basis; mod utils; pub use bezier::*; pub use subpath::*; -pub use symmetrical_basis::*; pub use utils::{Cap, Join, SubpathTValue, TValue, TValueType}; diff --git a/libraries/bezier-rs/src/symmetrical_basis.rs b/libraries/bezier-rs/src/symmetrical_basis.rs deleted file mode 100644 index 800de69d9d..0000000000 --- a/libraries/bezier-rs/src/symmetrical_basis.rs +++ /dev/null @@ -1,622 +0,0 @@ -/* - * Modifications and Rust port copyright (C) 2024 by 0Hypercube. - * - * Original version by lib2geom: - * - * The entirety of this file is specially licensed under MPL 1.1 terms: - * - * Original Authors: - * Nathan Hurst - * Michael Sloan - * Marco Cecchetti - * MenTaLguY - * Michael Sloan - * Nathan Hurst - * Krzysztof Kosiński - * And additional authors listed in the version control history of the following files: - * - https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/sbasis.h - * - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis.cpp - * - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis-to-bezier.cpp - * - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier.cpp - * - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/solve-bezier.cpp - * - https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/solve-bezier-one-d.cpp - * - * Copyright (C) 2006-2015 Original Authors - * - * 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 crate::{Bezier, BezierHandles}; -use glam::DVec2; - -impl std::ops::Index for Bezier { - type Output = DVec2; - fn index(&self, index: usize) -> &Self::Output { - match &self.handles { - BezierHandles::Linear => [&self.start, &self.end][index], - BezierHandles::Quadratic { handle } => [&self.start, handle, &self.end][index], - BezierHandles::Cubic { handle_start, handle_end } => [&self.start, handle_start, handle_end, &self.end][index], - } - } -} - -// Note that the built in signum cannot be used as it does not handle 0 the same way as in the C code. -fn sign(x: f64) -> i8 { - if x > 0. { - 1 - } else if x < 0. { - -1 - } else { - 0 - } -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/sbasis.h#L70 -#[derive(Debug, Clone)] -pub(crate) struct SymmetricalBasis(pub Vec); - -impl SymmetricalBasis { - // https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis.cpp#L323 - #[must_use] - fn derivative(&self) -> SymmetricalBasis { - let mut c = SymmetricalBasis(vec![DVec2::ZERO; self.len()]); - if self.iter().all(|x| x.abs_diff_eq(DVec2::ZERO, 1e-5)) { - return c; - } - for k in 0..(self.len() - 1) { - let d = (2. * k as f64 + 1.) * (self[k][1] - self[k][0]); - - c[k][0] = d + (k as f64 + 1.) * self[k + 1][0]; - c[k][1] = d - (k as f64 + 1.) * self[k + 1][1]; - } - let k = self.len() - 1; - let d = (2. * k as f64 + 1.) * (self[k][1] - self[k][0]); - if d == 0. && k > 0 { - c.pop(); - } else { - c[k][0] = d; - c[k][1] = d; - } - c - } - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis-to-bezier.cpp#L86 - #[must_use] - pub fn to_bezier1d(&self) -> Bezier1d { - let sb = self; - assert!(!sb.is_empty()); - - let n; - let even; - - let mut q = sb.len(); - if sb[q - 1][0] == sb[q - 1][1] { - even = true; - q -= 1; - n = 2 * q; - } else { - even = false; - n = 2 * q - 1; - } - - let mut bz = Bezier1d(vec![0.; n + 1]); - for k in 0..q { - let mut tjk = 1.; - for j in k..(n - k) { - // j <= n-k-1 - bz[j] += tjk * sb[k][0]; - bz[n - j] += tjk * sb[k][1]; // n-k <-> [k][1] - tjk = binomial_increment_k(tjk, n - 2 * k - 1, j - k); - } - } - if even { - bz[q] += sb[q][0]; - } - // the resulting coefficients are with respect to the scaled Bernstein - // basis so we need to divide them by (n, j) binomial coefficient - let mut bcj = n as f64; - for j in 1..n { - bz[j] /= bcj; - bcj = binomial_increment_k(bcj, n, j); - } - bz[0] = sb[0][0]; - bz[n] = sb[0][1]; - bz - } - - fn normalize(&mut self) { - while self.len() > 1 && self.last().is_some_and(|x| x.abs_diff_eq(DVec2::ZERO, 1e-5)) { - self.pop(); - } - } - - #[must_use] - pub(crate) fn roots(&self) -> Vec { - match self.len() { - 0 => Vec::new(), - 1 => { - let mut res = Vec::new(); - let d = self[0].x - self[0].y; - if d != 0. { - let r = self[0].x / d; - if (0. ..=1.).contains(&r) { - res.push(r); - } - } - res - } - _ => { - let mut bz = self.to_bezier1d(); - let mut solutions = Vec::new(); - if bz.len() == 0 || bz.iter().all(|&x| (x - bz[0]).abs() < 1e-5) { - return solutions; - } - while bz[0] == 0. { - bz = bz.deflate(); - solutions.push(0.); - } - // Linear - if bz.len() - 1 == 1 { - if sign(bz[0]) != sign(bz[1]) { - let d = bz[0] - bz[1]; - if d != 0. { - let r = bz[0] / d; - if (0. ..=1.).contains(&r) { - solutions.push(r); - } - } - } - return solutions; - } - bz.find_bernstein_roots(&mut solutions, 0, 0., 1.); - - solutions.sort_by(f64::total_cmp); - - solutions - } - } - } -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis.cpp#L228 -impl std::ops::Mul for &SymmetricalBasis { - type Output = SymmetricalBasis; - fn mul(self, b: Self) -> Self::Output { - let a = self; - if a.iter().all(|x| x.abs_diff_eq(DVec2::ZERO, 1e-5)) || b.iter().all(|x| x.abs_diff_eq(DVec2::ZERO, 1e-5)) { - return SymmetricalBasis(vec![DVec2::ZERO]); - } - let mut c = SymmetricalBasis(vec![DVec2::ZERO; a.len() + b.len()]); - - for j in 0..b.len() { - for i in j..(a.len() + j) { - let tri = (b[j][1] - b[j][0]) * (a[i - j][1] - a[i - j][0]); - c[i + 1] += DVec2::splat(-tri); - } - } - for j in 0..b.len() { - for i in j..(a.len() + j) { - for dim in 0..2 { - c[i][dim] += b[j][dim] * a[i - j][dim]; - } - } - } - c.normalize(); - c - } -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis.cpp#L88 -impl std::ops::Add for SymmetricalBasis { - type Output = SymmetricalBasis; - fn add(self, b: Self) -> Self::Output { - let a = self; - let out_size = a.len().max(b.len()); - let min_size = a.len().min(b.len()); - let mut result = SymmetricalBasis(vec![DVec2::ZERO; out_size]); - for i in 0..min_size { - result[i] = a[i] + b[i]; - } - for i in min_size..a.len() { - result[i] = a[i]; - } - for i in min_size..b.len() { - result[i] = b[i]; - } - result - } -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis.cpp#L110 -impl std::ops::Sub for SymmetricalBasis { - type Output = SymmetricalBasis; - fn sub(self, b: Self) -> Self::Output { - let a = self; - let out_size = a.len().max(b.len()); - let min_size = a.len().min(b.len()); - let mut result = SymmetricalBasis(vec![DVec2::ZERO; out_size]); - for i in 0..min_size { - result[i] = a[i] - b[i]; - } - for i in min_size..a.len() { - result[i] = a[i]; - } - for i in min_size..b.len() { - result[i] = -b[i]; - } - result - } -} - -impl std::ops::Deref for SymmetricalBasis { - type Target = Vec; - fn deref(&self) -> &Self::Target { - &self.0 - } -} -impl std::ops::DerefMut for SymmetricalBasis { - fn deref_mut(&mut self) -> &mut Self::Target { - &mut self.0 - } -} - -#[derive(Debug, Clone)] -pub(crate) struct SymmetricalBasisPair { - pub x: SymmetricalBasis, - pub y: SymmetricalBasis, -} - -impl SymmetricalBasisPair { - #[must_use] - pub fn derivative(&self) -> Self { - Self { - x: self.x.derivative(), - y: self.y.derivative(), - } - } - - #[must_use] - pub fn dot(&self, other: &Self) -> SymmetricalBasis { - (&self.x * &other.x) + (&self.y * &other.y) - } - - #[must_use] - pub fn cross(&self, rhs: &Self) -> SymmetricalBasis { - (&self.x * &rhs.y) - (&self.y * &rhs.x) - } -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/sbasis.h#L337 -impl std::ops::Sub for SymmetricalBasisPair { - type Output = SymmetricalBasisPair; - fn sub(self, rhs: DVec2) -> Self::Output { - let sub = |a: &SymmetricalBasis, b: f64| { - if a.iter().all(|x| x.abs_diff_eq(DVec2::ZERO, 1e-5)) { - return SymmetricalBasis(vec![DVec2::splat(-b)]); - } - let mut result = a.clone(); - result[0] -= DVec2::splat(b); - result - }; - - Self { - x: sub(&self.x, rhs.x), - y: sub(&self.y, rhs.y), - } - } -} - -#[derive(Debug, Clone)] -pub struct Bezier1d(pub Vec); - -impl std::ops::Deref for Bezier1d { - type Target = Vec; - fn deref(&self) -> &Self::Target { - &self.0 - } -} -impl std::ops::DerefMut for Bezier1d { - fn deref_mut(&mut self) -> &mut Self::Target { - &mut self.0 - } -} - -impl Bezier1d { - const MAX_DEPTH: u32 = 53; - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier.cpp#L176 - #[must_use] - fn deflate(&self) -> Self { - let bz = self; - if bz.is_empty() { - return Bezier1d(Vec::new()); - } - let n = bz.len() - 1; - let mut b = Bezier1d(vec![0.; n]); - for i in 0..n { - b[i] = (n as f64 * bz[i + 1]) / (i as f64 + 1.) - } - b - } - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/bezier.h#L55 - /// Compute the value of a Bernstein-Bezier polynomial using a Horner-like fast evaluation scheme. - #[must_use] - fn value_at(&self, t: f64) -> f64 { - let bz = self; - let order = bz.len() - 1; - let u = 1. - t; - let mut bc = 1.; - let mut tn = 1.; - let mut tmp = bz[0] * u; - for i in 1..order { - tn *= t; - bc = bc * (order as f64 - i as f64 + 1.) / i as f64; - tmp = (tmp + tn * bc * bz[i]) * u; - } - tmp + tn * t * bz[bz.len() - 1] - } - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/solve-bezier.cpp#L258 - #[must_use] - fn secant(&self) -> f64 { - let bz = self; - let mut s = 0.; - let mut t = 1.; - let e = 1e-14; - let mut side = 0; - let mut r = 0.; - let mut fs = bz[0]; - let mut ft = bz[bz.len() - 1]; - - for _n in 0..100 { - r = (fs * t - ft * s) / (fs - ft); - if (t - s).abs() < e * (t + s).abs() { - return r; - } - - let fr = self.value_at(r); - - if fr * ft > 0. { - t = r; - ft = fr; - if side == -1 { - fs /= 2.; - } - side = -1; - } else if fs * fr > 0. { - s = r; - fs = fr; - if side == 1 { - ft /= 2.; - } - side = 1; - } else { - break; - } - } - r - } - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/bezier.h#L78 - fn casteljau_subdivision(&self, t: f64) -> [Self; 2] { - let v = self; - let order = v.len() - 1; - let mut left = v.clone(); - let mut right = v.clone(); - - // The Horner-like scheme gives very slightly different results, but we need - // the result of subdivision to match exactly with Bezier's valueAt function. - let val = v.value_at(t); - for i in (1..=order).rev() { - left[i - 1] = right[0]; - for j in i..v.len() { - right[j - 1] = right[j - 1] + ((right[j] - right[j - 1]) * t); - } - } - right[0] = val; - left[order] = right[0]; - [left, right] - } - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/bezier.cpp#L282 - fn derivative(&self) -> Self { - let bz = self; - if bz.len() - 1 == 1 { - return Bezier1d(vec![bz[1] - bz[0]]); - } - let mut der = Bezier1d(vec![0.; bz.len() - 1]); - - for i in 0..(bz.len() - 1) { - der[i] = (bz.len() - 1) as f64 * (bz[i + 1] - bz[i]); - } - der - } - - // https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/solve-bezier-one-d.cpp#L76 - /// given an equation in Bernstein-Bernstein form, find all roots between left_t and right_t - fn find_bernstein_roots(&self, solutions: &mut Vec, depth: u32, left_t: f64, right_t: f64) { - let bz = self; - let mut n_crossings = 0; - - let mut old_sign = sign(bz[0]); - for i in 1..bz.len() { - let sign = sign(bz[i]); - if sign != 0 { - if sign != old_sign && old_sign != 0 { - n_crossings += 1; - } - old_sign = sign; - } - } - // if last control point is zero, that counts as crossing too - if sign(bz[bz.len() - 1]) == 0 { - n_crossings += 1; - } - // no solutions - if n_crossings == 0 { - return; - } - // Unique solution - if n_crossings == 1 { - // Stop recursion when the tree is deep enough - return 1 solution at midpoint - if depth > Self::MAX_DEPTH { - let ax = right_t - left_t; - let ay = bz.last().unwrap() - bz[0]; - - solutions.push(left_t - ax * bz[0] / ay); - return; - } - - let r = bz.secant(); - solutions.push(r * right_t + (1. - r) * left_t); - return; - } - // solve recursively after subdividing control polygon - let o = bz.len() - 1; - let mut left = Bezier1d(vec![0.; o + 1]); - let mut right = bz.clone(); - let mut split_t = (left_t + right_t) * 0.5; - - // If subdivision is working poorly, split around the leftmost root of the derivative - if depth > 2 { - let dbz = bz.derivative(); - - let mut d_solutions = Vec::new(); - dbz.find_bernstein_roots(&mut d_solutions, 0, left_t, right_t); - d_solutions.sort_by(f64::total_cmp); - - let mut d_split_t = 0.5; - if !d_solutions.is_empty() { - d_split_t = d_solutions[0]; - split_t = left_t + (right_t - left_t) * d_split_t; - } - - [left, right] = bz.casteljau_subdivision(d_split_t); - } else { - // split at midpoint, because it is cheap - left[0] = right[0]; - for i in 1..bz.len() { - for j in 0..(bz.len() - i) { - right[j] = (right[j] + right[j + 1]) * 0.5; - } - left[i] = right[0]; - } - } - // Solution is exactly on the subdivision point - left.reverse(); - while right.len() - 1 > 0 && (right[0]).abs() <= 1e-10 { - // Deflate - right = right.deflate(); - left = left.deflate(); - solutions.push(split_t); - } - left.reverse(); - if right.len() - 1 > 0 { - left.find_bernstein_roots(solutions, depth + 1, left_t, split_t); - right.find_bernstein_roots(solutions, depth + 1, split_t, right_t); - } - } -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/choose.h#L61 -/// Given a multiple of binomial(n, k), modify it to the same multiple of binomial(n, k + 1). -#[must_use] -fn binomial_increment_k(b: f64, n: usize, k: usize) -> f64 { - b * (n as f64 - k as f64) / (k + 1) as f64 -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/include/2geom/choose.h#L52 -/// Given a multiple of binomial(n, k), modify it to the same multiple of binomial(n - 1, k). -#[must_use] -fn binomial_decrement_n(b: f64, n: usize, k: usize) -> f64 { - b * (n as f64 - k as f64) / n as f64 -} - -// https://gitlab.com/inkscape/lib2geom/-/blob/master/src/2geom/sbasis-to-bezier.cpp#L86 -#[must_use] -pub(crate) fn to_symmetrical_basis_pair(bezier: Bezier) -> SymmetricalBasisPair { - let n = match bezier.handles { - BezierHandles::Linear => 1, - BezierHandles::Quadratic { .. } => 2, - BezierHandles::Cubic { .. } => 3, - }; - let q = (n + 1) / 2; - let even = n % 2 == 0; - let mut sb = SymmetricalBasisPair { - x: SymmetricalBasis(vec![DVec2::ZERO; q + even as usize]), - y: SymmetricalBasis(vec![DVec2::ZERO; q + even as usize]), - }; - - let mut nck = 1.; - for k in 0..q { - let mut tjk = nck; - for j in k..q { - sb.x[j][0] += tjk * bezier[k].x; - sb.x[j][1] += tjk * bezier[n - k].x; - sb.y[j][0] += tjk * bezier[k].y; - sb.y[j][1] += tjk * bezier[n - k].y; - tjk = binomial_increment_k(tjk, n - j - k, j - k); - tjk = binomial_decrement_n(tjk, n - j - k, j - k + 1); - tjk = -tjk; - } - tjk = -nck; - for j in (k + 1)..q { - sb.x[j][0] += tjk * bezier[n - k].x; - sb.x[j][1] += tjk * bezier[k].x; - sb.y[j][0] += tjk * bezier[n - k].y; - sb.y[j][1] += tjk * bezier[k].y; - tjk = binomial_increment_k(tjk, n - j - k - 1, j - k - 1); - tjk = binomial_decrement_n(tjk, n - j - k - 1, j - k); - tjk = -tjk; - } - nck = binomial_increment_k(nck, n, k); - } - if even { - let mut tjk = if q % 2 == 1 { -1. } else { 1. }; - for k in 0..q { - sb.x[q][0] += tjk * (bezier[k].x + bezier[n - k].x); - sb.y[q][0] += tjk * (bezier[k].y + bezier[n - k].y); - tjk = binomial_increment_k(tjk, n, k); - tjk = -tjk; - } - sb.x[q][0] += tjk * bezier[q].x; - sb.x[q][1] = sb.x[q][0]; - sb.y[q][0] += tjk * bezier[q].y; - sb.y[q][1] = sb.y[q][0]; - } - sb.x[0][0] = bezier[0].x; - sb.x[0][1] = bezier[n].x; - sb.y[0][0] = bezier[0].y; - sb.y[0][1] = bezier[n].y; - - sb -} - -#[cfg(test)] -mod tests { - use super::*; - - #[test] - fn find_bernstein_roots() { - let bz = Bezier1d(vec![50., -100., 170.]); - let mut solutions = Vec::new(); - bz.find_bernstein_roots(&mut solutions, 0, 0., 1.); - - solutions.sort_by(f64::total_cmp); - for &t in &solutions { - assert!(bz.value_at(t,).abs() < 1e-5, "roots should be roots {} {}", t, bz.value_at(t,)); - } - } -}