5 Examples

We show here two reconstruction examples using the above methodology. For the first example we have worked with a $256\times 256$pixel image of the core of galaxy NGC 4321 (M100). The data were generated in the following way: an image of NGC 4321 obtained by the refurbished Hubble Space Telescope (HST) was used as source (true image), shown in Fig. 1. The raw data for the reconstruction were then obtained by convolving the true image with the Point Spread Function (PSF) of the aberrated HST, computed by the STScI TinyTim Program and adding Poisson noise. We generated two data sets: A and B using two different seeds to generate two independent Poisson noise data sets. Figure 2 show the raw image A. This example combines a diffused object and several bright sources. All the figures corresponding to this example are shown in a logarithmic scale, except where noted.

  

Figure 1: True image of NGC 4321

  

Figure 2: Raw image of NGC 4321 to be reconstructed

For the segmentation step, we first carried out a MLE reconstruction stopped after 20 iterations. The stopping point was given by the cross-validation test for the whole image. Then, we carried out the wavelet decomposition into three planes plus the residual image. Figure 3 shows wavelet plane ${\bf w}_2$ and Fig. 4 shows the residual image ${\bf a}_r$.

  

Figure 3: Second wavelet plane of the MLE reconstruction after 20 iterations

  

Figure 4: Residual image of the wavelet decomposition

Wavelet plane ${\bf w}_2$ is used to identify the stars and other prominent features. Figure 5 shows the extracted objects obtained by setting the threshold to 22.8 counts. The threshold was computed using the automatic method described in Sect. 4.3. Finally, we carried out the segmentation using the self-organizing network described in Sect. 4.4. Figure 6 shows the final segmentation into 8 extended regions plus an additional region for each bright object totaling 50 regions.

  

Figure 5: Stars and other prominent objects extracted

  

Figure 6: Final segmentation

The regional hyperparameters were computed by the cross-validation technique described in Sect. 4.5.2. Figure 7 shows the cross-validation test curve for region No. 3 (corresponding to the spiral arms of the galaxy aproximately) from which an optimum value of $\Delta a_i = 65$ was obtained. In the 50 cross-validation tests carried out for this example, the following optimum values were obtained: for the eight extended regions $\Delta a_i = 20, 30, 70, 65, 60, 150, 70, 70$respectively; for the other 42 regions (stars and bright objects) the $\Delta a_i$ was set to the maximum of 300 in 29 cases while the other 13 objects received $\Delta a_i$ values between 50 and 200.

  

Figure 7: Cross-validation test for region No. 3

Using the above set of hyperparameters, we carried out the final Bayesian reconstruction of data sets A and B with the FMAPE algorithm given by Eqs. (26), (23). After adding both reconstructions we obtained the final result, as shown in Fig. 8. A smooth background, a well reconstructed nebula, and sharp images of the stars are obtained. Noise amplification in the background and in the nebula was suppressed by the algorithm, while images of the bright objects were reconstructed to near Maximum Likelihood. Finally, Fig. 9 shows the image of the residuals in a linear scale, which exhibits no visible structure or correlation with the solution.

  

Figure 8: Final reconstruction of NGC 4321 using the FMAPE algorithm with variable hyperparameters defined by the segmentation

  

Figure 9: Map of residuals of the reconstruction using the FMAPE algorithm with variable hyperparameters (linear scale)

In order to compare the FMAPE method with variable hyperparameter with the method with a single hyperparameter given by algorithm (27) and with the Maximum Likelihood method given by algorithm (28), we carried out also the reconstructions using those methods. To make the results comparable, we reconstructed both data sets A and B and added the results to obtain the final reconstructions. The single hyperparameter for Bayesian case (27), and the stopping point for the Maximum Likelihood method (28) were also determined by cross-validation for the entire image (Núñez & Llacer 1993a). We obtained a value of $\Delta a = 100$ for the Bayesian case and a stopping point of 30 iterations for the Maximum Likelihood method. Figure 10 and Fig. 11 show the results for the Bayesian with single hyperparameter and Maximum Likelihood cases respectively.

  

Figure 10: Reconstruction of NGC 4321 using the FMAPE algorithm with simple hyperparameter

  

Figure 11: Reconstruction of NGC 4321 using the MLE algorithm

If we compare the results of the three methods, we can observe than the reconstruction with a single hyperparameter (Fig. 10 left) is noisier both in the nebula and in the background than in the variable hyperparameter case. Also, the images of the stars are under-reconstructed. The MLE image (Fig. 10 right) is clearly under-reconstructed in the bright objects. If we increase the number of iterations to correct that effect, the noise in the lower intensity regions increases excessively, making the result unacceptable. The Bayesian method with variable hyperparameter appears to be the better method of reconstruction.

The second example is the reconstruction of a $512\times 512$ pixel real image of the planet Saturn obtained with the WF/PC camera of the Hubble Space Telescope, before the servicing mission. All images for this example are shown in a linear grey scale, except where noted. Figure 12 shows the raw data to be reconstructed. The Point Spread Function is shown in a logarithmic scale in Fig. 13.

Since the image of Saturn has no stars we segmented the image into 9 extended regions. For the segmentation step, we first carried out a MLE reconstruction stopped after 50 iterations. The stopping point was given by the feasibility test for the whole image. Then, as in the first example, we decomposed the image into three planes plus the residual image and carried out the segmentation by the self-organizing network described in Sect. 4.4. Figure 14 shows 9 regions resulting from the segmentation. In this example we do not have a second data set for cross-validation. In addition, the image was obtained from a CCD camera with readout noise, rendering the process of splitting the image into two halves by thinning questionable (Núñez & Llacer 1993a). Thus, in this example we have used the feasibility test described in Sect. 4.5.1 to determine the set of hyperparameters.

Figure 15 shows the result of the 9-region reconstruction. We have obtained a smooth background, a well reconstructed image of the planet and sharp divisions in the rings. By the use of the variable resolution FMAPE algorithm, we have avoided noise amplification in all regions of the image.

  

Figure 12: Raw image for the planet Saturn

  

Figure 13: Observed PSF used for the reconstruction

  

Figure 14: Segmentation of the image in 9 regions of Bayesian hyperparameters using the method described in the text

  

Figure 15: Reconstruction of the image of Saturn using the FMAPE algorithm with the 9 channels defined by the segmentation

Again, we compare the results of our new methodology with the constant $\Delta a$ approach and with the MLE method. For that purpose, we used the Saturn data for reconstructions by the FMAPE algorithm with constant values of $\Delta a = 200$ and 500, and a reconstruction by the Maximum Likelihood method at 100 iterations. The results show characteristics similar to the ones described in the previous example: the solution with variable hyperparameter is better reconstructed and noise amplification is well controlled. In order to characterize the solutions, we computed the mean normalized residuals for each of the 9 regions of the segmentation. Figure 16 shows the plot of the mean residual as a function of region number. An ideal reconstruction would give a straight line at a mean value of 1.0. The reconstruction obtained by MLE has mean residuals that are increasing from below to above the 1.0 value. The reconstruction with a single uniform $\Delta a = 200$ has mean residuals that are too high. With $\Delta a = 500$ most of the values are too low but some are still too high. However, the reconstruction with the variable $\Delta a_i$, although not perfect, is much improved, with most of the mean residuals close to the ideal value of 1.0.

  

Figure 16: Mean residuals for the 9 regions obtained with the FMAPE algorithm with variable hyperparameter in comparison with other methods. An ideal reconstruction would give a straight line at normalized mean residual 1.0