TrueBSplineInterpolation

Uniform B-Spline Interpolation passes through all data points.

★ 1 on GitHub

View on GitHub

Anyone may press Interesting, signed in or not. It helps others find this project.

Member rating

No ratings yet.

README

TrueBSplineInterpolation

English | 日本語

A method for generating control points (green) for a uniform B-spline curve (blue) that passes through all data points (red). This repository contains the full mathematical derivation of the interpolating filter, together with a Delphi/FireMonkey demo application that implements it.

利用ライブラリ

LUX ：The foundation math library of LUXOPHIA.
LUX.Draw ：A 2D drawing library.

1. B-spline Curve Is Not Interpolation

The uniform B-spline curve $\mathrm{BSC}_N(x)$ is a method for generating smooth and highly controllable piecewise polynomials by arranging uniform B-spline basis functions $\mathrm{CBS}_N(x)$ with different weights (control point positions) at equal intervals. By considering the sequence of control points $\mathrm{CPs}(i)$ as a discrete signal, this process can be interpreted as a form of signal processing using the B-spline basis function as a filter.

The degree of the B-spline basis function can be adjusted freely. A zero-degree (order 1) B-spline corresponds to a piecewise constant approximation using square (rectangular) basis functions, while a first-degree (order 2) B-spline corresponds to a piecewise linear approximation using triangular basis functions.

.svg)

CBS[ N_, X_ ] = BSplineBasis[ N-1, (X+N/2)/N ]; @ Wolfram Language

N: order = 1 + degree

math

\begin{aligned}
{\rm CBS}_1(x) &= \begin{dcases}
 1 & | x| \leq \frac{1}{2} \\
 0 & \text{otherwise} \\
\end{dcases} &\qquad \text{(1.1)} \\

{\rm CBS}_2(x) &= \begin{dcases}
 1-| x|  & | x| \leq 1 \\
 0 & \text{otherwise} \\
\end{dcases} &\qquad \text{(1.2)} \\

{\rm CBS}_3(x) &= \begin{dcases}
 \frac{1}{4} \left(3-4 | x| ^2\right) & | x| <\frac{1}{2} \\
 \frac{1}{8} (2 | x| -3)^2 & \frac{1}{2}\leq | x| \leq \frac{3}{2} \\
 0 & \text{otherwise} \\
\end{dcases} &\qquad \text{(1.3)} \\

{\rm CBS}_4(x) &= \begin{dcases}
 \frac{1}{6} \left(3 | x| ^3-6 | x| ^2+4\right) & | x| <1 \\
 -\frac{1}{6} (| x| -2)^3 & 1\leq | x| \leq 2 \\
 0 & \text{otherwise} \\
\end{dcases} &\qquad \text{(1.4)} \\

{\rm CBS}_5(x) &= \begin{dcases} 
 \frac{1}{192} \left(48 | x| ^4-120 | x| ^2+115\right) & | x| <\frac{1}{2} \\
 \frac{1}{96} \left(-16 | x| ^4+80 | x| ^3-120 | x| ^2+20 | x| +55\right) & \frac{1}{2}\leq | x| <\frac{3}{2} \\
 \frac{1}{384} (2 | x| -5)^4 & \frac{3}{2}\leq | x| \leq \frac{5}{2} \\
 0 & \text{otherwise} \\
\end{dcases} &\qquad \text{(1.5)} \\

{\rm CBS}_6(x) &= \begin{dcases}
 \frac{1}{60} \left(-5 | x| ^5+15 | x| ^4-30 | x| ^2+33\right) & | x| <1 \\
 \frac{1}{120} \left(5 | x| ^5-45 | x| ^4+150 | x| ^3-210 | x| ^2+75 | x| +51\right) & 1\leq | x| <2 \\
 -\frac{1}{120} (| x| -3)^5 & 2\leq | x| \leq 3 \\
 0 & \text{otherwise} \\
\end{dcases} &\qquad \text{(1.6)} \\
\end{aligned}

</br>

