A&A 444, 983-993 (2005)
DOI: 10.1051/0004-6361:20053480

Reconstruction of aperture-synthesis images from LBT LINC-NIRVANA data using the Richardson-Lucy and space-variant Building Block method

K.-H. Hofmann - T. Driebe - M. Heininger - D. Schertl - G. Weigelt

Max-Planck-Institut für Radioastronomie (MPIfR), Auf dem Hügel 69, 53121 Bonn, Germany

Received 20 May 2005 / Accepted 22 August 2005

Abstract
We present a new method, the regularized and space-variant Building Block method, which is able to reconstruct diffraction-limited aperture-synthesis images from Large Binocular Telescope (LBT) LINC-NIRVANA data. Images with the diffraction-limited resolution of a 22.8 m single-dish telescope can be derived if raw images are taken at several different hour angles. We simulated computer-generated and laboratory LBT interferograms that are similar to the data which can be obtained with the LINC-NIRVANA beam combiner instrument. From the simulated data, diffraction-limited images were reconstructed with the regularized Building Block method, which is an extension of the Building Block method (Hofmann & Weigelt 1993, A&A, 278, 328). We compare the Building Block reconstructions to images obtained with the Richardson-Lucy (RL) method (Richardson 1972, J. Opt. Soc. Am., 62, 55; Lucy 1974, AJ, 79, 745) and the Ordered Subsets Expectation Maximization (OSEM) method (Hudson & Larkin 1994, IEEE Trans. Med. Imag., 13, 601; Bertero & Boccacci 2000, A&AS, 144, 181). Our image reconstruction studies were performed with computer-simulated J-band and laboratory H-band raw data of a galaxy with simulated total magnitudes of $J=16^{\rm m}$ to 18$^{\rm m}$ and $H = 16^{\rm m}$ to 19$^{\rm m}$, respectively. One of the faintest structures in the images has a brightness of $J \sim 25$$^{\rm m}$. The simulated reference stars within the isoplanatic patch have magnitudes of $J = 20^{\rm m}{-}21^{\rm m}$ and $H = 19^{\rm m}$. All three methods are able to reconstruct diffraction-limited images with almost the same quality. Furthermore, raw data with space-variant point spread functions were simulated, and diffraction-limited images were reconstructed using the space-variant version of the Building Block method.

Key words: instrumentation: interferometers - instrumentation: high angular resolution - techniques: image processing - techniques: high angular resolution - techniques: miscellaneous

  
1 Introduction

The Large Binocular Telescope (LBT, Hill & Salinari 2003, 2004) with its two 8.4 m mirrors is currently built on Mt. Graham, Arizona. Due to the unique properties of the LBT facility, the raw images obtained with LINC-NIRVANA (Herbst et al. 2003, 2004), the optical/NIR beam combiner of the LBT, will have, in one direction, the spatial resolution equivalent to a diffraction-limited 8.4 m single-dish telescope, and in the orthogonal direction, that of a 22.8 m telescope. If exposures are taken at different hour angles in order to obtain a better uv-coverage, an aperture synthesis image reconstruction procedure can provide reconstructed images with a spatial resolution which reaches that of a 22.8 m single-dish telescope.

The reconstruction of high-resolution images from LBT raw data has already been discussed by several authors. First wide-field aperture synthesis studies were presented by Hege et al. (1995) and Angel et al. (1998). Reinheimer et al. (1997) have studied the reconstruction of diffraction-limited images from computer-simulated LBT speckle interferograms. The images were reconstructed with the Bispectrum Speckle Interferometry method and the Building Block method. Correia & Richichi (2000), Bertero & Boccacci (2000), Carbillet et al. (2002), Correia et al. (2002), and Anconelli et al. (2005) have investigated the deconvolution of Adaptive Optics (AO) corrected LBT raw images with the Richardson-Lucy algorithm and related methods.

In this paper, we will present an extension of the Building Block method and apply this new method to LINC-NIRVANA raw data obtained from computer and laboratory simulations. The main features of the extension are a regularization scheme to improve the quality of the reconstructed images and a method to reconstruct images in the case of space-variant point-spread functions (PSFs) in the large LINC-NIRVANA field-of-view (FOV) of $\sim$10 $\hbox{$^{\prime\prime}$ }$. Space-variant deconvolution is important for the optimal reconstruction of images obtained with the LINC-NIRVANA instrument.

The paper is organized as follows: in Sect. 2 we describe the Building Block image reconstruction algorithm (BB) in more detail than what was previously presented by Hofmann & Weigelt (1993). In addition, an extension of the BB method including regularization is discussed. In Sect. 3 we investigate the reconstruction of diffraction-limited images from computer-simultated LINC-NIRVANA LBT raw images (with partial AO corrections) with different image reconstruction algorithms: the Richardson-Lucy (RL) iterative algorithm (Richardson 1972; Lucy 1974), also known as the Expectation Maximation method (EM), the Ordered Subsets Expectation Maximization (OSEM) method (Hudson & Larkin 1994; Bertero & Boccacci 2000), which is an extension of the RL method for cases with more than one input image, and the Building Block (BB) method with regularization. Furthermore, in Sect. 4 we present the reconstruction of images from laboratory simulations of LINC-NIRVANA LBT raw images (called LBT raw images below) taken with a HAWAII array camera. Finally, in Sect. 5 we investigate the influence of space-variant PSFs on the quality of the reconstructed images. We present a modified version of the BB method that is able to handle space-variant PSFs.

  
2 The building block method

The basic principal of the Building Block method has been briefly outlined by Hofmann & Weigelt (1993). In this section, we discuss this method in more detail and present an extension of the BB method (including regularization) which we used for the image reconstruction studies described in Sects. 3 and 4.

