ZogyTask

class lsst.ip.diffim.ZogyTask(config=None, name=None, parentTask=None, log=None)

Bases: lsst.pipe.base.Task

Task to perform ZOGY proper image subtraction. See module-level documentation for additional details.

Methods Summary

calculateFourierDiffim(psf1, im1, varPlane1, …) Convolve and subtract two images in Fourier space.
calculateMaskPlane(mask1, mask2[, effPsf1, …]) Calculate the mask plane of the difference image.
checkCentroids(psfArr1, psfArr2) Check whether two PSF array centroids’ distance is within tolerance.
computePsfAtCenter(exposure) Computes the PSF image at the bbox center point.
emptyMetadata() Empty (clear) the metadata for this Task and all sub-Tasks.
estimateMatchingKernelSize(psf1, psf2) Estimate the image space size of the matching kernels.
finishResultExposures(diffExp[, scoreExp]) Perform final steps on the full difference exposure result.
generateGrid(imageBox, minEdgeDims, innerBoxDims) Generate a splitting grid for an image.
getAllSchemaCatalogs() Get schema catalogs for all tasks in the hierarchy, combining the results into a single dict.
getCentroid(A) Calculate the centroid coordinates of a 2D array.
getFullMetadata() Get metadata for all tasks.
getFullName() Get the task name as a hierarchical name including parent task names.
getName() Get the name of the task.
getSchemaCatalogs() Get the schemas generated by this task.
getTaskDict() Get a dictionary of all tasks as a shallow copy.
initializeSubImage(fullExp, innerBox, …[, …]) Initializes a sub image.
inverseFftAndCropImage(imgArr, origSize[, …]) Inverse FFT and crop padding from image array.
makeField(doc) Make a lsst.pex.config.ConfigurableField for this task.
makeKernelPsfFromArray(A) Create a non spatially varying PSF from a numpy.ndarray.
makeSpatialPsf(gridPsfs) Construct a CoaddPsf based on PSFs from individual sub image solutions.
makeSubtask(name, **keyArgs) Create a subtask as a new instance as the name attribute of this task.
padAndFftImage(imgArr) Prepare and forward FFT an image array.
padCenterOriginArray(A, newShape[, …]) Zero pad an image where the origin is at the center and replace the origin to the corner as required by the periodic input of FFT.
pasteSubDiffImg(ftDiff, diffExp[, scoreExp]) Paste sub image results back into result Exposure objects.
pixelSpaceSquare(D) Square the argument in pixel space.
prepareFullExposure(exposure1, exposure2[, …]) Performs calculations that apply to the full exposures once only.
prepareSubExposure(localCutout[, psf1, …]) Perform per-sub exposure preparations.
removeNonFinitePixels(imgArr) Replace non-finite pixel values in-place.
run(exposure1, exposure2[, calculateScore]) Task entry point to perform the zogy subtraction of exposure1-exposure2.
splitBorder(innerBox, outerBox) Split the border area around the inner box into 8 disjunct boxes.
subtractImageMean(image, mask, statsControl) In-place subtraction of sigma-clipped mean of the image.
timer(name[, logLevel]) Context manager to log performance data for an arbitrary block of code.

Methods Documentation

calculateFourierDiffim(psf1, im1, varPlane1, F1, varMean1, psf2, im2, varPlane2, F2, varMean2, calculateScore=True)

Convolve and subtract two images in Fourier space.

Calculate the ZOGY proper difference image, score image and their PSFs. All input and output arrays are in Fourier space.

Parameters:
psf1, psf2 : numpy.ndarray, (self.freqSpaceShape,)

Psf arrays. Must be already in Fourier space.

im1, im2 : numpy.ndarray, (self.freqSpaceShape,)

Image arrays. Must be already in Fourier space.

varPlane1, varPlane2 : numpy.ndarray, (self.freqSpaceShape,)

Variance plane arrays respectively. Must be already in Fourier space.

varMean1, varMean2 : numpy.float > 0.

Average per-pixel noise variance in im1, im2 respectively. Used as weighing of input images. Must be greater than zero.

F1, F2 : numpy.float > 0.

Photometric scaling of the images. See eqs. (5)–(9)

calculateScore : bool, optional

If True (default), calculate and return the detection significance (score) image. Otherwise, these return fields are None.
Returns:
result : pipe.base.Struct

All arrays are in Fourier space and have shape self.freqSpaceShape.

Fd

Photometric level of D (float).

D

The difference image (numpy.ndarray [numpy.complex]).

