From 0deda6cfd8da45720e930d281781f51b3c7cea40 Mon Sep 17 00:00:00 2001 From: Raph Levien Date: Fri, 1 Mar 2024 07:54:28 -0800 Subject: [PATCH] Add new hyperbezier math MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This commit adds a first implementation of the new hyperbezier curve. It doesn't yet offer much help in determining the raw parameters, but does produce nice Bézier paths through curve fitting. --- Cargo.toml | 4 +- README.md | 24 ++- examples/hyperbez.rs | 14 ++ src/hyperbezier.rs | 411 ++++++++++------------------------------- src/hyperbezier_old.rs | 349 ++++++++++++++++++++++++++++++++++ src/lib.rs | 5 +- src/simple_spline.rs | 2 +- src/spline.rs | 6 +- 8 files changed, 484 insertions(+), 331 deletions(-) create mode 100644 examples/hyperbez.rs create mode 100644 src/hyperbezier_old.rs diff --git a/Cargo.toml b/Cargo.toml index 302a258..59c1607 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,7 @@ name = "spline" version = "0.3.0" license = "MIT/Apache-2.0" authors = ["Raph Levien "] -edition = "2018" +edition = "2021" keywords = ["graphics", "curve", "curves", "bezier", "spline"] repository = "https://github.com/linebender/spline" description = "A spline for interactive 2D curve design" @@ -11,7 +11,7 @@ readme = "README.md" categories = ["graphics"] [dependencies] -kurbo = "0.9" +kurbo = "0.11" serde_ = { version = "1.0.117", package="serde", features = ["derive"], optional = true } [dev-dependencies] diff --git a/README.md b/README.md index 5c3c462..beaf97c 100644 --- a/README.md +++ b/README.md @@ -4,22 +4,32 @@ This crate implements a new spline designed and optimized for interactive design The work builds on previous iterations, notably the [Spiro] spline, and then another [research spline]. -## Hyperbeziers +## Status: in transition -The major innovation of this spline is the "hyperbezier" curve family. Like cubic Béziers and the Spiro curve, it is a four-parameter curve family. In fact, it's closely based on Spiro and there is significant overlap of the parameter space, including Euler spirals. +This crate is in transition from supporting the first draft of hyperbeziers to the second. The first draft had a reasonably clean mapping from Bézier control points to curve parameters (something still in progress for the second), but a number of shortcomings: -There is a significant difference, however. In the Spiro curve family, curvature is bounded, so it is not capable of cusp-like behavior. Rather, when "pushed," Spiro tends to wiggly, Shmoo-like shapes. Béziers are of course capable of high curvature regions, as are elastica when placed under very high tension. +* It was not closed under subdivision -A good way to parametrize the hyperbezier is by tangent angle and "tension," which correlates strongly with curvature at the endpoint. At low tension, the hyperbezier is equivalent to the Spiro curve. A natural tension value produces the Euler spiral (curvature is a linear function of arclength). But for higher tension values, a different function takes over, which approaches a cusp at the endpoint as tension increases. +* It was not capable of superellipse-like shapes -Unlike Béziers, the cusp happens *only* at the endpoint. Curvature maxima in the interior of a curve are ugly. With the hyperbezier, if the designer wants a sharp curvature maximum, simply place an on-curve point there. +* The shape when high tension on one side and low tension on the other was lumpy, with an undesirable curvature minimum -A particular strength of the hyperbezier is smooth (G2-continuous) transitions from straight to curved sections. The hyperbezier is capable (unlike a cubic Bézier) of an S-shaped curve with zero curvature at both ends. It's also capable of a wide range of Euler spiral like behavior where one end has zero curvature and the other is a nice rounded shape (in the general case a designer would use at least two Béziers to create this effect). +The new hyperbezier addresses all these problems. -The name "hyperbezier" clearly references its roots in the [cubic Bézier][A Primer on Bézier Curves], and the "hyper" part is a reference to the fact that the Euler spiral, an important section of its parameter space, is an instance of the [Hypergeometric function]. +## Hyperbeziers + +The major innovation of this spline is the "hyperbezier" curve family. Like cubic Béziers and the Spiro curve, it is a four-parameter curve family. In fact, it's based on Spiro and there is significant overlap of the parameter space, including Euler spirals. + +There is a significant difference, however. In the Spiro curve family, curvature is bounded, so it is not capable of cusp-like behavior. Rather, when "pushed," Spiro tends to wiggly, Shmoo-like shapes. Béziers are of course capable of high curvature regions, as are elastica when placed under very high tension. + +A good way to parametrize the hyperbezier is by tangent angle and "tension," which correlates strongly with curvature at the endpoint. At low tension, the hyperbezier has similar behavior to the Spiro curve. A natural tension value produces the Euler spiral (curvature is a linear function of arclength). But for higher tension values, curvature rises more sharply near the endpoints. In addition, the hyperbezier has a region of parameter space which produces curves resembling hyperbolas (also usable as a quadrant of a superellipse-like shape). + +The name "hyperbezier" references its roots in the [cubic Bézier][A Primer on Bézier Curves], and the "hyper" part is a reference to the fact that the Euler spiral, an important section of its parameter space, is an instance of the [Hypergeometric function]. Additionally, the "hyper" can be seen as a reference to hyperbola-like shapes, which it approximates and a cubic Bézier cannot. ## Focus on UX +Status: this section may be revised, as the new hyperbezier has not yet been hooked up to an interactive example. + A persistent challenge with spline-based curve design is getting the UX right. Bézier curves are not easy to master, but the [pen tool] has become highly refined over time, and is an extremely productive interface for designers. A major motivation for this work is to retain the good parts of the Bézier UX. In particular, the "control handle" maps to hyperbezier parameters in a natural, intuitive way. The tangent angle is obvious, and tension similarly dependent on the length of the control arm. So it's completely valid to use hyperbeziers simply as a drop-in replacement for Béziers. diff --git a/examples/hyperbez.rs b/examples/hyperbez.rs new file mode 100644 index 0000000..abf2408 --- /dev/null +++ b/examples/hyperbez.rs @@ -0,0 +1,14 @@ +use kurbo::{fit_to_bezpath, ParamCurve, Point}; +use spline::hyperbezier::{HyperbezParams, Hyperbezier}; + +fn main() { + let params = HyperbezParams::new(-10., 5., 2.0, -2.0); + let p0 = Point::new(100., 100.); + let p1 = Point::new(300., 200.); + let hb = Hyperbezier::from_points_params(params, p0, p1); + println!("{hb:?}"); + let p = fit_to_bezpath(&hb, 1.0); + println!("{}", p.to_svg()); + let p = fit_to_bezpath(&hb.subsegment(0.1 .. 0.8), 1.0); + println!("{}", p.to_svg()); +} \ No newline at end of file diff --git a/src/hyperbezier.rs b/src/hyperbezier.rs index 849eabe..cc4e5bb 100644 --- a/src/hyperbezier.rs +++ b/src/hyperbezier.rs @@ -1,349 +1,128 @@ -//! The math for the hyperbezier curve family. +//! Representation and computation of hyperbeziers. -use kurbo::common as coeffs; -use kurbo::{Affine, BezPath, PathEl, Point, Vec2}; +use kurbo::{common::GAUSS_LEGENDRE_COEFFS_32, CurveFitSample, ParamCurve, ParamCurveFit, Point, Vec2}; -use crate::util; +#[derive(Clone, Copy, Debug)] +pub struct HyperbezParams { + a: f64, + b: f64, + c: f64, + d: f64, -/// Parameters for a hyperbezier curve. -/// -/// A hyperbezier is a curve defined by curvature as a function of arclength. -/// It is similar to the Spiro curve in this way, but for some of the parameter -/// space the function is different. -/// -/// The parameter space is four dimensional. It is broken down symmetrically -/// into two parameters that predominantly affect one side of the curve, and -/// the curvature contributions are added: -/// -/// k(s) = k0 * f(bias0, 1 - s) + k1 * f(bias1, s) -/// -/// The "f" function takes a bias parameter, which can also be thought of as -/// tension. This value ranges from about -1 to exactly 2, with 2 representing -/// a cusp (infinitely high tension at the endpoint). For bias values less than -/// 1, it is defined thus: -/// -/// f(bias, s) = s + 6 * (1 - bias) * (s^2 - s^3 - s) -/// -/// For bias values greater than one, it is defined thus: -/// -/// f(bias, s) = c * s / (1 + (bias - 1) * s) ^ 2 -/// -/// Here, c is a normalization term chosen so that the integral of f from s=0 -/// to s=1 is 1. -/// -/// A few observation. If both bias values are 1, then the curve is an Euler -/// spiral. If both bias values are less than 1, then curvature is a cubic -/// polynomial as a function of arclength, so it is a Spiro curve. -#[derive(Copy, Clone, Debug)] -pub struct HyperBezier { - pub k0: f64, - pub bias0: f64, - pub k1: f64, - pub bias1: f64, + th_a: f64, + th_b: f64, } -/// An intermediate parametrization of the curve family. -/// -/// Here, angles are given relative to the chord, but the bias parameters -/// are the same as for `HyperBezier`. -#[derive(Copy, Clone, Debug)] -pub struct ThetaParams { - pub th0: f64, - pub bias0: f64, - pub th1: f64, - pub bias1: f64, +#[derive(Clone, Copy, Debug)] +pub struct Hyperbezier { + params: HyperbezParams, + p0: Point, + p1: Point, + scale_rot: Vec2, } -/// Result of measuring the curve. -/// -/// The `th0` and `th1` values are defined so that if they are have the -/// same sign, the curve is convex, but if they are opposite signs, it is -/// an "s" shape. -#[derive(Copy, Clone, Debug)] -pub struct HyperBezierResult { - /// Tangent angle from the chord to the curve at the start point. - pub th0: f64, - /// Tangent angle from the chord to the curve at the end point. - pub th1: f64, - /// Length of the chord assuming total arclength = 1. - pub chord: f64, - pub k0: f64, - pub k1: f64, -} - -impl HyperBezier { - /// Compute the angle for the given parameter. - /// - /// The argument is an arclength parametrization, ranging from 0 to 1. - /// - /// The returned angle is relative only, in other words there could be an - /// arbitrary rotation of the entire curve. - pub fn compute_theta(&self, s: f64) -> f64 { - self.k1 * integrate_basis(self.bias1, s) - self.k0 * integrate_basis(self.bias0, 1.0 - s) - } - - /// Compute the endpoint tangent angles and the chord length. - pub fn compute(&self) -> HyperBezierResult { - let integral = self.integrate(0.0, 1.0, 24); - let th_chord = integral.atan2(); - let chord = integral.hypot(); - let th0 = th_chord - self.compute_theta(0.0); - let th1 = self.compute_theta(1.0) - th_chord; - let k0 = chord * self.k0 * compute_k(self.bias0); - let k1 = chord * self.k1 * compute_k(self.bias1); - HyperBezierResult { - th0, - th1, - chord, - k0, - k1, - } +impl HyperbezParams { + /// Create a new hyperbezier with the given parameters. + pub fn new(a: f64, b: f64, c: f64, d: f64) -> Self { + // TODO: numerical stability issues when c is very small + // (which can happen for pure Euler spiral) or 4c-d^2 is + // very small (hyperbola-like with very high curvature). + let denom = 2. / (4. * c - d * d); + let dd = d * denom; + let th_a = (2. * b * c - d * a) * denom; + let th_b = b * dd - a * (1. + 0.5 * d * dd) / c; + HyperbezParams { a, b, c, d, th_a, th_b } } - fn integrate(&self, t0: f64, t1: f64, order: usize) -> Vec2 { - let c = match order { - 3 => coeffs::GAUSS_LEGENDRE_COEFFS_3, - 5 => coeffs::GAUSS_LEGENDRE_COEFFS_5, - 7 => coeffs::GAUSS_LEGENDRE_COEFFS_7, - 9 => coeffs::GAUSS_LEGENDRE_COEFFS_9, - 11 => coeffs::GAUSS_LEGENDRE_COEFFS_11, - 24 => coeffs::GAUSS_LEGENDRE_COEFFS_24, - _ => panic!("don't have coefficients for {}", order), - }; - let mut result = Vec2::ZERO; - let tm = 0.5 * (t1 + t0); - let dt = 0.5 * (t1 - t0); - for (wi, xi) in c { - let t = tm + dt * xi; - let th = self.compute_theta(t); - result += *wi * Vec2::from_angle(th); - } - dt * result - } - - /// Render to a [`BezPath`]. - pub fn render(&self, n: usize) -> BezPath { - self.render_elements(n).collect() - } - - /// Render to bezier elements. + /// Determine the angle for the given parameter. /// - /// The current algorithm just does a fixed subdivision based on arclength, - /// but should be adaptive in several ways; more subdivision for twistier - /// curves, and also more sophisticated parametrization (important as tension - /// increases). - pub fn render_elements<'a>(&'a self, n: usize) -> impl Iterator + 'a { - let order = 24; - let v = self.integrate(0.0, 1.0, order); - let a = Affine::new([v.x, v.y, -v.y, v.x, 0.0, 0.0]).inverse(); - let step = 1.0 / (n as f64); - fn calc_t(bias: f64) -> f64 { - if bias >= 1.0 { - (2.0 - bias).sqrt() * (1.0 / 3.0) - } else { - // Possibly this should increase for low tension curves, but that's not - // obvious. - 1.0 / 3.0 - } - } - let t1 = calc_t(self.bias0); - let t2 = 1.0 - calc_t(self.bias1); - let mut last_p = Point::ZERO; - let mut last_v = step * t1 * Vec2::from_angle(self.compute_theta(0.0)); - let mut i = 0; - let mut first = Some(PathEl::MoveTo(last_p)); - std::iter::from_fn(move || { - if let Some(first) = first.take() { - return Some(first); - } - i += 1; - if i <= n { - let u = (i as f64) * step; - let um = 1.0 - u; - let t = 3.0 * u * um * (um * t1 + u * t2) + u.powi(3); - let p = self.integrate(0.0, t, order).to_point(); - let p1 = last_p + last_v; - let dt = um * um * t1 + 2.0 * u * um * (t2 - t1) + u * u * (1.0 - t2); - let v = step * dt * Vec2::from_angle(self.compute_theta(t)); - let p2 = p - v; - let next = PathEl::CurveTo(a * p1, a * p2, a * p); - last_v = v; - last_p = p; - Some(next) - } else { - None - } - }) + /// This can be interpreted as a Whewell representation of the + /// curve. The `t` parameter ranges from 0 to 1, and the returned + /// value is 0 for `t = 0`. + fn theta(&self, t: f64) -> f64 { + let q = self.c * t * t + self.d * t + 1.; + (self.th_a * t + self.th_b) / q.sqrt() - self.th_b } - /// Suggest a number of subdivisions for rendering. + /// Evaluate the position of the raw curve. /// - /// This is a bit of a hacky heuristic. - pub fn render_subdivisions(&self) -> usize { - 2 + (self.k0.abs() + self.k1.abs()).floor() as usize - } - - /// Solve for curve params, given theta params. - pub fn solve_for_theta(params: &ThetaParams) -> HyperBezier { - let ThetaParams { - th0, - bias0, - th1, - bias1, - } = *params; - let mut dth = 0.0; - let mut lastxy: Option<(f64, f64)> = None; - const N: usize = 10; - for i in 0..N { - let params = HyperBezier { - k0: th0 + 0.5 * dth, - bias0, - k1: th1 - 0.5 * dth, - bias1, - }; - if i == N - 1 { - return params; - } - let result = params.compute(); - let th_err = util::mod_tau(th0 - th1 - (result.th0 - result.th1)); - if th_err.abs() < 1e-3 { - return params; - } - // Secant method - let nextxy = (dth, th_err); - let delta = if let Some(lastxy) = lastxy { - (nextxy.0 - lastxy.0) / (nextxy.1 - lastxy.1) - } else { - -0.5 - }; - dth -= delta * th_err; - lastxy = Some(nextxy); + /// This is simply the integral of the Whewell representation, + /// so that the total arc length is unit, and the initial tangent + /// is horizontal. + fn integrate(&self, t: f64) -> Vec2 { + // TODO: improve accuracy by subdividing in near-cusp cases + let mut xy = Vec2::ZERO; + let u0 = 0.5 * t; + for (wi, xi) in GAUSS_LEGENDRE_COEFFS_32 { + let u = u0 + u0 * xi; + xy += *wi * Vec2::from_angle(self.theta(u)); } - unreachable!() + u0 * xy } +} - /// Solve for curve params, given bezier control points. - /// - /// The points are given relative to p0 at (0, 0) and p3 at - /// (1, 0). - pub fn solve(p1: Point, p2: Point) -> HyperBezier { - let (th0, bias0) = Self::params_for_v(p1.to_vec2()); - let (th1, bias1) = Self::params_for_v(Point::new(1.0, 0.0) - p2); - // TODO: signs feel reversed here, but it all works out in the end. - let theta_params = ThetaParams { - th0: -th0, - bias0, - th1, - bias1, - }; - Self::solve_for_theta(&theta_params) +impl Hyperbezier { + /// Create a new hyperbezier curve with given parameters and end points. + pub fn from_points_params(params: HyperbezParams, p0: Point, p1: Point) -> Self{ + let uv = params.integrate(1.0); + let uv_scaled = uv / uv.length_squared(); + let d = p1 - p0; + let scale_rot = Vec2::new(uv_scaled.dot(d), uv_scaled.cross(d)); + Hyperbezier { params, p0, p1, scale_rot } } +} - /// Determine params for a control arm. - /// - /// Return values are theta and bias. - pub(crate) fn params_for_v(v: Vec2) -> (f64, f64) { - let th = v.atan2(); - // This formula ensures that bezier parameters approximating - // a circular arc map to a bias of 1.0. - let a = v.hypot() * 1.5 * (th.cos() + 1.0); - let bias = if a < 1.0 { - 2.0 - a * a +impl ParamCurve for Hyperbezier { + fn eval(&self, t: f64) -> Point { + if t == 1.0 { + self.p1 } else { - 1.0 + 2.0 * (0.5 * (1.0 - a)).tanh() - }; - (th, bias) + let s = self.scale_rot; + let uv = self.params.integrate(t); + self.p0 + Vec2::new(s.x * uv.x - s.y * uv.y, s.x * uv.y + s.y * uv.x) + } } - /// Determine control arm position from params. - /// - /// This function should be the inverse of `params_for_v` - pub(crate) fn v_for_params(th: f64, bias: f64) -> Vec2 { - let a = if bias >= 1.0 { - (2.0 - bias).sqrt() - } else { - // Bias may not be bounded from below, we probably want - // to rethink this... - 1.0 - 2.0 * (0.5 * (bias - 1.0)).atanh() - }; - let len = a / (1.5 * (th.cos() + 1.0)); - len * Vec2::from_angle(th) + fn start(&self) -> Point { + self.p0 } -} -const MAX_A: f64 = 1.0 - 1e-4; + fn end(&self) -> Point { + self.p1 + } -/// Compute integral of basis function. -/// -/// The integral of the basis function can be represented as a reasonably -/// simple closed-form analytic formula. -/// -/// Note: this is normalized so that f(1) - f(0) = 1. -/// -/// This is oriented for the rightmost control point. -fn integrate_basis(bias: f64, s: f64) -> f64 { - if bias <= 1.0 { - let iy0 = 4.0 * s.powi(3) - 3.0 * s.powi(4); - let iy1 = s.powi(2); - iy0 + bias * (iy1 - iy0) - } else if bias < 1.0002 { - // This is a more numerically robust approximation to the - // exact analytical formula in the next clause. - let b = (bias - 1.0) * (4.0 / 3.0); - (1.0 - b) * s.powi(2) + b * s.powi(3) - } else { - let a = (bias - 1.0).min(MAX_A); - let norm = 1.0 / (1.0 - a) + (1.0 - a).ln() - 1.0; - (1.0 / (1.0 - a * s) + (1.0 - a * s).ln() - 1.0) / norm + fn subsegment(&self, range: std::ops::Range) -> Self { + let (t0, t1) = (range.start, range.end); + let dt = t1 - t0; + let a = self.params.a * dt; + let b = self.params.b + self.params.a * t0; + let c = self.params.c * dt * dt; + let d = (self.params.d + 2. * self.params.c * t0) * dt; + let e = self.params.c * t0 * t0 + self.params.d * t0 + 1.; + let s = 1. / e; + let ps = dt * s * s.sqrt(); + let params = HyperbezParams::new(a * ps, b * ps, c * s, d * s); + let p0 = self.eval(t0); + let p1 = self.eval(t1); + Hyperbezier::from_points_params(params, p0, p1) } } -/// Compute curvature at endpoint. -fn compute_k(bias: f64) -> f64 { - if bias <= 1.0 { - bias * 2.0 - } else if bias < 1.0007 { - let a = bias - 1.0; - // A few terms of the Taylors series expansion of the formula below. - 2.0 + 4.0 / 3.0 * a + 11.0 / 9.0 * a * a - } else { - let a = (bias - 1.0).min(MAX_A); - // Reciprocal of integral - let sr = (a * a) / (1.0 / (1.0 - a) + (1.0 - a).ln() - 1.0); - sr / (1.0 - a).powi(2) +impl ParamCurveFit for Hyperbezier { + fn sample_pt_tangent(&self, t: f64, _sign: f64) -> CurveFitSample { + let (p, tangent) = self.sample_pt_deriv(t); + CurveFitSample { p, tangent } } -} -/// Compute bias given relative endpoint curvature. -/// -/// This is the inverse of compute_k. -pub(crate) fn compute_k_inv(k: f64) -> f64 { - if k <= 2.0 { - k * 0.5 - } else { - // We'll solve this by bisection for now, just for simplicity. - // I'm sure there are much better approximations. - let mut est_lo = 2.0 - 2.0 / k; - let mut est_hi = 2.0 - 1.0 / k; - const N: usize = 20; - for _ in 0..N { - let est = 0.5 * (est_lo + est_hi); - if compute_k(est) > k { - est_hi = est; - } else { - est_lo = est; - } - } - 0.5 * (est_lo + est_hi) + fn sample_pt_deriv(&self, t: f64) -> (Point, Vec2) { + let p = self.eval(t); + let uv = Vec2::from_angle(self.params.theta(t)); + let s = self.scale_rot; + let d = Vec2::new(s.x * uv.x - s.y * uv.y, s.x * uv.y + s.y * uv.x); + (p, d) } -} -#[test] -fn test_k() { - for k in &[0.0, 1.0, 2.0, 2.000001, 3.0, 5.0, 10.0, 20.0] { - let bias = compute_k_inv(*k); - let actual_k = compute_k(bias); - assert!((k - actual_k).abs() < 1e-5); - //println!("{}, {}, {}", k, bias, actual_k); + fn break_cusp(&self, _: std::ops::Range) -> Option { + None } } diff --git a/src/hyperbezier_old.rs b/src/hyperbezier_old.rs new file mode 100644 index 0000000..849eabe --- /dev/null +++ b/src/hyperbezier_old.rs @@ -0,0 +1,349 @@ +//! The math for the hyperbezier curve family. + +use kurbo::common as coeffs; +use kurbo::{Affine, BezPath, PathEl, Point, Vec2}; + +use crate::util; + +/// Parameters for a hyperbezier curve. +/// +/// A hyperbezier is a curve defined by curvature as a function of arclength. +/// It is similar to the Spiro curve in this way, but for some of the parameter +/// space the function is different. +/// +/// The parameter space is four dimensional. It is broken down symmetrically +/// into two parameters that predominantly affect one side of the curve, and +/// the curvature contributions are added: +/// +/// k(s) = k0 * f(bias0, 1 - s) + k1 * f(bias1, s) +/// +/// The "f" function takes a bias parameter, which can also be thought of as +/// tension. This value ranges from about -1 to exactly 2, with 2 representing +/// a cusp (infinitely high tension at the endpoint). For bias values less than +/// 1, it is defined thus: +/// +/// f(bias, s) = s + 6 * (1 - bias) * (s^2 - s^3 - s) +/// +/// For bias values greater than one, it is defined thus: +/// +/// f(bias, s) = c * s / (1 + (bias - 1) * s) ^ 2 +/// +/// Here, c is a normalization term chosen so that the integral of f from s=0 +/// to s=1 is 1. +/// +/// A few observation. If both bias values are 1, then the curve is an Euler +/// spiral. If both bias values are less than 1, then curvature is a cubic +/// polynomial as a function of arclength, so it is a Spiro curve. +#[derive(Copy, Clone, Debug)] +pub struct HyperBezier { + pub k0: f64, + pub bias0: f64, + pub k1: f64, + pub bias1: f64, +} + +/// An intermediate parametrization of the curve family. +/// +/// Here, angles are given relative to the chord, but the bias parameters +/// are the same as for `HyperBezier`. +#[derive(Copy, Clone, Debug)] +pub struct ThetaParams { + pub th0: f64, + pub bias0: f64, + pub th1: f64, + pub bias1: f64, +} + +/// Result of measuring the curve. +/// +/// The `th0` and `th1` values are defined so that if they are have the +/// same sign, the curve is convex, but if they are opposite signs, it is +/// an "s" shape. +#[derive(Copy, Clone, Debug)] +pub struct HyperBezierResult { + /// Tangent angle from the chord to the curve at the start point. + pub th0: f64, + /// Tangent angle from the chord to the curve at the end point. + pub th1: f64, + /// Length of the chord assuming total arclength = 1. + pub chord: f64, + pub k0: f64, + pub k1: f64, +} + +impl HyperBezier { + /// Compute the angle for the given parameter. + /// + /// The argument is an arclength parametrization, ranging from 0 to 1. + /// + /// The returned angle is relative only, in other words there could be an + /// arbitrary rotation of the entire curve. + pub fn compute_theta(&self, s: f64) -> f64 { + self.k1 * integrate_basis(self.bias1, s) - self.k0 * integrate_basis(self.bias0, 1.0 - s) + } + + /// Compute the endpoint tangent angles and the chord length. + pub fn compute(&self) -> HyperBezierResult { + let integral = self.integrate(0.0, 1.0, 24); + let th_chord = integral.atan2(); + let chord = integral.hypot(); + let th0 = th_chord - self.compute_theta(0.0); + let th1 = self.compute_theta(1.0) - th_chord; + let k0 = chord * self.k0 * compute_k(self.bias0); + let k1 = chord * self.k1 * compute_k(self.bias1); + HyperBezierResult { + th0, + th1, + chord, + k0, + k1, + } + } + + fn integrate(&self, t0: f64, t1: f64, order: usize) -> Vec2 { + let c = match order { + 3 => coeffs::GAUSS_LEGENDRE_COEFFS_3, + 5 => coeffs::GAUSS_LEGENDRE_COEFFS_5, + 7 => coeffs::GAUSS_LEGENDRE_COEFFS_7, + 9 => coeffs::GAUSS_LEGENDRE_COEFFS_9, + 11 => coeffs::GAUSS_LEGENDRE_COEFFS_11, + 24 => coeffs::GAUSS_LEGENDRE_COEFFS_24, + _ => panic!("don't have coefficients for {}", order), + }; + let mut result = Vec2::ZERO; + let tm = 0.5 * (t1 + t0); + let dt = 0.5 * (t1 - t0); + for (wi, xi) in c { + let t = tm + dt * xi; + let th = self.compute_theta(t); + result += *wi * Vec2::from_angle(th); + } + dt * result + } + + /// Render to a [`BezPath`]. + pub fn render(&self, n: usize) -> BezPath { + self.render_elements(n).collect() + } + + /// Render to bezier elements. + /// + /// The current algorithm just does a fixed subdivision based on arclength, + /// but should be adaptive in several ways; more subdivision for twistier + /// curves, and also more sophisticated parametrization (important as tension + /// increases). + pub fn render_elements<'a>(&'a self, n: usize) -> impl Iterator + 'a { + let order = 24; + let v = self.integrate(0.0, 1.0, order); + let a = Affine::new([v.x, v.y, -v.y, v.x, 0.0, 0.0]).inverse(); + let step = 1.0 / (n as f64); + fn calc_t(bias: f64) -> f64 { + if bias >= 1.0 { + (2.0 - bias).sqrt() * (1.0 / 3.0) + } else { + // Possibly this should increase for low tension curves, but that's not + // obvious. + 1.0 / 3.0 + } + } + let t1 = calc_t(self.bias0); + let t2 = 1.0 - calc_t(self.bias1); + let mut last_p = Point::ZERO; + let mut last_v = step * t1 * Vec2::from_angle(self.compute_theta(0.0)); + let mut i = 0; + let mut first = Some(PathEl::MoveTo(last_p)); + std::iter::from_fn(move || { + if let Some(first) = first.take() { + return Some(first); + } + i += 1; + if i <= n { + let u = (i as f64) * step; + let um = 1.0 - u; + let t = 3.0 * u * um * (um * t1 + u * t2) + u.powi(3); + let p = self.integrate(0.0, t, order).to_point(); + let p1 = last_p + last_v; + let dt = um * um * t1 + 2.0 * u * um * (t2 - t1) + u * u * (1.0 - t2); + let v = step * dt * Vec2::from_angle(self.compute_theta(t)); + let p2 = p - v; + let next = PathEl::CurveTo(a * p1, a * p2, a * p); + last_v = v; + last_p = p; + Some(next) + } else { + None + } + }) + } + + /// Suggest a number of subdivisions for rendering. + /// + /// This is a bit of a hacky heuristic. + pub fn render_subdivisions(&self) -> usize { + 2 + (self.k0.abs() + self.k1.abs()).floor() as usize + } + + /// Solve for curve params, given theta params. + pub fn solve_for_theta(params: &ThetaParams) -> HyperBezier { + let ThetaParams { + th0, + bias0, + th1, + bias1, + } = *params; + let mut dth = 0.0; + let mut lastxy: Option<(f64, f64)> = None; + const N: usize = 10; + for i in 0..N { + let params = HyperBezier { + k0: th0 + 0.5 * dth, + bias0, + k1: th1 - 0.5 * dth, + bias1, + }; + if i == N - 1 { + return params; + } + let result = params.compute(); + let th_err = util::mod_tau(th0 - th1 - (result.th0 - result.th1)); + if th_err.abs() < 1e-3 { + return params; + } + // Secant method + let nextxy = (dth, th_err); + let delta = if let Some(lastxy) = lastxy { + (nextxy.0 - lastxy.0) / (nextxy.1 - lastxy.1) + } else { + -0.5 + }; + dth -= delta * th_err; + lastxy = Some(nextxy); + } + unreachable!() + } + + /// Solve for curve params, given bezier control points. + /// + /// The points are given relative to p0 at (0, 0) and p3 at + /// (1, 0). + pub fn solve(p1: Point, p2: Point) -> HyperBezier { + let (th0, bias0) = Self::params_for_v(p1.to_vec2()); + let (th1, bias1) = Self::params_for_v(Point::new(1.0, 0.0) - p2); + // TODO: signs feel reversed here, but it all works out in the end. + let theta_params = ThetaParams { + th0: -th0, + bias0, + th1, + bias1, + }; + Self::solve_for_theta(&theta_params) + } + + /// Determine params for a control arm. + /// + /// Return values are theta and bias. + pub(crate) fn params_for_v(v: Vec2) -> (f64, f64) { + let th = v.atan2(); + // This formula ensures that bezier parameters approximating + // a circular arc map to a bias of 1.0. + let a = v.hypot() * 1.5 * (th.cos() + 1.0); + let bias = if a < 1.0 { + 2.0 - a * a + } else { + 1.0 + 2.0 * (0.5 * (1.0 - a)).tanh() + }; + (th, bias) + } + + /// Determine control arm position from params. + /// + /// This function should be the inverse of `params_for_v` + pub(crate) fn v_for_params(th: f64, bias: f64) -> Vec2 { + let a = if bias >= 1.0 { + (2.0 - bias).sqrt() + } else { + // Bias may not be bounded from below, we probably want + // to rethink this... + 1.0 - 2.0 * (0.5 * (bias - 1.0)).atanh() + }; + let len = a / (1.5 * (th.cos() + 1.0)); + len * Vec2::from_angle(th) + } +} + +const MAX_A: f64 = 1.0 - 1e-4; + +/// Compute integral of basis function. +/// +/// The integral of the basis function can be represented as a reasonably +/// simple closed-form analytic formula. +/// +/// Note: this is normalized so that f(1) - f(0) = 1. +/// +/// This is oriented for the rightmost control point. +fn integrate_basis(bias: f64, s: f64) -> f64 { + if bias <= 1.0 { + let iy0 = 4.0 * s.powi(3) - 3.0 * s.powi(4); + let iy1 = s.powi(2); + iy0 + bias * (iy1 - iy0) + } else if bias < 1.0002 { + // This is a more numerically robust approximation to the + // exact analytical formula in the next clause. + let b = (bias - 1.0) * (4.0 / 3.0); + (1.0 - b) * s.powi(2) + b * s.powi(3) + } else { + let a = (bias - 1.0).min(MAX_A); + let norm = 1.0 / (1.0 - a) + (1.0 - a).ln() - 1.0; + (1.0 / (1.0 - a * s) + (1.0 - a * s).ln() - 1.0) / norm + } +} + +/// Compute curvature at endpoint. +fn compute_k(bias: f64) -> f64 { + if bias <= 1.0 { + bias * 2.0 + } else if bias < 1.0007 { + let a = bias - 1.0; + // A few terms of the Taylors series expansion of the formula below. + 2.0 + 4.0 / 3.0 * a + 11.0 / 9.0 * a * a + } else { + let a = (bias - 1.0).min(MAX_A); + // Reciprocal of integral + let sr = (a * a) / (1.0 / (1.0 - a) + (1.0 - a).ln() - 1.0); + sr / (1.0 - a).powi(2) + } +} + +/// Compute bias given relative endpoint curvature. +/// +/// This is the inverse of compute_k. +pub(crate) fn compute_k_inv(k: f64) -> f64 { + if k <= 2.0 { + k * 0.5 + } else { + // We'll solve this by bisection for now, just for simplicity. + // I'm sure there are much better approximations. + let mut est_lo = 2.0 - 2.0 / k; + let mut est_hi = 2.0 - 1.0 / k; + const N: usize = 20; + for _ in 0..N { + let est = 0.5 * (est_lo + est_hi); + if compute_k(est) > k { + est_hi = est; + } else { + est_lo = est; + } + } + 0.5 * (est_lo + est_hi) + } +} + +#[test] +fn test_k() { + for k in &[0.0, 1.0, 2.0, 2.000001, 3.0, 5.0, 10.0, 20.0] { + let bias = compute_k_inv(*k); + let actual_k = compute_k(bias); + assert!((k - actual_k).abs() < 1e-5); + //println!("{}, {}, {}", k, bias, actual_k); + } +} diff --git a/src/lib.rs b/src/lib.rs index 62017a2..08f3fd9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -10,11 +10,12 @@ //! [Spiro]: https://github.com/raphlinus/spiro //! [research spline]: https://github.com/raphlinus/spline-research -mod hyperbezier; +mod hyperbezier_old; +pub mod hyperbezier; mod simple_spline; mod spline; mod util; pub use crate::spline::{Element, Segment, Spline, SplineSpec}; -pub use hyperbezier::{HyperBezier, ThetaParams}; +pub use hyperbezier_old::{HyperBezier, ThetaParams}; pub use simple_spline::SimpleSpline; diff --git a/src/simple_spline.rs b/src/simple_spline.rs index ca1f814..81caea5 100644 --- a/src/simple_spline.rs +++ b/src/simple_spline.rs @@ -4,7 +4,7 @@ use std::f64::consts::PI; use kurbo::{Affine, BezPath, Point, Vec2}; -use crate::hyperbezier::{HyperBezier, HyperBezierResult, ThetaParams}; +use crate::hyperbezier_old::{HyperBezier, HyperBezierResult, ThetaParams}; use crate::util; pub struct SimpleSpline { diff --git a/src/spline.rs b/src/spline.rs index 35f4921..76ed392 100644 --- a/src/spline.rs +++ b/src/spline.rs @@ -6,7 +6,7 @@ use kurbo::{Affine, BezPath, PathEl, Point, Vec2}; #[cfg(feature = "serde")] use serde_::{Deserialize, Serialize}; -use crate::hyperbezier::{self, HyperBezier, ThetaParams}; +use crate::hyperbezier_old::{self, HyperBezier, ThetaParams}; use crate::simple_spline; use crate::util; @@ -406,7 +406,7 @@ impl SplineSpec { let prev_seg = &self.segments[self.prev_ix(i) - 1]; let seg = &self.segments[i - 1]; let this_ch = seg.chord().hypot(); - let bias = hyperbezier::compute_k_inv(prev_seg.k1 * this_ch / (seg.hb.k0 * seg.ch)); + let bias = hyperbezier_old::compute_k_inv(prev_seg.k1 * this_ch / (seg.hb.k0 * seg.ch)); let bias = bias.max(MIN_BIAS); let bias = seg.hb.bias0 + scale * (bias - seg.hb.bias0); self.segments[i - 1].hb.bias0 = bias; @@ -417,7 +417,7 @@ impl SplineSpec { let next_seg = &self.segments[self.next_ix(i) - 1]; let seg = &self.segments[i - 1]; let this_ch = seg.chord().hypot(); - let bias = hyperbezier::compute_k_inv(next_seg.k0 * this_ch / (seg.hb.k1 * seg.ch)); + let bias = hyperbezier_old::compute_k_inv(next_seg.k0 * this_ch / (seg.hb.k1 * seg.ch)); let bias = bias.max(MIN_BIAS); let bias = seg.hb.bias1 + scale * (bias - seg.hb.bias1); self.segments[i - 1].hb.bias1 = bias;