PAINTER: a spatio-spectral image reconstruction algorithm for optical interferometry

06/29/2014 ∙ by Antony Schutz, et al. ∙ 0

Astronomical optical interferometers sample the Fourier transform of the intensity distribution of a source at the observation wavelength. Because of rapid perturbations caused by atmospheric turbulence, the phases of the complex Fourier samples (visibilities) cannot be directly exploited. Consequently, specific image reconstruction methods have been devised in the last few decades. Modern polychromatic optical interferometric instruments are now paving the way to multiwavelength imaging. This paper is devoted to the derivation of a spatio-spectral (3D) image reconstruction algorithm, coined PAINTER (Polychromatic opticAl INTErferometric Reconstruction software). The algorithm relies on an iterative process, which alternates estimation of polychromatic images and of complex visibilities. The complex visibilities are not only estimated from squared moduli and closure phases, but also differential phases, which helps to better constrain the polychromatic reconstruction. Simulations on synthetic data illustrate the efficiency of the algorithm and in particular the relevance of injecting a differential phases model in the reconstruction.



There are no comments yet.


page 9

page 10

page 11

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Current astronomical Optical Interferometers (OI) do not directly provide images, even if prospective studies exist in this direction lab1996 ; aime2012 ; 2012SPIE2 . Instead, OI data are related to the observed celestial scene at the observation wavelength by a two dimensional Fourier transform (FT) of the corresponding intensity distribution perpendicular to the line of sight. Ideally, the observables measured by an interferometer are the so-called complex visibilities

, which corresponds to complex samples of the Fourier spectrum. The sampling function is fully defined by the positions of the interfering telescopes and by the observation wavelength. Earth rotation provides additional samples for observations acquired at different epochs, as the telescopes’ positions are modified in time with respect to the line of sight (an effect called super synthesis). However, owing to the small number of telescopes involved (typically three or four), the sampling of the Fourier space that is achieved by OI is always sparse.

In essence, measuring moduli and phases of complex visibilities in OI amounts to measuring contrasts and phases of interference fringes (labeyrie75, ). Atmospheric turbulence randomly shifts these fringes on a timescale of ms. The moduli are thus obtained by estimating fringe contrasts in snapshot mode with short integration times that freeze the atmospheric turbulence. As for the phases, the unknown random phase shifts imply that optical interferometers cannot measure directly the phases (in absence of an artificial or real star, which would provide such a reference monnier2003 ). Astronomers thus use a turbulence independent observable related to the phases called closure phase, originally devised for radiointerferometers jennison58 , and from which the phase information can be partially extracted VisEst1984 ; CP1986 ; VisEst1994 . Note that the current situation in optical is however very different from radio, where the numbers of antennas or dipoles is several orders of magnitude greater than in optical (see e.g. MPvanHaarlem:2013gi for a recent example) and all phases can be estimated.

From an informational viewpoint, the image reconstruction problem posed by OI is highly underdetermined because of the sparse Fourier sampling and of additional missing phase information. This means that formally, infinitely many intensity distributions exist that are consistent with OI data. In this framework, a classical and well-understood strategy for image reconstruction is to adopt an inverse problem approach, where missing information is mitigated, and hopefully compensated for, by a priori knowledge thiebauteas . In this case, the image reconstruction algorithm aims at finding an intensity distribution that minimizes a cost function composed of a data fidelity term, which is related to the noise distribution, plus a regularization term and possibly other constraints, which are related to prior knowledge.

Following this path, various algorithms have blossomed in the last twenty years. Most of the proposed algorithms rely on gradient descent methods (WISARD meimon2005 ; Meimon:05 , BSMEM BSMEM , MiRA MIRA , BBM BBM , IRBis IRBis ). A different approach is used in MACIM MACIM and in its evolution SQUEEZE baron2012

, which rely on Markov Chain Monte Carlo (MCMC) method.

While these algorithms have proved very useful to astronomers, none of them is currently able to tackle polychromatic reconstruction for general sources having variations in their intensity distribution in wavelength. As such, they are monochromatic imagers and are not applicable for multiwavelength image reconstruction.