The frequency characteristics of the B-spline basis function are defined as the power of the $\mathrm{sinc}$ function. Due to its strong low-pass characteristics, filtering with the B-spline basis function significantly attenuates the high-frequency components inherent in the control point sequence. As a result, while the generated curve is very smooth because of the pronounced smoothing effect, it is an approximate curve that does not pass through any of the control points.

math

\widehat{{\rm CBS}}_N(\omega) = {\left[\frac{2 \sin \left(\frac{\omega }{2}\right)}{\omega }\right]}^N \qquad \text{(1.7)}

</br>

.svg)

2. B-spline Basis as a Discrete Filter

Although the B-spline basis functions are continuous filters, they can be regarded as discrete filters when focusing solely on the values at the control points. In other words, whether the curve passes through the control points depends on the frequency characteristics of the discrete B-spline basis function $\mathrm{DBS}_N(i)$.

.svg)

The frequency characteristics of the discrete B-spline basis function are periodic. For orders above 3, it exhibits low-pass characteristics with the greatest attenuation at the Nyquist frequency ($\pi$). Conversely, discrete B-spline bases of order 2 or lower do not attenuate the high-frequency components, so the generated curve always passes through the control points.

math

\begin{aligned}
\widehat{{\rm DBS}}_1(\omega)&=1 &\qquad \text{(2.1)} \\
\widehat{{\rm DBS}}_2(\omega)&=1 &\qquad \text{(2.2)} \\
\widehat{{\rm DBS}}_3(\omega)&=\frac{1}{4} (\cos (\omega )+3) &\qquad \text{(2.3)} \\
\widehat{{\rm DBS}}_4(\omega)&=\frac{1}{3} (\cos (\omega )+2) &\qquad \text{(2.4)} \\
\widehat{{\rm DBS}}_5(\omega)&=\frac{1}{192} (76 \cos (\omega )+\cos (2 \omega )+115) &\qquad \text{(2.5)} \\
\widehat{{\rm DBS}}_6(\omega)&=\frac{1}{60} (26 \cos (\omega )+\cos (2 \omega )+33) &\qquad \text{(2.6)} \\
\end{aligned}

</br>

.svg)

3. Discrete High-Enhancement Filter

To make the curve pass through the control points, it is necessary to restore the original high-frequency components of the control point sequence. This can be achieved by designing a discrete high-enhancement filter $\mathrm{DHE}_N(i)$ that cancels out the low-pass characteristics of the discrete B-spline basis function. In other words, its frequency characteristics are the reciprocal of that of the discrete B-spline basis function.

math

\begin{aligned}
\widehat{{\rm DHE}}_1(\omega)&=1 &\qquad \text{(3.1)} \\
\widehat{{\rm DHE}}_2(\omega)&=1 &\qquad \text{(3.2)} \\
\widehat{{\rm DHE}}_3(\omega)&=\frac{4}{\cos (\omega )+3} &\qquad \text{(3.3)} \\
\widehat{{\rm DHE}}_4(\omega)&=\frac{3}{\cos (\omega )+2} &\qquad \text{(3.4)} \\
\widehat{{\rm DHE}}_5(\omega)&=\frac{192}{76 \cos (\omega )+\cos (2 \omega )+115} &\qquad \text{(3.5)} \\
\widehat{{\rm DHE}}_6(\omega)&=\frac{60}{26 \cos (\omega )+\cos (2 \omega )+33} &\qquad \text{(3.6)} \\
\end{aligned}

</br>

.svg)

By performing an inverse Fourier transform on the frequency characteristics $\widehat{\mathrm{DHE}}_N(\omega)$, we can derive the impulse response in the time domain, i.e., the filter function $\mathrm{DHE}_N(i)$.

math

{\rm DHE}_N(i) = \frac{1}{2\pi} \int_{-\pi}^\pi \widehat{{\rm DHE}}_N(\omega) e^{j i \omega} \, d\omega \qquad \text{(3.7)}

Since $\widehat{\mathrm{DHE}}_N(\omega)$ is an even function, this is equivalent to the inverse Fourier cosine transform. Furthermore, because it is a periodic function, its impulse response is necessarily a discrete signal.

