Point Spread Function (PSF) models

This page covers PSF model construction and PSF fitting photometry. For object detection and aperture/optimal photometry see Object detection and measurement, and for photometric calibration see Photometric calibration.

STDPipe includes basic support for point spread function (PSF) construction and analysis through the interface to PSFEx code that is able to build supersampled PSF models from the object lists created using SExtractor. Please consider checking its documentation to better understand the concepts of PSFEx operation and its possible configuration options.

The stdpipe.psf.run_psfex() routine transparently calls the SExtractor on the image and then uses its output (it needs a special set of measured object parameters and postage stamps for them) to run PSFEx. The results are then parsed and returned as a structure representing the PSF model file as described here. This structure may later be used to visualize PSF shape at arbitrary position inside the image, as well as to generate artificial stars and inject them into the image.

stdpipe.psf.run_psfex(image, mask=None, thresh=2.0, aper=None, r0=0.0, gain=1, minarea=5, vignet_size=None, psf_size=None, order=0, sex_extra=None, checkimages=None, extra=None, psffile=None, get_obj=False, _workdir=None, _tmpdir=None, _exe=None, _sex_exe=None, verbose=False)

Wrapper around PSFEx to help extracting PSF models from images.

For the details of PSFEx operation we suggest to consult its documentation at https://psfex.readthedocs.io

Parameters:

imagenumpy.ndarray

Input image as a NumPy array.

masknumpy.ndarray, optional

Image mask as a boolean array (True values will be masked).

threshfloat, optional

Detection threshold in sigmas above local background, for running initial SExtractor object detection.

aperfloat, optional

Circular aperture radius in pixels, to be used for PSF normalization. Should contain most of object flux. If not specified, will be estimated as twice the FWHM.

r0float, optional

Smoothing kernel size (sigma) to be used for improving object detection in initial SExtractor call.

gainfloat, optional

Image gain.

minareaint, optional

Minimal number of pixels in the object to be considered a detection (DETECT_MINAREA parameter of SExtractor).

vignet_sizeint, optional

The size of postage stamps to be used for PSF model creation.

psf_sizeint, optional

The size of the supersampled PSF model.

orderint, optional

Spatial order of PSF model variance.

sex_extradict, optional

Dictionary of additional options to be passed to SExtractor for initial object detection (extra parameter of stdpipe.photometry.get_objects_sextractor()).

checkimageslist, optional

List of PSFEx checkimages to return along with PSF model.

extradict, optional

Dictionary of extra configuration parameters to be passed to PSFEx call, with keys as parameter names. See psfex -dd for the full list.

psffilestr, optional

If specified, PSF model file will also be stored under this file name, so that it may e.g. be re-used by SExtractor later.

get_objbool, optional

If set, also return the table with SExtractor detected objects.

_workdirstr, optional

If specified, all temporary files will be created in this directory, and will be kept intact after running SExtractor and PSFEx. May be used for debugging exact inputs and outputs of the executable.

_tmpdirstr, optional

If specified, all temporary files will be created in a dedicated directory (that will be deleted after running the executable) inside this path.

_exestr, optional

Full path to PSFEx executable. If not provided, the code tries to locate it automatically in your PATH.

_sex_exestr, optional

Full path to SExtractor executable. If not provided, the code tries to locate it automatically in your PATH.

verbosebool or callable, optional

Whether to show verbose messages during the run of the function or not.

Returns:

dict

PSF structure corresponding to the built PSFEx model.

The structure has at least the following fields:

• width, height - dimensions of supersampled PSF stamp

• fwhm - mean full width at half maximum (FWHM) of the images used for building the PSF model

• sampling - conversion factor between PSF stamp (supersampled) pixel size, and original image one (less than unity when supersampled resolution is finer than original image one)

• ncoeffs - number of coefficients pixel polynomials have

• degree - polynomial degree of a spatial variance of PSF model

• data - the data containing per-pixel polynomial coefficients for PSF model

• header - original FITS header of PSF model file, if get_header=True parameter was set

This structure corresponds to the contents of original PSFEx generated output file that is documented at https://psfex.readthedocs.io/en/latest/Appendices.html#psf-file-format-description

Attention

