stdpipe.sfft module

SFFT — Space-Frequency Fourier Transform image subtraction.

Implements spatially varying kernel fitting in a single global least-squares solve, following the approach of Hu et al. (2022, ApJ, 936, 157).

The model is:

science(x,y) = Σ_α Σ_β c_{α,β} · [P_β(x,y) · R_α(x,y)] + Σ_γ d_γ · P_γ(x,y)

where R_α is the reference shifted by kernel offset α, P_β are polynomial basis functions encoding spatial variation, and c_{α,β} / d_γ are scalar coefficients solved for globally.

The normal equations are assembled via batched BLAS matmuls over polynomial term pairs, avoiding materializing the full design matrix.

A soft kernel-sum constraint enforces that Σ_α a_α(x,y) = f(x,y) where f is a low-order polynomial modelling smooth flux-scale variation across the image.

class stdpipe.sfft.SFFTResult(diff: ndarray, model: ndarray, kernel_coeffs: ndarray, bg_coeffs: ndarray, kernel_shape: Tuple[int, int], kernel_poly_order: int, bg_poly_order: int, flux_poly_order: int, flux_poly_coeffs: ndarray, n_iter: int, rms: float, n_good: int, dmask: ndarray = None)

Bases: object

Result of SFFT image subtraction.

Attributes:

dmask

diff: ndarray

Difference image (science - model).

model: ndarray

Convolved reference + background model.

kernel_coeffs: ndarray

Kernel coefficients, shape (n_kernel, n_kpoly).

bg_coeffs: ndarray

Background coefficients, shape (n_bgpoly,).

kernel_shape: Tuple[int, int]

Kernel support size (ky, kx).

kernel_poly_order: int

Polynomial order for kernel spatial variation.

bg_poly_order: int

Polynomial order for differential background.

flux_poly_order: int

Polynomial order for the kernel-sum (flux scale) constraint.

flux_poly_coeffs: ndarray

Fitted flux-scale polynomial coefficients.

n_iter: int

Number of solver passes performed (sigma-clipping iterations plus the template-noise reweighting pass, if any).

rms: float

Final RMS of residuals over good (unmasked, unclipped) pixels.

n_good: int

Number of good (unmasked, unclipped) pixels in final iteration.

dmask: ndarray = None

Boolean mask (True = unreliable) of difference pixels whose model is built from missing template data: image edges and pixels within the kernel footprint of template defects (NaN/Inf or template_mask).

stdpipe.sfft.solve(image, template, mask=None, template_mask=None, err=None, template_err=None, kernel_shape=(7, 7), kernel_poly_order=2, bg_poly_order=2, flux_poly_order=1, flux_penalty=1000.0, ridge=1e-06, sigma_clip=3.0, max_iter=5, verbose=False)

SFFT image subtraction with spatially varying kernel.

Solves for a spatially varying convolution kernel and differential background in a single global least-squares problem. The kernel at each pixel is modelled as a delta-function basis with polynomial spatial variation. A soft constraint enforces that the kernel sum varies smoothly as a low-order polynomial (modelling flux-scale differences between science and template).

Iterative sigma-clipping rejects outlier pixels (transients, cosmic rays, artifacts) that would otherwise bias the kernel solution.

Since the template acts as the basis of the model, any template defect (NaN/Inf or template_mask) corrupts the model within a full kernel footprint around it. Such pixels are excluded from the fit, and flagged in the dmask attribute of the result — the difference image values there (and in the kernel half-width band along image edges) are unreliable and should not be searched for transients.

Parameters:

imagenumpy.ndarray

Science image as a 2-D NumPy array.

templatenumpy.ndarray

Template/reference image, same shape, aligned to science.

masknumpy.ndarray, optional

Boolean mask where True = bad pixels to exclude from the fit (science-side defects, or regions to ignore). Unlike template defects, these only affect their own pixel.

template_masknumpy.ndarray, optional

