by Yirong Shen and Rui Zhang

Image restoration is a vital part of many image processing applications. The purpose of image restoration is to recovery the original image from a degraded observation. In many instances, the degraded observation g(x,y) can be modeled as the two-dimensional convolution of the true image f(x,y) and the point-spread function (also called the blurring function) h(x,y) of a linear shift-invariant system plus some additive noise n(x,y). That is, g(x,y) = f(x,y)*h(x,y) + n(x,y)

In many situations, the point-spread function h(x,y) is known explicitly prior to the image restoration process. In these cases, the recovery of f(x,y) is known as the classical linear image restoration problem. This problem has been thoroughly studied and a long list of restoration methods for this situation includes numerous well-known techniques, such as inverse filtering, Wiener filtering, least-squares filtering, etc.

However, there are numerous situations in which the point-spread function is not explicitly known, and the true image f(x,y) must be identified directly from the observed image g(x,y) by using partial or no information about the true image and the point-spread function. In these cases, we have the more difficult problem of blind deconvolution. Several methods have been proposed for blind deconvolution and active research continues in this area. In our project, we seek to understand and implement two of the proposed algorithms for blind deconvolution: the iterative blind deconvolution algorithm and the simulated annealing algorithm.

The goal of the general blind deconvolution problem is to recover two convolved signals, f and h, when only a noisy version of their convolution is available along with some or no partial information about either signal. Some important characteristics of the problem of blind deconvolution is listed below:

1. The true image and the point-spread function must be irreducible for an unique solution to the problem. A signal is reducible if it can be expressed as the convolution of two other signals (neither of which is the two-dimensional delta function). If either the true image or the point-spread function is reducible, then there exists more than one way to perform the deconvolution.
2. Without the use of appropriate a priori information, only a scaled, shifted version of the original image can be obtained through blind deconvolution.
3. In the presence of additive noise, the exact blind deconvolution of the observed image g(x,y) is impossible and only an approximate solution can be found.

In practice, all blind deconvolution algorithms require some partial information to be known and some conditions to be satisfied. The iterative blind deconvolution algorithm and the simulated annealing algorithm are two of the most general algorithms with the least amount of required a priori information. Both methods require the true image and the point-spread function to be nonnegative with know finite support. These requirements are not overly constricting since image intensity is always nonnegative and support information can often be inferred in many image processing situations.

The iterative blind deconvolution algorithm was first proposed by Ayers and Dainty in 1988. It is one of the most popular and most famous algorithms in blind deconvolution. The basic structure of the algorithm is a simple iterative one and is illustrated in the figure below.

Starting with a random initial guess for the true image, the algorithm alternates between the image and Fourier domains, using known constraints in each domain. The image domain constraints are simply those imposed by the nonnegativity requirement and the known finite support information. Negative-valued pixels and nonzero pixels outside the region of support are set to zero. The Fourier domain constraint is a consequence of the fact that if the noise level is low, then the product of the Fourier transforms of the true image and the point-spread function should be approximately equal to the Fourier transform of the observed image. Using this fact, we can use inverse filtering to estimate the point-spread function (image) using the Fourier transforms of the degraded observation and the estimate of the true image (point-spread function). However, inverse filtering can result in the magnification of noise in regions where the function being inverted has low values. There are several different ways to deal with this problem. In our experiments, we found that the best performance is achieved when the following formula is used to generate a new estimate of H at the k-th iteration:

The formula for generating a new estimate of F is similar. The above formula resembles the Wiener filtering formula and it has the similar effect of providing robustness against noise. The constant beta represents the energy in the additive noise and must be chosen carefully for the algorithm to perform well.

The figures below show an example of the results that were obtained using the iterative blind deconvolution algorithm. Since scale factors cannot be recovered, all images are scaled so that the maximum pixel value is 1. The images are displayed using a gray scale in which 1 represents black and 0 represents white. The additive noise used in this case was white Gaussian noise with zero mean. The signal-to-noise ratio (SNR) was calculated to be 45dB. The results were obtained after 1000 iterations. The deconvolved images are clearly recognizable and appears to be slightly noisy versions of the originals.

The following figures contain the results of a different test run on the same observed image but using a different initial random estimate of the true image f(x,y). Although some features of the originals can still be recognized, the deconvolved images in this case are of much lower quality than the ones above.

This example illustrates some of the problems plaguing the iterative blind deconvolution algorithm. Although the IBD algorithm can be computed easily and quickly, it lacks reliability in its performance. The deconvolution results are sensitive to the initial image estimate, and the algorithm is sometimes unstable. In addition, the uniqueness and convergence properties of the IBD method are still unknown. The lack of theoretical understanding of the IBD algorithm has hindered its further maturation into a widely applicable algorithm in practice.

Simulated annealing is a well-known Monte-Carlo global minimization technique. Simulated annealing minimization is analogous to the annealing process of metals. When liquid metal is cooled, the random motions of the metal atoms are reduced. If the cooling is done sufficiently slowly, then the metal will reach the absolute minimum energy state related to the complete atomic ordering of the metal. If the cooling is done too quickly, then a suboptimal energy state may result.

