A&A 399, 807-811 (2003)
DOI: 10.1051/0004-6361:20021872

Restoration of astronomical images by an iterative superresolution algorithm

E. Boukouvala - A. H. Lettington

University of Reading, Department of Physics, JJ Thomson Physical Laboratory, PO Box 220, Reading RG6 6AF, UK

Received 21 April 2000 / Accepted 10 September 2002

Abstract
In this work a new fast method for super-resolving noisy images is described. The iterative technique is similar to the positive constraint algorithm but faster and more effective in terms of the resulting global noise in the image. The astronomical use of the technique is illustrated by applying it to the LMC 49-50 (04/01/98) image. Results from this method are compared with those of well-known algorithms such as the Wiener filter, the positive constraint, the maximum a posteriori probability, and the Lucy-Richardson method.

Key words: methods: data analysis - techniques: image processing


1 Introduction

In incoherent optical and radio-astronomical imaging the observed image illumination g(i) can be given by

 \begin{displaymath}% g(i) = h(i) \otimes f(i) + n(i) \end{displaymath} (1)

where f(i) is the illumination of an ideal image, h(i) is the point spread function (psf) of the system and n(i) represents the additive noise. The symbol $\otimes$ represents a convolution. If the psf is shift invariant and a circular convolution is allowed, (1) may be written in the Fourier domain as

 
G(k) = H(k) . F(k) + N(k) (2)

where the capital letters represent the Fourier transforms of the lower case counterparts and where . represents a scalar product. The quantity k is a specific spatial frequency within the M/2 sampled frequencies arising from the M sampled picture points. When an imager is band-limited the range of H, the optical transfer function, is less than M/2 while F is assumed to contain the full frequency spectrum.


  \begin{figure} \par\includegraphics[width=15.5cm,clip]{MS1883f1.eps}\end{figure} Figure 1: a) Original image, b) blurred image, c) linearly restored within the passband image, d) 500th iterations estimate (plotted as a function of position in the scene).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=15.5cm,clip]{MS1883f2.eps}\end{figure} Figure 2: Spectrum of a) original image, b) blurred image, c) linearly restored image and d) 500th iterations estimate.
Open with DEXTER

In linear restoration of band-limited images, spatial frequencies may be restored within the pass-band of the optics, through a process of linear deconvolution of the imager psf. This however introduces Gibbs ringing as well as enhancement of the high frequency components of the noise. Using a Wiener filter may reduce these effects but this additionally reduces the amount of restoration that can be achieved. Therefore, a number of non-linear algorithms have been introduced that attempt to suppress the Gibbs ringing through the addition of high spatial frequencies so achieving bandwidth extension or "super-resolution''.

Richardson (Richardson 1972) and Lucy (Lucy 1974) proposed the well-known Lucy-Richardson (LR) algorithm, a method of estimation based upon Bayes theorem. The same theorem formed the basis for the maximum a posteriori probability technique (MAP) proposed by Hunt (Hunt 1977).

The maximum-entropy method (ME) introduced by Jaynes (Jaynes 1968) and Frieden (Frieden 1972) has become widely used for astronomical image restoration. However, this technique is sensitive to noise and slow when compared with other methods. Later, researchers modified the technique in order to make the method more tolerant to noise and computationally simpler. In the present work, we have used the Gonsalves' and Kao's (Gonsalves & Kao 1987) ME algorithm when making comparisons.

The energy-error reduction algorithm proposed by Gerchberg (Gerchberg 1974) has a varying efficiency depending on the a priori knowledge of the extent of the finite object whose spectrum is to be continued. This technique applies constraints to the image in the object and the Fourier domains alternately. The diffraction-limited image of the localized object is set to zero outside the known spatial extent of the object. Then in the Fourier domain the calculated spectrum within the passband is replaced by the known low frequency spectrum. The corrected spectrum is then inverse transformed to the object domain to repeat the iterative process until a suitably restored image is obtained.

