peroxide/numerical/spline.rs
1//! Spline interpolations
2//!
3//! # Available splines
4//!
5//! * Cubic spline
6//! * Cubic Hermite spline
7//! * B-spline
8//!
9//! # `Spline<T>` trait
10//!
11//! ## Methods
12//!
13//! * `fn eval(&self, t: f64) -> T` : Evaluate the spline at t
14//! * For Cubic Spline, t means x (domain) and `T = f64`
15//! * For B-Spline, t means parameter of curve and `T = (f64, f64)`
16//! * `fn eval_vec(&self, v: &[f64]) -> Vec<T>` : Evaluate spline values for an array v
17//! * `fn eval_with_cond<F>(&self, t: f64, cond: F) -> T` : Evaluate the spline at t, with a condition
18//! * `fn eval_vec_with_cond<F>(&self, v: &[f64], cond: F) -> Vec<T>` : Evaluate spline values for an array v, with a condition
19//!
20//! # `PolynomialSpline` trait
21//!
22//! ## Methods
23//!
24//! * `fn polynomial_at(&self, x: f64) -> &Polynomial` : Get the polynomial at x
25//! * `fn number_of_polynomials(&self) -> usize` : Get the number of polynomials
26//! * `fn get_ranged_polynomials(&self) -> &Vec<(Range<f64>, Polynomial)>` : Get the polynomials
27//!
28//! # Low-level interface
29//!
30//! ## Members
31//!
32//! * `CubicSpline`: Structure for cubic spline
33//! * `fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result<Self>` : Create a cubic spline from nodes
34//! * `fn extend_with_nodes(&mut self, node_x: Vec<f64>, node_y: Vec<f64>) -> Result<()>` : Extend the spline with nodes
35//! * `CubicHermiteSpline`: Structure for cubic Hermite spline
36//! * `fn from_nodes_with_slopes(node_x: &[f64], node_y: &[f64], m: &[f64]) -> Result<Self>` : Create a Cubic Hermite spline from nodes with slopes
37//! * `fn from_nodes(node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod) -> Result<Self>` : Create a Cubic Hermite spline from nodes with slope estimation methods
38//! * `SlopeMethod`: Enum for slope estimation methods
39//! * `Akima`: Akima's method to estimate slopes ([Akima (1970)](https://dl.acm.org/doi/abs/10.1145/321607.321609))
40//! * `Quadratic`: Using quadratic interpolation to estimate slopes
41//! * `BSpline`: Structure for B-Spline
42//! * `fn open(degree: usize, knots: Vec<f64>, control_points: Vec<Vec<f64>>) -> Result<Self>` : Create an open B-Spline
43//! * `fn clamped(degree: usize, knots: Vec<f64>, control_points: Vec<Vec<f64>>) -> Result<Self>`
44//! : Create a clamped B-Spline
45//! * `fn cox_de_boor(t: f64, i: f64)` : Cox-de Boor recursion formula (Here, use iteration
46//! instead of recursion)
47//!
48//! ## Usage (Cubic Spline Family)
49//!
50//! ```rust
51//! use peroxide::fuga::*;
52//!
53//! fn main() -> Result<(), Box<dyn Error>> {
54//! let x = seq(0, 10, 1);
55//! let y = x.fmap(|t| t.sin());
56//!
57//! let cs = CubicSpline::from_nodes(&x, &y)?;
58//! let cs_akima = CubicHermiteSpline::from_nodes(&x, &y, Akima)?;
59//! let cs_quad = CubicHermiteSpline::from_nodes(&x, &y, Quadratic)?;
60//!
61//! cs.polynomial_at(0f64).print();
62//! cs_akima.polynomial_at(0f64).print();
63//! cs_quad.polynomial_at(0f64).print();
64//! // -0.1523x^3 + 0.9937x
65//! // 0.1259x^3 - 0.5127x^2 + 1.2283x
66//! // -0.0000x^3 - 0.3868x^2 + 1.2283x
67//!
68//! let new_x = seq(4, 6, 0.1);
69//! let new_y = new_x.fmap(|t| t.sin());
70//!
71//! let y_cs = cs.eval_vec(&new_x);
72//! let y_akima = cs_akima.eval_vec(&new_x);
73//! let y_quad = cs_quad.eval_vec(&new_x);
74//!
75//! let mut df = DataFrame::new(vec![]);
76//! df.push("x", Series::new(new_x));
77//! df.push("y", Series::new(new_y));
78//! df.push("y_cs", Series::new(y_cs));
79//! df.push("y_akima", Series::new(y_akima));
80//! df.push("y_quad", Series::new(y_quad));
81//!
82//! df.print();
83//! // x y y_cs y_akima y_quad
84//! // r[0] 5 -0.9589 -0.9589 -0.9589 -0.9589
85//! // r[1] 5.2 -0.8835 -0.8826 -0.8583 -0.8836
86//! // r[2] 5.4 -0.7728 -0.7706 -0.7360 -0.7629
87//! // r[3] 5.6 -0.6313 -0.6288 -0.5960 -0.6120
88//! // r[4] 5.8 -0.4646 -0.4631 -0.4424 -0.4459
89//! // r[5] 6 -0.2794 -0.2794 -0.2794 -0.2794
90//!
91//! Ok(())
92//! }
93//! ```
94//!
95//! ## Usage (B-Spline)
96//!
97//! ```rust
98//! use peroxide::fuga::*;
99//!
100//! fn main() -> Result<(), Box<dyn Error>> {
101//! let knots = vec![0f64, 1f64, 2f64, 3f64];
102//! let degree = 3;
103//! let control_points = vec![
104//! vec![0f64, 2f64],
105//! vec![0.2, -1f64],
106//! vec![0.4, 1f64],
107//! vec![0.6, -1f64],
108//! vec![0.8, 1f64],
109//! vec![1f64, 2f64],
110//! ];
111//!
112//!
113//! let spline = BSpline::clamped(degree, knots, control_points.clone())?;
114//! let t = linspace(0f64, 3f64, 200);
115//! let (x, y): (Vec<f64>, Vec<f64>) = spline.eval_vec(&t).into_iter().unzip();
116//!
117//! # #[cfg(feature = "plot")]
118//! # {
119//! let control_x = control_points.iter().map(|v| v[0]).collect::<Vec<f64>>();
120//! let control_y = control_points.iter().map(|v| v[1]).collect::<Vec<f64>>();
121//!
122//! let mut plt = Plot2D::new();
123//! plt
124//! .insert_pair((x.clone(), y.clone()))
125//! .insert_pair((control_x.clone(), control_y.clone()))
126//! .set_plot_type(vec![(1, PlotType::Scatter)])
127//! .set_color(vec![(0, "darkblue"), (1, "red")])
128//! .set_xlabel(r"$x$")
129//! .set_ylabel(r"$y$")
130//! .set_style(PlotStyle::Nature)
131//! .set_dpi(600)
132//! .set_path("example_data/b_spline_test.png")
133//! .savefig()?;
134//! # }
135//!
136//! let mut df = DataFrame::new(vec![]);
137//! df.push("t", Series::new(t));
138//! df.push("x", Series::new(x));
139//! df.push("y", Series::new(y));
140//! df.print();
141//!
142//! Ok(())
143//! }
144//! ```
145//!
146//! - Result for above code is:
147//! 
148//!
149//! # High-level interface
150//!
151//! ## Functions
152//!
153//! * `fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> CubicSpline` : Create a cubic spline from nodes
154//! * `fn cubic_hermite_spline(node_x: &[f64], node_y: &[f64], m: &[f64]) -> CubicHermiteSpline` : Create a cubic Hermite spline from nodes with slopes
155//!
156//! ## Usage
157//!
158//! ```rust
159//! use peroxide::fuga::*;
160//!
161//! fn main() -> Result<(), Box<dyn Error>> {
162//! let x = seq(0, 10, 1);
163//! let y = x.fmap(|t| t.sin());
164//!
165//! let cs = cubic_spline(&x, &y)?;
166//! let cs_akima = cubic_hermite_spline(&x, &y, Akima)?;
167//! let cs_quad = cubic_hermite_spline(&x, &y, Quadratic)?;
168//!
169//! cs.polynomial_at(0f64).print();
170//! cs_akima.polynomial_at(0f64).print();
171//! cs_quad.polynomial_at(0f64).print();
172//! // -0.1523x^3 + 0.9937x
173//! // 0.1259x^3 - 0.5127x^2 + 1.2283x
174//! // -0.0000x^3 - 0.3868x^2 + 1.2283x
175//!
176//! let new_x = seq(4, 6, 0.1);
177//! let new_y = new_x.fmap(|t| t.sin());
178//!
179//! let y_cs = cs.eval_vec(&new_x);
180//! let y_akima = cs_akima.eval_vec(&new_x);
181//! let y_quad = cs_quad.eval_vec(&new_x);
182//!
183//! let mut df = DataFrame::new(vec![]);
184//! df.push("x", Series::new(new_x));
185//! df.push("y", Series::new(new_y));
186//! df.push("y_cs", Series::new(y_cs));
187//! df.push("y_akima", Series::new(y_akima));
188//! df.push("y_quad", Series::new(y_quad));
189//!
190//! df.print();
191//! // x y y_cs y_akima y_quad
192//! // r[0] 5 -0.9589 -0.9589 -0.9589 -0.9589
193//! // r[1] 5.2 -0.8835 -0.8826 -0.8583 -0.8836
194//! // r[2] 5.4 -0.7728 -0.7706 -0.7360 -0.7629
195//! // r[3] 5.6 -0.6313 -0.6288 -0.5960 -0.6120
196//! // r[4] 5.8 -0.4646 -0.4631 -0.4424 -0.4459
197//! // r[5] 6 -0.2794 -0.2794 -0.2794 -0.2794
198//!
199//! Ok(())
200//! }
201//! ```
202//!
203//! # Calculus with polynomial splines
204//!
205//! ## Usage
206//!
207//! ```rust
208//! use peroxide::fuga::*;
209//! use std::f64::consts::PI;
210//!
211//! fn main() -> Result<(), Box<dyn Error>> {
212//! let x = seq(0, 10, 1);
213//! let y = x.fmap(|t| t.sin());
214//!
215//! let cs = cubic_spline(&x, &y)?;
216//! let cs_akima = cubic_hermite_spline(&x, &y, Akima)?;
217//! let cs_quad = cubic_hermite_spline(&x, &y, Quadratic)?;
218//!
219//! println!("============ Polynomial at x=0 ============");
220//!
221//! cs.polynomial_at(0f64).print();
222//! cs_akima.polynomial_at(0f64).print();
223//! cs_quad.polynomial_at(0f64).print();
224//!
225//! println!("============ Derivative at x=0 ============");
226//!
227//! cs.derivative().polynomial_at(0f64).print();
228//! cs_akima.derivative().polynomial_at(0f64).print();
229//! cs_quad.derivative().polynomial_at(0f64).print();
230//!
231//! println!("============ Integral at x=0 ============");
232//!
233//! cs.integral().polynomial_at(0f64).print();
234//! cs_akima.integral().polynomial_at(0f64).print();
235//! cs_quad.integral().polynomial_at(0f64).print();
236//!
237//! println!("============ Integrate from x=0 to x=pi ============");
238//!
239//! cs.integrate((0f64, PI)).print();
240//! cs_akima.integrate((0f64, PI)).print();
241//! cs_quad.integrate((0f64, PI)).print();
242//!
243//! // ============ Polynomial at x=0 ============
244//! // -0.1523x^3 + 0.9937x
245//! // 0.1259x^3 - 0.5127x^2 + 1.2283x
246//! // -0.0000x^3 - 0.3868x^2 + 1.2283x
247//! // ============ Derivative at x=0 ============
248//! // -0.4568x^2 + 0.9937
249//! // 0.3776x^2 - 1.0254x + 1.2283
250//! // -0.0000x^2 - 0.7736x + 1.2283
251//! // ============ Integral at x=0 ============
252//! // -0.0381x^4 + 0.4969x^2
253//! // 0.0315x^4 - 0.1709x^3 + 0.6141x^2
254//! // -0.0000x^4 - 0.1289x^3 + 0.6141x^2
255//! // ============ Integrate from x=0 to x=pi ============
256//! // 1.9961861265456702
257//! // 2.0049920614062775
258//! // 2.004327391790717
259//!
260//! Ok(())
261//! }
262//! ```
263//!
264//! # B-Spline utils
265//!
266//! - `UnitCubicBasis`: Single cubic B-Spline basis function
267//! - `CubicBSplineBases`: Uniform Cubic B-Spline basis functions
268//!
269//! # References
270//!
271//! - Gary D. Knott, *Interpolating Splines*, Birkhäuser Boston, MA, (2000).
272/// - [Wikipedia - Irwin-Hall distribution](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution#Special_cases)
273use self::SplineError::{NotEnoughNodes, NotEqualNodes, NotEqualSlopes, RedundantNodeX};
274#[allow(unused_imports)]
275use crate::structure::matrix::*;
276#[allow(unused_imports)]
277use crate::structure::polynomial::*;
278#[allow(unused_imports)]
279use crate::structure::vector::*;
280use crate::traits::matrix::LinearAlgebra;
281#[allow(unused_imports)]
282use crate::util::non_macro::*;
283use crate::util::useful::zip_range;
284use anyhow::{bail, Result};
285use peroxide_num::PowOps;
286#[cfg(feature = "serde")]
287use serde::{Deserialize, Serialize};
288use std::cmp::{max, min};
289use std::convert::From;
290use std::ops::{Index, Range};
291
292pub trait Spline<T> {
293 fn eval(&self, t: f64) -> T;
294 fn eval_vec(&self, v: &[f64]) -> Vec<T> {
295 v.iter().map(|&t| self.eval(t)).collect()
296 }
297 fn eval_with_cond<F: Fn(T) -> T>(&self, t: f64, cond: F) -> T {
298 cond(self.eval(t))
299 }
300 fn eval_vec_with_cond<F: Fn(T) -> T + Copy>(&self, t: &[f64], cond: F) -> Vec<T> {
301 t.iter().map(|&x| self.eval_with_cond(x, cond)).collect()
302 }
303}
304
305/// Trait for spline interpolation
306///
307/// # Available Splines
308///
309/// - `CubicSpline`
310/// - `CubicHermiteSpline`
311impl<P: PolynomialSpline> Spline<f64> for P {
312 fn eval(&self, x: f64) -> f64 {
313 self.polynomial_at(x).eval(x)
314 }
315}
316
317pub trait PolynomialSpline {
318 fn polynomial_at<T: Into<f64> + Copy>(&self, x: T) -> &Polynomial {
319 let x = x.into();
320
321 let poly = self.get_ranged_polynomials();
322
323 let index = match poly.binary_search_by(|(range, _)| {
324 if range.contains(&x) {
325 core::cmp::Ordering::Equal
326 } else if x < range.start {
327 core::cmp::Ordering::Greater
328 } else {
329 core::cmp::Ordering::Less
330 }
331 }) {
332 Ok(index) => index,
333 Err(index) => max(0, min(index, poly.len() - 1)),
334 };
335
336 &poly[index].1
337 }
338
339 fn number_of_polynomials(&self) -> usize {
340 self.get_ranged_polynomials().len()
341 }
342
343 fn get_ranged_polynomials(&self) -> &Vec<(Range<f64>, Polynomial)>;
344}
345
346// =============================================================================
347// High level functions
348// =============================================================================
349/// Cubic Spline (Natural)
350///
351/// # Description
352///
353/// Implement traits of Natural cubic splines, by Arne Morten Kvarving.
354///
355/// # Type
356/// `(&[f64], &[f64]) -> Cubic Spline`
357///
358/// # Examples
359/// ```
360/// #[macro_use]
361/// extern crate peroxide;
362/// use peroxide::fuga::*;
363///
364/// fn main() -> Result<(), Box<dyn Error>> {
365/// let x = c!(0.9, 1.3, 1.9, 2.1);
366/// let y = c!(1.3, 1.5, 1.85, 2.1);
367///
368/// let s = cubic_spline(&x, &y)?;
369///
370/// let new_x = c!(1, 1.5, 2.0);
371///
372/// // Generate Cubic polynomial
373/// for t in new_x.iter() {
374/// s.polynomial_at(*t).print();
375/// }
376/// // -0.2347x^3 + 0.6338x^2 - 0.0329x + 0.9873
377/// // 0.9096x^3 - 3.8292x^2 + 5.7691x - 1.5268
378/// // -2.2594x^3 + 14.2342x^2 - 28.5513x + 20.2094
379///
380/// // Evaluation
381/// for t in new_x.iter() {
382/// s.eval(*t).print();
383/// }
384///
385/// Ok(())
386/// }
387/// ```
388pub fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> Result<CubicSpline> {
389 CubicSpline::from_nodes(node_x, node_y)
390}
391
392pub fn cubic_hermite_spline(
393 node_x: &[f64],
394 node_y: &[f64],
395 slope_method: SlopeMethod,
396) -> Result<CubicHermiteSpline> {
397 CubicHermiteSpline::from_nodes(node_x, node_y, slope_method)
398}
399
400// =============================================================================
401// Cubic Spline
402// =============================================================================
403/// Cubic Spline (Natural)
404///
405/// # Description
406///
407/// Implement traits of Natural cubic splines, by Arne Morten Kvarving.
408///
409/// # Type
410/// `(&[f64], &[f64]) -> Cubic Spline`
411///
412/// # Examples
413/// ```
414/// #[macro_use]
415/// extern crate peroxide;
416/// use peroxide::fuga::*;
417///
418/// fn main() -> Result<(), Box<dyn Error>> {
419/// let x = c!(0.9, 1.3, 1.9, 2.1);
420/// let y = c!(1.3, 1.5, 1.85, 2.1);
421///
422/// let s = cubic_spline(&x, &y)?;
423///
424/// let new_x = c!(1, 1.5, 2.0);
425///
426/// // Generate Cubic polynomial
427/// for t in new_x.iter() {
428/// s.polynomial_at(*t).print();
429/// }
430/// // -0.2347x^3 + 0.6338x^2 - 0.0329x + 0.9873
431/// // 0.9096x^3 - 3.8292x^2 + 5.7691x - 1.5268
432/// // -2.2594x^3 + 14.2342x^2 - 28.5513x + 20.2094
433///
434/// // Evaluation
435/// for t in new_x.iter() {
436/// s.eval(*t).print();
437/// }
438///
439/// Ok(())
440/// }
441/// ```
442#[derive(Debug, Clone, Default)]
443#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
444#[cfg_attr(
445 feature = "rkyv",
446 derive(rkyv::Archive, rkyv::Deserialize, rkyv::Serialize)
447)]
448pub struct CubicSpline {
449 polynomials: Vec<(Range<f64>, Polynomial)>,
450}
451
452impl PolynomialSpline for CubicSpline {
453 fn get_ranged_polynomials(&self) -> &Vec<(Range<f64>, Polynomial)> {
454 &self.polynomials
455 }
456}
457
458#[derive(Debug, Copy, Clone)]
459pub enum SplineError {
460 NotEnoughNodes,
461 NotEqualNodes,
462 NotEqualSlopes,
463 RedundantNodeX,
464}
465
466impl std::fmt::Display for SplineError {
467 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
468 match self {
469 SplineError::NotEnoughNodes => write!(f, "node_x has less than 3 elements"),
470 SplineError::NotEqualNodes => write!(f, "node_x and node_y have different lengths"),
471 SplineError::NotEqualSlopes => write!(f, "nodes and slopes have different lengths"),
472 SplineError::RedundantNodeX => write!(f, "there are redundant nodes in node_x"),
473 }
474 }
475}
476
477impl CubicSpline {
478 /// # Examples
479 /// ```
480 /// #[macro_use]
481 /// extern crate peroxide;
482 /// use peroxide::fuga::*;
483 ///
484 /// fn main() -> Result<(), Box<dyn Error>> {
485 /// let x = c!(0.9, 1.3, 1.9, 2.1);
486 /// let y = c!(1.3, 1.5, 1.85, 2.1);
487 ///
488 /// let s = CubicSpline::from_nodes(&x, &y)?;
489 ///
490 /// for i in 0 .. 4 {
491 /// println!("{}", s.eval(i as f64 / 2.0));
492 /// }
493 ///
494 /// Ok(())
495 /// }
496 /// ```
497 pub fn from_nodes(node_x: &[f64], node_y: &[f64]) -> Result<Self> {
498 let polynomials = CubicSpline::cubic_spline(node_x, node_y)?;
499 Ok(CubicSpline {
500 polynomials: zip_range(node_x, &polynomials),
501 })
502 }
503
504 fn cubic_spline(node_x: &[f64], node_y: &[f64]) -> Result<Vec<Polynomial>> {
505 //! Pre calculated variables
506 //! node_x: n+1
507 //! node_y: n+1
508 //! h : n
509 //! b : n
510 //! v : n
511 //! u : n
512 //! z : n+1
513 let n = node_x.len() - 1;
514 if n < 2 {
515 bail!(NotEnoughNodes);
516 }
517 if n != node_y.len() - 1 {
518 bail!(NotEqualNodes);
519 }
520
521 // Pre-calculations
522 let mut h = vec![0f64; n];
523 let mut b = vec![0f64; n];
524 let mut v = vec![0f64; n];
525 let mut u = vec![0f64; n];
526 for i in 0..n {
527 if i == 0 {
528 h[i] = node_x[i + 1] - node_x[i];
529 b[i] = (node_y[i + 1] - node_y[i]) / h[i];
530 } else {
531 h[i] = node_x[i + 1] - node_x[i];
532 b[i] = (node_y[i + 1] - node_y[i]) / h[i];
533 v[i] = 2. * (h[i] + h[i - 1]);
534 u[i] = 6. * (b[i] - b[i - 1]);
535 }
536 }
537
538 // Tri-diagonal matrix
539 let mut m = matrix(vec![0f64; (n - 1) * (n - 1)], n - 1, n - 1, Col);
540 for i in 0..n - 2 {
541 m[(i, i)] = v[i + 1];
542 m[(i + 1, i)] = h[i + 1];
543 m[(i, i + 1)] = h[i + 1];
544 }
545 m[(n - 2, n - 2)] = v[n - 1];
546
547 // Calculate z
548 let z_inner = m.inv() * Vec::from(&u[1..]);
549 let mut z = vec![0f64];
550 z.extend(&z_inner);
551 z.push(0f64);
552
553 // Declare empty spline
554 let mut s: Vec<Polynomial> = Vec::new();
555
556 // Main process
557 for i in 0..n {
558 // Memoization
559 let t_i = node_x[i];
560 let t_i1 = node_x[i + 1];
561 let z_i = z[i];
562 let z_i1 = z[i + 1];
563 let h_i = h[i];
564 let y_i = node_y[i];
565 let y_i1 = node_y[i + 1];
566 let temp1 = poly(vec![1f64, -t_i]);
567 let temp2 = poly(vec![1f64, -t_i1]);
568
569 let term1 = temp1.powi(3) * (z_i1 / (6f64 * h_i));
570 let term2 = temp2.powi(3) * (-z_i / (6f64 * h_i));
571 let term3 = temp1 * (y_i1 / h_i - z_i1 * h_i / 6.);
572 let term4 = temp2 * (-y_i / h_i + h_i * z_i / 6.0);
573
574 s.push(term1 + term2 + term3 + term4);
575 }
576 Ok(s)
577 }
578
579 /// Extends the spline with the given nodes.
580 ///
581 /// # Description
582 ///
583 /// The method ensures that the transition between each polynomial is smooth and that the spline
584 /// interpolation of the new nodes is calculated around `x = 0` in order to avoid that
585 /// successive spline extensions with large x values become inaccurate.
586 pub fn extend_with_nodes(&mut self, node_x: Vec<f64>, node_y: Vec<f64>) -> Result<()> {
587 let mut ext_node_x = Vec::with_capacity(node_x.len() + 1);
588 let mut ext_node_y = Vec::with_capacity(node_x.len() + 1);
589
590 let (r, polynomial) = &self.polynomials[self.polynomials.len() - 1];
591 ext_node_x.push(0.0f64);
592 ext_node_y.push(polynomial.eval(r.end));
593
594 // translating the node towards x = 0 increases accuracy tremendously
595 let tx = r.end;
596 ext_node_x.extend(node_x.into_iter().map(|x| x - tx));
597 ext_node_y.extend(node_y);
598
599 let polynomials = zip_range(
600 &ext_node_x,
601 &CubicSpline::cubic_spline(&ext_node_x, &ext_node_y)?,
602 );
603
604 self.polynomials
605 .extend(polynomials.into_iter().map(|(r, p)| {
606 (
607 Range {
608 start: r.start + tx,
609 end: r.end + tx,
610 },
611 p.translate_x(tx),
612 )
613 }));
614
615 Ok(())
616 }
617}
618
619impl Into<Vec<Polynomial>> for CubicSpline {
620 fn into(self) -> Vec<Polynomial> {
621 self.polynomials
622 .into_iter()
623 .map(|(_, polynomial)| polynomial)
624 .collect()
625 }
626}
627
628impl From<Vec<(Range<f64>, Polynomial)>> for CubicSpline {
629 fn from(polynomials: Vec<(Range<f64>, Polynomial)>) -> Self {
630 CubicSpline { polynomials }
631 }
632}
633
634impl Into<Vec<(Range<f64>, Polynomial)>> for CubicSpline {
635 fn into(self) -> Vec<(Range<f64>, Polynomial)> {
636 self.polynomials
637 }
638}
639
640impl Index<usize> for CubicSpline {
641 type Output = (Range<f64>, Polynomial);
642
643 fn index(&self, index: usize) -> &Self::Output {
644 &self.polynomials[index]
645 }
646}
647
648impl Calculus for CubicSpline {
649 fn derivative(&self) -> Self {
650 let mut polynomials: Vec<(Range<f64>, Polynomial)> = self.clone().into();
651
652 polynomials = polynomials
653 .into_iter()
654 .map(|(r, poly)| (r, poly.derivative()))
655 .collect();
656
657 Self::from(polynomials)
658 }
659
660 fn integral(&self) -> Self {
661 let mut polynomials: Vec<(Range<f64>, Polynomial)> = self.clone().into();
662
663 polynomials = polynomials
664 .into_iter()
665 .map(|(r, poly)| (r, poly.integral()))
666 .collect();
667
668 Self::from(polynomials)
669 }
670
671 fn integrate<T: Into<f64> + Copy>(&self, interval: (T, T)) -> f64 {
672 let (a, b) = interval;
673 let a = a.into();
674 let b = b.into();
675
676 let mut s = 0f64;
677 for (r, p) in self.polynomials.iter() {
678 if r.start > b {
679 break;
680 } else if r.end < a {
681 continue;
682 } else {
683 // r.start <= b, r.end >= a
684 let x = r.start.max(a);
685 let y = r.end.min(b);
686 s += p.integrate((x, y));
687 }
688 }
689 s
690 }
691}
692
693// =============================================================================
694// Cubic Hermite Spline
695// =============================================================================
696#[derive(Debug, Clone, Default)]
697#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
698#[cfg_attr(
699 feature = "rkyv",
700 derive(rkyv::Archive, rkyv::Deserialize, rkyv::Serialize)
701)]
702pub struct CubicHermiteSpline {
703 polynomials: Vec<(Range<f64>, Polynomial)>,
704}
705
706impl PolynomialSpline for CubicHermiteSpline {
707 fn get_ranged_polynomials(&self) -> &Vec<(Range<f64>, Polynomial)> {
708 &self.polynomials
709 }
710}
711
712impl CubicHermiteSpline {
713 pub fn from_nodes_with_slopes(node_x: &[f64], node_y: &[f64], m: &[f64]) -> Result<Self> {
714 let n = node_x.len();
715 if n < 3 {
716 bail!(NotEnoughNodes);
717 }
718 if n != node_y.len() {
719 bail!(NotEqualNodes);
720 }
721 if n != m.len() {
722 bail!(NotEqualSlopes);
723 }
724
725 let mut r = vec![Range::default(); node_x.len() - 1];
726 let mut u = vec![Polynomial::default(); node_x.len() - 1];
727
728 for i in 0..node_x.len() - 1 {
729 let a_i = node_y[i];
730 let b_i = m[i];
731 let dx = node_x[i + 1] - node_x[i];
732 let dy = node_y[i + 1] - node_y[i];
733 let c_i = (3f64 * dy / dx - 2f64 * m[i] - m[i + 1]) / dx;
734 let d_i = (m[i] + m[i + 1] - 2f64 * dy / dx) / dx.powi(2);
735
736 let p = Polynomial::new(vec![1f64, -node_x[i]]);
737
738 r[i] = Range {
739 start: node_x[i],
740 end: node_x[i + 1],
741 };
742 u[i] = p.powi(3) * d_i + p.powi(2) * c_i + p.clone() * b_i;
743 u[i].coef[3] += a_i;
744 }
745
746 Ok(CubicHermiteSpline {
747 polynomials: r.into_iter().zip(u).collect(),
748 })
749 }
750
751 pub fn from_nodes(node_x: &[f64], node_y: &[f64], slope_method: SlopeMethod) -> Result<Self> {
752 match slope_method {
753 SlopeMethod::Akima => CubicHermiteSpline::from_nodes_with_slopes(
754 node_x,
755 node_y,
756 &akima_slopes(node_x, node_y)?,
757 ),
758 SlopeMethod::Quadratic => CubicHermiteSpline::from_nodes_with_slopes(
759 node_x,
760 node_y,
761 &quadratic_slopes(node_x, node_y)?,
762 ),
763 }
764 }
765}
766
767impl Into<Vec<Polynomial>> for CubicHermiteSpline {
768 fn into(self) -> Vec<Polynomial> {
769 self.polynomials
770 .into_iter()
771 .map(|(_, polynomial)| polynomial)
772 .collect()
773 }
774}
775
776impl From<Vec<(Range<f64>, Polynomial)>> for CubicHermiteSpline {
777 fn from(polynomials: Vec<(Range<f64>, Polynomial)>) -> Self {
778 CubicHermiteSpline { polynomials }
779 }
780}
781
782impl Into<Vec<(Range<f64>, Polynomial)>> for CubicHermiteSpline {
783 fn into(self) -> Vec<(Range<f64>, Polynomial)> {
784 self.polynomials
785 }
786}
787
788impl Index<usize> for CubicHermiteSpline {
789 type Output = (Range<f64>, Polynomial);
790
791 fn index(&self, index: usize) -> &Self::Output {
792 &self.polynomials[index]
793 }
794}
795
796impl Calculus for CubicHermiteSpline {
797 fn derivative(&self) -> Self {
798 let mut polynomials: Vec<(Range<f64>, Polynomial)> = self.clone().into();
799
800 polynomials = polynomials
801 .into_iter()
802 .map(|(r, poly)| (r, poly.derivative()))
803 .collect();
804
805 Self::from(polynomials)
806 }
807
808 fn integral(&self) -> Self {
809 let mut polynomials: Vec<(Range<f64>, Polynomial)> = self.clone().into();
810
811 polynomials = polynomials
812 .into_iter()
813 .map(|(r, poly)| (r, poly.integral()))
814 .collect();
815
816 Self::from(polynomials)
817 }
818
819 fn integrate<T: Into<f64> + Copy>(&self, interval: (T, T)) -> f64 {
820 let (a, b) = interval;
821 let a = a.into();
822 let b = b.into();
823
824 let mut s = 0f64;
825 for (r, p) in self.polynomials.iter() {
826 if r.start > b {
827 break;
828 } else if r.end < a {
829 continue;
830 } else {
831 // r.start <= b, r.end >= a
832 let x = r.start.max(a);
833 let y = r.end.min(b);
834 s += p.integrate((x, y));
835 }
836 }
837 s
838 }
839}
840
841// =============================================================================
842// Estimate Slopes
843// =============================================================================
844#[derive(Debug, Copy, Clone, PartialEq, Eq)]
845pub enum SlopeMethod {
846 Akima,
847 Quadratic,
848}
849
850fn akima_slopes(x: &[f64], y: &[f64]) -> Result<Vec<f64>> {
851 if x.len() < 3 {
852 bail!(NotEnoughNodes);
853 }
854
855 let mut m = vec![0f64; x.len()];
856 let mut s = vec![0f64; x.len() + 3]; // -2, -1, 0, ..., x.len()-1, x.len()
857
858 let l_i = lagrange_polynomial(x[0..3].to_vec(), y[0..3].to_vec());
859 let l_f = lagrange_polynomial(x[x.len() - 3..].to_vec(), y[y.len() - 3..].to_vec());
860
861 let x_i = x[0] - (x[1] - x[0]) / 10f64;
862 let x_ii = x_i - (x[1] - x[0]) / 10f64;
863 let x_f = x[x.len() - 1] + (x[x.len() - 1] - x[x.len() - 2]) / 10f64;
864 let x_ff = x_f + (x[x.len() - 1] - x[x.len() - 2]) / 10f64;
865
866 let y_i = l_i.eval(x_i);
867 let y_ii = l_i.eval(x_ii);
868 let y_f = l_f.eval(x_f);
869 let y_ff = l_f.eval(x_ff);
870
871 let new_x = concat(&concat(&[x_ii, x_i], x), &[x_f, x_ff]);
872 let new_y = concat(&concat(&[y_ii, y_i], y), &[y_f, y_ff]);
873
874 for i in 0..new_x.len() - 1 {
875 let dx = new_x[i + 1] - new_x[i];
876 if dx == 0f64 {
877 bail!(RedundantNodeX);
878 }
879 s[i] = (new_y[i + 1] - new_y[i]) / dx;
880 }
881
882 for (i, m_i) in m.iter_mut().enumerate() {
883 let j = i + 2;
884 let ds_f = (s[j + 1] - s[j]).abs();
885 let ds_i = (s[j - 1] - s[j - 2]).abs();
886
887 *m_i = if ds_f == 0f64 && ds_i == 0f64 {
888 (s[j - 1] + s[j]) / 2f64
889 } else {
890 (ds_f * s[j - 1] + ds_i * s[j]) / (ds_f + ds_i)
891 };
892 }
893 Ok(m)
894}
895
896fn quadratic_slopes(x: &[f64], y: &[f64]) -> Result<Vec<f64>> {
897 if x.len() < 3 {
898 bail!(NotEnoughNodes);
899 }
900 let mut m = vec![0f64; x.len()];
901 let q_i = lagrange_polynomial(x[0..3].to_vec(), y[0..3].to_vec());
902 let q_f = lagrange_polynomial(x[x.len() - 3..].to_vec(), y[y.len() - 3..].to_vec());
903
904 m[0] = q_i.derivative().eval(x[0]);
905 m[x.len() - 1] = q_f.derivative().eval(x[x.len() - 1]);
906
907 for i in 1..x.len() - 1 {
908 let q = lagrange_polynomial(x[i - 1..i + 2].to_vec(), y[i - 1..i + 2].to_vec());
909 m[i] = q.derivative().eval(x[i]);
910 }
911
912 Ok(m)
913}
914
915// =============================================================================
916// B-Spline
917// =============================================================================
918/// Unit Cubic Basis Function
919///
920/// # Description
921/// Unit cubic basis function from Irwin-Hall distribution (n=4).
922/// For general interval, we substitute t = 4 * (x - a) / (b - a).
923///
924/// # Reference
925/// [Wikipedia](https://en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution#Special_cases)
926#[derive(Debug, Copy, Clone)]
927#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
928#[cfg_attr(
929 feature = "rkyv",
930 derive(rkyv::Archive, rkyv::Deserialize, rkyv::Serialize)
931)]
932pub struct UnitCubicBasis {
933 pub x_min: f64,
934 pub x_max: f64,
935 pub scale: f64,
936}
937
938impl UnitCubicBasis {
939 pub fn new(x_min: f64, x_max: f64, scale: f64) -> Self {
940 Self {
941 x_min,
942 x_max,
943 scale,
944 }
945 }
946
947 pub fn eval(&self, x: f64) -> f64 {
948 let t = 4f64 * (x - self.x_min) / (self.x_max - self.x_min);
949
950 let result = if (0f64..1f64).contains(&t) {
951 t.powi(3) / 6f64
952 } else if (1f64..2f64).contains(&t) {
953 (-3f64 * t.powi(3) + 12f64 * t.powi(2) - 12f64 * t + 4f64) / 6f64
954 } else if (2f64..3f64).contains(&t) {
955 (3f64 * t.powi(3) - 24f64 * t.powi(2) + 60f64 * t - 44f64) / 6f64
956 } else if (3f64..4f64).contains(&t) {
957 (4f64 - t).powi(3) / 6f64
958 } else {
959 0f64
960 };
961
962 self.scale * result
963 }
964
965 pub fn eval_vec(&self, x: &[f64]) -> Vec<f64> {
966 x.iter().map(|x| self.eval(*x)).collect()
967 }
968}
969
970/// Uniform Cubic B-Spline basis functions
971///
972/// # Example
973///
974/// ```rust
975/// use peroxide::fuga::*;
976/// use core::ops::Range;
977///
978/// # #[allow(unused_variables)]
979/// fn main() -> anyhow::Result<()> {
980/// let cubic_b_spline = CubicBSplineBases::from_interval((0f64, 1f64), 5);
981/// let x = linspace(0f64, 1f64, 1000);
982/// let y = cubic_b_spline.eval_vec(&x);
983///
984/// # #[cfg(feature = "plot")] {
985/// let mut plt = Plot2D::new();
986/// plt.set_domain(x.clone());
987///
988/// for basis in &cubic_b_spline.bases {
989/// plt.insert_image(basis.eval_vec(&x));
990/// }
991///
992/// plt
993/// .insert_image(y)
994/// .set_xlabel(r"$x$")
995/// .set_ylabel(r"$y$")
996/// .set_style(PlotStyle::Nature)
997/// .tight_layout()
998/// .set_dpi(600)
999/// .set_path("example_data/cubic_b_spline.png")
1000/// .savefig()?;
1001/// # }
1002/// Ok(())
1003/// }
1004/// ```
1005#[derive(Debug, Clone)]
1006#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
1007#[cfg_attr(
1008 feature = "rkyv",
1009 derive(rkyv::Archive, rkyv::Deserialize, rkyv::Serialize)
1010)]
1011pub struct CubicBSplineBases {
1012 pub ranges: Vec<Range<f64>>,
1013 pub bases: Vec<UnitCubicBasis>,
1014}
1015
1016impl CubicBSplineBases {
1017 /// Create new Cubic B-Spline basis functions
1018 pub fn new(ranges: Vec<Range<f64>>, bases: Vec<UnitCubicBasis>) -> Self {
1019 Self { ranges, bases }
1020 }
1021
1022 /// Create new Cubic B-Spline basis functions for `[a, b]`
1023 pub fn from_interval((a, b): (f64, f64), num_bases: usize) -> Self {
1024 let nodes = linspace(a, b, num_bases + 4);
1025 let (ranges, bases) = nodes
1026 .iter()
1027 .zip(nodes.iter().skip(4))
1028 .map(|(a, b)| {
1029 (
1030 Range { start: *a, end: *b },
1031 UnitCubicBasis::new(*a, *b, 1f64),
1032 )
1033 })
1034 .unzip();
1035
1036 Self::new(ranges, bases)
1037 }
1038
1039 /// Rescale all basis functions
1040 ///
1041 /// # Arguments
1042 /// - `scale_vec` - scale vector
1043 pub fn rescale(&mut self, scale_vec: &[f64]) -> Result<()> {
1044 if scale_vec.len() != self.bases.len() {
1045 bail!("The number of scales should be equal to the number of basis functions");
1046 }
1047
1048 for (basis, scale) in self.bases.iter_mut().zip(scale_vec) {
1049 basis.scale = *scale;
1050 }
1051
1052 Ok(())
1053 }
1054
1055 pub fn eval(&self, x: f64) -> f64 {
1056 self.ranges
1057 .iter()
1058 .enumerate()
1059 .filter(|(_, range)| range.contains(&x))
1060 .fold(0f64, |acc, (i, _)| {
1061 let basis = &self.bases[i];
1062 acc + basis.eval(x)
1063 })
1064 }
1065
1066 pub fn eval_vec(&self, x: &[f64]) -> Vec<f64> {
1067 x.iter().map(|x| self.eval(*x)).collect()
1068 }
1069}
1070
1071/// B-Spline
1072///
1073/// # Description
1074/// : B-Spline is a generalization of Bezier curve.
1075///
1076/// # Caution
1077/// - Let K = the number of knots, C = the number of control points
1078/// - For open, K = C + degree + 1 (C = K - degree - 1)
1079/// - For clamped, K + 2 * degree = C + degree + 1 (C = K + degree - 1)
1080///
1081/// # Example
1082/// ```
1083/// use peroxide::fuga::*;
1084/// use anyhow::Result;
1085///
1086/// fn main() -> Result<()> {
1087/// let degree = 3;
1088/// let knots = vec![0f64, 0.5, 1f64];
1089/// let control_points = vec![
1090/// vec![0f64, 0f64],
1091/// vec![0.3, 0.5],
1092/// vec![0.5, 0.7],
1093/// vec![0.7, 0.3],
1094/// vec![1f64, 1f64],
1095/// ];
1096/// let spline = BSpline::clamped(degree, knots, control_points)?;
1097///
1098/// let t = linspace(0f64, 1f64, 100);
1099/// let (x, y): (Vec<f64>, Vec<f64>) = spline.eval_vec(&t).into_iter().unzip();
1100/// assert_eq!(x.len(), 100);
1101/// assert_eq!(y.len(), 100);
1102///
1103/// Ok(())
1104/// }
1105/// ```
1106#[derive(Debug, Clone, Default)]
1107#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
1108#[cfg_attr(
1109 feature = "rkyv",
1110 derive(rkyv::Archive, rkyv::Deserialize, rkyv::Serialize)
1111)]
1112pub struct BSpline {
1113 pub degree: usize,
1114 pub knots: Vec<f64>,
1115 pub control_points: Vec<Vec<f64>>,
1116}
1117
1118impl BSpline {
1119 /// Create new open B-Spline
1120 ///
1121 /// # Arguments
1122 /// - `degree` - Degree of B-Spline
1123 /// - `knots` - Knots (length = K)
1124 /// - `control_points` - Control points (length = C)
1125 ///
1126 /// # Caution
1127 /// - K = C + degree + 1 (C = K - degree - 1)
1128 pub fn open(degree: usize, knots: Vec<f64>, control_points: Vec<Vec<f64>>) -> Result<Self> {
1129 if knots.len() != control_points.len() + degree + 1 {
1130 bail!("The number of knots ({}) should be equal to the number of control points ({}) + degree ({}) + 1", knots.len(), control_points.len(), degree);
1131 }
1132
1133 Ok(Self {
1134 degree,
1135 knots,
1136 control_points,
1137 })
1138 }
1139
1140 /// Create new clamped B-Spline
1141 ///
1142 /// # Arguments
1143 /// - `degree` - Degree of B-Spline
1144 /// - `knots` - Knots (length = K)
1145 /// - `control_points` - Control points (length = C)
1146 ///
1147 /// # Caution
1148 /// - K + 2 * degree = C + degree + 1 (C = K + degree - 1)
1149 pub fn clamped(degree: usize, knots: Vec<f64>, control_points: Vec<Vec<f64>>) -> Result<Self> {
1150 let mut knots = knots.clone();
1151
1152 if knots.len() != control_points.len() + 1 - degree {
1153 bail!("For clamped, the number of knots ({}) should be equal to the number of control points ({}) + 1 - degree ({})", knots.len(), control_points.len(), degree);
1154 }
1155 for _ in 0..degree {
1156 knots.insert(0, knots[0]);
1157 knots.push(knots[knots.len() - 1]);
1158 }
1159 if knots.len() != control_points.len() + degree + 1 {
1160 bail!("The number of knots ({}) should be equal to the number of control points ({}) + degree ({}) + 1", knots.len(), control_points.len(), degree);
1161 }
1162
1163 Ok(Self {
1164 degree,
1165 knots,
1166 control_points,
1167 })
1168 }
1169
1170 /// Obtain basis function via Cox-de Boor algorithm
1171 #[allow(non_snake_case)]
1172 pub fn cox_de_boor(&self, t: f64, i: usize) -> f64 {
1173 let p = self.degree;
1174 let mut B = vec![vec![0f64; p + 1]; p + 1];
1175
1176 // Initialize 0th degree basis functions
1177 for (j, B_j) in B.iter_mut().enumerate() {
1178 if (self.knots[i + j] <= t && t < self.knots[i + j + 1])
1179 || (i + j == self.knots.len() - (p + 2) && t == self.knots[i + j + 1])
1180 {
1181 B_j[0] = 1f64;
1182 } else {
1183 B_j[0] = 0f64;
1184 }
1185 }
1186
1187 // Compute the basis functions for higher degree
1188 for k in 1..=p {
1189 for j in 0..=(p - k) {
1190 let a = if self.knots[i + j + k] == self.knots[i + j] {
1191 0f64
1192 } else {
1193 (t - self.knots[i + j]) / (self.knots[i + j + k] - self.knots[i + j])
1194 };
1195
1196 let b = if self.knots[i + j + k + 1] == self.knots[i + j + 1] {
1197 0f64
1198 } else {
1199 (self.knots[i + j + k + 1] - t)
1200 / (self.knots[i + j + k + 1] - self.knots[i + j + 1])
1201 };
1202
1203 B[j][k] = a * B[j][k - 1] + b * B[j + 1][k - 1];
1204 }
1205 }
1206
1207 B[0][p]
1208 }
1209
1210 /// Returns the derivative of the B-spline.
1211 /// The derivative of a B-spline of degree p is a B-spline of degree p-1.
1212 pub fn derivative(&self) -> Result<Self> {
1213 if self.degree == 0 {
1214 bail!("Cannot take the derivative of a B-spline with degree 0.");
1215 }
1216
1217 let p = self.degree;
1218 let c = self.control_points.len();
1219
1220 // The new control points (Q_i) are calculated from the old ones (P_i)
1221 // Q_i = p * (P_{i+1} - P_i) / (t_{i+p+1} - t_{i+1})
1222 let mut new_control_points = Vec::with_capacity(c - 1);
1223 for i in 0..c - 1 {
1224 let p_i = &self.control_points[i];
1225 let p_i1 = &self.control_points[i + 1];
1226 let denominator = self.knots[i + p + 1] - self.knots[i + 1];
1227 if denominator == 0.0 {
1228 new_control_points.push(vec![0.0; p_i.len()]);
1229 } else {
1230 let factor = p as f64 / denominator;
1231 let q_i = p_i1
1232 .iter()
1233 .zip(p_i)
1234 .map(|(a, b)| factor * (a - b))
1235 .collect();
1236 new_control_points.push(q_i);
1237 }
1238 }
1239
1240 // The new knot vector is the old one with the first and last knots removed.
1241 let new_knots = self.knots[1..self.knots.len() - 1].to_vec();
1242 let new_degree = self.degree - 1;
1243
1244 BSpline::open(new_degree, new_knots, new_control_points)
1245 }
1246
1247 /// Returns the integral of the B-spline.
1248 /// The integral of a B-spline of degree p is a B-spline of degree p+1.
1249 /// The first control point of the integral is set to zero (integration constant).
1250 pub fn integral(&self) -> Result<Self> {
1251 let p = self.degree;
1252 let c = self.control_points.len();
1253 let dim = self.control_points[0].len();
1254
1255 // The new knot vector is the old one with the first and last knots duplicated.
1256 let mut new_knots = self.knots.clone();
1257 new_knots.insert(0, self.knots[0]);
1258 new_knots.push(self.knots[self.knots.len() - 1]);
1259
1260 // The new control points (R_i) are calculated recursively.
1261 // R_j = R_{j-1} + P_{j-1} * (t_{j+p} - t_{j-1}) / (p+1)
1262 let mut new_control_points = vec![vec![0.0; dim]; c + 1]; // R_0 is the zero vector
1263
1264 for i in 1..=c {
1265 let factor = (self.knots[i + p] - self.knots[i - 1]) / ((p + 1) as f64);
1266 let prev_r = new_control_points[i - 1].clone();
1267 let p_im1 = &self.control_points[i - 1];
1268 let mut r_i = Vec::with_capacity(dim);
1269 for j in 0..dim {
1270 r_i.push(prev_r[j] + factor * p_im1[j]);
1271 }
1272 new_control_points[i] = r_i;
1273 }
1274
1275 let new_degree = self.degree + 1;
1276
1277 BSpline::open(new_degree, new_knots, new_control_points)
1278 }
1279}
1280
1281impl Spline<(f64, f64)> for BSpline {
1282 #[allow(non_snake_case)]
1283 fn eval(&self, t: f64) -> (f64, f64) {
1284 let n = self.control_points.len();
1285
1286 let mut x = 0f64;
1287 let mut y = 0f64;
1288 for i in 0..n {
1289 let B = self.cox_de_boor(t, i);
1290 x += B * self.control_points[i][0];
1291 y += B * self.control_points[i][1];
1292 }
1293
1294 (x, y)
1295 }
1296}