1. Introduction

The mathematical model for image formation is given by the linear operator equation

where d is the blurred, noisy image, f is the unknown true image, η is additive noise, and S is the blurring operator. In the case of spatially invariant blurs, Sf can be written as a convolution of the associated point spread function (PSF) s and the object f; that is,

However, if the blur is spatially variant, Sf cannot be written as a simple convolution operation. Space variant blur can occur from distortions due to telescope optics, such as was the case for the original Hubble Space Telescope Wide-Field/Planetary Camera, which had a large amount of spatial variation because of errors in the shaping mirrors [1]. Other important examples include objects moving at different velocities [2], and turbulence outside the telescope pupil [3]. Accurately modeling the spatially variant PSF in these applications can lead to substantially better reconstructions of the image. On the other hand, allowing for a fully spatially-varying PSF results in a computationally intractable problem [4]. Thus one typically makes the assumption that s is spatially-invariant on subregions of the image. That is, that the image d can be broken up into regions {Ω _n ${\}}_{n=1}^p$ such that the blur in each region is spatially-invariant. Then if we let c_n be the indicator function on Ω _n , i.e. c_n (u,v) = 1 for (u,v) ∈ Ω _n and 0 otherwise, and we define s_n to be the (approximately) spatially-invariant PSF on Ω _n , the PSF s will be given by

Substituting (1.2) into (1.1), we obtain

where f_n (ξ, η) := c_n (ξ, η)f(ξ, η).

Computational methods for the restoration of images arising from (1.3) have been sought by several researchers, though only in the case that the s_n ’s are known. Faisal et al. [5] apply the Richardson-Lucy algorithm, and discuss a parallel implementation. Boden et al. [6] also describe a parallel implementation of the Richardson-Lucy algorithm, and consider a model that allows for smooth transitions in the PSF between regions. Nagy and O’Leary [7] use a conjugate gradient (CG) algorithm with piecewise constant and linear interpolation, and also suggest a preconditioning scheme (for both interpolation methods) that can substantially improve rates of convergence.

In this paper, we present an algorithm for restoring images arising from model (1.3) when the s_n ’s are unknown. In our approach, the PSF in each region is estimated using the technique of phase diversity [8]. Once the PSF estimates have been obtained, we use a global restoration scheme of Nagy and O’Leary mentioned in the previous paragraph to obtain the restored image.

The paper is organized as follows. We begin in Section 2 with a description of the technique of phase diversity. Our computational method is then presented in Section 3. Results on some tests using simulated star field images taken through atmospheric turbulence are provided in Section 4, while conclusions and comments on the spatially-varying blur removal problem are given in Section 5.

2. PSF estimation and phase diversity

As discussed in, e.g. [9], with atmospheric turbulence the phase varies both with time and position in space. Adaptive optics systems use deformable mirrors to correct for these phase variations. Errors in this correction process arise from a variety of sources, e.g., errors in the measurement of phase, inability of the mirror to conform exactly to the phase shape, and lag time between phase measurement and mirror deformation. Thus, a spatially variant model can potentially produce better restorations.

In this section we describe a phase diversity-based scheme to approximate the PSF associated with each section of a segmented image. For a sufficiently fine segmentation, the PSF can be assumed to be essentially spatially invariant, and thus a blind deconvolution scheme can be applied. The mathematics of this phase recovery process was first described by Gonsalves [8], and has been applied extensively for imaging through atmospheric turbulence [4].

Assuming that light emanating from the object is incoherent, the dependence of the PSF on the phase is given by

where p denotes the pupil, or aperture, function, ι = √-1, and ℱ denotes the 2-D Fourier transform,

The pupil function p = p(x ₁,x ₂) is determined by the extent of the telescope’s primary mirror. In phase diversity-based blind deconvolution, the k^th diversity image is given by

where η_k represents noise in the data, f is the unknown object, s is the point spread function (PSF), ϕ is the unknown phase function, θ_k is the k^th phase diversity function.

In atmospheric optics [4], the phase ϕ (x ₁,x ₂) quantifies the deviation of the wave front from a reference planar wave front. This deviation is caused by variations in the index of refraction (wave speed) along light ray paths, and is strongly dependent on air temperature. Because of turbulence, the phase varies with time and position in space and is often modelled as a stochastic process.

