## Deconvolution of astronomical images using SOR with adaptive relaxation |

Optics Express, Vol. 19, Issue 14, pp. 13509-13524 (2011)

http://dx.doi.org/10.1364/OE.19.013509

Acrobat PDF (1867 KB)

### Abstract

We address the potential performance of the successive overrelaxation technique (SOR) in image deconvolution, focusing our attention on the restoration of astronomical images distorted by atmospheric turbulence. SOR is the classical Gauss-Seidel iteration, supplemented with relaxation. As indicated by earlier work, the convergence properties of SOR, and its ultimate performance in the deconvolution of blurred and noisy images, can be made competitive to other iterative techniques, including conjugate gradients, by a proper choice of the relaxation parameter. The question of how to choose the relaxation parameter, however, remained open, and in the practical work one had to rely on experimentation. In this paper, using constructive (rather than exact) arguments, we suggest a simple strategy for choosing the relaxation parameter and for updating its value in consecutive iterations to optimize the performance of the SOR algorithm (and its positivity-constrained version, +SOR) at finite iteration counts. We suggest an extension of the algorithm to the notoriously difficult problem of “blind” deconvolution, where both the true object and the point-spread function have to be recovered from the blurred image. We report the results of numerical inversions with artificial and real data, where the algorithm is compared with techniques based on conjugate gradients. In all of our experiments +SOR provides the highest quality results. In addition +SOR is found to be able to detect moderately small changes in the true object between separate data frames: an important quality for multi-frame blind deconvolution where stationarity of the object is a necesessity.

© 2011 OSA

## 1. Introduction

