bevy_math/cubic_splines/mod.rs
1//! Provides types for building cubic splines for rendering curves and use with animation easing.
2
3#[cfg(feature = "curve")]
4mod curve_impls;
5use crate::{
6 ops::{self, FloatPow},
7 Vec2, VectorSpace,
8};
9#[cfg(feature = "bevy_reflect")]
10use bevy_reflect::{std_traits::ReflectDefault, Reflect};
11use thiserror::Error;
12#[cfg(feature = "alloc")]
13use {alloc::vec, alloc::vec::Vec, core::iter::once, itertools::Itertools};
14
15/// A spline composed of a single cubic Bezier curve.
16///
17/// Useful for user-drawn curves with local control, or animation easing. See
18/// [`CubicSegment::new_bezier_easing`] for use in easing.
19///
20/// ### Interpolation
21///
22/// The curve only passes through the first and last control point in each set of four points. The curve
23/// is divided into "segments" by every fourth control point.
24///
25/// ### Tangency
26///
27/// Tangents are manually defined by the two intermediate control points within each set of four points.
28/// You can think of the control points the curve passes through as "anchors", and as the intermediate
29/// control points as the anchors displaced along their tangent vectors
30///
31/// ### Continuity
32///
33/// A Bezier curve is at minimum C0 continuous, meaning it has no holes or jumps. Each curve segment is
34/// C2, meaning the tangent vector changes smoothly between each set of four control points, but this
35/// doesn't hold at the control points between segments. Making the whole curve C1 or C2 requires moving
36/// the intermediate control points to align the tangent vectors between segments, and can result in a
37/// loss of local control.
38///
39/// ### Usage
40///
41/// ```
42/// # use bevy_math::{*, prelude::*};
43/// let points = [[
44/// vec2(-1.0, -20.0),
45/// vec2(3.0, 2.0),
46/// vec2(5.0, 3.0),
47/// vec2(9.0, 8.0),
48/// ]];
49/// let bezier = CubicBezier::new(points).to_curve().unwrap();
50/// let positions: Vec<_> = bezier.iter_positions(100).collect();
51/// ```
52#[derive(Clone, Debug)]
53#[cfg(feature = "alloc")]
54#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
55pub struct CubicBezier<P: VectorSpace> {
56 /// The control points of the Bezier curve.
57 pub control_points: Vec<[P; 4]>,
58}
59
60#[cfg(feature = "alloc")]
61impl<P: VectorSpace> CubicBezier<P> {
62 /// Create a new cubic Bezier curve from sets of control points.
63 pub fn new(control_points: impl IntoIterator<Item = [P; 4]>) -> Self {
64 Self {
65 control_points: control_points.into_iter().collect(),
66 }
67 }
68}
69
70#[cfg(feature = "alloc")]
71impl<P: VectorSpace> CubicGenerator<P> for CubicBezier<P> {
72 type Error = CubicBezierError;
73
74 #[inline]
75 fn to_curve(&self) -> Result<CubicCurve<P>, Self::Error> {
76 let segments = self
77 .control_points
78 .iter()
79 .map(|p| CubicSegment::new_bezier(*p))
80 .collect_vec();
81
82 if segments.is_empty() {
83 Err(CubicBezierError)
84 } else {
85 Ok(CubicCurve { segments })
86 }
87 }
88}
89/// An error returned during cubic curve generation for cubic Bezier curves indicating that a
90/// segment of control points was not present.
91#[derive(Clone, Debug, Error)]
92#[error("Unable to generate cubic curve: at least one set of control points is required")]
93pub struct CubicBezierError;
94
95/// A spline interpolated continuously between the nearest two control points, with the position and
96/// velocity of the curve specified at both control points. This curve passes through all control
97/// points, with the specified velocity which includes direction and parametric speed.
98///
99/// Useful for smooth interpolation when you know the position and velocity at two points in time,
100/// such as network prediction.
101///
102/// ### Interpolation
103///
104/// The curve passes through every control point.
105///
106/// ### Tangency
107///
108/// Tangents are explicitly defined at each control point.
109///
110/// ### Continuity
111///
112/// The curve is at minimum C1 continuous, meaning that it has no holes or jumps and the tangent vector also
113/// has no sudden jumps.
114///
115/// ### Parametrization
116///
117/// The first segment of the curve connects the first two control points, the second connects the second and
118/// third, and so on. This remains true when a cyclic curve is formed with [`to_curve_cyclic`], in which case
119/// the final curve segment connects the last control point to the first.
120///
121/// ### Usage
122///
123/// ```
124/// # use bevy_math::{*, prelude::*};
125/// let points = [
126/// vec2(-1.0, -20.0),
127/// vec2(3.0, 2.0),
128/// vec2(5.0, 3.0),
129/// vec2(9.0, 8.0),
130/// ];
131/// let tangents = [
132/// vec2(0.0, 1.0),
133/// vec2(0.0, 1.0),
134/// vec2(0.0, 1.0),
135/// vec2(0.0, 1.0),
136/// ];
137/// let hermite = CubicHermite::new(points, tangents).to_curve().unwrap();
138/// let positions: Vec<_> = hermite.iter_positions(100).collect();
139/// ```
140///
141/// [`to_curve_cyclic`]: CyclicCubicGenerator::to_curve_cyclic
142#[derive(Clone, Debug)]
143#[cfg(feature = "alloc")]
144#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
145pub struct CubicHermite<P: VectorSpace> {
146 /// The control points of the Hermite curve.
147 pub control_points: Vec<(P, P)>,
148}
149
150#[cfg(feature = "alloc")]
151impl<P: VectorSpace> CubicHermite<P> {
152 /// Create a new Hermite curve from sets of control points.
153 pub fn new(
154 control_points: impl IntoIterator<Item = P>,
155 tangents: impl IntoIterator<Item = P>,
156 ) -> Self {
157 Self {
158 control_points: control_points.into_iter().zip(tangents).collect(),
159 }
160 }
161
162 /// The characteristic matrix for this spline construction.
163 ///
164 /// Each row of this matrix expresses the coefficients of a [`CubicSegment`] as a linear
165 /// combination of `p_i`, `v_i`, `p_{i+1}`, and `v_{i+1}`, where `(p_i, v_i)` and
166 /// `(p_{i+1}, v_{i+1})` are consecutive control points with tangents.
167 #[inline]
168 fn char_matrix(&self) -> [[f32; 4]; 4] {
169 [
170 [1., 0., 0., 0.],
171 [0., 1., 0., 0.],
172 [-3., -2., 3., -1.],
173 [2., 1., -2., 1.],
174 ]
175 }
176}
177
178#[cfg(feature = "alloc")]
179impl<P: VectorSpace> CubicGenerator<P> for CubicHermite<P> {
180 type Error = InsufficientDataError;
181
182 #[inline]
183 fn to_curve(&self) -> Result<CubicCurve<P>, Self::Error> {
184 let segments = self
185 .control_points
186 .windows(2)
187 .map(|p| {
188 let (p0, v0, p1, v1) = (p[0].0, p[0].1, p[1].0, p[1].1);
189 CubicSegment::coefficients([p0, v0, p1, v1], self.char_matrix())
190 })
191 .collect_vec();
192
193 if segments.is_empty() {
194 Err(InsufficientDataError {
195 expected: 2,
196 given: self.control_points.len(),
197 })
198 } else {
199 Ok(CubicCurve { segments })
200 }
201 }
202}
203
204#[cfg(feature = "alloc")]
205impl<P: VectorSpace> CyclicCubicGenerator<P> for CubicHermite<P> {
206 type Error = InsufficientDataError;
207
208 #[inline]
209 fn to_curve_cyclic(&self) -> Result<CubicCurve<P>, Self::Error> {
210 let segments = self
211 .control_points
212 .iter()
213 .circular_tuple_windows()
214 .map(|(&j0, &j1)| {
215 let (p0, v0, p1, v1) = (j0.0, j0.1, j1.0, j1.1);
216 CubicSegment::coefficients([p0, v0, p1, v1], self.char_matrix())
217 })
218 .collect_vec();
219
220 if segments.is_empty() {
221 Err(InsufficientDataError {
222 expected: 2,
223 given: self.control_points.len(),
224 })
225 } else {
226 Ok(CubicCurve { segments })
227 }
228 }
229}
230
231/// A spline interpolated continuously across the nearest four control points, with the position of
232/// the curve specified at every control point and the tangents computed automatically. The associated [`CubicCurve`]
233/// has one segment between each pair of adjacent control points.
234///
235/// **Note** the Catmull-Rom spline is a special case of Cardinal spline where the tension is 0.5.
236///
237/// ### Interpolation
238///
239/// The curve passes through every control point.
240///
241/// ### Tangency
242///
243/// Tangents are automatically computed based on the positions of control points.
244///
245/// ### Continuity
246///
247/// The curve is at minimum C1, meaning that it is continuous (it has no holes or jumps), and its tangent
248/// vector is also well-defined everywhere, without sudden jumps.
249///
250/// ### Parametrization
251///
252/// The first segment of the curve connects the first two control points, the second connects the second and
253/// third, and so on. This remains true when a cyclic curve is formed with [`to_curve_cyclic`], in which case
254/// the final curve segment connects the last control point to the first.
255///
256/// ### Usage
257///
258/// ```
259/// # use bevy_math::{*, prelude::*};
260/// let points = [
261/// vec2(-1.0, -20.0),
262/// vec2(3.0, 2.0),
263/// vec2(5.0, 3.0),
264/// vec2(9.0, 8.0),
265/// ];
266/// let cardinal = CubicCardinalSpline::new(0.3, points).to_curve().unwrap();
267/// let positions: Vec<_> = cardinal.iter_positions(100).collect();
268/// ```
269///
270/// [`to_curve_cyclic`]: CyclicCubicGenerator::to_curve_cyclic
271#[derive(Clone, Debug)]
272#[cfg(feature = "alloc")]
273#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
274pub struct CubicCardinalSpline<P: VectorSpace> {
275 /// Tension
276 pub tension: f32,
277 /// The control points of the Cardinal spline
278 pub control_points: Vec<P>,
279}
280
281#[cfg(feature = "alloc")]
282impl<P: VectorSpace> CubicCardinalSpline<P> {
283 /// Build a new Cardinal spline.
284 pub fn new(tension: f32, control_points: impl IntoIterator<Item = P>) -> Self {
285 Self {
286 tension,
287 control_points: control_points.into_iter().collect(),
288 }
289 }
290
291 /// Build a new Catmull-Rom spline, the special case of a Cardinal spline where tension = 1/2.
292 pub fn new_catmull_rom(control_points: impl IntoIterator<Item = P>) -> Self {
293 Self {
294 tension: 0.5,
295 control_points: control_points.into_iter().collect(),
296 }
297 }
298
299 /// The characteristic matrix for this spline construction.
300 ///
301 /// Each row of this matrix expresses the coefficients of a [`CubicSegment`] as a linear
302 /// combination of four consecutive control points.
303 #[inline]
304 fn char_matrix(&self) -> [[f32; 4]; 4] {
305 let s = self.tension;
306 [
307 [0., 1., 0., 0.],
308 [-s, 0., s, 0.],
309 [2. * s, s - 3., 3. - 2. * s, -s],
310 [-s, 2. - s, s - 2., s],
311 ]
312 }
313}
314
315#[cfg(feature = "alloc")]
316impl<P: VectorSpace> CubicGenerator<P> for CubicCardinalSpline<P> {
317 type Error = InsufficientDataError;
318
319 #[inline]
320 fn to_curve(&self) -> Result<CubicCurve<P>, Self::Error> {
321 let length = self.control_points.len();
322
323 // Early return to avoid accessing an invalid index
324 if length < 2 {
325 return Err(InsufficientDataError {
326 expected: 2,
327 given: self.control_points.len(),
328 });
329 }
330
331 // Extend the list of control points by mirroring the last second-to-last control points on each end;
332 // this allows tangents for the endpoints to be provided, and the overall effect is that the tangent
333 // at an endpoint is proportional to twice the vector between it and its adjacent control point.
334 //
335 // The expression used here is P_{-1} := P_0 - (P_1 - P_0) = 2P_0 - P_1. (Analogously at the other end.)
336 let mirrored_first = self.control_points[0] * 2. - self.control_points[1];
337 let mirrored_last = self.control_points[length - 1] * 2. - self.control_points[length - 2];
338 let extended_control_points = once(&mirrored_first)
339 .chain(self.control_points.iter())
340 .chain(once(&mirrored_last));
341
342 let segments = extended_control_points
343 .tuple_windows()
344 .map(|(&p0, &p1, &p2, &p3)| {
345 CubicSegment::coefficients([p0, p1, p2, p3], self.char_matrix())
346 })
347 .collect_vec();
348
349 Ok(CubicCurve { segments })
350 }
351}
352
353#[cfg(feature = "alloc")]
354impl<P: VectorSpace> CyclicCubicGenerator<P> for CubicCardinalSpline<P> {
355 type Error = InsufficientDataError;
356
357 #[inline]
358 fn to_curve_cyclic(&self) -> Result<CubicCurve<P>, Self::Error> {
359 let len = self.control_points.len();
360
361 if len < 2 {
362 return Err(InsufficientDataError {
363 expected: 2,
364 given: self.control_points.len(),
365 });
366 }
367
368 // This would ordinarily be the last segment, but we pick it out so that we can make it first
369 // in order to get a desirable parametrization where the first segment connects the first two
370 // control points instead of the second and third.
371 let first_segment = {
372 // We take the indices mod `len` in case `len` is very small.
373 let p0 = self.control_points[len - 1];
374 let p1 = self.control_points[0];
375 let p2 = self.control_points[1 % len];
376 let p3 = self.control_points[2 % len];
377 CubicSegment::coefficients([p0, p1, p2, p3], self.char_matrix())
378 };
379
380 let later_segments = self
381 .control_points
382 .iter()
383 .circular_tuple_windows()
384 .map(|(&p0, &p1, &p2, &p3)| {
385 CubicSegment::coefficients([p0, p1, p2, p3], self.char_matrix())
386 })
387 .take(len - 1);
388
389 let mut segments = Vec::with_capacity(len);
390 segments.push(first_segment);
391 segments.extend(later_segments);
392
393 Ok(CubicCurve { segments })
394 }
395}
396
397/// A spline interpolated continuously across the nearest four control points. The curve does not
398/// necessarily pass through any of the control points.
399///
400/// ### Interpolation
401///
402/// The curve does not necessarily pass through its control points.
403///
404/// ### Tangency
405/// Tangents are automatically computed based on the positions of control points.
406///
407/// ### Continuity
408///
409/// The curve is C2 continuous, meaning it has no holes or jumps, the tangent vector changes smoothly along
410/// the entire curve, and the acceleration also varies continuously. The acceleration continuity of this
411/// spline makes it useful for camera paths.
412///
413/// ### Parametrization
414///
415/// Each curve segment is defined by a window of four control points taken in sequence. When [`to_curve_cyclic`]
416/// is used to form a cyclic curve, the three additional segments used to close the curve come last.
417///
418/// ### Usage
419///
420/// ```
421/// # use bevy_math::{*, prelude::*};
422/// let points = [
423/// vec2(-1.0, -20.0),
424/// vec2(3.0, 2.0),
425/// vec2(5.0, 3.0),
426/// vec2(9.0, 8.0),
427/// ];
428/// let b_spline = CubicBSpline::new(points).to_curve().unwrap();
429/// let positions: Vec<_> = b_spline.iter_positions(100).collect();
430/// ```
431///
432/// [`to_curve_cyclic`]: CyclicCubicGenerator::to_curve_cyclic
433#[derive(Clone, Debug)]
434#[cfg(feature = "alloc")]
435#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
436pub struct CubicBSpline<P: VectorSpace> {
437 /// The control points of the spline
438 pub control_points: Vec<P>,
439}
440#[cfg(feature = "alloc")]
441impl<P: VectorSpace> CubicBSpline<P> {
442 /// Build a new B-Spline.
443 pub fn new(control_points: impl IntoIterator<Item = P>) -> Self {
444 Self {
445 control_points: control_points.into_iter().collect(),
446 }
447 }
448
449 /// The characteristic matrix for this spline construction.
450 ///
451 /// Each row of this matrix expresses the coefficients of a [`CubicSegment`] as a linear
452 /// combination of four consecutive control points.
453 #[inline]
454 fn char_matrix(&self) -> [[f32; 4]; 4] {
455 // A derivation for this matrix can be found in "General Matrix Representations for B-splines" by Kaihuai Qin.
456 // <https://xiaoxingchen.github.io/2020/03/02/bspline_in_so3/general_matrix_representation_for_bsplines.pdf>
457 // See section 4.1 and equations 7 and 8.
458 let mut char_matrix = [
459 [1.0, 4.0, 1.0, 0.0],
460 [-3.0, 0.0, 3.0, 0.0],
461 [3.0, -6.0, 3.0, 0.0],
462 [-1.0, 3.0, -3.0, 1.0],
463 ];
464
465 char_matrix
466 .iter_mut()
467 .for_each(|r| r.iter_mut().for_each(|c| *c /= 6.0));
468
469 char_matrix
470 }
471}
472
473#[cfg(feature = "alloc")]
474impl<P: VectorSpace> CubicGenerator<P> for CubicBSpline<P> {
475 type Error = InsufficientDataError;
476
477 #[inline]
478 fn to_curve(&self) -> Result<CubicCurve<P>, Self::Error> {
479 let segments = self
480 .control_points
481 .windows(4)
482 .map(|p| CubicSegment::coefficients([p[0], p[1], p[2], p[3]], self.char_matrix()))
483 .collect_vec();
484
485 if segments.is_empty() {
486 Err(InsufficientDataError {
487 expected: 4,
488 given: self.control_points.len(),
489 })
490 } else {
491 Ok(CubicCurve { segments })
492 }
493 }
494}
495
496#[cfg(feature = "alloc")]
497impl<P: VectorSpace> CyclicCubicGenerator<P> for CubicBSpline<P> {
498 type Error = InsufficientDataError;
499
500 #[inline]
501 fn to_curve_cyclic(&self) -> Result<CubicCurve<P>, Self::Error> {
502 let segments = self
503 .control_points
504 .iter()
505 .circular_tuple_windows()
506 .map(|(&a, &b, &c, &d)| CubicSegment::coefficients([a, b, c, d], self.char_matrix()))
507 .collect_vec();
508
509 // Note that the parametrization is consistent with the one for `to_curve` but with
510 // the extra curve segments all tacked on at the end. This might be slightly counter-intuitive,
511 // since it means the first segment doesn't go "between" the first two control points, but
512 // between the second and third instead.
513
514 if segments.is_empty() {
515 Err(InsufficientDataError {
516 expected: 2,
517 given: self.control_points.len(),
518 })
519 } else {
520 Ok(CubicCurve { segments })
521 }
522 }
523}
524
525/// Error during construction of [`CubicNurbs`]
526#[derive(Clone, Debug, Error)]
527pub enum CubicNurbsError {
528 /// Provided the wrong number of knots.
529 #[error("Wrong number of knots: expected {expected}, provided {provided}")]
530 KnotsNumberMismatch {
531 /// Expected number of knots
532 expected: usize,
533 /// Provided number of knots
534 provided: usize,
535 },
536 /// The provided knots had a descending knot pair. Subsequent knots must
537 /// either increase or stay the same.
538 #[error("Invalid knots: contains descending knot pair")]
539 DescendingKnots,
540 /// The provided knots were all equal. Knots must contain at least one increasing pair.
541 #[error("Invalid knots: all knots are equal")]
542 ConstantKnots,
543 /// Provided a different number of weights and control points.
544 #[error("Incorrect number of weights: expected {expected}, provided {provided}")]
545 WeightsNumberMismatch {
546 /// Expected number of weights
547 expected: usize,
548 /// Provided number of weights
549 provided: usize,
550 },
551 /// The number of control points provided is less than 4.
552 #[error("Not enough control points, at least 4 are required, {provided} were provided")]
553 NotEnoughControlPoints {
554 /// The number of control points provided
555 provided: usize,
556 },
557}
558
559/// Non-uniform Rational B-Splines (NURBS) are a powerful generalization of the [`CubicBSpline`] which can
560/// represent a much more diverse class of curves (like perfect circles and ellipses).
561///
562/// ### Non-uniformity
563///
564/// The 'NU' part of NURBS stands for "Non-Uniform". This has to do with a parameter called 'knots'.
565/// The knots are a non-decreasing sequence of floating point numbers. The first and last three pairs of
566/// knots control the behavior of the curve as it approaches its endpoints. The intermediate pairs
567/// each control the length of one segment of the curve. Multiple repeated knot values are called
568/// "knot multiplicity". Knot multiplicity in the intermediate knots causes a "zero-length" segment,
569/// and can create sharp corners.
570///
571/// ### Rationality
572///
573/// The 'R' part of NURBS stands for "Rational". This has to do with NURBS allowing each control point to
574/// be assigned a weighting, which controls how much it affects the curve compared to the other points.
575///
576/// ### Interpolation
577///
578/// The curve will not pass through the control points except where a knot has multiplicity four.
579///
580/// ### Tangency
581///
582/// Tangents are automatically computed based on the position of control points.
583///
584/// ### Continuity
585///
586/// When there is no knot multiplicity, the curve is C2 continuous, meaning it has no holes or jumps and the
587/// tangent vector changes smoothly along the entire curve length. Like the [`CubicBSpline`], the acceleration
588/// continuity makes it useful for camera paths. Knot multiplicity of 2 in intermediate knots reduces the
589/// continuity to C1, and knot multiplicity of 3 reduces the continuity to C0. The curve is always at least
590/// C0, meaning it has no jumps or holes.
591///
592/// ### Usage
593///
594/// ```
595/// # use bevy_math::{*, prelude::*};
596/// let points = [
597/// vec2(-1.0, -20.0),
598/// vec2(3.0, 2.0),
599/// vec2(5.0, 3.0),
600/// vec2(9.0, 8.0),
601/// ];
602/// let weights = [1.0, 1.0, 2.0, 1.0];
603/// let knots = [0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 5.0];
604/// let nurbs = CubicNurbs::new(points, Some(weights), Some(knots))
605/// .expect("NURBS construction failed!")
606/// .to_curve()
607/// .unwrap();
608/// let positions: Vec<_> = nurbs.iter_positions(100).collect();
609/// ```
610#[derive(Clone, Debug)]
611#[cfg(feature = "alloc")]
612#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
613pub struct CubicNurbs<P: VectorSpace> {
614 /// The control points of the NURBS
615 pub control_points: Vec<P>,
616 /// Weights
617 pub weights: Vec<f32>,
618 /// Knots
619 pub knots: Vec<f32>,
620}
621
622#[cfg(feature = "alloc")]
623impl<P: VectorSpace> CubicNurbs<P> {
624 /// Build a Non-Uniform Rational B-Spline.
625 ///
626 /// If provided, weights must be the same length as the control points. Defaults to equal weights.
627 ///
628 /// If provided, the number of knots must be n + 4 elements, where n is the amount of control
629 /// points. Defaults to open uniform knots: [`Self::open_uniform_knots`]. Knots cannot
630 /// all be equal.
631 ///
632 /// At least 4 points must be provided, otherwise an error will be returned.
633 pub fn new(
634 control_points: impl IntoIterator<Item = P>,
635 weights: Option<impl IntoIterator<Item = f32>>,
636 knots: Option<impl IntoIterator<Item = f32>>,
637 ) -> Result<Self, CubicNurbsError> {
638 let mut control_points: Vec<P> = control_points.into_iter().collect();
639 let control_points_len = control_points.len();
640
641 if control_points_len < 4 {
642 return Err(CubicNurbsError::NotEnoughControlPoints {
643 provided: control_points_len,
644 });
645 }
646
647 let weights: Vec<f32> = weights
648 .map(|ws| ws.into_iter().collect())
649 .unwrap_or_else(|| vec![1.0; control_points_len]);
650
651 let mut knots: Vec<f32> = knots.map(|ks| ks.into_iter().collect()).unwrap_or_else(|| {
652 Self::open_uniform_knots(control_points_len)
653 .expect("The amount of control points was checked")
654 });
655
656 let expected_knots_len = Self::knots_len(control_points_len);
657
658 // Check the number of knots is correct
659 if knots.len() != expected_knots_len {
660 return Err(CubicNurbsError::KnotsNumberMismatch {
661 expected: expected_knots_len,
662 provided: knots.len(),
663 });
664 }
665
666 // Ensure the knots are non-descending (previous element is less than or equal
667 // to the next)
668 if knots.windows(2).any(|win| win[0] > win[1]) {
669 return Err(CubicNurbsError::DescendingKnots);
670 }
671
672 // Ensure the knots are non-constant
673 if knots.windows(2).all(|win| win[0] == win[1]) {
674 return Err(CubicNurbsError::ConstantKnots);
675 }
676
677 // Check that the number of weights equals the number of control points
678 if weights.len() != control_points_len {
679 return Err(CubicNurbsError::WeightsNumberMismatch {
680 expected: control_points_len,
681 provided: weights.len(),
682 });
683 }
684
685 // To align the evaluation behavior of nurbs with the other splines,
686 // make the intervals between knots form an exact cover of [0, N], where N is
687 // the number of segments of the final curve.
688 let curve_length = (control_points.len() - 3) as f32;
689 let min = *knots.first().unwrap();
690 let max = *knots.last().unwrap();
691 let knot_delta = max - min;
692 knots = knots
693 .into_iter()
694 .map(|k| k - min)
695 .map(|k| k * curve_length / knot_delta)
696 .collect();
697
698 control_points
699 .iter_mut()
700 .zip(weights.iter())
701 .for_each(|(p, w)| *p = *p * *w);
702
703 Ok(Self {
704 control_points,
705 weights,
706 knots,
707 })
708 }
709
710 /// Generates uniform knots that will generate the same curve as [`CubicBSpline`].
711 ///
712 /// "Uniform" means that the difference between two subsequent knots is the same.
713 ///
714 /// Will return `None` if there are less than 4 control points.
715 pub fn uniform_knots(control_points: usize) -> Option<Vec<f32>> {
716 if control_points < 4 {
717 return None;
718 }
719 Some(
720 (0..Self::knots_len(control_points))
721 .map(|v| v as f32)
722 .collect(),
723 )
724 }
725
726 /// Generates open uniform knots, which makes the ends of the curve pass through the
727 /// start and end points.
728 ///
729 /// The start and end knots have multiplicity 4, and intermediate knots have multiplicity 0 and
730 /// difference of 1.
731 ///
732 /// Will return `None` if there are less than 4 control points.
733 pub fn open_uniform_knots(control_points: usize) -> Option<Vec<f32>> {
734 if control_points < 4 {
735 return None;
736 }
737 let last_knots_value = control_points - 3;
738 Some(
739 core::iter::repeat_n(0.0, 4)
740 .chain((1..last_knots_value).map(|v| v as f32))
741 .chain(core::iter::repeat_n(last_knots_value as f32, 4))
742 .collect(),
743 )
744 }
745
746 #[inline(always)]
747 const fn knots_len(control_points_len: usize) -> usize {
748 control_points_len + 4
749 }
750
751 /// Generates a non-uniform B-spline characteristic matrix from a sequence of six knots. Each six
752 /// knots describe the relationship between four successive control points. For padding reasons,
753 /// this takes a vector of 8 knots, but only six are actually used.
754 fn generate_matrix(knots: &[f32; 8]) -> [[f32; 4]; 4] {
755 // A derivation for this matrix can be found in "General Matrix Representations for B-splines" by Kaihuai Qin.
756 // <https://xiaoxingchen.github.io/2020/03/02/bspline_in_so3/general_matrix_representation_for_bsplines.pdf>
757 // See section 3.1.
758
759 let t = knots;
760 // In the notation of the paper:
761 // t[1] := t_i-2
762 // t[2] := t_i-1
763 // t[3] := t_i (the lower extent of the current knot span)
764 // t[4] := t_i+1 (the upper extent of the current knot span)
765 // t[5] := t_i+2
766 // t[6] := t_i+3
767
768 let m00 = (t[4] - t[3]).squared() / ((t[4] - t[2]) * (t[4] - t[1]));
769 let m02 = (t[3] - t[2]).squared() / ((t[5] - t[2]) * (t[4] - t[2]));
770 let m12 = (3.0 * (t[4] - t[3]) * (t[3] - t[2])) / ((t[5] - t[2]) * (t[4] - t[2]));
771 let m22 = 3.0 * (t[4] - t[3]).squared() / ((t[5] - t[2]) * (t[4] - t[2]));
772 let m33 = (t[4] - t[3]).squared() / ((t[6] - t[3]) * (t[5] - t[3]));
773 let m32 = -m22 / 3.0 - m33 - (t[4] - t[3]).squared() / ((t[5] - t[3]) * (t[5] - t[2]));
774 [
775 [m00, 1.0 - m00 - m02, m02, 0.0],
776 [-3.0 * m00, 3.0 * m00 - m12, m12, 0.0],
777 [3.0 * m00, -3.0 * m00 - m22, m22, 0.0],
778 [-m00, m00 - m32 - m33, m32, m33],
779 ]
780 }
781}
782
783#[cfg(feature = "alloc")]
784impl<P: VectorSpace> RationalGenerator<P> for CubicNurbs<P> {
785 type Error = InsufficientDataError;
786
787 #[inline]
788 fn to_curve(&self) -> Result<RationalCurve<P>, Self::Error> {
789 let segments = self
790 .control_points
791 .windows(4)
792 .zip(self.weights.windows(4))
793 .zip(self.knots.windows(8))
794 .filter(|(_, knots)| knots[4] - knots[3] > 0.0)
795 .map(|((points, weights), knots)| {
796 // This is curve segment i. It uses control points P_i, P_i+2, P_i+2 and P_i+3,
797 // It is associated with knot span i+3 (which is the interval between knots i+3
798 // and i+4) and its characteristic matrix uses knots i+1 through i+6 (because
799 // those define the two knot spans on either side).
800 let span = knots[4] - knots[3];
801 let coefficient_knots = knots.try_into().expect("Knot windows are of length 6");
802 let matrix = Self::generate_matrix(coefficient_knots);
803 RationalSegment::coefficients(
804 points.try_into().expect("Point windows are of length 4"),
805 weights.try_into().expect("Weight windows are of length 4"),
806 span,
807 matrix,
808 )
809 })
810 .collect_vec();
811 if segments.is_empty() {
812 Err(InsufficientDataError {
813 expected: 4,
814 given: self.control_points.len(),
815 })
816 } else {
817 Ok(RationalCurve { segments })
818 }
819 }
820}
821
822/// A spline interpolated linearly between the nearest 2 points.
823///
824/// ### Interpolation
825///
826/// The curve passes through every control point.
827///
828/// ### Tangency
829///
830/// The curve is not generally differentiable at control points.
831///
832/// ### Continuity
833///
834/// The curve is C0 continuous, meaning it has no holes or jumps.
835///
836/// ### Parametrization
837///
838/// Each curve segment connects two adjacent control points in sequence. When a cyclic curve is
839/// formed with [`to_curve_cyclic`], the final segment connects the last control point with the first.
840///
841/// [`to_curve_cyclic`]: CyclicCubicGenerator::to_curve_cyclic
842#[derive(Clone, Debug)]
843#[cfg(feature = "alloc")]
844#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
845pub struct LinearSpline<P: VectorSpace> {
846 /// The control points of the linear spline.
847 pub points: Vec<P>,
848}
849
850#[cfg(feature = "alloc")]
851impl<P: VectorSpace> LinearSpline<P> {
852 /// Create a new linear spline from a list of points to be interpolated.
853 pub fn new(points: impl IntoIterator<Item = P>) -> Self {
854 Self {
855 points: points.into_iter().collect(),
856 }
857 }
858}
859
860#[cfg(feature = "alloc")]
861impl<P: VectorSpace> CubicGenerator<P> for LinearSpline<P> {
862 type Error = InsufficientDataError;
863
864 #[inline]
865 fn to_curve(&self) -> Result<CubicCurve<P>, Self::Error> {
866 let segments = self
867 .points
868 .windows(2)
869 .map(|points| {
870 let a = points[0];
871 let b = points[1];
872 CubicSegment {
873 coeff: [a, b - a, P::default(), P::default()],
874 }
875 })
876 .collect_vec();
877
878 if segments.is_empty() {
879 Err(InsufficientDataError {
880 expected: 2,
881 given: self.points.len(),
882 })
883 } else {
884 Ok(CubicCurve { segments })
885 }
886 }
887}
888
889#[cfg(feature = "alloc")]
890impl<P: VectorSpace> CyclicCubicGenerator<P> for LinearSpline<P> {
891 type Error = InsufficientDataError;
892
893 #[inline]
894 fn to_curve_cyclic(&self) -> Result<CubicCurve<P>, Self::Error> {
895 let segments = self
896 .points
897 .iter()
898 .circular_tuple_windows()
899 .map(|(&a, &b)| CubicSegment {
900 coeff: [a, b - a, P::default(), P::default()],
901 })
902 .collect_vec();
903
904 if segments.is_empty() {
905 Err(InsufficientDataError {
906 expected: 2,
907 given: self.points.len(),
908 })
909 } else {
910 Ok(CubicCurve { segments })
911 }
912 }
913}
914
915/// An error indicating that a spline construction didn't have enough control points to generate a curve.
916#[derive(Clone, Debug, Error)]
917#[error("Not enough data to build curve: needed at least {expected} control points but was only given {given}")]
918pub struct InsufficientDataError {
919 expected: usize,
920 given: usize,
921}
922
923/// Implement this on cubic splines that can generate a cubic curve from their spline parameters.
924#[cfg(feature = "alloc")]
925pub trait CubicGenerator<P: VectorSpace> {
926 /// An error type indicating why construction might fail.
927 type Error;
928
929 /// Build a [`CubicCurve`] by computing the interpolation coefficients for each curve segment.
930 fn to_curve(&self) -> Result<CubicCurve<P>, Self::Error>;
931}
932
933/// Implement this on cubic splines that can generate a cyclic cubic curve from their spline parameters.
934///
935/// This makes sense only when the control data can be interpreted cyclically.
936#[cfg(feature = "alloc")]
937pub trait CyclicCubicGenerator<P: VectorSpace> {
938 /// An error type indicating why construction might fail.
939 type Error;
940
941 /// Build a cyclic [`CubicCurve`] by computing the interpolation coefficients for each curve segment,
942 /// treating the control data as cyclic so that the result is a closed curve.
943 fn to_curve_cyclic(&self) -> Result<CubicCurve<P>, Self::Error>;
944}
945
946/// A segment of a cubic curve, used to hold precomputed coefficients for fast interpolation.
947/// It is a [`Curve`] with domain `[0, 1]`.
948///
949/// Segments can be chained together to form a longer [compound curve].
950///
951/// [compound curve]: CubicCurve
952/// [`Curve`]: crate::curve::Curve
953#[derive(Copy, Clone, Debug, Default, PartialEq)]
954#[cfg_attr(feature = "serialize", derive(serde::Serialize, serde::Deserialize))]
955#[cfg_attr(
956 feature = "bevy_reflect",
957 derive(Reflect),
958 reflect(Debug, Default, Clone)
959)]
960pub struct CubicSegment<P: VectorSpace> {
961 /// Polynomial coefficients for the segment.
962 pub coeff: [P; 4],
963}
964
965impl<P: VectorSpace> CubicSegment<P> {
966 /// Instantaneous position of a point at parametric value `t`.
967 #[inline]
968 pub fn position(&self, t: f32) -> P {
969 let [a, b, c, d] = self.coeff;
970 // Evaluate `a + bt + ct^2 + dt^3`, avoiding exponentiation
971 a + (b + (c + d * t) * t) * t
972 }
973
974 /// Instantaneous velocity of a point at parametric value `t`.
975 #[inline]
976 pub fn velocity(&self, t: f32) -> P {
977 let [_, b, c, d] = self.coeff;
978 // Evaluate the derivative, which is `b + 2ct + 3dt^2`, avoiding exponentiation
979 b + (c * 2.0 + d * 3.0 * t) * t
980 }
981
982 /// Instantaneous acceleration of a point at parametric value `t`.
983 #[inline]
984 pub fn acceleration(&self, t: f32) -> P {
985 let [_, _, c, d] = self.coeff;
986 // Evaluate the second derivative, which is `2c + 6dt`
987 c * 2.0 + d * 6.0 * t
988 }
989
990 /// Creates a cubic segment from four points, representing a Bezier curve.
991 pub fn new_bezier(points: [P; 4]) -> Self {
992 // A derivation for this matrix can be found in "General Matrix Representations for B-splines" by Kaihuai Qin.
993 // <https://xiaoxingchen.github.io/2020/03/02/bspline_in_so3/general_matrix_representation_for_bsplines.pdf>
994 // See section 4.2 and equation 11.
995 let char_matrix = [
996 [1., 0., 0., 0.],
997 [-3., 3., 0., 0.],
998 [3., -6., 3., 0.],
999 [-1., 3., -3., 1.],
1000 ];
1001 Self::coefficients(points, char_matrix)
1002 }
1003
1004 /// Calculate polynomial coefficients for the cubic curve using a characteristic matrix.
1005 #[inline]
1006 fn coefficients(p: [P; 4], char_matrix: [[f32; 4]; 4]) -> Self {
1007 let [c0, c1, c2, c3] = char_matrix;
1008 // These are the polynomial coefficients, computed by multiplying the characteristic
1009 // matrix by the point matrix.
1010 let coeff = [
1011 p[0] * c0[0] + p[1] * c0[1] + p[2] * c0[2] + p[3] * c0[3],
1012 p[0] * c1[0] + p[1] * c1[1] + p[2] * c1[2] + p[3] * c1[3],
1013 p[0] * c2[0] + p[1] * c2[1] + p[2] * c2[2] + p[3] * c2[3],
1014 p[0] * c3[0] + p[1] * c3[1] + p[2] * c3[2] + p[3] * c3[3],
1015 ];
1016 Self { coeff }
1017 }
1018
1019 /// A flexible iterator used to sample curves with arbitrary functions.
1020 ///
1021 /// This splits the curve into `subdivisions` of evenly spaced `t` values across the
1022 /// length of the curve from start (t = 0) to end (t = n), where `n = self.segment_count()`,
1023 /// returning an iterator evaluating the curve with the supplied `sample_function` at each `t`.
1024 ///
1025 /// For `subdivisions = 2`, this will split the curve into two lines, or three points, and
1026 /// return an iterator with 3 items, the three points, one at the start, middle, and end.
1027 #[inline]
1028 pub fn iter_samples<'a, 'b: 'a>(
1029 &'b self,
1030 subdivisions: usize,
1031 mut sample_function: impl FnMut(&Self, f32) -> P + 'a,
1032 ) -> impl Iterator<Item = P> + 'a {
1033 self.iter_uniformly(subdivisions)
1034 .map(move |t| sample_function(self, t))
1035 }
1036
1037 /// An iterator that returns values of `t` uniformly spaced over `0..=subdivisions`.
1038 #[inline]
1039 fn iter_uniformly(&self, subdivisions: usize) -> impl Iterator<Item = f32> {
1040 let step = 1.0 / subdivisions as f32;
1041 (0..=subdivisions).map(move |i| i as f32 * step)
1042 }
1043
1044 /// Iterate over the curve split into `subdivisions`, sampling the position at each step.
1045 pub fn iter_positions(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1046 self.iter_samples(subdivisions, Self::position)
1047 }
1048
1049 /// Iterate over the curve split into `subdivisions`, sampling the velocity at each step.
1050 pub fn iter_velocities(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1051 self.iter_samples(subdivisions, Self::velocity)
1052 }
1053
1054 /// Iterate over the curve split into `subdivisions`, sampling the acceleration at each step.
1055 pub fn iter_accelerations(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1056 self.iter_samples(subdivisions, Self::acceleration)
1057 }
1058}
1059
1060/// The `CubicSegment<Vec2>` can be used as a 2-dimensional easing curve for animation.
1061///
1062/// The x-axis of the curve is time, and the y-axis is the output value. This struct provides
1063/// methods for extremely fast solves for y given x.
1064impl CubicSegment<Vec2> {
1065 /// Construct a cubic Bezier curve for animation easing, with control points `p1` and `p2`. A
1066 /// cubic Bezier easing curve has control point `p0` at (0, 0) and `p3` at (1, 1), leaving only
1067 /// `p1` and `p2` as the remaining degrees of freedom. The first and last control points are
1068 /// fixed to ensure the animation begins at 0, and ends at 1.
1069 ///
1070 /// This is a very common tool for UI animations that accelerate and decelerate smoothly. For
1071 /// example, the ubiquitous "ease-in-out" is defined as `(0.25, 0.1), (0.25, 1.0)`.
1072 #[cfg(feature = "alloc")]
1073 pub fn new_bezier_easing(p1: impl Into<Vec2>, p2: impl Into<Vec2>) -> Self {
1074 let (p0, p3) = (Vec2::ZERO, Vec2::ONE);
1075 Self::new_bezier([p0, p1.into(), p2.into(), p3])
1076 }
1077
1078 /// Maximum allowable error for iterative Bezier solve
1079 const MAX_ERROR: f32 = 1e-5;
1080
1081 /// Maximum number of iterations during Bezier solve
1082 const MAX_ITERS: u8 = 8;
1083
1084 /// Given a `time` within `0..=1`, returns an eased value that follows the cubic curve instead
1085 /// of a straight line. This eased result may be outside the range `0..=1`, however it will
1086 /// always start at 0 and end at 1: `ease(0) = 0` and `ease(1) = 1`.
1087 ///
1088 /// ```
1089 /// # use bevy_math::prelude::*;
1090 /// # #[cfg(feature = "alloc")]
1091 /// # {
1092 /// let cubic_bezier = CubicSegment::new_bezier_easing((0.25, 0.1), (0.25, 1.0));
1093 /// assert_eq!(cubic_bezier.ease(0.0), 0.0);
1094 /// assert_eq!(cubic_bezier.ease(1.0), 1.0);
1095 /// # }
1096 /// ```
1097 ///
1098 /// # How cubic easing works
1099 ///
1100 /// Easing is generally accomplished with the help of "shaping functions". These are curves that
1101 /// start at (0,0) and end at (1,1). The x-axis of this plot is the current `time` of the
1102 /// animation, from 0 to 1. The y-axis is how far along the animation is, also from 0 to 1. You
1103 /// can imagine that if the shaping function is a straight line, there is a 1:1 mapping between
1104 /// the `time` and how far along your animation is. If the `time` = 0.5, the animation is
1105 /// halfway through. This is known as linear interpolation, and results in objects animating
1106 /// with a constant velocity, and no smooth acceleration or deceleration at the start or end.
1107 ///
1108 /// ```text
1109 /// y
1110 /// │ ●
1111 /// │ ⬈
1112 /// │ ⬈
1113 /// │ ⬈
1114 /// │ ⬈
1115 /// ●─────────── x (time)
1116 /// ```
1117 ///
1118 /// Using cubic Beziers, we have a curve that starts at (0,0), ends at (1,1), and follows a path
1119 /// determined by the two remaining control points (handles). These handles allow us to define a
1120 /// smooth curve. As `time` (x-axis) progresses, we now follow the curve, and use the `y` value
1121 /// to determine how far along the animation is.
1122 ///
1123 /// ```text
1124 /// y
1125 /// ⬈➔●
1126 /// │ ⬈
1127 /// │ ↑
1128 /// │ ↑
1129 /// │ ⬈
1130 /// ●➔⬈───────── x (time)
1131 /// ```
1132 ///
1133 /// To accomplish this, we need to be able to find the position `y` on a curve, given the `x`
1134 /// value. Cubic curves are implicit parametric functions like B(t) = (x,y). To find `y`, we
1135 /// first solve for `t` that corresponds to the given `x` (`time`). We use the Newton-Raphson
1136 /// root-finding method to quickly find a value of `t` that is very near the desired value of
1137 /// `x`. Once we have this we can easily plug that `t` into our curve's `position` function, to
1138 /// find the `y` component, which is how far along our animation should be. In other words:
1139 ///
1140 /// > Given `time` in `0..=1`
1141 ///
1142 /// > Use Newton's method to find a value of `t` that results in B(t) = (x,y) where `x == time`
1143 ///
1144 /// > Once a solution is found, use the resulting `y` value as the final result
1145 #[inline]
1146 pub fn ease(&self, time: f32) -> f32 {
1147 let x = time.clamp(0.0, 1.0);
1148 self.find_y_given_x(x)
1149 }
1150
1151 /// Find the `y` value of the curve at the given `x` value using the Newton-Raphson method.
1152 #[inline]
1153 fn find_y_given_x(&self, x: f32) -> f32 {
1154 let mut t_guess = x;
1155 let mut pos_guess = Vec2::ZERO;
1156 for _ in 0..Self::MAX_ITERS {
1157 pos_guess = self.position(t_guess);
1158 let error = pos_guess.x - x;
1159 if ops::abs(error) <= Self::MAX_ERROR {
1160 break;
1161 }
1162 // Using Newton's method, use the tangent line to estimate a better guess value.
1163 let slope = self.velocity(t_guess).x; // dx/dt
1164 t_guess -= error / slope;
1165 }
1166 pos_guess.y
1167 }
1168}
1169
1170/// A collection of [`CubicSegment`]s chained into a single parametric curve. It is a [`Curve`]
1171/// with domain `[0, N]`, where `N` is its number of segments.
1172///
1173/// Use any struct that implements the [`CubicGenerator`] trait to create a new curve, such as
1174/// [`CubicBezier`].
1175///
1176/// [`Curve`]: crate::curve::Curve
1177#[derive(Clone, Debug, PartialEq)]
1178#[cfg(feature = "alloc")]
1179#[cfg_attr(feature = "serialize", derive(serde::Serialize, serde::Deserialize))]
1180#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
1181pub struct CubicCurve<P: VectorSpace> {
1182 /// The segments comprising the curve. This must always be nonempty.
1183 segments: Vec<CubicSegment<P>>,
1184}
1185
1186#[cfg(feature = "alloc")]
1187impl<P: VectorSpace> CubicCurve<P> {
1188 /// Create a new curve from a collection of segments. If the collection of segments is empty,
1189 /// a curve cannot be built and `None` will be returned instead.
1190 pub fn from_segments(segments: impl IntoIterator<Item = CubicSegment<P>>) -> Option<Self> {
1191 let segments: Vec<_> = segments.into_iter().collect();
1192 if segments.is_empty() {
1193 None
1194 } else {
1195 Some(Self { segments })
1196 }
1197 }
1198
1199 /// Compute the position of a point on the cubic curve at the parametric value `t`.
1200 ///
1201 /// Note that `t` varies from `0..=(n_points - 3)`.
1202 #[inline]
1203 pub fn position(&self, t: f32) -> P {
1204 let (segment, t) = self.segment(t);
1205 segment.position(t)
1206 }
1207
1208 /// Compute the first derivative with respect to t at `t`. This is the instantaneous velocity of
1209 /// a point on the cubic curve at `t`.
1210 ///
1211 /// Note that `t` varies from `0..=(n_points - 3)`.
1212 #[inline]
1213 pub fn velocity(&self, t: f32) -> P {
1214 let (segment, t) = self.segment(t);
1215 segment.velocity(t)
1216 }
1217
1218 /// Compute the second derivative with respect to t at `t`. This is the instantaneous
1219 /// acceleration of a point on the cubic curve at `t`.
1220 ///
1221 /// Note that `t` varies from `0..=(n_points - 3)`.
1222 #[inline]
1223 pub fn acceleration(&self, t: f32) -> P {
1224 let (segment, t) = self.segment(t);
1225 segment.acceleration(t)
1226 }
1227
1228 /// A flexible iterator used to sample curves with arbitrary functions.
1229 ///
1230 /// This splits the curve into `subdivisions` of evenly spaced `t` values across the
1231 /// length of the curve from start (t = 0) to end (t = n), where `n = self.segment_count()`,
1232 /// returning an iterator evaluating the curve with the supplied `sample_function` at each `t`.
1233 ///
1234 /// For `subdivisions = 2`, this will split the curve into two lines, or three points, and
1235 /// return an iterator with 3 items, the three points, one at the start, middle, and end.
1236 #[inline]
1237 pub fn iter_samples<'a, 'b: 'a>(
1238 &'b self,
1239 subdivisions: usize,
1240 mut sample_function: impl FnMut(&Self, f32) -> P + 'a,
1241 ) -> impl Iterator<Item = P> + 'a {
1242 self.iter_uniformly(subdivisions)
1243 .map(move |t| sample_function(self, t))
1244 }
1245
1246 /// An iterator that returns values of `t` uniformly spaced over `0..=subdivisions`.
1247 #[inline]
1248 fn iter_uniformly(&self, subdivisions: usize) -> impl Iterator<Item = f32> {
1249 let segments = self.segments.len() as f32;
1250 let step = segments / subdivisions as f32;
1251 (0..=subdivisions).map(move |i| i as f32 * step)
1252 }
1253
1254 /// The list of segments contained in this `CubicCurve`.
1255 ///
1256 /// This spline's global `t` value is equal to how many segments it has.
1257 ///
1258 /// All method accepting `t` on `CubicCurve` depends on the global `t`.
1259 /// When sampling over the entire curve, you should either use one of the
1260 /// `iter_*` methods or account for the segment count using `curve.segments().len()`.
1261 #[inline]
1262 pub fn segments(&self) -> &[CubicSegment<P>] {
1263 &self.segments
1264 }
1265
1266 /// Iterate over the curve split into `subdivisions`, sampling the position at each step.
1267 pub fn iter_positions(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1268 self.iter_samples(subdivisions, Self::position)
1269 }
1270
1271 /// Iterate over the curve split into `subdivisions`, sampling the velocity at each step.
1272 pub fn iter_velocities(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1273 self.iter_samples(subdivisions, Self::velocity)
1274 }
1275
1276 /// Iterate over the curve split into `subdivisions`, sampling the acceleration at each step.
1277 pub fn iter_accelerations(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1278 self.iter_samples(subdivisions, Self::acceleration)
1279 }
1280
1281 #[inline]
1282 /// Adds a segment to the curve
1283 pub fn push_segment(&mut self, segment: CubicSegment<P>) {
1284 self.segments.push(segment);
1285 }
1286
1287 /// Returns the [`CubicSegment`] and local `t` value given a spline's global `t` value.
1288 #[inline]
1289 fn segment(&self, t: f32) -> (&CubicSegment<P>, f32) {
1290 if self.segments.len() == 1 {
1291 (&self.segments[0], t)
1292 } else {
1293 let i = (ops::floor(t) as usize).clamp(0, self.segments.len() - 1);
1294 (&self.segments[i], t - i as f32)
1295 }
1296 }
1297}
1298
1299#[cfg(feature = "alloc")]
1300impl<P: VectorSpace> Extend<CubicSegment<P>> for CubicCurve<P> {
1301 fn extend<T: IntoIterator<Item = CubicSegment<P>>>(&mut self, iter: T) {
1302 self.segments.extend(iter);
1303 }
1304}
1305
1306#[cfg(feature = "alloc")]
1307impl<P: VectorSpace> IntoIterator for CubicCurve<P> {
1308 type IntoIter = <Vec<CubicSegment<P>> as IntoIterator>::IntoIter;
1309
1310 type Item = CubicSegment<P>;
1311
1312 fn into_iter(self) -> Self::IntoIter {
1313 self.segments.into_iter()
1314 }
1315}
1316
1317/// Implement this on cubic splines that can generate a rational cubic curve from their spline parameters.
1318#[cfg(feature = "alloc")]
1319pub trait RationalGenerator<P: VectorSpace> {
1320 /// An error type indicating why construction might fail.
1321 type Error;
1322
1323 /// Build a [`RationalCurve`] by computing the interpolation coefficients for each curve segment.
1324 fn to_curve(&self) -> Result<RationalCurve<P>, Self::Error>;
1325}
1326
1327/// A segment of a rational cubic curve, used to hold precomputed coefficients for fast interpolation.
1328/// It is a [`Curve`] with domain `[0, 1]`.
1329///
1330/// Note that the `knot_span` is used only by [compound curves] constructed by chaining these
1331/// together.
1332///
1333/// [compound curves]: RationalCurve
1334/// [`Curve`]: crate::curve::Curve
1335#[derive(Copy, Clone, Debug, Default, PartialEq)]
1336#[cfg_attr(feature = "serialize", derive(serde::Serialize, serde::Deserialize))]
1337#[cfg_attr(
1338 feature = "bevy_reflect",
1339 derive(Reflect),
1340 reflect(Debug, Default, Clone)
1341)]
1342pub struct RationalSegment<P: VectorSpace> {
1343 /// The coefficients matrix of the cubic curve.
1344 pub coeff: [P; 4],
1345 /// The homogeneous weight coefficients.
1346 pub weight_coeff: [f32; 4],
1347 /// The width of the domain of this segment.
1348 pub knot_span: f32,
1349}
1350impl<P: VectorSpace> RationalSegment<P> {
1351 /// Instantaneous position of a point at parametric value `t` in `[0, 1]`.
1352 #[inline]
1353 pub fn position(&self, t: f32) -> P {
1354 let [a, b, c, d] = self.coeff;
1355 let [x, y, z, w] = self.weight_coeff;
1356 // Compute a cubic polynomial for the control points
1357 let numerator = a + (b + (c + d * t) * t) * t;
1358 // Compute a cubic polynomial for the weights
1359 let denominator = x + (y + (z + w * t) * t) * t;
1360 numerator / denominator
1361 }
1362
1363 /// Instantaneous velocity of a point at parametric value `t` in `[0, 1]`.
1364 #[inline]
1365 pub fn velocity(&self, t: f32) -> P {
1366 // A derivation for the following equations can be found in "Matrix representation for NURBS
1367 // curves and surfaces" by Choi et al. See equation 19.
1368
1369 let [a, b, c, d] = self.coeff;
1370 let [x, y, z, w] = self.weight_coeff;
1371 // Compute a cubic polynomial for the control points
1372 let numerator = a + (b + (c + d * t) * t) * t;
1373 // Compute a cubic polynomial for the weights
1374 let denominator = x + (y + (z + w * t) * t) * t;
1375
1376 // Compute the derivative of the control point polynomial
1377 let numerator_derivative = b + (c * 2.0 + d * 3.0 * t) * t;
1378 // Compute the derivative of the weight polynomial
1379 let denominator_derivative = y + (z * 2.0 + w * 3.0 * t) * t;
1380
1381 // Velocity is the first derivative (wrt to the parameter `t`)
1382 // Position = N/D therefore
1383 // Velocity = (N/D)' = N'/D - N * D'/D^2 = (N' * D - N * D')/D^2
1384 numerator_derivative / denominator
1385 - numerator * (denominator_derivative / denominator.squared())
1386 }
1387
1388 /// Instantaneous acceleration of a point at parametric value `t` in `[0, 1]`.
1389 #[inline]
1390 pub fn acceleration(&self, t: f32) -> P {
1391 // A derivation for the following equations can be found in "Matrix representation for NURBS
1392 // curves and surfaces" by Choi et al. See equation 20. Note: In come copies of this paper, equation 20
1393 // is printed with the following two errors:
1394 // + The first term has incorrect sign.
1395 // + The second term uses R when it should use the first derivative.
1396
1397 let [a, b, c, d] = self.coeff;
1398 let [x, y, z, w] = self.weight_coeff;
1399 // Compute a cubic polynomial for the control points
1400 let numerator = a + (b + (c + d * t) * t) * t;
1401 // Compute a cubic polynomial for the weights
1402 let denominator = x + (y + (z + w * t) * t) * t;
1403
1404 // Compute the derivative of the control point polynomial
1405 let numerator_derivative = b + (c * 2.0 + d * 3.0 * t) * t;
1406 // Compute the derivative of the weight polynomial
1407 let denominator_derivative = y + (z * 2.0 + w * 3.0 * t) * t;
1408
1409 // Compute the second derivative of the control point polynomial
1410 let numerator_second_derivative = c * 2.0 + d * 6.0 * t;
1411 // Compute the second derivative of the weight polynomial
1412 let denominator_second_derivative = z * 2.0 + w * 6.0 * t;
1413
1414 // Velocity is the first derivative (wrt to the parameter `t`)
1415 // Position = N/D therefore
1416 // Velocity = (N/D)' = N'/D - N * D'/D^2 = (N' * D - N * D')/D^2
1417 // Acceleration = (N/D)'' = ((N' * D - N * D')/D^2)' = N''/D + N' * (-2D'/D^2) + N * (-D''/D^2 + 2D'^2/D^3)
1418 numerator_second_derivative / denominator
1419 + numerator_derivative * (-2.0 * denominator_derivative / denominator.squared())
1420 + numerator
1421 * (-denominator_second_derivative / denominator.squared()
1422 + 2.0 * denominator_derivative.squared() / denominator.cubed())
1423 }
1424
1425 /// Calculate polynomial coefficients for the cubic polynomials using a characteristic matrix.
1426 #[cfg_attr(
1427 not(feature = "alloc"),
1428 expect(
1429 dead_code,
1430 reason = "Method only used when `alloc` feature is enabled."
1431 )
1432 )]
1433 #[inline]
1434 fn coefficients(
1435 control_points: [P; 4],
1436 weights: [f32; 4],
1437 knot_span: f32,
1438 char_matrix: [[f32; 4]; 4],
1439 ) -> Self {
1440 // An explanation of this use can be found in "Matrix representation for NURBS curves and surfaces"
1441 // by Choi et al. See section "Evaluation of NURB Curves and Surfaces", and equation 16.
1442
1443 let [c0, c1, c2, c3] = char_matrix;
1444 let p = control_points;
1445 let w = weights;
1446 // These are the control point polynomial coefficients, computed by multiplying the characteristic
1447 // matrix by the point matrix.
1448 let coeff = [
1449 p[0] * c0[0] + p[1] * c0[1] + p[2] * c0[2] + p[3] * c0[3],
1450 p[0] * c1[0] + p[1] * c1[1] + p[2] * c1[2] + p[3] * c1[3],
1451 p[0] * c2[0] + p[1] * c2[1] + p[2] * c2[2] + p[3] * c2[3],
1452 p[0] * c3[0] + p[1] * c3[1] + p[2] * c3[2] + p[3] * c3[3],
1453 ];
1454 // These are the weight polynomial coefficients, computed by multiplying the characteristic
1455 // matrix by the weight matrix.
1456 let weight_coeff = [
1457 w[0] * c0[0] + w[1] * c0[1] + w[2] * c0[2] + w[3] * c0[3],
1458 w[0] * c1[0] + w[1] * c1[1] + w[2] * c1[2] + w[3] * c1[3],
1459 w[0] * c2[0] + w[1] * c2[1] + w[2] * c2[2] + w[3] * c2[3],
1460 w[0] * c3[0] + w[1] * c3[1] + w[2] * c3[2] + w[3] * c3[3],
1461 ];
1462 Self {
1463 coeff,
1464 weight_coeff,
1465 knot_span,
1466 }
1467 }
1468}
1469
1470/// A collection of [`RationalSegment`]s chained into a single parametric curve. It is a [`Curve`]
1471/// with domain `[0, N]`, where `N` is the number of segments.
1472///
1473/// Use any struct that implements the [`RationalGenerator`] trait to create a new curve, such as
1474/// [`CubicNurbs`], or convert [`CubicCurve`] using `into/from`.
1475///
1476/// [`Curve`]: crate::curve::Curve
1477#[derive(Clone, Debug, PartialEq)]
1478#[cfg(feature = "alloc")]
1479#[cfg_attr(feature = "serialize", derive(serde::Serialize, serde::Deserialize))]
1480#[cfg_attr(feature = "bevy_reflect", derive(Reflect), reflect(Debug, Clone))]
1481pub struct RationalCurve<P: VectorSpace> {
1482 /// The segments comprising the curve. This must always be nonempty.
1483 segments: Vec<RationalSegment<P>>,
1484}
1485
1486#[cfg(feature = "alloc")]
1487impl<P: VectorSpace> RationalCurve<P> {
1488 /// Create a new curve from a collection of segments. If the collection of segments is empty,
1489 /// a curve cannot be built and `None` will be returned instead.
1490 pub fn from_segments(segments: impl IntoIterator<Item = RationalSegment<P>>) -> Option<Self> {
1491 let segments: Vec<_> = segments.into_iter().collect();
1492 if segments.is_empty() {
1493 None
1494 } else {
1495 Some(Self { segments })
1496 }
1497 }
1498
1499 /// Compute the position of a point on the curve at the parametric value `t`.
1500 ///
1501 /// Note that `t` varies from `0` to `self.length()`.
1502 #[inline]
1503 pub fn position(&self, t: f32) -> P {
1504 let (segment, t) = self.segment(t);
1505 segment.position(t)
1506 }
1507
1508 /// Compute the first derivative with respect to t at `t`. This is the instantaneous velocity of
1509 /// a point on the curve at `t`.
1510 ///
1511 /// Note that `t` varies from `0` to `self.length()`.
1512 #[inline]
1513 pub fn velocity(&self, t: f32) -> P {
1514 let (segment, t) = self.segment(t);
1515 segment.velocity(t)
1516 }
1517
1518 /// Compute the second derivative with respect to t at `t`. This is the instantaneous
1519 /// acceleration of a point on the curve at `t`.
1520 ///
1521 /// Note that `t` varies from `0` to `self.length()`.
1522 #[inline]
1523 pub fn acceleration(&self, t: f32) -> P {
1524 let (segment, t) = self.segment(t);
1525 segment.acceleration(t)
1526 }
1527
1528 /// A flexible iterator used to sample curves with arbitrary functions.
1529 ///
1530 /// This splits the curve into `subdivisions` of evenly spaced `t` values across the
1531 /// length of the curve from start (t = 0) to end (t = n), where `n = self.segment_count()`,
1532 /// returning an iterator evaluating the curve with the supplied `sample_function` at each `t`.
1533 ///
1534 /// For `subdivisions = 2`, this will split the curve into two lines, or three points, and
1535 /// return an iterator with 3 items, the three points, one at the start, middle, and end.
1536 #[inline]
1537 pub fn iter_samples<'a, 'b: 'a>(
1538 &'b self,
1539 subdivisions: usize,
1540 mut sample_function: impl FnMut(&Self, f32) -> P + 'a,
1541 ) -> impl Iterator<Item = P> + 'a {
1542 self.iter_uniformly(subdivisions)
1543 .map(move |t| sample_function(self, t))
1544 }
1545
1546 /// An iterator that returns values of `t` uniformly spaced over `0..=subdivisions`.
1547 #[inline]
1548 fn iter_uniformly(&self, subdivisions: usize) -> impl Iterator<Item = f32> {
1549 let length = self.length();
1550 let step = length / subdivisions as f32;
1551 (0..=subdivisions).map(move |i| i as f32 * step)
1552 }
1553
1554 /// The list of segments contained in this `RationalCurve`.
1555 ///
1556 /// This spline's global `t` value is equal to how many segments it has.
1557 ///
1558 /// All method accepting `t` on `RationalCurve` depends on the global `t`.
1559 /// When sampling over the entire curve, you should either use one of the
1560 /// `iter_*` methods or account for the segment count using `curve.segments().len()`.
1561 #[inline]
1562 pub fn segments(&self) -> &[RationalSegment<P>] {
1563 &self.segments
1564 }
1565
1566 /// Iterate over the curve split into `subdivisions`, sampling the position at each step.
1567 pub fn iter_positions(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1568 self.iter_samples(subdivisions, Self::position)
1569 }
1570
1571 /// Iterate over the curve split into `subdivisions`, sampling the velocity at each step.
1572 pub fn iter_velocities(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1573 self.iter_samples(subdivisions, Self::velocity)
1574 }
1575
1576 /// Iterate over the curve split into `subdivisions`, sampling the acceleration at each step.
1577 pub fn iter_accelerations(&self, subdivisions: usize) -> impl Iterator<Item = P> + '_ {
1578 self.iter_samples(subdivisions, Self::acceleration)
1579 }
1580
1581 /// Adds a segment to the curve.
1582 #[inline]
1583 pub fn push_segment(&mut self, segment: RationalSegment<P>) {
1584 self.segments.push(segment);
1585 }
1586
1587 /// Returns the [`RationalSegment`] and local `t` value given a spline's global `t` value.
1588 /// Input `t` will be clamped to the domain of the curve. Returned value will be in `[0, 1]`.
1589 #[inline]
1590 fn segment(&self, mut t: f32) -> (&RationalSegment<P>, f32) {
1591 if t <= 0.0 {
1592 (&self.segments[0], 0.0)
1593 } else if self.segments.len() == 1 {
1594 (&self.segments[0], t / self.segments[0].knot_span)
1595 } else {
1596 // Try to fit t into each segment domain
1597 for segment in self.segments.iter() {
1598 if t < segment.knot_span {
1599 // The division here makes t a normalized parameter in [0, 1] that can be properly
1600 // evaluated against a rational curve segment. See equations 6 & 16 from "Matrix representation
1601 // of NURBS curves and surfaces" by Choi et al. or equation 3 from "General Matrix
1602 // Representations for B-Splines" by Qin.
1603 return (segment, t / segment.knot_span);
1604 }
1605 t -= segment.knot_span;
1606 }
1607 return (self.segments.last().unwrap(), 1.0);
1608 }
1609 }
1610
1611 /// Returns the length of the domain of the parametric curve.
1612 #[inline]
1613 pub fn length(&self) -> f32 {
1614 self.segments.iter().map(|segment| segment.knot_span).sum()
1615 }
1616}
1617
1618#[cfg(feature = "alloc")]
1619impl<P: VectorSpace> Extend<RationalSegment<P>> for RationalCurve<P> {
1620 fn extend<T: IntoIterator<Item = RationalSegment<P>>>(&mut self, iter: T) {
1621 self.segments.extend(iter);
1622 }
1623}
1624
1625#[cfg(feature = "alloc")]
1626impl<P: VectorSpace> IntoIterator for RationalCurve<P> {
1627 type IntoIter = <Vec<RationalSegment<P>> as IntoIterator>::IntoIter;
1628
1629 type Item = RationalSegment<P>;
1630
1631 fn into_iter(self) -> Self::IntoIter {
1632 self.segments.into_iter()
1633 }
1634}
1635
1636impl<P: VectorSpace> From<CubicSegment<P>> for RationalSegment<P> {
1637 fn from(value: CubicSegment<P>) -> Self {
1638 Self {
1639 coeff: value.coeff,
1640 weight_coeff: [1.0, 0.0, 0.0, 0.0],
1641 knot_span: 1.0, // Cubic curves are uniform, so every segment has domain [0, 1).
1642 }
1643 }
1644}
1645
1646#[cfg(feature = "alloc")]
1647impl<P: VectorSpace> From<CubicCurve<P>> for RationalCurve<P> {
1648 fn from(value: CubicCurve<P>) -> Self {
1649 Self {
1650 segments: value.segments.into_iter().map(Into::into).collect(),
1651 }
1652 }
1653}
1654
1655#[cfg(feature = "alloc")]
1656#[cfg(test)]
1657mod tests {
1658 use crate::{
1659 cubic_splines::{
1660 CubicBSpline, CubicBezier, CubicGenerator, CubicNurbs, CubicSegment, RationalCurve,
1661 RationalGenerator,
1662 },
1663 ops::{self, FloatPow},
1664 };
1665 use alloc::vec::Vec;
1666 use glam::{vec2, Vec2};
1667
1668 /// How close two floats can be and still be considered equal
1669 const FLOAT_EQ: f32 = 1e-5;
1670
1671 /// Sweep along the full length of a 3D cubic Bezier, and manually check the position.
1672 #[test]
1673 fn cubic() {
1674 const N_SAMPLES: usize = 1000;
1675 let points = [[
1676 vec2(-1.0, -20.0),
1677 vec2(3.0, 2.0),
1678 vec2(5.0, 3.0),
1679 vec2(9.0, 8.0),
1680 ]];
1681 let bezier = CubicBezier::new(points).to_curve().unwrap();
1682 for i in 0..=N_SAMPLES {
1683 let t = i as f32 / N_SAMPLES as f32; // Check along entire length
1684 assert!(bezier.position(t).distance(cubic_manual(t, points[0])) <= FLOAT_EQ);
1685 }
1686 }
1687
1688 /// Manual, hardcoded function for computing the position along a cubic bezier.
1689 fn cubic_manual(t: f32, points: [Vec2; 4]) -> Vec2 {
1690 let p = points;
1691 p[0] * (1.0 - t).cubed()
1692 + 3.0 * p[1] * t * (1.0 - t).squared()
1693 + 3.0 * p[2] * t.squared() * (1.0 - t)
1694 + p[3] * t.cubed()
1695 }
1696
1697 /// Basic cubic Bezier easing test to verify the shape of the curve.
1698 #[test]
1699 fn easing_simple() {
1700 // A curve similar to ease-in-out, but symmetric
1701 let bezier = CubicSegment::new_bezier_easing([1.0, 0.0], [0.0, 1.0]);
1702 assert_eq!(bezier.ease(0.0), 0.0);
1703 assert!(bezier.ease(0.2) < 0.2); // tests curve
1704 assert_eq!(bezier.ease(0.5), 0.5); // true due to symmetry
1705 assert!(bezier.ease(0.8) > 0.8); // tests curve
1706 assert_eq!(bezier.ease(1.0), 1.0);
1707 }
1708
1709 /// A curve that forms an upside-down "U", that should extend below 0.0. Useful for animations
1710 /// that go beyond the start and end positions, e.g. bouncing.
1711 #[test]
1712 fn easing_overshoot() {
1713 // A curve that forms an upside-down "U", that should extend above 1.0
1714 let bezier = CubicSegment::new_bezier_easing([0.0, 2.0], [1.0, 2.0]);
1715 assert_eq!(bezier.ease(0.0), 0.0);
1716 assert!(bezier.ease(0.5) > 1.5);
1717 assert_eq!(bezier.ease(1.0), 1.0);
1718 }
1719
1720 /// A curve that forms a "U", that should extend below 0.0. Useful for animations that go beyond
1721 /// the start and end positions, e.g. bouncing.
1722 #[test]
1723 fn easing_undershoot() {
1724 let bezier = CubicSegment::new_bezier_easing([0.0, -2.0], [1.0, -2.0]);
1725 assert_eq!(bezier.ease(0.0), 0.0);
1726 assert!(bezier.ease(0.5) < -0.5);
1727 assert_eq!(bezier.ease(1.0), 1.0);
1728 }
1729
1730 /// Test that a simple cardinal spline passes through all of its control points with
1731 /// the correct tangents.
1732 #[test]
1733 fn cardinal_control_pts() {
1734 use super::CubicCardinalSpline;
1735
1736 let tension = 0.2;
1737 let [p0, p1, p2, p3] = [vec2(-1., -2.), vec2(0., 1.), vec2(1., 2.), vec2(-2., 1.)];
1738 let curve = CubicCardinalSpline::new(tension, [p0, p1, p2, p3])
1739 .to_curve()
1740 .unwrap();
1741
1742 // Positions at segment endpoints
1743 assert!(curve.position(0.).abs_diff_eq(p0, FLOAT_EQ));
1744 assert!(curve.position(1.).abs_diff_eq(p1, FLOAT_EQ));
1745 assert!(curve.position(2.).abs_diff_eq(p2, FLOAT_EQ));
1746 assert!(curve.position(3.).abs_diff_eq(p3, FLOAT_EQ));
1747
1748 // Tangents at segment endpoints
1749 assert!(curve
1750 .velocity(0.)
1751 .abs_diff_eq((p1 - p0) * tension * 2., FLOAT_EQ));
1752 assert!(curve
1753 .velocity(1.)
1754 .abs_diff_eq((p2 - p0) * tension, FLOAT_EQ));
1755 assert!(curve
1756 .velocity(2.)
1757 .abs_diff_eq((p3 - p1) * tension, FLOAT_EQ));
1758 assert!(curve
1759 .velocity(3.)
1760 .abs_diff_eq((p3 - p2) * tension * 2., FLOAT_EQ));
1761 }
1762
1763 /// Test that [`RationalCurve`] properly generalizes [`CubicCurve`]. A Cubic upgraded to a rational
1764 /// should produce pretty much the same output.
1765 #[test]
1766 fn cubic_to_rational() {
1767 const EPSILON: f32 = 0.00001;
1768
1769 let points = [
1770 vec2(0.0, 0.0),
1771 vec2(1.0, 1.0),
1772 vec2(1.0, 1.0),
1773 vec2(2.0, -1.0),
1774 vec2(3.0, 1.0),
1775 vec2(0.0, 0.0),
1776 ];
1777
1778 let b_spline = CubicBSpline::new(points).to_curve().unwrap();
1779 let rational_b_spline = RationalCurve::from(b_spline.clone());
1780
1781 /// Tests if two vectors of points are approximately the same
1782 fn compare_vectors(cubic_curve: Vec<Vec2>, rational_curve: Vec<Vec2>, name: &str) {
1783 assert_eq!(
1784 cubic_curve.len(),
1785 rational_curve.len(),
1786 "{name} vector lengths mismatch"
1787 );
1788 for (i, (a, b)) in cubic_curve.iter().zip(rational_curve.iter()).enumerate() {
1789 assert!(
1790 a.distance(*b) < EPSILON,
1791 "Mismatch at {name} value {i}. CubicCurve: {} Converted RationalCurve: {}",
1792 a,
1793 b
1794 );
1795 }
1796 }
1797
1798 // Both curves should yield the same values
1799 let cubic_positions: Vec<_> = b_spline.iter_positions(10).collect();
1800 let rational_positions: Vec<_> = rational_b_spline.iter_positions(10).collect();
1801 compare_vectors(cubic_positions, rational_positions, "position");
1802
1803 let cubic_velocities: Vec<_> = b_spline.iter_velocities(10).collect();
1804 let rational_velocities: Vec<_> = rational_b_spline.iter_velocities(10).collect();
1805 compare_vectors(cubic_velocities, rational_velocities, "velocity");
1806
1807 let cubic_accelerations: Vec<_> = b_spline.iter_accelerations(10).collect();
1808 let rational_accelerations: Vec<_> = rational_b_spline.iter_accelerations(10).collect();
1809 compare_vectors(cubic_accelerations, rational_accelerations, "acceleration");
1810 }
1811
1812 /// Test that a nurbs curve can approximate a portion of a circle.
1813 #[test]
1814 fn nurbs_circular_arc() {
1815 use core::f32::consts::FRAC_PI_2;
1816 const EPSILON: f32 = 0.0000001;
1817
1818 // The following NURBS parameters were determined by constraining the first two
1819 // points to the line y=1, the second two points to the line x=1, and the distance
1820 // between each pair of points to be equal. One can solve the weights by assuming the
1821 // first and last weights to be one, the intermediate weights to be equal, and
1822 // subjecting ones self to a lot of tedious matrix algebra.
1823
1824 let alpha = FRAC_PI_2;
1825 let leg = 2.0 * ops::sin(alpha / 2.0) / (1.0 + 2.0 * ops::cos(alpha / 2.0));
1826 let weight = (1.0 + 2.0 * ops::cos(alpha / 2.0)) / 3.0;
1827 let points = [
1828 vec2(1.0, 0.0),
1829 vec2(1.0, leg),
1830 vec2(leg, 1.0),
1831 vec2(0.0, 1.0),
1832 ];
1833 let weights = [1.0, weight, weight, 1.0];
1834 let knots = [0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0];
1835 let spline = CubicNurbs::new(points, Some(weights), Some(knots)).unwrap();
1836 let curve = spline.to_curve().unwrap();
1837 for (i, point) in curve.iter_positions(10).enumerate() {
1838 assert!(
1839 ops::abs(point.length() - 1.0) < EPSILON,
1840 "Point {i} is not on the unit circle: {point:?} has length {}",
1841 point.length()
1842 );
1843 }
1844 }
1845}