IMAGE ANALYSIS FOR COSMOLOGY: RESULTS FROM THE GREAT10 STAR CHALLENGE

In this paper we present the results from the GRavitational lEnsing Accuracy Testing 2010 (GREAT10) Star Challenge. GREAT10 was an image analysis challenge for cosmology that focused on the task of measuring the shapes of distant galaxies. Light from distant galaxies is deflected during its journey to us via gravitational lensing, and the images appear distorted into characteristic patterns (Hu 1999; Bartelmann & Schneider 2001). The amount of distortion depends on the intervening distribution of matter (including dark matter) and the geometry of spacetime (which is currently governed by dark energy). Such measurements thus probe directly the invisible dark sector and the fundamental nature of gravity—see reviews by Albrecht et al. (2006), Réfrégier (2003), Hoekstra & Jain (2008), Massey et al. (2010), and Weinberg et al. (2012).

All real imaging data are necessarily seen after convolution with (i.e., blurring by) a telescope's point-spread function (PSF). The PSF arises from the finite aperture of the telescope, charge diffusion within digital detectors, any imperfect elements along the optical path, and turbulence in the Earth's atmosphere (unless the telescope is in space). This increases the size of faint galaxies, and can spuriously change their ellipticity by an order of magnitude more than gravitational lensing (Bernstein & Jarvis 2002; Hoekstra 2004; Paulin-Henriksson et al. 2008, 2009; Massey et al. 2013). To recover the shape of the galaxy after only cosmological effects, it is necessary to (1) model the PSF and (2) somehow correct for its effect on the images of galaxies. The second half of this task has been widely addressed by teams analyzing individual surveys and as a vital community effort through the public Shear TEsting Programme (STEP) (Heymans et al. 2006; Massey et al. 2007), the GRavitational lEnsing Accuracy Testing (GREAT) galaxy challenges (Bridle et al. 2010; Kitching et al. 2012b), and the Mapping Dark Matter challenge (Kitching et al. 2012c). The first task (modeling the PSF) has so far only been investigated internally within teams (e.g., Bacon et al. 2003; Hoekstra et al. 2004; van Waerbeke et al. 2005; Rhodes et al. 2007; Schrabback et al. 2010; Hoekstra 2004; Rowe 2010; Jarvis & Jain 2005; Bergé et al. 2012). Here we present the results of the first blind, public trial of methods to model and interpolate the PSF of a typical astronomical telescope.

The PSF in an astronomical image can be measured from stars that happen to fall inside the field of view. Stars are so small that they are intrinsically point-like, and adopt the size and shape of the telescope's PSF. However, the PSF typically varies across the field of view, and stars only sparsely cover the extragalactic sky (Jarvis & Jain 2005; Jain et al. 2006; Heymans et al. 2012; Chang et al. 2012). It is therefore necessary to model the shapes of stars, then interpolate their shapes to the locations of galaxies (where there is necessarily not a bright star because otherwise the galaxy could not be seen). In practice the PSF also varies as a function of the wavelength of observed light, due to diffraction, reflection, and transmission effects in the telescope optics, filters, and CCDs and so must also be interpolated from the colors of the stars to the colors of the galaxy (Cypriano et al. 2010; Voigt et al. 2012; Plazas & Bernstein 2012). Color dependence is an important second order effect but in this paper we do not address this, focusing only on the primary changes in PSFs.

We simulated the spatial variation in the PSF of generic but realistic ground-, balloon-, and space-based telescopes (Kitching et al. 2012a, and see www.greatchallenges.info). We realized a large suite of sparse stellar fields in these different observing regimes, and publicly released most of the star images. Entrants were asked to then reconstruct the images of the missing stars on a pixel grid, at pre-defined locations. The performance of each entry was measured in real time using a single number, the "quality factor," which was designed to provide a crude ranking such that it could not be reverse-engineered to reveal the full solutions. In this paper, we analyze in detail the quantitative performance of 12 distinct algorithms submitted to model and interpolate the simulated PSFs. In particular we quantify how well the ellipticity and size of a spatially varying PSF can be reconstructed in a blind challenge.

