S. Correia^1,2 - M. Carbillet³ - P. Boccacci⁴ - M. Bertero⁴ - L. Fini³

1 - UMR 6525 Astrophysique, Université de Nice-Sophia Antipolis, Parc Valrose, 06108 Nice Cedex 2, France
2 - European Southern Observatory, Alonso de Cordova 3107, casilla 19001, Santiago 19, Chile
3 - Osservatorio Astrofisico di Arcetri, largo E. Fermi 5, 50125 Firenze, Italy
4 - INFM and DISI, Università di Genova, via Dodecaneso 33, 16146 Genova, Italy

Abstract
The Software Package AIRY (Astronomical Image Restoration in interferometrY) described in this paper is designed for simulation and/or multiple deconvolution of interferometric images. It was conceived for applications to the Large Binocular Telescope (LBT), but can also be used for future interferometers of the same class. AIRY is modular, IDL-based, and designed to be used together with the CAOS (Code for Adaptive Optics Systems) Application Builder. As concerns multiple deconvolution AIRY includes an implementation of the method OS-EM (Ordered Subsets-Expectation Maximization) which is an accelerated version of the Lucy-Richardson method. After a description of the structure and modules of the software package, its use is illustrated with a few applications. In particular the intrinsic performance of the implemented restoration method is explored in the case of binary stars of different angular separations and magnitude differences, as well as in the case of a diffuse object. The possibility of obtaining super-resolved images of unresolved binary stars is also demonstrated.

Key words: techniques: interferometric - techniques: image processing - methods: data analysis -
methods: numerical - stars: imaging

1 Introduction

The Large Binocular Telescope (LBT) is the first example of interferometers of a new conception, designed to provide high-resolution imaging of a wide field in the near-infrared and visible wavelengths domain. The basic feature of LBT and of future LBT-like interferometers (as for example the 20/20 telescope project) is the possibility of getting a good coverage of the u-v spatial frequency plane by means of a few observations at different parallactic angles. Indeed LBT will consist of two 8.4 m pupils placed 14.4 m apart on a common alt-azimuthal mount. It is currently under construction on the top of Mount Graham in Arizona and first light is planned for early 2004. As a permanent part of LBT, the foreseen high-level Adaptive Optics (AO) system, including adaptive secondary mirrors, pyramid wavefront sensors, and multi-conjugate AO techniques, is intended to achieve a high-quality correction (up to $\sim$ 90% Strehl ratio in the near-infrared) on a large part of the field (typically 2 arcmin). However, in order to be able to produce high-resolution, possibly deep and wide-field, imaging from interferometric data, appropriate restoration methods must be used.

The basic feature of the restoration problem for an interferometer like LBT is that a unique high-resolution image must be extracted from different interferometric observations of the same target, obtained with different orientations of the baseline with respect to the sky (i.e. different parallactic angles). Such a problem is essentially a problem of multiple image deconvolution (Piana & Bertero 1996), and is an extension of the usual image restoration problem where only one observation of the target is available. Well-know methods for the solution of the latter can be generalized to the case of multiple images. For instance Correia & Richichi (2000) provide the first application of the extension of the Lucy-Richardson (LR) method to LBT-like data while validation and comparison of the extensions of some iterative and non-iterative restoration methods is discussed in Bertero & Boccacci (2000b).

The main difficulty of the LR method, at least from the computational point of view, is the slow convergence of the iterations. However the analogy between the LBT restoration problem and computerized tomography suggests to investigate the applicability to LBT of methods introduced in that field for the acceleration of the Expectation Maximization (EM) method, which is the version of LR used in tomography (Shepp & Vardi 1982). As a result of this analysis a version of the method Ordered Subsets - Expectation Maximization (OS-EM), suitable to the LBT problem, is introduced in Bertero & Boccacci (2000a). The new method provides, when compared with LR, a reduction of the number of iterations by a factor p (where p is the number of observations), with approximately the same computational cost per iteration. Since the processing time is approximately independent of the number of observations, this method is well suited for a quick look in the restoration of images from interferometric data as it will be required by LBT. A practical implementation of this technique is included in the Software Package AIRY presented in this paper. However it is important to point out that, thanks to the particular modular structure of AIRY, it is easy to add to the package other modules implementing other restoration methods as those considered in Bertero & Boccacci (2000b), or blind deconvolution methods as those presented by Hege et al. (1999).

