2 The reconstruction method

In this work, we used a Lucy-Richardson (LR) iterative algorithm adapted for multiple deconvolution, i.e. an algorithm based on the LR method that allows to retrieve full angular resolution images from LBT interferograms. This adaptation was first developed and tested at the University of Genova (Bertero & Boccacci, in preparation) on simulated data made of one frame per orientation angle.

From this code, we have implemented the possibility to add the information of more than one frame per orientation angle, i.e. to do a simultaneous deconvolution of a set of frames per orientation angle. This approach is useful when dealing with noisy and atmospheric-degraded data, in the case of the LBT-like data we have collected. The method yields the common maximum-likelihood estimate object from the full data set. Details of this algorithm are outlined hereafter.

Let's denote f(x,y) the brightness distribution of the object, and h_i (x,y) the PSF corresponding to the $i^{{\rm th}}$ baseline position angle of the LBT interferometer. For each position angle, we observe the interferogram intensity g_i (x,y) defined by the object image convolution relationship:

$$\displaystyle g_i(x,y) = h_i(x,y) \ast f(x,y) + n(x,y)$$ (1)

where n(x,y) refers to a spatially variant noise process. The original LR algorithm is a non-linear algorithm derived from Bayesian considerations, and based on the knowledge of the PSF. The principle consist in multiplying the result of each iteration by a correction factor C^k(x,y), with k denoting the iteration number, that relates to the remaining fitting error. In standard LR image restoration, i.e. for i=1, the iterative relation is the following:

$$\displaystyle f^{k+1}(x,y) = ~f^{k}(x,y) ~C^{k}(x,y)$$ (2)

with

$$\displaystyle C^{k}(x,y) = ~h^{\ast}(-x,-y) ~\ast ~\left[\frac{g(x,y)}{h(x,y) \ast f^{k}(x,y)}\right]$$ (3)

where $h^{\ast}(-x,-y)$ is the conjugate of h(-x,-y). With a positive constant estimate f⁰(x,y), this algorithm leads to successive estimates f^k(x,y) which are implicitly positive. Since different values of this constant do not lead to significant consequences in term of convergence rate, f⁰(x,y) was fixed to unity. Another interesting characteristic of this algorithm is the conservation of the total energy. It is also demonstrated (Shepp & Vardi [1982]) that this relation leads to the maximum-likelihood object estimate, under the assumption of a Poisson process in image formation.

In the multiple deconvolution case, i.e. for i>1, we simply sum the contribution of each position angle at each iteration. The iterative scheme is identical to Eq.(2) only the correction factor is changed into:

$$\displaystyle C^{k}(x,y) = \frac{1}{H} ~ \sum_{i=1}^N ~ \lef... ...-x,-y) \ast ~\frac{g_i(x,y)}{ h_i(x,y) \ast f^{k}(x,y)}\right]$$ (4)

where N is the number of baseline position angles. 1/H is a normalized term where

$$\displaystyle H= \sum_{i=1}^N\,\left[\sum_{x,y}\,h_i(x,y)\right].$$ (5)

Because of the normalization of each h_i(x,y) to unity, H=N.

The same approach is used when a set of frames per orientation angle g_i,j(x,y) exists, with j denoting the position of the frame in the set. Adaptation of the algorithm to this leaves unchanged Eq.(2) and only modifies the correction factor as following:

$$\displaystyle C^{k}(x,y) {=} \frac{1}{H} \sum_{i{=}1}^N \sum... ...) \ast \frac{g_{i,j}(x,y)}{ h_i(x,y) \ast f^{k}(x,y)} \right]$$ (6)

where M is the number of frames per angle in the set. Note that the PSF may also be different, and therefore written as h_i,j(x,y). But, concerning our data, we were not able to obtain a good knowledge of the corresponding PSF for each frame and therefore preferred a common PSF estimate.

Figure 1: Numerical simulation of interferometric image synthesis with LBT at optical wavelength of a binary star object. a)simulated binary star image with main component magnitude mR=29 and $\Delta m_R=1$. b)one of the noise-free simulated interferograms of the binary (parallactic angle = $45\hbox {$^\circ $ }$). c)same interferogram noise-contaminated (10 linear contours levels are shown). d)reconstructed binary star after 1000 iterations of the algorithm (the z-scale is in electrons)

Figure 2: Same as Fig.1 for a binary star with main component magnitude mR=27.5 and $\Delta m_R=2.5$