This paper is organized as follows. In Section 2, we describe the simulations and competition in detail. In Section 3, we present results. We discuss and conclude in Section 4.

In this section we describe the simulations and the competition. For a full exposition of the background of the Star Challenge, see Kitching et al. (2012a).

2.1. Simulation Structure

2.2. Data Structure

2.3. Competition Structure

The competition started in 2010 December and ran for nine months until 2011 September; this was concurrent with the GREAT10 Galaxy Challenge (Kitching et al. 2012a, 2012b). As stated previously the total simulation size was ∼50 Gb and the total size of the uploaded submissions was ∼1 Gb (we allowed participants to tar, zip, or FITS-compress¹⁹ submissions to reduce size). Data and example code were provided online for participants.²⁰

The two parameters of the PSF that most directly impact the ability to interpret observations of galaxies are the ellipticity and the size of the PSF; any residual difference between the ellipticity or size of true PSF, and the respective quantities of the modeled PSF at any particular position, will result in errors and biases in parameters assigned to any galaxy at that position. Weak gravitational lensing is particularly sensitive to these types of error (Massey et al. 2013; Paulin-Henriksson et al. 2008, 2009). The ellipticity and size are defined here using the second order brightness moments of the image as

$$\begin{equation} q_{ij}=\frac{\sum _p w_p I_p (\theta _i - \bar{\theta }_i)(\theta _j - \bar{\theta }_j)}{\sum _p w_p I_p}, \quad i,j\in \lbrace 1,2\rbrace, \end{equation} \tag{ 1 }$$

where the sums are over pixels, I_p is the flux in the pth pixel, and θ is a pixel position (θ₁ = x_p and θ₂ = y_p). In order to make the results regular with regard to the impact of noise but not to constrain the interpretation to compact objects in the postage stamp, we include a weight function w_p chosen to be a broad Gaussian with a width of 24 pixels (we leave an investigation of how results vary as a function of weight for future work). These are almost unweighted quadrupole moments in this respect, and as a result, smooth analytical functions may be favored compared to models that try to reproduce details in the wings of the PSF. The weighted ellipticity (or technically the "polarizability") for a PSF in complex notation is defined as

$$\begin{eqnarray} e &=& \frac{q_{11} - q_{22} + 2{\rm i}q_{12}}{q_{11} + q_{22} + 2\left(q_{11}q_{22}-q^2_{12}\right)^{1/2}}, \end{eqnarray} \tag{ 2 }$$

where we have used a definition of ellipticity |e| = (1 − r)(1 + r)⁻¹, where r is the ratio of minor to major axes of the ellipse. For the weighted size we have a similar expression

$$\begin{equation} R^{2} = q_{11} + q_{22}. \end{equation} \tag{ 3 }$$

We can calculate the variance between the ellipticity of the model and true PSF $\sigma ^2(e)\equiv \langle (e - e_{\rm PSF}^{{\rm t}})^2\rangle$ and similarly for the size $\sigma ^2(R)\equiv \langle (R - R_{\rm PSF}^{{\rm t}})^2\rangle$. Submissions were scored in real-time on a leaderboard that displayed the metric P ≡ (1/1/2〈σ²(R) + σ²(e)〉) where the average was taken over images in each set only (the total number of asked star positions is N_Stars = N_images × N_{starsperimage}, the average was only taken over N_images), such that that a mean variance of 10⁻³ in both ellipticity and size would have P ∼ 1.0.

The P metric, while indicatively ranking the methods, does not offer any insight into the performance of a method on ellipticity and size reconstruction. In this paper we will present quantities that relate to the principal properties of the PSF more directly. These are the standard deviation of the mean of the residuals of the ellipticity σ(e) and size squared σ(R²)/R² over all asked stars, i.e., we compute the error on the mean of the residuals (the sample variance computed using centered second order moments). We assume that any mean bias could be removed through cross-validation, in this sense it is a generous analysis to those methods with a mean residual. We average these quantities over the 50 images in each set, but in fact for all methods we find that the fractional error between images in a set is ≲ 10%.

