© ESO, 2011

1. Introduction

Many astrophysical studies require milli-arcsecond (mas) resolution images at optical wavelengths (visible and infrared), for example, the understanding of the interplay between accretion and ejection in the inner part of the disks of young stellar objects, the expansion mechanisms in novae just a few hours or days after the explosion, and the nature of dust in active galactic nuclei. Information at such a high resolution at the optical wavelengths requires diffraction-limited images with pupil sizes of the order of tens to hundreds of meters that can only be achieved by interferometrically combining light from separate apertures. The Very Large Telescope Interferometer (VLTI) and the Center for High Angular Resolution Array (CHARA) are facilities that provide interferometric measurements that can be used to reconstruct images of stellar surfaces (e.g., Monnier et al. 2007; Haubois et al. 2009; Zhao et al. 2009), binaries (e.g., Zhao et al. 2008; Kraus et al. 2009), circumstellar shells around evolved stars (e.g., Le Bouquin et al. 2009), and the close environment of young stars (Renard et al. 2010; Kraus et al. 2010).

Owing to the sparse (u,v) coverage, the image reconstruction process is ill posed as there are more unknowns, e.g. the pixels of the image, than measurements. The data alone are insufficient to reconstruct an unambiguous image and some additional constraints, so-called the regularizations, are needed to converge to a unique and stable solution. Compared to radio interferometry, the data are much sparser in optical interferometry; hence, we expect the image reconstruction problem to be much more sensitive to the choice and the tuning of the regularization. Since the general study of regularization by Titterington (1985), many different methods have been proposed in the literature to adjust the regularization level. Since image reconstruction for optical interferometry is still in its infancy, it is fundamental to analyze the different types of regularization to find those that are the most suitable for the different astrophysical problems and to be able to tune the weight of the regularization. In this context, we carried out systematic tests using the image reconstruction algorithm devoted to optical interferometry data developed by Thiébaut (2008), called the Multi-aperture image Reconstruction Algorithm (MiRA). The analysis of these tests allow us to extract some general conclusions and establish practical rules for the users.

The mathematical principles of the image reconstruction technique is presented in Sect. 2. The parameters of the simulated data are presented together with the characteristics of the images and the strategy in Sect. 3. The results of the simulations are presented and discussed in Sect. 4, with an analysis of the role of the different terms and parameters in the image reconstruction. Finally, our conclusions are summarized in Sect. 5.

2. Principles of image reconstruction from optical interferometric data

We do not intend to provide a formal and precise description of image reconstruction in optical interferometry, but instead sufficient details to ensure that the paper is self-contained. Readers interested by a more detailed description of the method should refer to Thiébaut & Giovannelli (2010).

2.1. Data from optical interferometric observations

The principle of interferometry is to recombine coherently the beams from two or more independent telescopes and measure the so-called complex visibilities of the fringe patterns produced by the interferences. According to the van Cittert-Zernicke theorem, for an ideal interferometer, the complex visibility V_j1,j2(t) of the fringes produced by the interferences of the telescopes j₁ and j₂ at time t is proportional to the Fourier transform of the object brightness distribution Î(ν_j1,j2(t)) at spatial frequency , where λ is the wavelength, and the so-called baseline represents the separation between the two telescopes projected on a plane perpendicular to the line of sight (Lawson 2000; Malbet & Perrin 2007). Since the number of measurements is finite, to simplify the equations we introduce some notation for the mth measured complex visibility and the corresponding spatial frequency given by where j_1,m and j_2,m are the interfering telescopes and t_m is the time of observation.

2.2. Description of the image model

The final product of the image reconstruction is an image that can be treated as a grid of square pixels. In this context, the object brightness distribution as a function of the position θ can be approximated using the parametrization (3)where are the image parameters, e.g. the pixel values of the image, and is the chosen basis of functions, e.g. the response function of each pixel. The image reconstruction then consists of estimating the N parameters x that most closely fit the interferometric data. In this paper, we chose x_n to be proportional to the value of the nth pixel of a sampled image and b_n(θ) = b(θ − θ_n), where b(θ) is the pixel shape and θ_n the position of the nth pixel; thus (4)where α > 0 is a scaling factor such that x is normalized (this is required by the interferometric data format, cf. Pauls et al. 2005). With this model, the exact Fourier transform of the brightness distribution is given by (5)where and are the Fourier transforms of the basis functions. In our case, they correspond to the Fourier transform of the pixel response function, i.e. the pixel shape. Hence, the model of the mth complex visibility is given by (6)where A is a matrix with the complex coefficients (7)This matrix multiplication performs a linear transformation that contains the Fourier transform, the pixel shape, and the sparse sampling of the (u,v) plane.

The main problem in optical interferometry is the small number of telescopes (currently up to four or six), which leads to a sparse sampling of the spatial frequencies, the so-called (u,v) plane (see Fig. 2). Owing to random effects caused by the atmospheric turbulence, the visibility phase cannot be calibrated and the power spectrum and the closure phase are used. This results in a partial loss of the Fourier phase information (Thiébaut & Giovannelli 2010).

