MiRA image synthesis algorithm

MiRA is an algorithm for image reconstruction from data provided by optical interferometers. I wrote this algorithm to cope with problems specific to this type of data (sparseness of the measurements, missing of Fourier phase information, etc.) and not solved by existing methods, e.g., developed for radio-astronomy. MiRA won the 2008' Interferometric Imaging Beauty Contest [6] organized by IAU to compare the image synthesis algorithms designed for optical interferometry. An overview of the existing image synthesis algorithms in radio and optical interferometry is provided in [3]; more specifically, MiRA algorithm has been briefly described in [1, 5] and MiRA has been compared to Wisard method in [4].

In a nutshell, MiRA proceeds by direct minimization of a penalized likelihood. This penalty is the sum of two terms: a likelihood term (typically a χ²) which enforces agreement of the model with the data, plus a regularization term to account for priors. The priors are required to lever the many degeneracies due to the sparseness of the spatial frequency sampling. MiRA implements many different regularizations (quadratic or edge-preserving smoothness, total variation, maximum entropy, etc.) and let the user defines his own priors. The likelihood penalty is modular and designed to account for available data of any kind (complex visibilities, powerspectra and/or closure phase). One of the strength of MiRA is that it is purely based on an inverse problem approach and can therefore cope with incomplete data set; for instance, MiRA can build an image without any Fourier phase information.

MiRA is a research tool, it is however free software usable by others (though you need some basic knowledge in optical interferometry). For instance, MiRA has been used to teach image reconstruction in interferometry to the students of the 2008' VLTI Summer School held in Keszthely [3]. MiRA has been used to obtained the results published in a number of scientific journals [2, 7, 8] or conferences.

For now, input data must be in OI-FITS format, extension to other formats (e.g. UV-tables) is on the way. The algorithm is currently modified to achieve true 3-D reconstruction of multi-wavelength images and to account for other interferometers such as ALMA and next generation VLTI instruments (MATIS).

Installation

To use MiRA, you must have installed recent versions of Yorick (version ≥ 2.2 or CVS version), Yeti (version ≥ 6.3.1) and OptimPack (version ≥ 1.3.2).

Installing MiRA is as simple as unpacking the source archive anywhere. There is no other installation steps: the software has been packaged so that it is directly usable from its directory. You can however change the name of this directory to match your taste. See here in case of problems.

Launch Yorick (note that, in code samples below, $ is the shell prompt and > is Yorick prompt):
```
$ yorick
```
Load file mira.i (this should also load Yeti plugin, you may have to specify the full path to file mira.i):
```
> include, "mira.i";
```
Load OI-FITS data file (db will be your MiRA instance for this data file; in the data file, you must choose one with keyword eff_wave or choose a spectral range with keywords eff_wave and eff_band):
```
> db = mira_new("data/data1.oifits");
```
Configure data instance for image reconstruction parameters (keyword dim is the number of pixels along the width and height of the restored image; keyword pixelsize is the size of the pixel in radians; keyword xform is the name of the method to approximate the Fourier transform, can be "exact" or "fft", default is "exact"):
```
> mira_config, db, dim=50, pixelsize=0.5*MIRA_MILLIARCSECOND,
               xform="exact";
```
Choose a suitable regularization method, for instance:
```
> rgl = rgl_new("smoothness");
```
Attempt an image reconstruction (from scratch):
```
> img1 = mira_solve(db, maxeval=500, verb=1, xmin=0.0, normalization=1,
                    regul=rgl, mu=1e6);
```
where img1 is the output image, db is the data instance, maxeval is the maximum number of evaluations of the cost function, verb is set to 1 to display convergence information at every successful step (verb=10 to display this information every 10 steps and verb=0 or nil to display no information), xmin=0 to enforce a positivity constraint (xmin is the minimum allowed value in the restored image, you certainly do not to want to omit this option), regul set the regularization method, mu is the global weight of regularization (the higher its value the smoother the result), ftol controls the stopping criterion of OptimPack1 (which to see).
You can also try a reconstruction given an initial image img0:
```
> img2 = mira_solve(db, img0, maxeval=500, verb=1, xmin=0.0,
                    normalization=1, regul=rgl, mu=1e6);
```
Note that if img0 is not of size dim×dim it will be resampled to that size (assuming the field of view is the same).

With keyword zap_data, you can just build the default image as assumed by the regularization:

> img0 = mira_solve(db, maxeval=500, verb=1, xmin=0.0, normalization=1,
                    regul=rgl, mu=1e6, zap_data=1);

Troubleshooting