The intensity distribution of an object o(x) can be described as a sum of many building blocks (e.g. $\delta$-functions). The goal of the BB method is the calculation of a high-resolution aperture synthesis image that is equal to the convolution of o(x) with the diffraction-limited point-spread function (PSF) of a large single-dish telescope; for instance, a 22.8 m aperture in the case of the LBT. The BB method iteratively produces images o_k(x) (k=1, 2,...) by adding one or more building blocks $\delta(x-x_0)$ at each iteration step k. The goal is to find an image o_k(x) which minimizes the $\chi^2$ function

$$Q[o_k(x)] ~:=~ \int ~\left\vert \frac{o_k(x)\otimes p(x) - i(x)}{\sigma(x)} \right\vert^2 {\rm d}x,$$ (1)

where x is a 2-dimensional space coordinate, $\otimes$ denotes the convolution operator, p(x) is the sum LBT PSF discussed in Sect. 3.1 (i.e., the sum of the PSFs taken at different hour angles), and $\sigma(x)$ denotes the Poisson and detector noise in the image i(x) which is the sky background-subtracted sum LBT raw image (i.e., the sum of the individual raw images taken at different hour angles; see Sect. 3.1). The functions i(x), p(x), and o_k(x) are normalized to an integral intensity of 1.

The new building block added to o_k(x) to obtain o_k+1(x)is positioned at a particular coordinate x=x₀ in o_k(x). Since o_k(x) and o_k+1(x) are normalized to an integral intensity of 1, o_k+1(x) is given by

$$o_{k+1}(x) = \frac{n_k \cdot o_k(x) + \delta(x-x_0)}{n_k+1} \approx o_k(x) + \frac{1}{n_k}\cdot \delta(x-x_0),$$ (2)

where $\delta(x-x_0)$ is the new building block, and n_k denotes the number of building blocks of o_k(x). Since $n_k\gg1$ after a large number k of iteration steps, the right side of Eq. (2) holds. In this realization of the BB method all building blocks have the same amplitude. Another approach is to fit simultaneously the amplitude and position of the new building block as mentioned in Hofmann & Weigelt (1993). For the position x₀ of the new building block, a value is chosen that leads to an image which minimizes the $\chi^2$ function $Q[o_k(x) + 1/n_k \cdot \delta(x-x_0)]$. For a large number k of iterations, each contribution $1/n_k \cdot \delta(x-x_0)$ is very small. Therefore, $Q[o_k(x) + 1/n_k \cdot \delta(x-x_0)]$ can be approximated by a Taylor expansion:

 

$${Q[o_k(x) + 1/n_k \cdot \delta(x-x_0)] \approx Q[o_k(x)] + \frac{1}{n_k} \cdot \frac{\partial Q}{\partial o_k(x_0)} }$$

$$= Q[o_k(x)] + \frac{2}{n_k} \cdot \frac{o_k(x_0)\otimes p(x_0) - i(x_0)}{\sigma^2(x_0)} \otimes p(-x_0).$$ (3)

Equation (3) shows that the $\chi^2$ function Q is minimized by adding a positive (negative) building block at that position in o_k(x) where the function $\partial Q/\partial o_k(x_0)$ has its absolute minimum (maximum). Reconstructions are often improved by simultaneously adding positive and negative building blocks per iteration step. Negative building blocks are only added, if the positivity constraint is not violated. Adding more than one positive and one negative building block per iteration step, for instance by clipping the $\chi^2$ function in Eq. (3), accelerates the iteration process and can yield better reconstructions. This can be realized, for example, by adding positive (negative) building blocks at the positions where the values of $\partial Q/\partial o_k(x_0)$ are smaller (larger) than a certain fraction of the value of its absolute minimum (maximum).

One well-known technique that can lead to improved reconstructions is the regularization of the fitting procedure: instead of searching for the global minimum of the $\chi^2$ function Q, the minimum of Q must be estimated under certain conditions. One such condition is the smoothness of the reconstructed image, expressed, for example, by the entropy function (Wahl 1984)

$$H[o(x)] := \int o(x) \cdot \ln[o(x)] ~ {\rm d}x,$$ (4)

where $\ln$ is the natural logarithm. The goal is to find the image that minimizes the function

$$J[o(x)] := Q[o(x)] + \mu \cdot H[o(x)],$$ (5)

where Q[o(x)] is the $\chi^2$ function given in Eq. (1), and $\mu$ is a Lagrange multiplier which can be roughly estimated as described in Appendix A. To find the image that minimizes J[o(x)], the Taylor expansion of $J[o_k(x) + 1/n_k \cdot \delta(x-x_0)]$around o_k(x) is required:

 

$${J[o_k(x) + 1/n_k \cdot \delta(x-x_0)] \approx J[o_k(x)] + \frac{1}{n_k} \cdot \frac{\partial J}{\partial o_k(x_0)} }$$

$$= J[o_k(x)] + \frac{1}{n_k} \cdot \left\{\frac{\partial Q}{\partial o_k(x_0)} + \mu \cdot \frac{\partial H}{\partial o_k(x_0)}\right\},$$ (6)

where $\partial H/\partial o_k(x_0)$, according to Eq. (4), is given by

$$\frac{\partial H}{\partial o_k(x_0)} = \ln[o_k(x_0)] + 1.$$ (7)

Equation (6) shows that the function J is minimized by adding a positive (negative) building block at the position in o_k(x) where the function $\partial J/\partial o_k(x_0)$ has its absolute minimum (maximum). Applications of this BB method with regularization are presented in the next two sections.

  
3 Reconstruction of aperture-synthesis images from computer-simulated LBT data

In this section we investigate the reconstruction of diffraction-limited images from computer-simulated LBT raw images using the RL, OSEM, and BB methods.

  
3.1 Raw data of the computer simulations

The main goal is the reconstruction of diffraction-limited images from several LBT raw images recorded at different hour angles or pupil position angles. In other words, the goal is to obtain high-resolution images through aperture synthesis. The LBT raw images i'_j(x) are modeled as

$$i'_j(x) = o(x)\otimes p_j(x) + b_j(x) + r_j(x),$$ (8)