The acronym AIRY (Astronomical Image Restoration in interferometrY) was first introduced for describing the activity of a group of astronomers and mathematicians collaborating on the restoration of interferometric images (see the related WEB page: http://dirac.disi.unige.it). The first result of this collaboration is just the Software Package AIRY, which is based on IDL and designed to be used together with the CAOS (Code for Adaptive Optics Systems) Application Builder (Fini et al. 2001). It can be used for simulating optical and near-infrared interferometric observations, and to perform subsequent image restoration, with application to LBT or LBT-like interferometers.

The purpose of this paper is not only a description of the Software Package AIRY but also a demonstration of its use in the simulation and restoration of LBT images. Therefore the paper is organized as follows. In Sect. 2 we recall the restoration method we have implemented and we discuss the structure of AIRY. In Sect. 3 we give a detailed description of the modules of the present version of AIRY and of their tasks. In Sect. 4 we demonstrate the use of AIRY by a few applications to LBT. The topics considered in this last part concern the limits of the restoration method in evaluating astrometric and photometric parameters of binary stars with different values of angular separation and magnitude difference, the possibility of using OS-EM as a super-resolving method in the case of unresolved binary stars, and the imaging of a diffuse and extended object. In all these cases we assume perfectly co-phased and corrected point-spread functions (PSFs) as well as a declination of the object such that LBT can assure a complete coverage of the u-v plane. The effects on imaging due to partial angular coverage and to partial AO correction are investigated in the subsequent paper (Carbillet et al. 2002). We finally summarize the paper and give some conclusions and perspectives in Sect. 5.

2 Structure of the software package AIRY

We first give a brief description of the OS-EM method which is the restoration method implemented in the present version of AIRY. As already mentioned in the Introduction it is basically an accelerated version of the LR method, when applied to multiple images deconvolution, and therefore provides the same results as LR when the iterations are suitably stopped.

We denote by $\vec{f}$ the array $N \times N$ which represents the astronomical target to be recovered and by $\vec{g}_j$ , j=1, 2, ..., p, the arrays $N \times N$ corresponding to p observations of the target, performed with p different orientations of the baseline of the interferometer. The PSFs related to the observations $\vec{g}_j$ are denoted as $\vec{K}_j$ and, in this paper, are supposed to be space-invariant. Therefore we assume the following model for image formation:

where the symbol * denotes convolution, the quantities $\vec{b}_j$ are the average backgrounds due to sky emission and the quantities $\vec{w}_j$ are the noise terms. We point out that, even if we write the noise as an additive contribution to the images, this does not imply that we are only considering signal-independent noise. Indeed, in our simulations we introduce both read-out and Poisson noise.

This model of image formation is that used in the generation of LBT images by AIRY. As concerns the image restoration problem, it consists in estimating the astronomical target $\vec{f}$ from the knowledge of the p observations $\vec{g}_j$ , of the p corresponding PSFs $\vec{K}_j$ and of the p average backgrounds $\vec{b}_j$ . Therefore it is the problem of deconvolving multiple images of the same object in the presence of background and noise. As already mentioned in the Introduction, a few methods for the solution of this problem are proposed and compared in Bertero & Boccacci (2000b). As a result of this comparison, we have found convenient to implement only one method in the first version of AIRY, i.e. OS-EM.

If we denote by $\vec{K}_j^{(-)}$ the PSF obtained by reflecting $\vec{K}_j$ with respect to the center of the array then the proposed version of the OS-EM iterative method is as follows:

where, as usual, products and quotients of images have to be intended as products and quotients pixel by pixel and the index k characterizes the iterations. Each iteration consists of a cycle over the p observations. We also point out that this version of the method applies to the case where each PSF is normalized in such a way that the sum of all its values is one.

After these preliminary remarks, we describe now the structure of the Software Package AIRY which is modular and designed to simulate optical and near-infrared interferometric observations and/or to perform subsequent image restoration/deconvolution, using the model of image formation and the method of image deconvolution described above.

The Software Package AIRY can be applied to a wide range of imaging problems, including the post-processing of real data. It is written in IDL, and consists of a set of specific modules, each representing a particular task, which are described in the next section. AIRY is conceived to be used together with the CAOS Application Builder, from version 2.0 onwards (current version is 3.0). CAOS has been designed to provide a graphical programming environment where elementary building blocks can be assembled together to create complex simulation applications in a straightforward manner, so that the user can concentrate on the scientific aspects of her/his problem, while mundane coding problems are managed by some automatic tool. It basically allows a user to build a simulation program (a project) by putting together the needed elementary building blocks (the modules) and specifying the data flow between each block. When the project has been defined to the user's satisfaction the IDL code which implements the program is automatically generated. See Fini et al. (2001) for a more detailed description of the CAOS Application Builder.

As a CAOS-compatible software package, AIRY permits to build different applications by linking the modules together and setting their parameters through their own graphical user interfaces (GUIs). The concept of the Software Package AIRY is identical to that of the already existing Software Package CAOS (the original software package associated to the CAOS Application Builder - see Carbillet et al. 2001 - and devoted to adaptive optics systems simulations). In particular new modules can be easily implemented by the user for her/his particular needs. The Software Package AIRY also includes a library of ideal and AO-corrected LBT PSFs (produced using the Software Package CAOS), a library of useful routines (in addition to those related to each module), a hyper-text help which can be called from the individual GUI of each module, and a set of examples of typical simulation projects.

Table 1 shows a complete list, together with a very brief description, of the modules of the current version 1.0 of the Software Package AIRY. A more detailed description will be given in the next section. Nevertheless it may be important to remark that any module can have a number of inputs/outputs varying from zero to two and that only two types of input/output are defined within the present version of AIRY: the object type and the image type. Since AIRY is an image-oriented package, almost all modules have inputs and outputs of the type image, except the module OBJ (no input, and only one output of type object), and the module CNV (one input of type object, and one of type image, and only one output of type image). Finally the utility module DSP can display both types of input.

However, in order to clarify the role of each module, it is convenient to refer to a few typical applications. Two examples of projects are represented in Fig. 1. Both are essentially composed of three parts: a data modeling part, a data processing part, and a data analysis part.

Table 1: Descriptive list of the modules of the Software Package AIRY (version 1.0).
Module Purpose

Data simulation modules

OBJ - OBJect definition -to define the object characteristics

CNV - object-PSF CoNVolution -to perform convolution

ADN - ADd Noise to image -to add read-out and Poisson noise

Data processing modules

PRE - PRE-processing -to perform image pre-processing

DEC - DEConvolution process -to perform deconvolution (OS-EM method)

Data analysis modules

ANB - ANalysis of Binary -to analyse reconstructed images of binary objects

FSM - Find Star Module -to detect stars in the reconstructed images

Other modules and utilities

RFT - Read FiTs file format -to read FITS images

WFT - Write FiTs file format -to write FITS images

RSC - Restore im. Struct. Cubes -to restore image structure cubes (XDR or FITS format)

SIM - Save IMage struct. -to save image structure cubes (XDR or FITS format)

DSP - data DiSPlay -to display images

**Table 1:** Descriptive list of the modules of the Software Package AIRY (version 1.0).
Module	Purpose
Data simulation modules
`OBJ - OBJect definition`	-to define the object characteristics
`CNV - object-PSF CoNVolution`	-to perform convolution
`ADN - ADd Noise to image`	-to add read-out and Poisson noise
Data processing modules
`PRE - PRE-processing`	-to perform image pre-processing
`DEC - DEConvolution process`	-to perform deconvolution (OS-EM method)
Data analysis modules
`ANB - ANalysis of Binary`	-to analyse reconstructed images of binary objects
`FSM - Find Star Module`	-to detect stars in the reconstructed images
Other modules and utilities
`RFT - Read FiTs file format`	-to read FITS images
`WFT - Write FiTs file format`	-to write FITS images
`RSC - Restore im. Struct. Cubes`	-to restore image structure cubes (XDR or FITS format)
`SIM - Save IMage struct.`	-to save image structure cubes (XDR or FITS format)
`DSP - data DiSPlay`	-to display images

$\begin{figure} \par\includegraphics[width=17.2cm,clip]{H3418F1.eps}\end{figure}$	Figure 1: Two examples of projects using the Software Package AIRY. Left: restoration and data analysis in the case of real observations of binary stars. Right: simulation of star cluster observations, followed by image restoration and data analysis.
Open with DEXTER

The first project (left side of Fig. 1) is an example of treatment of real data in the case of the observation of a binary star. The data modeling part of this application is very simple because it is reduced to read the files containing the observed data and the corresponding PSFs obtained, for instance, from observations of a reference star. Therefore data modeling consists of two occurrences of the module RFT (one for the interferometric images and one for the PSFs), followed by a data processing part which in turns consists of a pre-processing stage (module PRE), and the subsequent multiple deconvolution (module DEC). Last part concerns the analysis of the deconvolved images, and in particular the detection and characterization of the retrieved binary star parameters with the help of the module ANB.

The second project (right side of Fig. 1) is the simulation of the observation, image restoration, and data analysis of a star cluster. Here the data modeling part is clearly more advanced because the object is first modeled within module OBJ, and then convolved with PSFs taken, for example, from the PSF library present in the package. To this purpose the modules RFT and CNV are used. Both the object data and reference data are then treated within module ADN in order to add the different kinds of noise contributions. The data processing part is then similar to that of the previous example (pre-processing of the data followed by multiple deconvolution), while the data analysis part consists of the analysis of the reconstructed star cluster by the module FSM.

In each application, utility modules allow to display and save the data at each stage of the data flow and at each iteration of the process.

As clarified by the previous examples, we use the following classification of the modules of AIRY: data modeling modules, data processing modules, data analysis modules, and utility modules. The next section describes in detail the tasks of the modules of these different classes.

It must also be stressed that, since the Software Package AIRY is fully modular, the possible applications are not limited to the previous two cases. The output of each module can be linked to the (compatible) input of another one, and ad hoc modules proper to the needs of the user can be easily added to AIRY, thanks to the template module given along with the package. Moreover, and thanks again to its modular structure, we believe the Software Package AIRY to be an ideal tool for fair comparison of different methods for image reconstruction, data pre-processing, etc.

3 Description of the AIRY modules

In this section we describe in detail the tasks of the modules of each one of the four classes indicated above.

3.1 Data modeling modules

This class contains the modules OBJ (OBJect definition), CNV (object-psf CoNVolution) and ADN (ADd Noise to image).

The module OBJ is designed to define the characteristics of an object chosen among several pre-defined types (binary object, open cluster, planetary nebula, Super Novae remnant, spiral galaxy, Young Stellar Object, stellar surface), or modeled by the user. For each object it is possible to adjust the corresponding relevant parameters. For instance planetary nebulae or Super Novae remnants modeling allows to enter emission-lines fluxes, which are subsequently converted into the flux within the observed bandwidth. Another example is provided by stellar surface modeling which allows to create maps of uniform elliptical disks with possibly the presence and characterization of a bright spot feature. Figure 2 shows the GUI corresponding to this module, with the main interface in the upper part, and the auxiliary interface corresponding to the particular case of a Young Stellar Object (YSO) in the bottom part. The whole set of available parameters is visible in this particular example, which is just the example of diffuse object considered in Sect. 4.

The module CNV performs an FFT-based convolution of the two inputs, typically the object map and a PSF cube (where usually each element of the cube represents the PSF at a given parallactic angle). The module CNV sets also the telescope- and detector-dependent parameters such as the overall efficiency, the integration time, and the telescope collecting surface. Moreover the module CNV rebins the object map with the pixel size of the PSF and puts each frame of these two inputs into an array with a dimension equal to the closest power of 2.

Finally the module ADN adds a defined sky background contribution to the outputs of the module CNV and provides a Poisson noise realization of the results. It can also add a computed contribution of the detector Read-Out Noise (RON).

3.2 Data processing modules

The module PRE performs all the operations to be done on the data before applying the deconvolution process. For instance it evaluates the value of the sky background, a fundamental parameter which must be accurately measured in order to fully exploit the positivity constraint assured by the implemented deconvolution method. This operation is performed using the so-called MMM algorithm taken from the astrolib library. The module makes also possible to enter the value of the background if this value has been estimated outside AIRY. Moreover the module PRE computes the peak Signal-to-Noise Ratio (SNR), defined as the ratio between the peak value and the rms of the noise contribution, and allows to re-center the frames of a data cube in order to correct for misalignment of the object from frame to frame. This operation is performed by means of an FFT-shift technique. The module PRE also includes the possibility of performing a border apodization which can be used in order to minimize "ripples" effects in the restored image. The filter used in this first version is a 2-D square function convolved with a Gaussian one of given FWHM, in order to smooth the borders.

The module DEC implements the OS-EM multiple deconvolution method, shortly described in Sect. 2. Its two inputs are the cube of the observations and that of the corresponding PSFs (in general different from those used in the module CNV). Both inputs include not only the observed images and PSFs but also an evaluation of the background contributions to them, obtained by means of the module PRE. As concerns the number of iterations its optimal value depends on the object and noise. Automatic criteria for stopping the iterations are under investigation. In the present version (1.0) of the Software Package AIRY the number of iterations has to be chosen by the user as a global parameter of the simulation project.

3.3 Data analysis modules

The module ANB performs the detection of a binary object, computes its photometric and astrometric parameters, and possibly compares them with entered binary model parameters. The detection of each component is achieved by an ad hoc routine, which detects in the input frame the peak location over a user-supplied threshold (expressed in rms background noise unit) and masks the pixels inside a disk of user-supplied detection radius, centered on the peak, before continuing. Centroid calculations are carried out with box size of user-supplied dimensions centered on the peak pixel of each detected star, yielding separation and position angle (PA) of the binary. Difference of magnitude is measured by using aperture photometry with a set of apertures of given diameters centered on the measured centroid location. The flux zero-point can be supplied by the user (when dealing with real data for instance) or can be computed.

The module FSM provides the detection of all the stars contained in an image, computing their location and aperture photometry. The module uses the astrolib routine findstar. Input parameters are the detection threshold (in rms background), the flux zero-points and the set of aperture photometry diameters.

3.4 Utility modules

4 Examples of application

The purpose of this section is to simulate high-resolution interferometric observations of astronomical objects of interest with LBT, and to evaluate the scientific parameters of these objects as retrieved by the deconvolution method implemented in AIRY. From this point of view the present work is to be considered as the continuation of preliminary studies already presented elsewhere (Correia et al. 2000; Carbillet et al. 2000; Correia et al. 2001).

As a relevant first example of application we treat the fundamental case of binary stars with different angular separations and magnitude differences, looking also at the possibility of achieving super-resolution with OS-EM. The second example of application concerns a diffuse stellar object. Indeed these two classes of objects have a different behaviour as concerns deconvolution because they correspond to two different spatial-frequency bands to be retrieved by the implemented method.

We restrict ourselves to the case of perfectly co-phased and AO-corrected PSFs, as well as to a complete coverage of the u-v plane. Therefore we first describe the set of the PSFs we are using in our simulations, assuming that they correspond to the observation of an object with a suitable declination.

4.1 The interferometric PSFs and other parameters

We consider a set of LBT interferometric PSFs corresponding to the observation of an object with a declination of +80 $\degr$ , assuming that observation is limited to airmass <2. More precisely we assume 6 observations equispaced in time in the allowed range of hour angles, that is between -5 h and +5 h. This corresponds to a range of parallactic angle of $\sim$ 200 $\degr$ , thus providing a complete u-v plane coverage. Indeed, due to the alt-az mount of the telescope, the projection of the baseline onto the observed object in the sky plane rotates and this position angle is what one calls the "parallactic angle'', which is a function of time (hour angle) and coordinates of the object on the sky (declination).

We assume an integration time per baseline of about 20 min. During this integration time and for the declination we are considering, the average rotation of the baseline is of about 5 $\degr$ in parallactic angle. This effect is taken into account in the computation of the PSFs by integrating over the interval of parallactic angles indicated above. The filter used is the K broad-band filter and we correct for the spectral width of the filter ( $\Delta\lambda =$ 400 nm) by integrating also over $\lambda$ .

The spatial sampling adopted for all the interferometric PSFs used in this paper (except when we mention explicitly a different one) is 3 pixels per Full-Width at Half-Maximum (FWHM) of the fringes, where the FWHM is equal to $\lambda /L$ , with L the total baseline length (22.65 m since the effective diameter of each mirror is 8.25 m due to the adaptive secondary mirror of LBT) and $\lambda$ the filter central wavelength (here 2.2 $\mu$ m).

Other parameters used in our simulations are: a total efficiency (mirror+optics+detector) of 25%, a sky background brightness of 12.5 mag/arcsec² and a read-out noise corresponding to the scientific imaging device of 2 ${\rm e}^-$ rms.

4.2 Intrinsic limits in the restoration of the parameters of binary stars

We investigate the intrinsic limitations of our method in the reconstruction of binary stars with different values of angular separation $\rho$ and difference in magnitude $\Delta m$ . The observed images are simulated by means of the modules CNV and ADN.

We consider five binary stars as observed with LBT and presenting increasing values of separation and magnitude difference. The five cases we consider here are the following:

The first example is likely to be considered as a limiting case for astrometry (if super-resolution is not attempted - see Sect. 4.3 ), and the last three examples for photometry. We place ourselves in a case of high SNR (with the assumed 20 min integration time per baseline and a main component magnitude K=10).

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H3418F3a.ps}\par\includegr... ...lip]{H3418F3b.ps}\par\includegraphics[width=8.8cm,clip]{H3418F3c.ps}\end{figure}$	Figure 3: Restoration of binary parameters in the case of ideal PSFs and complete coverage of the u-v plane: error plots of the angular separation $\rho$ , position angle PA, and magnitude difference $\Delta m$ , for the five cases defined in the text, versus the number of iterations of the OS-EM method. In the legends "res'' stands for $\lambda /L$ .
Open with DEXTER