Nowadays, modern OI is however polychromatic (see for instance AMBER AMBER2000SPIE , PIONIER PIONIER2011 , or VEGA Vega2009 ) and more powerful polychromatic instruments are in development (like MATISSE MATISSE and GRAVITY GRAVITY ). In such devices, interference fringes are simultaneously recorded in a number of wavelength channels that can reach several hundreds. Indeed, the fundamental justification of polychromatic observations is that the amount and distribution of electromagnetic radiation emitted by astrophysical sources may be very different from on wavelength channel to another, as light-matter interactions are highly variable in wavelength (emission and absorption lines for instance). Multiwavelength OI thus constitutes an extremely rich evolution over monochromatic OI. In order to fully exploit these instruments, the multiwavelength evolution of the mature but monochromatic OI reconstruction algorithms is mandatory.

Building on the monochromatic approaches mentioned above, some first steps have recently been undertaken in the direction of multiwavelength reconstruction. The work sparco implements a semi-parametric algorithm for the image reconstruction of chromatic objects, dedicated to the case of a central object surrounded by an extended structure such as a young star. The approach of admmthiebaut ; admmsoulez uses a sparsity-regularized approach, dedicated to the case where the observed scene is a collection of point-like sources. SelfCal uses MiRA as a key optimization engine combined to a Self Calibration approach and demonstrates the potential of using the differential phases in an image-reconstruction process.

This paper is devoted to the derivation of a multiwavelength or spatio–spectral images reconstruction (3D) algorithm named PAINTER (for Polychromatic opticAl INTErferometric Reconstruction software). This approach uses the absolute visibilities and closure phases, which are considered dependent of the wavelength. In addition, we also use the so-called differential phases, which are defined as the phases relatively to a reference channel and constitute an additional turbulence independent observable of the phases in multiwavelength observation mode.

The algorithm relies on the alternate estimation of the complex visibilities (from estimated phases and observed noisy moduli) and of the polychromatic intensity distribution (using spatio-spectral regularization and constraints). From a modeling viewpoint, the main originality of the approach is to estimate the unknown phases from both closure phases and differential phases. From an optimization viewpoint, the algorithm is based on ADMM methodology ADMMBoyd . PAINTER can be seen as an evolution of the MiRA–3D algorithm proposed in admmthiebaut ; admmsoulez . An implementation of PAINTER in matlab (octave compatible) with input data in OI-FITS format OIFITS , is available at

The paper is organized as follows: sections II and III introduce notations and data modeling. We derive here an extended model of phase differences which is specific to the 3D reconstruction. Section IV tackles the inverse problem approach. We introduce assumptions related to the data fidelity criterion (noise perturbations) as well as prior knowledge in the form of regularizations and constraints. Section V derives the resulting 3D image reconstruction algorithm. Performances of the algorithm are provided and analyzed in section VI.

Ii Notations

   phases of complex numbers

   (may result in a scalar, a vector or a matrix)

   moduli of complex numbers
   (may result in a scalar, a vector or a matrix)
   complex conjugate
   Fourier transform
   complex conjugate transpose
X    matrix
x    vector, with components

   Kronecker product (tensor product)

   direct sum:
   Hadamard product (dot product)

   identity matrix of size

   vector of length
   indicator function on positive orthant
   trace of X
   matrix vectorization i.e the columns of X
   are stacked into one column vector x
   diagonal matrix with x on its diagonal
   squared Frobenius norm
   update of x

Iii Data modeling

iii.1 Spatio-spectral model

In the absence of atmospheric turbulence the observable measured by an interferometer is the complex visibility thiebauteas . This observable is measured from the fringe pattern obtained by the interference of two beams collected from a pair of telescopes. The spatial position of each such pair defines one of the baselines of the telescope array. Hereafter, the notation for the baseline refers to the position vector of a telescope pair projected on a plane perpendicular to the line of sight.

In the considered case of polychromatic observations, an astrophysical source is described by an intensity distribution which is a function of wavelength. Because OI instrument always have limited fields of view, we assume that the distribution of interest accounts for the limited spatial response of the interferometer: it is an apodized version of the intensity distribution around the line of sight. The unkown distribution can be written as , which is a flux density at angular position of the sky and wavelength .

In absence of any perturbation and for purely monochromatic observations, a telescope pair of baseline provides a complex visibility defined by . Because OI instruments have limited angular and spectral resolutions (respectively set by the maximum distance between two telescopes and by the bandpass of the optical filters), a simple way to represent the unknown spatio-spectral distribution of the sources is to discretize over voxels. We consider here for simplicity the same discretization in both angular coordinates, resulting in parameters per image, and an instrument with wavelength channels of equal bandwidth , which is set to the spectral resolution. In this case, all voxels have the same size . If we further assume a unit transfer function in all channels, one voxel is simply defined as:


where refers to a pixel at angular position and to the reference wavelength of channel . The column vector collects all voxels and can be organized as the vectorization of a image of the astrophysical source at wavelength . In this setting, the goal of the multiwavelength reconstruction algorithm is to estimate the voxels, which represent the unknown parameters of the model.

Let be the complex visibility at the spatial frequency , and let be the column vector collecting the set of complex visibilities corresponding to all available baselines at wavelength . The complex visibilities can then be related in matrix form to the parameters by the direct model Thiebaut2010 ; admmthiebaut


where is obtained from a Non Uniform Discrete Fourier Transform (NuDFT) fessler2003 at the spatial frequencies imposed by the geometry of the telescope array and by the observation wavelength .

The previous expression describes the complex visibilities by wavelength. A compact notation including all wavelengths and baselines is

y (3)

where F is a block diagonal matrix with each block referring to the NuDFT at a particular wavelength. Vector y concatenates the complex visibility vectors ( of Eq. 2) for all wavelengths into a visibility vector, with associated moduli and phases given by

y (4)

In order to analyze the chromatic variation of the visibilities and of the images over the wavelengths, we also need to introduce the matrix Y and the matrix X defined by

Y (7)
X (8)

To clarify the use of a matrix notation note that the column of X corresponds to the vectorization of the image at the wavelength while the line is for the variation of the pixel along the wavelengths.

iii.2 Model of phase differences

In the presence of atmospheric turbulence, the beams received at each telescope are affected by random and different optical path differences, which corrupt the phases measurements of the complex visibilities. To overcome this difficulty, turbulence independent quantities need to be constructed.

iii.2.1 Closure phase

The first phase difference information used for image reconstruction in presence of turbulent measurements is the closure phase (operator ). It is defined as the phase of the bispectrum CP1986 , i.e., the Fourier transform of the triple correlation. For three baselines corresponding to a triplet of telescopes, the “atmospheric corrupted” instantaneous visibilities at a given wavelength can be modeled as


where are the quantities of interest (i.e., the uncorrupted phases) and are perturbation terms related to the corresponding telescopes. The closure phase associated with this triplet is defined as


where is as in Eq. 6 the vector containing all unperturbed phases for wavelength , and is a sparse row vector with only three non zeros entries which take values (as reflected by Eq. 12). Clearly, the phase closure allows to get rid of atmospheric effects for triplets of complex visibilities. If denotes the number of telescopes there are independent closure phases AdvancedImaging .

Let matrix concatenate in its rows the independent closure phases of the type (13) that can be obtained for the available triplets of telescopes at wavelength . For the simplicity of the presentation but without loss of generality, we assume here that all telescope pairs observe the same channels. In this case, the triplets involved in independent closure phases can be taken as the same in each wavelength channel. The global closure phase operator is then simply a block diagonal matrix that replicates , the closure phase matrix for


where is the vector of all closure phases and is the unknown unperturbed phase vector of Eq. 6.

iii.2.2 Differential Phases

A second phase difference information that is particularly interesting in polychromatic imaging is the differential phase . For one baseline , differential phases measure the phase evolution in wavelength with respect to a phase reference channel (see e.g. AMBER ). They are defined as

If the analyzed bandwidth is relatively narrow, which is generally the case, the phase turbulence terms on each telescope and are in first approximation independent of the wavelength DiffVis2012 . Thus the phase difference becomes


As closure phases, the differential phases in Eq. 16 are essentially not affected by the atmospheric perturbation.

The reference channel can be chosen as one of the available wavelengths. In this case differential phases are available per baseline. Without loss of generality, we denote below by the reference channel. Similarly to Eq. 13 we can write using definition of in Eq. 6:


where is a sparse row vector with only two non zeros entries which take values . Denoting by the differential phases vector which collects all wavelengths differences for all baselines, and by the matrix which concatenates all vectors in its rows, we have similarly to Eq. 14


Equations 14 and 18 can be merged in the single measurement equation


where models all phase difference information.

iii.3 Practical remarks

In practice, a specific phase reference channel is sometimes computed using all the visibilities associated with a baseline. One possibility is to set the phase reference as the angle of the empirical mean over the channels:

Another possibility, as for the AMBER instruments DiffVis2012 ; millour2006 is to define a channel dependent reference:

In these cases, the formalism of Eq. 16 is still valid assuming that all the reference channels for a given baseline are almost equal and will cancel in the difference of differential phases. Equation 16 can then be replaced by:


and the linear Eq. 18 is kept with the same matrix .

Consider the general case where telescopes configurations are used to observe at wavelengths. For a single exposure the number of spatial visibilities measured and so the number of phases to estimate is , where the number of bases . If the matrix H involved in the phases to the phases differences transformation is full ranked then there are no phases identifiability ambiguities. In practice, this matrix is the concatenation of two matrices, the first for the phase closures and the second for the differential phases.

A closure is said “closed” if the set of spatial frequencies involved in its construction are the vertex of a polygon. The minimal configuration to measure one closure phase is to use three telescopes, as the addition of several triplets of telescopes can describe al the possible polygon the information is clearly redundant. The rank of , the matrix which connects the phase closures to the polychromatic phases is given by the number of lines which ensure independent phase closures. For a single analyzed wavelength the rank of the matrix of section III.2.1 is meimon2005 . When polychromatic data are measured the augmented matrix follows the right term of Eq. 14, the phase closures for each wavelength are estimated from the same set of telescopes triplets. Consequently, the rank of is then increased up by a factor : . The ratio depends only on the number of telescopes.

For the differential phases, the construction of non redundant measures is detailed in section III.2.2 and is driven by . The rank of this matrix is given by the number of independent differential phases that may result in the transformation from phases and is . The ratio depends only on the number of analyzed wavelengths.

As H is the concatenation of and its rank is given by the number of independent relation between these matrices. Consider three telescopes and two wavelengths , is chosed to be the differential phases reference channel, a linear relation which connects the phase closures of Eq. 13 to the differential phases of Eq. 17 is:


the difference of lines of for a fixed wavelength only defines redundant phase closures. The combination over and of all the triplets of bases, as in Eq. 21, expresses all the differential phases and so the rank of H is at least the rank of . A more general phase closures differences can be expressed in term of , leads to redundant information and adds nothing. Consequently, the only remaining independent information is monochromatic and corresponds to all the phase closures related to the reference channel so . Finally, the ratio depends on the number of telescopes and analyzed wavelengths.

The ratio of the number of independent phases difference for the two models (ranks of and ) and the augmented model (rank of H) to the number of phases is shown in Fig .1. The ratio is drawn for telescopes which corresponds to actual interferometers and up to . Note as an example that for the high resolution mode of AMBER AMBER , and for the standard uses of PIONIER PIONIER2011 .

Figure 1: Ratio of the number of independent phases difference to the total number of phases as a function of the number of telescopes and analyzed wavelengths.

The use of all phases differences information largely increases the rank of the matrix which is always greater than for the phase closures alone.

Iv Inverse problem approach

iv.1 Image reconstruction problem formulation

According to Eq. 19 and notations defined in Eqs. 5 and 6, a model for measurements of difference phases and squared moduli can be written as


where and account for noise and modeling errors. Classical assumptions on their distributions are considered here. The noise is assumed to be jointly independent and Gaussian meimon2005 ; LITPRO2008 , and the noise is assumed to be jointly independent and marginally Von Mises distributed MIRA .

Writing the opposite logarithm of the likelihood related to the squared absolute value and phase differences model, we obtain:


where, and are relative weighting terms and:



is the variance of

. The constants are related to the variance of by

where is the modified Bessel function of order mardia . In practice is computed inverting numericaly the previous equation. The variance of the closure phase and differential phases, if not provided by the instrument pipeline, are estimated assuming independence of phase measurements.

Image reconstruction can be seen as an inverse problem IImgercons ; InversePBTarantola05 ; InversePBThiebaut . The model which connects object to the measured data involves a NuDFT, as described in Eq. 3 and Eqs. 22 and 23. This transformation leads to a poor coverage of the spatial frequency plan and makes the problem ill-conditioned which requires tackling the image reconstruction as a regularized optimization problem RenardReg . We will adopt here an objective function of the form:


where the image x can be constrained to be spatially limited in the support , and further constraints such as non negativity can be added in , which contains all the regularization terms. The support constraint is not included in for technical reasons related to the ADMM methodology described below.