While PSFEx may build PSF models depending on a variety of object parameters as described here, STDPipe currently allows reading and interpreting only the ones depending on the position on the image (e.g. representing optical distortions). On the other hand, the original PSFEx model file may be stored to the file using psffile parameter of stdpipe.psf.run_psfex(), and then supplied directly to SExtractor through psf parameter of stdpipe.photometry.get_objects_sextractor() to be used for SExtractor PSF photometry.

Using PSF models

The PSF model built by PSFEx may also be loaded from PSFEx model file. This model may then be used to construct PSF stamps both in supersampled PSF model pixels, or original image pixels.

stdpipe.psf.load_psf(filename, get_header=False, verbose=False)

Load PSF model from PSFEx file

The structure may be useful for inspection of PSF model with stdpipe.psf.get_supersampled_psf_stamp() and stdpipe.psf.get_psf_stamp(), as well as for injection of PSF instances (fake objects) into the image with stdpipe.psf.place_psf_stamp().

Parameters:

filenamestr

Name of a file containing PSF model built by PSFEx.

get_headerbool, optional

Whether to return the original FITS header of PSF model file or not. If set, the header will be stored in the header field of the returned structure.

verbosebool or callable, optional

Whether to show verbose messages during the run of the function or not.

Returns:

dict

PSF structure in the same format as returned from stdpipe.psf.run_psfex().

stdpipe.psf.get_supersampled_psf_stamp(psf, x=0, y=0, normalize=True)

Returns supersampled PSF model for a given position inside the image.

The returned stamp corresponds to PSF model evaluated at a given position inside the image, with its center always in the center of central stamp pixel. Every supersampled pixel of the stamp corresponds to psf['sampling'] pixels of the original image.

Parameters:

psfdict

Input PSF structure as returned by stdpipe.psf.run_psfex() or stdpipe.psf.load_psf().

xfloat, optional

x coordinate of the position inside the original image to evaluate the PSF model.

yfloat, optional

y coordinate of the position inside the original image to evaluate the PSF model.

normalizebool, optional

Whether to normalize the stamp to have flux exactly equal to unity or not.

Returns:

numpy.ndarray

Stamp of the PSF model evaluated at the given position inside the image.

stdpipe.psf.get_psf_stamp(psf, x=0, y=0, dx=None, dy=None, normalize=True)

Returns PSF stamp in original image pixel space with sub-pixel shift applied.

The PSF model is evaluated at the requested position inside the original image, and then downscaled from supersampled pixels of the PSF model to original image pixels, and then adjusted to accommodate for requested (dx, dy) sub-pixel shift.

Stamp is odd-sized, with PSF center at:

x0 = floor(width/2) + dx
y0 = floor(height/2) + dy

If dx=None or dy=None, they are computed directly from the floating point parts of the position x and y arguments:

dx = x - round(x)
dy = y - round(y)

The stamp should directly represent stellar shape at a given position (including sub-pixel center shift) inside the image.

Parameters:

psfdict

Input PSF structure as returned by stdpipe.psf.run_psfex() or stdpipe.psf.load_psf().

xfloat, optional

x coordinate of the position inside the original image to evaluate the PSF model.

yfloat, optional

y coordinate of the position inside the original image to evaluate the PSF model.

dxfloat, optional

Sub-pixel adjustment of PSF position in image space, x direction.

dyfloat, optional

Sub-pixel adjustment of PSF position in image space, y direction.

normalizebool, optional

Whether to normalize the stamp to have flux exactly equal to unity or not.

Returns:

numpy.ndarray

Stamp of the PSF model evaluated at the given position inside the image, in original image pixels.

PSF model utilities

Two small helpers complement the stamp-construction routines above:

• enclosed_psf_fraction() – fraction of normalised PSF flux inside a circular aperture at an arbitrary image position. Vectorised over radius (scalar in → scalar out, array in → array out). Useful for on-the-fly aperture corrections.

• select_psf_seeds() – pick bright, edge-clear sources distributed uniformly across the field for use as PSF seeds. Bins on a regular grid_x × grid_y spatial grid and keeps the max_per_tile brightest per cell. Configurable column names so it works directly on SExtractor catalogues (X_IMAGE / FLUX_AUTO).

from stdpipe import psf

# Aperture correction at one source's position, as a function of radius
radii = [3.0, 5.0, 7.0, 10.0]
frac = psf.enclosed_psf_fraction(psf_model, x=512.0, y=512.0, radius=radii)
print('aperture corrections:', dict(zip(radii, frac)))

