EPSFBuilder#

class photutils.psf.EPSFBuilder(*, oversampling=4, shape=None, smoothing_kernel='auto', alias_passband='auto', sigma_clip=<default: SigmaClip(sigma=3.0, maxiters=10)>, recentering_func=<function centroid_symmetry>, recentering_boxsize=(5, 5), recentering_maxiters=20, center_accuracy=0.001, converged_fraction=0.95, fitter=None, fit_shape='auto', fitter_maxiters=100, constrain_fluxes=True, maxiters=10, refinement_iters=5, wing_smoothing=True, progress_bar=True)[source]#

Bases: object

Class to build an effective PSF (ePSF).

The method is based on Anderson and King 2000 (PASP 112, 1360) and Anderson 2016 (WFC3 ISR 2016-12). It differs from them in several steps (see Notes).

Parameters:
oversamplingint or array_like (int), optional

The integer oversampling factor(s) of the output ePSF relative to the input stars along each axis. If oversampling is a scalar then it will be used for both axes. If oversampling has two elements, they must be in (y, x) order. The ePSF should have at least about four grid points per FWHM, so the value should be at least 4 / FWHM with the FWHM in pixels. The default of 4 is a good choice for a FWHM of about 1 pixel or more. A value above 4 / FWHM increases the run time and the memory use, and in tests a value of 8 was no more accurate than 4. For well-sampled data (a FWHM of 4 pixels or more), a value of 1 gives the same fitted star positions and fluxes in less time. See the guidelines in the ePSF building user guide for details.

shapeint, tuple of two ints, or None, optional

The (ny, nx) shape of the output ePSF. If the input shape is even along any axis, it will be made odd by adding one (with a warning). If the shape is None, it will be derived from the sizes of the input stars and the ePSF oversampling factor. The output ePSF will always have odd sizes along both axes to ensure a well-defined central pixel.

smoothing_kernelstr, 2D ndarray, or None, optional

The smoothing kernel to apply to the ePSF during each iteration step. A string must be 'auto', 'quartic', or 'quadratic'. If 'auto', a least-squares quartic polynomial kernel (see Notes) whose width is the largest odd number of oversampled grid points that is not larger than 0.7 times the FWHM of the current ePSF is used, and no smoothing is applied if that width is less than 5 grid points, i.e., for ePSFs with fewer than about 7 grid points per FWHM. The kernel is square and the FWHM is measured in each iteration along the narrowest axis of the ePSF, so with anisotropic oversampling the axis with the fewer grid points per FWHM sets the kernel size. If the FWHM cannot be measured, the 'quartic' kernel is used and a warning is emitted. The predefined 'quartic' and 'quadratic' kernels are 5x5 kernels (in oversampled grid points) derived from fourth and second degree polynomials, respectively. Alternatively, a custom 2D array can be input. If None then no smoothing will be performed. The kernels are applied on the oversampled grid, so the physical width of a fixed kernel depends on the oversampling factor. Structure that repeats with a period of about one input pixel or shorter is removed from the ePSF along oversampled axes independently of this parameter (see alias_passband).

alias_passband{‘auto’}, float, or None, optional

The end of the passband, in cycles per detector pixel, of the low-pass filter that is applied to the ePSF in each iteration along the axes with an oversampling factor greater than one. A spatial frequency of one cycle per pixel describes structure that repeats with a period of one detector pixel, and a frequency of 0.5 cycles per pixel describes structure that repeats every two pixels. The filter has unit gain up to alias_passband, a raised-cosine transition, and zero gain at and above one cycle per pixel. It removes the frequencies at which the star-pixel sampling lattice aliases onto the oversampled grid, which otherwise can grow into a checkerboard pattern (see Notes). The value must be greater than 0 and less than 1.

If 'auto' (default), the passband ends at 0.8 cycles per pixel, or at 0.7 cycles per pixel for an oversampling factor of 2 (70% of the Nyquist frequency of the oversampled grid).