math

{\rm DHE}_N(i) = \frac{1}{\pi} \int_{0}^\pi \widehat{{\rm DHE}}_N(\omega) \cos(i\omega) \, d\omega \qquad \text{(3.8)}

By introducing the following integral formula (3.9) [2]:

math

\int_0^\pi \frac{\cos(n\theta)}{a + b \cos \theta} \, d\theta = \frac{\pi}{\sqrt{a^2 - b^2}} \left( \frac{\sqrt{a^2 - b^2}-a}{b} \right)^n \quad |b| < |a|, \, 0 \leq n \qquad \text{(3.9)}

We can analytically derive $\mathrm{DHE}_N(i)$ for $N = 3, 4$:

math

\begin{aligned}
{\rm DHE}_3(i) &= \frac{1}{\pi} \int_{0}^\pi \widehat{{\rm DHE}}_3(\omega) \cos(i\omega) \, d\omega = \frac{4}{\pi} \int_{0}^\pi \frac{\cos(i\omega)}{3 + \cos \omega} \, d\omega \\
&= \frac{4}{\pi} \frac{\pi}{\sqrt{3^2 - 1}} \left( \frac{\sqrt{3^2 - 1}-3}{1} \right)^n
= \sqrt{2} \left( 2\sqrt{2}-3 \right)^n &\qquad \text{(3.10)} \\
{\rm DHE}_4(i) &= \frac{1}{\pi} \int_{0}^\pi \widehat{{\rm DHE}}_4(\omega) \cos(i\omega) \, d\omega = \frac{3}{\pi} \int_{0}^\pi \frac{\cos(i\omega)}{2 + \cos \omega} \, d\omega \\
&= \frac{3}{\pi} \frac{\pi}{\sqrt{2^2 - 1}} \left( \frac{\sqrt{2^2 - 1}-2}{1} \right)^n
= \sqrt{3} \left( \sqrt{3}-2 \right)^n &\qquad \text{(3.11)}
\end{aligned}

The discrete high-enhancement filter is an infinite impulse response (IIR) filter. Therefore, in practical implementations, it is necessary to truncate the coefficient sequence to a finite number of taps.

$`{\rm DHE}_3(i)`$	$`{\rm DHE}_4(i)`$
$`{\rm DHE}_3(\pm 0)=\sqrt{2}`$	$`{\rm DHE}_4(\pm 0)=\sqrt{3}`$
$`{\rm DHE}_3(\pm 1)=4-3 \sqrt{2}`$	$`{\rm DHE}_4(\pm 1)=3-2 \sqrt{3}`$
$`{\rm DHE}_3(\pm 2)=17 \sqrt{2}-24`$	$`{\rm DHE}_4(\pm 2)=7 \sqrt{3}-12`$
$`{\rm DHE}_3(\pm 3)=140-99 \sqrt{2}`$	$`{\rm DHE}_4(\pm 3)=45-26 \sqrt{3}`$
$`{\rm DHE}_3(\pm 4)=577 \sqrt{2}-816`$	$`{\rm DHE}_4(\pm 4)=97 \sqrt{3}-168`$
$`{\rm DHE}_3(\pm 5)=4756-3363 \sqrt{2}`$	$`{\rm DHE}_4(\pm 5)=627-362 \sqrt{3}`$
$`{\rm DHE}_3(\pm 6)=19601 \sqrt{2}-27720`$	$`{\rm DHE}_4(\pm 6)=1351 \sqrt{3}-2340`$
$`{\rm DHE}_3(\pm 7)=161564-114243 \sqrt{2}`$	$`{\rm DHE}_4(\pm 7)=8733-5042 \sqrt{3}`$
$`{\rm DHE}_3(\pm 8)=665857 \sqrt{2}-941664`$	$`{\rm DHE}_4(\pm 8)=18817 \sqrt{3}-32592`$

