CC BY
51Electronic Journal «Technical acoustics» http://webcenter.ru/~eeaa/ejta/
2003, 2
M. Djebbouri, D. Djebouri and R. Naoum
Telecommunications and Signal Processing Laboratory
Electronic Department, Djillali Liabes University, BP 89, Sidi Bel Abbes, Algeria 22000 e-mail: mdjebbouri@hotmail.com
Fourier And Wavelets For Blind Image Restoration
Received 28.10.2002, published 21.01.2003
This paper describes a technique for the blind deconvolution based on the wavelet domain deconvolution that comprises Fourier-domain followed by wavelet-domain noise suppression, in order to benefit from the advantages of each of them. The algorithm employs regularized Wiener filter, which allows it to operate even when the system is non-invertible. In fact, we model such image to be the result of a convolution of the original image with a point spread function (PSF). This PSF depends mainly on the image formation system. Unfortunately, it is often very difficult to model this PSF from the physical data, for this reason we consider the problem as a blind deconvolution. First, the identification of the blur is based on maximum likelihood and the solution is obtained iteratively by successive estimations of the PSF from the noisy blurred image. We propose a blind restoration by estimating the noise variance, the point spread function (PSF) and the original image from a blurred and noisy observation. Our method is based on regularized Wiener filter and RDWT (redundant discrete wavelet transform). We illustrate the results with simulations on some examples.
1. INTRODUCTION
The desire to remove distortion introduced by a spatially invariant point spread function (PSF) is a challenge faced in many imaging applications. The field of image restoration is generally concerned with the estimation of uncorrupted images from the noisy and blurred images acquired by imaging Systems. The image y, of an object x is given by twodimensional convolution.
y = h ® x + v , (1)
where x is the original undistorted image, y is the distorted noisy image, h is the PSF of the system, ® is the convolution operator, and v is the corrupting noise.
The diffraction of the optics of a camera produces a blur that is equivalent to a low pass filter that attenuates the highest image frequencies. Moreover, the electronics of the photoreceptors add a noise typically a white noise. The problem of estimating x given knowledge of y and h is referred as the inverse problem. The inverse problem is difficult to solve due to its ill-posed nature. Any noise in the measured imagery complicates the problem further. A solution to the inverse problem is usually defined as an element of set of feasible solutions consistent with the linear relationship (1). For simplicity, we assume periodic boundary conditions, which means that the convolution is circular and that the noise v is circular
stationary. One then seeks to incorporate prior knowledge in order to narrow the choices and produce better estimates. In most practical situations, the first step in restoring an image is to identify the kind of degradation the image has suffered. For example, practical satellite images are often blurred due to limitations such aperture effects of the camera, camera motion, or atmospheric turbulence.
There exist many classical and modern methods for solving linear inverse problem. The classical methods usually are deterministic in nature. The singular value decomposition (SVD) paradigm has been used, it has been shown to be mini-max best over certain classes of homogeneous functions [1]. However, as pointed out in [2], within such paradigm the bases are completely determined by the convolution operator and the underlying signal structure is totally ignored, which is very undesirable. One extreme is to perform no Fourier-domain regularization, this is equivalent to WVD (wavelet-vaguelette deconvolution) approach of Donoho [3] and the mirror-wavelet basis approach of Kalifa et al. [4]. The mirror-wavelet basis approach described in [4] (wavelet packet) adapts to the frequency response of the convolution operator H. Though the adapted basis improves upon the WVD. While multiscale/wavelet-based methods can avoid some of the drawbacks of traditional restoration methods, the existing method do have some limitations. First of all, some of the multiscale methods are analogous to classical Wiener filtering, reinterpreted in the wavelet-domain. These methods may have computational advantages, but do not overcome some of the pitfalls associated with linear filtering methods. However, Banham and Kastaggelos [5] apply multiscale Kalman filter to the deconvolution problem that circumvents this problem by switching the linear filter spatially based on an edge detection test. Their approach employs under-regularized constrained-least-squares prefilter to reduce the support of the state vectors in wavelet domain and a complicated prediction on edge and non edge quad-trees overcomplete wavelet basis. Recently a new method was developed by Neelamani and al [6] that combines classical Fourier domain regularization methods (Wiener filter) with wavelet domain thresholding/shrinkage. The advantage of such a formulation is the suppression of noise in singular or near singular operator situations. Most deconvolution algorithm are carried out using observed PSF, because theoretically computed PSF are usually not a very good match to the observations. The noise in observed PSF presents a problem, however when using a noisy PSF, the restoration algorithm ought to account for the fact that both the image and the PSF are noisy. The blind deconvolution approach allows one to construct an algorithm that takes as input a noisy image and a noisy PSF observation (this PSF is estimated and not known) and to construct as output estimates of both the PSF and the unblurred image. Blind deconvolution [8, 9, 12, 13] algorithms are predominately iterative and attempt to recover the image and PSF from a blurred observation using a variety of constraints such as non-negativity and object support.
This paper is structured as follows: in the next section, we review the methods related to our work. We present an algorithm BIWD (blind iterative Wiener deconvolution) scheme and an improved one, BWARD (blind wavelet-based regularized deconvolution) algorithm in section 3. Illustrative examples lie in section 4.
2. RELATED METHODS
We start with the imaging equation (1):
y(U J) = X h(P q)x(i - P’ J - q) + v(i,J) ,
(p,q ")£Sd
(2)
where we assume that the original image, denoted by x(i, J), i = 0,... N -1, J = 0,... N -1, and SD is the finite support region of the distortion filter.
Two degradations corrupt our observation y of the desired data x: convolution with a linear invariant system, having impulse response h and additive noise v of variance a2 (see Fig. 1).
observation
y
Fig. 1. Convolution model set-up
Our approach is based on two sections: blur operator identification and image estimation. Here we investigate its performances and theoretical properties in greater detail. Our main contribution in this paper is a combination of two complementary approaches: (i) blur estimation of Katsaggelos and Lay [8], (ii) wavelet based deconvolution of Neelamani et al. [6]. In [8] a blind deconvolution method without using the advantages offer by the wavelets is proposed, and in [6] a non blind deconvolution method is described. Basically we have combined the two methods by adding some modifications in that way the combined method (Katsaggelos and Neelamani) gives better performances.
Blur estimation
In this section, we give a brief derivation of the method of Katsaggelos and Lay. For more details, the reader is referred to [8]. The iterative algorithm developed by Katsaggelos and Lay arrived to the following equations:
Mïiy (m n)
H
( p)*
(m, n)S(Xp )(m, n)
H(p)(m, n) SXp)(m, n) + a.
2( p ) V
-Y (m, n),
(3)
SXP )(m, n)o2 p )
H(p)(m, n) SXP)(m, n) + a,
( p),
■ 2( p ) V
(4)
SX+1)(m, n) = S(xpl (m, n) + — M« y (m, n)
N2
1 X / y'
(5)
2
H(p+1)(m, n)
Y (m, n)Mxi*y (m, n)
N 2
(m,n) +
N 2
M(p)
1V1 X / j
(m, n)
(6)
1 N-1 N-1 2
P *1) = 77T ZLlH("1)(m-”)
N m=0 n=0 L
sX’K (m.«) - -fir
MX,%
(m, n)
+ -Nf (y(m, n)2 - 2Re[Y* (m, n)H(p+1) (m, n)MX^ (m, n)] .
(7)
In (7) Y(m, n) and M'xP/)y(m,n) are respectively the 2D DFT's (Discrete Fourier
Transform) of the observed image y(i, J) and restored image. Re[n], n and \tf\ denote Real part, complex conjugation and the magnitude of the complex numbern respectively. SX (p+1)(m, n) and H(p+1)(m, n) represent respectively the power spectral density of the input signal and the Fourier Transform of the PSF at the (p +1) iteration. Equation (3) shows that the restored image is the output of a Wiener filter, based on available estimate of the blur operator and the noise variance, with the observed image as input. The iterative algorithm developed in [8] (equations 3, 4, 5, 6 and 7) compute the restored image and the PSF. The Blackman-Tukey algorithm was used to compute an estimate of the power spectral density of y, which in turn was used as SX(0) , the initial estimate of the power spectral density of x. The 2D impulse ( h( x, y) = 1 for x = y = 0, and h( x, y) = 0 elsewhere) was used as h(0), the initial estimate of PSF. The method of Katsaggelos and Lay is computationally very intensive, as it consists of two nested iterative loops and, the algorithm takes a very large number of iterations to converge especially when the variance of the additive noise is important.
Image estimation
The method proposed in [6] consists of the following steps (see Fig. 2):
1. Pure inversion: similar to the method used in [3].
2. Fourier-domain noise attenuation: employ a small amount of Fourier domain shrinkage using weights to achieve partial noise attenuation in the estimate ~a.
3. Wavelet-domain signal estimation: Since the estimate in step 2 contains some residual noise, we shrink the wavelet coefficients of at each scale according to the noise variance at that scale to obtain the final estimate.
1
2
2
Fig. 2. Direct deconvolution in three steps (inversion, Fourier and Wavelet shrinkage)
Wavelet-domain signal estimation remains effective, since the noise corrupting the wavelet coefficients is not excessive. In this method, we have assumed knowledge of the convolution operator. However, in many cases such as most practical imaging systems, the convolution operator is also unknown. In such blind deconvolution problems, the convolution operator must be estimated from the observation after that, this algorithm can be used to perform the image restoration. If H is non-invertible the deconvolution method of Fig. 2 fails, and since Wiener filter is know to be optimal among linear filters, we used instead a deconvolution in two steps see (Fig. 3).
y
observation
x
estimate
Fig. 3. Direct deconvolution in two steps (Regularized Wiener filter and Wavelet shrinkage)
Fourier non-iterative image restoration techniques have as their main advantage the short computation time they need to produce a solution, when compared with more elaborate iterative restoration methods. Also, they are strictly linear by design. The reason behind that computational efficiency lies in the fact that the computations are made in the Fourier transform domain, where the 2-D matrix operations needed for inverting the imaging equation are expressible as 1-D vector-only operations. Also, the existence of very efficient FFT algorithms adds to the overall efficiency of the technique.
3. PROPOSED METHOD
In this section, we present two algorithms, the blind iterative Wiener deconvolution (BIWD) used as a PSF estimator and the blind wavelet-regularized deconvolution (BWARD) algorithm. The noise power can be estimated using the median estimate of [2] performed on the finest scale wavelet coefficients (where the signal energy is expected to be negligible).
Algorithm BIWD
We introduced a small modification on equations (3, 4, 5, 6) above to make them applicable for the blind deconvolution that we consider in this paper. The modification is the introduction of a regularization parameter a that conducted to a simplification of the algorithm described in [8] (we don't have to estimate the noise variance after each iteration). The equations (3, 4, 5, 6) become:
M ax I y (m n) =
H( p)* (m, n)S(xp\m, n)
H( p)(m, n)2 S(xp )(m, n) + aa
-Y(m, n).
(8)
(m, n)
S(XP) (m, n)aa,
H(p)(m, n) S(xp)(m, n) + aa,
(9)
(10)
(11)
We can use equation (11) to identify the PSF and then, we restore the image while using the equation (8).
To estimate SX, the power spectral density of the input signal, we use the equations (8, 9, 10) above. Pictures shown in Fig. 7(c) and 8(c) have been gotten by this method, (see test 1 and test 2).
Algorithm BWARD
In the previous section, we focused on the estimation of PSF. After that, Ward [6] (Wavelet-based regularized deconvolution) can be used to perform the deconvolution in the wavelet domain. BIWD (blind Wiener iterative deconvolution algorithm) generally converge after 2 or 3 iterations starting from a 2D impulse as initial value of the PSF, it is reasonable to think, that if the initial guess for the PSF was sufficiently good, then the algorithm converge faster to a more accurate PSF, which in turn leads to better restoration of the unblurred image. The PSF is estimated with the same method described above (BIWD), only the initial guess was changed (the number of iterations is fixed to one), and then restore the image by using the modified wavelet domain deconvolution algorithm. The implementation of the regularized inverse filter involves the estimation of the power spectrum of the original image in the spatial domain. The final algorithm is shown below, where 2D uniform blur is defined as h(x,y) = 1/d2 for (x,y) e SD, the support size of the PSF is dxd and h(x,y) = 0 elsewhere. This of course results in:
Input initial estimates for image x (the observed imagey), and PSF h (uniform 2D blur).
1. Estimate the noise variance.
2. Estimate the power spectral density of the input signal.
3. Estimate the near-optimal regularization parameter a.
4. Apply one iteration of the algorithm BIWD to estimate PSF.
5. Apply modified Ward algorithm (Fig. 5) to restore the image.
X h( x, y) =1
(12)
( x, y )ßSD
Our BWARD algorithm is shown in table 1:
Table 1. BWARD algorithm
One used a part of the BIWD algorithm to estimate the PSF, just the initial conditions (uniform blur) changed. We use the most simple Daubechies wavelets (Haar wavelet) and we apply the redundant wavelet decomposition and reconstruction algorithms of complexity O(lo n) for n given data (n=NxN pixels).
For the restoration of the image, we have modified the initial algorithm Ward shown in Fig. 4 by changing the position of the Fourier-domain shrinkage (see Fig. 5). The regularization parameter is chosen to be level dependent (same at each level). We have used 4 decompositions level in our simulations. To compare the results of the blur identification, the following figure of merit was used:
h - h
error ■
(13)
where h and h denotes the original and the estimated PSF, and SD their respective supports.
The noise added to the blurred image was evaluated by measuring the blurred signal to noise ratio denoted by
BSNR = 10log
x *
10
(14)
And the performance of the restoration algorithm was evaluated by measuring the improvement in signal to noise ratio defined in decibel (dB) by
ISNR = 10log
10
x - x
(15)
where x, x and y are respectively the original, restored and degraded images, ~ is the shifted version of the observed image that minimizes the ISNR.
Fig. 4. Wiener and Wavelet shrinkage with one regularization parameter. Blocks H0 and H1 denote the low-pass and the high-pass filters for the filter bank
2
2
2
2
4. NUMERICAL RESULTS
All of our code was written in MATLAB to take advantage of its visualization capabilities. Experiments will be given which show the performances of BIWD and BWARD algorithms. Our current implementation is not restricted to small images but larger images than 512x512 pixels take more time to compute the restored image and needs more memory. Many choices of the thresholding value is possible, usually it is chosen based on some estimate of the noise level av (standard deviation of the noise) in the data. In our case the shrinkage parameter is
fixed to -y/2 log(N2)av for simplicity, where N2 is the number of samples and hard
thresholding is performed to remove the residual noise. This choice is suggested in [2] as a probabilistic upper bound on the noise level (Donoho, Johnstone and collaborators have proven a lot of optimality for this choice this is the so-called universal threshold). Other choices include the SURE-threshold [2] and Bayesian estimates [14, 15]. Clearly one threshold cannot remove noise decently the amount of noise av depends on the resolution
level. Scale-dependent thresholds are more adaptive to signal characteristics (Ward algorithm is based on Fig. 4 and the BWARD is based on Fig. 5).
Comparison with the blind Lucy-Richardson algorithm
The blind Lucy-Richardson algorithm method is an extension of the well-known Lucy-Richardson method [11, 12]. The most common and efficient implementation makes use of the FFT to compute convolution. This implicitly imposes periodic boundary conditions on the image. The blind version is similar to the original method; each iteration alternately uses several iterations of the non-blind algorithm to estimate a new PSF and then a new image. It is generally more effective for images with fewer sharp edges, since convolution tends to smooth edge boundaries. For comparison, we included this method, which is widely used in image processing for deconvolution, especially in the astronomical domain.
Fig. 5. Wiener and Wavelet shrinkage with one regularization parameter for each wavelet scale. Blocks H0 and H1 denote the low-pass and the high-pass filters for the filter bank
Test 1
Our first test consists of an image of cameraman of size (256x256). The original PSF is a Gaussian blur with variance 10 truncated to a support of size (7x7). The blurred image was generated by convoluting the original image and the PSF (no noise is added in this test). The original and blurred images are shown in Fig. 6(a) and Fig. 6(b) respectively. An image (cameraman) is restored (see Fig. 6(e)) using the BIWR algorithm, the restored image shows significantly more detail than the blurred one (ISNR~12dB). In contrast to the Lucy Richardson method (Fig. 6(f)), the smooth region and most edges are well preserved in the image estimate. Because the Fourier basis used by the algorithm have support over the entire spatial domain, ripples in the restored image result. The restored PSF is shown in Fig. 6(d) and presents an error of 0.0200 it looks like the original one. When we take as an initial guess for the PSF a uniform blur, the method converge in one iteration.
Test 2
In the previous test (in test 1 no noise was added to the blurred image), our second test consists of an image boats of size (256x256). The blurred image was obtained by convoluting the original image and the PSF (the same as in test 1). The original and blurred images are shown in Fig. 7(a) and Fig. 7(b) respectively. Random, zero mean noise is added to the blurred image, resulting in a degraded image (BSNR = 40 dB). The restored PSF is shown in Fig. 7(d) and it is similar to the original one. An image (boats) is restored (see Fig. 7(f)) using our algorithm BWARD. The restored image shows significantly more detail than the blurred image and the edges are well preserved, (ISNR = 9.0 dB). It is slightly better than the one obtained by the non blind algorithm Ward proposed by Neelamani et al. [6] (Fig. 7(e) ISNR = 8.7 dB), although this last method uses the original PSF (this unexpected result may be explained by the fact that, the PSF act as a filter for white noise, so we have a colored noise which is better eliminated by an algorithm that takes some in consideration this point, see Fig. 4 and Fig. 5).
CONCLUSIONS
We have presented a blind deconvolution algorithm to enhance the resolution of images obtained by a system who produces a degraded image. The edges are preserved and the computational cost is reasonable. The proposed method is an improved version of the Ward algorithm by adding the estimation of the PSF. Our method (BWARD) can be potentially employed in a wide variety of applications such as satellite imagery, and seismic deconvolution to obtain enhanced deconvolution estimates. Several drawbacks still to overcome: (i) it is not trivial to use non-Gaussian noise, (ii) the investigation on the choice of the optimal regularization parameter.
Fig. 6. Image Cameraman, from top to bottom clockwise:
a) Original image
c) Original PSF (7x7 pixels)
e) BIWD method
b) Degraded image
d) Recovered PSF (error=0.0200)
f) Lucy-Richardson method
Fig. 7. Image Boats, from top to bottom clockwise:
a) Original image (256x256) b) Blurred and noisy image (BSNR = 40 dB)
c) BIWD method (notice the ripples) d) Recovered PSF (error = 0.0200)
e) Ward method f) Our method (BWARD)
REFERENCES
1. P. C. Hansen. Numerical aspects of deconvolutions, lectures notes. Department of Mathematical Modelling, Technical University of Denmark, 2000.
2. D. L. Donoho and I. M. Johnstone. Adapting to unknown smoothness ideal via wavelet shrinkage. J. Amer. Stat. Assoc., vol. 90, pp. 1200-1224, Dec. 1995.
3. D. L. Donoho. Nonlinear wavelet methods for recovery of signals, densities, and spectra from indirect and noisy data, in Different Perspectives on wavelets. Vol. 47 of Proc. Symp. Appl. Math., pp. 173-205, American Mathematical Society, 1993.
4. J. Kalifa, S. G. Mallat, and B. Rougé. Minimax solution of inverse problems and deconvolution by mirror wavelet thresholding. SPIE conference on Wavelet applications in signal and image processing VII, vol. 3813, (Denver, CO), pp. 42-57, Jul. 1999.
5. M. R. Banham and A. K. Katsaggelos. Spatially adaptive wavelet-based multiscale image restoration. IEEE Trans. Image Processing, vol. 5, pp. 619-634, April 1996.
6. R. Neelamani, H. Choi, and R. G. Baraniuk. Wavelet-based deconvolution for ill-conditioned systems. Dept of Rice University, Submitted to IEEE Trans on Image Processing February 2000.
7. S. Mallat. A Wavelet tour of signal processing. San Diego: Academic Press, 1998.
8. A. K. Katsaggelos and K. T. Lay. Maximum likelihood blur identification and image restoration using the EM algorithm. IEEE Transactions on Signal Processing, vol. 39, No. 3, March 1991.
9. A. K. Katsaggelos. Digital Image Restoration. New-York: Springer-Verlag, 1991.
10. M. Lang, H. Guo, J. E. Odegard, and C. S. Burrus. Noise reduction using an undecimated discrete wavelet transform. IEEE Sig. Proc. Letters, vol. 3, pp. 10-12, Jan. 1996.
11. D. C. Biggs and M. Andews. Acceleration of Iterative Restoration Algorithms. Applied Optics, vol. 36, No. 8, pp. 1766-1775, March 1997.
12. M. Djebbouri and D. Djrbbouri. A blind deconvolution for image restoration. Semémaire sur le traitement du signal et de l’image, Djelfa, Algérie, Novembre 1999.
13. A. Pruessber and D. P. O'Leary. Blind deconvolution using a regularized structured total least norm algorithm. Department of computer science and institute of advanced computer studies. University of Maryland, college Park, MD 20742, Sept. 2001.
14. M. Janseen and A. Bultheel. Multiple wavelet threshold estimation by generalized cross validation for images with correlated noise. IEEE transactions on image processing, 8(7), pp. 947-953, July 1999.
15. M. Janseen and A. Bultheel. Geometrical priors in a Bayesian approach to improve wavelet threshold procedures. Wavelet applications in signal and image processing VII, vol. 3813 of Spie proceedings, pp. 580-590, July 1999.