# Build a uniform PSF-seed sample for create_psf_model() / run_psfex()
seeds = psf.select_psf_seeds(
    obj, image_shape=image.shape,
    max_per_tile=20, grid=(8, 6), edge=20,
    accept_flags=0,
)

stdpipe.psf.enclosed_psf_fraction(psf, x=0, y=0, radius=None, subpixel=10)

Fraction of normalised PSF flux inside a circular aperture.

Evaluates the (possibly position-dependent) PSF model at (x, y) via get_psf_stamp(), normalises it to unit total flux, and sums the pixels lying within radius of the stamp centre. Useful for on-the-fly aperture corrections.

Pixels are weighted by the fractional coverage of the circular aperture using a sub-pixel grid (subpixel × subpixel samples per pixel) — set subpixel=1 to fall back to a hard pixel-centre mask, which is faster but has up to ~10 % radius-dependent jitter when the aperture edge cuts through individual pixels.

Parameters:

psfdict

PSF model from run_psfex(), load_psf(), or create_psf_model().

x, yfloat, optional

Position in image pixels at which to evaluate a position-dependent PSF model. The sub-pixel parts of x and y are folded into the stamp centre via get_psf_stamp().

radiusfloat or array-like

Aperture radius (or radii) in image pixels.

subpixelint, optional

Linear sub-pixel sampling per image pixel for partial-coverage weighting. Default 10 (100 sub-samples per pixel).

Returns:

float or ndarray

Enclosed flux fraction at each radius. Scalar when radius is a scalar, otherwise an array of the same shape.

stdpipe.psf.select_psf_seeds(obj, image_shape, *, max_per_tile=25, grid=6, edge=20, obj_col_x='x', obj_col_y='y', obj_col_flux='flux', obj_col_flags='flags', accept_flags=0)

Select bright, edge-clear sources spread uniformly across the field.

Bins obj on a regular grid_x × grid_y spatial grid and returns the max_per_tile brightest sources per cell that pass basic quality filters: finite x/y/flux, positive flux, optionally a flag mask, and inside the image after an edge-pixel margin.

The returned table is a strict subset of obj (same columns, sub-selected rows) and is the natural input for stdpipe.psf.create_psf_model() and similar PSF-modelling routines that benefit from a uniform spatial sampling.

Parameters:

objastropy.table.Table

Source catalogue.

image_shape(H, W) tuple

Image shape used for the spatial grid and edge mask.

max_per_tileint

Maximum number of seeds kept per spatial cell.

gridint or (nx, ny) tuple

Number of cells in the spatial grid. A scalar means a square grid × grid layout.

edgefloat

Edge margin in pixels; sources within edge of any image boundary are excluded.

obj_col_x, obj_col_ystr

Column names for source positions in pixel coordinates.

obj_col_fluxstr

Column name used to rank candidates within each cell (brightest first). Sources with non-positive values are dropped.

obj_col_flagsstr, optional

Column name for a SExtractor-style integer flag mask. Set to None to skip flag filtering. If the column is missing from obj, no flag filter is applied either.

accept_flagsint

Bitwise mask of acceptable flags. A source is kept only when flags & ~accept_flags == 0.

Returns:

astropy.table.Table

A subset of obj containing the selected seeds.

Placing artificial stars into the image

STDPipe also contains a couple of convenience functions to directly place a PSF model at a given image position - that would correspond to artificial point source injection into the image. Higher-level one, stdpipe.pipeline.place_random_stars(), inserts a number of randomly generated sources with realistically (log-uniform in fluxes) distributed brightness at random positions, as well as returns the catalogue of injected stars that may be used e.g. in deriving the object detection efficiency.

# Derive PSF model assuming that it does not change over the image
psf_model = psf.run_psfex(image, mask=mask, order=0, verbose=True)

# ..and now perform the detection efficiency analysis

sims = [] # will hold simulated stars

# We will repeatedly inject the stars (20 at a time), detect objects, and
# match them in order to see whether simulated stars are detectable

