The Large Binocular Telescope : Challenging Problems in Image Restoration

The Large Binocular Telescope (LBT), which is currently under construction on the top of Mount Graham in Arizona, will be a Fizeau interferometer with a maximum baseline of 22:8m. From the point of view of image reconstruction, the interesting feature of this instrument is that it will require routinely the use of multiple-images deconvolution methods if the aim is to produce a unique image with the resolution of a 22:8m telescope. Our group is working on this problem since a few years and has already produced methods and software for the data reduction problem. In this paper we provide a survey of our ongoing work.


INTRODUCTION
The Large Binocular Telescope (LBT) will consist of two Ñ mirrors on a common mount, with a spacing of ½ Ñ between their centres, so that a maximum baseline of ¾¾ Ñ will be available.
LBT will be equipped with a Fizeau interferometer, denoted as LINC/NIRVANA, (LN) which will be the result of a German-Italian collaboration.First light of LBT is scheduled for September 2004, second light for November 2005, and the first LN light is foreseen for July 2006.
The interferometric technique used in LN will provide direct imaging with the resolution of a ¾¾ Ñ telescope in the direction of the baseline and of a Ñ telescope in the orthogonal direction.
Since resolution is not uniform over the field, several images of the same scientific object must be acquired with different orientations of the baseline and they must be suitably processed in order to get a unique image with a uniform resolution over the field.Therefore imaging by LN will require routinely the use of multiple-images deconvolution methods if the goal is a unique image with the resolution of a ¾¾ Ñ telescope.
According to our experience in image reconstruction, two different kinds of methods will be required: ¯quick-look methods, computationally efficient even if not very accurate for any kind of astronomical object; they can be routinely used for a preliminary view of the target just after data acquisition.
¯ad hoc methods designed for particular classes of astronomical objects; each one of them must be accurate as far as possible for the objects in the corresponding class, even if it does not need to be computationally efficient.
Therefore the processing of LBT images will require, in general, a two-steps approach: the first step, based on a quick-look method, will be intended to understand qualitatively specific features of the particular object the astronomer is investigating; the second one will consist in the use of the ad hoc method which has been designed for the class of objects with those specific features.
Our group has already proposed possible quick-look methods (Bertero and Boccacci [1,2], Anconelli et al. [3]); they have been implemented in the Software Package AIRY (Astronomical Image Restoration in interferometrY), which is described in Correia et al. [4].AIRY is a set of modules for the simulation and/or analysis of interferometric images; it is written in IDL and is conceived to be used with the CAOS (Code for Adaptive Optics Systems) Application Builder (see ØØÔ »»ÛÛÛº Ö ØÖ º ×ØÖÓº Ø» Ó×).Validation of the implemented deconvolution methods is described in Correia et al. [4] and Carbillet et al [5].The investigation of ad hoc methods is in progress.One of them is proposed in Anconelli et al. [6].
The paper is organized as follows.In Section 2 we describe a general approach to the design of iterative methods for constrained deconvolution.In Section 3 we apply the approach to the derivation of fast RL-like methods (RL=Richardson-Lucy) which are candidates to be quick-look methods.The same approach is used in Section 4 for deriving three ad hoc methods we are presently investigating: the first one is for objects with high dynamic range; the second one is for faint objects if both photon and read-out noise must be taken into account; the third one is for the super-resolution of unresolved stellar objects.Finally in Section 5 we briefly discuss the state-of-the art of our activity.

