# Application of iterative blind deconvolution to the reconstruction of LBT LINC-NIRVANA images

G. Desiderá, B. Anconelli, M. Bertero, P. Boccacci, M. Carbillet

► To cite this version:

G. Desiderá, B. Anconelli, M. Bertero, P. Boccacci, M. Carbillet. Application of iterative blind deconvolution to the reconstruction of LBT LINC-NIRVANA images. Astronomy and Astrophysics - A&A, 2006, 452 (2), pp.727-734. 10.1051/0004-6361:20054481. hal-03473279

HAL Id: hal-03473279

https://hal.science/hal-03473279

Submitted on 13 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Application of iterative blind deconvolution to the reconstruction of LBT LINC-NIRVANA images

G. Desiderá<fnref n="1" />, B. Anconelli<fnref n="1" />, M. Bertero<fnref n="1" />, P. Boccacci<fnref n="1" />, and M. Carbillet<fnref n="2" />

Received 7 November 2005 / Accepted 11 February 2006

## ABSTRACT

**Context.** The paper is about methods for multiple image deconvolution and their application to the reconstruction of the images acquired by the Fizeau interferometer, denoted LINC-NIRVANA, under development for the Large Binocular Telescope (LBT). The multiple images of the same target are obtained with different orientations of the baseline.

**Aims.** To propose and develop a blind method for dealing with cases where no knowledge or very poor knowledge of the point spread functions (PSF) is available.

**Methods.** The approach is an iterative one where object and PSFs are alternately updated using deconvolution methods related to the standard Richardson-Lucy method. It is basically an extension, to the multiple image case, of iterative blind deconvolution methods proposed in the case of a single image.

**Results.** The method is applied to simulated LBT LINC-NIRVANA images and its limitations are investigated. The algorithm has been implemented in the module BLI of the software package AIRY (Astronomical Image Reconstruction in interferometry), available under request.

The preliminary results we have obtained are promising but an extensive simulation program is still necessary for a full understanding of the applicability of the method in the practice of the reconstruction of LINC-NIRVANA images.

**Key words.** methods: numerical – techniques: image processing

## 1. Introduction