The default is the best choice for most data. A different value can help in two cases, which are distinguished by the optical cutoff frequency of the telescope in cycles per pixel, cutoff = D * pixel_scale / wavelength, with the telescope diameter D and the mean wavelength of the bandpass in the same units and the pixel_scale in radians per pixel:

  • cutoff greater than about 1 (strongly undersampled, e.g., HST WFC3/IR F110W, JWST NIRCam F070W, or Roman WFI F062 and F106): the ePSF has real signal up to nearly one cycle per pixel. In tests the default recovered the peak of such ePSFs to within about 1 percent, except for the sharpest ones, whose peak was 3 percent low. A value of 0.85 with a maxiters of 15 to 20 recovered the peaks to within 0.3 percent and lowered the residuals in the core by 10 to 35 percent, except for the sharpest ePSF, which needed 0.9. Try 0.9 only if the ePSF built with 0.85 is still too broad, i.e., if the stars have positive residuals at their centers after the fitted ePSF is subtracted. A value above the default needs a large star sample (a few hundred stars), an oversampling of 4 or more, and more iterations. If the stars were observed in several dithered exposures, link them and keep the default instead (see constrain_fluxes). In tests, linking the stars of four exposures recovered the peaks to within 0.4 percent.

  • cutoff less than about 0.9 (e.g., JWST NIRCam F115W and redder, JWST MIRI, or most ground-based data): the ePSF has no signal to preserve near one cycle per pixel. A value of 0.7 rejects more noise (10 to 40 percent lower residuals in the core in tests) and converges in fewer iterations. It leaves the peak of a strongly undersampled ePSF low by 2 to 6 percent, so use it only when cutoff is known.

A value larger than needed makes the build converge more slowly, because a star sampled once per pixel constrains the frequencies near one cycle per pixel only weakly, and it is not recommended for small star samples or for an oversampling factor of 2.

If None, the filter is not applied. This is rarely appropriate. Without the filter, noise at the alias frequencies accumulates over the iterations, the build can stall before it converges, and heterogeneous or contaminated star samples can grow a checkerboard pattern. In tests with simulated HST, JWST, and Roman star fields, the unfiltered ePSF was less accurate than the filtered one in nearly every case, even for large, clean, and homogeneous star samples. The option is provided for experimentation, e.g., to check how much the filter changes a particular ePSF. Always compare the result with a filtered build. The filter is never applied along an axis with an oversampling factor of 1.

sigma_clipastropy.stats.SigmaClip instance, optional

A SigmaClip object that defines the sigma clipping parameters used to determine which pixels are ignored when stacking the ePSF residuals in each iteration step. The default is a SigmaClip with sigma=3.0 and maxiters=10. If None then no sigma clipping will be performed.

recentering_funccallable, optional

A callable object that is used to calculate the centroid of a 2D array. The callable must accept a 2D ndarray and have a mask keyword. The callable object must return a tuple of (x, y) centroids. The default is centroid_symmetry, the point about which the core of the ePSF is most symmetric, which is the center definition of Anderson 2016. With the default recentering_boxsize it measures the symmetry within about 1.5 detector pixels of the center. It needs a box of at least 5 oversampled pixels and an ePSF that is not constant. If it cannot be calculated (e.g., with a recentering_boxsize of 3 and an oversampling of 1), the center of mass (centroid_com) is used instead, and a warning is emitted if the box is too small. The center of mass is also used for a recentering step in which the ePSF is far from the center of the box (e.g., if the input star positions share a large offset), where the symmetry center cannot be found. This fallback applies only to centroid_symmetry itself and not to a function that wraps it. The center of mass of an asymmetric ePSF is pulled toward the asymmetric structure around its core, so the two definitions give centers that differ by a constant offset (up to about 0.1 pixel for undersampled space-telescope PSFs).

recentering_boxsizeint or tuple of two ints, optional

The size (in pixels) of the box used to calculate the centroid of the ePSF during each build iteration. The size is in the input star (i.e., undersampled) pixel space. It is automatically scaled by the oversampling factor when applied to the oversampled ePSF grid. If a single integer number is provided, then a square box will be used. If two values are provided, then they must be in (ny, nx) order. recentering_boxsize must have odd values and be greater than or equal to 3 for both axes.

recentering_maxitersint, optional

The maximum number of recentering iterations to perform during each ePSF build iteration.

center_accuracyfloat, optional

The desired accuracy for the centers of stars. The building iterations will stop when the centers of at least converged_fraction of the successfully fitted stars change by less than center_accuracy pixels between iterations (see converged_fraction).

converged_fractionfloat, optional

The fraction of the successfully fitted stars whose centers must change by less than center_accuracy pixels between iterations for the build to be considered converged. The change of each center is measured relative to the median change of all the centers, which must also be less than center_accuracy. A change that is common to all the stars is a shift of the center of the ePSF. It changes neither the shape of the ePSF nor the relative positions of the stars. The default of 0.95 allows a small number of stars (e.g., spurious detections or contaminated cutouts) whose centers never settle to not prevent convergence. Set to 1.0 to require all stars to converge. The fraction achieved in the final iteration is reported in the final_converged_fraction attribute of the returned EPSFBuildResults.

