An accurate centroid algorithm for PSF reconstruction
Tianhuan Lu Wentao Luo Jun Zhang Jiajun Zhang Hekun Li Fuyu Dong Yingke Li Dezi Liu Liping Fu Guoliang Li Zuhui Fan
astro-ph.IM








1 Department of Physics and Astronomy, Shanghai Jiao Tong University, 
  Shanghai, 200240, China; E-mail:
  mailto:njpipeorgan@sjtu.edu.cnnjpipeorgan@sjtu.edu.cn, 
  mailto:wentao.luo82@sjtu.edu.cnwentao.luo82@sjtu.edu.cn
  
2 Department of Astronomy, School of Physics, Peking University, 
  Beijing 100871, China

3 The Shanghai Key Lab for Astrophysics, Shanghai Normal University, 
  Shanghai 200234, China

4 Purple Mountain Observatory, Chinese Academy of Sciences, 
  Nanjing, 210000, China  

  


  In this work, we present a novel centroiding method based on Fourier space Phase Fitting (FPF) for Point
  Spread Function (PSF) reconstruction. We generate two sets of simulations to test our method.
  The first set is generated by GalSim with elliptical Moffat profile and strong 
  anisotropy which shifts the center of the PSF. The second set of simulation is drawn from CFHT i band stellar imaging data. We find non-negligible anisotropy
  from CFHT stellar images, which leads to ∼ 0.08 scatter in unit of pixels using polynomial fitting method <cit.>.  And we apply FPF method to 
  estimate the centroid in real space, this scatter reduces to ∼ 0.04 in SNR=200 CFHT-like sample.   In low SNR (50 and 100) CFHT-like samples, the background
  noise dominates the shifting of the centroid, therefore the scatter estimated from different methods are similar. We compare polynomial fitting and FPF using GalSim simulation with optical anisotropy. 
  We find that in all SNR (50, 100 and 200) samples, FPF performs better than polynomial fitting by a factor of ∼ 3. In general, we 
  suggest that in real observations there are anisotropy which shift the centroid, and FPF method is a better way to accurately locate it.
  









§ INTRODUCTION


Point Spread Function (PSF) is one of the major systematics in weak lensing
measurement. It introduces both multiplicative bias and additive bias. 
There are numerous methods in literature devoted to 
correcting PSF effects
<cit.>.  
  Lensfit <cit.> applies a
Bayesian based model-fitting approach; BFD (Bayesian Fourier Domain)
method <cit.> carries out Bayesian analysis in the
Fourier domain, using the distribution of unlensed galaxy moments as
a prior, and the Fourier_Quad method developed by <cit.> uses image moments in the Fourier Domain.

Many simulations are generated to test the accuracy of various methods, e.g. STEP (Shear TEsting Program)
<cit.>, Great08 <cit.>, Great 10
<cit.>, GREAT3 <cit.> or Kaggle – the
dark matter mapping competition[Supported by NASA & the
  Royal Astronomical Society.]. Other independent softwares, such as 
SHERA <cit.>, have also been designed
for specific surveys. But most of those simulations assume that the PSF
is perfectly known, which is not the case in reality. PSF at the position
of galaxy must be reconstructed using nearby star images.
  
In GREAT10 star challenge <cit.>, multiple PSF reconstruction methods has been 
tested, e.g. PSFEx <cit.>, PCA+Krigging <cit.>, Gaussianlets <cit.>,
B-slpline <cit.>, Inverse Distance Weighting (IDW), Radial Basis 
Function(RBF), and Krigging <cit.> etc. 
Especially, <cit.> tested various interpolation
methods to interpolate the PSF power sptectrum for Fourier_Quad shear 
estimator <cit.>, which achieves <1% level accuracy in GREAT3 simulation.

 As for estimating the centroid of
stellar images, a recent work by <cit.> claims that simple polynomial fitting works very well and close
to saturate the Cramér-Rao lower bound, while moment-based method 
does not deliver reliable centroid estimation. 
Our method, though, based on fitting the phase slope in Fourier space, not only
provide better centroid estimation in terms of scatter, but also automatically
shift the centroid to the center of a postage-stamp image after inverse Fourier
transformation.

We describe our method along with polynomial fitting,  in Sec. <ref>. In Sec. <ref>, we describe the simulations to 
test our method. The results are shown in Sec. <ref>. We summarize and conclude in
Sec. <ref>.



§ METHOD

 


In this section, we describe our centeroid measurement 
method along with polynomial fitting
method with a Gaussian smooth kernel.



 §.§ Fourier space Phase Fitting (FPF)



