libs/math/doc/html/math_toolkit/sf_poly/cardinal_b_splines.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Cardinal B-splines</title>
   5 <link rel="stylesheet" href="../../math.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
   7 <link rel="home" href="../../index.html" title="Math Toolkit 2.11.0">
   8 <link rel="up" href="../sf_poly.html" title="Polynomials">
   9 <link rel="prev" href="sph_harm.html" title="Spherical Harmonics">
  10 <link rel="next" href="gegenbauer.html" title="Gegenbauer Polynomials">
  11 </head>
  12 <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
  13 <table cellpadding="2" width="100%"><tr>
  14 <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../boost.png"></td>
  15 <td align="center"><a href="../../../../../../index.html">Home</a></td>
  16 <td align="center"><a href="../../../../../../libs/libraries.htm">Libraries</a></td>
  17 <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
  18 <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
  19 <td align="center"><a href="../../../../../../more/index.htm">More</a></td>
  20 </tr></table>
  21 <hr>
  22 <div class="spirit-nav">
  23 <a accesskey="p" href="sph_harm.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../sf_poly.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="gegenbauer.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
  24 </div>
  25 <div class="section">
  26 <div class="titlepage"><div><div><h3 class="title">
  27 <a name="math_toolkit.sf_poly.cardinal_b_splines"></a><a class="link" href="cardinal_b_splines.html" title="Cardinal B-splines">Cardinal B-splines</a>
  28 </h3></div></div></div>
  29 <h5>
  30 <a name="math_toolkit.sf_poly.cardinal_b_splines.h0"></a>
  31         <span class="phrase"><a name="math_toolkit.sf_poly.cardinal_b_splines.synopsis"></a></span><a class="link" href="cardinal_b_splines.html#math_toolkit.sf_poly.cardinal_b_splines.synopsis">Synopsis</a>
  32       </h5>
  33 <pre class="programlisting"><span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">cardinal_b_spline</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</span>
  34 </pre>
  35 <pre class="programlisting"><span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">{</span> <span class="keyword">namespace</span> <span class="identifier">math</span><span class="special">{</span>
  36
  37 <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">unsigned</span> <span class="identifier">n</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">&gt;</span>
  38 <span class="keyword">auto</span> <span class="identifier">cardinal_b_spline</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">);</span>
  39
  40 <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">unsigned</span> <span class="identifier">n</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">&gt;</span>
  41 <span class="keyword">auto</span> <span class="identifier">cardinal_b_spline_prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">);</span>
  42
  43 <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">unsigned</span> <span class="identifier">n</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">&gt;</span>
  44 <span class="keyword">auto</span> <span class="identifier">cardinal_b_spline_double_prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">);</span>
  45
  46 <span class="keyword">template</span><span class="special">&lt;</span><span class="keyword">unsigned</span> <span class="identifier">n</span><span class="special">,</span> <span class="keyword">typename</span> <span class="identifier">Real</span><span class="special">&gt;</span>
  47 <span class="identifier">Real</span> <span class="identifier">forward_cardinal_b_spline</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">);</span>
  48
  49 <span class="special">}}</span> <span class="comment">// namespaces</span>
  50 </pre>
  51 <p>
  52         Cardinal <span class="emphasis"><em>B</em></span>-splines are a family of compactly supported
  53         functions useful for the smooth interpolation of tables.
  54       </p>
  55 <p>
  56         The first <span class="emphasis"><em>B</em></span>-spline <span class="emphasis"><em>B</em></span><sub>0</sub> is simply
  57         a box function: It takes the value one inside the interval [-1/2, 1/2], and
  58         is zero elsewhere. <span class="emphasis"><em>B</em></span>-splines of higher smoothness are
  59         constructed by iterative convolution, namely, <span class="emphasis"><em>B</em></span><sub>1</sub> :=
  60         <span class="emphasis"><em>B</em></span><sub>0</sub> &#8727; <span class="emphasis"><em>B</em></span><sub>0</sub>, and <span class="emphasis"><em>B</em></span><sub><span class="emphasis"><em>n</em></span>+1</sub> :=
  61         <span class="emphasis"><em>B</em></span><sub><span class="emphasis"><em>n</em></span> </sub> &#8727; <span class="emphasis"><em>B</em></span><sub>0</sub>.
  62         For example, <span class="emphasis"><em>B</em></span><sub>1</sub>(<span class="emphasis"><em>x</em></span>) = 1 - |<span class="emphasis"><em>x</em></span>|
  63         for <span class="emphasis"><em>x</em></span> in [-1,1], and zero elsewhere, so it is a hat
  64         function.
  65       </p>
  66 <p>
  67         A basic usage is as follows:
  68       </p>
  69 <pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cardinal_b_spline</span><span class="special">;</span>
  70 <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cardinal_b_spline_prime</span><span class="special">;</span>
  71 <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cardinal_b_spline_double_prime</span><span class="special">;</span>
  72 <span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">0.5</span><span class="special">;</span>
  73 <span class="comment">// B&#8320;(x), the box function:</span>
  74 <span class="keyword">double</span> <span class="identifier">y</span> <span class="special">=</span> <span class="identifier">cardinal_b_spline</span><span class="special">&lt;</span><span class="number">0</span><span class="special">&gt;(</span><span class="identifier">x</span><span class="special">);</span>
  75 <span class="comment">// B&#8321;(x), the hat function:</span>
  76 <span class="identifier">y</span> <span class="special">=</span> <span class="identifier">cardinal_b_spline</span><span class="special">&lt;</span><span class="number">1</span><span class="special">&gt;(</span><span class="identifier">x</span><span class="special">);</span>
  77 <span class="comment">// First derivative of B&#8323;:</span>
  78 <span class="identifier">yp</span> <span class="special">=</span> <span class="identifier">cardinal_b_spline_prime</span><span class="special">&lt;</span><span class="number">3</span><span class="special">&gt;(</span><span class="identifier">x</span><span class="special">);</span>
  79 <span class="comment">// Second derivative of B&#8323;:</span>
  80 <span class="identifier">ypp</span> <span class="special">=</span> <span class="identifier">cardinal_b_spline_double_prime</span><span class="special">&lt;</span><span class="number">3</span><span class="special">&gt;(</span><span class="identifier">x</span><span class="special">);</span>
  81 </pre>
  82 <p>
  83         <span class="inlinemediaobject"><object type="image/svg+xml" data="../../../graphs/central_b_splines.svg"></object></span> <span class="inlinemediaobject"><object type="image/svg+xml" data="../../../graphs/central_b_spline_derivatives.svg"></object></span>
  84         <span class="inlinemediaobject"><object type="image/svg+xml" data="../../../graphs/central_b_spline_second_derivatives.svg"></object></span>
  85       </p>
  86 <h4>
  87 <a name="math_toolkit.sf_poly.cardinal_b_splines.h1"></a>
  88         <span class="phrase"><a name="math_toolkit.sf_poly.cardinal_b_splines.caveats"></a></span><a class="link" href="cardinal_b_splines.html#math_toolkit.sf_poly.cardinal_b_splines.caveats">Caveats</a>
  89       </h4>
  90 <p>
  91         Numerous notational conventions for <span class="emphasis"><em>B</em></span>-splines exist.
  92         Whereas Boost.Math (following Kress) zero indexes the <span class="emphasis"><em>B</em></span>-splines,
  93         other authors (such as Schoenberg and Massopust) use 1-based indexing. So
  94         (for example) Boost.Math's hat function <span class="emphasis"><em>B</em></span><sub>1</sub> is what Schoenberg
  95         calls <span class="emphasis"><em>M</em></span><sub>2</sub>. Mathematica, like Boost, uses the zero-indexing
  96         convention for its <code class="computeroutput"><span class="identifier">BSplineCurve</span></code>.
  97       </p>
  98 <p>
  99         Even the support of the splines is not agreed upon. Mathematica starts the
 100         support of the splines at zero and rescales the independent variable so that
 101         the support of every member is [0, 1]. Massopust as well as Unser puts the
 102         support of the <span class="emphasis"><em>B</em></span>-splines at [0, <span class="emphasis"><em>n</em></span>],
 103         whereas Kress centers them at zero. Schoenberg distinguishes between the
 104         two cases, called the splines starting at zero forward splines, and the ones
 105         symmetric about zero <span class="emphasis"><em>central</em></span>.
 106       </p>
 107 <p>
 108         The <span class="emphasis"><em>B</em></span>-splines of Boost.Math are central, with support
 109         support [-(<span class="emphasis"><em>n</em></span>+1)/2, (<span class="emphasis"><em>n</em></span>+1)/2]. If
 110         necessary, the forward splines can be evaluated by using <code class="computeroutput"><span class="identifier">forward_cardinal_b_spline</span></code>,
 111         whose support is [0, <span class="emphasis"><em>n</em></span>+1].
 112       </p>
 113 <h4>
 114 <a name="math_toolkit.sf_poly.cardinal_b_splines.h2"></a>
 115         <span class="phrase"><a name="math_toolkit.sf_poly.cardinal_b_splines.implementation"></a></span><a class="link" href="cardinal_b_splines.html#math_toolkit.sf_poly.cardinal_b_splines.implementation">Implementation</a>
 116       </h4>
 117 <p>
 118         The implementation follows Maurice Cox' 1972 paper 'The Numerical Evaluation
 119         of B-splines', and uses the triangular array of Algorithm 6.1 of the reference
 120         rather than the rhombohedral array of Algorithm 6.2. Benchmarks revealed
 121         that the time to calculate the indexes of the rhombohedral array exceed the
 122         time to simply add zeroes together, except for <span class="emphasis"><em>n</em></span> &gt;
 123         18. Since few people use <span class="emphasis"><em>B</em></span> splines of degree 18, the
 124         triangular array is used.
 125       </p>
 126 <h4>
 127 <a name="math_toolkit.sf_poly.cardinal_b_splines.h3"></a>
 128         <span class="phrase"><a name="math_toolkit.sf_poly.cardinal_b_splines.performance"></a></span><a class="link" href="cardinal_b_splines.html#math_toolkit.sf_poly.cardinal_b_splines.performance">Performance</a>
 129       </h4>
 130 <p>
 131         Double precision timing on a consumer x86 laptop is shown below:
 132       </p>
 133 <pre class="programlisting"><span class="identifier">Run</span> <span class="identifier">on</span> <span class="special">(</span><span class="number">16</span> <span class="identifier">X</span> <span class="number">4300</span> <span class="identifier">MHz</span> <span class="identifier">CPU</span> <span class="identifier">s</span><span class="special">)</span>
 134 <span class="identifier">CPU</span> <span class="identifier">Caches</span><span class="special">:</span>
 135   <span class="identifier">L1</span> <span class="identifier">Data</span> <span class="number">32</span><span class="identifier">K</span> <span class="special">(</span><span class="identifier">x8</span><span class="special">)</span>
 136   <span class="identifier">L1</span> <span class="identifier">Instruction</span> <span class="number">32</span><span class="identifier">K</span> <span class="special">(</span><span class="identifier">x8</span><span class="special">)</span>
 137   <span class="identifier">L2</span> <span class="identifier">Unified</span> <span class="number">1024</span><span class="identifier">K</span> <span class="special">(</span><span class="identifier">x8</span><span class="special">)</span>
 138   <span class="identifier">L3</span> <span class="identifier">Unified</span> <span class="number">11264</span><span class="identifier">K</span> <span class="special">(</span><span class="identifier">x1</span><span class="special">)</span>
 139 <span class="identifier">Load</span> <span class="identifier">Average</span><span class="special">:</span> <span class="number">0.21</span><span class="special">,</span> <span class="number">0.33</span><span class="special">,</span> <span class="number">0.29</span>
 140 <span class="special">-----------------------------------------</span>
 141 <span class="identifier">Benchmark</span>                            <span class="identifier">Time</span>
 142 <span class="special">-----------------------------------------</span>
 143 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">1</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">18.8</span> <span class="identifier">ns</span>
 144 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">2</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">25.3</span> <span class="identifier">ns</span>
 145 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">3</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">29.3</span> <span class="identifier">ns</span>
 146 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">4</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">33.8</span> <span class="identifier">ns</span>
 147 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">5</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">36.7</span> <span class="identifier">ns</span>
 148 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">6</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">39.1</span> <span class="identifier">ns</span>
 149 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">7</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">43.6</span> <span class="identifier">ns</span>
 150 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">8</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">62.8</span> <span class="identifier">ns</span>
 151 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">9</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">70.2</span> <span class="identifier">ns</span>
 152 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">10</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>       <span class="number">83.8</span> <span class="identifier">ns</span>
 153 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">11</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>       <span class="number">94.3</span> <span class="identifier">ns</span>
 154 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">12</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">108</span> <span class="identifier">ns</span>
 155 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">13</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">122</span> <span class="identifier">ns</span>
 156 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">14</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">138</span> <span class="identifier">ns</span>
 157 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">15</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">155</span> <span class="identifier">ns</span>
 158 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">16</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">170</span> <span class="identifier">ns</span>
 159 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">17</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">192</span> <span class="identifier">ns</span>
 160 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">18</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">174</span> <span class="identifier">ns</span>
 161 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">19</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">180</span> <span class="identifier">ns</span>
 162 <span class="identifier">CardinalBSpline</span><span class="special">&lt;</span><span class="number">20</span><span class="special">,</span> <span class="keyword">double</span><span class="special">&gt;</span>        <span class="number">194</span> <span class="identifier">ns</span>
 163 <span class="identifier">UniformReal</span><span class="special">&lt;</span><span class="keyword">double</span><span class="special">&gt;</span>               <span class="number">11.5</span> <span class="identifier">ns</span>
 164 </pre>
 165 <p>
 166         A uniformly distributed random number within the support of the spline is
 167         generated for the argument, so subtracting 11.5 ns from these gives a good
 168         idea of the performance of the calls.
 169       </p>
 170 <h4>
 171 <a name="math_toolkit.sf_poly.cardinal_b_splines.h4"></a>
 172         <span class="phrase"><a name="math_toolkit.sf_poly.cardinal_b_splines.accuracy"></a></span><a class="link" href="cardinal_b_splines.html#math_toolkit.sf_poly.cardinal_b_splines.accuracy">Accuracy</a>
 173       </h4>
 174 <p>
 175         Some representative ULP plots are shown below. The error grows linearly with
 176         <span class="emphasis"><em>n</em></span>, as expected from Cox equation 10.5.
 177       </p>
 178 <p>
 179         <span class="inlinemediaobject"><object type="image/svg+xml" data="../../../graphs/b_spline_ulp_3.svg"></object></span> <span class="inlinemediaobject"><object type="image/svg+xml" data="../../../graphs/b_spline_ulp_5.svg"></object></span>
 180         <span class="inlinemediaobject"><object type="image/svg+xml" data="../../../graphs/b_spline_ulp_9.svg"></object></span>
 181       </p>
 182 <h4>
 183 <a name="math_toolkit.sf_poly.cardinal_b_splines.h5"></a>
 184         <span class="phrase"><a name="math_toolkit.sf_poly.cardinal_b_splines.references"></a></span><a class="link" href="cardinal_b_splines.html#math_toolkit.sf_poly.cardinal_b_splines.references">References</a>
 185       </h4>
 186 <div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
 187 <li class="listitem">
 188             I.J. Schoenberg, <span class="emphasis"><em>Cardinal Spline Interpolation</em></span>,
 189             SIAM Volume 12, 1973
 190           </li>
 191 <li class="listitem">
 192             Rainer Kress, <span class="emphasis"><em>Numerical Analysis</em></span>, Springer, 1998
 193           </li>
 194 <li class="listitem">
 195             Peter Massopust, <span class="emphasis"><em>On Some Generalizations of B-splines</em></span>,
 196             arxiv preprint, 2019
 197           </li>
 198 <li class="listitem">
 199             Michael Unser and Thierry Blu, <span class="emphasis"><em>Fractional Splines and Wavelets</em></span>,
 200             SIAM Review 2000, Volume 42, No. 1
 201           </li>
 202 <li class="listitem">
 203             Cox, Maurice G. <span class="emphasis"><em>The numerical evaluation of B-splines.</em></span>,
 204             IMA Journal of Applied Mathematics 10.2 (1972): 134-149.
 205           </li>
 206 </ul></div>
 207 </div>
 208 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 209 <td align="left"></td>
 210 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
 211       Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
 212       Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
 213       R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
 214       Daryle Walker and Xiaogang Zhang<p>
 215         Distributed under the Boost Software License, Version 1.0. (See accompanying
 216         file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
 217       </p>
 218 </div></td>
 219 </tr></table>
 220 <hr>
 221 <div class="spirit-nav">
 222 <a accesskey="p" href="sph_harm.html"><img src="../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../sf_poly.html"><img src="../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="gegenbauer.html"><img src="../../../../../../doc/src/images/next.png" alt="Next"></a>
 223 </div>
 224 </body>
 225 </html>