where o(x) denotes the intensity distribution of the object; i'_j(x) and p_j(x) with j=1,..,m are the m LBT raw images and PSFs, respectively, recorded at m different pupil position angles; b_j(x) describes the sky background, and r_j(x) is the noise caused by the detector (mainly read-out noise) as well as the limited number of photons of the object and the sky background. The subtraction of the sky background yields the input images $i_j(x) = o(x)\otimes p_j(x) + r_j(x)$. The OSEM method uses the individual images i'_j(x) (including sky background) and the sky background-subtracted PSFs p_j(x)(without clipping negative values) as input. The RL method uses the sum i'(x) of all LBT raw images i'_j(x) (with sky background) and the sky background-subtracted sum PSF p(x) of the corresponding LBT raw PSFs p_j(x) as input data. The Building Block method uses the sum i(x) of all sky background-subtracted LBT raw images i_j(x), the sum p(x) of the corresponding sky background-subtracted LBT raw PSFs p_j(x), and the sum b(x) of all sky backgrounds b_j(x) as input data.

Figure 1: Reconstruction of aperture-synthesis images from computer-simulated LBT raw data (all images are shown on the same scale). a) test object (based on the HST image of NGC 3504); b) object from panel a) convolved with the PSF of a simulated diffraction-limited 22.8 m telescope; c), d) computer-simulated LBT raw images of the object (simulated total magnitude $J=17.0^{\rm m}$, sky background $J=16.0^{\rm m}$ per arcsec square, read-out noise of 11 e- rms, pupil position angles 0$^\circ$ and 30$^\circ$); e), f) LBT PSFs for position angles 0$^\circ$ and 30$^\circ$ These PSFs are derived from four unresolved stars in the vicinity of the target (see text and Fig. 2). The brightness of these reference stars ranges from $J=20.0^{\rm m}$ to $21.0^{\rm m}$; g), h), i) diffraction-limited images reconstructed from 6 LBT raw images (taken with pupil position angles of 0$^\circ$, 30$^\circ$, 60$^\circ$, 90$^\circ$, 120$^\circ$,and 150$^\circ$, two of which are shown in panels c) and d) and the corresponding 6 LBT PSFs (two of which are shown in panels  e) and f) using the RL g), the OSEM h), and the BB i) methods; the faint star-like structure in the upper left corner of panels a and b has a brightness of $J=25.2^{\rm m}$ (marked by an arrow in panel b). The restoration errors are 5.37% g), 5.61% h), and 5.20% i) (see text for details). The typical computing (PC) times for the RL, OSEM, and BB methods were 6 min, 3 min, and 12 min, respectively.

Open with DEXTER

Figure 2: a) One of six computer-simulated LBT raw images with the field-of-view (FOV) of the LINC-NIRVANA instrument of 10 arcsec corresponding to $2048\times 2048$ pixel (total magnitude of the galaxy $J=17^{\rm m}$, sky background $J=16^{\rm m}$ per arcsec square, J magnitudes of the four reference stars: 20.0 $^{\rm m}-21.0^{\rm m}$, distances of the four reference stars to the galaxy: 2.6 $\hbox {$^{\prime \prime }$ }{-} 5.8\hbox {$^{\prime \prime }$ }$; see Table 1 and Fig. 1). b) One of six LBT PSFs (pupil position angle 30$^\circ$, Strehl ratio $\sim$30%) derived from the computer-simulated LBT raw image a) by adding up the PSFs of the four re-centered reference stars (256 $\times$ 256 pixel FOV). c) One of six LBT raw images of the target used for the deconvolution process.

Open with DEXTER

Figure 1 presents the results of our image reconstruction experiments using the RL, OSEM, and BB methods. In these computer experiments, a computer-generated galaxy with a total J-band magnitude of 17$^{\rm m}$ was chosen. The faint, star-like structure in the upper left corner of Fig. 1b (marked by an arrow) has a simulated J magnitude of 25.2$^{\rm m}$. The magnitudes of the reference stars in the simulation are 20.0$^{\rm m}$ (top left in Fig. 2a), 20.5$^{\rm m}$ (top right), 20.5$^{\rm m}$ (bottom right), and 21.0$^{\rm m}$ (bottom left), respectively.

Figure 1 shows (a) the computer-generated reference object, (b) the reference object convolved with a theoretical PSF corresponding to a single-dish 22.8 m pupil, (c, d) two of the six raw images simulated for two different pupil position angles (0$^\circ$ and 30$^\circ$), and (g-i) the results of the image reconstruction experiments obtained with the three different deconvolution methods.

The raw images are obtained through convolution of the object intensity distribution (Fig. 1a) with the simulated LBT PSFs (for example, Figs. 1e, f), addition of sky background, and simulation of photon and detector read-out noise. LBT PSFs with partial AO corrections (Strehl ratio $\sim$30%) were simulated in the following way:

(1)

a phase screen was generated which produces speckle interferograms with one dominant diffraction-limited core surrounded by many faint speckles, which correspond to the speckle halo caused by partial AO corrections;

(2)

an intermediate short-exposure LBT PSF was obtained by calculating the squared Fourier transform of the complex amplitude transmission of the phase screen multiplied with the LBT pupil function;

(3)

the summation of a large number of intermediate short-exposure LBT PSFs, derived from many different phase screens, yielded the LBT PSF for the corresponding position angle and exposure time (see Figs. 1e, f).

For the computer experiments shown in Fig. 1, we simulated PSFs and raw images of the object corresponding to six different pupil position angles equally spaced between 0$^\circ$ and 180$^\circ$. We assumed an integration time of 600 s per position angle, a 30$\%$ efficiency (optics and detector), a read-out noise of 11 electrons rms, and a J-band sky background brightness of 16$^{\rm m}$ per arcsec square. The simulated photon noise of the six LBT raw images of the galaxy (two are shown Figs. 1c and d; and one in Fig. 2c) corresponds to a total magnitude of the whole galaxy of $J=17.0^{\rm m}$. The 6 LBT PSFs (two are displayed in Figs. 1e and f; and one is shown in Fig. 2b) were derived from four unresolved stars located close to the target and with distances to the center of the target ranging between 2.6 $\hbox{$^{\prime\prime}$ }$ and 5.8 $\hbox{$^{\prime\prime}$ }$. Figure 2a shows one of the six 2048$\times$2048 pixel computer-simulated LBT raw images with the galaxy and the four surrounding unresolved stars. The parameters of the computer simulation described here are listed in Table 1.

   
3.2 Image reconstruction results obtained using the RL, OSEM, and BB methods