Additional changes in the phase ϕ can occur after the light is collected by the primary mirror, e.g., when adaptive optics are applied. This involves mechanical corrections obtained with a deformable mirror to restore ϕ to planarity. By placing beam splitters in the light path and modifying the phase differently in each of the resulting paths, one can obtain more independent data. The phase diversity functions θ_k represent these deliberate phase modifications applied after light is collected by the primary mirror. The easiest to implement is defocus blur, modelled by a quadratic

where the parameters b_k are determined by defocus lengths. In practice, the number of diversity images is often quite small, e.g., K = 2 in the numerical simulations to follow. In addition, one of the images, which we will denote using index k = 1, is obtained with no deliberate phase distortion, i.e., θ ₁ = 0 in (2.6).

2.1. The minimization problem

To estimate the phase ϕ, and the object f from data (2.6), we consider the least squares fit-to-data functional

Here ∥ ∙ ∥ denotes the standard L ² norm. By the convolution theorem and the fact that Fourier transforms preserve the L ² norm, one can express this in terms of Fourier transforms, F = ℱ{f}, D_k = ℱ{d_k }, and S_k [ϕ} = ℱ{s[ϕ + θ_k ]},

Since deconvolution and phase retrieval are both ill-posed problems, any minimizer of J_data is unstable with respect to noise in the data. Hence we add regularization terms to obtain the full cost functional,

Here the regularization parameters γ and α are positive real numbers, and the regularization functionals J_object and J_phase provide stability and incorporate prior information.

Because of atmospheric turbulence, variations in the refractive index, and hence the phase itself, can be modelled as a random process [4]. We apply the von Karman turbulence model, which assumes this process is second order, wide sense stationary, and isotropic with zero mean; see, e.g., [4]. It can be characterized by its power spectral density,

where ω = (ω_x ,ω_y ) represents spatial frequency. Corresponding to this stochastic model for phase, we take the phase regularization functional

where 〈f,g〉 = ${\int }_{-\infty }^{\infty }$${\int }_{-\infty }^{\infty }$ f(ω)g ^*(ω)dω, and the superscript ^* denotes complex conjugate. For regularization of the object, we take the “minimal information prior”

Note that the object regularization functional (2.13) is quadratic and the dependence of the fit-to-data functional (2.8) on the object f is quadratic. Moreover, the Hessian with respect to the object of the full cost functional (2.10) is symmetric and positive definite with eigenvalues bounded below by γ. By setting the gradient with respect to the object equal to zero, one obtains a linear equation whose solution yields the object at a minimizer for J_full [10]. From (2.9)–(2.10) and (2.13) one obtains the Fourier representation for the minimizing object,

where

Substituting (2.14) back in (2.9)–(2.10), one obtains the cost functional that we will minimize,

where

See Appendix B of [9] for a detailed derivation.

Our objective is to minimize (2.16) in order to obtain a reconstruction of the phase ϕ. In order to do that, we need an effective computational method. Such a method is discussed in the next section.

2.2. The minimization algorithm

To minimize J in (2.16) we will use, as in [11], a quasi-Newton, or secant, method known as limited memory BFGS (L-BFGS) [12]. A generic quasi-Newton algorithm with line search globalization is presented below. We denote the gradient of J at ϕ by grad J(ϕ) and the Hessian of J at ϕ by Hess J(ϕ).

Quasi-Newton / Line Search Algorithm

In practice, the line search subproblem is solved inexactly [12].

The BFGS method [12] is one of the most popular quasi-Newton techniques. Given Bν, this method generates a new Hessian approximation B _ν+1 in terms of the differences in the successive approximates to the solution and its gradient,

L-BFGS is based on a recursion for the inverse of the B_ν ’s,

Given ν ∈ IRⁿ, computation of $B_{\mathrm{\nu }+1}^{-1}$ ν requires a sequence of inner products involving ν and the s_ν ’s and y_ν ’s, together with the application of $B_0^{-1}$. If B ₀ is symmetric positive definite (SPD) and the “curvature condition” $y_{\nu }^T$s_ν > 0 holds for each ν, then each of the B_ν ’s is also SPD, thereby guaranteeing that - $B_{\nu }^{-1}$ grad J(ϕ_ν ) is a descent direction. The curvature condition can be maintained by implementing the line search correctly [12].

