Physics-based model of the adaptive-optics corrected point-spread-function
Romain Fétick, Thierry Fusco, Benoit Neichel, Laurent Mugnier, Olivier Beltramo-Martin, Aurélie Bonnefois, Cyril Petit, Julien Milli, Joel Vernet, Sylvain Oberti, Roland Bacon
AAstronomy & Astrophysics manuscript no. aanda c (cid:13)
ESO 2019August 8, 2019
Physics-based model of the adaptive-optics-correctedpoint spread function (cid:63)
Applications to the SPHERE/ZIMPOL and MUSE instruments
R. JL. Fétick , , T. Fusco , , B. Neichel , L. M. Mugnier , O. Beltramo-Martin ,A. Bonnefois , C. Petit , J. Milli , J. Vernet , S. Oberti , and R. Bacon Aix Marseille Univ, CNRS, CNES, LAM, Marseille, Francee-mail: [email protected] ONERA, The French Aerospace Lab BP72, 29 avenue de la Division Leclerc, 92322 Chatillon Cedex, France European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura, Casilla 19001,Santiago, Chile ESO, European Southern Observatory, Karl-Schwarzschild Str. 2, 85748 Garching bei Muenchen, German CRAL, Observatoire de Lyon, CNRS, Université Lyon 1, 9 Avenue Ch. André, F-69561 Saint Genis Laval Cedex, FranceReceived May 3 rd th ABSTRACT
Context.
Adaptive optics (AO) systems greatly increase the resolution of large telescopes, but produce complex point spread function(PSF) shapes, varying in time and across the field of view. The PSF must be accurately known since it provides crucial informationabout optical systems for design, characterization, diagnostics, and image post-processing.
Aims.
We develop here a model of the AO long-exposure PSF, adapted to various seeing conditions and any AO system. This modelis made to match accurately both the core of the PSF and its turbulent halo.
Methods.
The PSF model we develop is based on a parsimonious parameterization of the phase power spectral density, with only fiveparameters to describe circularly symmetric PSFs and seven parameters for asymmetrical ones. Moreover, one of the parameters isthe Fried parameter r of the turbulence’s strength. This physical parameter is an asset in the PSF model since it can be correlatedwith external measurements of the r , such as phase slopes from the AO real time computer (RTC) or site seeing monitoring. Results.
We fit our model against end-to-end simulated PSFs using the OOMAO tool, and against on-sky PSFs from theSPHERE / ZIMPOL imager and the MUSE integral field spectrometer working in AO narrow-field mode. Our model matches theshape of the AO PSF both in the core and the halo, with a relative error smaller than 1% for simulated and experimental data. We alsoshow that we retrieve the r parameter with sub-centimeter precision on simulated data. For ZIMPOL data, we show a correlationof 97% between our r estimation and the RTC estimation. Finally, MUSE allows us to test the spectral dependency of the fitted r parameter. It follows the theoretical λ / evolution with a standard deviation of 0 . Key words.
Instrumentation: adaptive optics – Methods: analytical, observational – Atmospheric e ff ects – Telescopes
1. Introduction
Optical systems su ff er from aberrations and di ff raction e ff ectsthat limit their imaging performance. For ground-based observa-tions, the point spread function (PSF) is dramatically altered bythe atmospheric turbulence that distorts the incoming wavefront(Roddier 1981). The resolution under typical conditions of aseeing-limited telescope does not exceed the di ff raction limitof an ∼
12 cm aperture. Modern and future large telescopesthus include adaptive optics (AO) systems (Roddier 1999) thatcompensate for the atmospheric turbulence thanks to wavefrontsensors and deformable mirrors. The aberrated wavefrontis partially corrected and telescopes may operate near theirdi ff raction limited regime. Nevertheless the AO correction islimited by technical issues such as sensor noise, limited numberof actuators, or loop delay (Martin et al. 2017; Rigaut et al.1998). This results in a peculiar shape of the PSF made of asharp peak due to the partial AO correction, and a wide halo (cid:63) Our Python codes will be available on https://gitlab.lam.fr/lam-grd-public as soon as the paper is published. caused by the residual turbulence above the AO cuto ff frequency.The PSF thus provides critical information about an opti-cal system regarding its preliminary design, calibrations, test-ings, or diagnostic (Ascenso et al. 2015; Ragland et al. 2018).Image post-processing, such as deconvolution (Mugnier et al.2004), also requires knowledge of the PSF. Deconvolution oflong-exposure images using parametric PSFs has already beendemonstrated in Drummond (1998) and Fétick et al. (2019). Afine model of the PSF is necessary. The substantial advantageof parametric PSFs is to compress all the important informationof the physical PSF into a small number of parameters. The nu-merical values of these parameters might then be used for com-parisons, correlations, or any statistical analysis. Moreover if thePSF parameters are correlated to physical values (e.g. turbulencestrength, wind speed, AO residual phase variance), it is possibleto better constrain these parameters or better understand the AOresponse to given observing conditions. We state that an e ffi cientAO PSF model should fulfil the following requirements: Article number, page 1 of 11 a r X i v : . [ a s t r o - ph . I M ] A ug & A proofs: manuscript no. aanda – Accuracy . The model must represent accurately the shape ofthe AO-corrected PSF, especially the two areas correspond-ing to its central peak and to its wide turbulent halo. The re-quested accuracy depends on the application of the PSF (e.g.fitting, deconvolution, turbulence monitoring). – Versatility and robustness . The model must be used on dif-ferent AO system, with di ff erent AO correction levels, fordi ff erent turbulent strengths. – Simplicity . The model must have as few parameters as pos-sible without damaging its versatility or accuracy. – Physical parameters . Such parameters have a physicalmeaning related to the observing conditions. These param-eters have physical units.The literature already provides some models of AO-corrected PSF (Drummond 1998; Zieleniewski & Thatte 2013)often based on Gaussian, Lorentzian, and / or Mo ff at (1969)models. A trade-o ff is always drawn between a simple modelwith few parameters but imprecise, or a more precise but alsomore complex model. The di ffi culty often comes from thedescription of the turbulent halo with only a few parameters.Moreover, to the best of our knowledge, these PSF modelsrely only on mathematical parameters without direct physicalmeaning or units.We propose a long-exposure PSF model for AO-correctedtelescopes that describes accurately the shape of the PSF; thismodel is made of a small number of parameters with physicalmeaning whenever possible. Our method does not parameterizethe PSF directly in the focal plane, but rather from the phasepower spectral density (PSD). Indeed Goodman (1968) andRoddier (1981) have shown that the phase PSD containsall the necessary information to describe the long-exposureatmospheric PSF. Working in the PSD domain allows us toinclude physical parameters. Then Fourier transforms give theresulting PSF in the focal plane. Our PSF model also includespupil di ff raction e ff ects or any of the system static aberrations,provided they have been previously characterized.In Sect. 2 we first recall the expression of the Mo ff at functionand show its limits for AO PSF description. Then we developour PSF model, partially based on this Mo ff at function. Section3 validates the model by fitting PSFs from numerical simulationsand from observations made on two Very Large Telescope (VLT)instruments. Finally Sect. 4 concludes our work and discussesdirect and future applications for our PSF model.
2. Description of the PSF model
In the whole paper, we define ( x R , y R ) the reference coordinatesthat are respectively the detector horizontal and vertical coor-dinates. We also define ( x , y ) the proper PSF coordinates (e.g.along the major and minor PSF elongation axis), rotated by anangle θ R with respect to the reference frame. The reference frameto PSF frame transformation can be written as (cid:32) xy (cid:33) = (cid:32) cos θ R sin θ R − sin θ R cos θ R (cid:33) (cid:32) x R y R (cid:33) . (1)For the sake of simplicity we will use mainly the rotated coordi-nates, but it is important to keep in mind that θ R is a crucial PSFparameter. We also simplify the notations x = x ( x R , y R , θ R ) and y = y ( x R , y R , θ R ). In this section we first recall the usual Mo ff at PSF model,since it encompasses and generalizes Lorentzian and Gaussianmodels. We demonstrate the advantages of the Mo ff at model,but also its limitations. This motivates our search for a betterPSF model. However, the mathematical expression of the Mo ff atfunction will still be used inside our more complete PSF model. The AO-corrected PSF exhibits a sharp corrected peak, withwide wings extension. The Mo ff at (1969) model is often useddue to its good approximation of the AO PSF sharp peak (An-dersen et al. 2006; Müller Sánchez et al. 2006; Davies & Kasper2012; Orban de Xivry et al. 2015; Rusu et al. 2016). The Mo ff atfunction, of amplitude A , is written as M A ( x , y ) = A (1 + x /α x + y /α y ) β , (2)with α x , α y , and β strictly positive real numbers. Moreover thecondition β > β = β → + ∞ . The variable β parameter thus makesthe Mo ff at function a generalization of Lorentzian and Gaussianones. Since the PSF has a unit energy, demonstration in Ap-pendix A shows that the Mo ff at multiplicative constant is A = β − πα x α y , (3)so the PSF, called h , is made of only four free parameters α x , α y , β, and θ R . The Mo ff at PSF model thus is re-written as h ( x , y ) = β − πα x α y M ( x , y ) , (4)where the notation M must be understood as the Mo ff at M A with a multiplicative factor A = ff at function to mo-tivate our search for better functions. Indeed, as shown in Fig. 1,the Mo ff at function accurately fits the central peak of an actualPSF, but poorly describes the turbulent halo. Since this halo maycontain an important proportion of the PSF energy, depending onthe quality of the AO correction, it is necessary to model it ac-curately. Adding a constant background to the model artificiallyimproves the fitting (lower residuals). However, this method isnot suitable since it poorly describes the halo and mistaking thehalo for a background will yield the non-physical result of a PSFwith an infinite integral on an unlimited field of view. The mod-ulation transfer functions (MTF, bottom plot in Fig. 1), which isthe modulus of the PSF Fourier transform, also shows that theMo ff at does not match well the very low frequencies (halo) anddoes not model the telescope cuto ff frequency. Similarly, noneof the static aberrations of the telescope are taken into account.A more physical PSF model than a Mo ff at is thus required. Our PSF model is based on equations of image formation fromthe phase PSD to the focal plane. Indeed Roddier (1981) hasshown that the Fourier transform of the PSF, the optical transfer
Article number, page 2 of 11. JL. Fétick et al.: Physics based model of the AO-corrected PSF
Fig. 1.
Fitting of a SPHERE / ZIMPOL PSF (blue) using a Mo ff at modelwith background (green) and without background (dashed green). Top:PSF, Bottom: MTF. The insert plot is a zoom on the low spatial frequen-cies. function (OTF), can be written as the product of the telescopeaberrations OTF and the atmospheric turbulent OTF,˜h( ρ/λ ) = ˜h T ( ρ/λ ) · ˜h A ( ρ/λ ) , (5)where λ is the observation wavelength, ˜ h the total OTF, ˜h T thetelescope OTF, and ˜h A the atmospheric OTF. This OTF splittingequation is valid under the hypothesis of a spatially stationaryphase. This is the case for a purely turbulent phase, and a goodapproximation for an AO-corrected phase (Conan 1994). To es-tablish this result, Roddier also used the fact that the phase distri-bution follows a Gaussian process, as the sum of a large numberof independent turbulent layers. The telescope OTF is simplygiven by the autocorrelation of the pupil transmission function,whereas the atmospheric OTF is written as˜h A ( ρ/λ ) = e − B φ ( ) e B φ ( ρ ) , (6)with B φ the phase autocorrelation function defined as B φ ( ρ ) = (cid:104)(cid:104) φ ( r , t ) φ ( r + ρ, t ) (cid:105) t (cid:105) r . (7)The Wiener-Khintchine theorem states that the PSD is theFourier transform of the autocorrelation as W φ ( f ) = F { B φ ( ρ ) } , (8)where W φ denotes the phase PSD and F the Fourier operator; f and ρ are the Fourier conjugated variables. If we call u theangular variable conjugated to ρ/λ , the PSF is written as h ( u ) = F − (cid:110) ˜h T ( ρ/λ ) e − B φ ( ) e F − { W φ ( f ) } (cid:111) . (9) We note that B φ ( ) is the residual phase variance and is equal tothe integral of W φ on the whole frequency plane. This equationshows that only knowledge of the pupil and the static aberrationsof the term ˜h T , and the phase PSD W φ are necessary for thedescription of the long-exposure PSF. Di ff raction e ff ects – suchas finite aperture, central obstruction, and spiders – only dependon the pupil geometry and are known. Static aberrations aresecond-order e ff ects that can be either neglected (as we show inSects. 3.3 and 3.4), or measured (N’Diaye et al. 2013) and thenincluded in the ˜h T term for more accuracy. The term ˜h T beingfully determined, now we only have to parameterize the residualphase PSD to model the PSF. Actuators controlling the deformable mirror are separated by apitch that sets the maximal spatial frequency of the phase thatcan be corrected by the AO system. This is called the AO spa-tial cuto ff frequency, defined by f AO (cid:39) N act / D , where N act isthe linear number of actuators and D the pupil diameter. Thistechnical limitation induces the peculiar shape of the AO resid-ual phase PSD (and in fine a peculiar shape of the PSF). Theresidual phase PSD is thus separated into two distinct areas: – AO-corrected frequencies f ≤ f AO, – AO-uncorrected frequencies f > f AO.
The uncorrected area is not a ff ected by the AO system, and thephase PSD consequently follows the Kolmogorov law, W φ, Kolmo ( f ) = . r − / f − / , for f > f AO , (10)where r is the Fried parameter scaling the strength of theturbulence. The halo is thus set by the knowledge of only this r parameter.Regarding the AO-corrected area, it is di ffi cult to parameter-ize the phase residual PSD since it depends on the turbulence,the magnitude of the object, the AO loop delay, and the wave-front reconstruction algorithm. Our objective is not to build afull reconstruction of the phase PSD, but to only get a modelthat can match it. Racine et al. (1999) and Jolissaint & Veran(2002) have shown that in extreme AO correction (small resid-ual phase) the shape of the PSF is exactly the shape of the PSD.For partial AO correction, the shapes of the PSF and PSD are notexactly identical, but are still similar (Fétick et al. 2018). More-over Rigaut et al. (1998) have shown that the AO residual PSDis the sum of decreasing power laws of the spatial frequency. AMo ff at function used in the PSD domain would already describetwo regimes due to its shape, one regime for f ≤ α and oneregime for α < f < f AO . Adding a constant under the Mo ff at al-lows us to describe a third regime near the AO cuto ff frequencyat f (cid:46) f AO that is roughly similar to the shape of the aliasingPSD discussed by Rigaut et al. (1998). All the above pieces ofinformation suggest the possibility of using the Mo ff at functionfor a parsimonious parameterization of the AO-corrected PSD,rather than using it to directly parameterize the PSF in the focalplane. The full PSD model is written as W φ ( f ) = β − πα x α y M A ( f x , f y )1 − (cid:18) + f α x α y (cid:19) − β + C f ≤ f AO + (cid:104) W φ, Kolmo ( f ) (cid:105) f > f AO , (11) Article number, page 3 of 11 & A proofs: manuscript no. aanda where the Mo ff at normalization factor ensures a unit integral ofthe Mo ff at on the area f ≤ f AO (see Appendix A). Constant C is an AO-corrected phase PSD background. It is useful to modelthe residual AO PSD near the AO cuto ff , where the Mo ff at func-tion is close to zero. Thus the AO residual phase variance on thecircular domain below the AO cuto ff frequency is directly σ = A + C π f . (12)Parameter A , added to the C background contribution, has thephysical meaning of being the residual variance. An example ofour PSD model is given Fig. 2. We do not impose continuity atthe AO cuto ff frequency, so the PSD might be locally discon-tinuous. Indeed the transition area f (cid:39) f AO between correctedand uncorrected frequencies can lead to strong PSD gradients,which are modelled by an eventual PSD discontinuity. Fig. 2.
Three components of the PSD model: the Mo ff at (blue) and theconstant contribution (orange) below the AO cuto ff frequency, and theKolmogorov spectrum (green) above the AO cuto ff frequency. Discon-tinuity has been exaggerated by reducing C to show this degree of free-dom in our model. Plotting is in logarithmic-logarithmic scale. Our PSF model based on the PSD is made of the followingset of seven parameters: S = { α x , α y , β, θ R , C , r , A } in the asym-metric case. This reduces to five parameters in the symmetriccase (setting α y = α x and θ R = ffi cient in the majority of cases, asymmetries makeit possible to consider PSFs elongated due to strong wind e ff ectsor anisoplanetism. Once the parameters S are set, the PSF isthen computed from the AO PSD and static aberrations usingEq. (9).For the reader interested in deriving the Strehl ratio from ourmodel, we have to compute the integral of the Kolmogorov spec-trum above the AO cuto ff frequency, σ = . r − / π (cid:90) + ∞ f AO f − / f df = .
023 6 π r · f AO ) − / . (13)The Strehl ratio consequently is written asS R = exp (cid:104) − ( σ + σ ) (cid:105) = exp (cid:34) − A − C π f − .
023 6 π r · f AO ) − / (cid:35) . (14)
3. Validation
In this section we deal with images of PSFs (the data) that maycome from numerical simulations or observations of stars onVLT instruments. The fitting method consists in finding the PSFparameters so that the model PSF minimizes the square distanceto the data PSF: L ( S , γ, ζ, δ x , δ y ) = (cid:88) i , j w i , j (cid:104) γ · h i , j ( S , δ x , δ y ) + ζ − d i , j (cid:105) , (15)where h i , j is the discretized model of PSF on the pixels ( i , j ), S its set of parameters, and d i , j is the data PSF. Since the PSFmodel (given by Eqs. 9 and 11) has a unit flux, it is scaled by γ to match the flux of the data PSF, and ζ accounts for a possiblebackground. The shifts δ x and δ y centre the PSF with sub-pixelprecision on the data (by multiplication of the OTF with the cor-rect phasor). The weighting factor w i , j is the inverse of the noisevariance, which takes into account the photon noise and the de-tector read-out noise. As noted by Mugnier et al. (2004), for highfluxes (typically greater than ten photons per pixel), the Poissonphoton noise becomes nearly Gaussian and the weighting factoris written as w i , j = { d i , j , } + σ . (16)In this case, our approach can be seen from a statistical point ofview as maximizing the likelihood of the data d i , j corrupted byphoton and read-out noise. We thus minimize L , which is theneg-logarithm of the likelihood for a Gaussian process.Let us now note that the minimum of L has an analytic so-lution for γ and ζ (see Appendix B for a full demonstration).We actually do not need to numerically minimize over these twoparameters, the least-square criterion only relies on the PSF in-trinsic parameters S and the position parameters ( δ x , δ y ) as L (cid:48) ( S , δ x , δ y ) = L ( S , ˆ γ, ˆ ζ, δ x , δ y ) , (17)where ˆ γ and ˆ ζ are the analytic solutions for the flux and the back-ground, respectively. At each iteration of the minimization pro-cess, the minimizer evaluates our L (cid:48) criterion with a new setof parameters ( S , δ x , δ y ). The current PSF estimate h ( S , δ x , δ y )is computed, then the analytic solutions ˆ γ and ˆ ζ are computed.The quantity ˆ γ · h ( S , δ x , δ y ) + ˆ ζ is used to compute the residualswith the data d . Residuals are then provided to the minimizer toestimate a new set of parameters ( S , δ x , δ y ). The Object-Oriented Matlab Adaptive Optics (OOMAO) tool-box, presented by Conan & Correia (2014), provides end-to-endsimulations. For each time step, OOMAO generates a turbulentwavefront with a Von-Kármán spectrum defined as W φ, VK ( f ) = . r − / (cid:32) L (cid:33) + f − / , (18)where we have chosen the outer scale L =
30 m. Since1 / L (cid:28) f AO , the Von-Kármán spectrum is consistent with ourPSF model using the Kolmogorov spectrum above the AO cuto ff frequency. OOMAO then propagates the wavefront through Article number, page 4 of 11. JL. Fétick et al.: Physics based model of the AO-corrected PSF the telescope, simulates the wavefront sensor measurement,performs the wavefront reconstruction, and simulates the mirrorwavefront deformation. For each time step, we get a shortexposure PSF. Integration over time allows us to retrieve thelong-exposure PSF. It is important to notice that the methodto compute the PSF is then very di ff erent from our model,which directly uses the residual phase PSD in Eq. (9). We usedOOMAO to generate a set of PSFs, corresponding to di ff erentwavelengths from 0 . .
18 um, and seven di ff erent valuesof r from 7 . . r are given at 500 nm. Wetranslate them from the observation wavelength to the referencewavelength of 500 nm using the theoretical spectral depen-dency r ,λ / r ,λ = ( λ /λ ) / . For all the PSFs, the telescopeparameters are kept unchanged D = N act =
32, samplingat Shannon-Nyquist for all wavelengths. The phase screenconsists in one frozen flow (Taylor’s hypothesis) turbulent layertranslating at v =
10 m / s. Using these parameters, we generatedPSFs corresponding to exposure times of 0 .
1, 1, 10 , and 100seconds. For 0 . r estimation. For a 10 second exposure PSF,the random fluctuations of the phase are correctly averaged.This is confirmed by the 100 second exposure PSF, whichgives the same r estimation as the 10 second case. Since a 100second exposure is computationally demanding and does notsignificantly improve the results, we performed all our tests onthe 10 second exposure time.Parameter Values UnitDiameter D N act
32 –Wavelength λ r / sOuter scale L
30 mExposure 0 . , , ,
100 s
Table 1.
OOMAO parameters summary for our PSF simulations.
Param. Typical range Lower bound Guess Unit r −
30 eps 18 cm α x − − − eps 5 × − m − α y − − − eps 5 × − m − β . − + eps 1 . − θ − π − C − − − − rad m A − −
10 0 2 rad Table 2.
Typical range of PSF parameters for OOMAO simulations andSPHERE / ZIMPOL instrument, lower bounds and values used as ini-tial guess for the minimizer. Typical ranges are indicative and may varyaccording to the considered instrument. The value ’eps’ denotes the ma-chine precision. Parameters do not have any upper bound.
Each PSF is then fitted using the L (cid:48) ( S , δ x , δ y ) criteriongiven to an optimizer (e.g. Levenberg-Marquardt, trust region, orMarkov chain Monte Carlo). We used the Trust Region Reflec-tive algorithm, called ‘trf’, from the Python / SciPy (Jones et al.2001) library. This algorithm is gradient based, said to be robust, and allows bounds on the parameters. The robustness of this al-gorithm was verified for our applications of it to PSF fitting, eventhough the convexity of the problem is not demonstrated. So far,we have not encountered any local minimum and residuals arealways small. For all PSFs, the same initial conditions are pro-vided to the fitting algorithm (see Table 2), in particular we usedthe same value of r =
18 cm. Using the same initial parame-ters {S , δ x , δ y } init for all fits ensures that our model is suited forminimization procedures and that convergence is ensured even ifstarting far from the true values. Fitting results are presented onFig. 3. Our model fits well the OOMAO-generated PSF on boththe corrected and the uncorrected area, residuals being on aver-age one to two decades below the PSF. Let us define the relativeerror between fitted PSF and data PSF as (cid:15) h = (cid:113)(cid:80) i , j (cid:104) γ · h i , j ( S , δ x , δ y ) + ζ − d i , j (cid:105) (cid:80) i , j d i , j . (19)This error is the L2 norm of the di ff erences between fitting anddata, relative to the flux. Considering all the OOMAO fitting,we find an average relative error (cid:15) h = . × − with a standarddeviation of 2 . × − . For comparison, fitting with a Mo ff atmodel gives an average relative error of (cid:15) h = . × − with astandard deviation of 1 . × − . Using our model thus increasesthe fitting accuracy by a factor of approximately 5 with respectto the former Mo ff at model of Sect. 2.1. Regarding the OTFs,the fit is also accurate on the whole frequency range, from lowfrequencies (mainly the halo) to high frequencies (PSF peak andtelescope cuto ff ). Our model slightly over-estimates frequenciesjust below the telescope cuto ff frequency but this has never beenan issue in our applications, such as deconvolution, and is stillmuch better than the Mo ff at OTF (which has no telescope cuto ff frequency and give a poor estimation of the low frequencies).Regarding the flux, we consider the relative error betweenthe flux ˆ γ analytically estimated by our model fitting method,and the OOMAO flux that is directly the sum of the data on allthe pixels: (cid:15) γ = ˆ γ − (cid:80) i , j d i , j (cid:80) i , j d i , j . (20)On all our OOMAO simulations, we find an average relative er-ror of − . . − . , . r estimation As shown in Fig. 4, our r estimation is consistent with theOOMAO value of r . We find the best linear fit to be r , FIT = . r , OOMAO − . , (21)where values are given in centimetres. ThePearson correlation coe ffi cient is C Pearson = Cov( r , FIT , r , OOMAO ) / (cid:112) Var( r , FIT ) · Var( r , OOMAO ) = . r estimation with respect toOOMAO simulations with sub-centimetre precision. σ AO estimation Theoretically our model should also be able to retrieve the resid-ual variance σ AO on the corrected area and follow a λ − power Article number, page 5 of 11 & A proofs: manuscript no. aanda
Fig. 3.
OOMAO PSF fitting with our model. Left: Circular average forPSFs (given in photons). The vertical grey line corresponds to the AOcuto ff radius. Right: Corresponding circular average for OTFs (normal-ized to unity at the null frequency). From top to bottom, three wave-lengths are scanned from 500 nm to 1220 nm. Colours correspondto three values of the OOMAO required r . Solid curves: OOMAO.Dashed: fitting. Dotted: residuals. All PSFs, for all wavelengths, aresampled at Shannon-Nyquist. Fig. 4.
Fried parameter r estimated by fitting versus the r used inOOMAO to generate the PSF. All r are given at 500 nm. Here areshown results on 84 di ff erent PSFs, corresponding to seven values of r and 12 di ff erent wavelengths. The line is the linear fit between our r estimation and OOMAO r . Crosses show residuals | r , FIT − r , OOMAO | . Alog-log scale is used to show on the same graph both data and residuals. law. Figure 5 shows the fitting estimation of this variance versusthe wavelength. This data is then fitted with curves of equation a λ − . Except the two outliers for minimal r = . λ (cid:39) λ − power law is a good estima-tion of the σ evolution. This result gives confidence in theestimated parameter. Data from the real time computer (RTC)could be used in the future to provide the σ parameter for PSFestimation. The λ − spectral dependence is also an asset to shiftthe PSF from one wavelength to another. Fig. 5.
Estimation of the σ AO from PSF fitting versus the wavelength(dots). Colours correspond to the seven di ff erent values of OOMAO r .Curves of parametric equations σ = a λ − are fitted on the data. C estimation The C term in Eq. (11) accounts for multiple sources of residualPSD, including wavefront aliasing and other AO residual errors.Since this constant is dominated by the Mo ff at PSD in the core,it becomes more important near the AO cuto ff , where the alias-ing dominates. Since the aliasing scales in r − / (Rigaut et al.1998), we look for similar r dependencies for the PSF constant C . Figure 6 shows a clear decrease of C with r . Fitting the esti-mated C with a r − / power law shows a good match, with smallresiduals for nearly all r values. The power law is not exactly − / (cid:39) − .
67 but is closer to − .
46 for this OOMAO case.However, one can still think about normalizing the constant inEq. (11) by C = C (cid:48) r − / (22)and perform fitting over C (cid:48) instead of C . This would reduce thevariation range of this parameter. Reducing bounds or standarddeviation of a parameter is an asset for constraining the modeland improving minimization processes. The Spectro Polarimetric High-contrast Exoplanet REsearch(SPHERE) instrument (Beuzit et al. 2008; Beuzit et al. 2019) ofthe VLT includes the powerful SPHERE Adaptive optics for eX-oplanets Observation (SAXO) system described in Fusco et al.(2014) and Sauvage et al. (2010). The AO real time computer isbuilt on the ESO system called SPARTA (Fedrigo et al. 2006),which stands for Standard Platform for Adaptive optics RealTime Applications. In particular, for each observation, SPARTA
Article number, page 6 of 11. JL. Fétick et al.: Physics based model of the AO-corrected PSF
Fig. 6.
Estimation of the PSF constant C versus the r given at the ob-served wavelength (dots). A r − / fitting equation (solid line) is appliedon the data. Residuals between each data point and the r − / power laware also shown (crosses). is able to give an estimate of r from the mirror voltages andwavefront sensor slopes.The Zurich IMaging POLarimeter (ZIMPOL) instrument(Schmid et al. 2018) is mounted at the focal plane of SPHERE.ZIMPOL is also used as a very e ffi cient imager at visiblewavelengths. One of its applications in non-coronagraphic modeis the observation of asteroids (Vernazza et al. 2018; Viikinkoskiet al. 2018; Fétick et al. 2019) as part of an ESO Large Program(ID 199.C-0074, PI P. Vernazza). PSFs from stars were observedwith the ZIMPOL N_R filter (central wavelength 645.9 nm,width 56.7 nm) during the Large Program. When PSFs aresaved together with the SPARTA telemetry, we are able tocorrelate r given by our fitting and r given by SPARTA.In our sample, 28 PSFs were saved along with the SPARTAtelemetry. These PSFs were obtained during di ff erent nights, onstars of di ff erent magnitudes, with various seeing conditions.We fitted these 28 PSFs with our PSF model. Figure 7 showsthree of the 28 fittings, for the smallest r of the sample, themedian r , and the largest r , respectively. Our fitted PSFsmatch the shape of the core and halo. The average of relativeerror defined in Eq. (19) is (cid:15) h = . × − , with a standarddeviation of 1 . × − . We only consider di ff raction due to thetelescope 8 m aperture and its central obstruction in the staticOTF ˜h T (see Eq. (9)). Consequently non-circularly-symmetrice ff ects, such as the spiders or static aberrations visible onFig. 7, are not modelled. Even if the spiders could have beenincluded in our model, we have deliberately chosen to ignorethem in the ˜ h T term since they are negligible in comparison tothe other dominant e ff ects (AO residual core, turbulent halo,8 m aperture di ff raction). When performing PSF fitting, thecontribution of these e ff ects not taken into account in ˜ h T mightbias the atmospheric term ˜ h A during fitting procedure andslightly o ff set the estimation of the S parameters. Moreoverthese ZIMPOL images are field stabilized, meaning rotatingspiders, which are harder to model. Pupil stabilized imageswould make the description of the spider di ff raction e ff ect easier.The top graph on Fig. 8 shows that the values estimated bySPARTA (median =
22 cm) are greater than the values given byfitting (median =
13 cm). The best linear fit between r estimated Fig. 7.
Three ZIMPOL PSFs (top), model fittings (middle), residuals(bottom). Left: Minimal r of the sample. Middle: Median r . Right:Maximal r . The main di ff erences are due to some static aberrations nottaken into account in our model (only the pupil and its central obstruc-tion are taken into account). The hyperbolic arcsine of the intensity isshown to enhance details. The same intensity scale is used per column(data, model, residuals), but di ff ers between columns. by fitting and SPARTA is r , SPARTA = . r , FIT − . . (23)The Pearson correlation coe ffi cient between thetwo series r , FIT and r , SPARTA is C Pearson = Cov( r , FIT , r , SPARTA ) / (cid:112) Var( r , FIT ) · Var( r , SPARTA ) = . r , FIT and r , SPARTA are not identical, however they show a strong correlation. Wefurther investigated the di ff erence between the SPARTA and thefitting estimates thanks to the ESO atmospherical monitoringusing the Multi-Aperture Scintillation Sensor (MASS) com-bined with the Di ff erential Image Motion Monitor (DIMM).Since the MASS / DIMM instrument is located apart from thetelescopes, it does not see exactly the same turbulent volumeas the telescopes and does not su ff er the same dome e ff ect.There might be some uncertainties between MASS / DIMM r estimations and telescope r estimations (PSF fitting or RTC)due to the spatial evolution of the turbulence. Nevertheless thisinstrument is a valuable indicator of the Paranal atmosphericstatistics. For each PSF observation we retrieved the associatedMASS / DIMM seeing estimation within a delay of ± / DIMM is 0 . .
46” for SPARTA and 0 .
83” for PSF fitting.The over-estimation of the SPARTA r with respect to theMASS / DIMM r has been already discussed by Milli et al.(2017). The exact origin of the di ff erence between these threeestimations has not be found. Nevertheless we note that estima-tions with SPARTA are based on RTC measurements of the lowspatial frequencies of the phase (sensitive to the Von-Kármánouter scale L ), whereas our fitting method is based on the PSFhalo corresponding to the high spatial frequencies. Our PSF Article number, page 7 of 11 & A proofs: manuscript no. aanda fitting method might be sensitive to telescope internal wavefronterrors if they have not been previously calibrated and taken intoaccount in the PSF model.Even if there is still an uncertainty on the true value of r , thestrong correlation between SPARTA and fitting estimations issu ffi cient for many applications. Indeed it is still possible to get r , SPARTA from telemetry, use Eq. (23) to translate it into r , FIT ,and get an estimate of the PSF halo. This method constrains themodel for future PSF estimations without having access to theactual image of the PSF.
The Multi-Unit Spectroscopic Explorer MUSE (Bacon et al.2006, 2010) is an integral field spectrograph (IFS) workingmainly in the visible, from ∼
465 nm to ∼
930 nm. MUSE isequipped with the Ground Atmospheric Layer Adaptive Opticsfor Spectroscopic Imaging (GALACSI) adaptive optics system(Ströbele et al. 2012) to improve its spatial resolution in twodi ff erent modes, the so-called narrow-field mode (NFM) andwide-field mode (WFM), to correct di ff erent sizes of field ofview. The AO facility uses four laser guide stars (LGS) (Caliaet al. 2014) to perform a tomographic reconstruction of theturbulent phase. A 589 nm dichroic is present in the opticalpath to avoid light contamination from the sodium AO lasers,so no scientific information is available around this wavelength.Let us also note that MUSE is undersampled (sampling is 25mas in NFM) on the whole available spectrum. Our PSF modelmanages the undersampling issue by oversampling the PSD andthe OTF to safely perform numerical computations. The givenPSF is then spatially binned to retrieve the correct sampling.The shape of the PSF can be retrieved, but this method haslower precision on the parameters’ estimation inherent to under-sampled data. These di ff erences between the MUSE instrumentand SPHERE / ZIMPOL allow us to test the versatility of ourPSF model. Additionally the spectral resolution is an asset tovalidate our model at di ff erent wavelengths and to study thespectral evolution of our PSF parameters.During the May-June 2018 commissioning phase, MUSEobserved multiple targets in narrow-field mode. Among thesetargets, we have access to 40 PSFs observed on di ff erent stars,during di ff erent nights and at di ff erent seeing conditions. Theseselected PSFs have been spectrally binned into 92 bins of 5nmeach to increase the signal to noise ratio and reduce the numberof fittings. Then fitting is performed independently, spectralbin by spectral bin, without any spectral information on thetargets or the atmosphere. Figure 9 shows one MUSE datacubePSF fitting at three di ff erent wavelengths. The evolution of theAO correction radius is clearly visible in both the data and themodel. As for ZIMPOL, we did not take into account static PSF(except the occulted pupil di ff raction), which is the main visibledi ff erence between data and model. Fainter stars visible in thefield did not a ff ect the fitting and appear clearly in the residuals.For the 40 datacubes PSF, the relative error is (cid:15) h = . × − .This result is similar to the previous ones on OOMAO andSPHERE / ZIMPOL. Secondary stars in the field also count inthe residual error computation.The evolution of r with the wavelength is shown on Fig. 10for one datacube PSF. A least square fitting between our datapoints and the theoretical λ / evolution of the r gives a spectral Fig. 8.
Zenital r at 500 nm estimated on MUSE (filled circles, 40 datapoints) and SPHERE / ZIMPOL (empty circles, 27 data points) usingthree methods: PSF fitting, SPARTA, and MASS / DIMM. Top: r , SPARTA versus r , FIT . A di ff erent linear tendency is found for MUSE (plain line)and SPHERE / ZIMPOL (dashed line). Shaded areas show the standarddeviation between data points and the best linear fit. Middle: r , SPARTA and r , FIT versus r , MASS / DIMM . Two linear tendencies are identified forSPARTA, and only one for PSF fitting. Bottom: Histograms of seeingestimated with the three methods on both instruments. Dashed verticallines show median values of 0 .
46” (SPARTA), 0 .
69” (MASS / DIMM),and 0 .
83” (PSF fitting). averaged estimation of r = . r matches well the theory, with a standard deviation of 0 . Article number, page 8 of 11. JL. Fétick et al.: Physics based model of the AO-corrected PSF
Fig. 9.
MUSE PSFs (top), fittings (middle), and residuals (bottom) ofthe same star at three di ff erent wavelengths over the 92 spectral bins ac-tually fitted. The hyperbolic arcsine of the intensity is shown to enhancedetails. Fig. 10.
Estimation of the r ff erentwavelengths (blue dots), and comparison with a theoretical law in λ / (orange line). The best match between data and theory is achieved for r = . So far each spectral bin is fitted independently, however thespectral deterministic trend we recover is an asset for PSF de-termination. It makes possible the fitting of the whole datacubewith only one r parameter. The statistical contrast – ratio of thenumber of measurements over the number of unknowns – wouldbe increased. It would improve fitting robustness, especially forfaint stars where the halo is strongly a ff ected by noise.The results of fitting on the 40 datacubes give statistical in-formation on the r estimation. As for the SPHERE / ZIMPOLcase, we have access to SPARTA and MASS / DIMM data tocorrelate with our fitting estimations (see Fig. 8). The Pear-son correlation coe ffi cient between PSF fitting and SPARTA is C Pearson = .
96, which is similar to the SPHERE / ZIMPOL case.However, we get the linear relationship r , SPARTA = . r , FIT − . , (24)which is di ff erent from the SPHERE / ZIMPOL (Eq. (23)). Theexact origin of this di ff erent trend is unknown, nevertheless letus note that the actual implementation of the phase PSD esti-mation is slightly di ff erent for the SPHERE and for the MUSEinstruments. For SPHERE the Von-Kármán outer scale L is setto 25 m, whereas for MUSE the L is estimated jointly with r .This SPARTA double trend is corroborated by MASS / DIMMinformation (Fig. 8, middle graph). On the other hand, r es-timation using our PSF model gives similar results on bothSPHERE / ZIMPOL and MUSE instruments. This confirms therobustness of our PSF fitting method.
4. Conclusions
In this article we developed a parametric model of long-exposure AO-corrected PSF. The particularity of this modelis to parameterize the phase PSD using a Mo ff at core anda turbulent Kolmogorov halo. This model also incorporatesprior knowledge of the telescope, such as the optical cuto ff frequency, the obstruction and spider shapes, and even the staticaberrations if they are calibrated, for example by phase diversity(Mugnier et al. 2008). This model only requires five parametersfor circularly symmetrical PSFs, and seven for asymmetricalones. The sparsity of this PSF model makes it suitable fornumerical computation, such as minimization algorithms orleast-square fits. Tests on both simulated and real data validatedthe appropriateness of our model.One substantial advantage of our model over focal planemodels is to use physical parameters such as the Fried parameter r and the residual AO variance σ . Since these parametersare physical, their values in our PSF model can be correlated toexternal measurements. Tests on both OOMAO simulations andon-sky data (from the SPHERE / ZIMPOL and MUSE instru-ments) confirmed the physical meaning of the r parameter usedin our PSF. The ultimate goal would be to only use physicalparameters in the PSF description.Our model has already shown usability for di ff erent seeingsand di ff erent instruments, with di ff erent AO-correction quality.This shows the robustness and versatility of the model. We alsoplan to use it to parameterize the PSF for the future instrumentson bigger telescopes such as the Extremely Large Telescope(ELT).Finally, the small number of parameters makes this modelsuited for image post-processing techniques such as deconvolu-tion of long-exposure images. Deconvolution using parametricPSFs has already been demonstrated by Drummond (1998) andFétick et al. (2019). We plan to develop a myopic deconvolutionalgorithm estimating both the observed object and the PSF pa-rameters in a marginal approach similar to Blanco & Mugnier(2011). Acknowledgements.
This work was supported by the French Direction Généralede l’Armement (DGA) and Aix-Marseille Université (AMU). This work wassupported by the Action Spécifique Haute Résolution Angulaire (ASHRA) ofCNRS / INSU co-funded by CNES. This project has received funding from theEuropean Union’s Horizon 2020 research and innovation program under grantagreement No 730890. This material reflects only the authors views and the
Article number, page 9 of 11 & A proofs: manuscript no. aanda
Commission is not liable for any use that may be made of the information con-tained therein. This study has been partly funded by the French Aerospace Lab(ONERA) in the frame of the VASCO Research Project.The authors are thankful to Pierre Vernazza for providing data related to his ESOLarge Program ID 199.C-0074.
References
Andersen, D. R., Stoesz, J., Morris, S., et al. 2006, Publications of the Astro-nomical Society of the Pacific, 118, 1574Ascenso, J., Neichel, B., Silva, M., Fusco, T., & Garcia, P. 2015, in AdaptiveOptics for Extremely Large Telescopes IV (AO4ELT4), E1Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Ground-based and AirborneInstrumentation for Astronomy III, Vol. 7735, International Society for Opticsand Photonics, 773508Bacon, R., Bauer, S., Böhm, P., et al. 2006, in Ground-based and Airborne In-strumentation for Astronomy, Vol. 6269, International Society for Optics andPhotonics, 62690JBeuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, in Ground-based and airborneinstrumentation for astronomy II, Vol. 7014, International Society for Opticsand Photonics, 701418Beuzit, J.-L., Vigan, A., Mouillet, D., et al. 2019, arXiv e-prints[ arXiv:1902.04080 ]Blanco, L. & Mugnier, L. M. 2011, Optics Express, 19, 23227Calia, D. B., Hackenberg, W., Holzlöhner, R., Lewis, S., & Pfrommer, T. 2014,Advanced Optical Technologies, 3, 345Conan, J.-M. 1994, PhD thesis, Université Paris XI OrsayConan, R. & Correia, C. 2014, in Proc. SPIE, Vol. 9148, Adaptive Optics Sys-tems IV, 91486CDavies, R. & Kasper, M. 2012, Annual Review of Astronomy and Astrophysics,50, 305Drummond, J. D. 1998, in Proc. SPIE, Vol. 3353, Adaptive Optical System Tech-nologies, ed. D. Bonaccini & R. K. Tyson, 1030–1037Fedrigo, E., Donaldson, R., Soenke, C., et al. 2006, in Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series, Vol. 6272, 627210Fétick, R. J., Jorda, L., Vernazza, P., et al. 2019, A&A, 623, A6Fétick, R. J. L., Neichel, B., Mugnier, L. M., Montmerle-Bonnefois, A., & Fusco,T. 2018, MNRAS, 481, 5210Fusco, T., Sauvage, J.-F., Petit, C., et al. 2014, in Adaptive Optics Systems IV,Vol. 9148, International Society for Optics and Photonics, 91481UGoodman, J. W. 1968, Appendix B Sec. B, 3Jolissaint, L. & Veran, J.-P. 2002, in European Southern Observatory Conferenceand Workshop Proceedings, Vol. 58, European Southern Observatory Confer-ence and Workshop Proceedings, ed. E. Vernet, R. Ragazzoni, S. Esposito, &N. Hubin, 201Jones, E., Oliphant, T., Peterson, P., & Others. 2001, SciPy: Open source scien-tific tools for PythonMartin, O. A., Gendron, É., Rousset, G., et al. 2017, A&A, 598, A37Milli, J., Mouillet, D., Fusco, T., et al. 2017, arXiv e-prints[ arXiv:1710.05417 ]Mo ff at, A. 1969, Astronomy and Astrophysics, 3, 455Mugnier, L. M., Fusco, T., & Conan, J.-M. 2004, Journal of the Optical Societyof America A, 21, 1841Mugnier, L. M., Sauvage, J. F., Fusco, T., Cornia, A., & Dandy, S. 2008, OpticsExpress, 16, 18406Müller Sánchez, F., Davies, R. I., Eisenhauer, F., et al. 2006, A&A, 454, 481N’Diaye, M., Dohlen, K., Fusco, T., & Paul, B. 2013, Astronomy & Astro-physics, 555, A94Orban de Xivry, G., Rabien, S., Busoni, L., et al. 2015, in Adaptive Optics forExtremely Large Telescopes IV (AO4ELT4), E72Racine, R., Walker, G. A., Nadeau, D., Doyon, R., & Marois, C. 1999, Publica-tions of the Astronomical Society of the Pacific, 111, 587Ragland, S., Dupuy, T. J., Jolissaint, L., et al. 2018, in Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series, Vol. 10703, AdaptiveOptics Systems VI, 107031JRigaut, F. J., Veran, J.-P., & Lai, O. 1998, in Society of Photo-Optical Instru-mentation Engineers (SPIE) Conference Series, Vol. 3353, Adaptive OpticalSystem Technologies, ed. D. Bonaccini & R. K. Tyson, 1038–1048Roddier, F. 1981, Progress in Optics, 19, 281Roddier, F. 1999, Adaptive optics in astronomy (Cambridge university press)Rusu, C. E., Oguri, M., Minowa, Y., et al. 2016, MNRAS, 458, 2Sauvage, J.-F., Fusco, T., Petit, C., et al. 2010, in Adaptive Optics Systems II,Vol. 7736, International Society for Optics and Photonics, 77360FSchmid, H. M., Bazzon, A., Roelfsema, R., et al. 2018, A&A, 619, A9Ströbele, S., La Penna, P., Arsenault, R., et al. 2012, in Proc. SPIE, Vol. 8447,Adaptive Optics Systems III, 844737Vernazza, P., Brož, M., Drouard, A., et al. 2018, A&A, 618, A154Viikinkoski, M., Vernazza, P., Hanuš, J., et al. 2018, A&A, 619, L3Zieleniewski, S. & Thatte, N. 2013, in Proceedings of the Third AO4ELT Con-ference Appendix A: Integral of a truncated Moffat
A Mo ff at function as given in Eq. (2) shows elliptical con-tours (see Fig. A.1). In this appendix we compute the integralof the Mo ff at function inside one of these elliptical contours,called E ( R x , R y ), of semi-major axis R x and semi-minor axis R y = ( α y /α x ) R x . The integral to calculate is I ( R x , R y ) = (cid:34) E ( R x , R y ) M A ( x , y )dxdy . (A.1)Let us perform the change of variables φ : (cid:40) R ∗ + × [ − π, π [ −→ R × R \ (0 , r , θ ) −→ ( α x r cos θ, α y r sin θ ) . (A.2)The determinant of the Jacobian isdet J φ = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) α x cos θ − r α x sin θα y sin θ r α y cos θ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = α x α y r . (A.3)The integral of the Mo ff at in the ellipse is rewritten as I ( R x , R y ) = (cid:90) π − π (cid:90) R x /α x | det J φ | M A ( φ ( r , θ ))drd θ = π A α x α y (cid:90) R x /α x [1 + r ] − β = A πα x α y β − (cid:26) − (cid:104) + ( R x /α x ) (cid:105) − β (cid:27) , (A.4)where we assumed β (cid:44) α x (cid:39) α y and R x (cid:39) R y so theintegral over the ellipse E ( R x , R y ) is nearly equal to the energyin a disk of radius R x . This assumption is made for our modelof PSD to compute the residual variance below the AO cuto ff frequency in Eq. (11). Fig. A.1.
Visualization of an elongated Mo ff at (colour map), its ellipti-cal level curves (white), and the ellipse E ( R x , R y ) inside which the inte-gral is computed. Dashed circles of radius R x and R y help to visualizethe error done when computing the integral over the ellipse instead of acircle. A very elongated Mo ff at is shown here ( α x = . α y ). Moreover it is possible to calculate the integral of the Mo ff atfunction over the whole plane. The assumption β > R x → + ∞ and R y → + ∞ , one finds I ( + ∞ , + ∞ ) = A πα x α y β − . (A.5) Article number, page 10 of 11. JL. Fétick et al.: Physics based model of the AO-corrected PSF
Consequently, in order to get a Mo ff at of unit integral over thewhole plane, one must choose the Mo ff at amplitude factor as A = β − πα x α y . (A.6) Appendix B: Analytic solution for the flux andbackground
The minimum of L (given in Eq. 15) has an analytic solution for γ and ζ since nulling the partial derivative of L towards theseparameters gives ∂ L ∂γ = ⇐⇒ (cid:88) i , j w i , j h i , j (cid:104) γ · h i , j + ζ − d i , j (cid:105) = ⇐⇒ γ (cid:88) i , j w i , j h i , j + ζ (cid:88) i , j w i , j h i , j = (cid:88) i , j w i , j h i , j d i , j (B.1)and ∂ L ∂ζ = ⇐⇒ (cid:88) i , j w i , j (cid:104) γ · h i , j + ζ − d i , j (cid:105) = ⇐⇒ γ (cid:88) i , j w i , j h i , j + ζ (cid:88) i , j w i , j = (cid:88) i , j w i , j d i , j . (B.2)These two equations are linear in γ and ζ . They can be writtenwithin the matrix formalism A · (cid:32) ζγ (cid:33) = (cid:88) i , j w i , j d i , j (cid:32) h i , j (cid:33) , (B.3)where A is the 2 × A = (cid:88) i , j w i , j (cid:32) h i , j h i , j h i , j (cid:33) . (B.4)In order to invert A , we need to make sure that det( A ) (cid:44)
0. Thedeterminant of A isdet( A ) = (cid:88) i , j w i , j (cid:88) i , j w i , j h i , j − (cid:88) i , j w i , j h i , j = (cid:88) i , j √ w i , j (cid:88) i , j ( √ w i , j h i , j ) − (cid:88) i , j ( √ w i , j )( √ w i , j h i , j ) . (B.5)We can now define the two vectors √ W = { √ w i , j } and √ WH = { √ w i , j h i , j } . Using these notations the determinant is rewritten asdet( A ) = (cid:107) √ W (cid:107) (cid:107) √ WH (cid:107) − |(cid:104) √ W , √ WH (cid:105)| . (B.6)We recognize the Cauchy-Schwarz inequality. It follows thatdet( A ) ≥
0, with equality only if √ W and √ WH are colinearvectors. The colinearity is written as √ W // √ WH ⇐⇒ ∃ k ∈ R | ∀ ( i , j ) , k √ w i , j = √ w i , j h i , j ⇐⇒ ∃ k ∈ R | ∀ ( i , j ) where w i , j (cid:44) , k = h i , j . (B.7)This states that the determinant of A is null only if the PSFmodel h is constant on each point where w i , j (cid:44)
0. Since our PSFis not constant on the domain w i , j (cid:44)
0, we ensure that det( A ) (cid:44) γ and ζ is written as (cid:32) ζγ (cid:33) = A − (cid:88) i , j w i , j d i , j (cid:32) h i , j (cid:33) . (B.8)(B.8)