Another method similar to the Gerchberg technique is the positive constraint method (PC), proposed by Biraud (Biraud 1969), which may be applied to images where the intensity values are always positive. For images which contain intensities near zero, it is likely that negative intensities will be introduced, following inverse filtering, due to the lack of high spatial frequencies. Applying a positive constraint to such an image generates high frequencies and corrupts the low. These low frequencies are replaced by the already known original values and the process is iterated until the image is effectively restored.

In the present work a variant of the positive constraint method is proposed which works efficiently on astronomical images while also being fast and tolerant to noise.

2 Mathematical description of the algorithm

The process begins by constructing a Wiener filter solution $f_{{\rm Wiener}}(i)$ which is obtained from its Fourier transform $F_{{\rm Wiener}}(k)$ given by (3)

 \begin{displaymath}% F_{{\rm Wiener}}(k) = \frac{G(k)H^{\star}(k)}{ \vert H(k) \vert ^{2} + z}\cdot \end{displaymath} (3)

Here z corresponds to the ratio of the noise power spectrum to the image power spectrum and $H^{\star}(k)$ is the complex conjugate of the optical transfer function H(k). This Wiener filter solution provides the starting values for fn(i) and Fn(k) used in the following iterative procedure.

This procedure starts by computing the quantity abs(fn), which symbolizes the absolute value of fn. Next the Fourier transform Fn+1 of the updated image fn+1 is obtained from the following equation which includes replacing the known low spatial frequencies.

 
Fn+1(k)=Fn(k) . H(k) + (1-H(k)) . ABSfn(k) (4)

