Image subtraction

Any attempt of detection of variable or transient sources in a reasonably dense field, or on typical background for extragalactic transients (read - on the outskirts of galaxies, often bright enough) requires image subtraction methods.

Template images

Image subtraction requires template image that has to be astrometrically aligned with your science one. If you have your own set of deep enough images you may construct the template by the methods described in the Stacking the images section above.

STDPipe also has a couple of functions that may help you downloading template images from publicly available archives. All of them will deliver ready to use image that is already projected on the pixel grid of your science frame and has the same shape.

# Get r band image from PanSTARRS with science image original resolution and orientation
tmpl = templates.get_hips_image('PanSTARRS/DR1/r', wcs=wcs,
             width=image.shape[1], height=image.shape[0], get_header=False, verbose=True)

# Now mask some brighter stars in the template as we know they are saturated in Pan-STARRS
tmask = templates.mask_template(tmpl, cat, cat_col_mag='rmag', cat_saturation_mag=15,
             wcs=wcs, dilate=3, verbose=True)

# And now the same - using original Pan-STARRS images and masks
tmpl,tmask = templates.get_ps1_image_and_mask('r', wcs=wcs, width=, height=, verbose=True)

Using Hierarchical Progressive Survey (HiPS) images

First one, stdpipe.templates.get_hips_image(), acquires template images from any HiPS (Hierarchical Progressive Survey) formatted survey available on the net (see the full list - it is huge!). The routine uses CDS hips2fits service to do the job.

stdpipe.templates.get_hips_image(hips, ra=None, dec=None, width=None, height=None, fov=None, wcs=None, shape=None, header=None, asinh=None, normalize=True, upscale=False, get_header=True, verbose=False)[source]

Load an image from any HiPS survey using the CDS hips2fits service.

The pixel grid may be specified by center coordinates and field of view, or by passing a WCS solution or FITS header directly.

Parameters:
hipsstr

HiPS survey identifier. See https://aladin.u-strasbg.fr/hips/list for the full list.

rafloat, optional

Image center Right Ascension in degrees.

decfloat, optional

Image center Declination in degrees.

widthint, optional

Image width in pixels.

heightint, optional

Image height in pixels.

fovfloat, optional

Field of view (angular size) in degrees.

wcsastropy.wcs.WCS, optional

WCS defining the pixel grid; supersedes ra, dec, fov.

shapetuple of int, optional

Image shape (height, width); alternative to width / height.

headerastropy.io.fits.Header, optional

Header containing image dimensions and WCS; alternative to providing them explicitly.

asinhbool or None, optional

Whether the survey uses non-linear asinh flux scaling. Auto-detected for Pan-STARRS surveys; set to False to override.

normalizebool, optional

If True, pseudo-normalize the image so that background mean ≈ 100 and background RMS ≈ 10.

upscalebool or int, optional

If set, request an upscaled image then downscale before returning. Useful for Pan-STARRS to reduce asinh photometric errors at coarse pixel scales.

get_headerbool, optional

If True, also return the FITS header alongside the image.

verbosebool or callable, optional

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

Returns:
ndarray or tuple or None

Image projected onto the requested pixel grid, or (image, header) if get_header=True, or None on failure.

Attention

Depending on the survey, the results may wary. Be aware that masked pixels may either be represented as NaN in the image, or be silently interpolated, that may impact the analysis of these images.

Also, right now Pan-STARRS images are in the process of re-uploading to CDS that will supposedly fix the problems due to its orininal non-lineary flux scaling. So expect the routine to produce nonsensical output for it!

We have a convenience function that may help masking the pixels that are most probably unreliable in a given HiPS image (e.g. above survey saturation theshold) if they are somehow not masked (not set to NaN).

stdpipe.templates.mask_template(tmpl, cat=None, cat_saturation_mag=None, cat_col_mag='rmag', cat_col_mag_err='e_rmag', cat_col_ra='RAJ2000', cat_col_dec='DEJ2000', cat_sr=0.0002777777777777778, mask_nans=True, mask_masked=True, mask_photometric=False, aper=2, sn=5, wcs=None, dilate=5, verbose=False, _tmpdir=None)[source]

Apply masking heuristics to a template image.

The following masking steps are applied (where enabled):

  1. NaN pixels are masked if mask_nans=True.

  2. If cat is provided and wcs is set:

    • Stars brighter than cat_saturation_mag have their central pixel masked.

    • Stars with masked or NaN magnitudes are masked if mask_masked=True.

    • Photometrically saturated stars are masked if mask_photometric=True.

  3. The mask is optionally dilated by dilate pixels.

Parameters:
tmplndarray

Template image array.

catastropy.table.Table, optional