The results are shown in Fig. 3, where the errors in the retrieval of the three binary parameters (angular separation $\rho$ , position angle PA and magnitude difference $\Delta m$ ) are plotted as functions of the number of iterations of OS-EM. We comment on these results.

As expected it is much more difficult to retrieve the parameters whenever the companion is closer to and/or much fainter than the main component. More precisely, in the first two cases (same magnitude difference but different angular separation), the accuracy is higher when the angular separation is larger even if in both cases results are quite good and the correct values are reached after a rather small number of iterations (very small in the second case, but we recall that we are considering a case with a large SNR value).

As concerns the second and third case (same angular separation but different magnitude difference), the increase in $\Delta m$ implies that a higher number of iterations is required for distinguishing the weak companion from the remaining sparse energy of the main component. However, after detection, the retrieved binary parameters converge to rather accurate values. The effect increases for increasing magnitude difference. Indeed the detection of the right companion takes place after $\sim$ 50 iterations for the ( $\rho =$ 5 $\lambda /L$ , $\Delta m =9$ ) binary and after $\sim$ 200 iterations for the ( $\rho =$ 15 $\lambda /L$ , $\Delta m = 13$ ) binary; then the parameters converge to the right values, with the exception of the angular separation of the ( $\rho =$ 15 $\lambda /L$ , $\Delta m = 13$ ) case.

