DMTN-045: PSF Fitting: Literature Overview

  • Krzysztof Findeisen

Latest Revision: 2017-06-13

1   Overview

Understanding the point-spread function is important for several Alert Production tasks, including image differencing and photometry. The LSST point spread function is expected to vary in a complex way with chip position, telescope orientation and focus, and atmospheric conditions. The PSF must be fit independently for each image, and must be fit quickly enough to not form a bottleneck for Alert Production.

The requirements for PSF quality depend on downstream operations, but are in general poorly quantified. In particular:

  • the [AlardLupton98] image differencing algorithm requires that the FHWM of the PSF be known to within 10% (Reiss, priv. comm. 2017), but is otherwise self-contained
  • the ZOGY image differencing algorithm requires a “high quality” PSF (Reiss, priv. comm. 2017)
  • the limiting factor for our current DCR correction algorithm is noise in the images to be corrected, not errors in the PSF. Therefore, the only requirement DCR correction imposes is that the PSF not have color terms (Sullivan, priv. comm. 2017; [S+17])
  • PSF requirements specific to LSST photometry are not available, but PSF photometry is in general only moderately sensitive to PSF models ([AndersonBedinPiotto+06]), so photometry is unlikely to be the limiting application.
  • PSF requirements specific to LSST astrometry are not available

In addition, the Alert Pipeline PSF fitting code needs to be robust to crowded fields (Bosch, April 2017 PST Meeting)). Source crowding may be well in excess of source fill factor \(\Sigma n_\mathrm{eff} \sim 0.2\), where \(\Sigma\) is the source density and \(n_\mathrm{eff}\) is the PSF area defined in eq. 2 of [I+11]. The PSF fitting code may also need to support high-frequency atmospheric components ([ChangMarshallJernigan+12]), depending on the accuracy needed for image differencing.

This note reviews recent developments on PSF modeling, both in publicly available code and in the literature. Given the lack of general-purpose metrics for PSF fitting methods (new metrics will be presented in a future DMTN), the discussion is mostly qualitative. Section 2 covers recent software for PSF modeling of large astronomical datasets. Section 3 covers algorithms presented for ground-based observations of crowded fields. Section 4 covers algorithms designed for atmospheric variations and other high-frequency effects.

The following notation appears in discussions of an algorithm’s asymptotic complexity:

  • \(N_\mathrm{pixels}\): the number of pixels in the image
  • \(N_\mathrm{PSF}\): the number of pixels used to represent or fit the PSF
  • \(N_\mathrm{star}\): the number of stars used to fit the PSF (note: in general, \(N_\mathrm{star} N_\mathrm{PSF} \ll N_\mathrm{pixels}\))
  • \(N_\mathrm{cells}\): the number of cells used to model PSF spatial variations

2   General-Purpose Software

2.1   PSFEx

PSFEx is a generic PSF fitting tool last updated in 2014. While intended to be run as a standalone command-line program, the LSST Stack’s meas_extensions_psfex provides a wrapper that lets its core functionality be used as a PsfDeterminer.

When run as a command-line tool, PSFEx identifies PSF template sources by filtering an input catalog. The primary criterion is to find a locus of sources with constant half-light-radius as a function of magnitude, but the user may specify cuts on SNR, arbitrary flags, image ellipticity, or presence of zero-weight pixels. Crowded fields are handled by iteratively computing a PSF model, modeling each identified source, and removing 4-σ outliers (PSFEx Manual). It is not clear how high a source density is supported by this approach, but clearly there must be some isolated stars for the PSFEx approach to be viable.

The Stack’s version does not use the above system, instead providing an independent StarSelector that uses similar but manually specified filtering criteria, and no crowded-field handling.

PSFs are fit as a sum of basis functions, which may be pixel/nonparametric (as used by meas_extensions_psfex), Gauss-Laguerre, or a custom basis. The output PSF is always presented as a supersampled image; the developers concede that this approach does not handle PSF wings or diffraction spikes well.