*a priori*constraints on the solutions, like non-negativity which is particularly important in astronomical imaging. However, a problem which is common for the iterative techniques is their slow convergence, in particular when “blind” deconvolution is addressed. This provides the motivation for a search for new algorithms. In a recent paper [1

1. V. N. Strakhov and S. V. Vorontsov, “Digital image deblurring with SOR,” Inverse Probl. **24**, 025024 (2008) (Paper I). [CrossRef]

2. K. Sauer and C. Bouman, “A local update strategy for iterative reconstruction from projections,” IEEE Trans. Signal Process. **41**, 534–548 (1993). [CrossRef]

3. J. A. Fessler, “Penalized weighted least-squares image reconstruction for positron emission tomography,” IEEE Trans. Med. Imaging **13**, 290–300 (1994). [CrossRef] [PubMed]

*x*of size

*N*

^{2}represents the

*N*×

*N*true image, with image elements arranged by lexicographical ordering, the

*N*

^{2}×

*N*

^{2}matrix

*A*represents a blurring operator, and the “true” data

*f*is distorted by noise

*δf*. We assume that the

*N*×

*N*domain is large enough to accommodate both the true object and its blurred image, so that edge effects can be discarded and matrix

*A*can be considered as block-circulant with circulant blocks; the central column of

*A*represents lexicographically-ordered elements of the point-spread function (PSF).

*A*in the Eq. (3) is represented by the matrix sum where

^{T}A*L*is a strictly lower triangular matrix and

*D*is a diagonal matrix, and each consecutive iteration consists in solving the triangular linear system starting with an initial guess

*x*

^{(0)}, where

*τ*is the relaxation parameter, which we will later allow to change between the consecutive iterations. When

_{k}*τ*is set to 1, SOR reduces to the classical Gauss-Seidel iteration.

_{k}*A*is a convolution operator, matrix

*L*also represents a convolution (with a point-spread function which can be obtained from SPSF by setting to zero the elements of the central row starting from the central element and progressing to the right, together with all the elements below the central row). Matrix

*D*can be written as

*D*=

*dI*, where

*I*is the unit matrix and

*d*is the square of the Euclidean norm of the PSF. Therefore, the matrix-vector multiplications in Eq. (5) are convolutions, and computational expenses can be reduced by solving this equation in the Fourier domain.

1. V. N. Strakhov and S. V. Vorontsov, “Digital image deblurring with SOR,” Inverse Probl. **24**, 025024 (2008) (Paper I). [CrossRef]

4. A. Björck, *Numerical Methods for Least Squares Problems* (SIAM, 1996). [CrossRef]

*N*

^{2}minor steps where the merit function (Eq. (2)) is minimized separately with respect to each individual component

*x*of the vector

_{i}*x*in the consecutive order. At each minor step, an optimal variation of

*x*is scaled with the relaxation parameter

_{i}*τ*. This algorithm avoids the calculation of the SPSF, and allows us to work with a spatially-variant PSF, and/or with quadratic approximation for Poissonian likelihood in the merit function.

_{k}5. C. W. Cryer, “The solution of a quadratic programming problem using systematic overrelaxation,” SIAM J. Control **9**, 385–392 (1971). [CrossRef]

*a priori*information is the support-domain constraint. It is particularly important in “blind” deconvolution, where the true object needs to have a finite support for its deconvolution to be possible without prior knowledge of the point-spread function [6

6. R. G. Lane and R. H. T. Bates, “Automatic multidimensional deconvolution,” J. Opt. Soc. Am. **4**, 180–188 (1987). [CrossRef]

7. R. T. H. Bates, B. K. Quek, and C. R. Parker, “Some implications of zero sheets for blind deconvolution and phase retrieval,” J. Opt. Soc. Am. **7**, 468–479 (1990). [CrossRef]

## 2. Scheduling the relaxation parameter

*r*

^{(k)}|| = ||

*r*

^{(k−1)}||) when

*τ*∈ (0, 2], since

_{k}*D*is positive-definite. We now address the reduction of the residual vector in a different norm (in the

*AA*-norm). Taking the Euclidian norm of both sides of the equation it is straightforward to show that The matrices

^{T}*L*and

*L*represent convolutions, and hence

^{T}*LL*–

^{T}*L*= 0, when edge effects are discarded, since any two convolutions commute. Matrix

^{T}L*D*is

*D*=

*dI*, where

*I*is the unit matrix and

*d*is the central element of the SPSF (sum of squares of all the elements of the PSF). We thus have or, using Eq. (9), which proves that the residual is also reduced in the

*A*-norm, in the same range of

^{T}A*τ*. The last equation gives also an interesting hint. Suppose for a moment that the residual vector

_{k}*r*

^{(k)}can be reduced to zero (or almost zero) in one iteration. Setting

*r*

^{(k)}to zero in the Eq. (14) shows that the required value of the relaxation parameter

*τ*shall then be calculated from

_{k}*r*

^{(k−1)}using the simple equation

*τ*in two limiting situations, when the complete convergence can be achieved in one iteration. One is when the point-spread function differs from zero at one pixel only. We then have

_{k}*A*=

^{T}A*AA*=

^{T}*D*=

*dI*, and hence ||

*A*

^{T}r^{(k−1)}||

^{2}=

*d*||

*r*

^{(k−1)}||

^{2}, which gives

*τ*= 1. Since

_{k}*L*= 0, Eq. (7) gives

*A*Δ

^{T}A*x*=

*A*

^{T}r^{(k−1)}, and hence

*r*

^{(k)}= 0 from Eq. (11). Another limiting situation is when the PSF spreads over more than one pixel, but the residual image contains very low spatial frequencies only, so that the PSF works again like Dirac’s

*δ*-function. When the sum of the PSF’s elements is normalized to 1 (which is a standard practice, since a PSF is usually interpreted as a probability distribution function for a single photon), we have

*AA*

^{T}r^{(k−1)}=

*r*

^{(k−1)}, and the Eq. (15) suggests

*τ*= 2

_{k}*d/*(1 +

*d*). This value is the same as that suggested for the fastest restoration of the lowest frequencies by the convergence analysis of Paper I.

*ℱ*to designate a 2-D discrete Fourier transform (DFT), we have

*ℱ*(

*Lx*) =

*lℱ*(

*x*), where

*l*is DFT of a truncated point-spread function obtained from SPSF by setting to zero the elements of the central row starting with the central element and progressing to the right, together with all the elements below the central row. We also have

*ℱ*(

*L*) =

^{T}x*l*(

^{*}ℱ*x*), where

*l*is the complex conjugate of

^{*}*l*, and

*ℱ*(

*Dx*) =

*dℱ*(

*x*). Let also

*ℱ*(

*Ax*) =

*aℱ*(

*x*)

*, ℱ*(

*A*

^{T}*x*) =

*a*(

^{*}ℱ*x*), such that

*a*=

^{*}a*l*+

*d*+

*l*is DFT of the SPSF.

^{*}*designates the summation over all the spatial frequencies. The right-hand side of Eq. (16) is thus a weighted average of*

_{ω}*l*+

*d*+

*l*over the spatial frequencies; the weighting function is the Fourier power of the current residuals.

^{*}*l*) is imaginary part of

*l*.

*r*

^{(k−1)}contains only one Fourier component, or is dominated by the spectral components limited by a small frequency domain. The last equation shows that the convergence in one iteration is possible when Im(

*l*) = 0, if the relaxation parameter

*τ*is taken such that (2/

_{k}*τ*− 1)

_{k}*d*=

*l*+

*d*+

*l*, in agreement with the recipe suggested by Eq. (16). When Im(

^{*}*l*) differs from zero, the value suggested for

*τ*is nearly optimal (an optimal reduction of |

_{k}*W*(

*τ*)| is achieved with (2/

_{k}*τ*− 1)

_{k}*d*= |2

*l*+

*d*|).

*r*

^{(k−1)}spreads over a wide range of spectral components, the weighted average

*l*+

*d*+

*l*suggested by the Eq. (16) for (2/

^{*}*τ*− 1)

_{k}*d*will provide a value for the relaxation parameter which is more favourable for suppressing those spectral components in the residual which have bigger weights |

*ℱ*(

*r*

^{(k−1)}) |

^{2}, i.e. those which are bigger in power (note that

*d*does not depend on frequency). As iterations proceed, the power spectrum of the residuals will evolve (towards smaller values) with a general tendency of becoming more uniform in frequency. The flattening of the power spectrum of the approximation residuals during the iterative descent is the property which is most favourable for obtaining a high-quality restoration of the blurred and noisy image.

*τ*, close to its low-frequency limit 2

_{k}*d*/(1 +

*d*); this iteration will reduce the components with the lowest frequencies in the residuals, with almost no impact on the high frequencies. At further iterations, we expect

*τ*to increase gradually thus reducing the higher-frequency components, which now have larger relative weights in the residuals, with the net effect of flattening the Fourier power of the residuals towards the noise level. The iterations shall be terminated when the Fourier power of the residuals is made nearly flat, as further iterations will start to transfer the magnified high-frequency noise from the residuals into the solution. If the iterations are continued,

*τ*will continue growing, but we expect the rate of change of

_{k}*τ*with

_{k}*k*to slow down at high iterations. This is because at very high frequencies (it depends on the PSF what “very high” really means) the residuals can not be quickly suppressed (the effect of the imaginary part of

*l*in Eq. (20) becomes more significant).

*x–x*

_{true}||/||

*x*

_{true}|| versus iteration count, is shown by the solid line in Fig. 2. The adapted values of the relaxation parameter

*τ*are shown along this line. As expected, the iterations start with

_{k}*τ*close to the theoretical low-frequency limit 2

_{k}*d*/(1 +

*d*), which is 0.00039 in this experiment. As also expected, the relaxation parameter grows with iteration count

*k*, from the small values (deep underrelaxation) to

*τ*> 1 (overrelaxation), and the growth rate slows down at higher

_{k}*k*. The solution with minimum error (0.3294) is reached at iteration count

*k*= 13; this solution is shown in Figure 1(c). At higher iterations, the solution degrades because of the magnified high-frequency noise which is transferred into the solution.

*et al*[8

8. R. C. Puetter, T. R. Gosnell, and A. Yahil, “Digital image reconstruction: deblurring and denoising,” Ann. Rev. Astron. Astrophys. **43**, 139–194 (2005). [CrossRef]

*x*–

*x*

_{true}||/||

*x*

_{true}|| = 0.3299. We observe that the performance of the adaptive +SOR is much higher in this experiment: the minimum-error solution is more accurate, and the number of iterations required to achieve this solution is an order of magnitude smaller. The computational cost of a single iteration, however, is significantly larger: enforcing the positivity constraint at each sequential pixel update does not allow any reduction in the numerical expenses by excursions to the Fourier domain.

*x*–

*x*

_{true}||/||

*x*

_{true}|| = 0.3406 [9

9. D. Calvetti, G. Landi, L. Reichel, and F. Sgallari, “Non-negativity and iterative methods for ill-posed problems,” Inverse Probl. **20**, 1747–1758 (2004). [CrossRef]

10. P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. **34**, 561–580 (1992). [CrossRef]

11. C. R. Vogel, *Computational Methods for Inverse Problems* (SIAM, 2002). [CrossRef]

## 3. Blind deconvolution

12. C. L. Matson, K. Borelli, S. M. Jefferies, E. K. Hege, C. C. Beckner, and M. Lloyd-Hart, “A fast and optimal multiframe blind deconvolution algorithm for high-resolution, ground-based imaging of space objects,” Appl. Opt. **48**, A75–A92 (2009). [CrossRef]

^{−4}level (see Fig. 5) is hardly needed for working with astronomical data, since the noise level is usually much higher. However, it is important to make sure that the quality of the restored images is only limited by data quality, not by the efficiency of the inversion technique.

*δ*-function for the object, input data for the PSFs). The best restorations correspond, roughly, to the turning points in the convergence-history curves. At higher iterations, the restored images exhibit growing contamination by the high-frequency noise. In the +SOR results, the optimal iteration count can also be identified as near-flattening of the Fourier power of the approximation residuals (as in Fig. 3; we did not address the residuals of the +CG inversion).

## 4. Discussion

## Acknowledgments

## References and links

1. | V. N. Strakhov and S. V. Vorontsov, “Digital image deblurring with SOR,” Inverse Probl. |

2. | K. Sauer and C. Bouman, “A local update strategy for iterative reconstruction from projections,” IEEE Trans. Signal Process. |

3. | J. A. Fessler, “Penalized weighted least-squares image reconstruction for positron emission tomography,” IEEE Trans. Med. Imaging |

4. | A. Björck, |

5. | C. W. Cryer, “The solution of a quadratic programming problem using systematic overrelaxation,” SIAM J. Control |

6. | R. G. Lane and R. H. T. Bates, “Automatic multidimensional deconvolution,” J. Opt. Soc. Am. |

7. | R. T. H. Bates, B. K. Quek, and C. R. Parker, “Some implications of zero sheets for blind deconvolution and phase retrieval,” J. Opt. Soc. Am. |

8. | R. C. Puetter, T. R. Gosnell, and A. Yahil, “Digital image reconstruction: deblurring and denoising,” Ann. Rev. Astron. Astrophys. |

9. | D. Calvetti, G. Landi, L. Reichel, and F. Sgallari, “Non-negativity and iterative methods for ill-posed problems,” Inverse Probl. |

10. | P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. |

11. | C. R. Vogel, |

12. | C. L. Matson, K. Borelli, S. M. Jefferies, E. K. Hege, C. C. Beckner, and M. Lloyd-Hart, “A fast and optimal multiframe blind deconvolution algorithm for high-resolution, ground-based imaging of space objects,” Appl. Opt. |

**OCIS Codes**

(100.0100) Image processing : Image processing

(100.3020) Image processing : Image reconstruction-restoration

(100.3190) Image processing : Inverse problems

(100.1455) Image processing : Blind deconvolution

**ToC Category:**

Image Processing

**History**

Original Manuscript: April 22, 2011

Revised Manuscript: May 31, 2011

Manuscript Accepted: June 8, 2011

Published: June 28, 2011

**Citation**

S. V. Vorontsov, V. N. Strakhov, S. M. Jefferies, and K. J. Borelli, "Deconvolution of astronomical images using SOR with adaptive relaxation," Opt. Express **19**, 13509-13524 (2011)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-14-13509

Sort: Year | Journal | Reset

### References

- V. N. Strakhov and S. V. Vorontsov, “Digital image deblurring with SOR,” Inverse Probl. 24, 025024 (2008) (Paper I). [CrossRef]
- K. Sauer and C. Bouman, “A local update strategy for iterative reconstruction from projections,” IEEE Trans. Signal Process. 41, 534–548 (1993). [CrossRef]
- J. A. Fessler, “Penalized weighted least-squares image reconstruction for positron emission tomography,” IEEE Trans. Med. Imaging 13, 290–300 (1994). [CrossRef] [PubMed]
- A. Björck, Numerical Methods for Least Squares Problems (SIAM, 1996). [CrossRef]
- C. W. Cryer, “The solution of a quadratic programming problem using systematic overrelaxation,” SIAM J. Control 9, 385–392 (1971). [CrossRef]
- R. G. Lane and R. H. T. Bates, “Automatic multidimensional deconvolution,” J. Opt. Soc. Am. 4, 180–188 (1987). [CrossRef]
- R. T. H. Bates, B. K. Quek, and C. R. Parker, “Some implications of zero sheets for blind deconvolution and phase retrieval,” J. Opt. Soc. Am. 7, 468–479 (1990). [CrossRef]
- R. C. Puetter, T. R. Gosnell, and A. Yahil, “Digital image reconstruction: deblurring and denoising,” Ann. Rev. Astron. Astrophys. 43, 139–194 (2005). [CrossRef]
- D. Calvetti, G. Landi, L. Reichel, and F. Sgallari, “Non-negativity and iterative methods for ill-posed problems,” Inverse Probl. 20, 1747–1758 (2004). [CrossRef]
- P. C. Hansen, “Analysis of discrete ill-posed problems by means of the L-curve,” SIAM Rev. 34, 561–580 (1992). [CrossRef]
- C. R. Vogel, Computational Methods for Inverse Problems (SIAM, 2002). [CrossRef]
- C. L. Matson, K. Borelli, S. M. Jefferies, E. K. Hege, C. C. Beckner, and M. Lloyd-Hart, “A fast and optimal multiframe blind deconvolution algorithm for high-resolution, ground-based imaging of space objects,” Appl. Opt. 48, A75–A92 (2009). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.