$\begin{figure} \par\includegraphics[width=5.9cm,clip]{H3418F4a.ps}\hspace{1mm} ... ...F4b.ps}\hspace{1mm} \includegraphics[width=5.9cm,clip]{H3418F4c.ps}\end{figure}$	Figure 4: Illustration of the possibility of super-resolution. From left to right: the original object; the non-resolved restored image obtained with 1000 iterations of OS-EM; the resolved restored image obtained with the use of the mask defining the object domain as deduced from the first restoration.
Open with DEXTER

In conclusion, for all binaries considered here except one (the last case), the parameters converge to the correct values, although with different convergence rates.

4.3 Super-resolution

In astronomy as well as in microscopy, super-resolution is a word used for denoting a method, in general a processing of the data, which allows to reach a resolution beyond the diffraction limit. Super-resolution requires an extrapolation of the object spectrum outside the band of the instrument and is possible in a limited number of situations. In problems of image deconvolution it can be obtained by the use of specific constraints, satisfied by the object to be restored, which must be implemented in the restoration method (for a tutorial see Bertero & De Mol 1996). In the case of astronomical objects these constraints are essentially positivity and the extent of the object, which must not be much greater than the angular resolution of the telescope. The LR method, and hence the OS-EM method, can be nicely used for such a task. Indeed positivity is automatically assured. As concerns the information about the angular extent of the object it can be easily inserted in the method by an appropriate initialization. Indeed if the result of one iteration is zero in a certain pixel, then it will be zero at that pixel for all subsequent iterations.