Given an image I() and its centroid 𝐱_𝐜. The Fourier transformation
of the image is simply

    f()=1/2π∬ I() e^-i· ,

where f(𝐤) can be written as 

    f()=|f()| e^i∠ f(),

with ∠ (f) as the phase. Given that the centroid of a noise free PSF image I() is defined as 

    :=∬ I() /∬ I() ,

we find that the 
centroid  in real space corresponds to the slope of phase near =0, i.e.

    =-∇∠ f(0).

The proof is shown in the appendix.

In the case that the image is symmetric about the centroid, which is usually the case in the vicinity
of the centroid and the image value is real. Then according to the Fourier transformation properties
we can deduce that

    f()=f(-)=f^*().

As a result, f() is also a real function, meaning the imaginary part vanishes and
the phase function is zero. Any anisotropy can further introduce none
zero imaginary part, which can be reflected by the phase.









Fig. <ref> gives an example of how off-center affect the phase pattern in Fourier space.
The phase will be zero if the image is perfectly centered, while there will be stripe pattern
caused by off centering effect.

After discretization, the derivative becomes a minimization problem: 


    χ^2 = ∑_w()[∠(f() e^-i ϕ())]^2

where f()= ()+()  is the Fourier transform of the PSF (()) 
with noise (()), and ϕ() is the modeled phase pattern of f(), which depends on the PSF (f()), the centroid ().
The observed image is a matrix with real values, so its Fourier transformation 
satisfies the conditions that are f(-)=f^*() and 
ϕ(-)=-ϕ().

For a symmetric PSF. the phase is simply related with the centroids linearly, 

    ϕ=k_x x_c+k_y y_c.

Nevertheless, the higher order anisotropy that shifts centroids, should be fitted 
with higher order(here we apply 3rd order) polynomial
to capture this effect,

    ϕ=k_x x_c+k_y y_c+α_1 k_x^3+α_2 k_x^2 k_y+α_2 k_x k_y^2+α_4 k_y^3,
 
where {α_1, α_2, α_3, α_4} are free parameters. 

In the appendix, We show that the weights bearing the analytical form

    w()=
    2 |f()|^2/|()|^2, 
       |()|≪ |()|
    
    0,    |()|≈|()|

can out-weight the noise while preserve as many data points as possible for the fitting.





 §.§ Polynomial fitting method



We follow <cit.> and apply a fixed Gaussian kernel 

    k(x)=1/2π w^2exp(-x^2/2w^2).

where w is 1.2 pixels given that the PSF FWHM is 2.8 pixels, to smooth 
the image before fitting the centroids with 2D polynomial

    P(x,y)=a+bx+cy+dx^2+exy+fy^2.

Only the central 3× 3 patch around the brightest pixels are 
used to solve the coefficients X={a,b,c,d,e,f}. The design matrix
A can be constructed as

    A=
    [       1     x_1     y_1   x_1^2 x_1 y_1   y_1^2;       1     x_2     y_2   x_2^2 x_2 y_2   y_2^2;       ⋮       ⋮       ⋮       ⋮       ⋮       ⋮;       1     x_9     y_9   x_9^2 x_9 y_9   y_9^2 ].


Then the coefficients can be determined by solving

    X=(A^T A)^-1 A^T Z

where Z is the pixels of the flattened 3× 3 patch. 
The centroids can be then determined by those coefficients.

    [ x_c; y_c ]
    =D^-1[ -b; -c ],

where  D=[ 2d  e;  e 2f ], is the curvature matrix,





§ SIMULATIONS


We simulated two sets of images, one uses GalSim <cit.> and
the other is based on Principal Components (PCs) decomposed from CFHT w2 stellar images.
The GalSim simulation provides
a set of optical effects.
The CFHT simulation is for exploring how much shifts caused by the anisotropy in real
surveys.



 §.§ GalSim stellar image


We use GalSim to simulate 3 sets of star images with different signal to noise ratio: SNR∼50, SNR∼100 and SNR∼200. The SNR is not strictly 50, 100 or 200 due to the fact that
we simulate the images using Exposure Time Calculator (ETC) from a uniform magnitude distribution centered at 
each SNR.



 We apply Moffat profile for all the stellar images, and
uniform distribution for the r_d and β in 

    I(r)=(1+r^2/r_d^2)^-β,
 
