Denoise#

Overview#

The denoise module offers practical tools for removing noise and filling in missing values with fractional smoothing splines. 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 behaves like a Butterworth low-pass filter.

  • Fast cubic shortcut: A lightweight forward/backward IIR filter that approximates the cubic case and runs in a single pass—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 (\(\gamma=1\) gives the classic cubic penalty).

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),\qquad H(\omega)=\frac{1}{1+\lambda\,|\omega|^{2\gamma}},\]

where \(Y(\omega)\) is the DFT of the data and \(S(\omega)\) the DFT of the solution. The practical recipe is therefore

  1. FFT the data,

  2. multiply by \(H(\omega)\),

  3. inverse FFT to obtain the smoothed samples.

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

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}\), so the optimal filter is the radial version of the 1D one:

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

This looks and behaves like an order \(2\gamma\) Butterworth low-pass but now works the same in every direction. 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 cubic smoother#

When you only need the cubic case (\(\gamma = 1\)) the frequency response above simplifies so much that it can be implemented with two tiny first-order filters—one run forward, the other backward. The key quantity is the pole

\[z_1 \;=\;-\frac{\lambda}{1+\sqrt{\,1+4\lambda\,}}.\]

With that number in hand the algorithm is

  1. Causal pass: start at \(k = 0\) and accumulate \(c[k] = y[k] + z_1\,c[k-1]\).

  2. Anti-causal pass: start at the last sample and run backwards \(s[k] = c[k] + z_1\,s[k+1]\).

The two passes give the same zero-phase result you would obtain from the FFT method but at a cost that is strictly linear in the number of samples and with virtually no memory footprint. A detailed derivation appears in [1], 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.

Denoise examples#

References#