2.3. Inverting the problem of interferometric imaging

Because of the sparse (u,v) coverage and the possible lack of other information such as the phase, the reconstruction of an image obtained by the interferometric data alone is an ill-posed inverse problem. It needs additional a priori constraints to be recasted into a problem that has a unique and stable solution. A general prescription is to express the solution as one that minimizes a penalty function f under some strict constraints (Thiébaut 2005; Thiébaut & Giovannelli 2010) (8)with (9)where the so-called likelihood term f_data(x) measures the discrepancy between the model and the available data, while the so-called regularization term f_prior(x) measures the discrepancy with the prior information. In other words, minimizing the likelihood term f_data(x) enforces the fit with the actual data, while minimizing the regularization term f_prior(x) enforces the agreement with the priors. The so-called hyperparameter μ > 0 is used to adjust the relative weight of the constraints set by the measurements and the ones set by the priors. In Eq. (8), the positivity () and the normalization (∑ _nx_n = 1) of the brightness distribution are also included by default.

For the image reconstruction, we use MiRA (Thiébaut 2008) to find solutions of Eq. (8). MiRA can deal with the various kinds of data provided by an optical interferometer and implements a number of different regularizations (Thiébaut 2008; Thiébaut & Giovannelli 2010).

2.3.1. The likelihood term f_data and the data model

We focus on the choice of the priors and the tuning of their parameters. It is not within the scope of the paper to deal with global optimization issues and the search for a global minimum of the penalty function in Eq. (9). We therefore assume that the available measurements consist of complex visibilities, i.e. amplitude and phase, in order to have a convex likelihood term f_data(x). If the regularization term f_prior(x) is also convex, the global penalty function f(x) will be convex, which ensures that the solution of Eq. (8) is unique. Current optical interferometers only provide phase closures and power-spectrum data (i.e. the phase of the bispectrum and the squared amplitude of the complex visibilities), this means that our assumption will give somewhat optimistic results because some Fourier phase information is missing and because the likelihood term f_data(x) is non-convex when dealing with real data. However, as the number of simultaneous interfering telescopes increases, the number of missing phases becomes much less important and they can be reliably derived using self-calibration (Pearson & Readhead 1984) to achieve a situation similar to the case studied in our simulations. Moreover, new interferometers will make use of a phase reference source to directly measure the phase of the complex visibilities (Delplancke et al. 2003).

However, the OI-FITS standard (Pauls et al. 2005) imposes the use of complex visibilities in their polar representation with independent error bars. We therefore simulate each measured complex visibility as (10)where ϱ_m and ϕ_m are the measured amplitude and phase of the mth measure, V_m the corresponding complex visibility computed from the true object brightness distribution, and δϱ_m and δϕ_m are additive noise terms. In our simulations, the noise terms have independent Gaussian statistics such that (11)where  ⟨ ϱ_m ⟩  = |V_m| is the expected value of the amplitude that can be computed from the complex visibility of the true image. This particular choice follows approximately the model of Goodman (1985). In this context, to define the likelihood term we use the local approximation (Meimon et al. 2005) (12)where r_//,m(x) and r_⊥,m(x) are the two components of the complex residuals, r_m(x) = ϱ_m   e^i   ϕ_m − (A·x)_m, respectively, along and orthogonally to the measured complex visibility. Given the error bars σ_ϱ,m and σ_ϕ,m of the amplitude and the phase of the complex visibility, the standard deviations in the components of the complex residuals are (Pauls et al. 2005) (13)

2.3.2. The regularization term f_prior

In our simulations, we test 11 different regularization terms that are commonly used in image reconstruction methods and are implemented in MiRA (see Appendix A for detailed expressions).

• 1.

  Quadratic smoothness, which most closely describes a smoothimage and helps us to avoid unmeasured high frequencies.

• 2-3.

  Compactness, which describes compactness in the image plane and hence smoothness in the Fourier plane (Le Besnerais et al. 2008). Two different cases were studied in the simulations, with penalties of the second and third orders with respect to the distance of the center of the field of view (FOV).

• 4.

  Total variation (TV), which minimizes the total gradient of the image and helps us to describe uniform areas in the sought image with steep but localized changes (Strong & Chan 2003).

• 5.

  ℓ₁smoothness, which is useful for an extended object with sharp edges since it is linear for strong gradients.

• 6-8.

  ℓ_p-norm with p = 1.5, p = 2 and p = 3. For p > 1, the ℓ_p-norm regularization tends to produce a smooth image as it reduces the variance in the pixels.

• 9-11.

  Maximum entropy methods (MEM) aims to identify the least informative image consistent with the data (Gull & Skilling 1984; Narayan & Nityananda 1986). We try three types of entropies MEM-sqrt, MEM-log, and MEM-prior, respectively. The first two tend to reproduce an image with a flux spread across a minimum number of pixels. The last one is minimum when the image is as close as possible to a prior image. This prior image is the Gaussian that most closely reproduces the amplitude visibility data.