The frequency characteristics of a finite impulse response (FIR) filter, i.e., an approximate discrete high-enhancement filter $\mathrm{AHE}_N(i)$, closely match that of the ideal discrete high-enhancement filter $\mathrm{DHE}_N(i)$. Empirically, using 7 taps (±3) for the third-order case and 9 taps (±4) for the fourth-order case is sufficient to ensure accuracy. (The full derivations are also available as a Mathematica notebook [1].)

.svg)

4. Architecture

The demo is a Delphi FireMonkey (FMX) application. The interpolation itself lives in the LUX.Curve.BSpline unit of the bundled LUXOPHIA library: TBSInterp convolves the data points (Poins) with the truncated high-enhancement filter HEFilter4 — the closed form of equation (3.11) — to obtain the control points (Verts), and then evaluates the cubic B-spline through TBSCurve4. MargsN is the number of margin taps (±MargsN) used when truncating the IIR filter.

pascal

[ Application ??? object composition ]

???TForm1                          ????????? ( Main.pas ) main FMX form
  ??????TGenPoins                   ????????? ( GenPoins.pas ) random data points
  ??????TSingleBSInterp             ????????? ( LUX.Curve.BSpline ) interpolator
  ??????TCurveChart                 ????????? ( CurveChartFrame.pas ) chart display

[ Curve classes ??? inheritance ]

???TCurve<_TPoin_>                 ????????? ( LUX.Curve ) abstract curve base
  ??????TBSCurve4<_TPoin_>          ????????? ( LUX.Curve.BSpline ) cubic B-spline
  ???  ??????TSingleBSCurve4
  ???  ??????TDoubleBSCurve4
  ??????TBSInterp<_TPoin_,_TCurve_> ????????? ( LUX.Curve.BSpline ) interpolation
     ??????TSingleBSInterp
     ??????TDoubleBSInterp

[ Viewer classes ??? inheritance ]

???TDrawViewer                     ????????? ( LUX.Draw.Viewer ) 2D drawing frame
  ??????TCurveChart                 ????????? ( CurveChartFrame.pas )

[ Signal flow ??? nesting shows what is fed into what ]

???Poins                           ????????? ( TCurve ) data points ( red )
  ??????HEFilter4                   ????????? ( LUX.Curve.BSpline ) HE filter (3.11)
     ??????Verts                    ????????? ( TBSInterp ) control points (green)
        ??????TBSCurve4             ????????? ( LUX.Curve.BSpline ) curve ( blue )

pascal

???TrueBSplineInterpolation/
  ??????TrueBSplineInterpolation.dpr / .dproj ????????? project (FMX, Win32/Win64)
  ??????Main.pas / Main.fmx                   ????????? main form (TForm1)
  ??????CurveChartFrame.pas / .fmx            ????????? chart frame (TCurveChart)
  ??????GenPoins.pas                          ????????? data-point generator (TGenPoins)
  ??????_LIBRARY/LUXOPHIA/
  ???  ??????LUX/                              ????????? math primitives
  ???  ???  ??????Curve/                        ????????? LUX.Curve.BSpline etc.
  ???  ??????LUX.Draw/                         ????????? 2D drawing library
  ??????--------/
     ??????True B-Spline Interpolation.nb     ????????? Mathematica notebook [1]
     ??????_README/                           ????????? figures used in this document
     ??????_SCREENSHOT/                       ????????? screenshot

5. Building

Open TrueBSplineInterpolation.dproj in RAD Studio (Delphi, FireMonkey) and run. Target platforms: Win32 / Win64 (default: Win64).

Or from the command line, after running rsvars.bat:

pascal

msbuild TrueBSplineInterpolation.dproj /t:Build /p:Config=Release /p:Platform=Win64

All dependencies (the LUXOPHIA LUX / LUX.Draw libraries) are included under _LIBRARY; no external packages are required.

6. References

True B-Spline Interpolation.nb @ Mathematica
Table of Integrals, Series, and Products, Seventh Edition, 2007.
Bierens de Haan, D., Nouvelles tables d’intégrales définies, Amsterdam, 1867. (Reprint: G. E. Stechert & Co., New York, 1939).
Bierens de Haan, D., Tables d’intégrales définies, Verhandelingen der Koninklijke Akademie van Wetenschappen, Deel IV, 1858.