for _ in tqdm(range(100)):
    # We will operate on a copy of original image
    image1 = image.copy()

    # We will use this value as a saturation level
    saturation = 50000

    # Simulate 20 random stars
    sim = pipeline.place_random_stars(image1, psf_model, nstars=20, minflux=10, maxflux=1e6,
             wcs=wcs, gain=gain, saturation=saturation)

    # Initial metadata for injected stars
    sim['mag_calib'] = sim['mag'] + zero_fn(sim['x'], sim['y']) # Apply zero point from photometric calibration
    sim['detected'] = False
    sim['mag_measured'] = np.nan
    sim['magerr_measured'] = np.nan
    sim['flags_measured'] = np.nan

    # Mask corresponding to the saturation level
    mask1 = image1 >= saturation

    obj1 = photometry.get_objects_sextractor(image1, mask=mask|mask1, r0=1, aper=5.0,
             wcs=wcs, gain=gain, minarea=3, sn=5)
    # Apply zero point from photometric calibration
    obj1['mag_calib'] = obj1['mag'] + zero_fn(obj1['x'], obj1['y'])

    # Positional match within FWHM/2 radius
    oidx,sidx,dist = astrometry.spherical_match(obj1['ra'], obj1['dec'], sim['ra'], sim['dec'],
             pixscale*np.median(obj1['fwhm'])/2)

    # Mark matched stars
    sim['detected'][sidx] = True
    # Also store measured magnitude, its error and flags
    sim['mag_measured'][sidx] = obj1['mag_calib'][oidx]
    sim['magerr_measured'][sidx] = obj1['magerr'][oidx]
    sim['flags_measured'][sidx] = obj1['flags'][oidx]

    sims.append(sim)

# Now we may stack the data from all runs into a single table
from astropy.table import vstack
sims = vstack(sims)

# ..and plot it as a histograms of magnitudes
h0,b0,_ = plt.hist(sims['mag_calib'], range=[12,22], bins=50, alpha=0.3, label='All simulated stars');
h1,b1,_ = plt.hist(sims['mag_calib'][sims['detected']], range=[12,22], bins=50, alpha=0.3, label='Detected');
h2,b2,_ = plt.hist(sims['mag_calib'][idx], range=[12,22], bins=50, alpha=0.3, label='Detected and unflagged');

plt.legend()

plt.xlabel('Injected magnitude')
plt.show()

# ..or as a detection efficiency
plt.plot(0.5*(b0[1:]+b0[:-1]), h1/h0, 'o-', label='Detected')
plt.plot(0.5*(b0[1:]+b0[:-1]), h2/h0, 'o-', label='Detected and unflagged')

plt.legend()
plt.xlabel('Injected magnitude')

plt.title('Fraction of detected artificial stars')

The complete example of detection efficiency analysis may be seen in the corresponding notebook.

Attention

Note that the flux that you specify for the artificial star is a total flux - it will not necessarily correspond to the one measured inside an aperture (especially if the aperture size is quite small). The difference (aperture correction) is due to the fraction of the total flux that falls outside of your aperture, and it has to be somehow estimated and corrected if you wish to compare the generated fluxes of injected stars with measured values.

stdpipe.psf.place_psf_stamp(image, psf, x0, y0, flux=1, gain=None)

Places PSF stamp, scaled to a given flux, at a given position inside the image.

PSF stamp is evaluated at a given position, then adjusted to accommodate for required sub-pixel shift, and finally scaled to requested flux value. Thus, the routine corresponds to injection of an artificial point source into the image.

The stamp values are added on top of current content of the image. If gain value is set, the Poissonian noise is applied to the stamp.

The image is modified in-place.

Parameters:

imagenumpy.ndarray

The image where the artificial source will be injected (modified in-place).

psfdict

Input PSF structure as returned by stdpipe.psf.run_psfex() or stdpipe.psf.load_psf().

x0float

x coordinate of the position to inject the source.

y0float

y coordinate of the position to inject the source.

fluxfloat, optional

The source flux in ADU units.

gainfloat, optional

Image gain value. If set, used to apply Poissonian noise to the source.

stdpipe.pipeline.place_random_stars(image, psf_model, nstars=100, minflux=1, maxflux=100000, gain=None, saturation=65535, edge=0, wcs=None, verbose=False)

Randomly place artificial stars into an image in-place.

Stars are added on top of the existing image content. Poissonian noise is applied according to gain, and pixel values are clipped at saturation. Coordinates are uniformly distributed; fluxes are log-uniform.

Parameters:

imagendarray

2D image array; modified in-place.

psf_modeldict

PSF model structure as returned by stdpipe.psf.run_psfex().

nstarsint, optional

