Deconvolution

Usually, images acquired by a vision system suffer from degradation that can be modelled as a convolution. For example, some images present a camera shake effect (Fig. 101) or a blur due to poor focus (Fig. 102). The goal of deconvolution is to cancel the effect of a convolution.

Fig. 101 An example of motion blur (the parliament of Budapest shot by a camera).

Fig. 102 Hubble’s view of Ganymede in 1996.

The degradation phenomenon is modelled as in Fig. 103: The observed image \(y\) is degraded by the convolution with a PSF \(h\) and, possibly, by a noise \(b\) (considered to be additive).

\[ y = h*x + b \]

The deconvolution computes a deconvolved image \(\widehat{x}\) from the observation \(y\). We will consider only linear methods, thus deconvolution comes to filtering by \(g\):

\[ \widehat{x} = g*y \]

Fig. 103 Deconvolution model.

Deconvolution needs a degradation model, thus having knowledge about both \(h\) and \(b\).

• The PSF \(h\) can be estimated by observation, i.e. by finding in the image some factors to estimate \(h\). For example, a single point object in the image is \(h\). The PSF can also be estimated by experimentation by reproducing the observation conditions in a laboratory. So, the image of a pulse gives an estimate of \(h\). Finally, it is also possible to estimate the PSF from a mathematical model of the physics of the observation. Note also that some deconvolution methods estimate the PSF \(h\) at the same time as \(x\): these are called blind deconvolution methods (French: déconvolution myope).

• Models for the noise have already been presented in chapter denoising.

Inverse filter

The inverse filter is the simplest deconvolution method. Since the degradation is modelled \(y = h*x + b\), then this equation becomes in the Fourier domain:

\[ Y = HX + B \]

so we can write:

\[ X = \frac{Y-B}{H}. \]

We obtain \(x\) by calculating the inverse Fourier transform of the previous expression:

\[ x = \mathcal{F}^{-1} \left[ \frac{Y-B}{H} \right]. \]

As the noise (and therefore its spectrum \(B\)) is unknown, we can approximate the expression of \(x\) by cancelling \(B\) in the previous expression, and thus get the deconvolved image:

\[ \widehat{x} = \mathcal{F}^{-1} \left[ \frac{Y}{H} \right] \]

The result of the inverse filter applied on an image is given Fig. 104. The result is not usable, and yet the observed image is very little blurred with very little noise!

Fig. 104 Result of the inverse filter.

The catastrophic result of the inverse filter is due to the fact of having considered the noise to be zero. Indeed, according to the definition of \(\widehat{x}\) and considering \(Y = HX + B\), then:

\[ \widehat{x} = \mathcal{F}^{-1} \left[ \frac{Y}{H} \right] = \mathcal{F}^{-1} \left[ X + \frac{B}{H} \right] = x + \mathcal{F}^{-1} \left[ \frac{B}{H} \right] \]

Thus, the deconvolved image \(\widehat{x}\) corresponds to \(x\) with an additional term which is the inverse Fourier transform of \(B/H\). The PSF \(H\) is generally a low-pass filter, so the values of \(H(m,n)\) tend towards \(0\) for high frequencies \((m,n)\). Because \(H\) is in the denominator, this tends to drastically amplify the high frequencies of the noise, and then the term \(B/H\) quickly dominates \(X\). This explains the result of Fig. 104.

One solution consists in considering only the low frequencies of \(Y/H\). This is equivalent to truncating the result given by the inverse filter by cancelling the high frequencies before calculating the inverse Fourier transform. The result of the deconvolution is much more acceptable, as shown by Fig. 105, although the result is still far from perfect (there are many variations in intensity around objects, such as tree trunks)…

Fig. 105 Result of the truncated inverse filter with very small noise.

Wiener Filter

Wiener filter, denoted by \(g\) (with Fourier transform \(G\)), applies to the observation \(y\) such that:

\[ \widehat{x} = g * y \qquad\Leftrightarrow\qquad \widehat{X} = GY. \]

This filter is established in the statistical framework: the image \(x\) and the noise \(b\) are considered to be random variables. They are also assumed to be statistically independent. As a result, the observation \(y\) and the estimate \(\widehat{x}\) are also random variables.