THE GENERAL APPROACH
In the following we use bold letters for AE ¢ AE arrays.Therefore we denote by ½ ¾ Ô the Ô LN images of a given astronomical object, corresponding to Ô different orientations of the LBT baseline; moreover we often use the more compact notation for indicating the set of the Ô images.Each one of these images, acquired with a CCD camera, can be described by the following model (Snyder et al. [7]): where ´Ñ Òµ denotes the value of the detected -th image at pixel Ñ Ò; it is given by the sum of the number of photoelectrons due to both object radiation, Ó ´Ñ Òµ, and background radiation, ´Ñ Òµ, with the addition of the read-out noise due to the amplifier, Ö ´Ñ Òµ.The first two terms are realizations of independent Poisson processes, so that their sum is also a Poisson process, while the third term is the realization of a Gaussian process (white noise).The expected value of the -th image is given by: where: is the unknown object, an array formed by the average numbers of photons emitted at the pixels of the object domain and collected by the telescope; is the matrix describing the imaging process (optics plus atmosphere); is a background term, which can include a nonzero expected value of the read-out noise and, in general, can be assumed to be constant over the image domain.If we assume, for simplicity, space invariance, the matrix is given in terms of the point spread function (PSF) Ã of the instrument, corresponding to the -th orientation of the baseline: Ã £ (3) moreover we assume, as usual, that the PSFs are normalized in such a way that the sum of all their pixel values is one: The problem of image restoration is the problem of estimating the unknown array , being given the detected images , the backgrounds and the PSFs Ã .It is well known that it can be reduced to the minimization of suitable cost functionals in the framework of approaches such as regularization theory, maximum likelihood or Bayesian estimation (see, for instance, Bertero and Boccacci [8] for an introduction).In the literature there exists a great variety of cost functionals and of iterative algorithms for their minimization, with the possible addition of suitable constraints (for instance, non-negativity of the solution).All these methods have been designed for the deconvolution of a single image and in some cases their extension to the multiple-images case is not trivial.Certainly the software developed for these methods cannot be directly applied to multiple images.
In a series of recent papers Lanteri et al. [9,10,11] propose a general approach, denoted as SGM (Split Gradient Method), to the design of iterative methods for the minimization of functionals of the following type: with the additional constraints: In these equations Â ¼ ´ µ is the functional measuring the discrepancy between the computed images associated to and the detected images , and Â Ê ´ µ is a regularization functional while is the regularization parameter.As concerns the constraints, the first one is obvious while the second one is a constraint on the total flux of the object.The choice of will be discussed in Section 3.
The basic idea in the approach of Lanteri et al. relies on the following decomposition of the gradient of the functional Â ´ µ: ÖÂ ´ µ Í ´ µ Î ´ µ (7) where Í ´ µ and Î ´ µ are positive arrays depending on .Obviously such a decomposition always exists and is not unique.The applicability of the method requires explicit expressions for the dependence of these arrays on .Then the general structure of the iterative algorithm, as described in Lanteri et al. [9], is as follows: ¯choose an initial ´¼µ , satisfying the constraints of Eq. 6 (in general, a uniform image); ¯given ´ µ compute: ´ ¯set: ´ In Eq. 8 products and quotients of arrays are defined as pixel-by-pixel products and quotients while « and are relaxation parameters which, in principle, can be changed at each iteration.« is the step in the descent direction (which is modified by if ½), and it can be computed in order to guarantee both the non-negativity of the iterates and the convergence of the iterations.On the other hand is a parameter which can allow a reduction of the number of iterations.We will call it the Ð Ö Ø ÓÒ ÜÔÓÒ ÒØ.Indeed, as discussed by Lanteri et al. [9], an integer value of ½ can speed-up the convergence of the iterations.As already remarked, at the very first iterations modifies the descent direction; but, when the iterations are close to convergence, the algorithm with ½ is approximately equivalent to the algorithm with ½ and step-size « .
The algorithm takes a very simple form if we choose a unit step (i.e. « ½); in such a case we obtain: ¯given ´ µ compute: ´ •½µ ´Ñ Òµ (10) ¯set: ´ The convergence of this algorithm is not guaranteed in general, even if it can be proved in some particular cases; nevertheless it has been verified experimentally in all the applications we have considered.The most interesting feature of the algorithm is that non-negativity of the iterates is automatically satisfied, as one can easily verify.
In the case of this simplified algorithm the dependence on the regularization parameter can be made explicit by introducing decompositions of the gradients of the functionals Â ¼ ´ µ and Â Ê ´ µ similar to that of Eq. 7: then it turns out that the algorithm has the following structure: ¯given ´ µ compute: ´ ´ which shows the very simple dependence of the algorithm on the regularization parameter .
We apply now the previous approach to the case of multiple images.In such a case the functional Â ¼ ´ µ, which measures the discrepancy between the set of the detected images and the set of the computed images associated to , is, in general, the sum of the functionals corresponding to the different images.In a statistical approach this property follows from the statistical independence of the different images.Therefore we have: and the corresponding decompositions of the functions Í ¼ ´ µ and Î ¼ ´ µ follow from the equation: If we insert these expression in Eq. 14 and if we remark that the gradient of the regularization term does not depend on , so that it can be written as the sum of Ô identical contributions, we find that the complete algorithm takes the following form: ¯given ´ µ compute: ´ ´ The functions Í¼ ´ ´ µ µ and Î¼ ´ ´ µ µ for the two statistics considered in the text.We denote by ½ the array whose entries are all equal to 1.