Reference catalogue for catalogue-based masking.

cat_saturation_magfloat, optional

Stars brighter than this magnitude have their central pixel masked.

cat_col_magstr, optional

Catalogue column name for magnitude.

cat_col_mag_errstr, optional

Catalogue column name for magnitude error.

cat_col_rastr, optional

Catalogue column name for Right Ascension.

cat_col_decstr, optional

Catalogue column name for Declination.

cat_srfloat, optional

Catalogue matching radius in degrees.

mask_nansbool, optional

If True, mask NaN pixels in the template.

mask_maskedbool, optional

If True, mask catalogue stars whose magnitudes are masked or NaN.

mask_photometricbool, optional

If True, apply photometric-based masking to identify saturated stars.

aperfloat, optional

Aperture radius (pixels) for photometry in photometric masking.

snfloat, optional

Minimum S/N for star detection in photometric masking.

wcsastropy.wcs.WCS, optional

WCS of the template image.

dilateint, optional

Dilation kernel size in pixels for expanding masked regions.

verbosebool or callable, optional

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

_tmpdirstr, optional

Directory for temporary files created during photometric masking.

Returns:
ndarray of bool

Boolean mask for the template image (True = masked).

Using original Pan-STARRS or Legacy Survey images

STDPipe is also able to directly download the images from Pan-STARRS or Legacy Survey image archives, mosaic them and project onto requested pixel grid in order to produce the template. It is also able to simultaneously acquite the mask image so that you may properly exclude unreliable template pixels from the analysis.

stdpipe.templates.get_survey_image(band='r', ext='image', survey='ps1', wcs=None, shape=None, width=None, height=None, header=None, reproject='lanczos', extra=None, _cachedir=None, _cache_downscale=1, _tmpdir=None, _workdir=None, verbose=False, **kwargs)[source]

Download and reproject a Pan-STARRS or Legacy Survey image onto a WCS grid.

Pan-STARRS images are converted from asinh scaling to linear flux. Pan-STARRS mask bits: https://outerspace.stsci.edu/display/PANSTARRS/PS1+Pixel+flags+in+Image+Table+Data Legacy Survey mask bits: https://www.legacysurvey.org/dr10/bitmasks

Parameters:
bandstr, optional

Photometric band: 'g', 'r', 'i', 'z', or 'y'.

extstr, optional

Data type: 'image' or 'mask'.

surveystr, optional

Survey name: 'ps1' for Pan-STARRS or 'ls' for Legacy Survey.

wcsastropy.wcs.WCS, optional

Output WCS projection.

shapetuple of int, optional

Output image shape (height, width); alternative to width / height.

widthint, optional

Output image width in pixels.

heightint, optional

Output image height in pixels.

headerastropy.io.fits.Header, optional

Header containing WCS and image dimensions; alternative to providing them explicitly.

reprojectstr, optional

Reprojection method: 'lanczos' (default, pure-Python with flux conservation) or 'swarp' (external SWarp binary).

extradict, optional

Extra SWarp parameters (only used when reproject='swarp').

_cachedirstr, optional

Directory for caching downloaded sky cells between calls.

_cache_downscaleint, optional

Integer downscale factor for cached images.

_tmpdirstr, optional

Directory for temporary files.

_workdirstr, optional

If specified, temporary SWarp files are kept for debugging.

verbosebool or callable, optional

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

**kwargs

Additional keyword arguments passed to the reprojection function.

Returns:
ndarray or None

Image reprojected onto the requested pixel grid, or None on failure.

stdpipe.templates.get_survey_image_and_mask(band='r', **kwargs)[source]

Download both the image and mask from Pan-STARRS or Legacy Survey.

Convenience wrapper that calls get_survey_image() twice.

Parameters:
bandstr, optional

Photometric band: 'g', 'r', 'i', 'z', or 'y'.

**kwargs

All other keyword arguments are passed to get_survey_image().

Returns:
tuple of ndarray

(image, mask) reprojected onto the requested pixel grid.

stdpipe.templates.get_ps1_image(band='r', **kwargs)[source]
stdpipe.templates.get_ps1_image_and_mask(band='r', **kwargs)[source]

Running image subtraction

STDPipe provides several image subtraction methods:

All three wrappers follow the same interface conventions: they accept science and template images as NumPy arrays, optional masks and noise maps, and return the difference image along with optional noise maps, noise-normalized images, and convolved templates.

SFFT subtraction

stdpipe.subtraction.run_sfft() performs image subtraction using the SFFT algorithm (Hu et al. 2022, ApJ, 936, 157). It solves for a spatially varying convolution kernel and differential background in a single global least-squares problem, without requiring any external tools.