In total, 30 submissions were made from 9 teams. As a baseline benchmark, a method in which all stars were simply stacked in an image was created, where no spatial variation in the reconstructed stars was present. Several methods generated low scores due to misunderstanding of simulation details, resulting in scores below the benchmark, and in this paper we summarize only those not affected by these issues. In the following we choose the best performing submission, for size, for each of the 12 distinct method entries. All of the submitted methods are described in the Appendix. We show the results on the fiducial Airy set (set 1 in Table 1) and the fiducial Moffat set (set 8 in Table 1) in Tables 2 and 3, respectively. In Figures 1 and 2 we present general behaviors of methods over the sets as categories were changed, but for a quantitative presentation of each method we refer the reader to Figures 3–7 where we show pictographic tables of results. In the square root of the inverse variance, and referring to Figures 1 and 2, for example, the error between images in a set of 10% results in changes of Δ[1/σ(e)] ≈ ±200 and similarly for R², not being significant for most methods individually, but can be significant collectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall we find that the B-Splines, Inverse Distance Weighting (IDW), and Radial Basis Function (RBF) methods reconstruct the ellipticity and size most accurately (see Gentile et al. 2013), with σ(e) ≈ 2.5 × 10⁻⁴ and σ(R²)/R² ≈ 7.4 × 10⁻⁴ over all sets.²¹ We note, however, that this is a snapshot of performance and that further investigations into tunable aspects of code could result in improvements in all methods.

We summarize the behavior of the submissions below. In each test all other parameters are kept fixed except those discussed (with fiducial values of 1000 known star positions, no mask, and telescope parameters given in Section 2). We refer to Figures 1 and 2 which show the change in the square root of the inverse variance of the reconstructed PSFs over the fiducial sets (set 1 for Airy and set 8 for Moffat profiles; see Table 2) when each of the categories is varied. In Figures 3–7 we show pictographic tables of results.

• 1.  
  PSF type. For the best performing methods we find a trend that both ellipticity and size are estimated more accurately for the Airy PSF than for the Moffat PSF.
• 2.  
  Addition of Kolmogorov power. For each set combination where both Moffat and Moffat-plus-Kolmogorov power are available (e.g., the 4 arm + masks) we find evidence for methods performing less well with the addition of Kolmogorov power (see also Figures 3–5). In Figure 6 we also show the impact of adding a Kolmogorov power spectrum to a set that uses an Airy PSF profile. We find that the addition of this random component degrades the residual ellipticity reconstruction by a factor of ≳ 2–5, but has less impact on size reconstruction, as expected, since the power is in ellipticity only. These results will necessarily depend on the amplitude of the assumed power spectrum, which will vary for each ground-based telescope, and knowledge/information about this is improving (e.g., Heymans et al. 2012). In addition, atmospheric turbulence also changes the PSF size, but we do not simulate this here. It is possible that, depending on the site and weather, the impact of turbulence may be weaker or stronger than that simulated for this study.
• 3.  
  Masks. We show results for the mask variation in Figure 3. We find that for all methods the presence of diffraction spikes does not degrade the ability to measure the ellipticity of the PSF. For the Airy function the diffraction spikes act to increase the effective size of the PSF, which enables methods to measure the fractional error σ(R²)/R² more accurately; but note that for a fixed σ(R²) a large size will decrease the fractional error by definition. For the Moffat PSF the diffraction spikes impact the size estimation significantly. We note, however, that this was a simple addition of a mask with no commensurate change in the variation of ellipticity or size across a field of view. Also the diffraction spikes contained low flux (only observable with the eye if one stacked all stars); higher signal-to-noise stars would change this. We leave an investigation of these effects to future work.
• 4.  
  Number of stars. We find that all methods are only weakly dependent, or insensitive to the number of stars used to reconstruct the PSF in these simulations, except for those sets in which we include a Kolmogorov power spectrum where we find that a larger number of stars results in a better reconstruction for the best methods (see Figure 4). This indicates that PSFs with spatial power on smaller scales require more stars for a particular reconstruction accuracy than PSFs without power on small spatial scales.
