Iterative Image Restoration

For the algorithms discussed here, the restoration approaches have been formulated to handle the problem of an image degraded by a spatially invariant or spatially variant blur operator with additive Gaussian noise. This description represents the degradation present in the HST WFPC images. The image degradation process can be modeled according to

where a lexicographical ordering of the observed digital image, , the original image, , and the additive noise, , is used. represents the degradation operator of the imaging system. The image restoration problem calls for obtaining an estimate of given , , and some characterization of the noise process, .

Iterative approaches may be used to solve this inverse problem very effectively. The primary advantages of iterative techniques are (Schafer et al. 1981, Katsaggelos 1989): (i) there is no need to explicitly implement the inverse of an operator; (ii) knowledge about the solution may be directly incorporated into the restoration process; (iii) the process may be monitored as it progresses; (iv) the effect of noise may be controlled with certain constraints; and (v) parameters determining the solution can be updated as the iteration progresses.

The iterative techniques applied here to HST data are developed through a set theoretic approach. In this approach, prior constraints on the solution are imposed by (Katsaggelos et al. 1985, Katsaggelos 1989, Katsaggelos et al. 1991)

and

where the solution belongs to both ellipsoids described by Eqs. (2) and (3). In Eq. (2), the operator represents a highpass filter which bounds the high frequency energy of the restored image. If the bounds and are known, and the intersection of and is not empty, the solution may be found by solving

where , the regularization parameter, is equal to , which controls the trade-off between fidelity to the data, and smoothness of the solution.

Generalized Iteration Adaptive Algorithm

In order to solve for an optimal restored image and regularization parameter at each iteration, we use

as the functional to be minimized. The weighting matrices and are included to make the restoration algorithm spatially adaptive. The necessary condition for a minimum is that the gradient of with respect to be equal to zero. This results in

When the noise and the high pass filtered image are stationary, the weighting matrices and become symmetric, and thus the equation can be rewritten as

and it becomes

since with a proper choice of . When the noise and the high pass filtered image are uncorrelated with themselves, even though they are nonstationary, and become diagonal, and therefore we obtain the solution for Eq. (7). When the noise and the high pass filtered image are white, then the weighting matrices become constant identity matrices multiplied by the constant variances. The highly nonlinear term can be removed by the proper choice of the regularization functional, based on the global convexity of the smoothing functional. Since Eq. (8) is nonlinear, we can not solve for in a direct way, but we can use an iterative technique.

The regularization parameter is defined here as a function of the original image (but in practice becomes a function of an estimate of the original image). The form of the smoothing functional to be minimized is of great importance since it preserves convexity and exhibits only a global minimizer. After investigating the desirable properties for the regularization functional to satisfy, the following two forms have been shown to provide optimal solutions (Kang and Katsaggelos 1993)

where , and

where controls convergence and convexity.

Given an optimal choice for , we can solve Eq. (8) by the method of successive approximations with (Schafer et al. 1981, Katsaggelos et al. 1985, Katsaggelos 1989, Katsaggelos et al. 1991):

Using either choice of regularization functional (Eq. (9) or Eq. (10)), the iterative algorithm does not depend on the initial condition despite the nonlinearity of the iteration. This is due to the convexity of the functional and the convergence criteria satisfied by the globally optimal iteration. The algorithm is applicable to any type of degradation and stabilizing matrix (both of which may be spatially varying). So, there is no requirement that and be block circulant matrices. It is also important to note that no knowledge of the noise variance, or of the bound which determines the ellipsoid that expresses the smoothness of the image, is assumed.

Frequency Domain Iterative Adaptive Algorithm

The second iterative algorithm examined here is a nonlinear frequency domain algorithm in which the regularization parameter is frequency dependent, and updated at each iteration step. In this case, the algorithm considers that is a block circulant matrix representing a spatially invariant blur. Also, the set-theoretic formulation is constructed in a weighted space, such that

and

where and are both block circulant weighting matrices. These matrices are chosen to maximize the speed of convergence at every frequency component, and to compensate for the near-singular frequency components of the iteration. The solution which belongs to the intersection of the ellipsoids given by Eqs. (12) and (13) is given by

We define , and . The block circulant matrix is the spatial domain representation of the relaxation parameter which will be called in the frequency domain, and is the block circulant spatial domain representation of the regularization parameter which is shown next in the frequency domain as . Again, successive approximation may be used to solve Eq. (14). Since all of the matrices in this equation are block-circulant, the iteration may be written in the discrete frequency domain as

where represents a single 2-D frequency component. In this case,

Here, the term we use is defined by

where , and

and

The optimized parameter is given at each frequency component by (Strand 1974)

According to this iteration, since and are frequency dependent the convergence of the iteration can be accelerated, making this an attractive algorithm where speed is a concern.

Optimization for HST Data

In the course of testing our algorithms, we have made appropriate consideration of the optimization of these algorithms for HST data. In particular, we have included a positivity constraint at each step of the iteration. This imposes the condition that there should be no negative flux in our image source. The algorithms were developed with additive Gaussian noise as the observation noise in the model. So, these algorithms are well equipped to deal with the read-out noise problem of the WFPC.