¯set:
´ Eq. 18 shows the dependence on both the regularization parameter and the detected images.
The interesting feature of this algorithm is that the numerator and the denominator in Eq. 18 are given by a sum of terms which are just those appearing in the corresponding algorithm for singleimage deconvolution.This structure suggests a cyclic-iterative algorithm which is analogous to the OS-EM algorithm proposed by Hudson & Larkin [12] in the case of tomography: ´ • ½µ ÑÓ Ô and compute: ´ ´ We will call this algorithm the OS version of the multiple-images algorithm of Eq. 18.It is obvious that the computational cost of one of its iterations is just that of one iteration of the same algorithm in the case of a single image.Therefore, if the number of iterations of the algorithms of Eq. 18 and Eq.20 is approximately the same (a property which has been verified experimentally in several cases), it follows that the algorithm of Eq. 20 provides a computational gain by a factor which is approximately Ô. Finally it is obvious that the algorithm for single-image deconvolution is obtained by just taking the same image in each iteration.
We conclude this Section by showing the application of the previous approach to a few particular cases.In particular we consider two different kinds of noise statistics, leading to different discrepancy functionals: ¯Poisson Statistics (photon counting) -In this case the functionals Â ¼ ´ µ are given by: ¯Gauss statistics (white noise) -In this case the functionals Â ¼ ´ µ are given by: In Table 1 we give the expressions of the arrays Í ¼ ´ ´ µ µ and Î ¼ ´ ´ µ µ for the two cases.
In deriving these formulas we used the normalization of the PSF, as given in Eq. 4, which implies that Ì ½ ½ and Ì .
If we insert these expressions into Eq.18 with ¼ we obtain respectively the multi-image extensions of RL and ISRA (Iterative Space Reconstruction Algorithm), as proposed in Bertero & Boccacci [1,2].Moreover, by inserting the expressions for the Poisson case into Eq.20, with ½, we re-obtain the OS-EM method proposed in Bertero & Boccacci [1] for accelerated multiimages deconvolution.If ½ we obtain the further accelerated version of OS-EM investigated in Anconelli et al. [6].This point will be reconsidered in the next Section.
Let us emphasize the fact that the most interesting feature of the approach proposed in this paper is that one can easily insert regularization terms, obtaining new iterative algorithms which can be easily implemented.We mention, among others, Tikhonov regularization and Maximum Entropy both for Gauss and Poisson noise (with the addition of background terms, non-negativity and flux constraint).In Table 2 we give the expressions of the arrays Í Ê ´ µ and Î Ê ´ µ for the following regularization functionals: where ´Êµ ´Ñ Òµ denotes a reference object (in general, a constant array; very often, zero). ¯Smoothing where is a matrix with non-negative entries; such a form applies, for instance, to regularization in terms of discrete Laplacian, which can just be written in the form ¡ ´Á µ (Á being the identity matrix). ¯Entropy where, again, ´Êµ ´Ñ Òµ denotes a reference object.
The functions ÍÊ´ µ and ÎÊ´ µ in the case of Tikhonov, Smoothing and Entropy regularization.
In the case of Entropy the log of an array is defined pixel-by-pixel; moreover the flux normalization is taken into account, being the flux constant defined in Eq. 32.
If we combine the arrays of Table 1 and Table 2 we obtain several iterative algorithms for different kinds of regularization both in the Poisson and in the Gauss case.

