From 2ea00610188dce1eba1172a3ded8244570189230 Mon Sep 17 00:00:00 2001 From: Aevyrie Date: Fri, 3 Mar 2023 22:06:42 +0000 Subject: [PATCH] Add generic cubic splines to `bevy_math` (#7683) # Objective - Make cubic splines more flexible and more performant - Remove the existing spline implementation that is generic over many degrees - This is a potential performance footgun and adds type complexity for negligible gain. - Add implementations of: - Bezier splines - Cardinal splines (inc. Catmull-Rom) - B-Splines - Hermite splines https://user-images.githubusercontent.com/2632925/221780519-495d1b20-ab46-45b4-92a3-32c46da66034.mp4 https://user-images.githubusercontent.com/2632925/221780524-2b154016-699f-404f-9c18-02092f589b04.mp4 https://user-images.githubusercontent.com/2632925/221780525-f934f99d-9ad4-4999-bae2-75d675f5644f.mp4 ## Solution - Implements the concept that splines are curve generators (e.g. https://youtu.be/jvPPXbo87ds?t=3488) via the `CubicGenerator` trait. - Common splines are bespoke data types that implement this trait. This gives us flexibility to add custom spline-specific methods on these types, while ultimately all generating a `CubicCurve`. - All splines generate `CubicCurve`s, which are a chain of precomputed polynomial coefficients. This means that all splines have the same evaluation cost, as the calculations for determining position, velocity, and acceleration are all identical. In addition, `CubicCurve`s are simply a list of `CubicSegment`s, which are evaluated from t=0 to t=1. This also means cubic splines of different type can be chained together, as ultimately they all are simply a collection of `CubicSegment`s. - Because easing is an operation on a singe segment of a Bezier curve, we can simply implement easing on `Beziers` that use the `Vec2` type for points. Higher level crates such as `bevy_ui` can wrap this in a more ergonomic interface as needed. ### Performance Measured on a desktop i5 8600K (6-year-old CPU): - easing: 2.7x faster (19ns) - cubic vec2 position sample: 1.5x faster (1.8ns) - cubic vec3 position sample: 1.5x faster (2.6ns) - cubic vec3a position sample: 1.9x faster (1.4ns) On a laptop i7 11800H: - easing: 16ns - cubic vec2 position sample: 1.6ns - cubic vec3 position sample: 2.3ns - cubic vec3a position sample: 1.2ns --- ## Changelog - Added a generic cubic curve trait, and implementation for Cardinal splines (including Catmull-Rom), B-Splines, Beziers, and Hermite Splines. 2D cubic curve segments also implement easing functionality for animation. --- benches/benches/bevy_math/bezier.rs | 127 ++---- crates/bevy_math/src/bezier.rs | 457 ------------------- crates/bevy_math/src/cubic_splines.rs | 602 ++++++++++++++++++++++++++ crates/bevy_math/src/lib.rs | 12 +- 4 files changed, 650 insertions(+), 548 deletions(-) delete mode 100644 crates/bevy_math/src/bezier.rs create mode 100644 crates/bevy_math/src/cubic_splines.rs diff --git a/benches/benches/bevy_math/bezier.rs b/benches/benches/bevy_math/bezier.rs index e925976de22b0..8f9c866d32c50 100644 --- a/benches/benches/bevy_math/bezier.rs +++ b/benches/benches/bevy_math/bezier.rs @@ -1,128 +1,89 @@ use criterion::{black_box, criterion_group, criterion_main, Criterion}; -use bevy_math::*; +use bevy_math::{prelude::*, *}; fn easing(c: &mut Criterion) { - let cubic_bezier = CubicBezierEasing::new(vec2(0.25, 0.1), vec2(0.25, 1.0)); + let cubic_bezier = CubicSegment::new_bezier(vec2(0.25, 0.1), vec2(0.25, 1.0)); c.bench_function("easing_1000", |b| { b.iter(|| { (0..1000).map(|i| i as f32 / 1000.0).for_each(|t| { - cubic_bezier.ease(black_box(t)); + black_box(cubic_bezier.ease(black_box(t))); }) }); }); } -fn fifteen_degree(c: &mut Criterion) { - let bezier = Bezier::::new([ - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - ]); - c.bench_function("fifteen_degree_position", |b| { - b.iter(|| bezier.position(black_box(0.5))); - }); -} - -fn quadratic_2d(c: &mut Criterion) { - let bezier = QuadraticBezier2d::new([[0.0, 0.0], [0.0, 1.0], [1.0, 1.0]]); - c.bench_function("quadratic_position_Vec2", |b| { - b.iter(|| bezier.position(black_box(0.5))); - }); -} - -fn quadratic(c: &mut Criterion) { - let bezier = QuadraticBezier3d::new([[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 1.0]]); - c.bench_function("quadratic_position_Vec3A", |b| { - b.iter(|| bezier.position(black_box(0.5))); - }); -} - -fn quadratic_vec3(c: &mut Criterion) { - let bezier = Bezier::::new([[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 1.0]]); - c.bench_function("quadratic_position_Vec3", |b| { - b.iter(|| bezier.position(black_box(0.5))); - }); -} - fn cubic_2d(c: &mut Criterion) { - let bezier = CubicBezier2d::new([[0.0, 0.0], [0.0, 1.0], [1.0, 0.0], [1.0, 1.0]]); + let bezier = Bezier::new([[ + vec2(0.0, 0.0), + vec2(0.0, 1.0), + vec2(1.0, 0.0), + vec2(1.0, 1.0), + ]]) + .to_curve(); c.bench_function("cubic_position_Vec2", |b| { - b.iter(|| bezier.position(black_box(0.5))); + b.iter(|| black_box(bezier.position(black_box(0.5)))); }); } fn cubic(c: &mut Criterion) { - let bezier = CubicBezier3d::new([ - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - ]); + let bezier = Bezier::new([[ + vec3a(0.0, 0.0, 0.0), + vec3a(0.0, 1.0, 0.0), + vec3a(1.0, 0.0, 0.0), + vec3a(1.0, 1.0, 1.0), + ]]) + .to_curve(); c.bench_function("cubic_position_Vec3A", |b| { - b.iter(|| bezier.position(black_box(0.5))); + b.iter(|| black_box(bezier.position(black_box(0.5)))); }); } fn cubic_vec3(c: &mut Criterion) { - let bezier = Bezier::::new([ - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - ]); + let bezier = Bezier::new([[ + vec3(0.0, 0.0, 0.0), + vec3(0.0, 1.0, 0.0), + vec3(1.0, 0.0, 0.0), + vec3(1.0, 1.0, 1.0), + ]]) + .to_curve(); c.bench_function("cubic_position_Vec3", |b| { - b.iter(|| bezier.position(black_box(0.5))); + b.iter(|| black_box(bezier.position(black_box(0.5)))); }); } fn build_pos_cubic(c: &mut Criterion) { - let bezier = CubicBezier3d::new([ - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - ]); + let bezier = Bezier::new([[ + vec3a(0.0, 0.0, 0.0), + vec3a(0.0, 1.0, 0.0), + vec3a(1.0, 0.0, 0.0), + vec3a(1.0, 1.0, 1.0), + ]]) + .to_curve(); c.bench_function("build_pos_cubic_100_points", |b| { - b.iter(|| bezier.iter_positions(black_box(100)).collect::>()); + b.iter(|| black_box(bezier.iter_positions(black_box(100)).collect::>())); }); } fn build_accel_cubic(c: &mut Criterion) { - let bezier = CubicBezier3d::new([ - [0.0, 0.0, 0.0], - [0.0, 1.0, 0.0], - [1.0, 0.0, 0.0], - [1.0, 1.0, 1.0], - ]); + let bezier = Bezier::new([[ + vec3a(0.0, 0.0, 0.0), + vec3a(0.0, 1.0, 0.0), + vec3a(1.0, 0.0, 0.0), + vec3a(1.0, 1.0, 1.0), + ]]) + .to_curve(); c.bench_function("build_accel_cubic_100_points", |b| { - b.iter(|| bezier.iter_positions(black_box(100)).collect::>()); + b.iter(|| black_box(bezier.iter_positions(black_box(100)).collect::>())); }); } criterion_group!( benches, easing, - fifteen_degree, - quadratic_2d, - quadratic, - quadratic_vec3, cubic_2d, - cubic, cubic_vec3, + cubic, build_pos_cubic, build_accel_cubic, ); diff --git a/crates/bevy_math/src/bezier.rs b/crates/bevy_math/src/bezier.rs deleted file mode 100644 index 520fcebd0b3cf..0000000000000 --- a/crates/bevy_math/src/bezier.rs +++ /dev/null @@ -1,457 +0,0 @@ -use glam::{Vec2, Vec3, Vec3A}; - -use std::{ - fmt::Debug, - iter::Sum, - ops::{Add, Mul, Sub}, -}; - -/// A point in space of any dimension that supports the mathematical operations needed by -/// [`Bezier`]. -pub trait Point: - Mul - + Add - + Sub - + Sum - + Default - + Debug - + Clone - + PartialEq - + Copy -{ -} -impl Point for Vec3 {} // 3D -impl Point for Vec3A {} // 3D -impl Point for Vec2 {} // 2D -impl Point for f32 {} // 1D - -/// A cubic Bezier curve in 2D space -pub type CubicBezier2d = Bezier; -/// A cubic Bezier curve in 3D space -pub type CubicBezier3d = Bezier; -/// A quadratic Bezier curve in 2D space -pub type QuadraticBezier2d = Bezier; -/// A quadratic Bezier curve in 3D space -pub type QuadraticBezier3d = Bezier; - -/// A generic Bezier curve with `N` control points, and dimension defined by `P`. -/// -/// Consider the following type aliases for most common uses: -/// - [`CubicBezier2d`] -/// - [`CubicBezier3d`] -/// - [`QuadraticBezier2d`] -/// - [`QuadraticBezier3d`] -/// -/// The Bezier degree is equal to `N - 1`. For example, a cubic Bezier has 4 control points, and a -/// degree of 3. The time-complexity of evaluating a Bezier increases superlinearly with the number -/// of control points. As such, it is recommended to instead use a chain of quadratic or cubic -/// `Beziers` instead of a high-degree Bezier. -/// -/// ### About Bezier curves -/// -/// `Bezier` curves are parametric implicit functions; all that means is they take a parameter `t`, -/// and output a point in space, like: -/// -/// > B(t) = (x, y, z) -/// -/// So, all that is needed to find a point in space along a Bezier curve is the parameter `t`. -/// Additionally, the values of `t` are straightforward: `t` is 0 at the start of the curve (first -/// control point) and 1 at the end (last control point). -/// -/// ``` -/// # use bevy_math::{Bezier, Vec2, vec2}; -/// let p0 = vec2(0.25, 0.1); -/// let p1 = vec2(0.25, 1.0); -/// let bezier = Bezier::::new([p0, p1]); -/// assert_eq!(bezier.position(0.0), p0); -/// assert_eq!(bezier.position(1.0), p1); -/// ``` -/// -/// ### Plotting -/// -/// To plot a Bezier curve, simply plug in a series of values of `t` from zero to one. The functions -/// to do this are [`Bezier::position`] to sample the curve at a value of `t`, and -/// [`Bezier::iter_positions`] to iterate over the curve with a number of subdivisions. -/// -/// ### Velocity and Acceleration -/// -/// In addition to the position of a point on the Bezier curve, it is also useful to get information -/// about the curvature. Methods are provided to help with this: -/// -/// - [`Bezier::velocity`]: the instantaneous velocity vector with respect to `t`. This is a vector -/// that points in the direction a point is traveling when it is at point `t`. This vector is -/// tangent to the curve. -/// - [`Bezier::acceleration`]: the instantaneous acceleration vector with respect to `t`. This is a -/// vector that points in the direction a point is accelerating towards when it is at point -/// `t`. This vector will point to the inside of turns, the direction the point is being -/// pulled toward to change direction. -#[derive(Clone, Copy, Debug, PartialEq)] -pub struct Bezier(pub [P; N]); - -impl Default for Bezier { - fn default() -> Self { - Bezier([P::default(); N]) - } -} - -impl Bezier { - /// Construct a new Bezier curve. - pub fn new(control_points: [impl Into

; N]) -> Self { - let control_points = control_points.map(|v| v.into()); - Self(control_points) - } - - /// Compute the position of a point along the Bezier curve at the supplied parametric value `t`. - pub fn position(&self, t: f32) -> P { - generic::position(self.0, t) - } - - /// Compute the first derivative B'(t) of this Bezier at `t` with respect to t. This is the - /// instantaneous velocity of a point on the Bezier curve at `t`. - pub fn velocity(&self, t: f32) -> P { - generic::velocity(self.0, t) - } - - /// Compute the second derivative B''(t) of this Bezier at `t` with respect to t. This is the - /// instantaneous acceleration of a point on the Bezier curve at `t`. - pub fn acceleration(&self, t: f32) -> P { - generic::acceleration(self.0, t) - } - - /// A flexible iterator used to sample [`Bezier`] curves with arbitrary functions. - /// - /// This splits the Bezier into `subdivisions` of evenly spaced `t` values across the length of - /// the curve from start (t = 0) to end (t = 1), returning an iterator that evaluates the curve - /// with the supplied `sample_function` at each `t`. - /// - /// Given `subdivisions = 2`, this will split the curve into two lines, or three points, and - /// return an iterator over those three points, one at the start, middle, and end. - #[inline] - pub fn iter_samples( - &self, - subdivisions: usize, - sample_function: fn(&Self, f32) -> P, - ) -> impl Iterator + '_ { - (0..=subdivisions).map(move |i| { - let t = i as f32 / subdivisions as f32; - sample_function(self, t) - }) - } - - /// Iterate over the curve split into `subdivisions`, sampling the position at each step. - pub fn iter_positions(&self, subdivisions: usize) -> impl Iterator + '_ { - self.iter_samples(subdivisions, Self::position) - } - - /// Iterate over the curve split into `subdivisions`, sampling the velocity at each step. - pub fn iter_velocities(&self, subdivisions: usize) -> impl Iterator + '_ { - self.iter_samples(subdivisions, Self::velocity) - } - - /// Iterate over the curve split into `subdivisions`, sampling the acceleration at each step. - pub fn iter_accelerations(&self, subdivisions: usize) -> impl Iterator + '_ { - self.iter_samples(subdivisions, Self::acceleration) - } -} - -impl, P: Point, const N: usize> From<[T; N]> for Bezier { - fn from(control_points: [T; N]) -> Self { - Bezier::new(control_points) - } -} - -/// A 2-dimensional Bezier curve used for easing in animation. -/// -/// A cubic Bezier easing curve has control point `p0` at (0, 0) and `p3` at (1, 1), leaving only -/// `p1` and `p2` as the remaining degrees of freedom. The first and last control points are fixed -/// to ensure the animation begins at 0, and ends at 1. -#[derive(Default, Clone, Copy, Debug, PartialEq)] -pub struct CubicBezierEasing { - /// Control point P1 of the 2D cubic Bezier curve. Controls the start of the animation. - pub p1: Vec2, - /// Control point P2 of the 2D cubic Bezier curve. Controls the end of the animation. - pub p2: Vec2, -} - -impl CubicBezierEasing { - /// Construct a cubic Bezier curve for animation easing, with control points `p1` and `p2`. - /// These correspond to the two free "handles" of the Bezier curve. - /// - /// This is a very common tool for animations that accelerate and decelerate smoothly. For - /// example, the ubiquitous "ease-in-out" is defined as `(0.25, 0.1), (0.25, 1.0)`. - pub fn new(p1: impl Into, p2: impl Into) -> Self { - Self { - p1: p1.into(), - p2: p2.into(), - } - } - - /// Maximum allowable error for iterative Bezier solve - const MAX_ERROR: f32 = 1e-5; - - /// Maximum number of iterations during Bezier solve - const MAX_ITERS: u8 = 8; - - /// Given a `time` within `0..=1`, returns an eased value that follows the cubic Bezier curve - /// instead of a straight line. This eased result may be outside the range `0..=1`, however it - /// will always start at 0 and end at 1: `ease(0) = 0` and `ease(1) = 1`. - /// - /// ``` - /// # use bevy_math::CubicBezierEasing; - /// let cubic_bezier = CubicBezierEasing::new((0.25, 0.1), (0.25, 1.0)); - /// assert_eq!(cubic_bezier.ease(0.0), 0.0); - /// assert_eq!(cubic_bezier.ease(1.0), 1.0); - /// ``` - /// - /// # How cubic Bezier easing works - /// - /// Easing is generally accomplished with the help of "shaping functions". These are curves that - /// start at (0,0) and end at (1,1). The x-axis of this plot is the current `time` of the - /// animation, from 0 to 1. The y-axis is how far along the animation is, also from 0 to 1. You - /// can imagine that if the shaping function is a straight line, there is a 1:1 mapping between - /// the `time` and how far along your animation is. If the `time` = 0.5, the animation is - /// halfway through. This is known as linear interpolation, and results in objects animating - /// with a constant velocity, and no smooth acceleration or deceleration at the start or end. - /// - /// ```text - /// y - /// │ ● - /// │ ⬈ - /// │ ⬈ - /// │ ⬈ - /// │ ⬈ - /// ●─────────── x (time) - /// ``` - /// - /// Using cubic Beziers, we have a curve that starts at (0,0), ends at (1,1), and follows a path - /// determined by the two remaining control points (handles). These handles allow us to define a - /// smooth curve. As `time` (x-axis) progresses, we now follow the curve, and use the `y` value - /// to determine how far along the animation is. - /// - /// ```text - /// y - /// ⬈➔● - /// │ ⬈ - /// │ ↑ - /// │ ↑ - /// │ ⬈ - /// ●➔⬈───────── x (time) - /// ``` - /// - /// To accomplish this, we need to be able to find the position `y` on a Bezier curve, given the - /// `x` value. As discussed in the [`Bezier`] documentation, a Bezier curve is an implicit - /// parametric function like B(t) = (x,y). To find `y`, we first solve for `t` that corresponds - /// to the given `x` (`time`). We use the Newton-Raphson root-finding method to quickly find a - /// value of `t` that matches `x`. Once we have this we can easily plug that `t` into our - /// Bezier's `position` function, to find the `y` component, which is how far along our - /// animation should be. In other words: - /// - /// > Given `time` in `0..=1` - /// - /// > Use Newton's method to find a value of `t` that results in B(t) = (x,y) where `x == time` - /// - /// > Once a solution is found, use the resulting `y` value as the final result - /// - /// # Performance - /// - /// This operation can be used frequently without fear of performance issues. Benchmarks show - /// this operation taking on the order of 50 nanoseconds. - pub fn ease(&self, time: f32) -> f32 { - let x = time.clamp(0.0, 1.0); - let t = self.find_t_given_x(x); - self.evaluate_y_at(t) - } - - /// Compute the x-coordinate of the point along the Bezier curve at `t`. - #[inline] - fn evaluate_x_at(&self, t: f32) -> f32 { - generic::position([0.0, self.p1.x, self.p2.x, 1.0], t) - } - - /// Compute the y-coordinate of the point along the Bezier curve at `t`. - #[inline] - fn evaluate_y_at(&self, t: f32) -> f32 { - generic::position([0.0, self.p1.y, self.p2.y, 1.0], t) - } - - /// Compute the slope of the line at the given parametric value `t` with respect to t. - #[inline] - fn dx_dt(&self, t: f32) -> f32 { - generic::velocity([0.0, self.p1.x, self.p2.x, 1.0], t) - } - - /// Solve for the parametric value `t` that corresponds to the given value of `x` using the - /// Newton-Raphson method. See documentation on [`Self::ease`] for more details. - #[inline] - fn find_t_given_x(&self, x: f32) -> f32 { - // PERFORMANCE NOTE: - // - // I tried pre-solving and caching values along the curve at struct instantiation to give - // the solver a better starting guess. This ended up being slightly slower, possibly due to - // the increased size of the type. Another option would be to store the last `t`, and use - // that, however it's possible this could end up in a bad state where t is very far from the - // naive but generally safe guess of x, e.g. after an animation resets. - // - // Further optimization might not be needed however - benchmarks are showing it takes about - // 50 nanoseconds for an ease operation on my modern laptop, which seems sufficiently fast. - let mut t_guess = x; - for _ in 0..Self::MAX_ITERS { - let x_guess = self.evaluate_x_at(t_guess); - let error = x_guess - x; - if error.abs() <= Self::MAX_ERROR { - break; - } - // Using Newton's method, use the tangent line to estimate a better guess value. - let slope = self.dx_dt(t_guess); - t_guess -= error / slope; - } - t_guess.clamp(0.0, 1.0) - } -} - -impl> From<[P; 2]> for CubicBezierEasing { - fn from(points: [P; 2]) -> Self { - let [p0, p1] = points; - CubicBezierEasing::new(p0, p1) - } -} - -/// Generic implementations for sampling Bezier curves. Consider using the methods on -/// [`Bezier`](crate::Bezier) for more ergonomic use. -pub mod generic { - use super::Point; - - /// Compute the Bernstein basis polynomial `i` of degree `n`, at `t`. - /// - /// For more information, see . - #[inline] - pub fn bernstein_basis(n: usize, i: usize, t: f32) -> f32 { - (1. - t).powi((n - i) as i32) * t.powi(i as i32) - } - - /// Efficiently compute the binomial coefficient of `n` choose `k`. - #[inline] - fn binomial_coeff(n: usize, k: usize) -> usize { - let k = usize::min(k, n - k); - (0..k).fold(1, |val, i| val * (n - i) / (i + 1)) - } - - /// Evaluate the Bezier curve B(t) of degree `N-1` at the parametric value `t`. - #[inline] - pub fn position(control_points: [P; N], t: f32) -> P { - let p = control_points; - let degree = N - 1; - (0..=degree) - .map(|i| p[i] * binomial_coeff(degree, i) as f32 * bernstein_basis(degree, i, t)) - .sum() - } - - /// Compute the first derivative B'(t) of Bezier curve B(t) of degree `N-1` at the given - /// parametric value `t` with respect to t. Note that the first derivative of a Bezier is also a - /// Bezier, of degree `N-2`. - #[inline] - pub fn velocity(control_points: [P; N], t: f32) -> P { - if N <= 1 { - return P::default(); // Zero for numeric types - } - let p = control_points; - let degree = N - 1; - let degree_vel = N - 2; // the velocity Bezier is one degree lower than the position Bezier - (0..=degree_vel) - .map(|i| { - // Point on the velocity Bezier curve: - let p = (p[i + 1] - p[i]) * degree as f32; - p * binomial_coeff(degree_vel, i) as f32 * bernstein_basis(degree_vel, i, t) - }) - .sum() - } - - /// Compute the second derivative B''(t) of Bezier curve B(t) of degree `N-1` at the given - /// parametric value `t` with respect to t. Note that the second derivative of a Bezier is also - /// a Bezier, of degree `N-3`. - #[inline] - pub fn acceleration(control_points: [P; N], t: f32) -> P { - if N <= 2 { - return P::default(); // Zero for numeric types - } - let p = control_points; - let degree = N - 1; - let degree_vel = N - 2; // the velocity Bezier is one degree lower than the position Bezier - let degree_accel = N - 3; // the accel Bezier is one degree lower than the velocity Bezier - (0..degree_vel) - .map(|i| { - // Points on the velocity Bezier curve: - let p0 = (p[i + 1] - p[i]) * degree as f32; - let p1 = (p[i + 2] - p[i + 1]) * degree as f32; - // Point on the acceleration Bezier curve: - let p = (p1 - p0) * (degree_vel) as f32; - p * binomial_coeff(degree_accel, i) as f32 * bernstein_basis(degree_accel, i, t) - }) - .sum() - } -} - -#[cfg(test)] -mod tests { - use glam::Vec2; - - use crate::{CubicBezier2d, CubicBezierEasing}; - - /// How close two floats can be and still be considered equal - const FLOAT_EQ: f32 = 1e-5; - - /// Basic cubic Bezier easing test to verify the shape of the curve. - #[test] - fn easing_simple() { - // A curve similar to ease-in-out, but symmetric - let bezier = CubicBezierEasing::new([1.0, 0.0], [0.0, 1.0]); - assert_eq!(bezier.ease(0.0), 0.0); - assert!(bezier.ease(0.2) < 0.2); // tests curve - assert_eq!(bezier.ease(0.5), 0.5); // true due to symmetry - assert!(bezier.ease(0.8) > 0.8); // tests curve - assert_eq!(bezier.ease(1.0), 1.0); - } - - /// A curve that forms an upside-down "U", that should extend below 0.0. Useful for animations - /// that go beyond the start and end positions, e.g. bouncing. - #[test] - fn easing_overshoot() { - // A curve that forms an upside-down "U", that should extend above 1.0 - let bezier = CubicBezierEasing::new([0.0, 2.0], [1.0, 2.0]); - assert_eq!(bezier.ease(0.0), 0.0); - assert!(bezier.ease(0.5) > 1.5); - assert_eq!(bezier.ease(1.0), 1.0); - } - - /// A curve that forms a "U", that should extend below 0.0. Useful for animations that go beyond - /// the start and end positions, e.g. bouncing. - #[test] - fn easing_undershoot() { - let bezier = CubicBezierEasing::new([0.0, -2.0], [1.0, -2.0]); - assert_eq!(bezier.ease(0.0), 0.0); - assert!(bezier.ease(0.5) < -0.5); - assert_eq!(bezier.ease(1.0), 1.0); - } - - /// Sweep along the full length of a 3D cubic Bezier, and manually check the position. - #[test] - fn cubic() { - const N_SAMPLES: usize = 1000; - let bezier = CubicBezier2d::new([[-1.0, -20.0], [3.0, 2.0], [5.0, 3.0], [9.0, 8.0]]); - assert_eq!(bezier.position(0.0), bezier.0[0]); // 0 == Start - assert_eq!(bezier.position(1.0), bezier.0[3]); // 1 == End - for i in 0..=N_SAMPLES { - let t = i as f32 / N_SAMPLES as f32; // Check along entire length - assert!(bezier.position(t).distance(cubic_manual(t, bezier)) <= FLOAT_EQ); - } - } - - /// Manual, hardcoded function for computing the position along a cubic bezier. - fn cubic_manual(t: f32, bezier: CubicBezier2d) -> Vec2 { - let [p0, p1, p2, p3] = bezier.0; - p0 * (1.0 - t).powi(3) - + 3.0 * p1 * t * (1.0 - t).powi(2) - + 3.0 * p2 * t.powi(2) * (1.0 - t) - + p3 * t.powi(3) - } -} diff --git a/crates/bevy_math/src/cubic_splines.rs b/crates/bevy_math/src/cubic_splines.rs new file mode 100644 index 0000000000000..d6779f0e0ff80 --- /dev/null +++ b/crates/bevy_math/src/cubic_splines.rs @@ -0,0 +1,602 @@ +//! Provides types for building cubic splines for rendering curves and use with animation easing. + +use glam::{Vec2, Vec3, Vec3A}; + +use std::{ + fmt::Debug, + iter::Sum, + ops::{Add, Mul, Sub}, +}; + +/// A point in space of any dimension that supports the math ops needed for cubic spline +/// interpolation. +pub trait Point: + Mul + + Add + + Sub + + Add + + Sum + + Default + + Debug + + Clone + + PartialEq + + Copy +{ +} +impl Point for Vec3 {} +impl Point for Vec3A {} +impl Point for Vec2 {} +impl Point for f32 {} + +/// A spline composed of a series of cubic Bezier curves. +/// +/// Useful for user-drawn curves with local control, or animation easing. See +/// [`CubicSegment::new_bezier`] for use in easing. +/// +/// ### Interpolation +/// The curve only passes through the first and last control point in each set of four points. +/// +/// ### Tangency +/// Manually defined by the two intermediate control points within each set of four points. +/// +/// ### Continuity +/// At minimum C0 continuous, up to C2. Continuity greater than C0 can result in a loss of local +/// control over the spline due to the curvature constraints. +/// +/// ### Usage +/// +/// ``` +/// # use bevy_math::{*, prelude::*}; +/// let points = [[ +/// vec2(-1.0, -20.0), +/// vec2(3.0, 2.0), +/// vec2(5.0, 3.0), +/// vec2(9.0, 8.0), +/// ]]; +/// let bezier = Bezier::new(points).to_curve(); +/// let positions: Vec<_> = bezier.iter_positions(100).collect(); +/// ``` +pub struct Bezier { + control_points: Vec<[P; 4]>, +} + +impl Bezier

