1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
//! Tensor-product cubic B-spline interpolation.
//!
//! The regular-grid implementation stores one B-spline coefficient tensor with
//! the same shape and C-style memory order as the input values. Coefficients are
//! generated by applying a one-dimensional tridiagonal solve along each axis.
//!
//! For an interior regular-grid span, four neighboring B-spline coefficients
//! define the local cubic polynomial. If `t = (x - x[i]) / h` and the footprint
//! coefficients are
//!
//! ```text
//! c0 = c[i - 1]
//! c1 = c[i]
//! c2 = c[i + 1]
//! c3 = c[i + 2]
//! ```
//!
//! then evaluation uses the cardinal cubic B-spline basis:
//!
//! ```text
//! S(t) = c0*w0(t) + c1*w1(t) + c2*w2(t) + c3*w3(t)
//!
//! w0 = (1 - 3t + 3t^2 - t^3) / 6
//! w1 = (4 - 6t^2 + 3t^3) / 6
//! w2 = (1 + 3t + 3t^2 - 3t^3) / 6
//! w3 = t^3 / 6
//! ```
//!
//! Equivalently, in power-basis form `S(t) = a0 + a1*t + a2*t^2 + a3*t^3`,
//! the polynomial coefficients are
//!
//! ```text
//! a0 = (c0 + 4c1 + c2) / 6
//! a1 = (c2 - c0) / 2
//! a2 = (c0 - 2c1 + c2) / 2
//! a3 = (-c0 + 3c1 - 3c2 + c3) / 6
//! ```
//!
//! Boundary spans and rectilinear grids use the same local-footprint idea, but
//! with boundary ghost coefficients folded into the stored coefficients or with
//! nonuniform-grid basis weights.
//!
//! For a rectilinear-grid interior span, the polynomial mapping is
//! span-dependent because the local knot spacing affects the nonuniform
//! B-spline basis. For span `s`, define
//!
//! ```text
//! h = grid[s + 1] - grid[s]
//! t = (x - grid[s]) / h
//! ```
//!
//! and let the four footprint coefficients multiply the fixed-span nonuniform
//! basis weights:
//!
//! ```text
//! S(t) = c0*b0(t) + c1*b1(t) + c2*b2(t) + c3*b3(t)
//! ```
//!
//! where `b0..b3` are computed by the rectilinear `basis_span_weights` routine
//! for that span. To convert this span to a normalized power-basis polynomial,
//! sample the local spline at four normalized positions:
//!
//! ```text
//! t0 = 0
//! t1 = 1/3
//! t2 = 2/3
//! t3 = 1
//!
//! xr = grid[s] + tr*h
//! yr = b0(xr)*c0 + b1(xr)*c1 + b2(xr)*c2 + b3(xr)*c3
//! ```
//!
//! Then, for `S(t) = a0 + a1*t + a2*t^2 + a3*t^3`,
//!
//! ```text
//! a0 = y0
//! a1 = (-11*y0 + 18*y1 - 9*y2 + 2*y3) / 2
//! a2 = (18*y0 - 45*y1 + 36*y2 - 9*y3) / 2
//! a3 = (-9*y0 + 27*y1 - 27*y2 + 9*y3) / 2
//! ```
//!
//! In physical coordinates `u = x - grid[s]`, the same polynomial is
//!
//! ```text
//! S(u) = A0 + A1*u + A2*u^2 + A3*u^3
//!
//! A0 = a0
//! A1 = a1 / h
//! A2 = a2 / h^2
//! A3 = a3 / h^3
//! ```
//!
//! Boundary spans use the same conversion after replacing the raw nonuniform
//! basis weights with the low- or high-boundary weights that fold ghost
//! coefficients into the stored coefficient footprint.
pub use MultiBsplineRectilinear;
pub use MultiBsplineRegular;
// `usize::max` is not const-stable on this MSRV, so use a local helper in
// const scratch-sizing methods.
pub const
pub