The kernel at each pixel is modelled as a delta-function basis with polynomial spatial variation. A soft constraint enforces that the kernel sum (flux scale) varies smoothly as a low-order polynomial, accounting for differences in photometric calibration between science and template.

# Basic subtraction with automatic noise estimation
diff = subtraction.run_sfft(image, tmpl,
             mask=mask, template_mask=tmask,
             err=True, template_err=True,
             image_gain=gain, template_gain=1e6,
             verbose=True)

# Get all output planes (like HOTPANTS)
diff, conv, noise, scaled, sfft_result = subtraction.run_sfft(image, tmpl,
             mask=mask, template_mask=tmask,
             err=True, template_err=True,
             image_gain=gain, template_gain=1e6,
             get_convolved=True, get_noise=True,
             get_scaled=True, get_kernel=True,
             verbose=True)

# Now we have:
# - `diff` for the difference image
# - `conv` for the template convolved to match the science image
# - `noise` for the per-pixel noise map of the difference image
# - `scaled` for the noise-normalized difference image (diff / noise)
# - `sfft_result` for the full SFFTResult with kernel coefficients and metadata

# Inspect the spatially varying kernel at a specific position
from stdpipe import sfft
kernel_center = sfft.evaluate_kernel_at(sfft_result, 500, 500, image.shape)
flux_scale = sfft.evaluate_flux_scale(sfft_result, 500, 500, image.shape)
stdpipe.subtraction.run_sfft(image, template, mask=None, template_mask=None, err=None, template_err=None, image_gain=None, template_gain=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, obj=None, obj_size=None, get_convolved=False, get_scaled=False, get_noise=False, get_kernel=False, verbose=False)[source]

Image subtraction using SFFT with spatially varying kernel.

Wrapper around stdpipe.sfft.solve() that handles noise model construction, difference image noise propagation, and follows the same interface conventions as run_hotpants() and run_zogy().

The noise model for both science and template images may be provided directly (as arrays), or built automatically from the images using background estimation and gain values (same approach as HOTPANTS). If err or template_err is set to True, the routine builds a noise map from the image using SEP background estimation and the corresponding gain parameter.

Parameters:
imagenumpy.ndarray

Input science image as a NumPy array.

templatenumpy.ndarray

Input template image, same shape as science image.

masknumpy.ndarray, optional

Science image mask (True = bad).

template_masknumpy.ndarray, optional

Template image mask (True = bad).

errnumpy.ndarray or bool, optional

Science image noise map (per-pixel RMS). If True, built from the image and image_gain.

template_errnumpy.ndarray or bool, optional

Template image noise map. If True, built from the template and template_gain.

image_gainfloat, optional

Gain of science image in e-/ADU.

template_gainfloat, optional

Gain of template image in e-/ADU.

kernel_shapetuple of int, optional

(ky, kx) kernel size, must be odd. Default (7, 7).

kernel_poly_orderint, optional

Polynomial order for kernel spatial variation. Default 2.

bg_poly_orderint, optional

Polynomial order for differential background. Default 2.

flux_poly_orderint, optional

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

flux_penaltyfloat, optional

Penalty weight for kernel-sum constraint. Default 1e3.

ridgefloat, optional

Tikhonov regularization. Default 1e-6.

sigma_clipfloat, optional

Sigma threshold for outlier rejection. Default 3.0.

max_iterint, optional

Maximum sigma-clipping iterations. Default 5.

objastropy.table.Table, optional

List of detected objects with x, y columns. If provided, kernel fitting is restricted to rectangular regions around these objects, analogous to HOTPANTS stamp placement. This concentrates the fit on high-S/N regions and can improve kernel quality.

obj_sizeint, optional

Half-size in pixels of the fitting region around each object. If None (default), derived from the median fwhm column of obj (3 × median(fwhm)), or 15 pixels if fwhm is not available.

get_convolvedbool, optional

Whether to also return the convolved template.

get_scaledbool, optional

Whether to also return noise-normalized difference.

get_noisebool, optional

Whether to also return the difference image noise map.

get_kernelbool, optional

Whether to also return the SFFTResult object containing kernel coefficients and fit metadata.

verbosebool or callable, optional

Whether to show verbose messages.

Returns:
numpy.ndarray or list

The difference image, or a list [diff, ...] if any get_* option is set.

HOTPANTS subtraction

stdpipe.subtraction.run_hotpants() is a wrapper for the HOTPANTS image subtraction code. It requires HOTPANTS to be installed separately (see Installation). We recommend checking the HOTPANTS documentation to better understand the concepts and options for it.

# Run the subtraction getting back all possible image planes, assuming the template
# to be noise-less, and estimating image noise model from its statistics.
# And also pre-flatting the images before subtraction to get rid
# of possible background inhomogeneities.