The pixel-basis fitter constructs an image using “sinc interpolation” (essentially kernel smoothing by a Lanczos4 function) to resample the model to the pixel grid around each star. The fit is computed by minimizing the χ2 of the interpolated model across all sources.

PSFEx fits PSF variations simultaneously with the PSF itself. Variability may be taken with respect to source position, source parameters, or, for multi-image fits, image metadata. While PSFEx itself requires that this input be in SExtractor and FITS header format, respectively, meas_extensions_psfex does the appropriate translation from Sources and Exposures. PSF variations are represented by a low-order polynomial whose “coefficients” are themselves PSF images. There is no support for high-frequency components of the type discussed by [ChangMarshallJernigan+12]. The algorithm requires \(\mathcal{O}((N_\mathrm{PSF} N_\mathrm{coeff})^3 + N_\mathrm{PSF}^2 N_\mathrm{star} N_\mathrm{coeff})\) time, where \(N_\mathrm{coeff}\) is the number of terms in the polynomial.

PSFEx supports the use of PCA to get an optimized image basis. Presumably this is in the context of PSF variation fitting, but the documentation is very vague about what the principal component analysis produces. It can also produce PSF homogenization kernels on request.

2.2   PIFF

PIFF is a new PSF fitting tool that is still under development. As most of the current version consists of placeholders (e.g., the only PSF model supported is a single elliptical Gaussian), I will discuss this program only briefly.

PIFF is designed for wide-field telescopes and multi-chip cameras. It plans to offer fitting of PSFs as a sum of basis functions, which may be pixel/nonparametric, shapelet, or Gaussian. PSF variations will be supported (only?) with respect to position, and may be interpolated using polynomials, kriging, or PSFEnt. It will also allow the user to specify the optical component of the PSF while fitting for the atmospheric contribution.

It is not clear how PIFF will select PSF template sources, except that it will involve some kind of outlier rejection scheme.

3   PSF Fitting in Dense Fields

3.1   Iterative Derivation

[AndersonBedinPiotto+06] present a method to fit the PSF in dense fields for WFI on the ESO 2.2m. Beginning with crude source detections, centroided positions, and aperture fluxes, they reconstruct the PSF from visible stars, use the PSF to get new source positions and fluxes, use the improved source data to construct better PSFs, and so on. Template stars are selected to have high counts and no nearby neighbors; stars that are poor fits to the PSF can be rejected. The PSF is nonparametric, but with smoothness and centering constraints enforced at each iteration.

Each iteration of the algorithm requires \(\mathcal{O}(N_\mathrm{PSF} N_\mathrm{star} + N_\mathrm{PSF} N_\mathrm{smooth})\) time, where \(N_\mathrm{smooth}\) is the number of pixels in the PSF smoothing kernel. No stopping condition is presented, although the discussion suggests the algorithm converges quickly. Spatial variation is handled by solving an independent PSF for each section of a chip, then using bilinear interpolation to find the PSF at an arbitrary position (this adds an \(N_\mathrm{cells}\) factor to the complexity).

While [AndersonBedinPiotto+06] have tested their method in Baade’s Window (see their Fig. 2), they have no explicit support for cases where isolated stars do not exist. Forward modeling is used only to solve for a field’s astrometry and photometry once the final PSF model is available.

Their method can use saturated stars to fit PSF wings, but they admit the extended wings are not very accurate.

3.2   Blind Deconvolution

[SchreiberLaCameraPratoDiolaiti13] present a deconvolution-based nonparametric PSF estimator for adaptive optics observations of extremely dense, low-SNR fields. It does not require an explicit source selection. Starting from an initial PSF (which can be quite crude) and model image, they use scaled gradient projection to solve for a better image model, followed by using the refined image to improve the PSF, and so on. The algorithm is prone to overfitting, and their reconstructed PSF starts developing holes as they go to too many iterations (their Figures 3 and 4). Each iteration requires \(\mathcal{O}(N_\mathrm{PSF}^3 + N_\mathrm{pixels}^3)\) time. The authors say a “few hundred” iterations suffice to give a good PSF model, but give no explicit stopping condition.