where I(r) is the 2D brightness distribution. 
After we generate the stellar images, we shift
the centers using two uniform distributions from -0.5 to 0.5 to each image as input centroid value.  Then
we convolve the images with optical anisotropy effect coma as shown in Fig. 2. We exaggerate this effect for better illustration, in fact, it is hard to be noticed by visual inspection. Finally we add Gaussian noise using
GalSim to simulate difference SNR stellar images. Noise and randomly oriented optical 
anisotropy contribute more dispersion to the centroid of the image.





 §.§ CFHT w2 stellar image


 In real observations, atmospheric seeing dilutes the observed objects and
introduce extra ellipticity and shift on centroid. The background sky also shifts the 
centroid randomly. For high SNR images, the centroid-shifting mechanism is dominated by
atmospheric seeing. An accurate centroid estimation method should be designed to capture
this shifting and correct for it.




In this simulation, we focus on testing this high order centroid shifting effect
based on real data from CFHTLens survey field w2, which contains 7 exposures,
with 10 minutes exposure for each. 




The procedure of this simulation is described as follows:

  
  * Select all the stellar images with SNR>100 and without saturation from CFHT w2 area, there are
   ∼ 600,000 stars in total.
  
  * Extract components using PCA without centroiding.
  
  * Generate 10,000 stellar images using the first 16 PCs, the coefficients of the PCs are randomly drawn from the parent distribution of the original images.
  
  * Calculate the centroid using brightness weighted moments using 
   
    _0=∬ I() d/∬ I() d .
 to be the reference.
   This can be considered as the real input center when there is no noise.


The first 6 PCs decomposed from CFHT w2 stellar images are shown in Fig. <ref>. We preserve the dipoles by not centroiding
the postage stamp images so that the asymmetry can be directly reconstructed by the dipoles. 




We display the images from two simulations in Fig. <ref>.
 The right panel is the PSF from GalSim simulation and the left one is 
CFHT-like simulation. Despite of the anisotropy effect added to GalSim PSFs, we still
can not observe its existence by visual inspection.
 
 
 



§ RESULTS





We demonstrate the performance of three centroiding
methods Fig. <ref>, i.e. polynomial fitting, 1st order phase fitting and 3rd order phase fitting based on two simulations. The top three panels show the results from GalSim simulation, the bottom three are from CFHT-like simulation. From left to right, the SNR=50,100,200 respectively. 

In general, 1st order FPF is already better than polynomial fitting in both simulations, and
3rd FPF is significantly improves the centroid estimation. In CFHT-like SNR=50 simulation, 
all three methods perform similarly due to the fact that the sky background noise dominates the
scatter budget. The quantitative comparison is listed in Table. <ref>.




In GalSim simulation, where the optical anisotropy dominates
the scatter budget, the scatter from 3rd order FPF is ∼3 times smaller
than that from polynomial fitting in all SNR branches.

In CFHT-like simulation, the performance of three methods are similar in SNR=50 branch.
For SNR=100 branch, which are often used in real analysis, the scatter from 3rd order FPF is ∼1.5 smaller than that from polynomial fitting. This difference expands to ∼2.0 for SNR=200 branch. 



In the CFHT-like simulation, the performance of three methods are similar in SNR=50 branch. As we increase SNR to 100, which is closer to those in real measurements, the scatter from 3rd order FPF is ∼1.5 smaller than that from polynomial fitting. This difference expands to ∼2.0 for SNR=200 branch. As the effect of noise drops, 3rd order FPF has a better and better performances relative to polynomial fitting, due to its ability to accurately capture higher order anisotropy features. 

It can be further illustrated by the GalSim simulation, where higher order anisotropy dominates the scatter budget. In all SNR branches, the scatter from 3rd order FPF is ∼3 times smaller than that from polynomial fitting.




§ SUMMARY AND DISCUSSION



Centroiding is the first and important procedure for 
PSF reconstruction. An accurate PSF reconstruction
further affects shear measurement. We develop
our centroid estimation method in Fourier space – 3rd order FPF. 



In GalSim simulation, the centroiding shift is dominated by optical anisotropy. The scatter of 3rd order FPF are smaller than polynomial fitting
method by a factor of ∼2 to  ∼3,  from SNR=50 to SNR=200 branch. 

In CFHT-like simulation, where the higher order anisotropy 
is much smaller than GalSim simulation, we found that for the SNR=50
images, background noise dominate the scatter, while
for the SNR=200, optical anisotropy play a major role for this scatter. Therefore
3rd order FPF performs similarly in SNR=50 branch, but with half the scatter of polynomial
fitting method in SNR=200 branch.


Therefore, we conclude that 3rd order FPF method is so far the most accurate estimation in
centroiding. The scatter caused by noise can not well corrected for any methods, but for 
scatter introduced by optical anisotropies, 3rd order FPF can capture the shift precisely. 
This is very important for weak lensing measurements.