varplaneD

Variance plane of D (numpy.ndarray [numpy.complex]).

Pd

PSF of D (numpy.ndarray [numpy.complex]).

S

Significance (score) image (numpy.ndarray [numpy.complex] or None).

varplaneS

Variance plane of S ((numpy.ndarray [numpy.complex] or None).

Ps

PSF of S (numpy.ndarray [numpy.complex]).

Notes

All array inputs and outputs are Fourier-space images with shape of self.freqSpaceShape in this method.

varMean1, varMean2 quantities are part of the noise model and not to be confused with the variance of image frequency components or with varPlane1, varPlane2 that are the Fourier transform of the variance planes.

static calculateMaskPlane(mask1, mask2, effPsf1=None, effPsf2=None)

Calculate the mask plane of the difference image.

Parameters:
mask1, maks2 : lsst.afw.image.Mask

Mask planes of the two exposures.
Returns:
diffmask : lsst.afw.image.Mask

Mask plane for the subtraction result.

Notes

TODO DM-25174 : Specification of effPsf1, effPsf2 are not yet supported.

checkCentroids(psfArr1, psfArr2)

Check whether two PSF array centroids’ distance is within tolerance.

Parameters:
psfArr1, psfArr2 : numpy.ndarray of float

Input PSF arrays to check.
Returns:
None
Raises:
ValueError:

Centroid distance exceeds config.maxPsfCentroidDist pixels.
static computePsfAtCenter(exposure)

Computes the PSF image at the bbox center point.

This may be at a fractional pixel position.

Parameters:
exposure : lsst.afw.image.Exposure

Exposure with psf.
Returns:
psfImg : lsst.afw.image.Image

Calculated psf image.
emptyMetadata()

Empty (clear) the metadata for this Task and all sub-Tasks.

static estimateMatchingKernelSize(psf1, psf2)

Estimate the image space size of the matching kernels.

Return ten times the larger Gaussian sigma estimate but at least the largest of the original psf dimensions.

Parameters:
psf1, psf2 : lsst.afw.detection.Psf

The PSFs of the two input exposures.
Returns:
size : int

Conservative estimate for matching kernel size in pixels. This is the minimum padding around the inner region at each side.
finishResultExposures(diffExp, scoreExp=None)

Perform final steps on the full difference exposure result.

Set photometric calibration, psf properties of the exposures.

Parameters:
diffExp : lsst.afw.image.Exposure

The result difference image exposure to finalize.

scoreExp : lsst.afw.image.Exposure or None

The result score exposure to finalize.
Returns:
None.
static generateGrid(imageBox, minEdgeDims, innerBoxDims, minTotalDims=None, powerOfTwo=False)

Generate a splitting grid for an image.

The inner boxes cover the input image without overlap, the edges around the inner boxes do overlap and go beyond the image at the image edges.

Parameters:
imageBox : lsst.geom.Box2I

Bounding box of the exposure to split.

minEdgeDims : lsst.geom.Extent2I

Minimum edge width in (x,y) directions each side.

innerBoxDims : lsst.geom.Extent2I

Minimum requested inner box dimensions (x,y). The actual dimensions can be larger due to rounding.

minTotalDims: `lsst.geom.Extent2I`, optional

If provided, minimum total outer dimensions (x,y). The edge will be increased until satisfied.

powerOfTwo : bool, optional

If True, the outer box dimensions should be rounded up to a power of 2 by increasing the border size. This is up to 8192, above this size, rounding up is disabled.
Returns:
boxList : list of lsst.pipe.base.Struct

innerBox, outerBox : lsst.geom.Box2I, inner boxes and overlapping border around them.

Notes

Inner box dimensions are chosen to be as uniform as they can, remainder pixels at the edge of the input will be appended to the last column/row boxes.

See diffimTests/tickets/DM-28928_spatial_grid notebooks for demonstration of this code.

This method can be used for both PARENT and LOCAL bounding boxes.

The outerBox dimensions are always even.

getAllSchemaCatalogs()

Get schema catalogs for all tasks in the hierarchy, combining the results into a single dict.

Returns:
schemacatalogs : dict

Keys are butler dataset type, values are a empty catalog (an instance of the appropriate lsst.afw.table Catalog type) for all tasks in the hierarchy, from the top-level task down through all subtasks.

Notes

This method may be called on any task in the hierarchy; it will return the same answer, regardless.

The default implementation should always suffice. If your subtask uses schemas the override Task.getSchemaCatalogs, not this method.

static getCentroid(A)

Calculate the centroid coordinates of a 2D array.

Parameters:
A : 2D numpy.ndarray of float

The input array. Must not be all exact zero.
Returns:
ycen, xcen : tuple of float

Notes

Calculates the centroid as if the array represented a 2D geometrical shape with weights per cell, allowing for “negative” weights. If sum equals to exact (float) zero, calculates centroid of absolute value array.

The geometrical center is defined as (0,0), independently of the array shape. For an odd dimension, this is the center of the center pixel, for an even dimension, this is between the two center pixels.

getFullMetadata()

Get metadata for all tasks.

Returns:
metadata : lsst.daf.base.PropertySet or TaskMetadata

The keys are the full task name. Values are metadata for the top-level task and all subtasks, sub-subtasks, etc.

Notes

The returned metadata includes timing information (if @timer.timeMethod is used) and any metadata set by the task. The name of each item consists of the full task name with . replaced by :, followed by . and the name of the item, e.g.:

topLevelTaskName:subtaskName:subsubtaskName.itemName

using : in the full task name disambiguates the rare situation that a task has a subtask and a metadata item with the same name.

getFullName()

Get the task name as a hierarchical name including parent task names.

Returns:
fullName : str

The full name consists of the name of the parent task and each subtask separated by periods. For example:

  • The full name of top-level task “top” is simply “top”.
  • The full name of subtask “sub” of top-level task “top” is “top.sub”.
  • The full name of subtask “sub2” of subtask “sub” of top-level task “top” is “top.sub.sub2”.
getName()

Get the name of the task.

Returns:
taskName : str

Name of the task.

See also

getFullName

getSchemaCatalogs()

Get the schemas generated by this task.

Returns:
schemaCatalogs : dict

Keys are butler dataset type, values are an empty catalog (an instance of the appropriate lsst.afw.table Catalog type) for this task.

See also

Task.getAllSchemaCatalogs

Notes

Warning

Subclasses that use schemas must override this method. The default implementation returns an empty dict.

This method may be called at any time after the Task is constructed, which means that all task schemas should be computed at construction time, not when data is actually processed. This reflects the philosophy that the schema should not depend on the data.

Returning catalogs rather than just schemas allows us to save e.g. slots for SourceCatalog as well.

getTaskDict()

Get a dictionary of all tasks as a shallow copy.

Returns:
taskDict : dict

Dictionary containing full task name: task object for the top-level task and all subtasks, sub-subtasks, etc.
initializeSubImage(fullExp, innerBox, outerBox, noiseMeanVar, useNoise=True)

Initializes a sub image.

Parameters:
fullExp : lsst.afw.image.Exposure

The full exposure to cut sub image from.

innerBox : lsst.geom.Box2I

The useful area of the calculation up to the whole bounding box of fullExp. fullExp must contain this box.

outerBox : lsst.geom.Box2I

The overall cutting area. outerBox must be at least 1 pixel larger than inneBox in all directions and may not be fully contained by fullExp.

noiseMeanVar : float > 0.

The noise variance level to initialize variance plane and to generate white noise for the non-overlapping region.

useNoise : bool, optional

If True, generate white noise for non-overlapping region. Otherwise, zero padding will be used in the non-overlapping region.
Returns:
result : lsst.pipe.base.Struct
  • subImg, subVarImg : lsst.afw.image.ImageD The new sub image and its sub variance plane.

Notes

innerBox, outerBox must be in the PARENT system of fullExp.

Supports the non-grid option when innerBox equals to the bounding box of fullExp.

inverseFftAndCropImage(imgArr, origSize, filtInf=None, filtNaN=None, dtype=None)

Inverse FFT and crop padding from image array.

Parameters:
imgArr : numpy.ndarray of numpy.complex

Fourier space array representing a real image.

origSize : tuple of int

Original unpadded shape tuple of the image to be cropped to.

filtInf, filtNan : numpy.ndarray of bool or int, optional

If specified, they are used as index arrays for result to set values to numpy.inf and numpy.nan respectively at these positions.

dtype : numpy.dtype, optional

Dtype of result array to cast return values to implicitly. This is to spare one array copy operation at reducing double precision to single. If None result inherits dtype of imgArr.
Returns:
result : numpy.ndarray of dtype
classmethod makeField(doc)

Make a lsst.pex.config.ConfigurableField for this task.

Parameters:
doc : str

Help text for the field.
Returns:
configurableField : lsst.pex.config.ConfigurableField

A ConfigurableField for this task.

Examples

Provides a convenient way to specify this task is a subtask of another task.

Here is an example of use:

class OtherTaskConfig(lsst.pex.config.Config):
    aSubtask = ATaskClass.makeField("brief description of task")
static makeKernelPsfFromArray(A)

Create a non spatially varying PSF from a numpy.ndarray.

Parameters:
A : numpy.ndarray

2D array to use as the new psf image. The pixels are copied.
Returns:
psfNew : lsst.meas.algorithms.KernelPsf

The constructed PSF.
makeSpatialPsf(gridPsfs)

Construct a CoaddPsf based on PSFs from individual sub image solutions.

Parameters:
gridPsfs : iterable of lsst.pipe.base.Struct

Iterable of bounding boxes (bbox) and Psf solutions (psf).
Returns:
psf : lsst.meas.algorithms.CoaddPsf

A psf constructed from the PSFs of the individual subExposures.
makeSubtask(name, **keyArgs)

Create a subtask as a new instance as the name attribute of this task.

Parameters:
name : str

Brief name of the subtask.

keyArgs

Extra keyword arguments used to construct the task. The following arguments are automatically provided and cannot be overridden:

  • “config”.
  • “parentTask”.

Notes

The subtask must be defined by Task.config.name, an instance of ConfigurableField or RegistryField.

padAndFftImage(imgArr)

Prepare and forward FFT an image array.

Parameters:
imgArr : numpy.ndarray of float

Original array. In-place modified as numpy.nan and numpy.inf are replaced by array mean.
Returns:
result : lsst.pipe.base.Struct

Notes

Save location of non-finite values for restoration, and replace them with image mean values. Re-center and zero pad array by padCenterOriginArray.

static padCenterOriginArray(A, newShape, useInverse=False, dtype=None)

Zero pad an image where the origin is at the center and replace the origin to the corner as required by the periodic input of FFT.

Implement also the inverse operation, crop the padding and re-center data.

Parameters:
A : numpy.ndarray

An array to copy from.

newShape : tuple of int

The dimensions of the resulting array. For padding, the resulting array must be larger than A in each dimension. For the inverse operation this must be the original, before padding dimensions of the array.

useInverse : bool, optional

Selector of forward, add padding, operation (False) or its inverse, crop padding, operation (True).

dtype: `numpy.dtype`, optional

Dtype of output array. Values must be implicitly castable to this type. Use to get expected result type, e.g. single float (nympy.float32). If not specified, dtype is inherited from A.
Returns:
R : numpy.ndarray

The padded or unpadded array with shape of newShape and dtype of dtype.
Raises:
ValueError : newShape dimensions must be greater than or equal to the

dimensions of A for the forward operation and less than or equal to for the inverse operation.

Notes

For odd dimensions, the splitting is rounded to put the center pixel into the new corner origin (0,0). This is to be consistent e.g. for a dirac delta kernel that is originally located at the center pixel.

pasteSubDiffImg(ftDiff, diffExp, scoreExp=None)

Paste sub image results back into result Exposure objects.

Parameters:
ftDiff : lsst.pipe.base.Struct

Result struct by calculateFourierDiffim.

diffExp : lsst.afw.image.Exposure

The result exposure to paste into the sub image result. Must be dimensions and dtype of self.fullExp1.

scoreExp : lsst.afw.image.Exposure or None

The result score exposure to paste into the sub image result. Must be dimensions and dtype of self.fullExp1. If None, the score image results are disregarded.
Returns:
None

Notes

The PSF of the score image is just to make the score image resemble a regular exposure and to study the algorithm performance.

Add an entry to the self.gridPsfs list.

gridPsfs : list of lsst.pipe.base.Struct
static pixelSpaceSquare(D)

Square the argument in pixel space.

Parameters:
D : 2D numpy.ndarray of numpy.complex

Fourier transform of a real valued array.
Returns:
R : numpy.ndarray of numpy.complex

Notes

D is to be inverse Fourier transformed, squared and then forward Fourier transformed again, i.e. an autoconvolution in Fourier space. This operation is not distributive over multiplication. pixelSpaceSquare(A*B) != pixelSpaceSquare(A)*pixelSpaceSquare(B)

prepareFullExposure(exposure1, exposure2, correctBackground=False)

Performs calculations that apply to the full exposures once only.

Parameters:
exposure1, exposure2 : lsst.afw.image.Exposure

The input exposures. Copies are made for internal calculations.

correctBackground : bool, optional

If True, subtracts sigma-clipped mean of exposures. The algorithm assumes zero expectation value at background pixels.
Returns:
None
Raises:
ValueError : If photometric calibrations are not available while

config.scaleByCalibration equals True.

Notes

Set a number of instance fields with pre-calculated values.

prepareSubExposure(localCutout, psf1=None, psf2=None, sig1=None, sig2=None)

Perform per-sub exposure preparations.

Parameters:
sig1, sig2 : float, optional

For debug purposes only, copnsider that the image may already be rescaled by the photometric calibration.

localCutout : lsst.pipe.base.Struct
psf1, psf2 : lsst.afw.detection.Psf, optional

If specified, use given psf as the sub exposure psf. For debug purposes.

sig1, sig2 : float, optional

If specified, use value as the sub-exposures’ background noise sigma value.
Returns:
None
removeNonFinitePixels(imgArr)

Replace non-finite pixel values in-place.

Save the locations of non-finite values for restoration, and replace them with image mean values.

Parameters:
imgArr : numpy.ndarray of float

The image array. Non-finite values are replaced in-place in this array.
Returns:
result : lsst.pipe.base.Struct
  • filtInf, filtNaN : numpy.ndarray of bool

    The filter of the pixel values that were inf or nan.
run(exposure1, exposure2, calculateScore=True)

Task entry point to perform the zogy subtraction of exposure1-exposure2.

Parameters:
exposure1, exposure2 : lsst.afw.image.Exposure

Two exposures warped and matched into matching pixel dimensions.

calculateScore : bool, optional

If True (default), calculate the score image and return in scoreExp.
Returns:
resultName : lsst.pipe.base.Struct
  • diffExp : lsst.afw.image.Exposure

    The Zogy difference exposure (exposure1-exposure2).

  • scoreExp : lsst.afw.image.Exposure or None

    The Zogy significance or score (S) exposure if calculateScore==True.

  • ftDiff : lsst.pipe.base.Struct

    Lower level return struct by calculateFourierDiffim with added fields from the task instance. For debug purposes.

Notes

diffExp and scoreExp always inherit their metadata from exposure1 (e.g. dtype, bbox, wcs).

The score image (S) is defined in the ZOGY paper as the detection statistic value at each pixel. In the ZOGY image model, the input images have uniform variance noises and thus S has uniform per pixel variance (though it is not scaled to 1). In Section 3.3 of the paper, there are “corrections” defined to the score image to correct the significance values for some deviations from the image model. The first of these corrections is the calculation of the variance plane of S allowing for different per pixel variance values by following the overall convolution operation on the pixels of the input images. S scaled (divided) by its corrected per pixel noise is referred as Scorr in the paper.

In the current implementation, scoreExp contains S in its image plane and the calculated (non-uniform) variance plane of S in its variance plane. scoreExp can be used directly for source detection as a likelihood image by respecting its variance plane or can be divided by the square root of the variance plane to scale detection significance values into units of sigma. S should be interpreted as a detection likelihood directly on a per-pixel basis. The calculated PSF of S is merely an indication how much the input PSFs localize point sources.

TODO DM-23855 : Implement further correction tags to the variance of scoreExp. As of DM-25174 it is not determined how important these further correction tags are.

static splitBorder(innerBox, outerBox)

Split the border area around the inner box into 8 disjunct boxes.

Parameters:
innerBox : lsst.geom.Box2I

The inner box.

outerBox : lsst.geom.Box2I

The outer box. It must be at least 1 pixel larger in each direction than the inner box.
Returns:
resultBoxes : list of 8 boxes covering the edge around innerBox
Raises:
ValueError : If outerBox is not larger than innerBox.

Notes

The border boxes do not overlap. The border is covered counter clockwise starting from lower left corner.

static subtractImageMean(image, mask, statsControl)

In-place subtraction of sigma-clipped mean of the image.

Parameters:
image : lsst.afw.image.Image

Image to manipulate. Its sigma clipped mean is in-place subtracted.

mask : lsst.afw.image.Mask

Mask to use for ignoring pixels.

statsControl : lsst.afw.math.StatisticsControl

Config of sigma clipped mean statistics calculation.
Returns:
None
Raises:
ValueError : If image mean is nan.
timer(name, logLevel=10)

Context manager to log performance data for an arbitrary block of code.

Parameters:
name : str

Name of code being timed; data will be logged using item name: Start and End.

logLevel

A logging level constant.

See also

timer.logInfo

Examples

Creating a timer context:

with self.timer("someCodeToTime"):
    pass  # code to time