iv.2 Regularizations and constraints

OI images are by nature non negative and sometimes contain sources that are spatially localized. However, specifying the properties of the object parameters x only in terms of non negativity and spatial support is usually not constraining enough. It follows that the uses of regularization terms to emphasize some inherent a priori knowledge about the image structure is necessary.

Following the matrix notation for the 3D object as defined in Eq. 8, PAINTER in its current form accounts the following classical priors:

  • Ridge/Tikhonov regularization, motivated by the poor conditioning of the NuDFT operator.

  • Spatial or spectral Total Variation babacan09 .

  • Spatial or spectral smoothness admmthiebaut .

The support constraint is defined by the parameters space in Eq. 27 and the non-negativity constraint by the regularization term . In this description the regularization function in Eq. 27 writes

and are the matrices of finite difference, associated with the spatial and spectral regularizations respectively admmthiebaut . act on the column of X which are images processed independently while operates on the line of X to connects the pixels between wavelengths. and are matrix regularization terms that can be chosen as:


Finally, , and are hyper-parameters which control the weights of the associated regularization terms.

V 3D reconstruction algorithm

Owing to the unavoidable non convexity of the problem as defined by Eq. 27 (see e.g. in Meimon:05 ), the vast majority of image reconstruction algorithms rely on a descent optimization principle. So does PAINTER  by using the flexibility of the Alternate Direction Methods of Multipliers (ADMM) algorithm, which was already used in admmthiebaut ; admmsoulez to reconstruct stellar spectrum of point sources from complex visibilities. Within this framework the reconstruction algorithm will iterate as follow:

  1. Update the complex visibilities.

  2. Use the estimated complex visibilities and spatio–spectral regularization to reconstruct the polychromatic object (3D–images).

  3. Update the Lagrange multipliers.

Specific details are given in the rest of this section.

v.1 ADMM optimization algorithm

The ADMM framework allows to split complex problem into smaller and easier ones by introducing auxiliary variables. However, in this case additional terms have to be taken into account. The subproblems can be solved independently by means of proximal operators. Standard constraints are taken into account in PAINTER  and most proximal operators used in this section are known. The optimization problem of Eq. 27 where is given by Eqs. 2426 and the regularization term is given by Eq. IV.2 is equivalent to:

In this equation, the vector denotes the complex visibilities variables whose likelihood is associated with the measurements of visibilities squared absolute values, Eq. 25. The vector denotes the complex visibilities variables whose likelihood is associated with the measurements of the visibilities phases differences, Eq. 26.

Using the same approach for the three last regularization terms leads to the introduction of new auxiliary variables P, T and V and to replace each regularization term in Eq. V.1 by a constrained problem.

  • The non-negativity constraint and the support constraint in Eq. V.1 are replaced by

  • The total variation regularization in Eq. V.1 is replaced by

  • The spectral regularization in Eq. V.1 is replaced by


    the use of two auxiliary variables ensure spectral separability when solving for X.

The final optimization problem obtained by replacing Eqs. 31, 32 and 33 in Eq. V.1 is iteratively solved using the ADMM algorithm ADMMBoyd . Auxiliary variables related to the complex visibility: y, , have a proper Lagrange multiplier: (or , with columns , in matrix form), , and share the same augmented Lagrangian parameter . The auxiliary variables introduced by the regularization: P, T, V, S are associated with the Lagrange multipliers , , , and to the augmented Lagrangian parameters , and for V and S.

Denoting with a “+” superscript updated quantities, alternated minimization of the augmented Lagrangian gives:


with the definitions:


The update of the Lagrange multipliers are:


The proximal operators CombettesPesquet in Eqs. 34 and 35 are detailed in subsections V.2 and V.3 respectively and the proximal operators for the regularization terms in Eqs. 39, 40 and 42 are detailed in subsections V.4V.5.

v.2 Proximal operator for squared visibility

The proximal operator of Eq. 34 updates the estimation of the complex visibilities from the measured squared absolute visibilities. Replacing by its expression in Eq. 25, Eq. 34 separates on each component:


which separates again on the modulus and phase of as:


where and denote in this section the moduli and phases of .

The minimization of the previous fourth order polynomial for is obtained by computing the real roots of its derivative:


using Cardano’s method Jacobson09 . is the real positive root of Eq. 58 which minimizes the criterium Eq. 56. If Eq. 58 has no positive roots, then the polynomial to minimize is strictly increasing for and .

v.3 Proximal operator for phase differences

The proximal operator of Eq. 35 updates the estimation of the complex visibilities from the measured visibilities phase differences. This problem involves triplets and duets of visibility phases which, contrarily to the moduli, are not separable. Replacing given in Eq. 26 we obtain:


where and denote in this section the modulus and phases of . Minimization of Eq. 59 w.r.t. gives:


Eq. 60 is first replaced in Eq. 59 and the resulting function minimized w.r.t. . This minimization is carried out numerically using a gradient descent algorithm. A compact expression of the gradient of the function to be minimized is:


where and are diagonal matrices with elements and . Once is obtained, is given by . It is worthy to note that if at an iteration step of the descent algorithm , from Eq. 60 the corresponding term in the second part of Eq. 59 reduces to and the next estimate of will not depend on and . For this reason the initial estimate of is set at . The minimization in PAINTER  relies on the quasi-Newton algorithm implemented in minfunc , note that the non convexity of the problem leads to a local minimization to estimate .

v.4 Positivity and compactness

The solution of the constrained minimization in Eq. 39 ADMMBoyd ; CombettesPesquet is:


where is a binary mask matrix (composed of and ) which represents the support of the object (disk etc.), each column of is for a wavelength. If the object support is not constrained .

v.5 Spatial and spectral regularization

The proximal operator of Eqs. 40 and 42 is:

  • the soft thresholding ADMMBoyd if :

  • if :


where D can be T or V. and can be chosen according to a priori selected on image, if the desired spatial/spectral smoothness regularization must promotes brightness jump then otherwise . As an example, the proximal operator for the norm spectral smoothness regularization, as expressed in Eq. 42 with , is:


Note that according to the structure of matrix , see Eq. 47, the spatial proximal operator Eq. 63 will apply separately on images at different wavelengths. Similarly, according to the structure of matrix , the spectral proximal operators operates separately on the voxels. Finally, in order to scale properly the regularization parameters, the and of Eq. V.1 are divided by their number of non zero elements in and respectively. This normalization should makes the parameters independent of the size of the 3D–image.

Vi Simulations

Computer simulations are presented in order to illustrate the performances of PAINTER and the benefits of a combined use of phase closures and differentials phases. To do so, synthetic data set are used and the algorithm conforms to the complete chain needed in OI image reconstruction i.e going from files to 3D–image.

vi.1 Synthetic data

vi.1.1 Improvement related to phases differences

The first simulation aims at determining the improvement on the reconstruction due of the combination of phases differences model. The simulated noiseless object consists in two uniform disks with a decreasing radius with wavelengths. The intensity of each disk as a function of the wavelength is constant. The original cube is of size pixels channels. The instrumental configurations (geometrical position, wavelengths, acquisition time,…) have been generated with the ASPRO2 ASPRO software which simulates realistic interferometric observation and store the data into OIFITS files OIFITS . The simulation uses the configuration of the AMBER instrument at the VLTI AMBER with three telescopes. Height equispaced wavelengths in the range (high resolution) at five acquisition instants are analyzed. The spatial frequencies coverage, including the earth rotation effect, is shown in Fig. 4 for the spectral channel at . This results in the measurement of 129 squared absolute visibilities and 43 phase closures for each wavelength. The differential phases are calculated as in Eq. 20 and corresponds to measures. The object is shown in Fig. 4 in a channel per channel view.

Figure 2: Spatial frequencies coverage plan. Geometrical configuration the 2004 International Beauty Contest in Optical Interferometry.
Figure 3: Synthetic data, two disks with a decreasing radius along wavelength. Per channel view.
Figure 4: Left: A channel of the 3D image used in simulation. Right: Spatial frequencies coverage plan. Geometrical configuration the 2004 International Beauty Contest in Optical Interferometry.
Figure 2: Spatial frequencies coverage plan. Geometrical configuration the 2004 International Beauty Contest in Optical Interferometry.

vi.1.2 Resolved stars reconstruction

The second simulation deals with realistic noisy synthetic data with a given signal to noise ratio, SNR = 30 dB. More specifically, to generate the data we add zero-mean independent Gaussian random noise to each measurement. For each phase data

, the standard deviation of the noise is