Using simulated data, the authors show that they can get similar completeness and \(\sim 15\%\) worse photometry using the reconstructed PSF compared to using the true PSF. However, their background estimation method (involving an initial pass, creation of a source-subtracted smoothed image, and a second pass) does make their final catalog about 0.5 mag shallower than it would be with a perfectly estimated background. The effect may be less prominent for non-AO PSFs, which have smaller wings.

Unlike many modern PSF fitting algorithms (e.g., [AndersonBedinPiotto+06], [ChangMarshallJernigan+12]), the authors do not impose any kind of smoothness constraint on their PSF. It may be worth investigating whether, with such a constraint, the deconvolved PSF would be better-behaved, or whether the solution would diverge in some other way.

3.3   Forward Modeling with Image Differencing

Section 6.11.1 of [S+17] mentions a proposed algorithm by Lupton & Bosch to estimate PSFs in crowded LSST fields. Starting from an initial PSF (which can be quite crude, but should be narrower than the true PSF) and source list, a model image can be created, and [AlardLupton98] image differencing run on the original image and the model. The PSF is defined as the convolution of the previous PSF with the image matching kernel, and a new source list and model image are created. Each iteration requires \(\mathcal{O}(N_\mathrm{pixels}^2 + N_\mathrm{PSF} N_\mathrm{kernel})\) time, where \(N_\mathrm{kernel}\) is the number of pixels in the matching kernel. It is not clear how many iterations would be required or what the stopping criterion would be.

To work on the most crowded fields, this method requires a source identification algorithm that can deal with blended stars. Spatial PSF variations can be introduced using the [AlardLupton98] matching kernel.

4   PSF Fitting of Short Exposures

4.1   PSFEnt

PSFEnt [ChangMarshallJernigan+12] is a nonparametric PSF interpolation scheme to reconstruct small-scale PSF variations using maximum-entropy fitting. While it can reproduce arbitrary structure in the PSF as a function of detector position, it requires a parametric model for the PSF and (for best performance) prior knowledge of atmospheric turbulence properties.

PSFEnt models PSF variations as a grid of independent cells that are bilinearly interpolated to get the value at a specific point, much like the Stack does. The model is divided into seven “hidden” layers that are each forced to be smooth on a different spatial scale.

PSFEnt requires iterative maximization of a function whose complexity is \(\mathcal{O}(N_\mathrm{cells}^3 + N_\mathrm{layers} N_\mathrm{cells})\). It has been deemed too slow even for the Level 2 Data Reduction Pipeline (Bosch, priv. comm. 2017), so it is also too slow for Alert Production.

4.2   Compressive Sampling

[Suksmono14] proposes a method to reconstruct small-scale PSF variations from atmospheric turbulence by using properties of \(1/f\) random fields. The PSF model must be representable as a complex field; the authors, motivated by weak lensing work, use ellipticities.

The PSF variations are assumed to be sparse in the Fourier domain, and the sparsest solution consistent with the observations can be reconstructed by several optimization algorithms, which typically have \(\mathcal{O}(N_\mathrm{star} \log^2 N_\mathrm{star} + N_\mathrm{star} N_\mathrm{cells}^2 \log N_\mathrm{cells})\) complexity.

[Suksmono14] try their method on the [KitchingRoweGill+13] data set and find it can reconstruct ellipticity errors as well as other algorithms (and somewhat better than a polynomial fit to the PSF variations). This is encouraging, although ellipticity errors are not a metric relevant to the Alert Production pipeline.

While a single complex field is far too simple a model for LSST’s purposes, it may be possible to adapt their algorithm to fit an elliptical shear of a more general model. However, it is not clear how well this algorithm can handle a combination of atmospheric and optical PSF variation; while the [KitchingRoweGill+13] data does include terms for astigmatism, defocus, and coma as well as Kolmogorov turbulence, there is no theoretical basis for modeling the optics contributions using compressive sampling.

5   Future Work

