3 Determination of the level of regularization by comparison with the Wiener filter

As indicated in the previous sections, the reconstruction of the object must be made using regularization, in order to prevent unwanted amplification of the noise. We consider here that the regularization is obtained by stopping the iterations. The determination of the optimal iteration number (or of the regularization coefficient) is one of the main problems encountered in the reduction of image blurring. Several techniques have been suggested for that: feasibility tests and Cross-Validation tests (Nuñez & Llacer 1993), or Goodness of fit tests (Lucy 1994).

Figure 9: Top: raw results of the Wiener filter for the 104 photons/image (left) and the 106 photons/image (right). These images are then clipped to emphasize negative parts (bottom): about 45% of the pixels are negative (represented in black). Here also only the central parts of the image are given

We describe here the implementation of a new objective-stopping technique that is based on the comparison of the results of these algorithms (RLA or ISRA) with that of the optimum Wiener filter. The principle of the approach was recently proposed in an ESO/NOAO Conference (Lanteri et al. 1998). Fundamentally, it consists of selecting the iteration number k for ISRA (or RLA) that gives a result close to that of the Wiener filter. This technique could also be used to determine the optimal regularization factor when an explicit regularization is performed.

The idea is to take advantage of the two kinds of algorithms. The final reconstructed image is given by ISRA (or RLA) and therefore has an ensemble of interesting properties, such as it is non-negative. Closeness with the optimum Wiener solution ensures that the reconstructed image is not too much corrupted by noise, as wished.

The Wiener filter has the advantage of being defined from physical considerations on SNR in the Fourier plane. It is based on the simple fact that spatial frequency components are differently corrupted by noises. For correctly sampled data, low and intermediate spatial frequencies are likely to contain more signal than noise, and the highest spatial frequencies components tend to contain essentially pure noise. It is clear for example, that the spectral energy detected outside the cut-off frequency of the telescope used is entirely due to noise, or to algorithmic artifacts. The Wiener technique is well known (Brault & White 1971) and we just present here the main features.

The Wiener filter operates in the Fourier plane. Let u denote the two-dimensional spatial (or angular) frequency (u₁ u₂). Let Y(u), H(u) and $X_{\rm W}(u)$ denote the Fourier transforms of the observed image, point spread function and of the reconstructed Wiener image, respectively. In optics, the quantity H(u) is the optical transfer function.

We have:

$$X_{\rm W}(u)=\frac{Y(u)}{H(u)}\Phi(u),$$ (23)

where 1/H(u) is the raw inverse filter, and $\Phi(u)$ is the Wiener zero-phase filter defined as:

$$\Phi(u)=\frac{P_{Hx}(u)}{P_{Hx}(u)+P_N(u)}\cdot$$ (24)

The quantities P_Hx(u) and P_N(u) are the power spectra of the noiseless image and of the noise, respectively. Roughly, P_N(u) can be considered as an additive constant quantity. Of course, this may be true only in average. Since H(u) is given by the autocorrelation function of the phase-aberrated aperture, P_Hx(u) is limited in angular frequencies to $D/\lambda$, D being the telescope diameter and $\lambda$ the mean light wavelength. Spatial frequencies outside the telescope cut-off frequency are removed by the Wiener filter. Below that limit, the Wiener filter preserves the high SNR frequency components, while attenuating the low SNR frequency components. The resulting Wiener image x_W(r) is recovered by an inverse Fourier transform of X_W(u).

In our original proposal (Lanteri et al. 1998), we suggested to compare the overall effect of the Wiener procedure Abs $[\Phi(u/H(u))]$ on the modulus of the spatial frequencies to the equivalent quantity Abs [X^(k)(u)/Y(u)] obtained for ISRA or RLA during the iteration, and to choose the value of k for which these modulation transfer function best agree. In practice, depending on possible near-zero values of H(u) and Y(u), these quantities tend to present unwanted peaks that make difficult the computation of a distance between them.

In the present work, we propose as a measure of goodness-of-fit between the images given by the iterative algorithms and the Wiener filter, the Euclidean distances $E_{\rm W}(k)$ between modules of spatial frequencies of the images:

$$E_{\rm W}(k)=\vert\vert{\rm Abs}[X^k(u)]-{\rm Abs}[X_{\rm W}(u)]\vert\vert^2.$$ (25)