This work was supported
by the following programs;  NSFC
(Nos. 11503064), Shanghai Natural Science Foundation, Grant
No. 15ZR1446700. JZ  is  supported  by  the  NSFC  grants  (11673016, 11433001, 11621303) and the National Key Basic Research Program
of China (2015CB857001). L.P.F. acknowledges the support from NSFC grant 11333001 and
11673018, STCSM grants 13JC1404400 and 16R1424800, SHNU grant DYL201603.

This work was also supported by the High Performance Computing
Resource in the Core Facility for Advanced Research Computing at
Shanghai Astronomical Observatory.
















[Bacon & Taylor(2003)]Bacon2003
 Bacon, D. J., & Taylor, A. N. 2003, , 344, 1307




[Bergé et al.(2012)]Berge2012 Bergé, 
  J., Price, S., Amara, A., & Rhodes, J. 2012, , 419, 2356 




[Bertin(2011)]Bertin2011 Bertin, E. 2011, 
    Astronomical Data Analysis Software and Systems XX, 442, 435 


 [Bernstein & Jarvis(2002)]Bernstein2002
   Bernstein, G. M., & Jarvis, M. 2002, , 123, 583

 [Bernstein(2009)]Bernstein2009 Bernstein, G. M. 2009,
   , 695, 652

 [Bernstein & Armstrong(2014)]Bernstein2014
  Bernstein, G. M., & Armstrong, R. 2014, , 438, 1880



 
[Bertin & Arnouts(1996)]Bertin1996
 Bertin, E., & Arnouts, S. 1996, , 117, 393









 [Bridle et al.(2002)]Bridle2002 Bridle, S. L., Kneib,
J.-P., Bardeau, S.,
& Gull, S. F. 2002, The Shapes of Galaxies and their Dark Halos, 38

 [Bridle et al.(2009)]Bridle2009 Bridle, S.,
Shawe-Taylor, J., Amara, A., et al. 2009, Annals of Applied Statistics, 3,6












































  
[Gentile et al.(2013)]Gentile2013 Gentile, M., Courbin, F., 
   & Meylan, G. 2013, , 549, A1 













 [Heymans et al.(2004)]Heymans2004 Heymans, C., Brown, M.,
  Heavens, A., et al. 2004, , 347, 895

 [Heymans et al.(2005)]Heymans2005 Heymans, C., Brown,
  M. L., Barden, M., et al. 2005, , 361, 160

  [Heymans et al.(2006)]Heymans2006 Heymans, C., Van
   Waerbeke, L., Bacon, D., et al. 2006, , 368, 1323




 [Hirata & Seljak(2003)]Hirata2003
  Hirata, C., & Seljak, U. 2003, , 343, 459
  


















 [Kaiser et al.(1995)]Kaiser1995 Kaiser, N., Squires, G.,
   & Broadhurst, T. 1995, , 449, 460

  [Kaiser(2000)]Kaiser2000 Kaiser, N. 2000, , 537,
555




 [Kitching et al.(2008)]Kitching2008 Kitching, T. D.,
  Miller, L., Heymans, C. E., van Waerbeke, L.,
  & Heavens, A. F. 2008, , 390, 149

 [Kitching et al.(2010)]Kitching2010 Kitching, T., Balan,
S., Bernstein, G., et al. 2010, arXiv:1009.0779

 [Kilbinger et al.(2013)]Kilbinger2013 Kilbinger, M., Fu,
  L., Heymans, C., et al. 2013, , 430, 2200



























[Li, Xin, & Cui(2013)]LiXin2013 Li, G., Xin, B., 
  & Cui, W. 2013, Astrophysics from Antarctica, 288, 306 










  




[Lu et al.(2016)]Lu2016 Lu, T., Zhang, J., Dong, F., et al. 2016, arXiv:1610.09828 





[Luo et al.(2017)]Luo2017 Luo, W., Yang, X.,
& Zhang, J. 2017, , 836, 1








[Kitching et al.(2013)]Kitching2013 Kitching, T. D., 
 Rowe, B., Gill, M., et al. 2013, , 205, 12 






















[Mandelbaum et al.(2012)]Mandelbaum2012 Mandelbaum, R.,
Hirata, C. M., Leauthaud, A., Massey, R. J.,
& Rhodes, J. 2012, , 420, 1518




 [Mandelbaum et al.(2014)]Mandelbaum2014 Mandelbaum, R.,
  Rowe, B., Bosch, J., et al. 2014, , 212, 5