“Limited memory” means that at most N vector pairs {(s_ν ,y_ν ),…, (s _ν-N+1,y _ν-N+1)} are stored and at most N steps of the recursion are taken, i.e., if ν ≥ N, apply the recursion (2.20) for ν, ν - 1,…, ν - N, and set $B_{\nu -N}^{-1}$ equal to an SPD matrix $M_{\nu }^{-1}$. We will refer to M_ν as the preconditioning matrix. In standard implementations, M_v is taken to be a multiple of the identity [12]. For our application we choose an M_ν which has the operator representation on an operand ψ, given by

where Φ is defined in (2.11).

Under mild assumptions, the local convergence rate for the BFGS method is superlinear [12]. In the limited memory case, this rate can become linear.

Before continuing, we make the observation that the above computational approach provides estimates of both the phase ϕ and object f (via (2.14)). In our application of removing spatially-varying blur, the main interest is in using this approach for obtaining the PSF as a function of the phase, as given in (2.4). The object f is then reconstructed in a second stage, which we now discuss.

3. Sectioning the image and applying a global restoration scheme

For many spatially variant blurs, in small regions of the image the blur can be well approximated by a spatially invariant PSF. This property has motivated several types of sectioning methods [2, 13, 14, 15] that partition the image, restoring each local region using its corresponding spatially invariant PSF. The results are then sewn together to obtain the restored image. To reduce blocking artifacts at the region boundaries, larger, overlapping regions are used, and then the restored sections are extracted from their centers. Trussell and Hunt [2] proposed using the Landweber iteration for the local deblurring, and suggested a complicated stopping criteria based on a combination of local and global convergence constraints. Fish, Grochmalicki and Pike [14] use a truncated singular value decomposition (TSVD) to obtain the local restorations.

An alternative scheme can be derived by partitioning the image into subregions on which the blur is assumed to be spatially invariant; but, rather than deblurring the individual subregions locally and then sewing the individual results together, we first sew (interpolate) the individual PSFs, and restore the image globally. A global reconstruction avoids the problem of boundary artifacts that can occur in traditional sectioning methods [2, 13, 14, 15]. In algebraic terms, the blurring matrix S, given in (1.1), can be written as

where S_ij a matrix representing the spatially invariant PSF in region (i, j), and D_ij is a nonnegative diagonal matrix satisfying ∑∑D_ij = I. For example, if piecewise constant interpolation is used, then the lth diagonal entry of D_ij is 1 if the lth point is in region (i, j), and 0 otherwise. In principle, any partitioning scheme can be used. The most obvious approach is to partition the image into rectangular subregions, or to use concentric annular regions in the case that the PSFs vary radially from the on-axis PSF.

Faisal et al. [5] use this formulation of the spatially variant PSF, apply the Richardson-Lucy algorithm with piecewise constant interpolation of the PSFs, and discuss a parallel implementation. Boden et al. [6] also describe a parallel implementation of the Richardson-Lucy algorithm, and consider piecewise constant as well as piecewise linear interpolation. Nagy and O’Leary [7] use a conjugate gradient algorithm with piecewise constant and linear interpolation, and also suggest a preconditioning scheme (for both interpolation methods) that can substantially improve the rate of convergence. Furthermore, it is shown in [16] that by generalizing overlap-add and overlap-save convolution techniques, very efficient FFT-based algorithms can be obtained when the image is partitioned into rectangular subregions. Although other partitioning schemes can be used, the computational cost will increase. If the PSFs vary smoothly, then a rectangular partitioning should provide a good approximation of the spatially variant blur. For these reasons, we assume throughout the paper that a rectangular partitioning of the image is used.

Note that, using (3.21), the image formation model can be written as:

where s_ij is the PSF for region (i, j). This model assumes the PSFs are known, so the question is how to use this for blind deconvolution. The algorithm proceeds follows:

4. Numerical examples

In this section we describe some numerical experiments to illustrate the potential advantages of the space variant approach discussed in this paper. All tests were done on a simulated star field image, shown in Fig. 1, which was obtained from the Space Telescope Science Institute (www.stsci.edu).

 

Fig. 1. Simulated image of a star field.

Download Full Size | PDF

To simulate an image blurred by a spatially variant point spread function, we begin by generating 1024 different PSFs. The PSFs were created by generating a single pupil function and a moving phase screen, each on a 32 × 32 grid. Figure 2 shows the pupil function, two selected phase screens, and their corresponding PSFs.

