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