Number of artificial stars to inject.

minfluxfloat, optional

Minimum star flux in ADU.

maxfluxfloat, optional

Maximum star flux in ADU.

gainfloat, optional

Detector gain in e-/ADU; used to add Poissonian noise to each injected star.

saturationfloat, optional

Saturation level in ADU; injected pixels above this are clipped.

edgeint, optional

Minimum distance to image edges for star placement.

wcsastropy.wcs.WCS, optional

If provided, sky coordinates ra and dec are added to the catalogue.

verbosebool or callable, optional

Whether to show verbose messages. May be boolean or a print-like callable.

Returns:

astropy.table.Table

Catalogue of injected stars as returned by make_random_stars(), with columns x, y, flux, mag, and optionally ra, dec.

PSF photometry with photutils

STDPipe provides PSF fitting photometry using the photutils library as an alternative to SExtractor. This approach is pure Python, more flexible, and works with various PSF models including PSFEx, Gaussian, and empirical ePSF.

The stdpipe.photometry_psf.measure_objects_psf() function performs PSF photometry at the positions of already detected objects. It supports:

• Multiple PSF model types (PSFEx, Gaussian, empirical ePSF)

• Position-dependent PSF for wide-field images

• Grouped fitting for crowded fields

• Comprehensive quality metrics

from stdpipe import photometry, photometry_psf, psf

# Detect objects
obj = photometry.get_objects_sep(image, mask=mask, thresh=3.0)

# Build PSF model
psf_model = psf.run_psfex(image, mask=mask, order=1, verbose=True)

# Perform PSF photometry
obj_psf = photometry_psf.measure_objects_psf(
    obj, image,
    psf=psf_model,
    mask=mask,
    verbose=True
)

The output table includes fitted positions (x_psf, y_psf), fluxes (flux, fluxerr), magnitudes (mag, magerr), and quality metrics:

• qfit_psf - Fit quality (0 = good, higher values indicate poor fits)

• cfit_psf - Central pixel fit quality

• flags_psf - Photutils fit flags

• npix_psf - Number of unmasked pixels used in fit

• reduced_chi2_psf - Reduced chi-squared (photutils >= 2.3.0)

stdpipe.photometry_psf.measure_objects_psf(obj, image, psf=None, psf_size=None, fwhm=None, mask=None, bg=None, err=None, gain=None, bg_size=64, bkgann=None, bkg_order=1, sn=None, fit_shape='circular', fit_size=None, maxiters=3, recentroid=True, keep_negative=True, get_bg=False, use_position_dependent_psf=False, group_sources=True, grouper_radius=None, verbose=False)

PSF photometry at the positions of already detected objects using photutils.

Performs PSF fitting photometry which is more accurate than aperture photometry, especially for point sources in crowded fields or when accurate flux measurement of PSF wings is important.

This function will estimate and subtract the background unless external background estimation (bg) is provided, and use user-provided noise map (err) if requested.

If a PSF model is not provided, a simple Gaussian PSF will be constructed based on the fwhm parameter or estimated from the data.

Parameters:

obj~astropy.table.Table

Table with initial object detections to be measured. Must have ‘x’ and ‘y’ columns.

image~numpy.ndarray

Input image as a 2D NumPy array.

psfphotutils PSF model, dict, or None, optional

PSF model to use. Can be a photutils PSF model (e.g., IntegratedGaussianPRF, FittableImageModel), a PSFEx PSF structure from stdpipe.psf.run_psfex(), or None (will create Gaussian PSF based on fwhm).

psf_sizeint or None, optional

Size of the PSF model in pixels. If None, will be estimated from PSF or set to 5*fwhm.

fwhmfloat or None, optional

Full width at half maximum in pixels. Used if PSF model is not provided, or to estimate psf_size. If None, will be estimated from obj[‘fwhm’] if available.

mask~numpy.ndarray or None, optional

Image mask as a boolean array (True values will be masked).

bg~numpy.ndarray or None, optional

If provided, use this background (same shape as input image) instead of automatically computed one.

err~numpy.ndarray or None, optional

Image noise map to be used instead of automatically computed one.

gainfloat or None, optional

Image gain in e-/ADU, used to build image noise model.

bg_sizeint, optional

Background grid size in pixels.

bkgannlist of float or None, optional