The blurred image, shown in Fig. 3, was then created by convolving each of the 1024 PSFs with successive 32 × 32 pixels of the true image. By overlapping these regions, we obtain successive 8×8 pixel regions of the blurred image with virtually no blocking (boundary) artifacts; see Fig. 3.

4.1. Advantages of space variant model

To illustrate the potential advantages of using a space variant model, we construct certain average PSFs as follows:

The conjugate gradient (CG) method was used to reconstruct the image with the various approximations of the PSF as described above. Efficient implementation of the matrix vector multiplications for the space variant cases (4 PSFs, 16 PSFs, and 64 PSFs) was done using the approach outlined in Section 3. Goodness of fit is measured using the relative error,

where f_true is the (diffraction limited) true image, and f_k is the (diffraction limited) solution at the kth CG iteration.

A plot of the relative errors for the special case of noise free data is shown in Fig. 4. Note that when more PSFs are used, a better approximation of the space variant blur is obtained. This is confirmed by the results in Fig. 4. The computed reconstructions at iterations corresponding to the minimum relative errors for each case are shown in Fig. 5.

Similar results occur if noise is added to the blurred image, however the ill-conditioning of the problem means that the reconstructions will be more contaminated with noise. For example, with Poisson and Gaussian noise (mean 0, standard deviation 2) added to the blurred image, the relative errors for CG are shown in Fig. 6 and the corresponding reconstructions are shown in Fig. 7.

 

Fig. 4. Relative errors at each iteration of the conjugate gradient method, using increasingly more accurate approximations of the spatially variant blur. These results were computed using a noise free blurred image.

Download Full Size | PDF

4.2. Practical results

The results from the previous section clearly show that if we are able to get accurate approximations of the PSFs, then the space variant model will provide better reconstructions than by using a space invariant model. In this section we show results when the PSFs are reconstructed using a blind deconvolution algorithm. In particular, we use phase diversity, with one frame and two channels. The blind deconvolution (BD) algorithm we use is described in Section (2.2). Our regularization parameters were taken to be α= 1×10^-1 and γ= 1×10^-4. We remark that the problem of choosing appropriate regularization parameters is one of the most difficult issues in the numerical solution of ill-posed problems. There are techniques for estimating parameters for simple noise models; see, for example [17, 18, 19]. In many cases good heuristics based on computational experience with the algorithm and the particular application can be useful. The parameters we use are similar to those reported in the literature; see [11].

 

Fig. 5. Reconstructions computed by the conjugate gradient algorithm, using accurate approximations of the PSF, and with a noise free blurred image. For comparison, we also show the true image and the true image convolved with the diffraction-limited PSF.

Download Full Size | PDF

 

Fig. 6. Relative errors at each iteration of the conjugate gradient method, using increasingly more accurate approximations of the spatially variant blur. These results were computed using a noisy (Poisson and Gaussian) blurred image.

Download Full Size | PDF

 

Fig. 7. Reconstructions computed by the conjugate gradient algorithm, using accurate approximations of the PSF, with Poisson and Gaussian noise added to the blurred image. For comparison, we also show the true image and the true image convolved with the diffraction-limited PSF.

Download Full Size | PDF

In the phase estimation step, regularization is applied before computing a minimum of the least squares functional. Thus it is reasonable to seek estimates of the phase and object that yield a gradient norm near zero. We therefore stop L-BFGS iterations once the norm of the gradient had been reduced by six orders of magnitude in all cases. Note that no noise amplification occurs provided the regularization parameters are properly chosen, and, as previously mentioned, the choice of these parameters is dependent on the data.

Once the PSFs are computed from the BD algorithm, we then use the CG algorithm to reconstruct the image, as was done in the previous subsection. Figure 9 shows plots of the relative errors when no noise is added to the data. At first glance these appear to be disappointing results; 4 PSFs produce lower relative errors than a single PSF, but additional PSFs actually increase the relative error. It is important to note, however, that to obtain more PSFs, the blind deconvolution algorithm must be implemented on small regions of the image. If the small regions do not contain significant object information, then we cannot hope to obtain good approximations of the PSFs. Consequently, given the fact that the object that we are using in our experiment has significant black background, it is not surprising that several such small regions occur. Thus the poor results when using 64 PSFs.