fitterFitter or EPSFFitter, optional

A Fitter object used to fit the ePSF to stars. If None, then the default TRFLSQFitter will be used.

Deprecated since version 3.0: Passing an EPSFFitter instance is deprecated. Use the fitter, fit_shape, and fitter_maxiters parameters instead.

fit_shape{‘auto’}, int, tuple of int, or None, optional

The size (in detector pixels) of the box centered on the star to be used for ePSF fitting. This allows using only a small number of central pixels of the star (i.e., where the star is brightest) for fitting. If 'auto', a square box of twice the FWHM of the current ePSF (measured in each iteration along its narrowest axis, in detector pixels), with a minimum of 5 pixels and a maximum of the smallest star cutout size, is used. For an elongated ePSF, the narrower axis sets the box size along both axes. If the FWHM cannot be measured, a 5x5 box is used and a warning is emitted. If fit_shape is a scalar then a square box of size fit_shape will be used. If fit_shape has two elements, they must be in (ny, nx) order. fit_shape must have odd values and be greater than or equal to 3 for both axes. If None, the fitter will use the entire star image.

fitter_maxitersint, optional

The maximum number of iterations in which the fitter is called for each star. The value can be increased if the fit is not converging. This parameter is passed to the fitter if it supports the maxiter parameter and ignored otherwise.

constrain_fluxesbool, optional

Whether to constrain the fluxes of the stars within each LinkedEPSFStar (i.e., the same physical star observed in multiple dithered images) to their mean value after each fitting iteration, in addition to constraining their centers to the same sky coordinate. This breaks the degeneracy between the flux of a star and its subpixel position caused by intra-pixel sensitivity variations, which would otherwise be absorbed into the ePSF. It assumes that the linked images have the same flux scale (e.g., the same exposure time and throughput). Set to False if the linked images have different flux scales. This parameter has no effect on stars that are not linked.

Linking the stars of dithered exposures is the recommended way to build the ePSF of an undersampled detector. In tests with simulated star fields, the peaks of strongly undersampled ePSFs built from four linked exposures were recovered to within 0.4 percent with this constraint and to within about 1 percent without it, compared to as much as 3 percent low for one exposure.

maxitersint, optional

The maximum number of ePSF building iterations to perform.

refinement_itersint, optional

The number of refinement iterations to perform after the building iterations. The alias low-pass filter of the building iterations (see alias_passband) removes some real signal of an undersampled ePSF just below one cycle per pixel. The filter acts separately along each axis, so the missing signal shows as a ripple pattern with a period of about one pixel along the row and the column through the center of the ePSF. The refinement iterations restore that signal. In each refinement iteration, the ePSF is updated five times from the star residuals with the star centers and fluxes held fixed, and the stars are then refit with the updated ePSF. The ePSF is recentered in each update, as in the building iterations, which keeps the ePSF and the star centers from drifting together. The low-pass filter of these updates has unit gain up to 1.1 cycles per pixel and zero gain at and above 1.33 cycles per pixel, so it does not remove signal near one cycle per pixel. It removes only the frequencies that the star residuals do not constrain. Such a wide filter cannot be used from the start of the build, because the build then converges slowly and is more sensitive to the initial star centers. The refinement is not performed if refinement_iters is 0, if alias_passband is None, or if the oversampling factor is less than 4 along both axes. It roughly doubles the run time of a build.

The converged attribute of the results describes the last refinement iteration. The first refit of the refinement moves some stars by more than center_accuracy, so with refinement_iters=1 a build whose building iterations converged can report converged=False. The star centers settle again within a few refinement iterations. With different oversampling factors along the two axes they can need more refinement iterations than the default (about 10 in tests).

wing_smoothingbool, optional

Whether to smooth the wings of the final ePSF more strongly than its core. Far from the center the ePSF is faint and varies slowly, so its noise can be averaged over a larger area than in the core. If True, each value of the final ePSF beyond 3.5 FWHM from its center is blended into the value of a least-squares quadratic fit to the values in a box 1.25 FWHM wide around it, and beyond 5 FWHM into the fit in a box 1.75 FWHM wide. The ePSF within 3.5 FWHM of its center is not changed, except by the renormalization of the smoothed ePSF (less than 0.03 percent in tests). The smoothing is applied once, after the last iteration, so it does not affect the star fits or the convergence of the build. The fluxes of the returned stars were therefore fit before that renormalization. It is modeled on the approach of Anderson 2016, which smooths the wings of HST ePSFs more strongly than their cores. In tests with a few hundred stars it lowered the residuals of the wings of undersampled ePSFs by up to about 50 percent. This matters when the wings are used, e.g., to subtract bright stars, to make model images, or to measure encircled energies.