The LBT raw images and LBT PSFs were preprocessed in the following way:

(1) The PSFs for the deconvolution process were obtained by adding up the PSFs of the four re-centered, unresolved reference stars surrounding the target (see Fig. 2).

(2) The sky background was estimated and subtracted from the LBT PSFs (without clipping negative values); for the BB method, the sky background was subtracted from the LBT raw images of the target (without clipping negative values);

(3) In the Fourier transform of each individual LBT raw image, all elements outside the uv-coverage corresponding to a given pupil function were set to zero in order to reduce the noise. For illustration, Fig. 3a shows the modulus of the Fourier transform of a single LBT raw image of a point source with photon noise. In the case of an ideal, i.e. noise-free PSF, the area with Fourier modulus values larger than zero is the instantaneous uv-coverage. In Fig. 3a the noise outside the uv-coverage can be seen. Figure 3b shows the Fourier transform of the same raw image, but with all Fourier elements outside the uv-coverage set to zero. In Appendix B we will show that diffraction-limited images reconstructed from Fourier-masked LBT raw data have lower restoration errors (see definition in Eq. (9)) than reconstructions obtained without Fourier-masking.

(4) Co-adding of the re-centered LBT raw images recorded at all six pupil position angles and preprocessed as described in (2) and (3) yields the input data for the RL and BB methods. The LBT raw images were re-centered using the shift-vectors determined from the re-centering procedure of the surrounding reference stars for the derivation of the LBT-PSFs (1). OSEM uses the six individual preprocessed LBT raw images and LBT PSFs as input data.

The reconstruction methods RL, OSEM, and BB iteratively reconstruct diffraction-limited images from the LBT raw data. For the RL and OSEM deconvolution experiments presented in this paper, the publicly available AIRY package (version 2.0) was used (Correia et al. 2002) with the estimated sky background level as an additional input value. The BB image reconstruction experiments presented in this paper were performed with the algorithm and the regularization scheme outlined in Sect. 2.

Table 1: Parameters of the computer-simulated LBT raw images (Fig. 1).

Figures 1g-i show the images derived from six LBT raw images (see examples in Figs. 1c and d) and the corresponding six LBT PSFs (two are shown in Figs. 1e and f) by using the three image reconstruction methods RL, OSEM, and BB, respectively. In our simulations, the restoration error, which is a measure of the quality of the image reconstruction, is defined as

$$\rho := \frac{\sqrt{\int \vert o_k(x)\otimes p'(x) - o(x)\oti... ...rm d}x}}{\sqrt{\int \vert o(x)\otimes p'(x)\vert^2 {\rm d}x}},$$ (9)

where $o(x)\otimes p'(x)$ is the computer object convolved with the theoretical diffraction-limited PSF p'(x) of a hypothetical 22.8 m single-dish telescope (Fig. 1b), and $o_k(x)\otimes p'(x)$ (Figs. 1g-i) is the reconstruction convolved with the same PSF p'(x). Both images are normalized according to

$\int o_k(x)\otimes p'(x) {\rm d}x = \int o(x)\otimes p'(x) {\rm d}x$. In the definition of $\rho$, the reconstruction o_k(x) and the reference object o(x) are convolved with p'(x)since the goal of interferometric LBT imaging is usually to obtain an image with the diffraction-limited resolution of a hypothetical 22.8 m single-dish telescope.

The average restoration errors of the three reconstructions were determined to be $5.41\pm0.10$% (RL), $5.58\pm0.08$% (OSEM), and $5.30\pm0.10$% (BB). These errors were calculated from 10 raw data sets (see parameters in Table 1) that were statistically independent with respect to photon and detector read-out noise and the phase screens simulating partial AO. Figure 4 shows intensity cuts along a line through the nucleus and the faint star-like object (upper left corner in Fig. 1b) for the test object (Fig. 1b; red), the RL reconstruction (Fig. 1g; green), the OSEM reconstruction (Fig. 1h; blue), and the BB reconstruction (Fig. 1i; pink).

In order to study the dependence of the quality of the reconstructions on the brightness of the astronomical target, we performed two additional experiments with the computer-generated galaxy shown in Fig. 1a and total J-band magnitudes of 16$^{\rm m}$ and 18$^{\rm m}$. All other parameters of the computer-simulated LBT raw images were the same as for the experiment shown in Fig. 1 (see parameters in Table 1). The average restoration errors of the reconstructions obtained with the RL, OSEM, and BB methods are listed in Table 2. As in the case of object brightness $J=17^{\rm m}$, all three methods yield reconstructions of similar quality and restoration errors.

	Figure 3: a) Modulus of the Fourier transform of a single computer-simulated LBT raw image of a point source (pupil position angle 30$\circ$). b) Same as a), but with all Fourier elements outside the uv-coverage set to zero (see text).
Open with DEXTER

	Figure 4: Intensity cuts along a line through the nucleus and the faint star-like object (upper left corner in Fig. 1b) for the test object (Fig. 1b; red), the RL reconstruction (Fig. 1g; green), the OSEM reconstruction (Fig. 1h; blue), and the BB reconstruction (Fig. 1i; pink).
Open with DEXTER

   
4 Image reconstruction from LINC-NIRVANA data obtained with a laboratory LBT simulator

4.1 Laboratory simulated LINC-NIRVANA raw data

Simulated LBT LINC-NIRVANA data of a galaxy (see Fig. 5a) were recorded with a laboratory LBT simulator and a HAWAII-1 array camera. The simulator consists of the lenses of an imaging setup, a slide which contains the laboratory object and the reference star, a rotatable pupil mask which simulates the LBT pupil, and a rotating phase screen that allows the simulation of partial AO correction.