On the other hand, some of the small regions may contain enough significant object information so that the blind deconvolution algorithm can reconstruct good PSFs. This motivates us to consider an approach where we refine the partitioning in each subregion only if further refinement produces good reconstructions of the PSFs. For example:

Thus, this adaptive approach uses a “mix” of PSFs computed by the various partitionings of the image.

 

Fig. 8. Example of adaptive refinement of region partitions. The left plot shows an initial partitioning of an image, the middle shows a refinement of the partitioning, and the right plot shows what might happen in an adaptive approach where only certain subregions are refined further.

Download Full Size | PDF

To automate the process of determining when a PSF is reconstructed sufficiently accurately, an appropriate measure must be used. Such a measure will be dependent on the application, and the type of images being reconstructed. In order to test the potential of this adaptive scheme, we used a systematic approach for choosing a good PSF by comparing the computed PSFs with the average PSFs described in the previous subsection. Of course in a realistic problem the average PSFs will not be available, and some other mechanism must be used to decide on the quality of the computed PSFs. In any case, our systematic approach shows that if the best computed PSFs can be found, then substantially better reconstructions can be obtained; see the bottom curve in Fig. 9. The reconstructed images are shown in Fig. 10.

Similar results are obtained with noisy data. For example, with Poisson and Gaussian noise (mean 0, standard deviation 2) added to the blurred image, the relative errors using the CG iterative method are shown in Fig. 11. The corresponding reconstructions are shown in Fig. 12, where we also include arrows to indicate regions in which the space variant approach produces significantly better reconstructions (compared to the diffraction limited true image) than the space invariant (1 PSF) model.

 

Fig. 9. Relative errors at each iteration of the conjugate gradient method, using increasingly more accurate approximations of the spatially variant blur. The PSFs used to approximate the space variant blur were computed using phase diversity-based blind deconvolution. These results were computed using noise free blurred image data.

Download Full Size | PDF

5. Conclusions

We have presented a computational approach for solving image restoration problems in which the blur in the image is both unknown and spatially varying. The approach has three stages. The first stage involves sectioning the image into regions in which the blur is believed to be approximately spatially invariant. In the second stage, phase diversity-based blind deconvolution via the L-BFGS optimization algorithm is implemented in order to obtain an estimate of the phase in each region. From these reconstructed phases, the corresponding PSF in each region can be computed. In the final stage, with these PSFs in hand, the object can be reconstructed globally via the algorithm of Nagy and O’Leary.

Our numerical experiments show, first, that using a spatially varying model when a spatially varying blur is present does indeed provide more accurate results. Secondly, we find that in regions with little object information, the phase, and hence, PSF, reconstructions can be inaccurate. This motivates a “PSF mixing” scheme, in which the object is divided into further subregions only in areas in which there is enough object information. A conclusion of particular importance that follows from our numerical experiments is that in the presence of an unknown, spatially varying blur, our approach is much more effective than is standard one-PSF phase-diversity.

We remark that, in principle, our space variant approach should be applicable with other blind deconvolution methods, but as evidenced from our numerical results, the quality of the reconstructed PSFs is important. Thus we would expect that additional diversity information should improve the results. Similarly, we would expect a multi-frame scheme, coupled with our space variant technique, to outperform its single frame counterpart.

 

Fig. 10. Reconstructions computed by the conjugate gradient algorithm, using PSFs computed from a phase diversity blind deconvolution algorithm. The blurred image data in this case is noise free.

Download Full Size | PDF

 

Fig. 11. Relative errors at each iteration of the conjugate gradient method, using increasingly more accurate approximations of the spatially variant blur, and with a noisy (Poisson and Gaussian) blurred image.

Download Full Size | PDF

 

Fig. 12. Reconstructions computed by the conjugate gradient algorithm, using PSFs computed from a phase diversity-based blind deconvolution algorithm, and with a noisy (Poisson and Gaussian) blurred image.

Download Full Size | PDF

The problem of optimal sectioning of the image, and efficient implementation details, is worth further exploration. One might base the partitioning on a-priori knowledge of how the PSF varies; for example, one could partition into concentric annular regions for a radially varying PSF. Another approach is to base the partitioning on the amount of information (e.g., signal to noise ratio) in various regions. In this case, one could use an approach where regions are repartitiioned into smaller and smaller subregions until, for example, the mean intensity of the subregion falls below a predetermined threshold (which might be obtained through extensive simulations).