...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#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 B0 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, B1 := B0 ∗ B0, and Bn+1 := Bn ∗ B0. For example, B1(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);
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 B1 is what Schoenberg
calls M2. 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].
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.
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.
Some representative ULP plots are shown below. The error grows linearly with n, as expected from Cox equation 10.5.