Table 2: Restoration errors of the reconstructions using the galaxy test object shown in Fig. 1 (see also Table 1) with total magnitudes of $J=16.0^{\rm m}$, $J=17.0^{\rm m}$, and $J=18.0^{\rm m}$. The restoration errors listed are averages of 10 statistically independent LBT data sets (see text).

Figure 5: Laboratory simulation of LINC-NIRVANA and image reconstruction from laboratory data using the RL, OSEM, and BB methods. a) Test object (galaxy) used for the laboratory experiment (simulated resolution corresponding to a diffraction-limited 22.8 m single-dish telescope). b), c) Laboratory LBT raw images and PSFs for pupil position angles 0$^\circ$ and 60$^\circ$. d), e), f) Diffraction-limited reconstructions derived from 6 LBT raw images (total magnitude $H = 16.2^{\rm m}$) using the RL d), OSEM e), and BB methods f). The restoration errors are 6.85% d), 6.76% e), and 6.23% f). The typical computing (PC) times for the RL, OSEM, and BB methods were 6 min, 4 min, and 3 min, respectively.

Open with DEXTER

Table 3: Parameters of our laboratory simulation.

Table 4: Restoration errors of the reconstructions obtained with the RL, OSEM, and BB methods from the Fourier-masked laboratory LBT raw images.

The laboratory data was taken in the H band. Images of the object and one unresolved star were simultaneously recorded (see Figs. 5b and c). For each image reconstruction experiment, images were taken with 6 different pupil position angles (see Table 3). Data sets with simulated total galaxy magnitudes of $H = 16.2^{\rm m}$, 17.1$^{\rm m}$, 17.8$^{\rm m}$, and 19.1$^{\rm m}$ were recorded. For each of the 4 object brightnesses, 10 statistically independent data sets were recorded, which allowed the determination of the error of the restoration errors. Sky background was simulated as thermal background by heating the LBT pupil mask. The object and the PSF were recorded with a simulated full 22.8 m pupil to be able to compare the reconstructions of the LINC-NIRVANA simulations with full-pupil images. A summary of the parameters of the laboratory experiments is given in Table 3.

The simulated PSFs consist of a dominant diffraction-limited 3-fringe core and an extended AO speckle halo (see Figs. 5b and c). The speckle contrast in this halo is high if the exposure times are short compared to the simulated speckle life time and low (as in our experiments) if the exposure times are long compared to the speckle life time. The Strehl ratio simulated in this experiment was 0.10.

The raw laboratory images were preprocessed in the following way:

(1) The images were flat-fielded, and pixel bias and bad pixel effects were corrected.

(2) The sky background was estimated by fitting the measured sky background to the actual background in the data.

(3) To reduce the noise in the individual data corresponding to a certain hour angle, all elements in the Fourier transform outside the uv-coverage of the simulated LBT pupil were set to zero. As discussed in Appendix B, this Fourier-masking reduces the restoration error.

(4) For the RL and BB methods, all 6 images corresponding to the different hour angles were co-added. For OSEM, the individual images were used.

   
4.2 Image reconstruction results obtained using the RL, OSEM, and BB methods

Figure 5 summarizes the results of one of the laboratory LINC-NIRVANA image reconstruction experiments. In this experiment the following parameters were simulated: total magnitude of the galaxy $H = 16.2^{\rm m}$, magnitude of the reference star $H = 19.1^{\rm m}$, and sky background $H = 14.0^{\rm m}$ per square arcsec. Figure 5a shows an image of the laboratory object taken with a simulated 22.8 m single dish pupil. Figures 5b and c show two of the six simulated raw images (pupil position angle = 0$^\circ$ and 60$^\circ$). Figures. 5d-f show the images reconstructed with the RL, the OSEM, and the BB methods. The restoration errors obtained for the four different object magnitudes $H = 16.2^{\rm m}$, 17.1$^{\rm m}$, 17.8$^{\rm m}$, and 19.1$^{\rm m}$ are summarized in Table 4.

	Figure 6: Intensity cuts along a horizontal line through the nucleus for the test object (Fig. 5a; red), the RL reconstruction (Fig. 5d; green), the OSEM reconstruction (Fig. 5e; blue), and the BB reconstruction (Fig. 5f; pink).
Open with DEXTER

As in the case of the computer simulations discussed in the previous section, the reconstructed images obtained with the three methods RL, OSEM, and BB are of similar quality. For a given object brightness all restoration errors agree within the error bars. As expected, in general the restoration quality becomes worse with decreasing object brightness (from $\sim$5% for $H = 16^{\rm m}$ to $\sim$13% for $H = 19^{\rm m}$). Figure 6 shows intensity cuts along a horizontal line through the nucleus for the test object (Fig. 5a; red), the RL reconstruction (Fig. 5d; green), the OSEM reconstruction (Fig. 5e; blue), and the BB reconstruction (Fig. 5f; pink).

  
5 Space-variant Building Block method: image reconstruction from LBT raw data in the case of space-variant PSFs

Within the large field-of-view of the LINC-NIRVANA instrument, the PSF will be space-variant. In this section we investigate the influence of space-variant PSFs on the quality of reconstructions obtained (a) with the conventional BB method described above, and (b) with a modified version of the BB method that is able to deal with space-variant PSFs.

5.1 Modification of the BB method for space-variant PSFs

If the PSF is space-variant, the convolution operations in the approximation of the $\chi^2$ function Q (see Eq. (3)) have to be replaced by space-variant convolution operations. The space-variant version of Eq. (3) is given by (see Appendix C)

 

$${Q[o_k(x) + 1/n_k \cdot \delta(x-x_0)] \approx ~~ Q[o_k(x)] }$$

$$\quad+ \frac{2}{n_k} \cdot \left\{\int o_k(\xi)\cdot f(x_0;\xi)~{\rm d}\xi - \int \frac{i(\xi)}{\sigma^2(\xi)}\cdot p(-x_0;\xi)~{\rm d}\xi\right\}$$