The different regularization terms are expected to behave as follows. The positivity and the normalization imposed in all the reconstructions are an ℓ₁ norm and lead to the sparsity of the solution, i.e. to a minimum number of bright pixels to explain the data. As most astrophysical images are smooth and compact, we expect that the regularizations that can describe these images will behave well, i.e. smoothness (quadratic or ℓ₁), compactness, and TV. The ℓ_p-norm regularization with p = 2 (and by extension for p > 1 as the regularization has the same behavior) has the tendency to force to zero the spatial frequencies that have not been measured (according to the Parseval theorem). Since the regularization has to interpolate correctly between the data, which is closer to the reality, it is not expected to give good results. Finally, the MEM-prior is expected to yield more reliable results than the other type of entropy because our choice of the a priori image can closely describe a compact object.

3. Description of the simulations

We now describe all the various simulated data that we compute for different objects, (u,v) coverages, and signal-to-noise ratios (SNR), as well as the parameters used for the image reconstructions.

3.1. Simulated data

Our simulated data sets are saved as OI-Fits files (Pauls et al. 2005) and depend on several setting of the object type, the (u,v) coverage, and the SNR. The 90 data files are available on the JMMC website¹.

3.1.1. Astrophysical objects

For our simulations, we consider ten astrophysical objects (see Fig. 1) that differ in term of their morphology and the typical length scales of their structures.

• 1.

  LkHα: the model of LkHα describes a compact object with a peak of intensity and a smooth envelope. The model comes from the 2004 International Imaging Beauty Contest organized by Lawson for the IAU (Lawson et al. 2004).

• 2.

  A stellar surface: this is a model of the supergiant α Ori that presents some convective cells, producing small scales on a smooth background. This model was produced by Chiavassa for the 2009 Workshop on Interferometry Imaging (WII09) organized by Berger and Malbet (Berger et al., in prep.).

• 3.

  A stellar cluster: the model consists of a hundred point sources. The position and the brightness of the sources follow a normal law.

• 4.

  Eta Carinae: the image of Eta Carinae presents many different scales and structures, such as the extended gas and the stars. This image was retrieved from the Hubble Space Telescope’s website². Some treatments have been applied to the image, i.e. a mean of the three different color channels to produce a grayscale image and a cut of the low intensities to produce a null background.

• 5.

  A protoplanetary disk: the model represents an Herbig Ae/Be star with a point source (the star) and an extended structure (the disk). This model was computed by J.-P. Berger for the phase A science case of the VLTI-Spectro Imager instrument (Filho et al. 2008).

• 6.

  A limb-darkened star: we used the power-law model of Hestroffer (1997) with an exponent α = 0.3. The image has a very smooth core with steep edges.

• 7.

  The galaxy system M 51: this image of M 51 consists of as many different scales and structures as Eta Carinae (gas, stars, spiral arms). This image was retrieved from the Hubble Space Telescope’s website and was processed in a way similar to the image of Eta Carinae.

  Fig. 2 (u,v) coverage. From left to right: rich (245 sampled frequencies), medium (88 sampled frequencies), and poor (31 sampled frequencies).

• 8.

  The AGN M 87: this AGN has a jet that consists of a narrow structure surrounded by a smooth background due to the gas. The image was retrieved from the Hubble Space Telescope’s website and the same treatments as the Eta Carinae’s image have been applied.

• 9.

  A gravitational microlensing image: gravitational microlensing is an astronomical phenomenon due to the gravitational lens effect. When a distant star or quasar becomes sufficiently aligned with a massive compact foreground object, the bending of light due to its gravitational field leads to two distorted unresolved images resulting in an observable magnification. The image shows four very compact structures. This model was developed by J. Surdej for the phase A science case of the VLTI-Spectro Imager instrument (Filho et al. 2008).

• 10.

  A rapid rotator: the rapid rotation of a star affects the stellar shape and the local emitted flux. We used the model of Domiciano de Souza et al. (2002), with parameters D = 0.78 and T_p = 35   000 K. The resulting light distribution was projected onto a plane with an inclination of 45°.

Since the goal of our tests is to determine the influence of the object’s type on the image reconstruction, all the considered objects have a similar angular size of  ~34 mas, which is consistent with the typical size of the FOV of an optical interferometer such as Amber/VLTI. To generalize our results, we expect the most important figure for a given object and instrument to be the number of resolved elements, which is equal to the ratio of the angular size of the object support to the effective resolution of the imaging system. Estimated by this ratio, the complexity of the objects that we have considered is in the range of 200–600 resolved elements.

3.1.2. (u,v) coverage

