1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 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.12.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"><</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">></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"><</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">></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"><</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">></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"><</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">></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"><</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">></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> ∗ <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> ∗ <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₀(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"><</span><span class="number">0</span><span class="special">>(</span><span class="identifier">x</span><span class="special">);</span> 75<span class="comment">// B₁(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"><</span><span class="number">1</span><span class="special">>(</span><span class="identifier">x</span><span class="special">);</span> 77<span class="comment">// First derivative of B₃:</span> 78<span class="identifier">yp</span> <span class="special">=</span> <span class="identifier">cardinal_b_spline_prime</span><span class="special"><</span><span class="number">3</span><span class="special">>(</span><span class="identifier">x</span><span class="special">);</span> 79<span class="comment">// Second derivative of B₃:</span> 80<span class="identifier">ypp</span> <span class="special">=</span> <span class="identifier">cardinal_b_spline_double_prime</span><span class="special"><</span><span class="number">3</span><span class="special">>(</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> > 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"><</span><span class="number">1</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">18.8</span> <span class="identifier">ns</span> 144<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">2</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">25.3</span> <span class="identifier">ns</span> 145<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">3</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">29.3</span> <span class="identifier">ns</span> 146<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">4</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">33.8</span> <span class="identifier">ns</span> 147<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">5</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">36.7</span> <span class="identifier">ns</span> 148<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">6</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">39.1</span> <span class="identifier">ns</span> 149<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">7</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">43.6</span> <span class="identifier">ns</span> 150<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">8</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">62.8</span> <span class="identifier">ns</span> 151<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">9</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">70.2</span> <span class="identifier">ns</span> 152<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">10</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">83.8</span> <span class="identifier">ns</span> 153<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">11</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">94.3</span> <span class="identifier">ns</span> 154<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">12</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">108</span> <span class="identifier">ns</span> 155<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">13</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">122</span> <span class="identifier">ns</span> 156<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">14</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">138</span> <span class="identifier">ns</span> 157<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">15</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">155</span> <span class="identifier">ns</span> 158<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">16</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">170</span> <span class="identifier">ns</span> 159<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">17</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">192</span> <span class="identifier">ns</span> 160<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">18</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">174</span> <span class="identifier">ns</span> 161<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">19</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">180</span> <span class="identifier">ns</span> 162<span class="identifier">CardinalBSpline</span><span class="special"><</span><span class="number">20</span><span class="special">,</span> <span class="keyword">double</span><span class="special">></span> <span class="number">194</span> <span class="identifier">ns</span> 163<span class="identifier">UniformReal</span><span class="special"><</span><span class="keyword">double</span><span class="special">></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 © 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å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> 226