In previous papers (Bertero & Boccacci 2000a,b; Correia et al. 2002; Carbillet et al. 2002; Anconelli et al. 2005a) we developed methods and software for the deconvolution of multiple interferometric images of the same astronomical target. Moreover, our group has produced the software package AIRY (see http://dirac.disi.unige.it), which can be used in conjunction with the software package CAOS (Code for Adaptive Optics Systems; Carbillet et al. 2005), within the common CAOS "problem solving environment" (see http://www.arcetri.astro.it/caos or http://www-astro.unice.fr/caos). This tool can be applied to Fizeau interferometers such as LINC-NIRVANA (Lbt INterferometric Camera and Near-InfraRed/Visible Adaptive INterferometer for Astronomy), the German-Italian beam combiner for LBT (in the following LINC-NIRVANA will be denoted by LN). We recall that LBT consists of two 8.4 m mirrors on a common mount, with a spacing of 14.4 m between their centres, so that a maximum baseline of 22.8 m will be available. The first primary mirror has been installed in March 2004, and "first light" has been achieved on October 12, 2005. The second primary mirror has been aluminized on board in January 2006, so that LBT is now a true binocular telescope with two coated primary mirrors.

The methods we have implemented and validated are based on the extension of the well-known Richardson-Lucy method (RLM) to the problem of multiple image deconvolution. Indeed, LN will provide different images of the same astronomical

target, corresponding to different orientations of the baseline. In order to get a unique high-resolution reconstruction of the target, one needs the knowledge of the interferometric PSFs corresponding to the different images. These PSFs cannot be modeled and computed since they depend on the configuration of the atmosphere during the exposures; therefore they must be obtained from the images of reference stars. These images are corrupted by several factors such as noise and backgrounds (we assume calibrated images, already corrected for effects such as bad pixels, flat field, etc.); moreover, as a consequence of non-uniform adaptive optics (AO) correction, they may not correspond to the PSFs in the region of the scientific object. In other words, one is faced with the problem of improving the quality of the measured PSFs. Sometimes the knowledge may be so poor that the PSFs, in practice, are not known.

In the case of single-image deconvolution, methods of *blind deconvolution* (BD) have been developed for extracting both the object and the PSF from the recorded image. These methods can be used also when a poor knowledge of the PSF is available, as discussed above. In such a case the term *myopic deconvolution* is sometimes used.

BD methods can be divided into two classes. In the first class BD is attempted by minimizing an error function designed to optimize both the reconstructed object and the PSF, with the addition of suitable constraints (non-negativity, bounded domain in the angular or *u*, *v* variables, etc.). The underlying object and the PSF are obtained simultaneously. An example of these methods is that proposed by Jefferies & Christou (1993). The second class consists of methods that restore the object and the PSF separately

<footnote n="1">¹ DISI, Università di Genova, via Dodecaneso 35, 16146 Genova, Italy
e-mail: bertero@disi.unige.it</footnote>

<footnote n="2">² Laboratoire Universitaire d'Astrophysique de Nice, UMR 6525, Parc Valrose, 06108 Nice Cedex 02, France</footnote>
in an iterative form: within each iteration, which we will call a cycle, either the object or the PSF is kept fixed while the other is updated, using one of the usual deconvolution techniques (remember that the convolution product is symmetric with respect to the exchange of object and PSF). Therefore the output of each cycle updates both the object and the PSF, as provided by the previous one. These methods are referred to as iterative blind deconvolution (IBD).

IBD was first proposed with the Wiener filter as deconvolution method (Ayers & Dainty 1988). However, most of the IBD methods are based on RLM (Holmes 1992; Tsumuraya et al. 1994; Fish et al. 1995; Biggs & Andrews 1998). It is obvious that one can use different deconvolution methods for the object and the PSF. For instance, an IBD approach based on the projected Landweber method, with different constraints for the object and the PSF, is proposed in Bertero et al. (1998).

In this paper we extend IBD to the case of multiple-image deconvolution and we apply the method to the reconstruction of LN images. In Sect. 2 we describe the general structure of IBD in terms of cycles: if *p* is the number of detected images, each cycle provides an update both of the object and of the *p* PSFs and is described in terms of an Object box and of *p* PSF boxes. As suggested by the names, the Object box updates the object while the PSFs boxes update the *p* PSFs. We describe the input/output relationship between the boxes and the initialization of the cycles. In Sect. 3 we give the algorithms implemented within the two boxes. In the Object box we use the method denoted ordered subset expectation maximization (OSEM), introduced by Hudson & Larkin (1994) in emission tomography, and extended to LBT imaging in our previous papers (Bertero & Boccacci 2000a; Anconelli et al. 2005a); in the *p* PSF boxes we use single image RLM, with additional constraints on the PSFs. In Sect. 4 we describe the module BLI of the software package AIRY, implementing the IBD method we have developed, and, finally, in Sect. 5 we give the results of our numerical experiments. In Sect. 6 we discuss the limitations of the method, as derived from our experience.

## 2. The structure of the method

We use bold letters for denoting $N \times N$ arrays, whose pixels are indexed by a multi-index $n = \{n_1, n_2\}$. Moreover, we assume to have *p* images acquired with LN, corresponding to *p* different orientations of the baseline and denoted by $g_1, g_2, ..., g_p$. Then, according to the model proposed by Snyder et al. (1993) for images acquired with a CCD camera, the value of one of these images at pixel *n* is given by:

$$g_j(\mathbf{n}) = g_{\text{obj},j}(\mathbf{n}) + g_{\text{back},j}(\mathbf{n}) + r_j(\mathbf{n}), \quad (1)$$

where: $g_{\text{obj},j}(\mathbf{n})$ is the number of photoelectrons due to radiation from the object; $g_{\text{back},j}(\mathbf{n})$ is the number of photoelectrons due to external and internal background, dark current, etc.; $r_j(\mathbf{n})$ is the read-out noise due to the amplifier. The first two terms are realizations of independent Poisson processes (photon noise), so that their sum is also a Poisson process; its expected value is given by:

$$E\{g_{\text{obj},j}(\mathbf{n}) + g_{\text{back},j}(\mathbf{n})\} = (K_j * f)(\mathbf{n}) + b_j(\mathbf{n}), \quad (2)$$

where: $K_j$ is the point spread function (PSF), corresponding to the *j*th orientation of the baseline; *f* is the object array; $b_j$ is the expected value of the (constant) background. We assume, as usual, that the PSFs are normalized to unit volume:

$$\sum_{n} K_j(n) = 1. \quad (3)$$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Pictorial representation of a blind cycle.</figcaption></figure>

The data of the problem are the *p* detected images and backgrounds, {$g_j, b_j$} ($j = 1, ..., p$), and the goal is an estimate of the unknown PSFs $K_j$ ($j = 1, ..., p$), and object *f*.

For the convenience of the reader, in Fig. 1 we give a pictorial representation of one IBD cycle, which consists of one Object box, followed by *p* PSF boxes, one for each PSF. We introduce an index *k* characterizing the IBD cycles. If {$K_j^{(k-1)}, j = 1, ..., p; f^{(k-1)}$} is the output of the cycle *k* − 1 (or the initial estimate in the case *k* = 1), then the functions of the boxes in the cycle *k* are the following.

- **Object box**: this box performs a multiple-image deconvolution. The input consists of {$g_j, b_j$}, ($j = 1, ..., p$) and of the outputs of the previous cycle, {$K_j^{(k-1)}, j = 1, ..., p; f^{(k-1)}$}, while the output will be an update $f^{(k)}$ of the object.

- **PSF boxes**: each box performs a single-image deconvolution for estimating one of the PSFs. For a given *j*, the input of the *j*th PSF box consists again of {$g_j, b_j$} and, in addition, of the estimate of the *j*th PSF, $K_j^{(k-1)}$, provided by the previous cycle, as well as of the estimate of the object, $f^{(k)}$, provided by the previous Object box. The output will be an update of the *j*th PSF, $K_j^{(k)}$.

Since the input data {$g_j, b_j$}, ($j = 1, ..., p$) are common to all boxes, in Fig. 1 they are indicated as inputs of all IBD cycles. Moreover, it is clear that the initialization of IBD is just the initialization of the first Object box which, in Fig. 1, is indicated by {$K_j^{(0)}, f^{(0)}$}.

The deconvolution algorithms implemented in the different boxes, and derived from OSEM and RLM, will be described in detail in the next Section. Therefore we conclude by discussing possible initializations of the first IBD cycle. The initialization requires initial estimates both of the object and of the PSFs.

A possible choice of $f^{(0)}$ is provided by the standard initialization of the RLM algorithm, namely a constant array. For the PSFs the following choices can be considered:

- if the object consists mainly of stars, then, as suggested by Biggs & Andrews (1998), the initial PSFs can be obtained
- Set:

$$ \mathbf{K}_j^{(k,l+1)} = \tilde{\mathbf{K}}_{P,j}^{(k,l+1)} / \sum_n \tilde{\mathbf{K}}_{P,n}^{(k,l+1)}(n). \quad (9) $$

The output of the jth PSF box is just $\mathbf{K}_j^{(k)} = \mathbf{K}_j^{(k,l_{psf})}$, the update of the jth PSF provided by the cycle k.

About the choice of the low-pass filter, we observe that, for any given telescope, the band $\mathcal{B}$ of the instrument, namely the domain of the $\{u, v\}$-plane where the modulus transfer function (MTF) of the telescope is different from zero, is known, since it is determined by the pupil of the instrument. In the particular case of LN, which is the main subject of this paper, $\mathcal{B}$ consists of the union of three discs, which are replicas of the band of a 8.4 m mirror. More precisely, we have different bands $\mathcal{B}_j$ corresponding to the different orientations of the baseline. Then, with reference to the LN-case, given $\mathcal{B}_j$, the most simple choice of the low-pass filter consists in taking $\mathbf{F}_j$ as the inverse Fourier transform of the mask of $\mathcal{B}_j$, namely the array $\mathbf{M}_j$ which is one over $\mathcal{B}_j$ and zero outside:

$$ \hat{\mathbf{F}}_j = \mathbf{M}_j. \quad (10) $$

However, as a consequence of the discontinuity at the boundary of the band, this filter array contains strong side lobes of the central peak, producing strong artifacts.

In order to circumvent this difficulty one can proceed in the following way. First, define an effective band $\mathcal{B}_{eff,j}$ as the domain in the $\{u, v\}$-plane where the MTF is greater than some thresholding value $\sigma$ (i.e. $|\hat{\mathbf{K}}_j(\mathbf{n})| > \sigma$); next, convolve the mask of this domain with some suitable smoothing function such as a Gaussian or a Butterworth function.

## 4. The module BLI of AIRY

On the basis of the ideas presented in the previous Sections, we have developed a module BLI (BLInd deconvolution) of the software package AIRY. As already mentioned, a description of AIRY is given in Correia et al. (2002).

The input of the module BLI is a structure containing the images and additional information (estimate of the backgrounds, pixel size, etc.), while the output consists of two structures: one containing the reconstructed object and the other containing the reconstructed PSFs. In the case of a simulation study, the data are generated inside an AIRY project, while in the case of data reduction application they must be uploaded.

The user has first to create, in the worksheet of CAOS (Correia et al. 2002), an AIRY project linking together a set of modules including BLI. In the worksheet he has to set the global number of iterations which, in this application, coincides with the maximum number of cycles, $k_{max}$, as defined in Sect. 2. In Fig. 2 we show the graphical user interface (GUI) of this module which allows the user to set the parameters of a blind deconvolution application, as described in the previous sections.

About the choice of the PSFs we observe that the first item in the enumeration of Sect. 2 concerns PSFs which are computed by the module while the second and third item must be intended as PSFs provided by the user. Moreover, other parameters which must be selected by the user are: the maximum number of iterations inside the Object box, $l_{obj}$, denoted as *Object iterations* in the GUI; the maximum number of iterations inside the PSFs boxes, $l_{psf}$, denoted as *PSFs iterations* in the GUI. The choice of these parameters will be discussed in the next section.

The choice of the filter is an important issue and can be done through a specific GUI (see Fig. 3), which allows to create the

<figure><img src="image_3.png" /><figcaption>Fig. 2. GUI of the module BLI of AIRY.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. GUI for the choice of the filter.</figcaption></figure>

different types of filters described in Sect. 3 or to upload a user defined filter.
# 5. Numerical experiments

IBD contains several parameters which must be selected by the user, such as the number of iterations within the Object and PSFs boxes, $l_{obj}$ and $l_{psf}$ respectively, and the total number of cycles (or global iterations) $k_{max}$. Following Biggs & Andrews (1998), we will call symmetric IBD the case $l_{obj} = l_{psf}$ and asymmetric IBD the case $l_{obj} \neq l_{psf}$.

Other important degrees of freedom are the initialization of the IBD cycles and the choice of the filter in the PSFs boxes. About the last point, after a number of numerical experiments, we concluded that, in the case of LN, the best choice is provided by a superposition, clipped to one, of three circular Butterworth filters of order 10; the three filters are centered on the three discs forming the band, and have a width slightly smaller than the radius of the discs.

In our numerical experiments we consider objects represented by 256×256 arrays and we assume that they are observed in K-band by means of the LN detector; hence the pixel size is 5.12 mas. Moreover we assume that each observation of a given object provides a set of three images obtained with different orientations of the baseline, corresponding to relative parallactic angles of 0°, 60° and 120°. Therefore, for simulating one observation we need a set of three PSFs. Since we are interested to investigate the dependence of the results on the Strehl ratio (SR) of the PSFs, we consider three different sets, with different SR, each consisting of three PSFs with the orientations of the baseline indicated above. The first is a set of ideal interferometric PSFs of LN, the second and third one are sets of AO-corrected PSFs, with SR of about 70% and 26%, respectively. These PSFs have been obtained by means of the software package CAOS, according to a model already described in previous papers (Anconelli et al. 2004, 2005a).

For a given observation of a given object, the three images are generated first by convolving the object with the selected set of PSFs. Next, sky background emission in K band (12.5 mag/arcsec²) is added to the results, which are also perturbed with Poisson and read-out Gaussian noise (10 e⁻ rms). Moreover the following parameters are used: a total efficiency of 25%, an integration time of 20 min and a telescope surface of 104 m².

With these observational parameters, the average number of background photons per pixel is about $1.4 \times 10^4$; moreover the total number of object photons is about $5.4 \times 10^9$ for an object of magnitude 10, and $2.1 \times 10^9$ for an object of magnitude 11. We remark that the assumed read-out noise, that corresponds to the expected value for the LN detector, is negligible with respect to background and photon noise.

In the following subsections we describe the results provided by the implemented IBD method for two different kinds of objects: binary systems and diffuse objects.

## 5.1. Binary systems

We consider three different binary systems, with different angular separations: 22 mas, 36 mas and 195 mas. The first is just at the limit of the angular resolution of LN in K-band (20 mas), while the two others are well resolved. Moreover we assume that the two stars have magnitude 10 and 11 in K-band, so that, by taking into account background and photon noise, the average SNR per pixel is about 320.

In the following the three binaries will be denoted as Binary 1, 2 and 3, respectively. For each of them we consider three sets of images, corresponding to the three sets of PSFs

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Reconstruction of the binary with angular separation of 36 mas, using exact PSFs (inverse crime). Behaviour of the flux of the two components of the binary as a function of the number of OSEM iterations: ideal PSFs (dotted line); AO-corrected PSFs with SR = 70% (dashed line) and SR = 26% (dash-dotted line). The horizontal full lines give the exact values of the fluxes.</figcaption></figure>

indicated above. Hence, each set simulates an observation with a given AO-correction.

The best accuracy achievable in the reconstruction of the binaries is estimated by an inverse crime experiment, namely the set of PSFs used for generating the images is also used for reconstructing the binaries by means of the OSEM method implemented in AIRY. The accuracy is tested by computing the fluxes of the two components at each iteration, using a square of 3 × 3 pixel centered on each one of the two stars. The results in the case of Binary 2 are shown in Fig. 4. We use a maximum number of 100 iterations and in all cases the accuracy is satisfactory, even if the fluxes are a bit underestimated in the case of Binary 1 and of the PSFs with the smallest SR. Probably one needs a larger number of iterations in this case. Indeed, the main difference between the three cases consists in the convergence rate: convergence is faster for higher values of SR and for wider angular separation of the two stars.

Since we are considering point objects, we first use a symmetric IBD with $l_{obj} = l_{psf} = 1$, as suggested by Tsumuraya et al. (1994) and Biggs & Andrews (1998). In all cases the initialization of the PSFs is provided by the autocorrelations of the images. In the case of ideal PSFs the results are quite good. The rate of convergence depends on the angular separation of the binary but one can conclude that, for all binaries, after about 50 IBD cycles, the fluxes of the two stars, as well as the PSFs, are reconstructed with a few percent error. The correct value of $\Delta m$ is reached after about 25 iterations (errors of the order of 1%).

In the case of AO-corrected PSFs, results are acceptable for Binary 3. For both sets of PSFs, the fluxes of the two stars are underestimated (approximately one half of the correct value) but the difference of magnitude is correctly reproduced (error smaller than 1% after about 30 IBD cycles). Finally the AO-corrected PSFs are reconstructed with good accuracy. In Fig. 5 we show the result we have obtained in the case SR = 26%, for the PSF with an orientation angle of 0°. Important details of the AO-corrected PSF are recovered; the reconstruction error is of the order of 15%.

The results are not so good for the two other binaries. Both the fluxes and the difference of magnitude are underestimated. In the attempt of circumventing this difficulty we consider asymmetric IBD with $l_{obj} > l_{psf}$. In the case SR = 70% we do not find an improvement with respect to the symmetric case. Therefore, for both binaries, the fluxes and the $\Delta m$ are underestimated with an error of about 20%. On the other hand, in the case SR = 26% the best results are obtained with the asymmetry (18 : 1), namely 18 OSEM iterations in the Object box and 1 iteration in the PSF
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** *Upper panels:* the AO-corrected PSF with SR = 26% and 0° (left) and the corresponding image of the binary with angular separation of 195 mas (right). *Lower panels:* autocorrelation of the previous image (left) and reconstruction of the PSF provided by IBD (right).</figcaption></figure>

boxes. The fluxes are underestimated with errors ranging between 30 and 40%, but the $\Delta m$ is overestimated with an error of 20% for Binary 1 and an error of 10% for Binary 2.

It must be observed that the lack of accuracy discussed above concerns photometry; astrometry is quite good in all cases. The reconstructed binaries consist mainly of two bright pixels, corresponding to the correct positions of the two stars, and of a number of weaker scattered pixels (which can be visualized by a representation of the reconstruction in logarithmic scale), accounting for the under-estimation of the flux of the binary. A way for improving photometry could consist in estimating the positions of the two stars and the PSFs by means of IBD and then estimate the magnitudes by a least-square fitting of the detected images, using the procedure proposed in Anconelli et al. (2005b).

## 5.2. Diffuse objects

We first consider a diffuse object already used in previous simulations (Correia et al. 2002), and precisely the image of a young stellar object (YSO) derived from HST near-infrared observation of IRAS 04302+2247. The integrated magnitude of the object is set to 11 in K-band, so that the average SNR per pixel is about 150 (again only background and photon noise is taken into account). Moreover, three sets of LN images are produced by means of the three sets of PSFs mentioned above.

Also in this case we first perform an *inverse crime* experiment in order to estimate the best accuracy which can be achieved by means of the different sets of PSFs. At each OSEM iteration the accuracy is evaluated by computing the relative rms error, defined as the Euclidean norm of the difference between the iterate and the true object, divided by the Euclidean norm of the object. We call this quantity the *reconstruction error*, and in Fig. 6 we plot the results of the first 4000 OSEM iterations for the three sets of PSFs. In all cases the convergence is very slow. Indeed, in the case of ideal PSFs, the minimum reconstruction

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Behaviour of the reconstruction error as a function of the number of OSEM iterations in the case of exact PSFs (inverse crime). The full line corresponds to ideal PSFs, the dotted and dashed lines to AO-corrected PSFs, respectively with SR = 70% and 26%.</figcaption></figure>

error is about 3.7%, and is reached after 7000 iterations, while in the cases SR = 70% and SR = 26%, the minimum reconstruction error is about 4% and 4.8%, respectively, and is reached after 9000 and 20 000 iterations. Therefore, if the quality of the PSFs decreases, both the minimum reconstruction error and the number of iterations increase, even if the quality of the reconstruction is not significantly degraded.

We also checked how good is the data fitting provided by the reconstructed objects. To this purpose, we compute at each iteration the Euclidean norm of $g_j - b - K_j * f^{(k)}$, divided by the Euclidean norm of $g_j - b$, a quantity that can be called *relative discrepancy*. In all cases we find that this quantity decreases monotonically, as a function of the number of iterations, reaching a value of 0.2% in the case of ideal PSFs, and values of 0.62% and 0.65% in the two other cases. In all cases, data fitting is quite good.

The slow convergence in the inverse crime experiment indicates the need of a large number of IBD cycles. Since we are working with an extended object, we use asymmetric IBD with $l_{obj} < l_{psf}$, as suggested by Biggs & Andrews (1998), and $l_{obj} = 1$, to reduce the number of free parameters. Moreover, the PSFs boxes are initialized with the autocorrelations of the ideal PSFs, as proposed in Sect. 2.

In the case of ideal PSFs and AO-corrected PSFs with SR = 70%, the asymmetry $\{1:7\}$ is optimal, since an increase of the asymmetry does not produce an improvement of convergence. In the ideal case, after 500 IBD cycles, the reconstruction error of the object is about 8% (against 3.7% in the case of inverse crime) while the reconstruction error of the PSFs is about 24%. In the case SR = 70%, after 1000 IBD cycles, the reconstruction error of the object is about 13.5% (against 4.0%) while the reconstruction error of the PSFs depends on the PSF and is ranging between 22 and 24%. In Fig. 7 we show the results obtained in this case at the end of the IBD cycles. The reconstruction of the object is quite satisfactory also from the visual point of view and the improvement in the PSF, with respect to the initial guess provided by the autocorrelation of the ideal PSF, is evident. On the other hand, in the case SR = 26% results are very poor: the error in the reconstruction of the object is of the order of 30% (against 4.8%), while the error in the reconstruction of the PSFs is ranging between 40 and 50%.

The failure of IBD in the case of PSFs with low SR may be due to the difficulty of reconstructing an object such as that of Fig. 7, which is characterized by sharp edges in its central part. Therefore we considered a smoother object, namely the image of the galaxy NGC 1288, as extracted from a HST image. Again the integrated magnitude in K band is set to 11 and three LN images
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** *Upper panels:* the object (**left**) and the AO-corrected PSF with SR = 70% and 0° (**right**). *Middle panels:* the image corresponding to the previous PSF (**left**) and the autocorrelation of the ideal PSF with 0° (**right**). *Lower panels:* the reconstruction of the object provided by IBD (**left**) and the reconstruction of the PSF with 0° (**right**).</figcaption></figure>

are generated by convolving the object with the three PSFs with SR = 26%. The same observational parameters of the previous simulations are used.

The inverse crime experiment produces an error of 8% after 3000 iterations. Therefore a large number of IBD cycles must be expected. We find that, in this case, the optimal asymmetry is {1 : 3}; the best reconstruction is obtained with 2000 IBD cycles and the reconstruction error is of the order of 15%, even if the reconstruction error of the PSFs is still quite large, of the order of 40%. In Fig. 8 we show the original object, the simulated LN image with 0°, the result provided by the inverse crime and that provided by IBD. The result looks satisfactory, even if some artifacts due to the “hexagonal” structure of the covering in the *u*, *v* plane are evident. We recall that we consider only three equi-spaced orientations of the baseline; therefore, these artifacts could be reduced by a bit larger number of observations.

**6. Concluding remarks**

In this paper we describe the implementation of IBD in the case of multiple image deconvolution and its application to LBT LINC-NIRVANA interferometer. The method has been inserted in a module of the software package AIRY, version 3.0, which can be downloaded from http://dirac.disi.unige.it or

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** *Upper panels:* the object (**left**) and the image corresponding to the PSF with 0° and SR = 26%. *Lower panels:* the inverse-crime reconstruction (**left**) and the IBD reconstruction (**right**).</figcaption></figure>

from http://www-astro.unice.fr/caos. The results obtained with a limited number of numerical experiments indicate that it is possible to obtain satisfactory results both in the case of stellar objects and in the case of diffuse objects if the SR of the PSFs is sufficiently high. However, it is obvious that the use of IBD in the practice of the reconstruction of the LN images will require a large amount of simulations to optimize the parameters involved in the method for different classes of objects and different observation conditions.

In this first implementation of IBD the methods used in the Object and PSFs boxes are respectively OSEM and RLM. The convergence of these methods is notoriously slow and, as a consequence, a very large number of IBD cycles (global iterations) is typically required before reaching reliable estimates. Therefore it is crucial to investigate the applicability of acceleration methods such as those proposed by Biggs & Andrews (1997, 1998). These acceleration methods must be applied to the global iterations and not to the internal iterations of the Object and PSFs boxes.

We also stress that our software is very flexible, so that other methods, in addition to OSEM and RLM, can be easily inserted in the different boxes. For instance, reduction of noise propagation and artifacts can be obtained by means of Bayesian approaches, leading to the regularization of functionals such as the least-squares functional or the Csiszár I-divergence. It is worth mentioning that iterative methods for the minimization of these functional can be derived from the so-called *split gradient method* (SGM), proposed by Lanteri et al. (2001, 2002) in the case of single image deconvolution, and extended to multiple image deconvolution by Anconelli et al. (2004). Lastly, reduction of boundary effects can be obtained by implementing methods recently proposed by Bertero & Boccacci (2005) and Anconelli et al. (2006) respectively for single image and multiple image deconvolution.
References

Anconelli, B., Bertero, M., Boccacci, P., et al. 2004, Proc. SPIE, 5491, 932

Anconelli, B., Bertero, M., Boccacci, P., Carbillet, M., & Lanteri, H. 2005a,
A&A, 430, 731

Anconelli, B., Bertero, M., Boccacci, P., & Carbillet, M. 2005b, A&A, 431, 747

Anconelli, B., Bertero, M., Boccacci, P., Carbillet, M., & Lanteri, H. 2006, A&A, in press

Ayers, G. R., & Dainty, J. C. 1988, Opt. Lett., 13, 547

Bertero, M., & Boccacci, P. 2000a, A&AS, 144, 181

Bertero, M., & Boccacci, P. 2000b, A&AS, 147, 323

Bertero, M., & Boccacci, P. 2005, A&A, 437, 369

Bertero, M., Bindi, D., Boccacci, P., et al. 1998, Inverse Problems, 14, 815

Biggs, D. S. C., & Andrews, M. 1997, App. Opt., 36, 1766

Biggs, D. S. C., & Andrews, M. 1998, Proc. SPIE, 3461, 33

Carbillet, M., Correia, S., Boccacci, P., & Bertero, M. 2002, A&A, 387, 744

Carbillet, M., Vérinaud, C., Femeniá, B., et al. 2005, MNRAS, 356, 1263

Correia, S., Carbillet, M., Boccacci, P., Bertero, M., & Fini, L. 2002, A&A, 387, 733

Fish, D. A., Brinicombe, A. M., & Pike, E. R. 1995, JOSA, A12, 58

Holmes, T. J. 1992, JOSA, A9, 1052

Hudson, H. M., & Larkin, R. S. 1994, IEEE Trans. Med. Im., 13, 601

Jefferies, S. M., & Christou, J. C. 1993, ApJ, 415, 862

Lanteri, H., Roche, M., Cuevas, O., & Aime, C. 2001, Signal Processing, 81, 945

Lanteri, H., Roche, M., & Aime, C. 2002, Inverse Problems, 18, 1397

Shepp, L. A., & Vardi, T. 1982, IEEE Trans. Med. Im., 1, 113

Snyder, D. L., Hammoud, A. M., & White, R. L. 1993, JOSA, A10, 1014

Tsumuraya, F., Miura, N., & Baba, N. 1994, A&A, 282, 699