To study the influence of the instrumental configuration and analyze the capability of the regularization to interpolate the available data and thus fill the voids in the (u,v) plane, we consider several sparse (u,v) coverages (see Fig. 2). To remain general, our different (u,v) coverages were constructed to be uniformly spread. Each (u,v) coverage consists of concentric rings modulated by a sinusoid along the ring and with phases of the modulations chosen to maximize the distance between the points of two adjacent rings. This also avoids symmetries in the (u,v) coverage. The concentration of the rings is more important on small baselines than on the largest ones to insure a good sampling of low spatial frequencies. In this paper, we consider three (u,v) coverages which differ in terms of the number of samples, called hereafter rich (245 sampled frequencies), medium (88 sampled frequencies), and poor (31 sampled frequencies) coverages. The chosen (u,v) configurations could be considered as very good by the standards of existing data, but the goal of the paper is not to cover all possible cases but to show the general trend.

3.1.3. Signal-to-noise ratio

To investigate the effects of varying the SNR, we use a SNR factor γ in the standard deviations given by (14)Therefore with these settings, the error bars become (15)To analyze the influence of the SNR on the reconstructed images, three values of γ are tested: high SNR (1%) , intermediate SNR (5%), and low SNR (10%). We limit our study to these SNR values to maintain a small amount of results. Therefore, there is no systematic attempt to search for the limit of SNR and we realize that our worst case (10%) might be considered moderate noise in real data.

3.2. Parameters of the synthesized image

The resolution and the FOV of the reconstructed image have to be chosen with care. On the one hand, the pixel size Δθ must be small enough to properly account for the highest spatial frequencies available in the data. Using Shannon’s rule, we find that , where B_max is the maximum length of the observed baselines and λ the wavelength. In our reconstructions, the pixel size is one third of the limit set by Shannon’s rule. This oversampling allows us to check whether image reconstruction can achieve some level of super-resolution. For instance, this yields Δθ = 0.4 mas at λ = 2.2   μm with B_max = 190 m. On the other hand, the FOV has to be large enough to avoid field aliasing and therefore we take an image three times larger than the object itself. As all our objects have a size of 34 mas and thus approximately fit into 85 × 85 pixels, the reconstructed images have 256 × 256 pixels.

3.3. Reconstruction strategy

Given the data (determined by the object, the (u,v) coverage, and the noise realization) and the regularization, we conduct a sequence of image reconstructions for different values of the regularization weight μ. Since the problems solved are convex (as explained in Sect. 2.3.1), their solutions do not depend on the starting point. To reduce the calculation time, we therefore try to use a starting point that is as close as possible to the final solution. We begin the sequences of reconstructions with the highest value of μ and use the true image as the starting point for this first reconstruction. We then reduce μ and use the image previously obtained as the starting point. This latter step is repeated until we reach the lowest value of μ. Each reconstruction is an iterative process that is stopped when the norm of the gradient of the penalty function f(x) is below a preset threshold. This threshold is derived from the true image according to (16)where η > 0 is a small value and f_μ(x_ref) the penalty computed for the true image and for each value of μ. The ℓ₂ norm of the true image is normalized and we assume that η = 10^-5.

3.4. Image quality criterion: the mean-squared error

To assess the quality of the reconstructed images, we consider the mean-squared error (MSE) of the reconstructed images. In our simulations, we use normalized input images with the same pixel size as for the reconstructed image. Hence, to compare the reconstructed image x^rec and the true image x^ref, we can simply use the MSE defined as (17)where N is the total number of pixels. We note that, since we use complex visibilities and our priors, apart from the FOV one, are shift-invariant, the reconstructed image is correctly centered and there is no need to compensate for registration errors.

4. Results and discussion

We performed a total of  ~24 000 simulations corresponding to the reconstruction of all cases described in Sect. 3 with different values of the hyperparameter μ. We present here the results of these simulations and analyze the consequences for the image reconstruction process in order to draw general conclusions.

We begin by discussing the optimal value of the hyperparameter and determining whether the MSE is a good quality criterion. We then discuss the effects of the following parameters:

• the regularization: what are the good and bad regularizations?

• the limits: what are the minimal (u,v) coverage and SNR value?

• the hyperparameter μ: what is the optimal value? With which parameters does it vary?

• the likelihood term: how does it vary? Can it be used to tune the regularization term instead of the hyperparameter μ?

• the effective resolution: what degree of super-resolution can be achieved? How does it vary?

4.1. Optimal regularization weight μ ⁺ 

We first investigate whether there is an optimal regularization weight μ for a given situation. Therefore, for each object, configuration, SNR level, and regularization, we reconstruct an image for different values of the hyperparameter μ.

The top row of the right panel of Fig. 3 shows the effects of different values of μ, whereas the left panel displays the obtained MSE (see Sect. 3.4) for each values of μ. For too small a value of μ, the under-regularized image (labeled with an A) has plenty of artifacts. In contrast, for too large a μ, the over-regularized image (labeled with a C) is blurred, and many fine features are lost. These two extreme situations correspond to high values of the MSE (A and C), but there is a situation where the MSE reaches a minimum and the image appears to have far fewer artifacts (B). Visual comparison of the A and C images with the one obtained for B confirms that the MSE is a good criterion to correctly set the regularization weight.

We conclude that there is indeed an optimal value of μ, the one that gives the smallest MSE (18)where is the image reconstructed with a regularization weight set to μ. As the minimum of the curve is quite flat, the optimal value of μ is not precisely defined but may vary by a factor as large as either two or three with a negligible influence on the result.

