An iterative method for restoring noisy blurred images
This paper introduces a new iterative image restoration method which is capable of restoring noisy blurred images by incorporating a priori knowledge about the image and noise statistics into the iterative procedure. The iteration equation consists of a prediction part which is based on a noncausal image model description and an innovation part which is weighted by a gain factor. The gain is computed using a linear MSE optimization procedure and is updated at each step of the iteration. The convergence of the algorithm, the resolution of some convergence difficulties by using "reblurring," and methods for the introduction of physical constraints will be discussed. This image restoration scheme can be interpreted as an iterative procedure with a statistical constraint on the image data. Results of several experiments with noisy blurred data are presented to demonstrate the feasibility of this approach. 1. Introduct ion Images are produced to provide useful informat ion about some phenomenon o f interest. Unfor tunate ly , since physical imaging systems are imperfect , a recorded image invariably represents a degraded version o f an original image or scene. For example, in aerial pho tography and remote sensing, the images can be degraded by atmospheric turbulence, optical aberrations, and relative mot ion between the camera and the scene. Medical images are typically o f low resolution and low contrast and electron micrographs are degraded by the spherical aberration o f the electron lens. The c o m m o n problem confronting researchers in these fields is that o f restoring the image data to improve quality. The problem is complicated by the r andom noise that is inevitably * Received September 21, 1983; revised December 15, 1983. This work was supported in part by the Joint Services Electronics Program under contract DAAG29-81-K-0024 and in part by the Netherlands Organization for the Advancement of Pure Research (ZWO). School of Electrical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, USA. 2 Information Theory Group, Department of Electrical Engineering, Delft University of Technology, Delft, The Netherlands. 140 KATSAGGELOS, BIEMOND, MERSEREAU, AND SCHAFER mixed with the data. This may originate with the image formation process, the transmission medium, the recording process, or any combination of these. Such noise represents a significant limitation to the perfect restoration of images. In many practical situations, the image degradation can be adequately modelled by a linear blur (motion, defocussing, atmospheric turbulence) and an additive white Gaussian noise process [1]. If the observed image is represented by an M x N array of real picture elements t y ( i , j ) : 1 <_ i <_ M, 1 _< j < N], then in the spatially invariant case, it can be described by the following two-dimensional convolution summation: y( i , j ) = ~ ~ b ( m , n ) x ( i m , j n ) + w ( i , j ) , (I) (m,n) E Wb where y ( i, j ) is the degraded image, x ( i, j ) is the original image, w ( i, j ) is the additive observation noise, which is assumed to be uncorrelated with the image data, and b ( m, n ) is the impulse response or point-spread function ~SF) of the imaging system that is introducing the blur. It is assumed that the support of the PSF, Wb, is much smaller than the size of the image. With this model, the problem of image restoration is represented as the problem of operating on the degraded image y (i, j ) in order to get an improved image 2 ( i, j ) which is as close to the original image x ( i, j ) as possible, subject to a suitable optimality criterion and given some prior knowledge about the PSF of the blur, the image, and the noise statistics. Inverse filtering techniques, which aim at perfect restoration of the image by using the convolutional inverse of the blur, become poor restoration techniques when noise is present [2]. This can also be the case with the class of iterative restoration algorithms which were introduced in [3], since these are essentially iterative implementations of the inverse or pseudoinverse filter. The iterative techniques, however, do have advantages when compared with the inverse filter, such as the possibility of incorporating physical constraints on the data [3], man-machine interaction [4], and the ability to deal with nonlinear or shift-varying blurs [3]. Therefore, considerable effort has been expended in trying to diminish the high noise sensitivity of the iterative procedures, while still producing reasonably sharp images. One such procedure called for lowpass filtering the observed data prior to applying a constrained iterative procedure [51. With another approach, which we will call "reblurring," the observed image is purposely degraded prior to the restoration. The purpose of reblurring is to ameliorate some convergence difficulties, but this procedure also prevents extreme noise amplifications [3,6]. One can also reduce the effects of noise by limiting the number of iterations or by imposing some type of stopping rule based on the error residual and the variance of the observation noise [7,8]. In all of these approaches, however, no attempt has been made to incorporate statistical knowledge of the image and the noise directly into the restoration scheme, as is common practice with Wiener and Kalman filters [9]. AN ITERATIVE METHOD FOR RESTORING NoisY BLURRED IMAGES 141 The purpose of this paper is to show that by incorporating a priori knowledge about the image and noise statistics directly into an iterative restoration procedure and by using an optimization criterion, we are able to obtain restored images with good resolution and without noise amplification, while still retaining the attractive features of iterative algorithms. This new restoration scheme can be interpreted as an iterative procedure with a statistical constraint on the image data. In Section 2, we discuss a noncausal autoregressive model for the undistorted image and procedures for identifying the model parameters. In Section 3, we introduce our new iterative algorithm, which consbts of a predictive part based on the noncausal image model and an innovation part weighted by a fiker-galn factor. This gain is computed using a linear MSE optimization and is updated at every step in the iteration. The convergence of the algorithm, the solution of some convergence problems by using reblurring, and the introduction of physical constraiflts into the iteration are also discussed in this section. Several experiments on noisy blurred data will be presented in Section 4. 2. Image modelling We will begin by assuming that the original undistorted image can be interpreted as a sample from a discrete homogeneous random field { x ( i , j ) : 1 _< i _< M, 1 _< j _< iV}, which in turn has a spatially invariant mean, #, and a translation invariant autocovariance function r (k , / ) . For convenience we will assume that the mean of the image has been estimated and subtracted. For a discussion of autoregressive models for images with nonzero means, see [lO]. Under these assumptions the original image x ( i , j ) can be modeled by the twodimensional autoregressive equation x(i, j) = E E a ( p ' q ) x ( i p ' J q ) +u( i , j ) , (2) (p,q) 6S where x ( i, j ) represents the image-intensity value at spatial coordinates ( i, j ) and where u (i, j ) can be viewed as either the input process or as the error in generating x (i, j ) . The support consists of a set of index pairs (p , q) , which are independent of ( i , j ) and which do not include the origin. The shape of S defines the type of model which is being used. Let us investigate the choice of the set S further. In general, the shape of S specifies that some samples of x (i, j ) must be computed before others. Causality, however, is a time-domain concept which has no meaning when talking about the spatial variables in an image. Imposing an ordering relation on the samples is unnatural and should be avoided if possible. Therefore, we restrict ourselves to the class of noncausal image models for which S --{ ( p , q ) : p2 + q~ <_ p 2 ' ( p , q ) ~ (0,0)} . (3) 142 KATSAGGELOS, BIEMOND~ MERSEREAU, AND ~CHAFER In this case, the image can be seen to be, in general, a sample from an anisotropic homogeneous random field that satisfies a noncausal stochastic difference equation [11,12,13]. Note that with this description we make no assumptions about the separability or the exponential form of the autocovariance function. There are different approaches that can be followed to estimate the parameters a ( p , q ) of Equation (2). If we have a noise-free unblurred prototype image, we could measure its autocovariance function under the assumption that the image data is homogeneous. Then the noncausal model could be fitted to the measured autocovariance function for different values of o by using a linear MSE fitting procedure. Which model ultimately serves to describe the image data would depend upon some model quality criterion [14]. It should be noted that unlike causal minimum variance models, noncausal minimum variance models are not driven by white noise [15,161. In the case where we have noisy image data, we could make use of the model parameter identification procedure described by Kaufman et aI. [17] and Hoogterp and Loh [18]. 3. A new iterative restoration algorithm In this section we introduce a new iterative algorithm that is capable of restoring noisy blurred images by incorporating a priori knowledge about image and noise statistics into the iterative restoration procedure. Furthermore, we derive a formula for the iterative filter-gain update procedure and analyze the convergence of the resulting algorithm. Finally, the possibility of incorporating physical constraints will be discussed. FORMULATION. As we have already discussed, exact image restoration is not possible and iterative filtering schemes which are based on the convolutional inverse of the blur become po