import photutils

bg = photutils.Background2D(image, 128, mask=mask, exclude_percentile=30).background
tbg = photutils.Background2D(tmpl, 128, mask=tmask, exclude_percentile=30).background

diff,conv,sdiff,ediff = subtraction.run_hotpants(image-bg, tmpl-tbg,
             mask=mask, template_mask=tmask,
             get_convolved=True, get_scaled=True, get_noise=True,
             image_fwhm=fwhm, template_fwhm=1.5,
             image_gain=gain, template_gain=1e6,
             err=True, verbose=True)

# Now we have:
# - `diff` for the difference image
# - `conv` for the template convolved to match the original image
# - `sdiff` for noise-normalized difference image - ideal for quickly seeing significant deviations!
# - `ediff` for the difference image noise model - you may use it to weight object detection to reject the subtraction artefacts e.g. around brighter objects
stdpipe.subtraction.run_hotpants(image, template, mask=None, template_mask=None, err=None, template_err=None, extra=None, image_fwhm=None, template_fwhm=None, image_gain=None, template_gain=1000, rel_r=3, rel_rss=4, nx=1, ny=None, obj=None, get_convolved=False, get_scaled=False, get_noise=False, get_kernel=False, get_header=False, _tmpdir=None, _workdir=None, _exe=None, verbose=False)[source]

Wrapper for running HOTPANTS to subtract a template from science image.

To better understand the logic and available subtraction options you may check HOTPANTS documentation at https://github.com/acbecker/hotpants

The routine tries to be clever and select reasonable defaults for running HOTPANTS, but you may always directly specify any HOTPANTS option manually through extra parameter.

The noise model for both science and template images may optionally be provided as an external parameters (err and template_err). Moreover, if these parameters are set to True, the routine will try to build noise models itself. To do so, it will estimate the image background and background RMS. Then, this background RMS map is used as a primary component of the noise map. On top of that, the contribution of the sources is added by subtracting the background from original image, smoothing it a bit to mitigate pixel-level noise a bit, and then converting everything above zero to corresponding Poissonian noise contribution by dividing with gain value and taking a square root.

rel_r and rel_rss define the scaling of corresponding HOTPANTS parameters (which are the convolution kernel half width and half-width of the sub-stamps used for kernel fitting) with image FWHM.

The routine also uses “officially recommedned” logic for defining HOTPANTS ng parameter, that is to set it to [3, 6, 0.5*sigma_match, 4, 1.0*sigma_match, 2, 2.0*sigma_match] where sigma_match = np.sqrt(image_fwhm**2 - template_fwhm**2) / 2.35. This approach assumes that image_fwhm > template_fwhm, and of course needs image_fwhm and template_fwhm to be provided.

Parameters:
imagenumpy.ndarray

Input science image as a NumPy array.

templatenumpy.ndarray

Input template image, should have the same shape as the science image.

masknumpy.ndarray, optional

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

template_masknumpy.ndarray, optional

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

errnumpy.ndarray or bool, optional

Science image error map (expected RMS of every pixel). If set to True, the code will try to build the noise map directly from the image and image_gain parameter.

template_errnumpy.ndarray or bool, optional

Template image error map. If set to True, the code will try to build the noise map directly from the template and template_gain parameter.

extradict, optional

Extra parameters to be passed to HOTPANTS executable, with parameter names as keys.

image_fwhmfloat, optional

FWHM of science image in pixels.

template_fwhmfloat, optional

FWHM of template image in pixels.

image_gainfloat, optional

Gain of science image.

template_gainfloat, optional

Gain of template image.

rel_rfloat, optional

If specified, HOTPANTS r parameters will be set to image_fwhm*rel_r.

rel_rssfloat, optional

If specified, HOTPANTS rss parameters will be set to image_fwhm*rel_rss.

nxint, optional

Number of image sub-regions in x direction.

nyint, optional

Number of image sub-regions in y direction.

objastropy.table.Table, optional

List of objects detected on science image. If provided, it will be used to better place the sub-stamps used to derive the convolution kernel.

get_convolvedbool, optional

Whether to also return the convolved template image.

get_scaledbool, optional

Whether to also return noise-normalized difference image.

get_noisebool, optional

Whether to also return difference image noise model.

get_kernelbool, optional

Whether to also return convolution kernel used in the subtraction.

get_headerbool, optional

Whether to also return the FITS header from HOTPANTS output file.

_workdirstr, optional

If specified, all temporary files will be created in this directory, and will be kept intact after running HOTPANTS. May be used for debugging.

_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 HOTPANTS 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:
numpy.ndarray or list

The difference image and, optionally, other kinds of images as requested by the get_* options above.