This procedure cannot be used in practice because the true image is obviously unknown. Nevertheless, this procedure allows us to define the most accurate image that can be reconstructed given the data and the type of regularization. In the following analysis, the reconstructed images are always obtained with μ = μ ⁺ .

4.2. Dependence of μ ⁺  on the MSE quality criterion

We wish to determine the type of relationship that links the optimized regularization factor μ ⁺  to the MSE in order to detect any trends. The scatter plot in Fig. 4 reports the values of μ ⁺  and MSE obtained for each simulation. In the left panel, the different colors and symbols correspond to different objects. In the right panel, the different colors and symbols correspond to different regularizations. In the bottom row, we present our computed histograms of MSE (left) and μ ⁺  (right). As in the rest of the paper, the plotted histograms are approximations of the probability density functions (PDF) of our results. These curves were computed from our samples following the optimal method described in Scott (1992).

In the left part of Fig. 4, the different colors representing the objects seem to be aligned vertically. We therefore conclude that the MSE depends mostly on the structure of the observed object.

More precisely, two classes of objects can be distinguished in the distribution of MSE in the left panel of Fig. 4. We therefore grouped together the objects with similar behaviors in the top panels of Fig. 5: (i) the objects with very compact sources, i.e. the star cluster, the protoplanetary disk and the microlensing (curve in blue labeled B); and (ii) the other objects with extended structures (curve in red labeled C). The MSE is systematically higher for objects of the first class.

Following these observations, we try to find a way of renormalizing the MSE as a function of the object. To join the curves together, we define a new MSE criterion, called MSE⁺, by dividing each MSE by the smallest MSE for each object separately. As expected, this normalization cancels the different object classes, as seen in the right part of Fig. 5. Two peaks clearly appear on the graph, distinguishing the good reconstructions (left peak) from the bad ones (right peak). An example of each case is shown on the bottom part of Fig. 5. The low quality reconstructions are caused by bad configurations (not enough (u,v) points, low SNR) and/or bad regularizations as it can be seen in the right part of Fig. 6 and will be described in the next sections. To be sure that the peaks represent the good and bad reconstructions and do not come from different object’s types, we verify that each object is represented in each peak, as shown in the left part of Fig. 6.

By visual inspection, we assess that the value of the MSE⁺ leads to a correct ordering of the images reconstructed for a given object when the other settings change (data quality, type of regularization, etc.): the lower the MSE⁺, the higher the quality of the image is. Other attempts at the renormalization of the MSE are explained in Appendix B.

Now that a good quality criterion has been defined and that the optimal value of μ have been obtained, we can study the other parameters.

4.3. Limits due to the (u,v) coverage and the SNR

In this section, we classify the observational configurations ((u,v) coverage and SNR) on the basis of the MSE⁺. For each data set (unique combination of object and regularization), we order the pair [(u,v) coverage, SNR] according to the value of MSE⁺ they reach, giving them a rank from 1 for the best configuration (lowest MSE⁺) to 9 for the worst (highest MSE⁺). In the left panel of Fig. 7, we display the cumulative distributions of the ranks reached by every configuration, determining how many times a given configuration reaches at least the first rank, the second rank, etc. The highest quality configurations are the ones in the upper-left part of the plot.

The poor (u,v) coverage combined with any value of SNR is clearly too sparse to reconstruct good images. While acceptable reconstruction is possible for any considered SNR when there are enough samples in the (u,v) coverage. We deduce that there is a minimal (u,v) coverage needed to perform image reconstruction, whereas there is no such limit set by the tested SNR. However, when the (u,v) coverage is sufficient, the higher the SNR or the more filled the (u,v) coverage, the higher the quality of the reconstructed image is. The bottom row of the right panel of Fig. 3 shows how the visual quality of the optimal image depends on the SNR. As expected, the higher the SNR, the better the reconstructed image is.

Figure 7 (right) shows the histogram of the MSE⁺ for different configurations of (u,v) coverage and SNR. It indicates that the success image reconstruction is more influenced by the amount of data than by the SNR: the MSE⁺ is lower for a rich (u,v) coverage with a poor SNR than for a medium (u,v) coverage with an intermediate SNR. We note that all the tested (u,v) coverage are uniform, but we expect that the amount of data has to be sufficient and also homogeneously spread in the (u,v) plan.

After the removal of the sparsest coverage, the MSE⁺ distribution is shown in Fig. 6 in dashed line. There is still a little bump of bad MSE⁺ caused by bad regularizations, as discussed in the next section.

4.4. Quality of the regularizations

Using the same principles as in the previous section, we classify the regularizations as a function of the MSE⁺ for different realizations. Figure 8 shows the corresponding cumulative distributions. Isotropic TV appears to be the most successful regularization for the two classes of objects. The compactness prior with is the second highest quality choice for the very compact objects. The worst regularizations are the ℓ_p-norm, the MEM-sqrt, and the MEM-log, as expected and explained in Sect. 2.3.2. Reconstructions for good and bad regularizations are illustrated in Fig. 9.

