Smoothing Splines#

Overview#

The smoothing splines module in SplineOps provides tools to fit fractional smoothing splines to noisy data. Think of these splines as flexible low-pass filters whose sharpness can be tuned continuously, making them effective for signals and images that exhibit repeating, self-similar patterns.

You will find:

  • Exact 1D routine: Works on a 1D array and returns the mathematically exact fractional-spline result.

  • Isotropic N-D routine: Extends the idea to 2D pictures or 3D volumes through one FFT; internally it uses a Butterworth low-pass filter.

  • Fast linear shortcut: A lightweight forward/backward IIR filter that implements the first-order (piecewise-linear) smoothing spline in linear time—handy for real-time streams.

  • Extra helpers to generate test data (fractional Brownian motion) and to compute spline autocorrelations.

Core Idea in One Dimension#

Given noisy samples \(y[k]\) at integer positions, we look for a smooth curve \(s(t)\) that minimises

\[\sum_{k}\lvert y[k]-s(k)\rvert^{2} \;+\; \lambda\,\lVert\partial^{\gamma}s\rVert_{L^{2}}^{2},\]

where

  • the first term measures closeness to the data,

  • the second term penalises roughness,

  • \(\lambda\) balances the two,

  • \(\partial^{\gamma}\) is a fractional derivative: the optimal solution is a fractional spline of degree \(2\gamma - 1\).

Taking the discrete Fourier transform (DFT) of both sides turns the problem into a simple, frequency-by-frequency scaling

\[S(\omega) \;=\; H(\omega)\,Y(\omega),\]

where \(Y(\omega)\) is the DFT of the data and \(S(\omega)\) the DFT of the solution. \(H(\omega)\) has been computed and has a closed form [1]. A convenient Butterworth-like approximation of \(H(\omega)\) is \(H_{2\gamma}(\omega)\), with

\[H_{2\gamma}(\omega)=\frac{1}{1+\lambda\,|\omega|^{2\gamma}}.\]

Note

Here \(\omega\) denotes the (angular) frequency variable. This form is equivalent to the standard Butterworth parameterization \(1/(1+|\omega/\omega_0|^{2\gamma})\) with \(\omega_0=\lambda^{-1/(2\gamma)}\).

The practical recipe is therefore

  1. FFT the data,

  2. multiply by \(H_{2\gamma}(\omega)\),

  3. inverse FFT to obtain the smoothed samples.

A full derivation of this result can be found in [1], [2] and [3].

The next figure, from 1D Fractional Brownian Motion, illustrates these ideas on a noisy fractional Brownian motion: the original process, noisy measurements, and the smoothed spline estimates at the original and oversampled grids.

../_images/sphx_glr_01_1d_fractional_brownian_motion_001.png

Core Idea in Higher Dimensions#

For a 2D image or a 3D volume we replace the one-dimensional fractional derivative with the fractional Laplacian \((-\Delta)^{\gamma/2}\). The variational cost therefore becomes

\[\sum_{\mathbf k}\bigl|\,y[\mathbf k]-s(\mathbf k)\bigr|^{2} \;+\; \lambda\,\bigl\lVert(-\Delta)^{\gamma/2}s\bigr\rVert_{L^{2}}^{2}.\]

In the Fourier domain the Laplacian turns into \(\|\boldsymbol\omega\|^{2}\). A good approximation of the optimal filter is the radial counterpart of the 1D Butterworth filter:

\[S(\boldsymbol\omega) \;=\; \frac{1}{1+\lambda\,\lVert\boldsymbol\omega\rVert^{2\gamma}}\, Y(\boldsymbol\omega).\]

The practical algorithm is identical to the 1D case:

  1. Run an n-dimensional FFT to obtain \(Y(\boldsymbol\omega)\).

  2. Multiply by the gain above.

  3. Apply the inverse FFT to get the smoothed image or volume.

Because the filter is applied element-wise in the frequency domain, the computation still needs just one forward FFT and one inverse FFT, whatever the data dimension.

Fast Recursive Linear Smoother#

When you only need the first-order case (\(\gamma = 1\)), the corresponding discrete frequency response simplifies to

\[H(\omega)=\frac{1}{1+4\lambda\,\sin^{2}(\omega/2)}.\]

This symmetric all-pole response can be factorized into a causal and an anti-causal first-order filter with pole \(z_1\in(0,1)\):

\[H(\omega)=\frac{(1-z_1)^2}{1-2z_1\cos\omega+z_1^2}, \qquad z_1=\frac{\sqrt{1+4\lambda}-1}{\sqrt{1+4\lambda}+1}.\]

With that number in hand the algorithm is

  1. Causal pass (forward), with steady-state initialization:

    \[c[0]=\frac{y[0]}{1-z_1},\qquad c[k]=y[k]+z_1\,c[k-1].\]

  2. Anti-causal pass (backward), with steady-state initialization:

    \[s[K-1]=\frac{c[K-1]}{1-z_1},\qquad s[k]=c[k]+z_1\,s[k+1].\]

  3. Normalization:

    \[s[k]\leftarrow (1-z_1)^2\,s[k].\]

The two passes yield the same zero-phase result as applying \(H(\omega)\) in the Fourier domain, but at a cost that is strictly linear in the number of samples and with constant memory. A detailed derivation (including boundary handling for finite-length signals) appears in [4], Sections II-B and II-D.

Note

The cubic smoothing spline case requires a higher-order recursion (complex-conjugate poles) and does not reduce to a single real pole; see [4], Section IV-B.

Choosing the Parameters#

  • gamma: controls how steeply the filter rolls off (larger values ⇒ steeper transition). Typical range: \(0.5 \le \gamma \le 3\).

  • lambda: moves the cut-off frequency (small values keep more detail, large values smooth harder). For most images, \(10^{-3} \le \lambda \le 10^{-1}\) is a good starting interval.

Smoothing Splines Examples#

References#