Technical Report - WISARD software documentation
Contents
1. Moreover with Eqs 4 12 and 4 11 we have ata A 1 res 2 1 2 J x a t t 2 Wrad ij t a t t g Utanjj t Wien ij 2 a t t 1 res igdata 1 Yos _jgdata Wrad ij t yee Ui x a t t e g 5 tan iji m yi x a t t e 95 H We gather I e a0 t O Re if a a t ij 2 F e e a a t t Re ui a t t 4 13 W22 ij t res with for clarity we omit here the j and t indexes Wil Wrad cos gu SE Wtan sin W12 Wrad Wtan sin cos 6 4 14 99 Wrad sin oe Wtan COS giis It is from this expression that the criterion and its gradient w r t the object are computed in WISARD L MUGNIER amp S MEIMON UNCLASSIFIED 22 SANS MENTION MAY 2007 DE PROTECTION APPENDIX I THE CLOSURE AND BASELINE OPERATORS C AND B Let N be the number of telescopes of the interferometric array We have the following definitions 1 1 1 1 A 1 4 1 Id SA Id By Bo 1 3 Cy Bn Id L4 ct oT oct L 5 for N gt 3 It 15 easy to see that C B 0 ONERA EMEN THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 23 SANS MENTION MAY 2007 DE PROTECTION APPENDIX II SQUARE ROOT OF A GAUSSIAN DISTRIBUTION Let us suppose we measure the squared value s of a positive value a with an additive noise system gdata a graise 3106 being 0 mean2. ra t r t u t r t a t us4 t ra t rs t or in a matrix formulation u t B r t 4 3 where the N x N operator B is the baseline operator see appendix I page 22 4 3 1 Ideal interferometric data We consider a monochromatic source with a wavelength and denote z the brightness distribution being angular coordinates For each baseline it is possible to access to the following short exposure data aphase t a modulus a t The complex vector y t adet 1 ei O is called the complex visibility vector According to the Van Cittert Zernike theorem 12 complex visibilities are ideally linked to the Fourier Transform FT of x at the 2D spatial frequency A ult u t PX 4 4 through yiz t FT z vij t i e l py pyl t arg FT z 6 wlt m a t si s t FT 2 v t or in a matrix formulation olw t T ad t a x t In what follows we will only consider a discretized version z of The Fourier Transform then corresponds to a matrix product by a Fourier operator H y z t H t x 4 7 For convenience we also define a t arg H t z a z t H t z 4 8 s x t H t x The operators and arg are point to point operators ONERA ee r E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 17 SANS MENTION MAY 2007 DE PROTECTION
3. Ref 1 is available on the French national multidisciplinary thesis server TEL at the follow ing address http tel archives ouvertes fr docs 00 05 44 98 PDF these finale pdf The two other ref erences are also available on line at http laurent mugnier free fr publis Meimon JOSAA 05 pdf and http laurent mugnier free fr publis Meimon OL 05 pdf respectively ONERA ee m E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON MAY 2007 UNCLASSIFIED SANS MENTION DE PROTECTION 2 INSTALLATION OF THE WISARD SOFTWARE 2 1 Requirements The WISARD software should run on all platforms such that a IDL or GDL is installed and b the OptimPack optimization package of ric Thi baut see Section 2 3 is installed This means that WISARD can be installed and run on at least all Unix like platforms 2 2 Unpacking WISARD is distributed in a compressed archive called Wisard date tar gz Unpacking it is straightforward and will create a WISARD directory under which is the distribution From the shell prompt type tar xvfz Wisard date tar gz The directories created under W1 doc lib lib wisardlib lib extralib lib oneralib lib optimpacklib optimpack pro inputdata SARD and the nature of their contents are the following This documentation of WI
4. WISARD lib PATH If you have unpacked the archive elsewhere simply replace with the appropriate directory in the command above See the example of Section 3 3 for using a relative path instead of an absolute one 2 5 Acknowledgements The authors of WISARD hereby express their gratitude to the following people as WISARD uses some routines of theirs Eric Thi baut from CRAL WISARD heavily uses his neat OptimPack software as its minimization engine Fr d ric Cassaing and Jean Marc Conan WISARD uses a few routines by these ONERA scientists These routines and some others by Laurent Mugnier are distributed along with WISARD with these authors permissions in the 1ib oneralib directory John D Monnier for his OIFITS reading routines as well as the authors and maintainers of the NASA astronomical library for their routines used by John s ones The OIFITS reading routines are avail able at http www astro lsa umich edu monnier while NASA s astronomical library can be found at http Adlastro gsfc nasa gov For convenience their routines are included in the 1ib extralib di rectory of the WISARD distribution ONERA ee r E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 8 SANS MENTION MAY 2007 DE PROTECTION 3 WISARD S USER MANUAL 3 1 General presentation The WISARD reconst
5. 4 3 2 Data model The long exposure observables considered in this report are squared moduli sdata z and closure phases 39 1 These observable are affected by a noise and we suppose in this report that only second degree statistics are provided and that the closure phase measurement set and the squared modulus measurement set are two uncorrelated sets of variables The noise distribution we can make out of these statistics is then gaussian s t a ag t s t s 0 N 0 Raw Bn C z t Bee noise t A N 0 Row The algebraic structure of the closure operator C is described in appendix I page 22 4 9 4 3 3 Input Data format in WISARD The data format is a vector of structures one structure for each time of measurement with the following fields VIS2 the squared visibilities for each baseline VIS2ERR the standard deviation on the squared visibilities CLOT the closure phases for each triplet of telescopes involving telescope 1 CLOTERR the standard deviation on the closure phases FREQS 0 first coordinate of the spatial frequencies involving telescope 1 FREQS V second coordinate of the spatial frequencies involving telescope 1 For example at each moment t for a 4 telescope array the fields contain field notation VIS2 sif 0 sig 0 a menja el spet VIS2ERR G so D C si3 D Osalt O s s t Cann Ossa t CLOT
6. I l I l 0 2 0107 4040 6040 8010 1 0410 1 240 frequency Figure 3 2 Typical plot displayed by WISARD at convergence see text ONERA ee m E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 4 SANS MENTION MAY 2007 DE PROTECTION 4 WISARD DESIGN REPORT 4 1 Introduction Optical interferometry allows one to reach the angular resolution that a hundred meter telescope would provide using several ten meter telescopes The interferograms of current instruments are affected by turbulence which corrupts the recorded object phases so one is led to form quantities that are turbulence independent such as squared visibilities and phase closures In order to cope with the missing phase information we introduce phase calibration parameters to be estimated jointly with the observed object and we propose WISARD for Weak phase Interferometric Sample Alternating Reconstruction Device to perform this estimation This algorithm combines within a Bayesian framework an alternating estimation of the object and phase parameters in the spirit of self calibration algo rithms proposed by radio astronomers 11 a recently developped noise model suited to optical interferometry data 2 and an edge preserving regularization 9 to deal with the sparsity of the data typical of optical inter ferometry 4 2 Structure of WISARD WISARD is made of the following major blocks
7. 1 Introduction c2 ao om om mo ok RU R xk Ee 14 4 2 Structure of WISARD llle 14 4 3 Interferometricobservables 15 4 3 1 Ideal interferometric data 16 4 3 2 Datamodel 17 4 3 3 Input Data format in WISARD 00002 eee 17 4 4 From input data toa myopic model 17 4 4 1 Visibility modulus pseudo data ll 18 4 4 2 Visibility phase pseudo data 18 4 5 From a myopic model to a myopic convexified model 19 4 6 The criterion minimized gs ks ck ka 93593 KD ORS OK SE KR EK xx 20 4 6 1 Expression for the aberration step 20 4 6 2 Expression for the object step 21 APPENDIX I THE CLOSURE AND BASELINE OPERATORS C AND B 22 APPENDIX II SQUARE ROOT OF A GAUSSIAN DISTRIBUTION 23 APPENDIX III CARTESIAN GAUSSIAN APPROXIMATION TO A POLAR GAUSSIAN DISTRIBUTION OR REO 49 CORRI Roe LR es Eg 25 HII General expression oom o R9 m RR xx RA s 25 ONERA we r S THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 3 SANS MENTION MAY 2007 DE PROTECTION IIL2 Gaussian Approximatio
8. 1 e 2 data Oa I lt Utanjj t oan 20 a 45 t i l e P d 1 e IM data 4 2 C ais t I Qij These equations are explained in appendix III 4 6 The criterion minimized We define EIUS 2 Re CER i hs Vz C Cf eq 4 10 Ztanaj t Sm Vz C Cf eq 4 10 4 11 res A data iBa y z a t t y t y a t e Po and propose to use the following data likelihood criterion gre a a n a t px a t t 1 res 2 5 Utan j t aa 4 a a t t 4 12 A IE a alt t waaut alt NI 4 6 1 Expression for the aberration step With Eqs 4 11 we have y z a t t y t y z t e P9 Before starting the aberration step we compute all the elements useful for the minimzation independant of a ie y z t a z t and d x t To compute the gradient in we also compute two other quantities called wx and wx2 in the codes All these quantities are arranged in the structure weights x which does not depend on the a ONERA ee m E THE FRENCH AEROSPACE LAB 2 L MUGNIER amp S MEIMON UNCLASSIFIED SANS MENTION MAY 2007 DE PROTECTION 4 6 2 Expression for the object step With Eqs 4 11 and 4 8 we gather y z a t t y t H t aDiag es y t Diag ee eu H t x PH t
9. 3 us t Boi t Bisa t CLOTERR 3 O Bia t O Biza t O Pisa t FREQS U 3 uio t ua t FREQS V 3 Vy12 t vas t vaa t The covariance matrices R and Agi are implicitely supposed diagonal 4 4 From input data to a myopic model This section describes the code wisard data2mdata pro We want to recast the data model in a myopic one with visibility phase and modulus pseudo data The corresponding format is described below The myopic data format is a vector of structures one structure for each time of measurement with the following fields ONERA ee m E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 18 SANS MENTION MAY 2007 DE PROTECTION ISAMP the visibility moduli for each baseline ISAMPERR the standard deviation on the visibility moduli SPHI the visibility phases for each baseline SPHIERR the standard deviation on the visibility phases FREQS 0 first coordinate of every spatial frequency FREQS V second coordinate of every spatial frequency V V V V For example at each moment t for a 4 telescope array the fields contain field notation VISAMP 6 and ds aalt quie Goat one VISAMPERR 2 6 Galt Tais t Caralt Cont s Cut Cau VISEHI 6 ix t di3 t Ai 0 5 03 t VISPHIERR 6 O 12 t V oia t dia t 0 os t O 24 t 2 osa t FRE
10. AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 7 SANS MENTION MAY 2007 DE PROTECTION 2 3 Dependencies WISARD heavily uses ric Thi baut s neat OptimPack software as its minimization engine So WISARD needs the OptimPack library for IDL Opt imPack_IDL so and its IDL frontend op pro files to be installed Although OptimPack is not part of WISARD we have included the OptimPack IDL frontend and the pre compiled OptimPack library for IDL for Linux and Solaris both 32 and 64 bit architectures in the 11b optimpacklib directory of the WISARD distribution with their author s permission So if your computer is a Sun running Solaris or a PC running Linux you should have nothing to do Should these pre compiled librairies not work for you we have also included all the necessary sources c handMakefile underthe opt impack directory A make command from the opt impack idl directory should build the library for your architecture Then simply copy the resulting Opt imPack_IDL so files in the 1ib optimpacklib directory 2 4 Installing WISARD Once you have unpacked the archive Sect 2 2 and once you have compiled the OptimPack library for IDL if necessary Sect 2 3 the only installation to be done is to make all routines under the 1ib directory known to IDL Assuming that you unpacked the archive in your home directory this is achieved by typing the following at the JDL prompt IPATH expand path
11. a factor 10 or more i e assume a stronger object If you set PSD 100D psd for instance you will get an under regularized solution very similar to that obtained without a regularization quite noisy but with the emergence of the central peak present in the original object Let s now try the white L1 L2 spike preserving regularization Set WHITE 1 and choose SCALE of the order of the object s mean value this is easy because the reconstructed object is of unit sum For DELTA 1 one obtains a quadratic L2 regularization and a reconstruction similar to a PSD based one For the desired spike preserving regularization you can set white 1B scale 1D NP min 2 for white L1 L2 scale avg object level 1 NP 2 delta 1D Run it x lll2white WISARD data FOV fov NP MIN np min GUESS 0 THRESHOLD threshold SCALE scale DELTA delta WHITE white AUX OUTPUT aux output lll2white DISPLAY The result is much smoother than an under regularized solution and has the central peak of the object as visible as the nose in the middle of one s face French expression Figure 3 1 shows three reconstructions an under regularized PSD based reconstruction with the above PSD multipied by 100 the PSD based reconstruction run above and the white L1 L2 reconstruction just obtained Figure 3 1 Three reconstructions obtained with WISARD under r
12. regularization of the reconstruction MEAN 0 15 used both for PSD regularization and for white L1 L2 regularization 1 e when DELTA SCALE and WHITE are set scalar factor for L1 L2 regularization used to set the thresh old between quadratic L2 and linear L1 regularization See routine J PRIOR L1L2 and the example file for some more details See 9 for a complete description of the Ll1 L2 regularization The PDF of this paper is on line at http laurent mugnier free fr publis Mugnier JOSA A 04 pdf scalar SCALE factor for L1 L2 regularization should be of the order of the average object value if WHITE is set and of the or der of the RMS object s gradient value if WHITE is not set See PRIOR 1112 and example file for some more details See the paper cited above for the whole story flag to switch between edge preserving regularization WHITE O de fault and spike preserving regularization WHITE 1 In the lat ter case the regularization is performed independently on each pixel value hence the flag name full path of the OptimPack library see FMIN OP for details if nec essary Under Unix systems the OptimPack library should be found automatically structure containing various optional AUXiliary outputs for diagnos tic purposes The details of this structure is given in the on line docu mentation for the wisard routine display the reconstructed object and the fit to the
13. visibilities along the way prints version number before execution prints the on line documentation and exits prints information about copyright and exits An example of use of WISARD for several regularizations types and levels is the wisard_batch pro file located in the pro directory of the distribution Studying this example running it and modifying its pa rameters should be a good starting point for building your own reconstruction know how The data used here from the Imaging Beauty Contest 2004 organized by Peter Lawson for the IAU 4 With the software distributed here you will obtain reconstructions notably better than the ones we obtained with WISARD at the time 10 and at least as good as the best one of 4 Below are described the essential portions of this file First change directory to the pro directory of the WISARD distribution and run the IDL interpreter ONERA ee r E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 1 SANS MENTION MAY 2007 DE PROTECTION If for instance WISARD has been unpacked in your home directory it is under the WISARD directory the commands to be issued at the shell prompt are cd WISARD pro idl Then at the IDL prompt add the libraries needed by WISARD to IDL s search PATH IPATH expand path lib PATH Read the license dummy wisard copyright If you accept
14. AY 2007 DE PROTECTION Elliptic gaussian approximation Noise statistics O Figure III 1 Polar gaussian distribution contour plot and its cartesian gaussian approximation ONERA ee THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 28 SANS MENTION MAY 2007 DE PROTECTION 1 2 3 4 5 6 7 8 BIBLIOGRAPHY S Meimon Reconstruction d images astronomiques en interf rom trie optique PhD thesis Universit Paris Sud 2005 1 S Meimon L M Mugnier and G Le Besnerais A convex approximation of the likelihood in optical interferometry J Opt Soc Am A November 2005 1 4 1 III 1 III 2 S Meimon L M Mugnier and G Le Besnerais Reconstruction method for weak phase optical interferometry Opt Lett 30 14 pp 1809 1811 July 2005 1 P R Lawson W D Cotton C A Hummel J D Monnier M Zhao J S Young H Thorsteinsson S C Meimon L Mugnier G Le Besnerais E Thi baut and P G Tuthill An interferometric imaging beauty contest In New frontiers in stellar interferometry edited by W A Traub vol 5491 pp 886 899 Proc Soc Photo Opt Instrum Eng 2004 Conference date June 2004 Glasgow UK 113 3 G Demoment Image Reconstruction and Restoration Overview of Common Estimation Structures and Problems IEEE Trans Acoust Speech Sig
15. CH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 25 SANS MENTION MAY 2007 DE PROTECTION APPENDIX III CARTESIAN GAUSSIAN APPROXIMATION TO A POLAR GAUSSIAN DISTRIBUTION HII General expression With consider a polar distribution of a gaussian vector y of modulus a and phase go a oe 11 1 qata aras 2 where and a are 0 mean real gaussian vectors of covariance matrices Ra and Rg the vectors and a are supposed uncorrelated With the definitions gt aexpid noise A data P noise i x i j 113 Gm Lyme Yrad n Utan y lI Se n Yrad ll n Utan noise gt we gather Yad la a cos a I4 Yon la De qose sin pre A complex vector is gaussian if and only if each of its components is gaussian A complex is gaussian if and only if in any cartesian basis its two components are gaussian So y is gaussian if and only if y is gaussian which is not the case 2 In what follows we show how to optimally approximate the distribution of by a gaussian distribution IIL 2 Gaussian Approximation We caracterize our Cartesian additive gaussian approximation i e its mean yr and covariance Rise by minimizing the Kullback Leibler distance between the two noise distributions which gives 2 noise Yrad Uraa F a Uis Uian 311 n n n T FU poise E Grad B
16. ONERA THE FRENCH AEROSPACE LAB DEPARTEMENT OPTIQUE THEORIQUE ET APPLIQUEE Technical Report WISARD software documentation L MUGNIER amp S MEIMON May 2007 L MUGNIER amp S MEIMON MAY 2007 UNCLASSIFIED SANS MENTION DE PROTECTION CONTENTS ABBREVIATIONS AND ACRONYMS les 4 1 INTRODUCTION SESE RSEESERS CER OES PRE HES ex ES 5 2 INSTALLATION OF THE WISARD SOFTWARE 6 2 1 Re uirements ode ee eh ek Oe eS KO OSE apos EUH e 6 22 Unpacking 6 2 3 Dep ndenti s 2 om s 369 q Q aeg ep a ee ep de d Reap de densi les EEK ats 7 2 4 Installing WISARD 24446 46 48 489 ee s ee ss s es s RR s 7 2 5 AcknowledgemenfS _____ __ _ aaa e 7 3 WISARD SUSERMANUAL 8 3 1 General presentation RR RR ERE RR E A A 8 3 2 List of parameters and accepted keywords 8 3 3 Example of use of WISARD 10 4 WISARD DESIGN REPORT 14 4
17. PhD thesis work of S Meimon 1 and in refer ences 213 The WISARD software described here is a complete rewrite from scratch by S Meimon and L Mug nier of earlier programs of the aforementioned thesis This rewrite has bought us a dramatic gain in speed and in code maintainability as well as the elimination of quite a few bugs Yet there may of course remain some This code and its documentation are distributed in the hope that it will be useful but on an as is basis without any warranty of any kind More precisely WISARD is a free software licensed under the CeCILL B license version 1 which can be found at http www cecill info WISARD is written in the commercial language IDL of ITT Visual Information Solutions http www ittvis com idl It should thus also run with GDL the GNU Data Language available at http gnudatalanguage sourceforge net Chapter 2 contains the instruction for installing WISARD on your computer Chapter 3 documents the use of WISARD and contains an example of reconstruction from interferomet ric data The data comes from the Imaging Beauty Contest 2004 organized by Peter Lawson for the IAU 4 The example batch file is part of the distribution and can be used and modified for self study Chapter 4 contains an in depth description of the method implemented by the WISARD code For more details on the method the interested reader should consult 1 2 3
18. QS U 6 Ur2 t U13 t 14 t U23 t U24 t 34 t FREQS V 6 valt vis t ui t voa t voa t Vga t 4 4 1 Visibility modulus pseudo data Each ageta t Casj p i is computed from input data sc t Tsy acording to in every case 4 sde e t os if s D gt then o Ree ij M sig t e See appendix II for justification of these formulae else oat 4 4 2 Visibility phase pseudo data data t containing the fata t and the vector og t containing nni i t iio For each instant t the vector the o4 is computed from the vector 8 containing the input data c5 acording to EN p t C gasta q e t 3 CiDiagIoa t Ch Id containing the input data 71 and the vector o g t ONERA ee r E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 190 SANS MENTION MAY 2007 DE PROTECTION where denotes a term to term product and C is defined in appendix I Let R be a diagonal matrix the vector of the diagonal components being r For any matrix M with compatible dimensions we have DURM OH SK m ij x R A uMi Y k So we can write c t 3 C C og t This is how the computation of a t is coded in wisard data2mdata pro 4 5 From a myopic model to a myopic convexified model This section describes the code wisard mdata2cmdata pro The
19. SARD WISARD its sub routines and routines used by them the WISARD main routine and its sub routines These routines are under the CeCILL B license Several extra routines quick and dirty code to input some OIFITS data wisard oifits2data pro which needs some work anyone interested plus routines used by the former code to read the OIFITS format by J D Mon nier see http www astro Isa umich edu monnier plus some routines from NASA s astro library used in turn by the former rou tines available from http idlastro gsfc nasa gov ONERA routines necessary to WISARD These routines are nof part of WISARD but used and distributed by permission of their respective authors These routines are under the CeCILL C license OptimPack library for IDL Opt imPack_IDL so and its IDL frontend routines op pro by Thi baut This library is not part of WISARD but is used by WISARD and distributed here by special permission of the author OptimPack sources c h and Makefile see Section 2 3 plus a copy of the IDL frontend routines op pro Once again this library is not part of WISARD but is used by WISARD and distributed here by special permission of the author These routines are under the GPL example batch file for WISARD A good starting point for self study cf Sec tion 3 3 input data for the example batch file ONERA ee r E THE FRENCH
20. Yrad Yrad Mi Yrad Ytan Ytan Ytan Ytan 115 ONERA lt THE FRENCH AEROSPACE LAB 26 L MUGNIER amp S MEIMON UNCLASSIFIED SANS MENTION MAY 2007 DE PROTECTION and we define A Rrad rad ssa tan Renoise T i IU usn Rian tan For a 0 mean gaussian vector of covariance matrix Rg E sin dp 0 Ro E noise cos 0 exp E sin sin 0j sinh R exp 1II 6 Rg Ro E cos 9j cos 69 cosh Rpp exp ae E cos gore sin 0 By combining Eqs 1II 5 III 3 III 4 and 1II 6 we obtain Roi B yha oa e 1 E 7 0 Rrad rad lai cosh Rg 1 Ra cosh Rg e Rrad tan 0 Fete teen ia Ga Rez sinh Bs ig 1J Rg Rg Re 2 IIL3 The scalar case Now we make the additionnal asumption that both and an ise are decorrelated i e Ra Diag 102 Rg Diag 10 We obtain Fara Diag o uud Riantan Diag wat Rbrad tan 0 with 72 2 2 El 1 nut 1 2 2 III 8 2 1 32 Cai 1 2o2 Cun bi pi duco et mue In this case we can plot for one complex visibility the true noise distribution i e a gaussian noise in phase and modulus and our gaussian approximation see fig III 1 L MUGNIER amp S MEIMON UNCLASSIFIED 27 SANS MENTION M
21. afirst block wisard data2mdata pro computes so called myopic complex pseudo data from the input data and error bars on these pseudo data These complex data have the following properties the phases measurements and error bars are computed from and compatible with the measured phase closures and their error bars the amplitudes measurements and error bars are computed from and compatible with the measured squared visibilities and their error bars aconvexification block wisard mdata2cmdata pro computes a gaussian approximation of the pseudo visibility data model This approximation is optimal in the sense of a Kullback Leibler distance see III athird block wisard set regul pro proposes an adapted prior term to be chosen among a list of regularizations such as positivity Power Spectral Density quadratic regularization linear quadratic also called L L2 regularization etc The algorithmic structure of WISARD make it easy to add other kinds of prior terms to the code a self calibration block performs the minimization It alternates a minimization of the criterion w r t the object for given aberrations which consists essentially in a minimization of a convex criterion under positivity constraint anoptimization of the aberrations for the current object It is accelerated by performing in parallel several optimizations of a subset of the aberrations instead of one global optimization The pattern of WISARD is
22. convexified myopic data format is a vector of structures one structure for each time of measure ment with the following fields VIS the complex visibilities for each baseline W RAD radial weight on the visibility moduli to be used in criterion see sect 4 6 W TAN tangential weight on the visibility moduli to be used in criterion see sect 4 6 FREQS U first coordinate of every spatial frequency FREQS V second coordinate of every spatial frequency For example at each moment t for a 4 telescope array the fields contain field notation Ves 6 yi3 t vis t via t uos t vba E vs 6 W RAD 6 Wrad 12 t Wrad 13 t Wrad 14 t Wrad 23 t Wrad 24 t Wrad 34 t W_TAN 6 Wtan 12 z Wtan 13 t tan 14 t Wtan 23 t Wtan 24 t Wtan 34 t FREQS U 6 uj t uia t 114 t 123 t 1124 t ua4 t FREQS V 6 vi t U13 t vi4 t U23 t 094 t v3a t Each Tuna Wean ig e is computed from myopic data ede t C p a agPte t ONERA ee r E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 20 SANS MENTION MAY 2007 DE PROTECTION acording to ata A ata ata ata age a t 10 t 4 10 data ata M yt asp t in o2 data data data 229g gt yj t aj t exp oi t 2 e 2 1 Wrad ij t loa g2 2 205 0 1 e suo
23. d For a disk like object you should take a PSD with a 1 f dependence on the spatial frequency as the square of the FT of a disk has an envelope that goes to zero with this dependence You should threshold this PSD so that it remains finite at 0 and so that its dynamic range remains reasonable as a large dynamic range for the PSD may hinder the minimization distance double shift dist NP min NP min 2 NP min 2 PSD 1D distance 3 gt 1D lt 1d6 Run WISARD displaying both the object being reconstructed in window 0 and the fit of the recon structed visibilities to the data in window x psd WISARD data FOV fov NP MIN np min GUESS guess THRESHOLD threshold PSD psd AUX OUTPUT aux_output_psd DISPLAY ONERA lt eee THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 12 SANS MENTION MAY 2007 DE PROTECTION The obtained reconstruction shows the correct overall shape of the reconstructed object on a dark and quite noiseless background but without the central peak of the original object see 4 A global factor multiplying the PSD acts as the inverse of a regularization parameter In other words for more regularization to get a smoother object you should divide the PSD by e g a factor 10 or more i e assume a weaker object and for less regularization to get an object with more details you should multiply the PSD by e g
24. described in Fig 4 1 ONERA ee r E THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 15 SANS MENTION MAY 2007 DE PROTECTION Raw data y Recasting i Myopic pseudo data Convexification i Myopic approx data Initialization prior object guess a 0 Self calibration i Object step T Aberration step Reconstruction Figure 4 1 WISARD algorithm loop 4 3 Interferometric observables Consider a model 2 telescope interferometer which two identical apertures 1 T gt are located at three space positions On and OTa We denote P the plane normal to the pointing direction and r and r the projection of OT and OF onto P The baseline wu is defined by the projection of the displacement Ti onto P A Uiz T2 Ti Because of the Earth rotation the pointing direction changes during an observing night so the baseline is time dependant too ui2 t ra t r t 4 1 Let us consider now a N telescope interferometer There are as many baselines as ways to choose 2 telescopes among Ni N N N 1 Ny a lt 4 2 b 2 5 4 2 ONERA a THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 16 SANS MENTION MAY 2007 DE PROTECTION For N 4 the N 6 baselines are uis t ro t r t uia t ra t ri t uj4 t r4 t ri t u 3 t
25. ed jdata is the current value of the likelihood term and jprior is the current value of the likelihood term Lastly we recall that as for any properly written IDL program the on line documentation on WISARD or any of its sub routines can be obtained by typing doc library routine name at the IDL prompt where routine name is for instance wisard 3 2 List of parameters and accepted keywords The only positional parameter is data All the other parameters are passed as keywords The notation KEYWORD in the following table means KEYWORD 1 as is customary with IDL ONERA we S THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON MAY 2007 data FOV NP_MIN OVERSAMPLI GUESS NBITER THRESHOLI POSITIV IY NG input input input optional input optional input optional input optional input optional input UNCLASSIFIED rs SANS MENTION DE PROTECTION interferometric input data for the reconstruction data is a vector of structures one structure per time of measurement each containing the following fields VIS2 the squared visibilities for each baseline VIS2ERR the standard deviation on the squared visibilities CLOT the closure phases for each triplet of telescopes involving telescope 1 CLOTERR the standard deviation
26. egularized PSD based left correctly regularized PSD based center and white L1 L2 right The false color table is given on the left side Figure 3 2 shows the plots of WISARD displayed during the white L1 L2 reconstruction These plots aim at showing the quality of the fit of the reconstruction to the data the red crosses show the modulus of ONERA ee S THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 13 SANS MENTION MAY 2007 DE PROTECTION the reconstructed visibilities FT of the object at the measured frequencies while the green squares show the measured visibilities i e the data The blue line shows the difference of the two normalized by 10 times the standard deviation of the visibilities In other words when the blue line is at a level of 0 1 it means that the fit between the reconstructed and the measured visibilities 15 at one standard deviation which means the reconstruction fits the data well 1 07 7 Abs Reconstructed Vis x j Abs Measured Vis 5 Abs Difference 10 x stddev 0 8 0 6 22 0 4 E L T 3 T x i L x d L gt Rpg z 0 0 l I I 3 I I
27. gaussian with the variance o2 Let be the estimator of a from sdata defined by Vsa if sdate gt 0 I 0 else l Although is not gaussian vector we will approximate its distribution by a gaussian one In many cases this approximation is valid as shown in fig II 1 1 0 i Er 7 7 N E x gaussian x gaussian Figure II 1 Comparison between a gaussian distribution of x and a gaussian distribution of x However if c is greater than a fiew a this estimator is biased We have studied the behavior of the mean lt gt and standard deviation o of this estimator in function of 0 a see figs II 2 and II 3 mean a 3 sqrt o 6 N N gt o gt 5 e 5 2 Figure II 2 Mean of the estimator in function of a a We then propose the data model for a a date a q noise ONERA ee m E THE FRENCH AEROSPACE LAB 24 L MUGNIER amp S MEIMON UNCLASSIFIED SANS MENTION MAY 2007 DE PROTECTION T get Figure II 3 Standard deviation of the estimator in function of a a with 0 if sdata lt o o 6 if 0 lt sta lt c Ig Vaita if sqata gt o 6 f vs 2 abate c Oa I if g data gt Os q data EN We also decide to discard the data such that s o lt 0 ONERA gu om NS THE FREN
28. ise smooth prior is obtained by a linear quadratic also called L L gt regularization introduced in 8 9 This prior is called edge preserving as it allows sharp edges in the object if the data is compatible with them contrarily to a quadratic regularization A variant of this regularization has been designed to grant the solution with some smoothness while allowing spikes it is a pixel independent or white linear quadratic regularization called spike preserving in short Both the edge preserving and the spike preserving priors are implemented in j prior 1112 pro The main parameters for WISARD are the data of course the Field Of View FOV for the reconstruction the minimum number of pixel NP MIN the user wants in this FOV and keywords specific to the desired regularization For a quadratic regularization use keyword PSD and possibly keyword MEAN 0 For an edge preserving regularization use keywords DELTA and SCALE For a spike preserving regularization use keywords DELTA SCALE and WHITE 1 See 3 2 for the full list of keywords and details on how to set them Additionally some information is displayed while the iterative minimization is in progress under the following form ITER NBITER Convergence x Criterion jtotal jdata jprior where ITER is the number of iterations performed so far NBITER is the maximum allowed number of iterations x is the value of the convergence test jt otal is the current value of the criterion being minimiz
29. it you can go on and pre process the data file which is in OIFITS format into a structure suitable for input into WISARD datafilename inputdata datalmagingBeautyContest2004 0ifits data wisard oifits2data datafilename Choose the reconstructed Field Of View FOV and grid size The default FOV is the inverse of the minimum spatial frequency present in the data which is often too small if the data lacks low frequencies For a 32 x 32 reconstructed object and a FOV of 18 milli arcseconds for instance type NP min 32L onemas 1d 3 DPi 180D 3600D one mas in radian fov 18 0 onemas Choose the convergence threshold used to stop the iterations By default one waits until the criterion no longer evolves the default threshold being given by machine precision and of the order of 1077 For experimenting with the code you may choose a larger value such as 10 5 threshold 1d 6 You may also give an initial guess if you have performed other reconstructions or if you want to see how stable the reconstruction is with respect to the initial guess By default if the guess is O or absent the so called dirty map is computed by WISARD on the appropriate grid thresholded to positive values and used as an initial guess guess 0 For a first reconstruction you can use a Gaussian or PSD based regularization The PSD must be a 2 D map in Fourier space containing the assumed PSD of the object to be reconstructe
30. n Re 25 II 3 Thescalarcase 0 26 5 BIBLIOGRAPHY uan 3 Re eRe Pe ORG RS ESR ESE RS 28 ONERA EL THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 4 SANS MENTION MAY 2007 DE PROTECTION ABBREVIATIONS AND ACRONYMS CeCILL Cela C nrs I nria L ogiciel L ibre free software license agreement EII European Interferometry Initiative FITS Flexible Image Transport System FP 6 Sixth Framework Programme of the European Union GPL GNU General Public License IAU International Astronomical Union JRA 4 Join Research Activity 4 NA Not Applicable OIFITS Optical Interferometry FITS VLTI Very Large Telescope Interferometer WP Work Package ONERA lt ee THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 5 SANS MENTION MAY 2007 DE PROTECTION 1 INTRODUCTION This report documents the image reconstruction software called WISARD developed by ONERA in WP 2 5 of the JRA 4 within the framework of the FP 6 Together with the software archive it constitutes ONERA s contribution to the outputs of WP 2 5 Image Reconstruction Tools WISARD stands for Weak phase Interferometric Sample Alternating Reconstruction Device It is a software for the reconstruction of images from interferometric data such as the ones that can be recorded by e g the European VLTI It is based on and quite thoroughly described in the
31. nal Process 37 12 pp 2024 2036 December 1989 3 1 J Idier editor Approche bay sienne pour les probl mes inverses Herm s Paris 2001 3 1 J M Conan L M Mugnier T Fusco V Michau and G Rousset Myopic Deconvolution of Adaptive Optics Images by use of Object and Point Spread Function Power Spectra Appl Opt 37 21 pp 4614 4622 July 1998 3 1 L M Mugnier C Robert J M Conan V Michau and S Salem Myopic deconvolution from wavefront sensing J Opt Soc Am A 18 pp 862 872 April 2001 3 1 ONERA we ee THE FRENCH AEROSPACE LAB L MUGNIER amp S MEIMON UNCLASSIFIED 29 SANS MENTION MAY 2007 DE PROTECTION 9 10 11 12 L M Mugnier T Fusco and J M Conan MISTRAL a Myopic Edge Preserving Image Restoration Method with Application to Astronomical Adaptive Optics Corrected Long Exposure Images J Opt Soc Am A 21 10 pp 1841 1854 October 2004 3 1 3 2 4 1 S C Meimon L M Mugnier and G Le Besnerais A novel method of reconstruction for weak phase optical interferometry In New frontiers in stellar interferometry edited by W A Traub vol 5491 pp 909 919 Proc Soc Photo Opt Instrum Eng 2004 Conference date June 2004 Glasgow UK 3 3 T J Cornwell and P N Wilkinson A new method for making maps with unstable radio interferomete
32. on the closure phases FREQS U first coordinate of the spatial frequencies involving telescope 1 FREQS V second coordinate of the spatial frequencies involving telescope 1 The zero frequency must not be present in the data the data must be normalized as are OIFITS files so that the value of the data at this frequency would be 1 In other words the reconstructed object has a unit sum Field Of View of the reconstructed image in units consistent with the data More precisely the unit for FOV must be the inverse of the unit of the arrays of frequencies FREQS_U and FREOS V of the data Usually FREQS_U and FREQS_ V are in rd 1 so FOV should be in rd radians MINimum width Number of Points of the reconstructed image The Number of Points of the reconstructed object may be greater depend ing on FOV OVERSAMPLING factor and frequencies present in the data See routine WISARD MAKE 11 for details oversampling factor for the reconstructed image By default OVERSAMPLING 1 and the maximum spatial frequency of the re constructed object is the maximum frequency of the data See routine WISARD MAKE H for details initial guess for the reconstructed image or dirty map if not present This guess is massaged in the following way resampled if necessary to the correct size thresholded to positive values and normalized to a unit sum This massaged guess is available on exit as a field of AUX OUTPUT ma
33. rs Mon Not R Astr Soc 196 pp 1067 1086 1981 4 1 J W Goodman Statistical optics John Wiley amp Sons New York 1985 4 3 1 ONERA weer ee THE FRENCH AEROSPACE LAB
34. ruction method is a regularized reconstruction 5 6 where the solution recon structed object is defined as the one that minimizes a compound criterion or metric with two terms One term called the likelihood term measures the fit of the reconstruction to the data and the other term regu larization or penalty term measures the compatibility of the reconstruction with our prior knowledge on the true object The main routine of WISARD is wisard pro The input data for wisard pro is typically recorded by an optical interferometer measuring in time the simultaneous interferences of at least three tele scopes It consists in squared visibilities and closure phases the error bars on these quantities and the spatial frequencies corresponding to each datum The structure of the data is given in Section 3 2 and more precisions are available in chapter 4 The prior knowledge is on the one hand the hard constraint of positivity of the reconstruction and on the other hand a somewhat fuzzy knowledge that the object to be reconstructed has some kind of smoothness or piece wise smoothness or global smoothness apart from some spikes etc In WISARD several regular ization terms are available to embody this prior knowledge on the solution A smooth solution is obtained by a quadratic regularization which uses the object s PSD Power Spectral Density as input This regularization is described in 7 and implemented in j prior gauss pro A piecew
35. ximum NumBer of ITERations for the reconstruction 500 by de fault For a better control of the reconstruction one should rather use THRESHOLD below and not lower this value convergence THRESHOLD to be used as a stopping criterion for the iterations By default set to the machine precision in simple precision approximately 1 19e 07 although computations are done in double precision For a rather quick look result set to a smaller value e g 1076 POSITIVITY constraint for the reconstruction It 15 set to true 1 by default Set it explicitly to 0 if by misplaced curiosity you do not want to use the positivity constraint in the reconstruction ONERA ee s THE FRENCH AEROSPACE LAB L MUGNIER amp S MAY 2007 PSD MEAN O DELTA SCALE WHI TE LIBRARY AUX OUTPUT ISPLAY VERSI HELP COPYRI ON GHT MEIMON optional input optional input optional input optional input optional input optional input optional output optional input optional input optional input optional input 3 3 Example of use of WISARD UNCLASSIFIED 10 SANS MENTION DE PROTECTION 2 D map of size NP MINXNP MIN containing the PSD for the quadratic regularization of the reconstruction See routine J PRIOR GAUSS and the example file for details 2 D map of size NP MINXNP MIN containing the MEAN Object to be used for
Download Pdf Manuals
Related Search