A possible application to the case of unresolved binary stars is the following. By means of OS-EM, initialized with a constant image, one first estimates the domain of the unresolved object, by taking for instance all pixels where the values of the restored object are greater than a certain threshold. Next the deconvolution based on OS-EM is repeated by initializing the method with a mask which is one over the estimated domain and zero elsewhere.

We present here the results of a simulation carried out in order to demonstrate the possibility of super-resolution. To this purpose we consider a binary star whose angular separation is about one half the resolution limit of LBT in K-band ( $\rho = \lambda/2L$ ) while the magnitude difference is the same as in the first two cases of the previous section ( $\Delta m =1$ ). Since we are looking for super-resolution we assume that the sampling of the PSFs defined in Sect. 4.1 is increased by a factor two, so that a frequency which is twice the cut-off frequency of the telescope can be sampled. It is outside the scope of this paper to discuss how such an oversampling can be reached in practice. With a K=20 primary and $\Delta m =1$ the SNR of one baseline observation is $\sim$ 60.

To this example we apply the procedure outlined above by taking a threshold value of $10 \%$ of the peak value of the primary. The results are shown in Fig. 4. We can easily see that, as expected, the binary is not resolved by the first restoration while it is resolved by the second one when the object domain is provided as an input to the algorithm. The binary parameters of the super-resolved restoration are in a satisfactory agreement with the true ones (see Table 2), although they are obviously not so precise as in the case of a standard resolved binary.