Mask of template defects (True = bad). Grown by the kernel footprint, since the model at every pixel within a kernel half-width of a template defect depends on the missing data. Pass template-side defects here rather than merging them into mask.

errnumpy.ndarray, optional

Per-pixel error (standard deviation) map for inverse-variance weighting. If None, uniform weights are used. A constant error map (or True) is equivalent to uniform weighting and does not change the solution (unless template_err is also given).

template_errnumpy.ndarray, optional

Per-pixel error map of the template. If given, one reweighting pass is performed: after an initial solve, the template variance is propagated through the fitted kernel and the weights are rebuilt as 1 / (err² + Σ_α a_α² · template_err²), followed by a re-solve. This properly de-weights pixels dominated by template noise and matters when the template is not much deeper than the science image. Costs one extra solver pass.

kernel_shapetuple of int, optional

(ky, kx) size of the convolution kernel, must be odd. Default (7, 7). Larger kernels handle bigger PSF differences but are slower.

kernel_poly_orderint, optional

Polynomial order for spatial variation of each kernel coefficient. Default 2 (quadratic). Higher orders capture more complex PSF variation but need more pixels.

bg_poly_orderint, optional

Polynomial order for the differential background model. Default 2.

flux_poly_orderint, optional

Polynomial order for the kernel-sum constraint (flux scale variation). Default 1 (linear gradient). Set to 0 for constant flux scale.

flux_penaltyfloat, optional

Dimensionless penalty weight for the kernel-sum constraint, relative to the median diagonal of the normal matrix (so its effect does not depend on image size or flux units). Default 1e3, which enforces the constraint nearly strictly. Set to 0 to disable the constraint entirely.

ridgefloat, optional

Dimensionless Tikhonov regularization, relative to the median diagonal of the normal matrix. Default 1e-6.

sigma_clipfloat, optional

Sigma threshold for iterative outlier rejection. Default 3.0. Set to None or 0 to disable clipping.

max_iterint, optional

Maximum number of sigma-clipping iterations. Default 5. Must be at least 1. The reweighting pass triggered by template_err is performed in addition to these.

verbosebool or callable, optional

If True, print progress. If callable, use as log function.

Returns:

SFFTResult

Result object with difference image and all fit metadata. Pixels flagged in its dmask attribute are unreliable in the difference.

stdpipe.sfft.evaluate_kernel_at(result, x, y, image_shape)

Evaluate the spatially varying kernel at a single image position.

The returned stamp is the convolution kernel in the usual sense: model(p) = Σ k[dy, dx] · template(p - (dy, dx)), with the element k[hy + dy, hx + dx] at array index (hy + dy, hx + dx) for half-sizes hy, hx.

Parameters:

resultSFFTResult

SFFTResult from solve().

xfloat

X pixel coordinate.

yfloat

Y pixel coordinate.

image_shapetuple of int

(ny, nx) of the original image.

Returns:

numpy.ndarray

2-D array of shape result.kernel_shape.

stdpipe.sfft.evaluate_flux_scale(result, x, y, image_shape)

Evaluate the flux-scale polynomial at image position(s).

Parameters:

resultSFFTResult

SFFTResult from solve().

xfloat or array-like

X coordinate(s), scalar or array.

yfloat or array-like

Y coordinate(s), scalar or array.

image_shapetuple of int

(ny, nx) of the original image.

Returns:

float or numpy.ndarray

Flux scale value(s), same shape as x/y.

stdpipe.sfft.evaluate_background(result, x, y, image_shape)

Evaluate the differential background model at image position(s).

Parameters:

resultSFFTResult

SFFTResult from solve().

xfloat or array-like

X coordinate(s), scalar or array.

yfloat or array-like

Y coordinate(s), scalar or array.

image_shapetuple of int

(ny, nx) of the original image.

Returns:

float or numpy.ndarray

Background value(s), same shape as x/y.