• 5.  
  Size of PSF. For the Airy profile we find that the larger the PSF the more accurately its size can be measured; for the Moffat we find a weak dependence with size. This is understandable because a larger PSF is better sampled and hence the size is easier to measure. This statement applies for PSFs with sizes 1.5 and 3 pixels in this study, but is certainly not true in all generality: once the sampling is better than critical, other factors (e.g., noise) take over (which our results also support). We also stress that an increase of the size of the PSF relative to the apparent size of galaxies will cause the galaxies to be less well-resolved, losing information and placing greater demands on shape measurement (Paulin-Henriksson et al. 2008). Also with the simulations presented, the impact of sampling on weak lensing shape measurement was not tested, only the performances of the PSF interpolation methods. We show results for the PSF size variation in Figure 5. When trading requirements of PSF model residuals against requirements for resolution (i.e., the absolute size of the ellipticity and PSF) such behavior should be noted.
• 6.  
  Telescope parameters. We show results for the PSF size variation in Figure 7. In varying the telescope parameters in the Jarvis et al. (2008) model we change the fiducial parameters respectively (a₀ = 0.014, d₀ = −0.006, c₀ = −0.010) to a₀ = −0.011, d₀ = 0.009, and c₀ = −0.011, i.e., an opposite astigmatism, a positive de-focus, and a 10% increased coma. We find that methods in this experiment were not affected by the change in de-focus, but performed better with the change in these astigmatism and coma parameters.

We discuss each method individually in the Appendix.

This paper presents the first blind simulation challenge aimed at testing optical PSF reconstruction methods. Simulations were generated in which participants were presented with a spatially varying PSF, sparsely sampled by stars, and asked to reconstruct the PSF at non-star positions. The competition, the GREAT10 Star Challenge, attracted 30 submissions from 9 teams; several of these teams were from non-astronomy backgrounds. The simulation presented participants with 27,500 stars over 1300 images subdivided into 26 sets, where in each set a category change was made in the type or spatial variation of the PSF. The simulations were intentionally simplistic so as to present the problem in an approachable way; in particular, the spatial variation of the PSF and the form of the PSF use simple analytic functions. In addition, only spatial variation, not temporal variation, was tested; hence these results should not be used to make specific statements about any particular experiment but should provide a benchmark with which methods can be tested and improved.²²

In this paper we analyze the submissions by testing how well each one can measure the ellipticity and size of the PSF. We quantify this as the inverse variance in the modeled PSF in each image for ellipticity and sized squared—defined using weighted quadrupole moments. This study was motivated by a desire to find methods that will be of use for weak gravitational lensing, where the PSF must be reconstructed to high accuracy (Paulin-Henriksson et al. 2008, 2009) at galaxy positions, but these results should also be of more general interest for any science case that analyzes galaxy images with optical data.

The submissions, and this paper, present a snapshot of any methods' ability to model the PSF. Due to the nature of the competitive blind submissions, post-challenge tuning of methods, which may yield significant improvements for any given method over the results presented here (see Gentile et al. 2013, for example), were not investigated. Each method submitted is summarized in the Appendix. We can, however, make some general statements about regimes in which methods tend to perform well or poorly when run in a blind way.

The functional form of the PSF was either a Moffat function or an Airy function, the spatial variation of the PSF was modeled using the analytic function given in Jarvis et al. (2008). In addition we optionally included diffraction spikes (+ or * forms), changed the PSF size (from 3.0 pixels to 1.5 or 6.0 pixels), changed the number of stars (from 1000 to 500 or 2000), and added an atmospheric turbulence pattern in ellipticity (with a Kolmogorov power spectrum). To summarize the conclusions we find the following.

• 1.  
  The best methods can reconstruct the PSF with an accuracy of σ(e) ≈ 2.5 × 10⁻⁴ and σ(R²)/R² ≈ 7.4 × 10⁻⁴ over all sets.
• 2.  
  Methods that performed poorly did so in part because the functional form of the PSF was not modeled correctly (in particular the Airy function).
• 3.  
  Smaller PSFs were more difficult to model than larger PSFs for the Airy function, but we add a caution that this does not mean larger PSFs are better for weak lensing because information on a target object is lost; instead this means that well sampled PSFs are better for weak lensing.