Table 2: Measured binary parameters of the restored super-resolved image compared with the true parameters. Measurements are performed by means of the module ANB.
$\rho$ (mas) $\theta (\degr)$ $\Delta m$

model 10.5 18.4 1.0

restoration 9.2 16.9 0.3

**Table 2:** Measured binary parameters of the restored super-resolved image compared with the true parameters. Measurements are performed by means of the module ANB.
	$\rho$ (mas)	$\theta (\degr)$	$\Delta m$
model	10.5	18.4	1.0
restoration	9.2	16.9	0.3

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H3418F5.ps}\end{figure}$

Figure 5: Cuts of the power spectra along the binary axis. The dashed line corresponds to the original object, the dotted line to the result of the first restoration while the solid line to the result of the second restoration. The super-resolution effect is evident (the telescope cut-off is located at 128/6 $\simeq$ 20 pixels), but it is also evident that the magnitude difference is not correctly retrieved.

Open with DEXTER

These results are preliminary and are given to demonstrate the feasibility of super-resolution. It is obvious that much more refined simulations are required and much wider experimentation is needed in order to explore the limits of the method.

4.4 Restoration of a diffuse object

$\begin{figure} \par\includegraphics[width=5.9cm,clip]{H3418F6a.ps}\hspace{1mm} ... ...F6b.ps}\hspace{1mm} \includegraphics[width=5.9cm,clip]{H3418F6c.ps}\end{figure}$	Figure 6: From left to right: the original object; the image corresponding to a parallactic angle of $\sim$ 60 $\degr$ ; the restored image. The image extension is about $0\farcs7$ . The grey scale is square root.
Open with DEXTER

