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/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 320e772d9f..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.
@@ -105,10 +106,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.
@@ -121,10 +133,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.
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,));
- }
- }
-}