{ + /// Create a new Bezier curve from sets of control points. + pub fn new(control_points: impl Into>) -> Self { + Self { + control_points: control_points.into(), + } + } +} +impl CubicGenerator

for Bezier

{ + #[inline] + fn to_curve(&self) -> CubicCurve

{ + let char_matrix = [ + [1., 0., 0., 0.], + [-3., 3., 0., 0.], + [3., -6., 3., 0.], + [-1., 3., -3., 1.], + ]; + + let segments = self + .control_points + .iter() + .map(|p| CubicCurve::coefficients(*p, 1.0, char_matrix)) + .collect(); + + CubicCurve { segments } + } +} + +/// A spline interpolated continuously between the nearest two control points, with the position and +/// velocity of the curve specified at both control points. This curve passes through all control +/// points, with the specified velocity which includes direction and parametric speed. +/// +/// Useful for smooth interpolation when you know the position and velocity at two points in time, +/// such as network prediction. +/// +/// ### Interpolation +/// The curve passes through every control point. +/// +/// ### Tangency +/// Explicitly defined at each control point. +/// +/// ### Continuity +/// At minimum C0 continuous, up to C1. +/// +/// ### Usage +/// +/// ``` +/// # use bevy_math::{*, prelude::*}; +/// let points = [ +/// vec2(-1.0, -20.0), +/// vec2(3.0, 2.0), +/// vec2(5.0, 3.0), +/// vec2(9.0, 8.0), +/// ]; +/// let tangents = [ +/// vec2(0.0, 1.0), +/// vec2(0.0, 1.0), +/// vec2(0.0, 1.0), +/// vec2(0.0, 1.0), +/// ]; +/// let hermite = Hermite::new(points, tangents).to_curve(); +/// let positions: Vec<_> = hermite.iter_positions(100).collect(); +/// ``` +pub struct Hermite { + control_points: Vec<(P, P)>, +} +impl Hermite

{ + /// Create a new Hermite curve from sets of control points. + pub fn new( + control_points: impl IntoIterator, + tangents: impl IntoIterator, + ) -> Self { + Self { + control_points: control_points + .into_iter() + .zip(tangents.into_iter()) + .collect(), + } + } +} +impl CubicGenerator

for Hermite

{ + #[inline] + fn to_curve(&self) -> CubicCurve

{ + let char_matrix = [ + [1., 0., 0., 0.], + [0., 1., 0., 0.], + [-3., -2., 3., -1.], + [2., 1., -2., 1.], + ]; + + let segments = self + .control_points + .windows(2) + .map(|p| { + let (p0, v0, p1, v1) = (p[0].0, p[0].1, p[1].0, p[1].1); + CubicCurve::coefficients([p0, v0, p1, v1], 1.0, char_matrix) + }) + .collect(); + + CubicCurve { segments } + } +} + +/// A spline interpolated continuously across the nearest four control points, with the position of +/// the curve specified at every control point and the tangents computed automatically. +/// +/// **Note** the Catmull-Rom spline is a special case of Cardinal spline where the tension is 0.5. +/// +/// ### Interpolation +/// The curve passes through every control point. +/// +/// ### Tangency +/// Automatically defined at each control point. +/// +/// ### Continuity +/// C1 continuous. +/// +/// ### Usage +/// +/// ``` +/// # use bevy_math::{*, prelude::*}; +/// let points = [ +/// vec2(-1.0, -20.0), +/// vec2(3.0, 2.0), +/// vec2(5.0, 3.0), +/// vec2(9.0, 8.0), +/// ]; +/// let cardinal = CardinalSpline::new(0.3, points).to_curve(); +/// let positions: Vec<_> = cardinal.iter_positions(100).collect(); +/// ``` +pub struct CardinalSpline { + tension: f32, + control_points: Vec

, +} + +impl CardinalSpline

{ + /// Build a new Cardinal spline. + pub fn new(tension: f32, control_points: impl Into>) -> Self { + Self { + tension, + control_points: control_points.into(), + } + } + + /// Build a new Catmull-Rom spline, the special case of a Cardinal spline where tension = 1/2. + pub fn new_catmull_rom(control_points: impl Into>) -> Self { + Self { + tension: 0.5, + control_points: control_points.into(), + } + } +} +impl CubicGenerator

for CardinalSpline

{ + #[inline] + fn to_curve(&self) -> CubicCurve

{ + let s = self.tension; + let char_matrix = [ + [0., 1., 0., 0.], + [-s, 0., s, 0.], + [2. * s, s - 3., 3. - 2. * s, -s], + [-s, 2. - s, s - 2., s], + ]; + + let segments = self + .control_points + .windows(4) + .map(|p| CubicCurve::coefficients([p[0], p[1], p[2], p[3]], 1.0, char_matrix)) + .collect(); + + CubicCurve { segments } + } +} + +/// A spline interpolated continuously across the nearest four control points. The curve does not +/// pass through any of the control points. +/// +/// ### Interpolation +/// The curve does not pass through control points. +/// +/// ### Tangency +/// Automatically computed based on the position of control points. +/// +/// ### Continuity +/// C2 continuous! The acceleration continuity of this spline makes it useful for camera paths. +/// +/// ### Usage +/// +/// ``` +/// # use bevy_math::{*, prelude::*}; +/// let points = [ +/// vec2(-1.0, -20.0), +/// vec2(3.0, 2.0), +/// vec2(5.0, 3.0), +/// vec2(9.0, 8.0), +/// ]; +/// let b_spline = BSpline::new(points).to_curve(); +/// let positions: Vec<_> = b_spline.iter_positions(100).collect(); +/// ``` +pub struct BSpline { + control_points: Vec

, +} +impl BSpline

{ + /// Build a new Cardinal spline. + pub fn new(control_points: impl Into>) -> Self { + Self { + control_points: control_points.into(), + } + } +} +impl CubicGenerator

for BSpline

{ + #[inline] + fn to_curve(&self) -> CubicCurve

{ + let char_matrix = [ + [1., 4., 1., 0.], + [-3., 0., 3., 0.], + [3., -6., 3., 0.], + [-1., 3., -3., 1.], + ]; + + let segments = self + .control_points + .windows(4) + .map(|p| CubicCurve::coefficients([p[0], p[1], p[2], p[3]], 1.0 / 6.0, char_matrix)) + .collect(); + + CubicCurve { segments } + } +} + +/// Implement this on cubic splines that can generate a curve from their spline parameters. +pub trait CubicGenerator { + /// Build a [`CubicCurve`] by computing the interpolation coefficients for each curve segment. + fn to_curve(&self) -> CubicCurve

; +} + +/// A segment of a cubic curve, used to hold precomputed coefficients for fast interpolation. +/// +/// Segments can be chained together to form a longer compound curve. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct CubicSegment { + coeff: [P; 4], +} + +impl CubicSegment

{ + /// Instantaneous position of a point at parametric value `t`. + #[inline] + pub fn position(&self, t: f32) -> P { + let [a, b, c, d] = self.coeff; + a + b * t + c * t.powi(2) + d * t.powi(3) + } + + /// Instantaneous velocity of a point at parametric value `t`. + #[inline] + pub fn velocity(&self, t: f32) -> P { + let [_, b, c, d] = self.coeff; + b + c * 2.0 * t + d * 3.0 * t.powi(2) + } + + /// Instantaneous acceleration of a point at parametric value `t`. + #[inline] + pub fn acceleration(&self, t: f32) -> P { + let [_, _, c, d] = self.coeff; + c * 2.0 + d * 6.0 * t + } +} + +/// The `CubicSegment` can be used as a 2-dimensional easing curve for animation. +/// +/// The x-axis of the curve is time, and the y-axis is the output value. This struct provides +/// methods for extremely fast solves for y given x. +impl CubicSegment { + /// Construct a cubic Bezier curve for animation easing, with control points `p1` and `p2`. A + /// cubic Bezier easing curve has control point `p0` at (0, 0) and `p3` at (1, 1), leaving only + /// `p1` and `p2` as the remaining degrees of freedom. The first and last control points are + /// fixed to ensure the animation begins at 0, and ends at 1. + /// + /// This is a very common tool for UI animations that accelerate and decelerate smoothly. For + /// example, the ubiquitous "ease-in-out" is defined as `(0.25, 0.1), (0.25, 1.0)`. + pub fn new_bezier(p1: impl Into, p2: impl Into) -> Self { + let (p0, p3) = (Vec2::ZERO, Vec2::ONE); + let bezier = Bezier::new([[p0, p1.into(), p2.into(), p3]]).to_curve(); + bezier.segments[0].clone() + } + + /// Maximum allowable error for iterative Bezier solve + const MAX_ERROR: f32 = 1e-5; + + /// Maximum number of iterations during Bezier solve + const MAX_ITERS: u8 = 8; + + /// Given a `time` within `0..=1`, returns an eased value that follows the cubic curve instead + /// of a straight line. This eased result may be outside the range `0..=1`, however it will + /// always start at 0 and end at 1: `ease(0) = 0` and `ease(1) = 1`. + /// + /// ``` + /// # use bevy_math::prelude::*; + /// let cubic_bezier = CubicSegment::new_bezier((0.25, 0.1), (0.25, 1.0)); + /// assert_eq!(cubic_bezier.ease(0.0), 0.0); + /// assert_eq!(cubic_bezier.ease(1.0), 1.0); + /// ``` + /// + /// # How cubic easing works + /// + /// Easing is generally accomplished with the help of "shaping functions". These are curves that + /// start at (0,0) and end at (1,1). The x-axis of this plot is the current `time` of the + /// animation, from 0 to 1. The y-axis is how far along the animation is, also from 0 to 1. You + /// can imagine that if the shaping function is a straight line, there is a 1:1 mapping between + /// the `time` and how far along your animation is. If the `time` = 0.5, the animation is + /// halfway through. This is known as linear interpolation, and results in objects animating + /// with a constant velocity, and no smooth acceleration or deceleration at the start or end. + /// + /// ```text + /// y + /// │ ● + /// │ ⬈ + /// │ ⬈ + /// │ ⬈ + /// │ ⬈ + /// ●─────────── x (time) + /// ``` + /// + /// Using cubic Beziers, we have a curve that starts at (0,0), ends at (1,1), and follows a path + /// determined by the two remaining control points (handles). These handles allow us to define a + /// smooth curve. As `time` (x-axis) progresses, we now follow the curve, and use the `y` value + /// to determine how far along the animation is. + /// + /// ```text + /// y + /// ⬈➔● + /// │ ⬈ + /// │ ↑ + /// │ ↑ + /// │ ⬈ + /// ●➔⬈───────── x (time) + /// ``` + /// + /// To accomplish this, we need to be able to find the position `y` on a curve, given the `x` + /// value. Cubic curves are implicit parametric functions like B(t) = (x,y). To find `y`, we + /// first solve for `t` that corresponds to the given `x` (`time`). We use the Newton-Raphson + /// root-finding method to quickly find a value of `t` that is very near the desired value of + /// `x`. Once we have this we can easily plug that `t` into our curve's `position` function, to + /// find the `y` component, which is how far along our animation should be. In other words: + /// + /// > Given `time` in `0..=1` + /// + /// > Use Newton's method to find a value of `t` that results in B(t) = (x,y) where `x == time` + /// + /// > Once a solution is found, use the resulting `y` value as the final result + #[inline] + pub fn ease(&self, time: f32) -> f32 { + let x = time.clamp(0.0, 1.0); + self.find_y_given_x(x) + } + + /// Find the `y` value of the curve at the given `x` value using the Newton-Raphson method. + #[inline] + fn find_y_given_x(&self, x: f32) -> f32 { + let mut t_guess = x; + let mut pos_guess = Vec2::ZERO; + for _ in 0..Self::MAX_ITERS { + pos_guess = self.position(t_guess); + let error = pos_guess.x - x; + if error.abs() <= Self::MAX_ERROR { + break; + } + // Using Newton's method, use the tangent line to estimate a better guess value. + let slope = self.velocity(t_guess).x; // dx/dt + t_guess -= error / slope; + } + pos_guess.y + } +} + +/// A collection of [`CubicSegment`]s chained into a curve. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct CubicCurve { + segments: Vec>, +} + +impl CubicCurve

{ + /// Compute the position of a point on the cubic curve at the parametric value `t`. + /// + /// Note that `t` varies from `0..=(n_points - 3)`. + #[inline] + pub fn position(&self, t: f32) -> P { + let (segment, t) = self.segment(t); + segment.position(t) + } + + /// Compute the first derivative with respect to t at `t`. This is the instantaneous velocity of + /// a point on the cubic curve at `t`. + /// + /// Note that `t` varies from `0..=(n_points - 3)`. + #[inline] + pub fn velocity(&self, t: f32) -> P { + let (segment, t) = self.segment(t); + segment.velocity(t) + } + + /// Compute the second derivative with respect to t at `t`. This is the instantaneous + /// acceleration of a point on the cubic curve at `t`. + /// + /// Note that `t` varies from `0..=(n_points - 3)`. + #[inline] + pub fn acceleration(&self, t: f32) -> P { + let (segment, t) = self.segment(t); + segment.acceleration(t) + } + + /// A flexible iterator used to sample curves with arbitrary functions. + /// + /// This splits the curve into `subdivisions` of evenly spaced `t` values across the + /// length of the curve from start (t = 0) to end (t = 1), returning an iterator that evaluates + /// the curve with the supplied `sample_function` at each `t`. + /// + /// Given `subdivisions = 2`, this will split the curve into two lines, or three points, and + /// return an iterator over those three points, one at the start, middle, and end. + #[inline] + pub fn iter_samples( + &self, + subdivisions: usize, + sample_function: fn(&Self, f32) -> P, + ) -> impl Iterator + '_ { + (0..subdivisions).map(move |i| { + let segments = self.segments.len() as f32; + let t = i as f32 / subdivisions as f32 * segments; + sample_function(self, t) + }) + } + + /// Iterate over the curve split into `subdivisions`, sampling the position at each step. + pub fn iter_positions(&self, subdivisions: usize) -> impl Iterator + '_ { + self.iter_samples(subdivisions, Self::position) + } + + /// Iterate over the curve split into `subdivisions`, sampling the velocity at each step. + pub fn iter_velocities(&self, subdivisions: usize) -> impl Iterator + '_ { + self.iter_samples(subdivisions, Self::velocity) + } + + /// Iterate over the curve split into `subdivisions`, sampling the acceleration at each step. + pub fn iter_accelerations(&self, subdivisions: usize) -> impl Iterator + '_ { + self.iter_samples(subdivisions, Self::acceleration) + } + + /// Returns the [`CubicSegment`] and local `t` value given a spline's global `t` value. + #[inline] + fn segment(&self, t: f32) -> (&CubicSegment

, f32) { + if self.segments.len() == 1 { + (&self.segments[0], t) + } else { + let i = (t.floor() as usize).clamp(0, self.segments.len() - 1); + (&self.segments[i], t - i as f32) + } + } + + #[inline] + fn coefficients(p: [P; 4], multiplier: f32, char_matrix: [[f32; 4]; 4]) -> CubicSegment

{ + let [c0, c1, c2, c3] = char_matrix; + // These are the polynomial coefficients, computed by multiplying the characteristic + // matrix by the point matrix. + let mut coeff = [ + p[0] * c0[0] + p[1] * c0[1] + p[2] * c0[2] + p[3] * c0[3], + p[0] * c1[0] + p[1] * c1[1] + p[2] * c1[2] + p[3] * c1[3], + p[0] * c2[0] + p[1] * c2[1] + p[2] * c2[2] + p[3] * c2[3], + p[0] * c3[0] + p[1] * c3[1] + p[2] * c3[2] + p[3] * c3[3], + ]; + coeff.iter_mut().for_each(|c| *c = *c * multiplier); + CubicSegment { coeff } + } +} + +#[cfg(test)] +mod tests { + use glam::{vec2, Vec2}; + + use crate::cubic_splines::{Bezier, CubicGenerator, CubicSegment}; + + /// How close two floats can be and still be considered equal + const FLOAT_EQ: f32 = 1e-5; + + /// Sweep along the full length of a 3D cubic Bezier, and manually check the position. + #[test] + fn cubic() { + const N_SAMPLES: usize = 1000; + let points = [[ + vec2(-1.0, -20.0), + vec2(3.0, 2.0), + vec2(5.0, 3.0), + vec2(9.0, 8.0), + ]]; + let bezier = Bezier::new(points).to_curve(); + for i in 0..=N_SAMPLES { + let t = i as f32 / N_SAMPLES as f32; // Check along entire length + assert!(bezier.position(t).distance(cubic_manual(t, points[0])) <= FLOAT_EQ); + } + } + + /// Manual, hardcoded function for computing the position along a cubic bezier. + fn cubic_manual(t: f32, points: [Vec2; 4]) -> Vec2 { + let p = points; + p[0] * (1.0 - t).powi(3) + + 3.0 * p[1] * t * (1.0 - t).powi(2) + + 3.0 * p[2] * t.powi(2) * (1.0 - t) + + p[3] * t.powi(3) + } + + /// Basic cubic Bezier easing test to verify the shape of the curve. + #[test] + fn easing_simple() { + // A curve similar to ease-in-out, but symmetric + let bezier = CubicSegment::new_bezier([1.0, 0.0], [0.0, 1.0]); + assert_eq!(bezier.ease(0.0), 0.0); + assert!(bezier.ease(0.2) < 0.2); // tests curve + assert_eq!(bezier.ease(0.5), 0.5); // true due to symmetry + assert!(bezier.ease(0.8) > 0.8); // tests curve + assert_eq!(bezier.ease(1.0), 1.0); + } + + /// A curve that forms an upside-down "U", that should extend below 0.0. Useful for animations + /// that go beyond the start and end positions, e.g. bouncing. + #[test] + fn easing_overshoot() { + // A curve that forms an upside-down "U", that should extend above 1.0 + let bezier = CubicSegment::new_bezier([0.0, 2.0], [1.0, 2.0]); + assert_eq!(bezier.ease(0.0), 0.0); + assert!(bezier.ease(0.5) > 1.5); + assert_eq!(bezier.ease(1.0), 1.0); + } + + /// A curve that forms a "U", that should extend below 0.0. Useful for animations that go beyond + /// the start and end positions, e.g. bouncing. + #[test] + fn easing_undershoot() { + let bezier = CubicSegment::new_bezier([0.0, -2.0], [1.0, -2.0]); + assert_eq!(bezier.ease(0.0), 0.0); + assert!(bezier.ease(0.5) < -0.5); + assert_eq!(bezier.ease(1.0), 1.0); + } +} diff --git a/crates/bevy_math/src/lib.rs b/crates/bevy_math/src/lib.rs index e054475f16e4c..f9786ac792d35 100644 --- a/crates/bevy_math/src/lib.rs +++ b/crates/bevy_math/src/lib.rs @@ -6,14 +6,10 @@ #![warn(missing_docs)] -mod bezier; +pub mod cubic_splines; mod ray; mod rect; -pub use bezier::{ - generic as generic_bezier, Bezier, CubicBezier2d, CubicBezier3d, CubicBezierEasing, - QuadraticBezier2d, QuadraticBezier3d, -}; pub use ray::Ray; pub use rect::Rect; @@ -21,9 +17,9 @@ pub use rect::Rect; pub mod prelude { #[doc(hidden)] pub use crate::{ - BVec2, BVec3, BVec4, Bezier, CubicBezier2d, CubicBezier3d, CubicBezierEasing, EulerRot, - IVec2, IVec3, IVec4, Mat2, Mat3, Mat4, QuadraticBezier2d, QuadraticBezier3d, Quat, Ray, - Rect, UVec2, UVec3, UVec4, Vec2, Vec3, Vec4, + cubic_splines::{BSpline, Bezier, CardinalSpline, CubicGenerator, CubicSegment, Hermite}, + BVec2, BVec3, BVec4, EulerRot, IVec2, IVec3, IVec4, Mat2, Mat3, Mat4, Quat, Ray, Rect, + UVec2, UVec3, UVec4, Vec2, Vec3, Vec4, }; }