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 χ² 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.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.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.