QUICK-LOOK METHODS
As it is well-known the RL method is one of the most frequently used methods for image deconvolution in Astronomy.It has many interesting features: ¯it is an iterative method which maximizes the likelihood function in the case of Poisson statistics; ¯it provides automatically nonnegative objects; ¯in the case of zero background it provides objects whose total flux coincides with that of the detected image.
The main difficulties are generated by the slow convergence of the algorithm and to the fact that, in some cases, iterations must not be pushed to convergence.For instance, in the case of a diffuse object, such as a nebula or a galaxy, the iterations first provide better and better approximations of the object but, after a certain point, the so-called checkerboard effect appears, which is basically due to amplified noise propagation.More precisely, in the case of numerical simulations, let us define, at each iteration , the following relative restoration error ´ µ Ö Ð : ´ µ Ö Ð ´ µ where the numerator is the r.m.s.value of the difference between the -th iterate and the true object while the denominator is the r.m.s.value of the true object.Then, if we look at the behaviour of this restoration error as a function of the number of iterations, we find that it first decreases, goes trough a minimum and then increases, showing a typical behaviour which is known as semiconvergence (see, for instance, Bertero and Boccacci [8]).On the other hand, in the case of pointwise objects such as a binary system, the algorithm seems to converge.An example is given in Fig. 1.In the left panel we plot, in the case of the diffuse object of Fig. 2, the restoration error of Eq. 27 as a function of the number of iterations; the minimum of the restoration error is evident even if it is very flat.It is obvious that one can use about 100 iterations without a significant loss of accuracy.On the other hand. in the right panel we plot the behaviour of the reconstructed flux of a binary system, with Ñ ¾, again as a function of the number of iterations.It is evident that, at least as concerns photometry, we have a convergent behaviour and, in such a case, the convergence is very fast since we need about 20 iterations to get an accurate value.
If we apply Eq. 18 to the case of Poisson noise with ¼, we obtain the following algorithm for multi-images deconvolution: ADA III ¯given ´ µ compute: ´ ´ As we already remarked in the previous Section, in the case ½ this is precisely the multiimages extension of the RL-method proposed in Bertero and Boccacci [1].In a similar way from Eq. 20 we obtain the following algorithm: ¯given ´ µ , set ´ • ½µ ÑÓ Ô and compute: ´ ´ As concerns flux normalization, which is taken into account by all our algorithms, the constant can be selected in the following way.In the case of zero background, as shown in Bertero & Boccacci [1], the total flux of each iterate of the multi-image RL algorithm coincides with the arithmetic mean of the total fluxes of the images .As a consequence of the introduction of the background terms, the iterates of the corresponding RL-algorithm and of its accelerated version OS-EM, both given in Bertero & Boccacci [2], do not have this property.Therefore, it is convenient to constrain the solution to satisfy the flux condition of the RL-method with zero background; it follows that the constant can be estimated as the arithmetic mean of the total fluxes of the images , after subtraction of the background terms: In Fig. 2 we give an example of reconstruction of a diffuse object obtained with the OS-EM method.
The object is the HST image of the Southern Crab Nebula.In the simulation we have used three images obtained by convolving the object with three ideal interferometric PSFs corresponding to equally spaced orientations of the baseline and by adding Poisson and Gaussian noise.More details about our simulations, which are performed by means of the Software Package AIRY, are given in Anconelli et al. [3].We have considered several different diffuse objects with different integrated magnitudes and different sets of PSFs, including different levels of adaptive optics (AO) correction.
The plot of the restoration error corresponding to the object of Fig. 2 is just that given in the left panel of Fig. 1.The number of iterations required for getting the best restoration is not terrible but is not small.Therefore it is interesting to investigate the effect of the acceleration exponent .As concerns this point we can summarize the results of the numerical simulations described in Anconelli et al. [3].The tests have been performed mainly on diffuse objects by looking at the number of iterations needed for reaching the minimum of the restoration error.As a result the choice ¾ always leads to a converging algorithm, with a reduction by a factor 2 in the number of iterations, hence a gain of a factor 2 in computation time.Moreover by changing the magnitude of the object we have verified that, when the magnitude increases, one can also increase the acceleration exponent.In some cases we have obtained converging algorithms with and a corresponding reduction of the number of iterations by a factor 8. Therefore such an approach seems quite promising for producing quick-look methods.

METHODS FOR PARTICULAR CLASSES OF OBJECTS
As we stated in the Introduction we expect that, in general, an observation of a complex object with LBT will require two steps: the first one is a reconstruction with a quick-look method which can be used for identifying the object and understanding qualitatively its main features; the second one is a reconstruction based on a specialized method which has been designed for objects of that type.In this Section we give a few examples which may help to understand our point of view.