= 1 rd/SNR (in radians). For each amplitude data , the standard deviation of the noise is = SNR. The simulation consists in two resolved stars of different spectral types: G8V and M3V for the large and small stellar disks respectively. The spectra and their chromatic limb darkening laws are based on Kurucz and van Hamme models Kurucz ; VH ; refId0 . A channel of the original 3D-image is shown in Fig. 4, . The two disks have a diameters of and milli-arcsecond (mas) and individual brightness distribution is varying among wavelengths. The original cube with a pixel resolution of mas is of size pixels channels, the field of view (FOV) is mas. The instrumental configuration was used in the 2004 International Beauty Contest in Optical Interferometry IBC04 : Thirty equi-spaced wavelengths in the range (high resolution) at 13 acquisition instants are analyzed. The spatial frequencies coverage, including the earth rotation effect, is shown in Fig. 4. This results in the measurement of 195 squared absolute visibilities and 130 phase closures for each wavelength. The differential phases are calculated as in Eq. 20 and corresponds to measures and simulation data are stored into OIFITS files OIFITS .

vi.2 Reconstruction parameters

The initial estimate was a Dirac function centered on the image for each wavelength and the algorithm was stopped after iterations.

Figure 5: Estimated 3–D object with closure phases only.

The type of regularization was chosen according to the a priori on the object: and to enforce spatial continuity with sharp edges and spectral smoothness. The regularization weights and , were empirically tuned to have a visually acceptable solution (though not the best one). The associated augmented parameters, and , which drive the convergence of the spatial/spectral regularization sub-problems have been tuned to be not too small (to avoid divergence) and not too large (to not slow down the convergence). Without loss of generality the problem can be dimensionless in and thus we took . We tuned the weights of the two terms of Eqs. 34 and 35 so that, at convergences, these terms are of the same order. This yields and . Finally, we took . The reconstruct cube is of size pixels per channel ( and for the first and the second simulation respectively), leading to parameters to estimate.

vi.3 Simulation results

vi.3.1 Improvement related to phases differences

In order to compare the improvement due to the combinations of phase differences we compare the results obtained using the squared visibility with: the closure phases, the differential phases and both phase differences.

Figure 6: Estimated 3–D object with differentials phases only.

The classical observable of phase differences used in OI images reconstructions is the closure phases. The first estimation uses only this phase information. The result is shown in Fig. 5. Due to the lack of information: phase closures for baselines per wavelength, the reconstruction performances are quite poor. The narrow objects at high wavelength are well reconstructed. The object structure (two disks) is not correctly found at lower wavelengths where artifacts are present.

The second simulation uses the differentials phases as the only phase differences measure. It leads to the estimated objects shown in Fig. 6. The object structures are better reconstructed. However some artifacts are present and the object surface are not as smooth as they should be.

Finally, the last simulation uses all the available phase information: phase closures and differentials phases. The result in Fig. 7 shows that the object is extremely well reconstructed at all wavelengths. The Relative Mean Squared Error (RMSE) of the reconstructed object () at wavelength defined as has been computed for each wavelength. The RMSEs in dB are presented in Fig. 5, 6 and 7 to quantify the improvement.

Figure 7: Synthetic data. Estimated 3–D object with closure phases and differentials phases.

vi.3.2 Resolved stars reconstruction

Images estimated per channel are shown in Fig. 9, values lower than of the maximum of the cube are thresholded. The shape and the size of the two stars are correctly restored. The polychromatic brightness distributions of the two estimated disks are compared to the solution in Fig. 9. Even if the estimated images are not as smooth as the original ones the estimated brightness follows correctly the true distribution.

Vii Acknowledgments

The present work was funded by the french ANR project POLCA (ANR-2010-BLAN-0511-02). One aim of POLCA is to elaborate dedicated algorithms for model-fitting and image reconstruction using polychromatic interferometric observations. The authors thank R. Petrov for proposing the use of the differential phase as a key item of this image reconstruction process.

Figure 8: Estimated Images per channels, the units are in arc second.
Figure 9: Original and estimated Brightness distribution overs channels for the two disks.
Figure 8: Estimated Images per channels, the units are in arc second.

Viii Conclusions