Background annulus for local background estimation, [inner_radius, outer_radius] in pixels. If None, no local background subtraction is performed (relies only on global Background2D subtraction). If set, uses gradient-aware local background fitting to handle non-uniform backgrounds. Note: radii are NOT scaled by FWHM (unlike measure_objects).

bkg_orderint, optional

Polynomial order for local background fitting. 0 = constant (mean), 1 = plane (linear gradient, recommended), 2 = quadratic surface. Only used if bkgann is set.

snfloat or None, optional

Minimal S/N ratio for the object to be considered good. If set, all measurements with magnitude errors exceeding 1/sn will be discarded.

fit_shapestr, optional

Shape of fitting region. Options: ‘circular’ (default), ‘square’. Determines the aperture used for PSF fitting.

fit_sizeint or None, optional

Size of fitting region in pixels. If None, defaults to psf_size.

maxitersint, optional

Maximum number of iterations for PSF fitting.

recentroidbool, optional

If True, allow PSF position to vary during fitting (recommended).

keep_negativebool, optional

If False, measurements with negative fluxes will be discarded.

get_bgbool, optional

If True, the routine will also return estimated background and background noise images.

use_position_dependent_psfbool, optional

If True and PSF is a PSFEx model, use polynomial evaluation for position-dependent PSF (evaluates PSF at each source position).

group_sourcesbool, optional

If True, use grouped PSF fitting for overlapping sources. Fits nearby sources simultaneously for better accuracy in crowded fields.

grouper_radiusfloat or None, optional

Radius in pixels for grouping nearby sources. If None, defaults to 2*psf_size. Only used if group_sources is True.

verbosebool or callable, optional

Whether to show verbose messages during the run. May be either boolean, or a print-like function.

Returns:

result~astropy.table.Table or tuple

Copy of original table with flux, fluxerr, mag, magerr, x_psf, y_psf columns from PSF fitting. Also includes quality of fit columns: qfit_psf (fit quality, 0=good), cfit_psf (central pixel fit quality), flags_psf (photutils fit flags), npix_psf (number of unmasked pixels used in fit), and reduced_chi2_psf (reduced chi-squared, available in photutils >= 2.3.0). If get_bg is True, returns a tuple of (table, background, background_error).

Position-dependent PSF

For wide-field images where the PSF varies across the field, enable position-dependent evaluation of PSFEx polynomial models:

# Build PSF model with spatial variation (order > 0)
psf_model = psf.run_psfex(image, mask=mask, order=2, verbose=True)

# Use position-dependent PSF
obj_psf = photometry_psf.measure_objects_psf(
    obj, image,
    psf=psf_model,
    use_position_dependent_psf=True,  # Evaluate PSF at each position
    verbose=True
)

This evaluates the PSFEx polynomial (constant, linear, quadratic terms) at each source position for more accurate photometry in fields with PSF gradients.

Grouped PSF fitting for crowded fields

In crowded fields with overlapping sources, grouped fitting simultaneously fits nearby sources to reduce contamination:

obj_psf = photometry_psf.measure_objects_psf(
    obj, image,
    psf=psf_model,
    group_sources=True,      # Enable grouped fitting
    grouper_radius=10.0,     # Group sources within 10 pixels
    verbose=True
)

# Check grouping results
print(f"Grouped {len(np.unique(obj_psf['group_id']))} groups")

This is slower but more accurate when sources are close together.

Building empirical PSF models

When PSFEx is not available or you want a purely empirical PSF, use stdpipe.psf.create_psf_model() to build an ePSF from stars in your image:

from stdpipe import psf

# Create empirical PSF from stars in the image
epsf_model = psf.create_psf_model(
    image,
    mask=mask,
    fwhm=3.0,        # Approximate FWHM
    size=25,         # Cutout size
    oversampling=2,  # ePSF oversampling factor
    verbose=True
)

# Use it just like a PSFEx model
obj_psf = photometry_psf.measure_objects_psf(
    obj, image,
    psf=epsf_model,
    verbose=True
)

The ePSF model is built by combining postage stamps of isolated stars and returns a PSFEx-compatible dictionary structure that works with all PSF functions.

stdpipe.psf.create_psf_model(image, obj=None, fwhm=None, size=None, mask=None, oversampling=None, degree=0, regularization=1e-06, subtract_neighbors=True, subtract_background=False, isolation=2.0, get_raw=False, verbose=False)

