123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139 |
- [/
- 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]
|