In simulated annealing minimization, a cost function Q (analogous to the energy of the metal) is minimized through the random perturbation of the minimization parameters. As the iterations progress, the temperature parameter is gradually reduced and the amount of random changes is also decreased. Just as in the case of the annealing of metals, the temperature must be reduced sufficiently slowly every iteration in order for the method to reach the global minimum of the cost function. In the limit of infinite precision and infinitely many iterations, simulated annealing is guaranteed to reach the global minimum of the cost function.

The algorithm starts with random estimates, f₀ and h₀, of the true image and point-spread function respectively, and proceeds as follows:

1. Set cycle counter n_c = 1, set scan counter n_s = 1.
2. Calculate values for T (temperature) and a (scale of perturbations).
3. Scale the estimates of the true image and the point-spread function by b and 1/b respectively, where b is chosen such that the scaled functions have the same root mean square value. (Note: the convolution of the two estimates is not changed after scaling.)
4. The following steps are performed for each pixel in the support of the estimates of the true image.
1. Perturb the current pixel with a random value uniformly distributed in the range [-a , a ].
2. If new pixel value is negative, then set it to zero. (Nonnegativity constraint)
3. Calculate dQ the change in the cost function.
4. If dQ <= 0, then the perturbation is accepted. If dQ > 0, then the perturbation is accepted with probability exp(-dQ/T).
5. If n_s < N_s, then increment n_s and return to step 4.
6. If n_c < N_c, then increment n_c and return to step 2, otherwise stop.

The actual details of step 2 is left ambiguous because it can be implemented in several ways. Special care is needed here because incorrect choices can lead to the algorithm being trapped in a local minimum or a excessively slow convergence. In step 4, some perturbations that increase the value of the cost function are accepted to prevent the algorithm from being trapped in a local minimum (This only works when temperature is reduced slowly)

Results from two examples are shown below. For both examples, the additive noise was uniformly distributed with zero mean and the SNR was 20dB. The contamination level level c is defined as c=En{n}/En{f*h}=1/SNR. In our examples, c=0.01. In general, the value of the cost function Q cannot be reduced to 0 when the convolution is contaminated with additive noise. The best that can be done is to find solutions with Q approximately equal to (on the same order of magnitude as) c. For both examples, the temperature at the k-th cycle, T_k, is varied as T_k = 0.8T_k-1.

In the first example, the algorithm is run with N_s = 50, N_c = 40, and T₁ = Q(f_0, h₀, g)/10. The value of a at the k-th iteration, denoted by a _k is varied as a _k = 50(T_k^1/2).

In the seond example, the parameters N_s = 25, N_c = 30, and T₁ = Q(f_0, h₀, g)/100 were used. The value of a was varied according to the ratio of the total number of accepted perturbations to the total number of rejected perturbations over the k-th cycle. If this ratio is less than 0.8, then a is scaled by a factor of 0.8; if the ratio is greater than 1.2, then a is scaled by a factor of 1/0.8; if the ratio is between 0.8 and 1.2, then a is not changed.

The results of the deconvolution in both examples are of high subjective quality and can be easily recognized in comparison to the originals. The final values of the cost function Q when the algorithm stopped were 0.006 and 0.005 respectively, which are both slightly smaller than the contamination level. Considering the small sizes of the images used, the examples took a relatively long period of time to run. The second example required approximately 30 minutes on a 500 MHz DEC Alpha workstation.

In this project, we have studied two proposed algorithms for the blind deconvolution of images, iterative blind deconvolution and simulated annealing. We have demonstrated that both algorithms can be applied successfully to the blind deconvolution problem and produce high quality results using only minimal partial information. However, problems exist with both algorithms that prevent them from being widely applied in practice. The iterative blind deconvolution algorithm is sufficiently fast, but it lacks reliability and has unknown convergence properties. The simulated annealing algorithm has guaranteed reliability and convergence in theory, but it is too computationally intensive and would require an excess amount of time to run on realistically sized images. The two algorithms exhibit a poor tradeoff between computational complexity and convergence properties. This poor tradeoff leaves ample room for possible future research in this area.

1. G. R. Ayers and J. C. Dainty, "Iterative blind deconvolution method and its applications," Optics Letters, vol. 13(7), pp. 547-549, July 1988.
2. B. C. McCallum, "Blind deconvolution by simulated annealing," Optics Communication, vol. 75(2), pp. 101-105, Feb. 1990.
3. D. Kundur and D. Hatzinakos, "Blind image deconvolution," IEEE Signal Processing Magazine, vol. 13, pp. 43-64, May 1996.
4. R. H. T. Bates and H. Jiang, "Blind deconvolution --recovering the seemingly irrecoverable!," in International Trends in Optics (J. W. Goodman ed.), pp. 423-437, 1991.