The main priority for improving LSST PSF estimation at the time of writing is robust handling of crowded fields, which PSFEx handles poorly (Reiss, priv. comm. 2017). Early attempts to improve crowded field handling will likely involve implementing Lupton’s image differencing. While this algorithm is appealing because of its quadratic running time (assuming a constant bound on the number of iterations) and reuse of existing code, it is untested and may fail if it cannot find suitable PSF template stars.

The best approach may involve combining the strengths of recent algorithms. For example, blind deconvolution is appealing because it does not require source identification, but the published algorithm requires hand tuning to avoid overfitting. A regularization scheme like that used by [AndersonBedinPiotto+06] may make it more stable without significantly increasing its asymptotic complexity.

Another possibility is to use multiple algorithms in different contexts. For example, we may find that Lupton’s algorithm performs well on all but the most crowded fields, where source identification fails catastrophically. If so, its high speed would make it the preferred algorithm for LSST data reduction, and we could fall back to a much slower algorithm like blind deconvolution in dense fields (blind deconvolution has nearly as extreme computational requirements as PSFEnt, so while it’s promising from a reliability standpoint it is likely to be too expensive to run on all LSST images). So long as fewer than 2% of all images require an expensive algorithm, such a strategy can comply with LSST system requirements ([C+16], LSR-REQ-0025).

6   References

[C+16]Charles Claver and others. LSST System Requirements. LPM-29,, August 2016.
[I+11]Željko Ivezić and others. LSST Science Requirements Document. LPM-17,, July 2011.
[S+17](1, 2) John Swinbank and others. LSST Data Management Science Pipelines Design. LDM-151,, May 2017.
[AlardLupton98](1, 2, 3) C. Alard and R. H. Lupton. A Method for Optimal Image Subtraction. ApJ, 503:325–331, August 1998. arXiv:astro-ph/9712287, doi:10.1086/305984.
[AndersonBedinPiotto+06](1, 2, 3, 4, 5) J. Anderson, L. R. Bedin, G. Piotto, R. S. Yadav, and A. Bellini. Ground-based CCD astrometry with wide field imagers. I. Observations just a few years apart allow decontamination of field objects from members in two globular clusters. A&A, 454:1029–1045, August 2006. arXiv:astro-ph/0604541, doi:10.1051/0004-6361:20065004.
[ChangMarshallJernigan+12](1, 2, 3, 4) C. Chang, P. J. Marshall, J. G. Jernigan, J. R. Peterson, S. M. Kahn, S. F. Gull, Y. AlSayyad, Z. Ahmad, J. Bankert, D. Bard, A. Connolly, R. R. Gibson, K. Gilmore, E. Grace, M. Hannel, M. A. Hodge, L. Jones, S. Krughoff, S. Lorenz, S. Marshall, A. Meert, S. Nagarajan, E. Peng, A. P. Rasmussen, M. Shmakova, N. Sylvestre, N. Todd, and M. Young. Atmospheric point spread function interpolation for weak lensing in short exposure imaging data. MNRAS, 427:2572–2587, December 2012. arXiv:1206.1383, doi:10.1111/j.1365-2966.2012.22134.x.
[KitchingRoweGill+13](1, 2) T. D. Kitching, B. Rowe, M. Gill, C. Heymans, R. Massey, D. Witherick, F. Courbin, K. Georgatzis, M. Gentile, D. Gruen, M. Kilbinger, G. L. Li, A. P. Mariglis, G. Meylan, A. Storkey, and B. Xin. Image Analysis for Cosmology: Results from the GREAT10 Star Challenge. ApJS, 205:12, April 2013. arXiv:1210.1979, doi:10.1088/0067-0049/205/2/12.
[SchreiberLaCameraPratoDiolaiti13]L. Schreiber, A. La Camera, M. Prato, and E. Diolaiti. Point Spread Function extraction in crowded fields using blind deconvolution. In S. Esposito and L. Fini, editors, Proceedings of the Third AO4ELT Conference, 78. December 2013. doi:10.12839/AO4ELT3.13358.
[Suksmono14](1, 2) A. B. Suksmono. Interpolation of PSF based on compressive sampling and its application in weak lensing survey. MNRAS, 443:919–926, September 2014. doi:10.1093/mnras/stu1169.