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:

  • Fractional 1D routine: Implements the periodized fractional-spline frequency response on a real, uniformly sampled signal.

  • Isotropic N-D approximation: Applies a radial Butterworth-like frequency response to periodic real data. It is useful, but is not the same discrete operator as the exact 1D fractional-spline construction.

  • 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.

For repeated arrays of the same shape and parameters, splineops.smoothing_splines.SmoothingSplinePlan constructs the radial response once and stores only the non-redundant real-FFT half spectrum:

import numpy as np
from splineops.smoothing_splines import SmoothingSplinePlan

plan = SmoothingSplinePlan(image.shape, lamb=0.02, gamma=1.5)
output = np.empty(image.shape, dtype=np.float64)
plan.apply(image, out=output)
next_output = plan(next_image)

retained_bytes reports the response storage and configuration describes its fixed shape and parameters. The plan accepts finite real arrays whose selected axes match the construction shape. Without axes, every input dimension remains spatial for backward compatibility. With explicit axes, remaining batch/channel slices are transformed independently by one batched real FFT:

plan = SmoothingSplinePlan(batch.shape[-2:], lamb=0.02, gamma=1.5)
smooth = plan(batch, axes=(-2, -1))

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 realize this response with a steady-state end treatment at a cost linear in the number of samples. It is not identical to circular FFT filtering on finite signals because their boundary assumptions differ. A detailed derivation appears in [4], Sections II-B and II-D.

Current contract#

The public smoothing functions accept finite, real-valued arrays. The 1D fractional and N-D FFT methods assume periodic extension; the N-D path and its plan use rfftn/irfftn because the public contract is real-valued. The recursive method uses steady-state causal and anticausal initialization. These routines remain experimental because the three formulations intentionally implement different operators. Current independent checks cover identity at zero regularization, the analytical response of periodic N-D cosine modes, a dense linear-system reference for the recursive boundary equations, and constant preservation. A fixed multitone fixture now evaluates the published fractional transfer function through an explicit dense DFT, independently of the implementation’s SciPy FFT. These checks support the stated contracts but do not make the three methods interchangeable.

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#