140 lines
5.4 KiB
Plaintext
140 lines
5.4 KiB
Plaintext
[/
|
|
Copyright 2019, Nick Thompson
|
|
Distributed under the Boost Software License, Version 1.0.
|
|
(See accompanying file LICENSE_1_0.txt or copy at
|
|
http://www.boost.org/LICENSE_1_0.txt).
|
|
]
|
|
|
|
[section:cardinal_b_splines Cardinal B-splines]
|
|
|
|
[h4 Synopsis]
|
|
|
|
``
|
|
#include <boost/math/special_functions/cardinal_b_spline.hpp>
|
|
``
|
|
|
|
namespace boost{ namespace math{
|
|
|
|
template<unsigned n, typename Real>
|
|
auto cardinal_b_spline(Real x);
|
|
|
|
template<unsigned n, typename Real>
|
|
auto cardinal_b_spline_prime(Real x);
|
|
|
|
template<unsigned n, typename Real>
|
|
auto cardinal_b_spline_double_prime(Real x);
|
|
|
|
template<unsigned n, typename Real>
|
|
Real forward_cardinal_b_spline(Real x);
|
|
|
|
}} // namespaces
|
|
|
|
Cardinal /B/-splines are a family of compactly supported functions useful for the smooth interpolation of tables.
|
|
|
|
The first /B/-spline /B/[sub 0] is simply a box function:
|
|
It takes the value one inside the interval \[-1/2, 1/2\], and is zero elsewhere.
|
|
/B/-splines of higher smoothness are constructed by iterative convolution, namely, /B/[sub 1] := /B/[sub 0] \u2217 /B/[sub 0],
|
|
and /B/[sub /n/+1] := /B/[sub /n/ ] \u2217 /B/[sub 0].
|
|
For example, /B/[sub 1](/x/) = 1 - |/x/| for /x/ in \[-1,1\], and zero elsewhere, so it is a hat function.
|
|
|
|
A basic usage is as follows:
|
|
|
|
using boost::math::cardinal_b_spline;
|
|
using boost::math::cardinal_b_spline_prime;
|
|
using boost::math::cardinal_b_spline_double_prime;
|
|
double x = 0.5;
|
|
// B₀(x), the box function:
|
|
double y = cardinal_b_spline<0>(x);
|
|
// B₁(x), the hat function:
|
|
y = cardinal_b_spline<1>(x);
|
|
// First derivative of B₃:
|
|
yp = cardinal_b_spline_prime<3>(x);
|
|
// Second derivative of B₃:
|
|
ypp = cardinal_b_spline_double_prime<3>(x);
|
|
|
|
|
|
[$../graphs/central_b_splines.svg]
|
|
[$../graphs/central_b_spline_derivatives.svg]
|
|
[$../graphs/central_b_spline_second_derivatives.svg]
|
|
|
|
[h3 Caveats]
|
|
|
|
Numerous notational conventions for /B/-splines exist.
|
|
Whereas Boost.Math (following Kress) zero indexes the /B/-splines,
|
|
other authors (such as Schoenberg and Massopust) use 1-based indexing.
|
|
So (for example) Boost.Math's hat function /B/[sub 1] is what Schoenberg calls /M/[sub 2].
|
|
Mathematica, like Boost, uses the zero-indexing convention for its `BSplineCurve`.
|
|
|
|
Even the support of the splines is not agreed upon.
|
|
Mathematica starts the support of the splines at zero and rescales the independent variable so that the support of every member is \[0, 1\].
|
|
Massopust as well as Unser puts the support of the /B/-splines at \[0, /n/\], whereas Kress centers them at zero.
|
|
Schoenberg distinguishes between the two cases, called the splines starting at zero forward splines, and the ones symmetric about zero /central/.
|
|
|
|
The /B/-splines of Boost.Math are central, with support support \[-(/n/+1)\/2, (/n/+1)\/2\].
|
|
If necessary, the forward splines can be evaluated by using `forward_cardinal_b_spline`, whose support is \[0, /n/+1\].
|
|
|
|
[h3 Implementation]
|
|
|
|
The implementation follows Maurice Cox' 1972 paper 'The Numerical Evaluation of B-splines', and uses the triangular array of Algorithm 6.1 of the reference rather than the rhombohedral array of Algorithm 6.2.
|
|
Benchmarks revealed that the time to calculate the indexes of the rhombohedral array exceed the time to simply add zeroes together, except for /n/ > 18.
|
|
Since few people use /B/ splines of degree 18, the triangular array is used.
|
|
|
|
[h3 Performance]
|
|
|
|
Double precision timing on a consumer x86 laptop is shown below:
|
|
|
|
``
|
|
Run on (16 X 4300 MHz CPU s)
|
|
CPU Caches:
|
|
L1 Data 32K (x8)
|
|
L1 Instruction 32K (x8)
|
|
L2 Unified 1024K (x8)
|
|
L3 Unified 11264K (x1)
|
|
Load Average: 0.21, 0.33, 0.29
|
|
-----------------------------------------
|
|
Benchmark Time
|
|
-----------------------------------------
|
|
CardinalBSpline<1, double> 18.8 ns
|
|
CardinalBSpline<2, double> 25.3 ns
|
|
CardinalBSpline<3, double> 29.3 ns
|
|
CardinalBSpline<4, double> 33.8 ns
|
|
CardinalBSpline<5, double> 36.7 ns
|
|
CardinalBSpline<6, double> 39.1 ns
|
|
CardinalBSpline<7, double> 43.6 ns
|
|
CardinalBSpline<8, double> 62.8 ns
|
|
CardinalBSpline<9, double> 70.2 ns
|
|
CardinalBSpline<10, double> 83.8 ns
|
|
CardinalBSpline<11, double> 94.3 ns
|
|
CardinalBSpline<12, double> 108 ns
|
|
CardinalBSpline<13, double> 122 ns
|
|
CardinalBSpline<14, double> 138 ns
|
|
CardinalBSpline<15, double> 155 ns
|
|
CardinalBSpline<16, double> 170 ns
|
|
CardinalBSpline<17, double> 192 ns
|
|
CardinalBSpline<18, double> 174 ns
|
|
CardinalBSpline<19, double> 180 ns
|
|
CardinalBSpline<20, double> 194 ns
|
|
UniformReal<double> 11.5 ns
|
|
``
|
|
|
|
A uniformly distributed random number within the support of the spline is generated for the argument, so subtracting 11.5 ns from these gives a good idea of the performance of the calls.
|
|
|
|
[h3 Accuracy]
|
|
|
|
Some representative ULP plots are shown below.
|
|
The error grows linearly with /n/, as expected from Cox equation 10.5.
|
|
|
|
[$../graphs/b_spline_ulp_3.svg]
|
|
[$../graphs/b_spline_ulp_5.svg]
|
|
[$../graphs/b_spline_ulp_9.svg]
|
|
|
|
[h3 References]
|
|
|
|
* I.J. Schoenberg, ['Cardinal Spline Interpolation], SIAM Volume 12, 1973
|
|
* Rainer Kress, ['Numerical Analysis], Springer, 1998
|
|
* Peter Massopust, ['On Some Generalizations of B-splines], arxiv preprint, 2019
|
|
* Michael Unser and Thierry Blu, ['Fractional Splines and Wavelets], SIAM Review 2000, Volume 42, No. 1
|
|
* Cox, Maurice G. ['The numerical evaluation of B-splines.], IMA Journal of Applied Mathematics 10.2 (1972): 134-149.
|
|
|
|
[endsect]
|