The smoothing also removes real structure in the wings that is finer than about two FWHM, such as diffraction rings and spikes. Applied to noise-free ePSFs, it changed the wings by 2 to 8 percent of their mean value for JWST and Roman ePSFs and by 11 to 14 percent for HST WFC3/IR ePSFs. With a few hundred stars the noise that it removes is larger than this in every case that was tested. For a large star sample of high signal-to-noise (thousands of stars), the noise in the wings can be smaller than this change, and the wings are then more accurate without the smoothing.

The boxes are at least 5 oversampled grid points wide. For an ePSF with a FWHM of less than 4 grid points they are therefore wider than given above, and they remove more of the real structure. This matters most for an oversampling factor of 1, where the boxes of an undersampled ePSF (a FWHM of about 1.3 pixels) are nearly 4 FWHM wide. In tests with an oversampling factor of 1, the smoothing made the wings of the most undersampled HST and JWST ePSFs less accurate, by up to a factor of about 2, and those of most other ePSFs slightly more accurate. With an oversampling factor of 2 or larger it made the wings more accurate or left them unchanged in every case. For an oversampling factor of 1, compare the ePSFs built with and without the smoothing.

Set to False to keep the wings as built.

progress_barbool, optional

Whether to print the progress bar during the build iterations. The progress bar requires that the tqdm optional dependency be installed.

Notes

In each build iteration, the residual between each star and the current ePSF model is deposited on the oversampled grid points within 0.375 detector pixel (and at least one grid spacing) of each star pixel center along each axis. Each grid point is therefore estimated from the star pixels in a box three quarters of a pixel wide around it, and every star contributes to both parities of the grid regardless of its subpixel phase. Anderson and King (2000) combined the pixels within 0.25 pixel of each grid point of an ePSF with an oversampling factor of 4. That narrower box gives a slightly sharper ePSF for clean star samples, but it grows a checkerboard pattern for heterogeneous or contaminated ones. After the residuals are combined and the ePSF is smoothed, power is removed along each oversampled axis with a low-pass filter that by default has unit gain up to 0.8 cycles per input pixel (0.7 for an oversampling factor of 2) and zero gain at and above one cycle per input pixel (see alias_passband). A pixel-integrated PSF has essentially no power at one cycle per input pixel, but that is the frequency at which the star-pixel sampling lattice aliases onto the oversampled grid. Without these two measures, noise in the ePSF grid from heterogeneous or contaminated stars can bias the fitted star centers toward particular subpixel phases and grow into a checkerboard pattern in the ePSF. A warning is emitted if the subpixel phases of the fitted star centers are strongly nonuniform at the end of the build.

The default smoothing_kernel='auto' and fit_shape='auto' scale the smoothing kernel and the fitting box with the FWHM of the ePSF. The 5x5 'quartic' kernel (in oversampled grid points) of Anderson and King 2000 and the 5-pixel fitting box (in detector pixels) of Anderson 2016 were designed for HST images with an oversampling factor of 4, where they correspond to about 0.7 and 2.5 FWHM. A fixed kernel oversmooths heavily undersampled ePSFs, and a fixed 5-pixel fitting box applied to a well-sampled star uses only its flat core, which biases the fitted centers and can prevent convergence. The polynomial kernels replace each grid value by the value at the center of a least-squares polynomial fit to the surrounding grid values. The kernel and fitting box chosen in the final iteration are reported in the smoothing_kernel and fit_shape attributes of the returned EPSFBuildResults.