• 4.  
  Diffraction spikes caused the size of Moffat PSFs to be modeled less accurately, but Airy PSFs more accurately, due to the increase in the effective size.
• 5.  
  The addition of atmospheric Kolmogorov power (equivalent to short exposure PSFs; see Heymans et al. 2012) made ellipticity and size reconstruction less accurate by a factor of ≳ 2–5 for all methods. We add the caveat that the temporal nature of varying PSFs was not investigated, therefore methods such as cross-correlation between sequential images, which could potentially improve modeling, were not investigated.

For subsequent blind PSF modeling challenges the realism of the temporal and wavelength dependent nature of PSF variation could be included, and the simulations could be tailored to specific experiments.

Modeling the PSF is of critical importance in efforts to understand the nature of dark energy and dark matter using weak gravitational lensing where any inaccuracy in the modeled PSF can cause biases, and increased errors of cosmological parameters of interest. To address this crucial open problem this initial presentation of a blind PSF reconstruction challenge will hopefully provide a benchmark upon which methods can continue to be refined and tested.

T.D.K. is supported by a Royal Society University Research Fellowship, and was supported by a Royal Astronomical Society 2010 Fellowship and the University of Edinburgh for some of this work. B.R. and C.H. acknowledge support from the European Research Council in the form of a Starting Grant with numbers 24067 (B.R.) and 240185 (C.H.). R.M. is supported by a Royal Society University Research Fellowship. D.G. was supported by SFB-Transregio 33 "The Dark Universe" by the Deutsche Forschungsgemeinschaft (DFG) and the DFG cluster of excellence "Origin and Structure of the Universe" and thanks Gary Bernstein and Stella Seitz for helpful discussions. M.Ge., G.C., and G.M. are supported by the Swiss National Science Foundation (SNSF). G.L. thanks Wei Cui for useful discussions. G.L. and B.X. were supported in part by the U.S. Department of Energy through Grant DE-FG02-91ER4068 and G.L. is also supported by the one-hundred talents program of the Chinese Academy of Sciences (CAS). M.K. thanks Liping Fu. This work was funded by a EU FP7 PASCAL 2 Challenge Grant. Workshops for the GREAT10 Challenge were funded by the eScience STFC Theme and by NASA JPL, and hosted at the eScience Institute Edinburgh and by IPAC Caltech, Pasadena. We thank Francesca Ziolkowska, Harry Teplitz, and Helene Seibly for local organization of the workshops. We thank the GREAT10 Advisory team, co-authors of the GREAT10 handbook (Kitching et al. 2012a), for discussions before and after the challenge.

Contributions: All authors contributed to the writing and analysis presented. T.D.K. was PI of GREAT10, defined and created the simulations, and lead the analysis. T.D.K., B.R., M.G., C.H., R.M. were active members of the GREAT10 team during 2010 December to 2011 September and after the challenge. B.R. created the FITS image simulation code. M.Ge., F.C., G.M., D.G., M.K., K.G., A.S., A.M., G.L., B.X. submitted entries to the GREAT10 Star Challenge. D.W. maintained the GREAT10 leaderboard and processed submissions with T.D.K. during the challenge.

Here we include a brief description and references for each of the methods submitted to the challenge.

Several methods use the name "Kriging," which is in fact the same method as Gaussian process regression, the method submitted for the methods MoffatGP; Kriging is a different term which has been in use in the geostatistics field but all are types of Gaussian processes.

A.1. PSFEx (Gruen)

A.2. PCA+Kriging (Li, Xin)

A.3. Gaussianlets (Li, Xin)

A.4. B-Splines (Gentile, Courbin, Meylan)

A.5. Inverse Distance Weighting (IDW) (Gentile, Courbin, Meylan)

A.6. Radial Basis Function (RBF and RBF-thin) (Gentile, Courbin, Meylan)

A.7. Kriging (Gentile, Courbin, Meylan)

A.8. IDWStk (Gentile, Courbin, Meylan)

A.9. MoffatGP (Georgatzis, Mariglis, Storkey)