We test here our method in the case of an object with an extended and diffuse shape. To this purpose we consider an embedded young stellar object, also known as Class I object, derived from HST near-infrared observations of IRAS 04302 +2247. This is a remarkable object showing a bipolar scattered light nebulae and a totally opaque dust lane, probably featuring a large optically thick circumstellar disk seen edge-on (Padgett et al. 1999). We scale the image to have a $0\farcs85 \times 0\farcs85$ extension, so that our object is more distant by a factor $\sim$ 10 than the original one, i.e. about 1.4 kpc. With such a distance and the K=15 mag that we set here, instead of the true K=10.9 (Lucas & Roche 1997), our simulated object is intrinsically brighter by roughly 1 mag. With this brightness and the integration time considered in all our simulations (20 min), the single baseline image peak SNR is $\sim$ 180. Figure 6 shows the image of the object, as well as one of the observed images and the result of the restoration.

The accuracy of the restoration provided by the kth iteration can be measured by computing the relative rms error defined as follows:

where $\vec{f}^{(k)}$ is the result of the kth iteration (see Eq. (2)) and $\vec{f}$ is the original model (see Eq. (1)). The norm is the Euclidean one, i.e. the square root of the sum of the squares of the pixel values. In Fig. 7 we plot the restoration error R^(k) as a function of the number of iterations k. The curve has a minimum at k = 263, with a minimum value of $5.6 \%$ (the restored image of Fig. 6 corresponds to this minimum). We observe however that in this example the minimum is rather flat so that after 70 iterations we already have a restoration error of about $7 \%$ . Such a restoration is not visually distinguishable from that given in Fig. 6.

The behaviour of the restoration error shown in Fig. 7 is typical of the regularizing iterative methods and is known as semiconvergence - see, for instance, Bertero & Boccacci (1998). The first iterations improve the accuracy of the approximation of the unknown object but after a certain number of iterations the amplification of the noise increases and the quality of the restoration is degraded. In some cases the effect is much stronger than in our example; in these cases the choice of the number of iterations may be critical.