$$\qquad= Q[o_k(x)] + \frac{2}{n_k} \cdot A_k(x_0) - \frac{2}{n_k} \cdot B(x_0),$$ (10)

with the space-variant function

$$f(x_0;\xi) := \int \frac{p(\xi';\xi)}{\sigma^2(\xi')} \cdot p(-x_0;\xi')~{\rm d}\xi'.$$ (11)

For the determination of the position of the building block to be set in the next iteration step, we need only the second and third term in Eq. (10). The third term B(x₀) consists of the space-variant convolution of the measured LBT raw image i(x) with the space-variant LBT PSF p(x₀;x), weighted with the inverse of $\sigma^2(x)$. This time-consuming operation is independent of the iteration process and, therefore, can be carried out before starting the iteration process. The second term A_k(x₀) in Eq. (10) is the space-variant convolution of the iterated image o_k(x) at iteration step k with the space-variant function f(x₀;x)defined in Eq. (11). This term must be calculated at each iteration step. In order to save time, A_k(x) is, like o_k(x), built up iteratively by adding the space-variant function f(x₀;x) at position x₀ of the next building block.

Figure 7: Image reconstruction from computer-simulated LBT data with space-variant PSFs using the space-variant BB method. a) Computer-generated object convolved with the PSF corresponding to a single-dish 22.8 m pupil: the object is a cluster of 7 stars (5 stars have magnitude J = 22.21$^{\rm m}$, the 2 fainter stars have J = 22.97$^{\rm m}$). b) One of the simulated LBT raw images (position angle 60$^\circ$): Poisson noise corresponding to an average number of $8\times 10^5$ photons, detector noise of 11 e- rms, and a sky background brightness of J = 16$^{\rm m}$ per arcsec square. The three reference stars are located at the top left, bottom left and at the right-hand side; c)-e) PSFs of the 3 reference stars magnified by a factor of $\sim$ 2.5 and displayed without photon noise to visualize their space-variant Strehl ratios of 92$\%$, 60$\%$, and 17$\%$ (from left to right); f), g) diffraction-limited images reconstructed from 3 LBT raw images recorded with pupil position angles 0$^\circ$, 60$^\circ$, and 120$^\circ$ obtained with the space-variant f) and space-invariant g) version of the Building Block method. The restoration errors of the two reconstructions are 5.55% f); space-variant BB method) and 21.67% g); space-invariant BB method). The typical computing (PC) times for the space-variant and space-invariant BB methods were 5 min and 5 min, respectively.

Open with DEXTER

5.2 Image reconstruction from computer-simulated LBT raw data with space-variant PSFs using the space-variant BB method

Figure 7 shows a computer experiment of image reconstruction with the space-variant BB method. The target chosen for this simulation is a star cluster which consists of 7 stars and is surrounded by 3 unresolved reference stars (see Figs. 7a, b). Five of the target stars have identical magnitudes (J = 22.21$^{\rm m}$) and two stars are $\sim$2 times fainter (J = 22.97$^{\rm m}$). The 3 additional unresolved stars located in the vicinity of the target (see Fig. 7b) serve as reference stars. They all have J = 22.21$^{\rm m}$. In Fig. 7b one of three simulated LBT raw images with space-variant PSF is shown. Each simulated space-variant PSF is the weighted sum of a diffraction-limited LBT PSF (i.e. a fringed Airy disk) plus a halo resulting from simulated partial AO correction. The space-variant weight w(x) of the fringed Airy disk and the weight 1-w(x) of the halo were assumed to be a linear function of the two cartesian coordinates of the 2-dimensional space vector x. The space-variant PSFs of the three reference stars shown in Figs. 7c, d, e (magnified by a factor of 2.5) have Strehl ratios of 92$\%$, 60$\%$, and 17$\%$, respectively. Poisson noise of the target corresponding to $8\times 10^5$ photons, a sky background brightness of J = 16$^{\rm m}$, and detector read-out noise of 11 electrons rms were simulated. Figures 7f and g show BB reconstructions obtained from 3 simulated LBT raw images (position angles 0$^\circ$, 60$^\circ$, and 120$^\circ$) with space-variant PSFs. For comparison, the reconstructions were obtained with both a space-variant and a space-invariant version of the BB method discussed below.

In the space-invariant version of the BB method applied in this computer experiment (Fig. 7), the $\chi^2$ function discussed in Sect. 2 was approximated by

$$Q[o_k(x) + 1/n_k \cdot \delta(x-x_0)] \approx ~~ Q[o_k(x)] + \frac{2}{n_k} \cdot \{o_k(x_0)\otimes p(x_0) - i(x_0)\}.$$ (12)

The space-variant version of the BB method, used in this experiment (Fig. 7), is an extension of Eq. (12) for the case of space-variant PSFs. The approximation of the $\chi^2$ function in the space-variant version is given by

 

$$Q[o_k(x) + 1/n_k \cdot \delta(x-x_0)] \approx Q[o_k(x)] + \frac{2}{n_k} \cdot \left\{\int o_k(\xi)\cdot p(x_0;\xi)~d\xi - i(x_0)\right\}$$

$$Q[o_k(x)] + \frac{2}{n_k} \cdot \{A(x_0) - B(x_0)\}.$$ = (13)

For the space-invariant BB method, the required input PSF was the sum of the PSFs of all reference stars. For the space-variant BB method, the PSF at the positions of the individual new building blocks were derived from the PSFs of the 3 reference stars through linear interpolation.

In Table 5 the theoretical and reconstructed magnitudes of the seven stars are summarized. The comparison shows that the space-variant BB method yields reconstructions with roughly 4 times smaller magnitude errors than the conventional space-invariant BB method. The restoration errors (see definition in Eq. (9)) were determined to be 5.55% for the space-variant BB reconstruction (Fig. 7f), and 21.67% for the space-invariant BB reconstruction (Fig. 7g).

  
6 Conclusions