A polychromatic 3–D image reconstruction for optical interferometry namely PAINTER  is presented. The main contribution of PAINTER  is to uses both phase closure and differential phases information to estimates the original phases of the observed scene. Compared to the use of phases difference independently, the combination leads to use a transformation matrix with a “number of wavelengths/telescopes” dependent rank. The constraint optimization follows the ADMM methodology with spatial and spectral regularization. PAINTER  is able to deal with the complete chain of astronomical standard data (OIFITS file to data). Realistic simulation on noisy data are encouraging and shows the potential of PAINTER.


  • (1) C. Aime, H. Lantéri, M. Diet, and A. Carlotti. “Strategies for the deconvolution of hypertelescope images,” A&A, 543, A42, 1–9, 2012.
  • (2) A. Labeyrie. “Resolved imaging of extra-solar planets with future 10-100km optical interferometric arrays,” A&A Sup. Series, 118, 517–524, 1996.
  • (3) A. Labeyrie, D. Mourard, F. Allouche, R. Chakraborthy, J. Dejonghe, A. Surya, Y. Bresson, C. Aime, D. Mary and A. Carlotti. “Concept study of an “extremely Large hyper telescope” (eLhyt) with 1200m sparse aperture for direct imaging at 100 micro-arcsecond resolution,” In Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, July 2012.
  • (4) A. Labeyrie. “Interference fringes obtained on VEGA with two optical telescopes,” ApJ, 196:L71–L75, March 1975.
  • (5) J.D. Monnier. “Optical interferometry in astronomy,” Reports on Progress in Physics, 66(789-857), 2003.
  • (6) R.C. Jennison. “A phase sensitive interferometer technique for the measurement of the Fourier transforms of spatial brightness distribution of small angular extent,” MNRAS, 118(276-284), 1958.
  • (7) D. Mourard, I. Tallon-Bosc, F. Rigal, F. Vakili, D. Bonneau, F. Morand and P. Stee. “Estimation of visibility amplitude by optical long-baseline Michelson interferometry with large apertures,” Astron. Astrophy., 288, 675–682, August 1994.
  • (8) F. Roddier. “Triple correlation as a phase closure technique.” Optics Communications, 60:145–148, November 1986.
  • (9) F. Roddier and P. Lena. “Long-baseline Michelson interferometry with large ground-based telescopes operating at optical wavelengths. I - General formalism: Interferometry at visible wavelengths.” Journal of Optics, 15:171–182, August 1984.
  • (10) M. P. van Haarlem, M. W. Wise, A. W. Gunst, G. Heald, J. P. McKean, J. W. T. Hessels, A. G. de Bruyn, R. Nijboer, J. Swinbank, R. Fallows, M. Brentjens, A. Nelles, R. Beck, H. Falcke, R. Fender, J. Hörandel, L. V. E. Koopmans, G. Mann, G. Miley, H. Röttgering, B. W. Stappers, R. A. M. J. Wijers, S. Zaroubi, M. van den Akker, A. Alexov, J. Anderson, K. Anderson, A. van Ardenne, M. Arts, A. Asgekar, I. M. Avruch, F. Batejat, L. Bähren, M. E. Bell, M. R. Bell, I. van Bemmel, P. Bennema, M. J. Bentum, G. Bernardi, P. Best, L. Bïrzan, A. Bonafede, A.-J. Boonstra, R. Braun, J. Bregman, F. Breitling, R. H. van de Brink, J. Broderick, P. C. Broekema, W. N. Brouw, M. Br ggen, H. R. Butcher, W. van Cappellen, B. Ciardi, T. Coenen, J. Conway, A. Coolen, A. Corstanje, S. Damstra, O. Davies, A. T. Deller, R.-J. Dettmar, G. van Diepen, K. Dijkstra, P. Donker, A. Doorduin, J. Dromer, M. Drost, A. van Duin, J. Eislöffel, J. van Enst, C. Ferrari, W. Frieswijk, H. Gankema, M. A. Garrett, F. de Gasperin, M. Gerbers, E. de Geus, J.-M. Grie meier, T. Grit, P. Gruppen, J. P. Hamaker, T. Hassall, M. Hoeft, H. A. Holties, A. Horneffer, A. van der Horst, A. van Houwelingen, A. Huijgen, M. Iacobelli, H. Intema, N. Jackson, V. Jelic, A. de Jong, E. Juette, D. Kant, A. Karastergiou, A. Koers, H. Kollen, V. I. Kondratiev, E. Kooistra, Y. Koopman, A. Koster, M. Kuniyosh