Another way for checking the quality of the restored image may consist in comparing its power spectrum with that of the object. This is shown in Fig. 8, where the profiles of the azimuthaly averaged power spectra of both images are plotted. The two curves match fairly well up to the telescope cut-off frequency located at about 40 pixels.

5 Summary and concluding remarks

In this paper we have presented the first release of the image restoration and simulation package AIRY (Astronomical Image Restoration in interferometrY). After a description of the restoration method OS-EM implemented in the package, we discuss the modular structure of AIRY and the tasks of the modules available in the present version 1.0. Next we investigate the intrinsic capabilities of the implemented restoration method in the particular case of LBT. To this purpose we consider the ideal situation of perfectly co-phased and AO-corrected PSFs and we apply the method both to the restoration of binary objects and to the restoration of diffuse objects. In the first case resolution and detection are analysed both for angularly close objects and for high brightness ratios.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H3418F7.ps}\end{figure}$	Figure 7: Restoration error R^(k) as a function of the number of iterations. The convergence occurs at iteration number 263, with $R^{(263)} = 5.6\%$ .
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H3418F8.ps}\end{figure}$	Figure 8: Profiles of the azimuthaly averaged power spectra. Solid line corresponds to the restored image (k = 263), dotted-line to the object. The telescope cut-off frequency is located at 128/3 $\simeq$ 40 pixels.
Open with DEXTER

Our simulations indicate that the fundamental limit of the OS-EM method, at least as concerns the restoration of binary stars, is essentially a limit in the estimation of the brightness ratio which can be larger in the case of wider angular separation. Indeed the difference of magnitude can be as high as $\Delta m =$ 5 for a separation of only two resolution limits and $\Delta m =$ 13 for 15 resolution limits. The upper limit in brightness ratio for more separated objects is essentially given by the detection limit of the companion star alone. This is similar to the limitation imposed by the uncorrelated speckle noise, present in the halo of the PSFs in classical AO, for the detection of faint companions (see e.g. Racine et al. 1999). We have also found that the higher is the brightness ratio between the companion and the primary, the larger is the number of iterations required to actually detect it. But, even if the convergence of the parameters is slower, the bias in their estimation tends to decrease systematically with the number of iterations. Such a conclusion obviously holds true only in the case of high SNR, assuming that it is possible to calibrate exactly the PSFs of the system (in time and space). In real life it would hardly be the case, except for very good seeing conditions and/or with the help of techniques such as multi-conjugate AO which should be able to approach this limit.

We have also demonstrated the possibility of obtaining super-resolution for this kind of objects, at least in the case of a reasonable magnitude difference. Such a result, which might open very promising future developments, can be obtained without modifying the deconvolution algorithm implemented in AIRY; it is only necessary to use a suitable initialization of the iteration process. Of course such a remarkable property holds true not only in the case of OS-EM but also in the case of LR and, more generally, of any multiplicative iterative method.

In the case of a realistic diffuse object we have shown that it is possible to obtain quite accurate restorations up to the cut-off frequency of the telescope, if the amplification of the noise is not limiting.

All the results presented in this paper have been obtained by assuming ideal PSFs and complete coverage of the spatial frequency plane. The study of the effects of partial angular coverage and partial AO correction is the subject of the subsequent paper (Carbillet et al. 2002). However we believe that we have already demonstrated that OS-EM can be conveniently used at least for a quick look of the observations with LBT.

Obviously if some specific object, with a structure more complex than that of the examples considered in this paper, requires a more refined analysis, that can be done outside AIRY with methods more sophisticated than OS-EM or also using future modules of AIRY implementing those methods.

Indeed new modules can be easily implemented by any user, depending on his/her specific need, and can be included and distributed in future versions of AIRY. From our part, we are currently studying the possibility of implementing improved versions of existing deconvolution methods, including criteria for stopping the iterations or for choosing regularization parameters, as well as blind deconvolution methods.

The Software Package AIRY, version 1.0, is obtainable under request and subscription to its dedicated mailing-list by visiting the web site http://dirac.disi.it/airy. It is delivered together with the CAOS Application Builder (version 3.0): the IDL-based graphical environment within which it was developed.

Restoration of interferometric images

I. The software package AIRY