A deconvolution example

To introduce Jest we will discuss what probably is the Hello World program of image processing: image restoration by deconvolution. Many physical transmission systems blur their input signal. For example, atmospheric turbulences blur satellite and telescope images. The broadband wavelet in a seismic experiment blurs the reflectivity series of the Earth. Instruments themselves cause a degradation of their input which often is well described by a blur.

Mathematically, a blurred image y(t) is usually modeled as the convolution of the original image f(t) by a blurring filter b(t):

$$y(t) = b(t) \ast x(t) = \int b(t-t') x(t') dt'.$$

If the image is digitized, convolution can be expressed discretely as

$$y_i = \sum_j b_{i-j} x_j.$$

The discrete version can be formulated as a matrix multiplication ${\bf y} = {\bf B}~{\bf x}$:

$$\left[ \begin{array} {c} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ ... ...n{array} {c} x_1 \\ x_2 \\ x_3 \\ x_4 \\ x_5 \end{array}\right]$$

which in this case implements transient boundary conditions.

In our example, we deconvolve the time series ${\bf y}$ recorded by a seismograph with known blur filter ${\bf b}$. To estimate the unblurred image ${\bf x}$, we seek the unconstrained minimum of

$$\min_{\bf x} \Vert{\bf y} - {\bf B}~{\bf {\hat x}}\Vert$$

which yields the least square solution

$${\bf {\hat x}} = ({\bf B}^T {\bf B})^{-1} {\bf B}^T {\bf y}.$$

The routine that deconvolves the time series combines Jest vectors that represent ${\bf y}$, ${\bf b}$ and ${\bf {\hat x}}$, a Jest linear-operator-with-adjoint to implement the convolution ${\bf B}$, and finally Jest's conjugate-gradient solver to estimate the solution vector ${\bf {\hat x}}$. Figure 1 shows the known filter, the blurred input time series, and the deconvolved result.

 

tcai
Figure 1 Seismogram deconvolution. The top signal is the known filter. The second signal is the received, blurred signal. The third signal is the deblurred signal estimated from the signals above. The bottom signal is the correct answer.