This class follows the ePSF building procedure of Anderson and King 2000 and Anderson 2016, with several modifications. The main differences are listed below. See the ePSF building user guide (Differences from the Anderson and King Algorithm) for the complete comparison and the reasons for each difference.

  • A single ePSF is built with any integer oversampling factor. Anderson builds a 3x3 array of ePSFs across each detector with an oversampling factor of 4.

  • The star cutouts must be background subtracted by the user. Anderson measures the background of each star in an annulus around it.

  • The residuals are combined in a box of 0.375 pixel half width with a sigma-clipped median (3 sigma by default). Anderson uses a half width of 0.25 pixel and a mean with rejection at 2.5 sigma.

  • The smoothing kernel and the wing smoothing scale with the FWHM of the ePSF, and the wings are smoothed only once, after the last iteration. Anderson uses a 5x5 quartic kernel in the core and stronger smoothing beyond fixed radii of 3 to 5 pixels, in every iteration.

  • A low-pass Fourier filter is applied in every iteration for oversampling factors greater than 1 (see alias_passband). Anderson applies none.

  • The ePSF is centered on its point of maximal symmetry within a radius of about 1.5 pixels by default, as in Anderson (2016). Anderson (2000) requires equal values half a pixel on either side of the center.

  • The ePSF is normalized so that its values sum to the product of the oversampling factors over the whole grid. Fitted fluxes are therefore the fluxes within the area of the grid. Anderson normalizes the ePSF to unit flux in the central 5x5 pixels (2000) or within a radius of 5.5 pixels (2016).

  • The stars are fit after every ePSF update until their centers converge, and the ePSF is then refined with five updates per fit (see refinement_iters). Anderson uses five updates per fit throughout and iterates until the fitted positions and fluxes show no trend with pixel phase.

  • The stars are fit in a box of twice the FWHM with a nonlinear least-squares fitter. Anderson fits the pixels within about 2 pixels of the center (2000) or the central 5x5 pixels (2016) with Poisson weights.

  • Dithered exposures are optional (see LinkedEPSFStar). They are central to Anderson’s method, which averages the positions and fluxes of each star over the exposures after every fit.

  • The ePSF is evaluated with a single bicubic spline over the whole grid. Anderson uses a bicubic spline within 4 pixels of the center and bilinear interpolation farther out (2016).

This class stores per-call state on the instance (e.g., the automatic smoothing kernel), so a single instance must not be called concurrently from multiple threads. Create one instance per thread for concurrent use. Sharing a single Astropy fitter instance across concurrently-used objects is also unsafe because Astropy fitters store fit_info on themselves.

Attributes Summary

coord_transformer

Methods Summary

__call__(stars)

Build an ePSF from input stars.

build_epsf(stars, *[, epsf])

Build iteratively an ePSF from star cutouts.

Attributes Documentation

coord_transformer#

Deprecated since version 3.1: The coord_transformer attribute is deprecated and may be removed in a future version.

Methods Documentation

__call__(stars)[source]#

Build an ePSF from input stars.

Parameters:
starsEPSFStars

The stars used to build the ePSF.

Returns:
resultEPSFBuildResults

The result of the ePSF building process.

build_epsf(stars, *, epsf=None)[source]#

Build iteratively an ePSF from star cutouts.

This method builds an ePSF from an initial model when epsf is provided, or from scratch when epsf is None. In the latter case, it is equivalent to invoking an EPSFBuilder instance on the input stars.

Parameters:
starsEPSFStars object

The stars used to build the ePSF.

epsfImagePSF object, optional

The initial ePSF model. Its oversampling must match the oversampling of the builder, and its shape must be at least the largest star cutout shape times the oversampling (plus one along an axis where that product is even). The output ePSF has the shape of this model. If None, then the ePSF will be built from scratch.

Returns:
resultEPSFBuildResults

The ePSF building results, with detailed information about the building process. For backward compatibility, the result can be unpacked as a tuple: (epsf, fitted_stars) = epsf_builder(stars).

Raises:
ValueError

If the oversampling or shape of the input epsf is incompatible with the builder or the stars, if the shape of the builder is too small for the stars, if the stars are invalid, or if the ePSF fit fails for all stars.

Notes

The structured result object contains:

  • epsf: The final constructed ePSF

  • fitted_stars: Stars with updated centers/fluxes

  • iterations: Number of iterations performed

  • converged: Whether convergence was achieved

  • final_center_accuracy: Final center movement accuracy

  • final_converged_fraction: Fraction of stars that met the accuracy

  • n_excluded_stars: Number of stars excluded due to fit failures

  • excluded_star_indices: Indices of excluded stars

  • smoothing_kernel: Smoothing kernel of the final iteration

  • fit_shape: Fitting box of the final iteration

  • iteration_epsfs: ePSF image after each iteration

  • iteration_info: Table of the convergence statistics of each iteration

  • initial_epsf: Image of the input ePSF, if any

  • center_asymmetry: Offset of the center of mass of the final ePSF from its symmetry center