Create an empirical PSF (ePSF) model from stars in the image.

For degree=0 (default), builds a position-invariant ePSF using photutils EPSFBuilder (iterative recentering and stacking).

For degree > 0, builds a position-dependent PSF model by fitting per-pixel polynomial coefficients to resampled star stamps, following the same approach as PSFEx. The polynomial model is:

\[PSF(i,j; x,y) = \sum_k c_k(i,j) \cdot dx^{p1_k} \cdot dy^{p2_k}\]

where (dx, dy) are normalized image coordinates and (p1_k, p2_k) are polynomial exponents with p1_k + p2_k <= degree.

The returned dictionary structure is compatible with PSFEx output from stdpipe.psf.run_psfex() and can be used with the same evaluation functions like stdpipe.psf.get_psf_stamp().

Parameters:

imagenumpy.ndarray

Input image as a NumPy array, must be background subtracted.

objastropy.table.Table, optional

Table of star positions. If None, stars will be detected automatically. Should have ‘x’, ‘y’ columns and optionally ‘flux’.

fwhmfloat, optional

Approximate FWHM of stars in pixels. If None, will be estimated.

sizeint, optional

Size of cutouts to extract around stars (should be odd). If None, automatically determined from FWHM as max(15, round_up_to_odd(5 * fwhm)).

masknumpy.ndarray, optional

Image mask as a boolean array (True values will be masked).

oversamplingint, optional

Oversampling factor for the ePSF. If None (default), it is auto-selected from the FWHM: 1 when fwhm >= 2.5 image pixels (well-sampled PSF) and 2 otherwise (under-sampled PSF). Pass an explicit integer to override.

degreeint, optional

Polynomial degree for spatial PSF variation (default: 0 = constant). Degree 1 = linear (3 coefficients), degree 2 = quadratic (6 coefficients), etc.

regularizationfloat, optional

Tikhonov regularization parameter for polynomial fitting (default: 1e-6). Only used when degree > 0. Set to 0 for unregularized least-squares.

subtract_neighborsbool, optional

If True (default), subtract estimated flux from neighboring stars before extracting cutouts. Reduces contamination in crowded fields. Only used when degree > 0.

subtract_backgroundbool or {‘none’, ‘median’, ‘plane’}, optional

Local background handling for each training stamp before normalization. False or 'none' leaves the current image values unchanged, True or 'median' subtracts the median of the stamp border, and 'plane' fits and subtracts a tilted background plane from the border pixels. Only used when degree > 0. When building from the raw image instead of a background-subtracted one, 'plane' is usually the most robust choice.

isolationfloat, optional

Minimum nearest-neighbor distance in FWHM units for selecting stars for ePSF building (default: 2.0). Stars with a neighbor closer than isolation * fwhm are excluded. Combined with subtract_neighbors=True, the lower default value works in both sparse and dense fields. Set to 0 or None to disable isolation filtering.

get_rawbool, optional

If True and degree=0, returns raw photutils EPSFModel object. Ignored when degree > 0.

verbosebool or callable, optional

Whether to show verbose messages.

Returns:

dict

Dictionary with PSFEx-compatible structure containing the PSF model.

PSF photometry with DAOPHOT (IRAF)

As an alternative to the photutils backend, STDPipe also provides PSF photometry through the classic IRAF DAOPHOT workflow (phot → psf → allstar) via stdpipe.photometry_iraf.measure_objects_psf(). This requires PyRAF and IRAF to be installed.

DAOPHOT automatically selects bright, isolated stars to build a PSF model, then fits all sources using the allstar task. It supports multiple analytical PSF functions (Gaussian, Moffat, Lorentzian, Penny) and spatially variable PSF models.

from stdpipe import photometry_iraf

# Basic DAOPHOT PSF photometry (automatic PSF star selection)
result = photometry_iraf.measure_objects_psf(obj, image, fwhm=fwhm, gain=gain)

# Custom PSF function and fitting parameters
result = photometry_iraf.measure_objects_psf(obj, image, fwhm=fwhm, gain=gain,
             psf_function='moffat25', psfrad=12.0, fitrad=4.0, sn=5)

# Provide specific PSF stars
bright = np.argsort(obj['flux'])[-20:]
result = photometry_iraf.measure_objects_psf(obj, image, fwhm=fwhm,
             psf_stars=bright, verbose=True)

