From 6a10bcc3fc752d8604989eb2e18d55082021f161 Mon Sep 17 00:00:00 2001 From: 0HyperCube <78500760+0HyperCube@users.noreply.github.com> Date: Sun, 7 Jan 2024 14:23:04 +0000 Subject: [PATCH] Bezier-rs: add Tangents/Normals to Point functions (#1547) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add tangent/normal to point to bézier-rs * Add license information and general cleanup * Bump version * Add iframe demos --------- Co-authored-by: Keavon Chambers --- Cargo.lock | 2 +- libraries/bezier-rs/Cargo.toml | 6 +- libraries/bezier-rs/src/bezier/solvers.rs | 79 +++ libraries/bezier-rs/src/lib.rs | 2 + libraries/bezier-rs/src/symmetrical_basis.rs | 612 ++++++++++++++++++ .../src/components/BezierDemo.ts | 6 + .../src/features/bezier-features.ts | 13 + .../other/bezier-rs-demos/wasm/src/bezier.rs | 25 + website/other/bezier-rs-demos/wasm/src/lib.rs | 12 + 9 files changed, 753 insertions(+), 4 deletions(-) create mode 100644 libraries/bezier-rs/src/symmetrical_basis.rs diff --git a/Cargo.lock b/Cargo.lock index 2399b21a..d509c759 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -559,7 +559,7 @@ checksum = "35636a1494ede3b646cc98f74f8e62c773a38a659ebc777a2cf26b9b74171df9" [[package]] name = "bezier-rs" -version = "0.2.0" +version = "0.3.0" dependencies = [ "dyn-any", "glam", diff --git a/libraries/bezier-rs/Cargo.toml b/libraries/bezier-rs/Cargo.toml index 2497d0d5..2cbcb31f 100644 --- a/libraries/bezier-rs/Cargo.toml +++ b/libraries/bezier-rs/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "bezier-rs" -version = "0.2.0" +version = "0.3.0" rust-version = "1.66.0" edition = "2021" authors = ["Graphite Authors "] @@ -16,5 +16,5 @@ documentation = "https://graphite.rs/libraries/bezier-rs/" [dependencies] glam = { version = "0.24", features = ["serde"] } -dyn-any = { workspace = true, optional = true } -serde = { workspace = true, optional = true } +dyn-any = { version = "0.3.0", path = "../dyn-any", optional = true } +serde = { version = "1.0", workspace = true, optional = true } diff --git a/libraries/bezier-rs/src/bezier/solvers.rs b/libraries/bezier-rs/src/bezier/solvers.rs index 326d24a0..4e8d7edb 100644 --- a/libraries/bezier-rs/src/bezier/solvers.rs +++ b/libraries/bezier-rs/src/bezier/solvers.rs @@ -1,11 +1,19 @@ use super::*; use crate::utils::{solve_cubic, solve_quadratic, TValue}; +use crate::{to_symmetrical_basis_pair, SymmetricalBasis}; use glam::DMat2; use std::ops::Range; /// Functionality that solve for various curve information such as derivative, tangent, intersect, etc. 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()] + } + /// 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. /// The `i`th element of the list represents the set of points in the `i`th iteration. /// More information on the algorithm can be found in the [De Casteljau section](https://pomax.github.io/bezierinfo/#decasteljau) in Pomax's primer. @@ -72,12 +80,32 @@ impl Bezier { } } + /// Find the `t`-value(s) such that the tangent(s) at `t` pass through the specified point. + /// + #[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) + } + /// Returns a normalized unit vector representing the direction of the normal at the point `t` along the curve. /// pub fn normal(&self, t: TValue) -> DVec2 { self.tangent(t).perp() } + /// Find the `t`-value(s) such that the normal(s) at `t` pass through the specified point. + /// + #[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) + } + /// Returns the curvature, a scalar value for the derivative at the point `t` along the curve. /// Curvature is 1 over the radius of a circle with an equivalent derivative. /// @@ -582,6 +610,27 @@ mod tests { assert_eq!(cubic.tangent(TValue::Parametric(1.)), DVec2::new(30., 120.).normalize()); } + #[test] + fn tangent_at_point() { + let validate = |bz: Bezier, p: DVec2| { + let solutions = bz.tangents_to_point(p); + assert_ne!(solutions.len(), 0); + for t in solutions { + let pos = bz.evaluate(TValue::Parametric(t)); + let expected_tangent = (pos - p).normalize(); + let tangent = bz.tangent(TValue::Parametric(t)); + assert!(expected_tangent.perp_dot(tangent).abs() < 0.2, "Expected tangent {expected_tangent} found {tangent} pos {pos}") + } + }; + let bz = Bezier::from_quadratic_coordinates(55., 50., 165., 30., 185., 170.); + let p = DVec2::new(193., 83.); + validate(bz, p); + + let bz = Bezier::from_cubic_coordinates(55., 30., 18., 139., 175., 30., 185., 160.); + let p = DVec2::new(127., 121.); + validate(bz, p); + } + #[test] fn test_normal() { // Test normals at start and end points of each Bezier curve type @@ -604,6 +653,36 @@ mod tests { assert_eq!(cubic.normal(TValue::Parametric(1.)), DVec2::new(-120., 30.).normalize()); } + #[test] + fn normal_at_point() { + let validate = |bz: Bezier, p: DVec2| { + let solutions = bz.normals_to_point(p); + assert_ne!(solutions.len(), 0); + for t in solutions { + let pos = bz.evaluate(TValue::Parametric(t)); + let expected_normal = (pos - p).normalize(); + let normal = bz.normal(TValue::Parametric(t)); + assert!(expected_normal.perp_dot(normal).abs() < 0.2, "Expected normal {expected_normal} found {normal} pos {pos}") + } + }; + + let bz = Bezier::from_linear_coordinates(50., 50., 100., 100.); + let p = DVec2::new(100., 50.); + validate(bz, p); + + let bz = Bezier::from_quadratic_coordinates(55., 50., 165., 30., 185., 170.); + let p = DVec2::new(193., 83.); + validate(bz, p); + + let bz = Bezier::from_cubic_coordinates(55., 30., 18., 139., 175., 30., 185., 160.); + let p = DVec2::new(127., 121.); + validate(bz, p); + + let bz = Bezier::from_cubic_coordinates(55.0, 30.0, 85.0, 140.0, 175.0, 30.0, 185.0, 160.0); + let p = DVec2::new(17., 172.); + validate(bz, p); + } + #[test] fn test_curvature() { let p1 = DVec2::new(10., 10.); diff --git a/libraries/bezier-rs/src/lib.rs b/libraries/bezier-rs/src/lib.rs index bc8a8a4f..afe6c85c 100644 --- a/libraries/bezier-rs/src/lib.rs +++ b/libraries/bezier-rs/src/lib.rs @@ -5,8 +5,10 @@ pub(crate) mod compare; mod bezier; mod consts; 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 new file mode 100644 index 00000000..8b11f514 --- /dev/null +++ b/libraries/bezier-rs/src/symmetrical_basis.rs @@ -0,0 +1,612 @@ +/* + * 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 LGPL or 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], + } + } +} + +// 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 bz[0].signum() != bz[1].signum() { + 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<'a> std::ops::Mul for &'a 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.0 - 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 = bz[0].signum(); + for i in 1..bz.len() { + let sign = bz[i].signum(); + 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 bz[bz.len() - 1].signum() == 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.0, -100.0, 170.0]); + 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,)); + } + } +} diff --git a/website/other/bezier-rs-demos/src/components/BezierDemo.ts b/website/other/bezier-rs-demos/src/components/BezierDemo.ts index 3a48faa4..f5e4772d 100644 --- a/website/other/bezier-rs-demos/src/components/BezierDemo.ts +++ b/website/other/bezier-rs-demos/src/components/BezierDemo.ts @@ -39,6 +39,9 @@ class BezierDemo extends HTMLElement implements Demo { sliderUnits!: Record; + // Avoids "recursive use of an object detected which would lead to unsafe aliasing in rust" error when moving mouse fast. + locked!: boolean; + async connectedCallback() { this.title = this.getAttribute("title") || ""; this.points = JSON.parse(this.getAttribute("points") || "[]"); @@ -85,6 +88,8 @@ class BezierDemo extends HTMLElement implements Demo { } onMouseMove(event: MouseEvent) { + if (this.locked) return; + this.locked = true; const mx = event.offsetX; const my = event.offsetY; const figure = event.currentTarget as HTMLElement; @@ -96,6 +101,7 @@ class BezierDemo extends HTMLElement implements Demo { } else if (this.triggerOnMouseMove) { this.drawDemo(figure, [mx, my]); } + this.locked = false; } getSliderUnit(sliderValue: number, variable: string): string { diff --git a/website/other/bezier-rs-demos/src/features/bezier-features.ts b/website/other/bezier-rs-demos/src/features/bezier-features.ts index 7b907262..fb56b9e9 100644 --- a/website/other/bezier-rs-demos/src/features/bezier-features.ts +++ b/website/other/bezier-rs-demos/src/features/bezier-features.ts @@ -126,6 +126,13 @@ const bezierFeatures = { }, }, }, + "tangents-to-point": { + name: "Tangents To Point", + callback: (bezier: WasmBezierInstance, _: Record, mouseLocation?: [number, number]): string => + mouseLocation ? bezier.tangents_to_point(mouseLocation[0], mouseLocation[1]) : bezier.to_svg(), + triggerOnMouseMove: true, + demoOptions: { Linear: { disabled: true } }, + }, normal: { name: "Normal", callback: (bezier: WasmBezierInstance, options: Record, _: undefined): string => bezier.normal(options.t, BEZIER_T_VALUE_VARIANTS[options.TVariant]), @@ -135,6 +142,12 @@ const bezierFeatures = { }, }, }, + "normals-to-point": { + name: "Normals To Point", + callback: (bezier: WasmBezierInstance, _: Record, mouseLocation?: [number, number]): string => + mouseLocation ? bezier.normals_to_point(mouseLocation[0], mouseLocation[1]) : bezier.to_svg(), + triggerOnMouseMove: true, + }, curvature: { name: "Curvature", callback: (bezier: WasmBezierInstance, options: Record, _: undefined): string => bezier.curvature(options.t, BEZIER_T_VALUE_VARIANTS[options.TVariant]), diff --git a/website/other/bezier-rs-demos/wasm/src/bezier.rs b/website/other/bezier-rs-demos/wasm/src/bezier.rs index 7bf022b4..9c1b13ec 100644 --- a/website/other/bezier-rs-demos/wasm/src/bezier.rs +++ b/website/other/bezier-rs-demos/wasm/src/bezier.rs @@ -300,6 +300,31 @@ impl WasmBezier { let content = format!("{bezier}{}", draw_line(projected_point.x, projected_point.y, x, y, RED, 1.),); wrap_svg_tag(content) } + pub fn tangents_to_point(&self, x: f64, y: f64) -> String { + let bezier = self.get_bezier_path(); + let mut content = String::new(); + for t in self.0.tangents_to_point(DVec2::new(x, y)) { + let point = self.0.evaluate(TValue::Parametric(t)); + content += &draw_line(x, y, point.x, point.y, RED, 1.); + } + use std::fmt::Write; + write!(content, "{bezier}").unwrap(); + + wrap_svg_tag(content) + } + + pub fn normals_to_point(&self, x: f64, y: f64) -> String { + let bezier = self.get_bezier_path(); + let mut content = String::new(); + for t in self.0.normals_to_point(DVec2::new(x, y)) { + let point = self.0.evaluate(TValue::Parametric(t)); + content += &draw_line(x, y, point.x, point.y, RED, 1.); + } + use std::fmt::Write; + write!(content, "{bezier}").unwrap(); + + wrap_svg_tag(content) + } pub fn local_extrema(&self) -> String { let local_extrema: [Vec; 2] = self.0.local_extrema(); diff --git a/website/other/bezier-rs-demos/wasm/src/lib.rs b/website/other/bezier-rs-demos/wasm/src/lib.rs index 3363017b..7b8d849c 100644 --- a/website/other/bezier-rs-demos/wasm/src/lib.rs +++ b/website/other/bezier-rs-demos/wasm/src/lib.rs @@ -4,7 +4,9 @@ mod svg_drawing; mod utils; use wasm_bindgen::prelude::*; + pub static LOGGER: WasmLog = WasmLog; +thread_local! { pub static HAS_CRASHED: std::cell::RefCell = std::cell::RefCell::new(false); } /// Initialize the backend #[wasm_bindgen(start)] @@ -12,6 +14,16 @@ pub fn init() { // Set up the logger with a default level of debug log::set_logger(&LOGGER).expect("Failed to set logger"); log::set_max_level(log::LevelFilter::Trace); + + fn panic_hook(info: &core::panic::PanicInfo) { + // Skip if we have already panicked + if HAS_CRASHED.with(|cell| cell.replace(true)) { + return; + } + log::error!("{}", info); + } + + std::panic::set_hook(Box::new(panic_hook)); } /// Logging to the JS console