Deconvolution

Blurring is a convolution, so in the frequency domain it is a multiplication, and dividing should undo it. It does, until you add a tiny amount of noise. Blur an image with a known kernel, restore it, and find out what it takes to keep the noise under control.

Original ff
Degraded g=h∗f+ng = h \ast f + n what the camera records
Restored f^=F−1{R⋅G}\hat f = \mathcal{F}^{-1}\{R \cdot G\}
Error f^−f\hat f - f amplified 4×, gray = 0
Filters along the uu axis, v=0v = 0 blur ∣H∣|H| restoration ∣R∣|R| together ∣RH∣|R H| shaded: R=0R = 0

Blur as a multiplication

A blurred image is the sharp image convolved with a kernel hh, the point spread function, plus some noise nn. By the convolution theorem, the Fourier transform turns the convolution into a product of spectra:

g=h∗f+n⟺G(u,v)=H(u,v) F(u,v)+N(u,v)g = h \ast f + n \qquad \Longleftrightarrow \qquad G(u, v) = H(u, v) \, F(u, v) + N(u, v)

A Gaussian blur with σ\sigma pixels has the transfer function H(u,v)=e−2π2σ2(u2+v2)/N2H(u, v) = e^{-2\pi^2 \sigma^2 (u^2 + v^2) / N^2}: it keeps the low frequencies and damps the high ones. Horizontal motion blur, a box of LL pixels along xx, has a transfer function with ripples that pass close to zero near u=N/L,2N/L,…u = N / L, 2N / L, \dots. Here the kernel is known; estimating it from the blurred image as well is called blind deconvolution.

The inverse filter

If the blur is a multiplication, dividing by HH undoes it:

F^=GH=F+NH\hat F = \frac{G}{H} = F + \frac{N}{H}

Without noise this is exact. With noise, the second term is the problem. Noise such as sensor noise or rounding is spread evenly over all frequencies, while HH becomes tiny at high frequencies: 10−910^{-9} at the Nyquist frequency for σ=2\sigma = 2. Dividing multiplies the noise there by a billion. Where ∣HF∣|H F| is smaller than ∣N∣|N|, the recorded coefficient is mostly noise; the information about FF at that frequency is lost and no filter can bring it back. Even rounding the blurred image to 8 bits is enough to destroy the result.

Taming the noise

The truncated inverse (pseudo-inverse) only divides where ∣H∣≥ε|H| \ge \varepsilon and sets all other frequencies to zero. It gives up on frequencies that are too weak to recover, so the result stays somewhat blurred, and the abrupt cut-off causes ringing next to edges, like an ideal low-pass filter.

Regularization trades fidelity to the data against the size of the result. Minimizing ∥h∗f^−g∥2+K∥f^∥2\|h \ast \hat f - g\|^2 + K \|\hat f\|^2 gives, frequency by frequency,

F^=H∗∣H∣2+K G,R=HH2+K for a real H.\hat F = \frac{H^*}{|H|^2 + K} \, G, \qquad R = \frac{H}{H^2 + K} \text{ for a real } H.

Where H2≫KH^2 \gg K, R≈1/HR \approx 1 / H as before; where H2≪KH^2 \ll K, R≈H/KR \approx H / K goes smoothly to zero instead of exploding. The largest gain is 1/(2K)1 / (2\sqrt{K}). This is the Wiener filter for a constant noise-to-signal ratio KK; the full Wiener filter uses the ratio of the noise power to the signal power at each frequency, which is optimal in the mean-squared-error sense. More noise calls for a larger KK.

The blur in this demo wraps around the image borders, so that the model G=HFG = H F holds exactly. In real photos, the content beyond the border is unknown; this mismatch causes ringing at the borders unless the image is padded or its borders are tapered first.

Try this