libs/math/doc/sf/cardinal_b_splines.qbk

   1 [/
   2   Copyright 2019, Nick Thompson
   3   Distributed under the Boost Software License, Version 1.0.
   4   (See accompanying file LICENSE_1_0.txt or copy at
   5   http://www.boost.org/LICENSE_1_0.txt).
   6 ]
   7
   8 [section:cardinal_b_splines Cardinal B-splines]
   9
  10 [h4 Synopsis]
  11
  12 ``
  13 #include <boost/math/special_functions/cardinal_b_spline.hpp>
  14 ``
  15
  16    namespace boost{ namespace math{
  17
  18    template<unsigned n, typename Real>
  19    auto cardinal_b_spline(Real x);
  20
  21    template<unsigned n, typename Real>
  22    auto cardinal_b_spline_prime(Real x);
  23
  24    template<unsigned n, typename Real>
  25    auto cardinal_b_spline_double_prime(Real x);
  26
  27    template<unsigned n, typename Real>
  28    Real forward_cardinal_b_spline(Real x);
  29
  30    }} // namespaces
  31
  32 Cardinal /B/-splines are a family of compactly supported functions useful for the smooth interpolation of tables.
  33
  34 The first /B/-spline /B/[sub 0] is simply a box function:
  35 It takes the value one inside the interval \[-1/2, 1/2\], and is zero elsewhere.
  36 /B/-splines of higher smoothness are constructed by iterative convolution, namely, /B/[sub 1] := /B/[sub 0] \u2217 /B/[sub 0],
  37 and /B/[sub /n/+1] := /B/[sub /n/ ] \u2217 /B/[sub 0].
  38 For example, /B/[sub 1](/x/) = 1 - |/x/| for /x/ in \[-1,1\], and zero elsewhere, so it is a hat function.
  39
  40 A basic usage is as follows:
  41
  42     using boost::math::cardinal_b_spline;
  43     using boost::math::cardinal_b_spline_prime;
  44     using boost::math::cardinal_b_spline_double_prime;
  45     double x = 0.5;
  46     // B₀(x), the box function:
  47     double y = cardinal_b_spline<0>(x);
  48     // B₁(x), the hat function:
  49     y = cardinal_b_spline<1>(x);
  50     // First derivative of B₃:
  51     yp = cardinal_b_spline_prime<3>(x);
  52     // Second derivative of B₃:
  53     ypp = cardinal_b_spline_double_prime<3>(x);
  54
  55
  56 [$../graphs/central_b_splines.svg]
  57 [$../graphs/central_b_spline_derivatives.svg]
  58 [$../graphs/central_b_spline_second_derivatives.svg]
  59
  60 [h3 Caveats]
  61
  62 Numerous notational conventions for /B/-splines exist.
  63 Whereas Boost.Math (following Kress) zero indexes the /B/-splines,
  64 other authors (such as Schoenberg and Massopust) use 1-based indexing.
  65 So (for example) Boost.Math's hat function /B/[sub 1] is what Schoenberg calls /M/[sub 2].
  66 Mathematica, like Boost, uses the zero-indexing convention for its `BSplineCurve`.
  67
  68 Even the support of the splines is not agreed upon.
  69 Mathematica starts the support of the splines at zero and rescales the independent variable so that the support of every member is \[0, 1\].
  70 Massopust as well as Unser puts the support of the /B/-splines at \[0, /n/\], whereas Kress centers them at zero.
  71 Schoenberg distinguishes between the two cases, called the splines starting at zero forward splines, and the ones symmetric about zero /central/.
  72
  73 The /B/-splines of Boost.Math are central, with support support \[-(/n/+1)\/2, (/n/+1)\/2\].
  74 If necessary, the forward splines can be evaluated by using `forward_cardinal_b_spline`, whose support is \[0, /n/+1\].
  75
  76 [h3 Implementation]
  77
  78 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.
  79 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.
  80 Since few people use /B/ splines of degree 18, the triangular array is used.
  81
  82 [h3 Performance]
  83
  84 Double precision timing on a consumer x86 laptop is shown below:
  85
  86 ``
  87 Run on (16 X 4300 MHz CPU s)
  88 CPU Caches:
  89   L1 Data 32K (x8)
  90   L1 Instruction 32K (x8)
  91   L2 Unified 1024K (x8)
  92   L3 Unified 11264K (x1)
  93 Load Average: 0.21, 0.33, 0.29
  94 -----------------------------------------
  95 Benchmark                            Time
  96 -----------------------------------------
  97 CardinalBSpline<1, double>        18.8 ns
  98 CardinalBSpline<2, double>        25.3 ns
  99 CardinalBSpline<3, double>        29.3 ns
 100 CardinalBSpline<4, double>        33.8 ns
 101 CardinalBSpline<5, double>        36.7 ns
 102 CardinalBSpline<6, double>        39.1 ns
 103 CardinalBSpline<7, double>        43.6 ns
 104 CardinalBSpline<8, double>        62.8 ns
 105 CardinalBSpline<9, double>        70.2 ns
 106 CardinalBSpline<10, double>       83.8 ns
 107 CardinalBSpline<11, double>       94.3 ns
 108 CardinalBSpline<12, double>        108 ns
 109 CardinalBSpline<13, double>        122 ns
 110 CardinalBSpline<14, double>        138 ns
 111 CardinalBSpline<15, double>        155 ns
 112 CardinalBSpline<16, double>        170 ns
 113 CardinalBSpline<17, double>        192 ns
 114 CardinalBSpline<18, double>        174 ns
 115 CardinalBSpline<19, double>        180 ns
 116 CardinalBSpline<20, double>        194 ns
 117 UniformReal<double>               11.5 ns
 118 ``
 119
 120 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.
 121
 122 [h3 Accuracy]
 123
 124 Some representative ULP plots are shown below.
 125 The error grows linearly with /n/, as expected from Cox equation 10.5.
 126
 127 [$../graphs/b_spline_ulp_3.svg]
 128 [$../graphs/b_spline_ulp_5.svg]
 129 [$../graphs/b_spline_ulp_9.svg]
 130
 131 [h3 References]
 132
 133 * I.J. Schoenberg, ['Cardinal Spline Interpolation], SIAM Volume 12, 1973
 134 * Rainer Kress, ['Numerical Analysis], Springer, 1998
 135 * Peter Massopust, ['On Some Generalizations of B-splines], arxiv preprint, 2019
 136 * Michael Unser and Thierry Blu, ['Fractional Splines and Wavelets], SIAM Review 2000, Volume 42, No. 1
 137 * Cox, Maurice G. ['The numerical evaluation of B-splines.], IMA Journal of Applied Mathematics 10.2 (1972): 134-149.
 138
 139 [endsect]