In this paper we presented an extended version of the Building Block algorithm for the reconstruction of diffraction-limited images from data obtained with the LBT LINC-NIRVANA instrument. The main features of this extension, a regularization scheme and a method to handle space-variant PSFs, are outlined. We carried out a large number of image deconvolution experiments to study the quality of images reconstructed from simulated LINC-NIRVANA data with the regularized Building Block method. We compared the Building Block reconstructions to images obtained with the Richardson-Lucy method and the Ordered Subsets Expectation Maximization method, and determined the restoration error as a function of the object's brightness. We find that all three methods are able to reconstruct diffraction-limited images with almost the same quality. Our image reconstruction studies were performed with computer-simulated and laboratory raw data of a galaxy with simulated total magnitudes of $J=16^{\rm m}$ to 18$^{\rm m}$, and $H = 16^{\rm m}$ to 19$^{\rm m}$, respectively. One of the faintest structures in the images has a brightness of $J \sim 25$$^{\rm m}$. The simulated reference stars have magnitudes in the ranges of $J = 20^{\rm m}$ to 21$^{\rm m}$ and $H = 19^{\rm m}$.

Furthermore, simulated LINC-NIRVANA data of a compact star cluster with space-variant point spread functions in the FOV (caused by partial AO correction) were simulated, and diffraction-limited images were reconstructed using the space-variant version of the Building Block method. This space-variant Building Block method allowed us to obtain much higher photometric accuracy than with the space-invariant method. Therefore, the simulation illustrates the importance of space-variant deconvolution methods.

Table 5: Computer simulation with space-variant PSFs (Fig. 7): magnitudes of the stars reconstructed with the space-variant BB (sv. BB) and the conventional BB (conv. BB) methods, and the differences ($\Delta$) between the reconstructed and original star magnitudes. The values show that the deviations from the true magnitudes are on average $\sim$4 times smaller in the reconstruction obtained with the space-variant BB method.

Acknowledgements

We would like to thank the referee for the various valuable comments which helped to widely improve the manuscript.

References

• Anconelli, B., Bertero, M., Boccacci, P., Carbillet, M., & Lanteri, H. 2005, A&A, 430, 731 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Angel, J. R. P., Strittmatter, P. A., Salinari, P., & Weigelt, G. 1998, Astronomical Interferometry, ed. R. D. Reasenberg, Proc. SPIE, 3350, 881 (In the text)
• Bertero, M., & Boccacci, P. 2000, A&AS, 144, 181 [EDP Sciences] [NASA ADS] (In the text)
• Carbillet, M., Correia, S., Boccacci, P., & Bertero, M. 2002, A&A, 387, 744 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Correia, S., & Richichi, A. 2000, A&AS 141, 301 (In the text)
• Correia, S., Carbillet, M., Boccacci, P., Bertero, M., & Fini, L. 2002, A&A, 387, 733 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Hege, E. K., Angel, J. R. P., Cheselka, M., & Lloyd-Hart, M. 1995, Optical Aperture Synthesis, ed. J. D. Gonglewski, Proc. SPIE, 2566B, 144 (In the text)
• Herbst, T., Raggazzoni, R., Andersen, D., et al. 2003, Interferometry for Optical Astronomy II, ed. W. A. Traub, Proc. SPIE, 4838, 456 (In the text)
• Herbst, T. M., Ragazzoni, R., Eckart, A., & Weigelt, G. 2004, Ground-based Instrumentation for Astronomy, ed. A. F. Moorwood, & M. Iye, Proc. SPIE, 5492, 1045 (In the text)
• Hill, J. M., & Salinari, P. 2003, Large Ground-based Telescopes, ed. J. M. Oschmann, & L. M. Stepp, Proc. SPIE, 4837, 140 (In the text)
• Hill, J. M., & Salinari, P. 2004, Ground-based Telescopes, ed. J. M. Oschmann, Jr., Proc. SPIE, 5489, 603 (In the text)
• Hofmann, K.-H., & Weigelt, G. 1993, A&A, 278, 328 [NASA ADS] (In the text)
• Hudson, H. M., & Larkin, R. S. 1994, IEEE Trans. Med. Imag., 13, 601 (In the text)
• Lucy, L. B. 1974, AJ, 79, 745 [NASA ADS] [CrossRef] (In the text)
• Reinheimer, T., Hofmann, K.-H., Schöller, M., & Weigelt, G. 1997, A&AS, 121, 191 [EDP Sciences] [NASA ADS] (In the text)
• Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55 [NASA ADS] (In the text)
• Wahl, F. M. 1984, Digitale Bildverarbeitung (Heidelberg-New-York-Tokio: Springer Verlag) (In the text)

 

  
Online Material

  
Appendix A: The Lagrange multiplier $\mu$ in the regularized Building Block method

In this appendix we derive an expression for the Lagrange multiplier $\mu$ in J[o(x)] (Eq. (5)).

The best-fitting image (the desired high-resolution reconstruction) is obtained when the distance function J (see Eq. (5)) has reached its minimum value. Therefore, the derivative $\partial J/\partial o(x_0)$vanishes:

$$\frac{\partial J}{\partial o(x_0)} = \frac{\partial Q}{\partial o(x_0)} + \mu \cdot \frac{\partial H}{\partial o(x_0)} = 0.$$ (A.1)

Thus, the integral of $\partial J/\partial o(x_0)$ multiplied with o(x₀)must also vanish:

$$\int \frac{\partial Q}{\partial o(x_0)} \cdot o(x_0) \, {\rm ... ...c{\partial H}{\partial o(x_0)} \cdot o(x_0) \, {\rm d}x_0 = 0.$$ (A.2)

Next we calculate the integral $\int \partial Q/\partial o(x_0) \cdot o(x_0) \, {\rm d}x_0$ using Eq. (3):

$$\int \frac{\partial Q}{\partial o(x_0)} \cdot o(x_0) \, {\rm d}x_0 = \nonumber \\$$

$$\qquad 2 \int \frac{o(x_0)\otimes p(x_0) - i(x_0)}{\sigma^2(x_0)} \otimes p(-x_0) \cdot o(x_0) \, {\rm d}x_0 \nonumber \\$$

$$\qquad=2 \int \int \frac{o(x)\otimes p(x) - i(x)}{\sigma^2(x)... ...t p(x-x_0) \cdot o(x_0) \, {\rm d}x\,{\rm d}x_0 \nonumber \\$$

$$\qquad= 2 \int \frac{o(x)\otimes p(x) - i(x)}{\sigma^2(x)} \cdot o(x)\otimes p(x) \, {\rm d}x \nonumber \\$$

$$\qquad= 2 \int \frac{\vert o(x)\otimes p(x)\vert^2 - i(x)\cdot o(x)\otimes p(x)}{\sigma^2(x)} \, {\rm d}x.$$ (A.3)

Because of the noise in the raw data i(x), we can get many different reconstructions which are in agreement with the measurement i(x) within the error $\sigma(x)$ of the measurement. Therefore, we can consider o(x) and also $o(x)\otimes p(x)$ as a random variable with expectation value

$$\langle o(x)\otimes p(x)\rangle = i(x),$$ (A.4)

and second moment

$$\langle \vert o(x)\otimes p(x)\vert^2\rangle = \vert i(x)\vert^2 + \sigma^2(x).$$ (A.5)

With these assumptions the average of the integral $\int \partial Q/\partial o(x_0) \cdot o(x_0) \, {\rm d}x_0$ is calculated to be

$$\bigg\langle \int \frac{\partial Q}{\partial o(x_0)} \cdot o(x_0) \, {\rm d}x_0 \rangle= \nonumber \\$$

$$\qquad 2 \int \frac{\langle \vert o(x)\otimes p(x)\vert^2\ran... ...x)\otimes p(x)\rangle }{\sigma^2(x)} \, {\rm d}x \nonumber \\$$

$$\qquad= 2 \int \frac{\vert i(x)\vert^2 + \sigma^2(x) - \vert i(x)\vert^2}{\sigma^2(x)} \,{\rm d}x \nonumber \\$$

$$\qquad= 2 \int \frac{\sigma^2(x)}{\sigma^2(x)} \, {\rm d}x = 2 \int 1 \, {\rm d}x = 2\cdot \langle Q\rangle,$$ (A.6)

where $\langle Q\rangle$ ist the mean of the $\chi^2$ function defined in Eq. (1). Using Eqs. (A.2) and (A.6) the Lagrange multiplier $\mu$ can be expressed by

$$\mu = - \frac{2\cdot \langle Q\rangle }{\int \partial H/\partial o(x_0) \cdot o(x_0) \, {\rm d}x_0}\cdot$$ (A.7)

For the calculation of $\mu$ from this equation a reasonable estimate of the true object is required. Our experiments show, that the measured image i(x) is such a good estimate.

  
Appendix B: The advantage of Fourier-masked LBT raw images

In this section we investigate the benefit of the preprocessing step (3) in Sect. 3.2, i.e., setting of all Fourier elements of an LBT raw image outside its uv-coverage to zero (Fourier-masking of the LBT raw data). The reconstruction results listed in Table 2 were obtained from computer-generated LBT raw images (see simulation parameters in Table 1) after application of this Fourier-masking step. The restoration errors of the reconstructions obtained from the same LBT raw images, but without application of the Fourier-masking step, are listed in Table B.1. These computer simulations show that all images reconstructed with the RL, OSEM, and BB methods have nearly 10%-20% larger restoration errors if Fourier-masking is not applied. A similar result is also obtained for laboratory LBT raw images (see Fig. 5, and Table 3). The restoration errors of the reconstructions obtained from the laboratory LBT raw images with Fourier-masking are listed in Table 4, whereas the reconstruction results without Fourier-masking are summarized in Table B.2. As the comparison shows, in the case of the laboratory data the restoration errors of the reconstructions without Fourier-masking of the LBT raw data are typically 10% larger than the errors obtained using Fourier-masked raw data.

Table B.1: Restauration errors of the reconstructions obtained with the RL, OSEM, and BB methods from the computer-simulated LBT raw images of the galaxy (see Table 1) without Fourier-masking (see Sect. 3.2).

Table B.2: Restauration errors of the reconstructions obtained with the RL, OSEM, and BB methods from the laboratory LBT raw images without Fourier-masking (see Sect. 4.2).

  
Appendix C: Space-variant BB method

The convolution between the object intensity distribution o_k(x) and the space-variant PSF $p(x;\xi)$ can be described by

$$o_k(x_0)*p(x_0;\xi) = \int o_k(\xi)\cdot p(x_0;\xi) \, {\rm d}\xi,$$ (C.1)

where * denotes space-variant convolution. The space-variant version of $\partial Q/\partial o_k(x_0)$ in Eq. (3) is given by

$$\frac{o_k(x_0) * p(x_0;\xi)}{\sigma^2(x_0)} * p(-x_0;\xi') - \frac{i(x_0)}{\sigma^2(x_0)} * p(-x_0;\xi')=$$

$$\qquad \int \int \frac{o_k(\xi)\cdot p(\xi';\xi)}{\sigma^2(\xi')} \cdot p(-x_0;\xi') \, {\rm d}\xi' {\rm d}\xi -$$

$$\qquad- \int \frac{i(\xi)}{\sigma^2(\xi)} \cdot p(-x_0;\xi) \, {\rm d}\xi$$

$$\qquad= \int o_k(\xi)\cdot f(x_0;\xi) \, {\rm d}\xi - \int \frac{i(\xi)}{\sigma^2(\xi)} \cdot p(-x_0;\xi) \, {\rm d}\xi,$$ (C.2)

where $f(x_0;\xi)$ is

$$f(x_0;\xi) := \int \frac{p(\xi';\xi)}{\sigma^2(\xi')} \cdot p(-x_0;\xi') \, {\rm d}\xi'.$$ (C.3)