In the MSE⁺ distribution after the elimination of the bad regularizations (see Fig. 6 in dot-dashed line), there is still much evidence of lower quality MSE⁺. We conclude that the bad MSE⁺ are mainly caused by the sparsest (u,v) coverage. However, both have to be eliminated to obtain a cleaned sample of reconstructed images.

In what follows, we made a selection and only retained the regularizations (quadratic and ℓ₁ smoothness, compactness, isotropic TV and MEM with a Gaussian prior) and the (u,v) coverages (rich and medium) that lead to correct reconstructed images. We kept data for all values of SNR.

4.5. Predetermined value of the hyperparameter μ ⁺ ?

In a fully Bayesian framework, the data noise and the variability of the object are independent. We therefore expect that μf_prior does not depend on the data (Tarantola 2005). As a result, for a given type of object, the optimal value for μ should be the same regardless of the SNR or the (u,v) coverage. However, we are not in a truly Bayesian framework because the regularizations are derived from general considerations that must help us to solve for the degeneracies of the problem and are not really derived from the statistics of the brightness distributions of the observed objects. Since the degeneracies of the problem are mostly due to the sparsity of the (u,v) coverage, this observational parameter may have some influence on the tuning of the regularization weight. Our expectation is thus that, for a given type of object, a quasi optimal value of μ ⁺  can be derived from simulations and from simple considerations to scale this parameter if different image reconstruction settings are used.

The hyperparameter μ ⁺  depends mostly on the regularizations, as seen in the right part of Fig. 4, where the colors, representing the regularizations, appear aligned horizontally and the pre-eminent peaks of the histograms are separated. Moreover, the hyperparameter μ appears to be quasi independent of the SNR and the (u,v) coverage, as expected. As seen in Fig. 10 for the TV regularization, if there is still a variation in the value of μ, it depends on the object morphology but not on the (u,v) plane and the SNR value.

This finding allows us to link each regularization to a mean value of the hyperparameter μ. This mean value gives a useful start point for the user (see Table 1). The equations to rescale this value in the case of different image settings are given in Appendix C.

4.6. How different noise realizations affect the MSE and the optimal μ?

In all simulations, in order not to influence the classification of the results, the same random seed is used to compute the noisy data. Therefore, we test 100 noise realizations for each regularization in the case of the galaxy with a medium (u,v) coverage and an intermediate SNR to study its impact on the curves computed from a single realization. From Fig. 11, we conclude that the MSE is not very different, thus the image reconstruction does not depend on the particular noise realization (as shown for example in Fig. 12). Moreover, the optimal value of μ varies by less than an order of magnitude.

4.7. The effective spatial resolution

To investigate whether super-resolution is achievable and to quantify the amount of trustable super-resolution, we estimate the effective resolution of the reconstructed images. We define the effective resolution as the full width at half maximum (hereafter FWHM) of the Gaussian that yields the smallest value of MSE between the reconstructed image and the true image smoothed by that Gaussian (19)where G(w) is a linear smoothing filter that convolves its argument by a Gaussian of FWHM equal to w.

As shown in the left part of Fig. 13, the effective resolution varies with both the amount of data and the SNR: the more extensively the (u,v) parameter space is filled and the higher the SNR, the higher the effective resolution is. As for the MSE (cf. Sect. 4.3), the effective resolution is more influenced by the (u,v) coverage than by the SNR. Super-resolution can be achieved in only the best cases.

The bottom-right part of Fig. 13 shows the expected behavior of the effective resolution: the higher the hyperparameter μ, the lower the effective resolution is. In order words, the more the image is regularized, the fewer the details that are visible.

An interesting question is why this super-resolution exists. Although we do not have a definitive answer, we can speculate that it is due to the role of the positivity in the image reconstruction (Biraud 1969). This is confirmed by our simulations (cf. the upper right part of Fig. 13): without the positivity constraints, the distribution of the FWHMs has a peak higher than when using the positivity constraints.

4.8. Other methods for tuning the regularization?

For fully Gaussian statistics (i.e. all the likelihood and prior penalties are quadratic and there are no constraints such as the normalization and the positivity), the expected values of the total penalty function f(x ⁺ ) and the likelihood term f_data(x ⁺ ) at the optimal solution x ⁺  is given by (Tarantola 2005) where N_data is the number of data points (for both visibility amplitudes and phases), and edf is the number of equivalent degrees of freedom. Despite our being unable to use fully Gaussian statistics, the left panel of Fig. 14 shows that the distribution of f(x ⁺ ) peaks approximatively at the value of N_data. The spread of this distribution however prevents us to be able to tune the regularization level according to the criterion that f(x ⁺ ) ≈ N_data.

