Use poly-cool as the root finder in bezier-rs

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

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 }

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.
 	/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#bezier/project/solo" title="Project Demo"></iframe>
 	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.)))));

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<f64>; 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 {
 	/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#bezier/tangents-to-point/solo" title="Tangents to Point Demo"></iframe>
 	#[must_use]
 	pub fn tangents_to_point(self, point: DVec2) -> Vec<f64> {
-		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 {
 	/// <iframe frameBorder="0" width="100%" height="300px" src="https://graphite.rs/libraries/bezier-rs#bezier/normals-to-point/solo" title="Normals to Point Demo"></iframe>
 	#[must_use]
 	pub fn normals_to_point(self, point: DVec2) -> Vec<f64> {
-		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.

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};