OBJECTS WITH HIGH DYNAMIC RANGE
An important feature of Astronomical images, which is not present in the natural ones, is that, in many important cases, they can have a very high dynamic range, i. e., in the same image one can find objects with very different magnitudes, in general associated with different scales.One example is provided by a model of young binary star with scattered light emission coming from a dusty circumbinary ring.This model is proposed in Carbillet et al. [5] and is produced using parameters taken from observations of the close binary of the quadruple GG Tau system in the near infrared.The reconstruction of this object is a difficult problem because the ratio between the flux of one of the two stars and the maximum value of the flux emitted by the pixels of the ring is of the order of ½¼ .
The object is shown in the upper-left panel of Fig. 3.The main component has K-magnitude Ñ=19 and the difference of magnitude is ¡Ñ ½ while the integrated magnitude of the ring is about 15.The corresponding image is represented in the upper-right panel using a logarithmic scale.This image has been obtained by convolving the object with an Airy function and by corrupting the result with Poisson and Gaussian noise.It is evident that the rings of the images of the two stars are superimposed to the ring and that they mask its structure.Finally in the lower-left panel we show the best reconstruction obtained with the RL method, in the sense that it corresponds to the minimum of the reconstruction error defined in Eq. 27 and computed on the region of the ring.Several artifacts are evident; if iterations are further pushed, then the image of the ring is corrupted by the typical checkerboard effect of the RL method when applied to the reconstruction of diffuse objects.
It is obvious that the RL-reconstruction already provides important information which cannot be directly deduced from the original image.On the other hand it is also obvious that the reconstruction of the ring is not completely satisfactory for a quantitative analysis.Moreover, by  looking at the behaviour of the iterations, it seems clear that a method acting in a different way in the region of the binary and in that of the disc should be desirable.Such a remark suggest to use a kind of adaptive regularization which can be obtained, for instance, by adding to the discrepancy of the photon noise, as given in Eq. 22, the following regularization functional: This functional, with replaced by the gradient of , was introduced by Geman and Mc Clure [13] in order to get a different regularization in regions with very different intensity variations.It is not convex but is bounded from below and its Hessian has the same property.It has a unique minimum at ¼ and is convex in a neighborhood of this point.The parameter is basically a thresholding parameter which separates the low-flux regions from the high-flux ones.Indeed it is easy to recognize that, if ´Ñ Òµ is much smaller than , then the functional Â Ê ´ µ practically coincides with the Tikhonov regularization functional of Eq. 24; on the other hand, if ´Ñ Òµ is much greater than , then the functional is constant and therefore no regularization of the RL method is introduced.
If we have in mind the object discussed above, we easily understand that, even if a preliminary reconstruction of the object by means of the RL method, does not provide an accurate reconstruction of the ring, however it provides sufficient information for attempting an estimation of the parameter .
Our problem is now to find an iterative algorithm for the minimization of the functional of Eq. 5, with the functional Â ¼ ´ µ given by Eqs.16 and 22, and the functional Â Ê ´ µ given by Eq. 33.
In order to simplify our equations, we write the regularization parameter in the form: then, by applying SGM to this functional we obtain a very simple modification of the algorithm of Eq. 28 (with ´ • ½µ ÑÓ Ô and compute: ´ µ ´ µ ´ In the case of a single image we have applied this algorithm to the object of Fig. 3.We have found that the best restoration is obtained with ¿ ¢ ½¼ and ¢ ½¼ ½¼ .The result is shown in the lower-right panel of Fig. 3; the r.m.s error on the ring is of the order of ±.However, in order to make evident the accuracy of the reconstruction, in Fig. 4 we compare the reconstructed object with the original one by looking at the level curves.Except for a small region the accuracy is quite good.Work is in progress for investigating the applicability of this approach to LBT images.The results published in Carbillet et al. [5] indicate that the result of the first step, based on the OS-EM method, is not so reliable as that obtained in the case of a single image and shown in Fig. 3.

FAINT OBJECTS
In the observation of faint objects such that only a few photons per pixel are available, it may happen that the read-out noise is comparable with the photon noise.In such a situation the approximation of the likelihood function based on the assumption of a dominating photon noise is not valid and one must use the more complete expression which takes into account both Poisson and Gaussian noise, as given by Snyder et al. [7].If we consider the negative ÐÓ of this function and we modify the terms independent of in order to get an expression closer to that of the Poisson case, we obtain the following form of the discrepancy functional for the th image: È ÖÓÒ ´ µ being the probability density of the read-out noise, which, in the case of white Gaussian noise with expected value Ö and variance ¾ , is given by: The SGM can be applied to this functional and, as shown by Lanteri and Theys [14], the expressions of the U, V arrays are the following: where: By inserting these expressions into Eq.18 or Eq. 20, both with ¼, we obtain the iterative algorithm for the multiple-images case and its OS version.The latter is just the algorithm for the single-image case when only one image is available.
We point out that we do not consider the case where a regularization functional is added to the discrepancy functional because, as shown in Lanteri and Theys [14], the addition of standard regularization functionals such as Tikhonov or Entropy does not produce a significant improvement of the restored image.Moreover the convergence of the algorithm is very fast since very few iterations (between 5 and 10) are needed for reaching the best image.Therefore an interactive choice of the best reconstruction, based on the visualization of the result of each iteration, is possible As an example we show here the result of a numerical simulation described in Lanteri and Theys [14].The object is a ½¾ ¢ ½¾ HST image of a sun-like star nearing the end of its life.The image has been convolved with a PSF including atmospheric turbulence effects and a total flux of ½¼ photons has been assumed (average number of photons per pixel smaller than 1); moreover a Gaussian white noise with ½ has been added.The image has been reconstructed using both the RL-method and the algorithm outlined above.The results are shown in Fig. 5 where we show the object, the blurred and noisy image as well as the reconstruction provided by the RL-method and that provided by the new method.The visual improvement obtained with the new method is evident, even if the improvement in the r.m.s.error is not very large.Indeed the minimum restoration error of the RL-method is 3.9 % and is reached after 8 iterations, while the minimum restoration error of the new method is 3.5 % and is reached after 6 iterations.

SUPER-RESOLUTION OF STELLAR SYSTEMS
As another example of a two-steps approach for an optimal reconstruction of LBT images we briefly describe our recent work on the super-resolution of unresolved objects.
In Astronomy and Microscopy super-resolution is a term which, in general, is used to indicate any method able to improve resolution beyond the diffraction limit, and therefore it implies the extrapolation of the Fourier transform of the object outside the band of the optical instrument.As remarked in a paper by Wolter [15], this is possible if the object has a finite extent.However, since the extrapolation problem is very sensitive to noise propagation from the data to the solution, unlimited super-resolution is impossible in the presence of noise.A first attempt of estimating the amount of super-resolution achievable in practice is developed in Bertero and Pike [16], where it is shown that it is controlled by two parameters: the space-bandwidth product (SBP) and the signal-to-noise ratio (SNR).In the case of a telescope the SBP is approximately given by the ratio between the angular diameter of the object and the diffraction limit; then the results derived in the paper mentioned above imply that super-resolution is feasible when this parameter is not much greater than one, namely in the case of compact and unresolved objects.On the other hand the amount of super-resolution increases with increasing SNR, as shown by Lucy [17].
According to the previous remarks, the development of a super-resolving method requires an estimate of the domain of the object, which must be inserted as a constraint in the restoration algorithm so that the total flux of the restored object is concentrated within .However it has been claimed by many authors that other constraints may provide super-resolution such as nonnegativity and a constraint on the total flux of the object.For instance, it has been proved that maximum entropy provides super-resolution in the case of the so-called nearly-black objects, i.e. objects which are essentially zero in the vast majority of the samples (Donoho et al. [18]) and that, in general, non-negativity works well when combined with the sparsity of the object (Donoho [19]).For sparse objects the constraint on the total flux can provide improved reconstructions.
The RL-like methods introduced in the previous Sections implement both the constraint of nonnegativity and that on the total flux of the image.Moreover, they can implement also the constraint on the domain of the object if we use a suitable initialization of the iterations.Indeed, this property derives from the fact that the result of iteration • ½ contains as a factor the result of iteration ; therefore if ´¼µ ¼ at a given pixel then all the iterates will be zero at that pixel.We call such a property the localization property of the RL-like methods, in the sense that we can constrain the object to be localized in a given domain by a suitable initialization of the iterations.Therefore RL-like methods provide algorithms implementing in a very simple way all the relevant constraints and they look quite suitable for estimating the amount of super-resolution achievable in practical situations.
Our approach to super-resolution applies to the case where the astronomical target consists of unresolved objects (angular size of the order of the diffraction limit) separated by an angle considerably greater than the diffraction limit (sparsity condition), so that there is no overlapping of their images.For simplicity we assume that there is only one of these unresolved objects, surrounded by empty sky.Then the approach consists of the following steps.
Step 1 -We apply the RL-like method initialized with a constant image (correctly normalized) and we use a number of iterations such that the restored object is sufficiently well localized.If the angular separation is not much smaller than the diffraction limit and the magnitude difference is small, it may happen that the stars of the object are already resolved in this first step, as a consequence of the super-resolving properties of the RL-like methods.In such a case one can estimate their positions by computing the centroids and go directly to Step 3 for a more accurate estimate of their magnitudes.
If the stars are not resolved, then one can go to Step 2.
Step 2 -We define the domain of the object by identifying the pixels where the flux of the result of Step 1 is greater than a selected threshold (for instance some percent of its maximum value); in alternative we can take a disc with diameter equal to the diffraction limit and containing most of the flux of the reconstructed object.Next we apply again the RL-like method to the detected images, but now initialized with the mask of the domain, namely a function which is constant over the domain and zero elsewhere, The constant is evaluated in such a way that the normalization condition is satisfied.Each iterate is localized in the selected domain; if, after a number of iterations, the stars are separated, their positions are obtained by computing their centroids.
In order to illustrate this approach, in Fig. 6 we show an example of its application.The object consists of two unresolved binaries with an angular separation of their centres which is about on-half of the diffraction limit (here we consider the case of LBT).Also for this object we have considered three images obtained by convolving the object with three ideal PSFs of LBT, corresponding to uniformly spaced directions of the baseline; the images have been corrupted with photon and read-out noise.The result of 1000 OS-EM iterations allows to resolve the two binaries but not the stars within the binaries.Next a mask is obtained by identifying the pixels where the flux is greater than 50 % of the maximum value of the restored image.This mask is also shown in Fig. 6.Finally OS-EM is reused, but now initialized with the mask derived from the first step, and the result of 1000 iterations clearly shows now the structure of the two binaries.
A third step may be required when the unresolved object consists of stars with different magnitudes.Indeed, in such a case, the first two steps may not provide the correct photometry.In order to overcome this difficulty we propose the following third step.
Step 3 -Let us assume that the restored object, provided by Step 1 or Step 2, consists of Õ stars localized at the pixels È ½ È ¾ È Õ (in the application of this paper Õ ¾).Then the positions of the stars are deduced from these restorations while their intensities are obtained by solving a linear least-squares problem.
In Anconelli et al. [6] this approach is applied to the case of binary systems and, by means of a careful analysis, it is shown that a considerable amount of super-resolution is achievable both in the case of a single mirror and in the case of a Fizeau interferometer such as LBT.

CONCLUDING REMARKS
In this paper we have presented the state-of-the art of our research about the reconstruction of LBT images as well as the many problems this unique instrument is suggesting.All the methods we have described are implemented in the software package AIRY which is now available in the version 2.1 and can be required at the site ØØÔ »» Ö º × ºÙÒ º Ø. Work is in progress about the validation of the methods described in the Sections 4.1 and 4.2.

FIGURE 1 :
FIGURE 1:Illustrating the convergence of the RL-algorithm: to the left, the plot of the behaviour of the reconstruction error, defined in Eq. 27, as a function of the number of iterations, in the case of a diffuse object; to the right, the plot of the behaviour of the flux of the reconstructed image in the case of a binary system.

FIGURE 2 :
FIGURE 2: An example of diffuse object used in our numerical simulations: in the left panel, the object, a HST image of the Southern Crab Nebula; in the central panel, one of the three LBT images used for the reconstruction; in the right panel, the reconstructed object

If½,
we get the version of the OS-EM algorithm proposed in Bertero and Boccacci[1] while, if ½, we get the accelerated version of the same algorithm proposed and investigated in Anconelli et al.[3].All the algorithms described above are implemented in the present version of the Software Package AIRY.

FIGURE 3 :
FIGURE 3: Upper-left: the model of young binary system used as a test-object with high dynamic range; upper-right: the image o/btained by convolving the object with an Airy function and by corrupting the result with Poisson an Gaussian noise; lower-left: the best reconstruction obtained with the RL-method; lower-right: the best reconstruction obtained with the adaptive regularization proposed in this paper.

FIGURE 4 :
FIGURE 4: Left panel: level curves of the object of Fig.3; right panel: level curves of the reconstruction obtained by means of the method of Eq. 34.

FIGURE 6 :
FIGURE 6: Upper-left: the object consisting of a pair of binaries (the diameter of the circle is just ); upperright: the OS-EM reconstruction (1000 iterations); lower-left: the mask; lower-right: the OS-EM reconstruction initialized with the previous mask (1000 iterations).