The value of f_data(x ⁺ )/N_data is around 0.1, which means that edf/N_data ≈ 90% of the data information is resolved. The image reconstruction is able to estimate almost as many parameters as data points. Since the f_data(x ⁺ )/N_data ratio has a smaller range (from  ~0.3 to  ~0.003) than the possible value of the hyperparameter μ (from 100 to 10¹³), it may be easier to tune the balance between the terms in the penalty function thanks to the f_data(x ⁺ )/N_data value instead of μ. However, this criterion may be really variable with the noise statistics and a study of the variation in f_data(x ⁺ )/N_data with different noise statistics should be done before using it as a tune factor (Gull 1988; Pichon & Thiebaut 1998).

Another way of determining the value of μ is the so-called L-curve. The L-curve is a log-log plot of the regularization term f_prior versus the likelihood term f_data for a range of the hyperparameter μ. In the L-curve criterion, the regularization parameter μ is such that the corresponding point on the L-curve lies in the corner. This choice is motivated by the corner being the separation between the flat part where the solution is dominated by regularization errors and the vertical part where it is dominated by the perturbation errors (Hansen 2000). The correct behavior of the L-curve is confirmed by the simulations, as seen in Fig. 15 (right): since the curves are plotted as a function of the highest quality image, they should cross in their corner, which is globally the case. We note that the shape of the L-curve seems more complicated as there are at least two corners and not only one (see Fig. 15 left). The L-curve appears to be an appropriate tool for finding the optimal value of μ but a more general study has to be done to confirm its suitability (see comparison with GCV) and its practical implementation in an image reconstruction algorithm.

5. Conclusions

Thanks to the use of a flexible algorithm devoted to image reconstruction in optical interferometry, we have performed a detailed study of the regularization terms. This study is the first one to compare such a number of regularizations on an equal footing, i.e. with the same algorithm and using the same data type. Performing these systematic tests has allowed us to discuss the different parameters and terms in the image reconstruction and to extract some practical rules, which are summarized in the following:

• 1.

  A minimal (u,v) coverage is necessary to reconstruct an image.Even if the image quality improves with the SNR, such a limit doesnot exist for the SNR. In other words, in the image reconstructiontechnique, increasing the number of telescopes is moreinteresting than constructing larger ones. The homogeneity ofthe (u,v) coverage is probably also critical but has not been tested.

• 2.

  Some regularizations are suitable for the optical image reconstruction and others not, regardless of the object being targeted. The holes in the (u,v) plane are the major issue in optical interferometry and the main role of the regularization is the correct interpolation of the missing data, regardless of the object. The highest quality regularization among the tested ones is the isotropic total variation.

• 3.

  The hyperparameter μ does not depend on the (u,v) coverage and the SNR as theoretically expected and depends mostly on the type of regularization. An optimal value for each regularization tested in this paper is given in Table 1. This optimal value may vary by a factor of 2–3 without there being any major changes in the images. A slight dependence on both the object structures and pixel size is also discernible and the equation to rescale the optimal values are computed. It should be interesting to implement regularizations independent of the pixel size.

• 4.

  Super-resolution can be achieved in the image reconstruction process and its level rises with the (u,v) coverage filling.

• 5.

  There are various possible ways of tuning the regularization level:

  • The visual tuning is enough as the μ value can slightly vary without causing any large changes in the reconstructed image.

  • Setting the likelihood term f_data seems to be a more effective way of fixing the balance between the regularization and the likelihood terms. However, the variation in the likelihood term with noise statistics needs to be investigated.

  • The L-curve criterion could give correct results.

---

Acknowledgments

