2 The LR/EM method for LBT imaging

We assume that p images of the same astronomical object, corresponding to p different orientations of the baseline, are available and that each image consists of $N\times N$ pixels, characterized by the indices m,n. We denote by f(m,n) and g_j(m,n) (m,n = 0,1,...,N-1; j = 1,2,...,p) the values, at the pixel m,n, of the brightness distribution of the unknown object and of its j-th image, respectively. Moreover we denote by $\vec{f}$ and $\vec{g}_j$ the arrays formed by these values.

In the case of space-invariant PSF's the mathematical model of image formation is given by

 

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

+ w_j(m,n)

where K_j(m,n) is the value, at the pixel m,n, of the PSF corresponding to the j-th orientation of the baseline and w_j(m,n) is the value of the noise (Poisson noise, read-out noise etc.) corrupting the j-th image. We will write this equation in the following synthetic form

$$\vec{g}_j =\vec{A}_j \vec{f} +\vec{w}_j$$ (2)

$\vec{A}_j$ being, in general, a block-Toeplitz matrix which, as usual, can be approximated by a block-circulant matrix in order to easily compute matrix-vector products by means of FFT. This approximation is sufficiently good in many circumstances, and precisely when the size of the main lobe of the PSF is much smaller than the size of the image.

When the approximation of a periodic PSF is used, by computing the FFT of both sides of Eq. (2) we obtain the usual relationship

$$\hat g_j(m,n)=\hat K_j(m,n)\hat f(m,n)+\hat w_j(m,n).$$ (3)

The FFT of the PSF K_j(m,n), i.e. $\hat {K}_j(m,n)$, is the transfer function of the system, corresponding to the j-th orientation of the baseline.

In order to formulate the LR/EM method for the restoration of LBT images in a compact form, let us define the product of two arrays $\vec g$ and $\vec h$, denoted by $\vec g \vec h$, as the array which is obtained by multiplying each component of $\vec g$by the corresponding component of $\vec h$, i.e. $(\vec {gh})(m,n) = g(m,n)h(m,n)$. Analogously the quotient $\vec{g}/\vec{h}$ is defined by $(\vec{g}/\vec{h})(m,n) = g(m,n)/h(m,n)$, provided that h(m,n) is different from zero everywhere.

LR/EM is an iterative algorithm for approximating the solutions of the maximum likelihood problem in the case where images are corrupted by Poisson noise (Shepp & Vardi 1982). In such a case, if we assume that the images $\vec{g}_j$ are statistically independent and that their values at distinct pixels are also statistically independent Poisson processes, then the log-likelihood function is given by (except for a constant term depending on the images $\vec{g}_j$)

 

$$l(\vec{f}) \sum_{j=1}^p\sum_{m,n=0}^{N-1}\{ g_j(m,n) \ln[(\vec{A}_j\vec{f})(m,n)]-$$ =

$$(\vec{A}_j\vec{f})(m,n) \}.$$ - (4)

By means of the arguments developed by Shepp & Vardi (1982) it is easy to derive the LR/EM algorithm in the present case. The parameters related to the normalization of the columns of the imaging matrix, in the case of block-circulant matrices, are given by

$$\eta_j=\sum_{m,n=0}^{N-1} K_j(m,n)=\hat K_j(0,0).$$ (5)

so that, if we put

$$\eta=\sum_{j=1}^{p} \hat K_j(0,0),$$ (6)

the LR/EM algorithm is as follows

$$\eta\vec{f}^{(k+1)}=\vec{f}^{(k)}\sum_{j=1}^{p} \vec{A}_j^T{{\vec{g}_j}\over{\vec{A}_j\vec{f}^{(k)}}} .$$ (7)

The iterations must be initialized with a positive image $\vec{f}^{(0)}$, the most simple choice being that of a uniform image.

It may be convenient to normalize the PSF's in such a way that $\hat K_j(0,0)=1$. Then we get $\eta=p$. As follows from Eq. (3), this normalization implies that all quantities $\hat g_j(0,0)$are equal, except for noise fluctuations, so that it can be consistently used if all images have essentially the same number of counts.

The usual properties of the EM algorithm can be easily extended to the present case: the result $\vec{f}^{(k)}$ of the k-th iteration is a non-negative estimate of the unknown object; in the case $\eta=p$ it satisfies the condition

$$\sum_{m,n=0}^{N-1} f^{(k)}(m,n)={1\over p} \sum^{p}_{j=1}\sum_{m,n=0}^{N-1} g_j(m,n),$$ (8)

i.e. the number of counts of the restored image is always equal to the arithmetic mean of the numbers of counts of the images $\vec{g}_j$; $l(\vec{f}^{(k)})$ is an increasing function of k, so that, thanks to an argument due to Vardi et al. (1985), the iterates $\vec{f}^{(k)}$ converge to a maximum point of $l(\vec{f})$.

However it is important to remark that, as it is well-known, it is not convenient to push the algorithm to full convergence: too many iterations generate artifacts, mainly due to noise amplification, which reduce the quality of the restoration. In other words the maximum points of $l(\vec{f})$ do not provide sensible solutions. These can be obtained by a suitable stopping of the iterations since, as it may be shown by means numerical simulations, the number of iterations acts as a regularization parameter, in the sense that the method as a property which is often denoted as semiconvergence: the RMS error representing the discrepancy between $\vec{f}^{(k)}$ and $\vec{f}$(the known object in the case of simulated images) first decreases for increasing k, remains approximately constant for a certain number of iterations and then grows again (Bertero & Boccacci 1998). An example of this behaviour will be given in Fig. 1.

This property implies that, for any given image, there exists an optimum number of iterations. In the case of real images it may be difficult to find this optimum number and this is the main shortcoming of the LR/EM method. This problem is even more important for the method we present in the next section because convergence (or semiconvergence) is faster. Indeed, in the case of LR/EM the minimum of the RMS error is, in general, rather flat and this property implies that, in the case of real images, one can stop the iterations when no significant modification is observed in the restoration. Another stopping rule may be provided by the so called discrepancy principle: the iterations are stopped when the RMS error representing the discrepancy between $\vec{A}_j\vec{f}^{(k)}$ and the observed image $\vec{g}_j$ becomes smaller than some estimated noise level.