where

 \begin{displaymath}% ABSf_{n}(k) = \mathcal{F} (abs(f_{n}(i)) \end{displaymath} (5)

$\mathcal{F}$ symbolizing the Fourier transform and ABSfn the Fourier transform of abs(fn).

Then taking the inverse Fourier transform of Fn+1

 \begin{displaymath}% f_{n+1}(i)= \mathcal{F} ^{-1}(F_{n+1}(k)) \end{displaymath} (6)

provides the new estimated image and the whole iterative procedure starts again. Because of its simplicity the iterative procedure is fast and offers the prospect of non-linear super-resolution processing in real time to images where a positive constraint may be applied.

3 Application of the technique to simulated images

To evaluate the performance of the proposed algorithm, it has been compared with results obtained from other more established techniques. A one-dimensional simulated image consisting of two blurred, closely spaced spikes was restored using the present method and the PC (Biraud 1969), LR (Lucy 1974; Richardson 1972), the MAP (Hunt 1977) methods.

The synthetic image was a linear array of 128 picture points containing two spikes of length 100 at positions 69 and 72 in the array. This array was then convolved with a normalized Gaussian function having a standard deviation of 1.5 picture points. The separation of the objects is then 0.85 of the full width at half maximum (FWHM) of the psf which satisfies the criterion suggested by Lucy (1992a,b), that the separation of the stars should be less than the FWHM of the psf. In the present case 5% random white noise, having a mean value of zero, was then added to this blurred image. The original and blurred images are illustrated as Figs. 1a and 1b respectively. The two spikes are unresolved in the blurred image which was the starting point for each of the restoration methods.

In the present method this blurred image was first restored using a Wiener filter with a z value of 0.0001, as illustrated in Fig. 1c. The results of using the present algorithm are illustrated in Fig. 1d after 500 iterations.

The spectrum of the restored image is given in Fig. 2d. It can be seen that previously lost high frequency information has been effectively recovered after 500 iterations, even in the presence of 5% noise. The result in Fig. 2d had converged to within 99% after 500 iterations. Moreover the suggested algorithm contains no adjustable parameter which needs optimization while converging, although the initial background level needs to be estimated precisely, a problem discussed by Bijaoui (1980).

The variation in the signal-to-noise ratio (SNR) as a function of the number of iterations is illustrated in Fig. 3 together with the results calculated for the other methods. It can be seen that the present method gave by far the best results. Values of the global root-mean-square (rms) noise levels contained in the images are given in Table 1 after 200 and 500 iterations, for the various methods discussed. It can be seen that the present method compares favourably with the others.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{MS1883f3.eps}\end{figure} Figure 3: SNR as a function of the number of iterations for present algorithm (* star line), PC (solid line), MAP (. dotted line) and LR (o circle line).
Open with DEXTER


 

 
Table 1: Results on the global noise (rms value of noise).
Number of it. Present algor. PC MAP LR
200 0.64 0.83 0.69 0.82
500 0.45 0.58 0.48 0.55



 

 
Table 2: Saddle-to-peak intensity ratios (for 500 iter.).
Results on R1 and R2 Present algor. PC MAP LR
R1 0.003 0.004 0.122 0.586
R2 0.996 0.910 0.984 0.988


Another parameter compared in Table 2 is the saddle-to-peak intensity ratio in the resulting images. According to Lucy's work on the resolution criteria for the deconvolved images (Lucy 1992b), an algorithm is considered to have resolved the double star when R1[= 2I2/(I1+I3)]<0.735 and R2[= I1/I3]>2/3. Here $I_{1}({\rm max})$, $I_{2}({\rm min})$, $I_{3}({\rm max})$ (with I1<I3), are the double star relative intensities.

Values of R1 and R2 obtained using the various algorithms are also included in Table 2. It can be seen that Lucy's criteria of resolution are satisfied in all cases but that the present algorithm gave the best results. It was found that the MAP and the LR were very sensitive to the actual noise included in the image and the results varied from run to run. Similar results were observed with the ME method which performed like the LR but not as well as the MAP algorithm. The results for ME were omitted in the interest of clarity.

A similar trend was observed in the restoration of two-dimensional images as described in the following section.

4 Application of the technique to astronomical images

In order to test our algorithm we used an LMC 49-50 image recorded by the Danish 1.54 m telescope on 4th January 1998 using a yellow filter. The DFOSC camera was used as the detector array and the exposure time was 1200 s. Before processing, the image was flat-field corrected and bias corrected using the procedures of the IRAF package. The image contained 256 by 256 picture points and was in bitmap (bmp) format.

The psf of the image was estimated by applying a least-square fitting method to a number of suitable stars. The resulting profile was best fitted with a two-dimensional Gaussian function having a standard deviation of 2.5 picture points.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{MS1883f4.eps}\end{figure} Figure 4: 3-dimensional plot of the original image where intensities are presented as a function of position.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{MS1883f5.eps}\end{figure} Figure 5: 3-dimensional plot of image resulting from application of the proposed algorithm after 3 iter.
Open with DEXTER

Figure 4 shows a 3-dimensional plot of the original image where the intensities of the stars are presented as a function of position. The application of the present super-resolving algorithm is illustrated in Fig. 5 after 3 iterations. It can be seen that even after so few iterations, the present algorithm is effective in sharpening the quasi-point sources as well as in revealing data previously hidden in noise. These results converged, after 30 iterations, with a further improvement of the image quality.

5 Discussion and conclusions

The proposed algorithm starts with a Wiener filter solution in which ringing occurs due to boosting the high frequency content of the noise and to the lack of higher spatial frequencies in the image beyond the pass band. The effect of taking an absolute value of this image doubles the frequency of the ringing and of the noise at each iteration, ultimately converting it to dc where it is removed in the low frequency substitution.

The present algorithm has been successful in recovering high spatial frequency information and is very powerful on the basis of its resolving power. It is similar to the positive constraint technique of Biraud but takes an absolute value instead of cutting negative-going intensities to zero. It is very tolerant of noise and can therefore be applied successfully to astronomical images and to other positive images with well-defined baselines. In those applications where the image to be restored is superimposed on a non-zero background then other techniques must be applied to subtract this background before the present method and many other methods can be successfully applied. It is very effective and faster in comparison with many other super-resolution algorithms and offers the prospect of non-linear super-resolution image processing in real-time.

A copy of the present algorithm in a form of the Matlab program is available on request from the authors.

Acknowledgements
The authors wish to thank Dr. A. Dapergolas and Dr. I. Georgantopoulos, researchers in the National Observatory of Athens for supplying the LMC 49-50 image and for the useful discussions. Also the first author wishes to express her particular thanks to Dr. M. Kontizas for supervision and the provision of resources.

References

 

Copyright ESO 2003