`123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140` ``````MIRA: Multi-aperture Image Reconstruction Algorithm --------------------------------------------------- MIRA (Multi-aperture Image Reconstruction Algorithm) is an algorithm for image reconstruction from interferometric data. Copyright (c) 2001-2009, Eric Thi�baut Quick instructions for using MIRA --------------------------------- You must have Yorick (at least version 2.1), Yeti (version 6.2.0 http://www-obs.univ-lyon1.fr/labo/perso/eric.thiebaut/yeti.html) and OptimPack (version 1.2) installed. 1. launch Yorick (note that '\$' is the shell prompt and '>' is Yorick prompt): \$ yorick 2. load "mira.i" (this should also load Yeti plugin): > include, "mira.i"; 3. 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"); 4. 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"; 5. Choose a suitable regularization method: > rgl = rgl_new("smoothness"); 6. 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 one 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). 6. 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). 7. 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); Notes: a. *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). b. 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 strictly positive value for XMIN (choose a very small value, for instance: 1e-50*NRM/NPIX where NRM is the normalization level and NPIX the total number of pixels). c. Likewise, NORMALIZATION=1 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). d. The cost function 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); > ... e. It is usually better to work with a purposely too high regularization level. And to lower the value of MU as the image reconstruction converges. Guidelines ---------- The pixel size must not be too small: pixel size < (1/2) wavelength / max(baseline) The field of view must not be too small otherwise: - oversmoothing of the Fourier spectrum - risk of field aliasing field of view must not be too large otherwise: - undersmoothing of the Fourier spectrum w.r.t. to the (u,v) sampling and, for real data, to the (u,v) smoothing due to the finite size of the telescopes and to the tapering (spatial filtering by fiber optics). Typically: FOV ~ a few times (wavelength / telescope diameter) ``````