2 Formulation of the problem

If we assume that p observations, corresponding to p different orientations of the LBT baseline, have been performed and that for each observation the PSF is space-invariant, then the mathematical model describing the formation of the p images g_j(x,y), (j = 1,..., p)is the following

 

$$\int\int K_j(x-x^\prime,y-y^\prime)f_0(x^\prime,y^\prime){\rm d}x^\prime {\rm d}y^\prime$$ g_j(x,y) =

$$+ \,\, b_j(x,y) + w_j(x,y)$$ (1)

where f₀(x,y) is the brightness distribution of the unknown astronomical object, K_j(x,y) is the space-invariant PSF of the LBT interferometer, b_j(x,y) is a function describing the average background due to sky emission and, finally, w_j(x,y) represents the noise (typically Poisson noise and read-out noise) corrupting the j-th image.

The image restoration problem consists in evaluating an estimate f(x,y) of f₀(x,y), being given the detected images g_j(x,y), the PSFs K_j(x,y) as well as estimates of the average backgrounds b_j(x,y). The latter estimates are important in many circumstances. Indeed, if the restoration method includes the positivity constraint, such a constraint is not active without background subtraction. In general one can assume that the functions b_j(x,y) are constant over the field of view of the telescope.

In practice the LBT images are discrete and consist of $N\times N$ pixels, which we characterize by the indices m,n = 0, 1,..., N-1. If we denote by g_j(m,n), f₀(m,n) etc. the value of g_j, f₀ etc. at the pixel m,n, then the discrete version of Eq. (1) is given by

 

$$\sum_{m^\prime,n^\prime=0}^{N-1} K_j(m-m^\prime,n-n^\prime)f_0(m^\prime,n^\prime)$$ g_j(m,n) =

$$+ \,\, b_j(m,n)+ w_j(m,n).$$ (2)

In many circumstances, and precisely when the size of the central peak of the PSF is much smaller than the size of the image, a sufficiently good approximation is obtained if the convolution product of Eq. (2) is intended as a cyclic convolution (periodic PSF), so that, by taking the Discrete Fourier Transform (DFT) of both sides of this equation, we obtain the usual relationship

$$\hat g_j(m,n)=\hat K_j(m,n)\hat f_0(m,n)+\hat b_j(m,n)+ \hat w_j(m,n),$$ (3)

where the hat denotes the DFTs. The DFTs of the discrete PSFs K_j(m,n) are the discrete transfer functions (TF) of the system corresponding to the different orientations of the baseline.

Since the PSFs of LBT are bandlimited, if they are correctly sampled then the discrete TFs should be zero outside the set of pixels corresponding to the band. However, in practical applications the PSFs of Eq. (2) are measured (for example by imaging a guide star). In such a case the DFTs of these PSFs are, in general, different from zero everywhere, as a consequence of experimental errors and noise. Then the band of the j-th PSF can be defined as the set of the pairs of indices m,n where $\hat {K}_j(m,n)$ is greater than some threshold value related to the experimental errors. We denote this set by ${\cal B}_j$ and we assume that $\hat {K}_j(m,n)$ is set to zero outside ${\cal B}_j$.

In this section we formulate the image restoration problem for LBT as a least-squares problem. This is equivalent to a Maximum Likelihood approach in the case of white Gaussian noise.

In order to formulate this problem and the related equations in a concise form, we introduce a vector-matrix notation and we write Eq. (2) as follows

$$\vec{g}_j =\vec{A}_j \vec{f}_0 +\vec{b}_j + \vec{w}_j$$ (4)

where $\vec{g}_j, \vec{f}_0$ etc. are the arrays formed by g_j(m,n), f₀(m,n) etc. and $\vec{A}_j$ is the block-circulant matrix defined by the cyclic convolution with the j-th PSF. Then, if we denote by $\vec{b}_j$ an estimate of the j-th background contribution (obtained, for instance, by means of some suitable preprocessing of the detected images $\vec{g}_j$), a least-square estimate of $\vec{f}_0$ is any array $\vec{f}$ which minimizes the discrepancy functional:

$$\varepsilon^2(\vec{f})= \sum_{j=1}^p \Vert\vec{A}_j\vec{f}+\vec{b}_j - \vec{g}_j\Vert^2$$ (5)

where $\Vert .\Vert$ denotes the usual Euclidean norm. This equation implies that it is convenient to introduce subtracted images $\vec{g}_{s,j}$ defined as follows

$$\vec{g}_{s,j}=\vec{g}_j-\vec{b}_j.$$ (6)

By means of standard variational procedures (see I), it is easy to prove that to minimize the functional of Eq. (5) is equivalent to solve the following normal equation

$$\sum_{j=1}^p \vec{A}_j^{\rm T} \vec{A}_j \vec{f}= \sum_{j=1}^p \vec{A}_j^{\rm T} \vec{g}_{s,j},$$ (7)

where $\vec{A}_j^{\rm T}$ denotes the transposed of the block-circulant matrix $\vec{A}_j$. Then, by taking the DFT of both sides of this equation one easily gets (the ^* denotes complex conjugation)

$$\mid \hat K(m,n)\mid ^2 \hat f(m,n)=\sum_{j=1}^p \hat{K}^\ast_j(m,n) \hat g_{s,j}(m,n)$$ (8)

where

$$\mid \hat K(m,n)\mid ^2=\sum_{j=1}^p \mid \hat {K}_j(m,n)\mid ^2 .$$ (9)

We point out that the band of $\vert\hat K(m,n)\vert$ is the set $\cal B$ which is the union of the sets ${\cal B}_j$ associated with the PSFs corresponding to the different orientations of the baseline. This set represents the coverage of the u-v plane provided by the observations we are considering.

Equation (8) implies that $\hat f(m,n)$ is uniquely determined inside but not determined outside the band $\cal B$ defined above. In other words the solution of the equation is not unique because the images $\vec{g}_{s,j}$ do not provide information about the DFT of the object outside $\cal B$. In such a case it is usual to introduce the least-squares solution with minimal Euclidean norm, the so-called generalized solution (Golub & van Loan 1983), denoted by $\vec{f}^\dagger$, which is obtained by setting $\hat f(m,n) = 0$ outside $\cal B$. Therefore the DFT of $\vec{f}^\dagger$ is given by

$$\hat f^\dagger(m,n)=\sum_{j=1}^p {{\hat K_j^\ast(m,n)}\over{\mid \hat K(m,n) \mid^2}} \hat g_{s,j}(m,n)$$ (10)

when the pair m,n belongs to $\cal B$ and is zero elsewhere. By taking the inverse DFT of this array, one gets $\vec{f}^\dagger$.

However the generalized solution is not a satisfactory solution of the restoration problem. In the case of a single image it coincides with the solution provided by the so-called inverse filter (Frieden 1975) which, as it is well-known, is completely corrupted by noise propagation as a consequence of the ill-posedness of the image restoration problem (see I for a general discussion). The same result holds true also in the case of multiple images (Piana & Bertero 1996), hence the need of introducing the so-called regularization methods (see again I for an introduction).