The calculations are done in the Fourier domain for simplicity (since convolutions become multiplications). The goal of Wiener filter is to find the image \(\widehat{X} = \mathcal{F}[\widehat{x}]\) closest to \(X = \mathcal{F}[x]\), in the sense of the mean squared error \(\mathrm{MSE} = \mathbb{E}\left[(\widehat{X}-X)^2\right]\). Thereby :

\[\begin{split} \mathrm{MSE} &= \mathbb{E}\left[(\widehat{X}-X)^2\right] \\ &= \mathbb{E}\left[(GY-X)^2\right] \\ &= \mathbb{E}\left[\big(G(HX+B)-X\big)^2\right] \\ &= \mathbb{E}\left[\big((GH-I)X+GB\big)^2\right] \end{split}\]

where \(I\) is an image filled with 1. So:

\[ \mathrm{MSE} = \mathbb{E}\Big[ (GH-I)^*(GH-I)X^*X + (GH-I)^*GX^*B + G^*(GH-I)B^*X + G^*GB^*B \Big] \]

where \(\cdot^*\) denotes the conjugate of the variables. Since the expectation \(\mathbb{E}\) is linear and only \(X\) and \(B\) are random variables, we can decompose the previous expression into four terms:

\[\begin{split} \mathrm{MSE} = & (GH-I)^*(GH-I) \,\mathbb{E}\big[X^*X\big]\\ & + (GH-I)^*G \,\mathbb{E}\big[X^*B\big]\\ & + G^*(GH-I) \,\mathbb{E}\big[B^*X\big]\\ & + G^*G \,\mathbb{E}\big[B^*B\big]. \end{split}\]

Since \(X\) and \(B\) are independent, then the covariances \(\mathbb{E}\big[X^*B\big]\) and \(\mathbb{E}\big[B^*X\big] \) are zeros. Moreover, \(\mathbb{E}\big[X^*X\big]\) and \(\mathbb{E}\big[B^*B\big]\) are the power spectral densities denoted as \(S_x\) and \(S_b\) (the power spectral density is the expectation of the square of the modulus of the Fourier transform). So the mean squared error simplifies into:

\[ \mathrm{MSE} = (GH-1)^*(GH-1) S_x + G^*G S_b \]

We look for the filter \(G\) that minimizes the MSE, or equivalently, that cancels the derivative of MSE:

\[\begin{split} & \frac{\partial \mathrm{MSE}}{\partial G} = (GH-1)^*H S_x + G^* S_b = 0 \\ \Leftrightarrow\quad &G^*H^*H S_x - H S_x + G^* S_b = 0 \\ \Leftrightarrow\quad &G^* ( H^*H S_x + S_b) = H S_x \\ \Leftrightarrow\quad &G^* = \frac{H S_x}{H^*H S_x + S_b} \\ \Leftrightarrow\quad &G = \frac{H^* S_x}{H^*H S_x + S_b} \\ \Leftrightarrow\quad &G = \frac{H^* S_x}{|H|^2 S_x + S_b} \end{split}\]

Here we are, we get the expression of the Wiener filter \(G\)! 🥳 Finally, the deconvolved image is the inverse Fourier transform of \(GY\):

\[ \widehat{x} = \mathcal{F}^{-1} \Bigg[\frac{H^* S_x}{|H|^2 S_x + S_b} Y \Bigg] \]

We can consider that the power spectral densities \(S_x\) and \(S_b\) are constant (for \(S_b\), it is necessary to assume white noise). Thus, the expression of the Wiener filter can be written

\[ \widehat{x} = \mathcal{F}^{-1} \Bigg[ \frac{H^*}{|H|^2 + S_b/S_x} Y \Bigg] \]

and the term \(S_b/S_x\) is replaced by a constant \(K\), which becomes the parameter of the method, to be set by the user.

Two remarks:

• where \(H\) vanishes (typically in high frequencies), the problem of noise increase is no longer observed as with the inverse filter, since the inverse filter tends towards 0,

• moreover, if the noise in the image is zero, then \(S_b = 0\) and Wiener filter comes back to the inverse filter:

  \[ \widehat{x} = \mathcal{F}^{-1} \Bigg[ \frac{H^*}{|H|^2} Y \Bigg] = \mathcal{F}^{-1} \Bigg[ \frac{Y}{H} \Bigg] \]

The result of Wiener filter is presented Fig. 106: it is clearly much better than the inverse filter, even its truncated version!

Fig. 106 Result of Wiener filter (\(\lambda\) is set so that the estimation is the best in terms of MSE).