In MiRA's top directory you'll find different Yorick source files, among others some utilities that I wrote but which may be distributed elsewhere. I try to always distribute the last versions of these files with MiRA so beware of possible conflicts if you have installed older versions in one of the paths searched by Yorick.
All units are in SI (i.e. angles are in radians, wavelengths in meters, etc). A very common error in parameter settings is to use completely out of range values because you assume the wrong units. To make things a little bit easier, some constants are pre-defined by MiRA package, and you can use for instance: 3*MIRA_MILLIARCSECOND (instead of 1.45444e-08 radians) or 5*MIRA_MICRON (instead of 5e-6 meters).
Positivity is a must, you cannot expect a good image reconstruction (at least with any practical u-v coverage) without option xmin=0 in mira_solve. If you use certain regularizations such as the entropy, you must specify a non-negative value xmin (you may also set the epsilon attribute of the regularizer to a very small value, for instance: 1e-20*NRM/NPIX where NRM is the normalization level and NPIX the total number of pixels).
Likewise, normalization=1.0 must not be omitted if your data set obeys OI-FITS standard (i.e. visibilities are normalized) and has no explicit measurement at frequency (0,0).

The cost function which is minimized by MiRA is highly non-quadratic and may cause difficulties for OptimPack to converge. To overcome this, it is sometimes very effective to simply restart mira_solve with the an initial image provided by a previous reconstruction (possibly after recentering by mira_recenter). For instance:

> img2 = mira_solve(db, img0, maxeval=500, verb=1, xmin=0.0,
                    normalization=1, regul=rgl, mu=1e6);
> img2 = mira_recenter(img2);
> img2 = mira_solve(db, img2, maxeval=500, verb=1, xmin=0.0,
                    normalization=1, regul=rgl, mu=1e6);
> img2 = mira_recenter(img2);
> img2 = mira_solve(db, img2, maxeval=500, verb=1, xmin=0.0,
                    normalization=1, regul=rgl, mu=1e6);
> ...

It is usually better to start with a regularization level purposely too high; and to lower the value of mu as the image reconstruction converges.

References

[1]	E. Thiébaut, P.J.V. Garcia & R. Foy, "Imaging with Amber/VLTI: the case of microjets," Astrophysics and Space Science 286, pp. 171—176, 2003.
[2]	S. Lacour, S. Meimon, É. Thiébaut, G. Perrin, T. Verhoelst, E. Pedretti, P.A. Schuller, L. Mugnier, J. Monnier, J.P. Berger, X. Haubois, A. Poncelet, G. Le Besnerais, K. Eriksson, R. Millan-Gabet, M. Lacasse, & W. Traub, "The limb darkened Arcturus Imaging with the IOTA/IONIC interferometer," Astronomy & Astrophysics 485, 561—570, 2008.
[3]	É. Thiébaut, "Image Reconstruction with Optical Interferometers," in VLTI Summer School: "Astrometry and imaging with the very large telescope interferometer," Keszthely (Hungary), 2008.
[4]	G. Le Besnerais, S. Lacour, L. M. Mugnier, É. Thiébaut, G. Perrin & S. Meimon, "Advanced imaging methods for long-baseline optical interferometry," IEEE Journal of Selected Topics in Signal Processing 2, 767—780, 2008.
[5]	É. Thiébaut, "MiRA: an effective imaging algorithm for optical interferometry," in Astronomical Telescopes and Instrumentation, Proc. SPIE 7013, 70131I (Marseille, France), 2008.
[6]	W. Cotton, J. Monnier, F. Baron, K.-H. Hofmann, S. Kraus, G. Weigelt, S. Rengaswamy, É. Thiébaut, P. Lawson, W. Jaffe, C. Hummel, T. Pauls, H. Schmitt, P. Tuthill & J. Young, "2008 imaging beauty contest," in Astronomical Telescopes and Instrumentation, Proc. SPIE 7013, 70131N (Marseille, France), 2008.
[7]	J.-B. Le Bouquin, S. Lacour, S. Renard, É. Thiébaut & A. Merand, 2009, "Pre-maximum spectro-imaging of the Mira star T Lep with AMBER/VLTI," Astronomy & Astrophysics 496, L1—L4, 2009.
[8]	F. Millour, O. Chesneau, M. B. Fernandes, A. Meilland, G. Mars, C. Benoist, É. Thiébaut, P. Stee, F. Baron, J. Young, A. Carciofi, A. D. de Souza, P. Kervella, R. G. Petrov, F. Vakili, K.-H. Hofmann, P. Bendjoya, L. Waters & G. Weigelt, 2009, "A binary engine fueling HD 87643's complex circumstellar environment using AMBER/VLTI imaging," to be published in Astronomy & Astrophysics.

MiRA – a Multi-aperture Image Reconstruction Algorithm

Links

Downloads

Installation

Troubleshooting

References