7. Appendix

7.1 Derivation of Integral Formula (3.9)

math

\int_0^\pi \frac{\cos(n\theta)}{a + b \cos \theta} \, d\theta = \frac{\pi}{\sqrt{a^2 - b^2}} \left( \frac{\sqrt{a^2 - b^2}-a}{b} \right)^n \quad |b| < |a|, \, 0 \leq n \qquad \text{(3.9)}

An integral formula that has already been proven is given below. It is listed as equation (3.613-1) on page 391 of the Table of Integrals, Series, and Products (7th edition) [2] (see also the classical tables of definite integrals [3][4]).

math

\int_0^\pi \frac{\cos(n\theta)}{1 + \alpha \cos(\theta)}d\theta = \frac{\pi}{\sqrt{1 - \alpha^2}} \left( \frac{\sqrt{1 - \alpha^2} - 1}{\alpha} \right)^n, \quad \alpha^2 < 1, \, 0 \leq n \qquad \text{(7.1)}

By performing the following variable transformation, we can derive integral formula (3.9):

math

\begin{aligned}
\int_0^\pi \frac{\cos(n \theta)}{a + b \cos \theta} \, d\theta
&= \frac{1}{a} \int_0^\pi \frac{\cos(n \theta)}{1 + \dfrac{b}{a} \cos \theta} \, d\theta
= \frac{1}{a} \int_0^\pi \frac{\cos(n\theta)}{1 + \alpha \cos(\theta)}d\theta &\qquad \text{(7.2)} \\
&= \frac{1}{a} \cdot \frac{\pi}{\sqrt{1 - \left( \dfrac{b}{a} \right)^2}} \left( \frac{ \sqrt{1 - \left( \dfrac{b}{a} \right)^2 } - 1 }{ \dfrac{b}{a} } \right)^n &\qquad \text{(7.3)} \\
&= \frac{\pi}{a} \cdot \dfrac{a}{\sqrt{a^2 - b^2}} \left[ \dfrac{a}{b} \left( \dfrac{\sqrt{a^2 - b^2}}{a} - 1 \right) \right]^n &\qquad \text{(7.4)} \\
&= \frac{\pi}{\sqrt{a^2 - b^2}} \left( \frac{ \sqrt{a^2 - b^2} - a }{ b } \right)^n &\qquad \text{(7.5)} \\
\end{aligned}

The condition for the validity of the formula is as follows:

math

\alpha^2 < 1
\  \Leftrightarrow \  \left( \dfrac{b}{a} \right)^2 < 1
\  \Leftrightarrow \  b^2 < a^2
\  \Leftrightarrow \  |b| < |a| \qquad \text{(7.6)}

💖 Embarcadero Delphi

Integrated Development Environment (IDE) for Creating Native Cross-Platform Apps.

Live chat for Free Pascal, Lazarus, Delphi and every Pascal programmer.

TrueBSplineInterpolation

Member rating

README

TrueBSplineInterpolation

利用ライブラリ

1. B-spline Curve Is Not Interpolation

2. B-spline Basis as a Discrete Filter

3. Discrete High-Enhancement Filter

4. Architecture

5. Building

6. References

7. Appendix

7.1 Derivation of Integral Formula (3.9)

💖 Embarcadero Delphi

Free Download: Delphi Community Edition

Live chat for Free Pascal, Lazarus, Delphi and every Pascal programmer.

TTrueBSplineInterpolation

Member rating

README

TrueBSplineInterpolation

利用ライブラリ

1. B-spline Curve Is Not Interpolation

2. B-spline Basis as a Discrete Filter

3. Discrete High-Enhancement Filter

4. Architecture

5. Building

6. References

7. Appendix

7.1 Derivation of Integral Formula (3.9)

💖 Embarcadero Delphi

Free Download: Delphi Community Edition

TrueBSplineInterpolation