This quantity gives better results than the simple Euclidean distance between the reconstructed images of ISRA (or RLA) and Wiener, because of possible negative values in the Wiener result. Indeed, the drawback of the Wiener technique is that the non-negativity of $x_{\rm W}(r)$ is not assured. Moreover, the Wiener technique is developed in the Fourier plane and it seems reasonable to make the comparison there.

A numerical simulation is made using the data material presented in Sect. 2; for the sake of conciseness, we limit the presentation to the case of Poisson noise. The Wiener filter is constructed taking for P_Hx(u) the squared modulus of the Fourier transform of the noiseless blurred image. The noise power spectrum P_N(u) is determined as the average noise value outside the frequency cut-off of the telescope.

Wiener images corresponding to the deconvolution of Fig. 2 are shown in Fig. 9. For the two reconstructed images in the Wiener procedure, about 45% of the reconstructed pixels are negative. The total negative part corresponds to about 7% of the positive parts for the 10⁴ photons image, and 4% for the 10⁶ photons image.

The modulus of the Fourier transform of these images is used to compute the distance $E_{\rm W}(k)$ defined by Eq. (25). It is represented in Figs. 10 for ISRA and RLA, and compared to the Euclidean distance $E_{\rm O}(k)$ between x^k and the object. For the 10⁴ photons image, the minimum values of $E_{\rm W}(k)$ and $E_{\rm O}(k)$ perfectly coincides for both algorithms (iteration # 7 ISRA, iteration # 6 for RLA). For the 10⁶ photons image, the minimum values of $E_{\rm W}(k)$ is obtained earlier than expected (iteration # 92 instead of # 164 for ISRA, iteration # 81 instead of # 134 for RLA). In this later case the use of the minimum value of $E_{\rm W}(k)$ tends to stop the iteration prematurely, but still give a valuable order for stopping the iterations.

Figure 10: Comparative representation of the distances $E_{\rm W}(k)$, defined by Eq. (25), and E0(k)as function of the iteration number "k''. Left frames are for the 104 photons image, right ones for the 106 photons image. Top frames ISRA, bottom ones: RLA. For all the representations, $E_{\rm W}(k)$ is lower than E0(k). For both algorithms, the comparison with Wiener gives the correct iteration number for the 104 photons image. For the 106 photons image, the Wiener filter leads to a value of k smaller than the optimum, then to slightly oversmoothed images

Figure 11: Comparison of frequency "equivalent filters" of ISRA and RLA with the Wiener filter. The circular average of ${\rm Abs}\lfloor X^{(k)}(u)/Y(u)\rfloor$ obtained in the iterative deconvolution is drawn as a function of the spatial frequency "u", for the 3 values of k of Figs. 4 and 5. They are compared to the circular average of ${\rm Abs}\lfloor X_{\rm W}(u)/Y(u)\rfloor$ given by the Wiener filter (individual black squares). For the iterative algorithms, the data material of Figs. 11A and 11B is that of Figs. 4 and 5, respectively. For the 104 photoelectrons image (Fig. 11A), the mean frequency gain is of the order of two; it is of the order of six for the 106 photoelectrons image (Fig. 11B). There is a strong similarity between ISRA and RLA results. In the very low range of frequencies, all curves tend to fit the inverse filter. In the intermediate range of frequencies, there is a rather good agreement between the best results of the iterative algorithms and the Wiener filter. In the highest frequencies, the iterative algorithms are very soon much higher than the Wiener filter. A good agreement between ISRA or RLA and Wiener in all range of frequencies wouldrequire an explicit regularization

In Figs. 11A and 11B, we give a comparison of the overall effect produced on the spatial frequencies of the images. For each spatial frequency "u'', the ratios ${\rm Abs}\lfloor X^{(k)}(u)/Y(u)\rfloor$ for ISRA and RLA are compared with the equivalent quantity given by the Wiener filter ${\rm Abs}[X_{\rm W}(u)/Y(u)]$. To obtain smoother results, we have represented, for each modulus of frequency | u|, the ratio of the circular average of Abs [X^(k)(u)] and Abs $[X_{\rm W}(u)]$ to the equivalent circular average of Abs[Y(u)]. The comparison is limited to the results obtained in the Poisson statistics.

For a given image, the equivalent filters obtained for ISRA and RLA are very similar. In the very low range of frequencies, all curves tend to fit the inverse filter. In the intermediate range of frequencies, there is a rather good agreement between the best results of the iterative algorithms and the Wiener filter. In the highest frequencies, the iterative algorithms are very soon much higher than the Wiener filter. A good agreement between ISRA or RLA and Wiener in all range of frequencies would require an explicit regularization.