This research has made use of Yorick, a free data processing language written by D. Munro (http://yorick.sourceforge.net) and the Service de Calcul Intensif de l’Observatoire de Grenoble.

References

Appendix A: The regularization expressions

In our simulations, we test 11 different regularization terms commonly used in image reconstruction methods and implemented in MiRA:

Quadratic smoothness:(A.1)where S is a smoothing operator implemented via finite differences.

Compactness (Le Besnerais et al. 2008): (A.2)which is a separable quadratic regularization. To enforce compactness, the weights have to increase with the distance from the center of the image. Two different cases were studied in the simulations: and , where  ∥ θ_n/Δθ ∥  is the distance in pixels of the nth pixel from the center of the FOV.

Total variation (Strong & Chan 2003): (A.3)where is the squared magnitude of the spatial gradient in the image, ϵ > 0 is a small number inserted to avoid the discontinuity in zero, and (n₁,n₂) ~ n are the two dimensional indices of the nth pixel. In our reconstructions, ϵ has always been chosen so as to be negligible compared to the gradient of significant structures.

ℓ₂ − ℓ₁smoothness (Charbonnier et al. 1997): (A.4)where ψ(z) = z − log (1 + z) is a ℓ₂-ℓ₁ norm, τ > 0 is a threshold level, and D is a finite difference operator approximating the qth spatial derivatives of its argument. When  ∥ (D·x)_n ∥  is much smaller than the threshold, τ²   ψ(∥ (D·x)_n ∥ /τ) ≈ 1/2    ∥ (D·x)_n ∥ ², while, for  ∥ (D·x)_n ∥  it is much larger than the threshold, τ²   ψ(∥ (D·x)_n ∥ /τ) ≈ τ    ∥ (D·x)_n ∥ . This regularization attempts to strongly smooth the weak spatial gradients and slightly smooth the strong gradients. In our simulations, we take D to approximate the spatial Laplacian and choose a threshold small enough such that the regularization behaved mostly like a linear smoothness.

ℓ_p-norm: (A.5)where ϵ > 0 is a small value introduced to avoid the singularity in zero when p ≤ 1. For p > 1, the ℓ_p-norm regularization tends to produce a smooth image as it reduces the variance of the pixels. For p < 1, the ℓ_p-norm regularization tends to promote sparsity in the image. This is interesting mostly for objects consisting of point sources. True sparsity constraints would be obtained for p = 0, although when p < 1 the regularization is no longer convex and the optimization problem becomes extremely difficult to solve as p gets closer to 0. Results in compress sensing (Candes et al. 2006) have proven that choosing p = 1 yields the most sparse solution, like p = 0 for a large class of problems. However, taking p = 1 yields a non-smooth but convex problem that is much easier to solve than the combinatorial problem resulting from the choice p = 0. With positivity and normalization constraints, the ℓ₁-norm of x is constant. Hence, taking p = 1 is meaningless in our framework and we consider only p = 1.5, p = 2 and p = 3.

Maximum entropy methods (Narayan & Nityananda 1986): (A.6)Here the prior is to assume that the image is drawn towards a prior model according to a non-quadratic potential h, called the entropy. We try three entropies: MEM was first introduced in radioastronomy and is useful for images made of bright point-like sources on a smooth background. In our simulations, we took the prior image to be the isotropic 2D Gaussian that most accurately fits the amplitude visibility data.

Appendix B: Renormalization of MSE

We attempted to normalize the MSE with two additional methods:

• 1.

  Since the squared difference between the real image and asmoothed version of the real image is higher for images with sharpor point-like structures, we computed a first normalized MSE as(B.1)where S is a smoothing operator. The distribution of this normalized MSE is shown in the right panel of Fig. B.1: the distribution is narrower but the two classes remain, despite the normalization.

• 2.

  In the second normalization, the MSE is compared to the norm of the reference image (B.2)This normalization is more effective than the previous one, as it joins the curves together. However, the distribution is unimodal and does not enable us to distinguish the good and the bad reconstructions. It is thus not a useful normalization.

Appendix C: Scaling the hyperparameter μ

In a Bayesian framework, the prior penalty μ   f_prior(x) should only depend on both the sought after brightness distribution I(θ) and the image parametrization. In this appendix and from these simple principles, we derive a method to adapt the value of the hyperparameter μ when the image parameters such as the pixel size are modified.

Using the sampled image model in Eq. (4), the normalization of x implies that where we have used the Riemann approximation of integrals, Δθ is the pixel size, and L = 2 is the number of dimensions in θ. Without any loss of generality, we can assume that the sought after brightness distribution is normalized such that ^∫I(θ)d^Lθ = 1; hence (C.1)and Eq. (4) becomes (C.2)

C.1. Separable ℓ_p norm

Using Eq. (4) and then, the Riemann approximation, the prior penalty for a separable ℓ_p norm regularization, after substituting in Eq. (A.5), becomes (C.3)which shows that for a regularization by a separable ℓ_p norm (C.4)Hence, with the normalization constraint, the optimal value of μ should be the same regardless of the pixel size for a regularization given by a separable ℓ₁ norm.

C.2. ℓ_p norm on the gradient

Using 1D notation to simplify the equations, the prior penalty for a regularization by the ℓ_p norm on the gradient is given by where ∂_θI(θ) denotes the partial derivative of the brightness distribution along the angular direction. In L dimensions and using the Riemann approximation, this gives which shows that (C.5)Applying this result for a regularization by quadratic smoothness in 2D, e.g. Eq. (A.1), we found that, with a normalization constraint, μ ∝ Δθ^-4.

C.3. Total variation

The preceding result, with p = 1, readily applies to regularization by the total variation, that is (C.6)We can also deduce that, if a relaxed version of TV is used, as in Eq. (A.3), the relaxation parameter ϵ must scale as the pixel size Δθ to have the prior penalty approximately insensitive to the pixel size.

We note that, with our particular choice of the threshold τ for the ℓ₂ − ℓ₁ smoothness regularization defined in Eq. (A.4), we expect this regularization to behave mostly like TV.

C.4. Quadratic compactness

The quadratic compactness we used in MiRA is given by Eq. (A.2) with q = 2 or 3. From this last approximation, we derive the scaling of μ with the pixel size (C.7)Hence, in 2D (L = 2) and for a normalized image, μ does not depend on the pixel size for q = 2, while it scales as Δθ for q = 3.

C.5. Maximum entropy

For maximum entropy methods, we have and

with as in Eq. (A.9). from which we deduce that (C.8)Hence, for a normalized image, μ does not depend on the pixel size in a MEM-prior regularization.

All Tables

All Figures