[Mandelbaum et al.(2015)]Mandelbaum2015 Mandelbaum, R., 
Rowe, B., Armstrong, R., et al. 2015, , 450, 2963 


 [Maoli et al.(2000)]Maoli2000 Maoli, R., Mellier, Y.,
  van Waerbeke, L., et al. 2000, The Messenger, 101, 10




 [Massey et al.(2007a)]Massey2007 Massey, R., Rhodes, J.,
  Leauthaud, A., et al. 2007a, , 172, 239










[Miller et al.(2007)]Miller2007 Miller, L., Kitching,
T. D., Heymans, C., Heavens, A. F., & van Waerbeke, L. 2007, , 382, 315

 [Miller et al.(2013)]Miller2013 Miller, L., Heymans, C.,
  Kitching, T. D., et al. 2013, , 429, 2858


























 [Refregier(2003)]Refregier2003 Refregier, A. 2003, , 41, 645
 
 [Refregier et al.(2010)]Refregier2010 Refregier, A., Amara, 
A., Kitching, T. D., et al. 2010, arXiv:1001.0061 

 [Rhodes et al.(2000)]Rhodes2000 Rhodes, J., Refregier,
A., & Groth, E. J. 2000, , 536, 79

 [Rhodes et al.(2007)]Rhodes2007 Rhodes, J. D., Massey,
  R. J., Albert, J., et al. 2007, , 172, 203

 [Rowe et al.(2015)]Rowe2015 Rowe, B. T. P.,
 Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121







































[Vakili & Hogg(2016)]Vakili2016 Vakili, M., 
& Hogg, D. W. 2016, arXiv:1610.05873 




[van Waerbeke(2001)]vanWaerbeke2001 van Waerbeke, L. 2001,
Cosmological Physics with Gravitational Lensing, 165









































  [Zhang (2010)]Zhang2010 Zhang, J. 2010, , 403,
 673

 [Zhang (2011)]Zhang2011 Zhang, J. 2011, JCAP, 11, 041


 [Zhang et al.(2015)]Zhang2015 Zhang, J., Luo, W.,
  & Foucaud, S. 2015, JCAP, 1, 024

[Zhang (2016)]Zhang2016 Zhang, J  2016, National Science Review, 3, 159-164













§ CENTROID OF A PSF IMAGE IN FOURIER SPACE


First, we calculate the derivative of f() with respect to k_x at =0: 

    .∂/∂ k_x|_=0 f()
       =   .∂/∂ k_x|_=0{1/2π∬ I() e^-i· }
       =   .∂/∂ k_x|_k_x=0{1/2π∬ I(x,y) e^-i k_x x  x y}
       =   1/2π∬ I(x,y) .∂/∂ k_x|_k_x=0 e^-i k_x x  x y 
       =   -i/2π∬ I(x,y) x  x y 
       =   -i x_c 1/2π∬ I(x,y)  x y, 
       =   -i f(0) x_c.

Next, we related the derivative of f() with respect to k_x with that of ∠ f(): 

    .∂/∂ k_x|_=0 f()
       =   .∂/∂ k_x|_=0{|f()| e^i∠ f()}
       =   |f(0)| e^i∠ f(0)(i.∂/∂ k_x|_=0∠ f())+∂ f(0)/∂ k_x e^i∠ f(0)
       =   i f(0).∂/∂ k_x|_=0∠ f()

Compare Eqn.(<ref>)  with Eqn.(<ref>), we find 

    x_c=-.∂/∂ k_x|_=0∠ f(),

and similarly, 

    y_c=-.∂/∂ k_y|_=0∠ f().

Thus, we get 

    =-∇∠ f(0).




§ WEIGHTS IN FOURIER PHASE FITTING

Suppose the Fourier transform of the noise free PSF is (), and the noise has magnitude |()| and random phase. Given f=+, we can deduce the scatter of ∠ f in terms of its variance. 

When ||≪|| (without loss of generality, we let ∈ℝ),  

    Var(∠ f)
       =   ⟨(∠ f-⟨∠ f⟩)^2⟩
       =   ⟨(∠(+)-∠)^2⟩
       =   ⟨(Im().∇_Im∠ f|_f=)^2⟩
       =   1/2||^2 (. f^-1/ f|_f=)^2 
       =   ||^2/2 ||^2,

where ∇_Im denotes the directional derivative along the imaginary part. 


The weights should be inversely proportional to the variances. We take

    w()=Var(∠ f)^-1=2 ||^2/||^2≈2 |f|^2/||^2.


When ||≈||, the scatters of ∠ f are so large that such pixels do not provide useful information of ∠. Thus we take 

    w()=0.