stdpipe.photometry_iraf.measure_objects_psf(obj, image, fwhm=None, psf_stars=None, n_psf_stars=20, mask=None, bg=None, err=None, gain=None, bg_size=64, sn=None, psfrad=None, fitrad=None, psf_function='auto', psf_varorder=0, keep_negative=True, get_bg=False, _workdir=None, _tmpdir=None, verbose=False)

PSF photometry using IRAF DAOPHOT (psf + allstar tasks).

This function performs PSF photometry using the full DAOPHOT workflow: 1. Select PSF stars (or use provided psf_stars) 2. Build PSF model using DAOPHOT psf task 3. Fit all sources using DAOPHOT allstar task 4. Return results in stdpipe format with quality metrics

The PSF model is built automatically from bright, isolated stars in the image. If psf_stars is provided, those specific stars will be used instead of automatic selection.

Note: Requires PyRAF and IRAF to be installed and properly configured.

Parameters:

objastropy.table.Table

Initial object detections. Must have ‘x’ and ‘y’ columns. Should also have ‘flux’ or ‘mag’ for PSF star selection.

imagenumpy.ndarray

Input image as a 2D NumPy array.

fwhmfloat

FWHM in pixels (required for PSF building and fitting).

psf_starsarray-like or None

Indices of objects to use as PSF stars. If None, stars are selected automatically.

n_psf_starsint

Number of PSF stars to select automatically.

masknumpy.ndarray or None

Boolean image mask (True values will be masked).

bgnumpy.ndarray or None

Background image to subtract instead of automatically computed one.

errnumpy.ndarray or None

Image noise map.

gainfloat or None

Image gain in e/ADU.

bg_sizeint

Background grid size in pixels.

snfloat or None

Minimal S/N ratio for output filtering.

psfradfloat or None

PSF radius in pixels (default: 3*fwhm). Sets how far the PSF extends.

fitradfloat or None

Fitting radius in pixels (default: fwhm). Sets how many pixels to use for fitting each star.

psf_functionstr

PSF function type: ‘auto’, ‘gauss’, ‘moffat15’, ‘moffat25’, ‘lorentz’, ‘penny1’, or ‘penny2’.

psf_varorderint

PSF variability order: 0=constant, 1=linear, 2=quadratic.

keep_negativebool

If False, discard measurements with negative fluxes.

get_bgbool

If True, also return background and error images.

_workdirstr or None

Working directory for temporary files.

_tmpdirstr or None

Parent directory for temporary directory creation.

verbosebool or callable

Verbose output.

Returns:

astropy.table.Table

Table with PSF photometry results including flux, mag, x_psf, y_psf, qfit_psf, cfit_psf, flags_psf columns. If get_bg=True, returns a tuple of (table, background, background_error).

Raises:

ImportError

If PyRAF is not available.

ValueError

If required parameters are missing.

RuntimeError

If DAOPHOT tasks fail.

PSF photometry in SExtractor (alternative)

As an alternative to the photutils-based approach described above, derived PSF models may also be used for performing PSF photometry in SExtractor.

As currently STDPipe lacks the routine for saving PSF models back to FITS files, in order to use PSF photometry you will need to pass psffile argument to stdpipe.psf.run_psfex() in order to directly store the produced PSF model into the file. Then, this file may be used in SExtractor by passing its path as a psf argument into stdpipe.photometry.get_objects_sextractor(). It will then produce the list of detected objects having a set of additional measured parameters including flux_psf, fluxerr_psf, mag_psf, magerr_psf - they correspond to the flux derived from fitting the objects with the PSF model.

# We probably do not have enough stars to study PSF spatial variance, so we use order=0 here
psf_model = psf.run_psfex(image, mask=mask, order=0, gain=gain, psffile='/tmp/psf.psf', verbose=True)

obj_psf = photometry.get_objects_sextractor(image, mask=mask, aper=3.0, edge=10, gain=gain,
             wcs=wcs, psf='/tmp/psf.psf', verbose=False)

# Now we may calibrate the photometry using PSF fluxes
m_psf = pipeline.calibrate_photometry(obj_psf, cat, sr=1/3600,
             obj_col_mag='mag_psf', obj_col_mag_err='magerr_psf',
             cat_col_mag='rmag', cat_col_mag1='gmag', cat_col_mag2='rmag',
             order=0, verbose=True)