Skip to content

DM-51368: Move kernel candidate residual metric to utils #408

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
isullivan merged 1 commit into from
Jun 13, 2025
Merged

DM-51368: Move kernel candidate residual metric to utils #408

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 37 additions & 78 deletions python/lsst/ip/diffim/detectAndMeasure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@
import lsst.afw.table as afwTable
import lsst.daf.base as dafBase
import lsst.geom
from lsst.ip.diffim.utils import evaluateMaskFraction
from lsst.ip.diffim.utils import evaluateMaskFraction, computeDifferenceImageMetrics
from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask
from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
import lsst.meas.deblender
Expand Down Expand Up @@ -68,6 +68,12 @@ class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections,
storageClass="ExposureF",
name="{fakesType}{coaddName}Diff_differenceTempExp",
)
kernelSources = pipeBase.connectionTypes.Input(
doc="Final selection of sources used for psf matching.",
dimensions=("instrument", "visit", "detector"),
storageClass="SourceCatalog",
name="{fakesType}{coaddName}Diff_psfMatchSources"
)
outputSchema = pipeBase.connectionTypes.InitOutput(
doc="Schema (as an example catalog) for output DIASource catalog.",
storageClass="SourceCatalog",
Expand Down Expand Up @@ -390,7 +396,7 @@ def runQuantum(self, butlerQC: pipeBase.QuantumContext,
butlerQC.put(outputs, outputRefs)

@timeMethod
def run(self, science, matchedTemplate, difference,
def run(self, science, matchedTemplate, difference, kernelSources,
Copy link
Contributor

@fred3m fred3m Jun 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't forget to add kernelSources to the docstring.

idFactory=None):
"""Detect and measure sources on a difference image.

Expand All @@ -411,6 +417,8 @@ def run(self, science, matchedTemplate, difference,
difference image.
difference : `lsst.afw.image.ExposureF`
Result of subtracting template from the science image.
kernelSources : `lsst.afw.table.SourceCatalog`
Final selection of sources that was used for psf matching.
idFactory : `lsst.afw.table.IdFactory`, optional
Generator object used to assign ids to detected sources in the
difference image. Ids from this generator are not set until after
Expand Down Expand Up @@ -473,7 +481,7 @@ def run(self, science, matchedTemplate, difference,
results.negative)

result = self.processResults(science, matchedTemplate, difference,
sources, idFactory,
sources, idFactory, kernelSources,
positives=positives,
negatives=negatives)
result.differenceBackground = background
Expand All @@ -483,7 +491,7 @@ def run(self, science, matchedTemplate, difference,
negatives = afwTable.SourceCatalog(self.schema)
results.negative.makeSources(negatives)
result = self.processResults(science, matchedTemplate, difference,
results.sources, idFactory,
results.sources, idFactory, kernelSources,
positives=positives,
negatives=negatives)
result.differenceBackground = background
Expand Down Expand Up @@ -514,7 +522,7 @@ def _prepareInputs(self, difference):
mask.addMaskPlane(mp)
mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)

def processResults(self, science, matchedTemplate, difference, sources, idFactory,
def processResults(self, science, matchedTemplate, difference, sources, idFactory, kernelSources,
Copy link
Contributor

@fred3m fred3m Jun 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please update this docstring too.

positives=None, negatives=None,):
"""Measure and process the results of source detection.

Expand All @@ -532,6 +540,8 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
idFactory : `lsst.afw.table.IdFactory`
Generator object used to assign ids to detected sources in the
difference image.
kernelSources : `lsst.afw.table.SourceCatalog`
Final selection of sources that was used for psf matching.
positives : `lsst.afw.table.SourceCatalog`, optional
Positive polarity footprints.
negatives : `lsst.afw.table.SourceCatalog`, optional
Expand Down Expand Up @@ -597,7 +607,7 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
if self.config.doForcedMeasurement:
self.measureForcedSources(diaSources, science, difference.getWcs())

self.calculateMetrics(science, difference, diaSources)
self.calculateMetrics(science, difference, diaSources, kernelSources)

measurementResults = pipeBase.Struct(
subtractedMeasuredExposure=difference,
Expand Down Expand Up @@ -787,20 +797,22 @@ def measureForcedSources(self, diaSources, science, wcs):
for diaSource, forcedSource in zip(diaSources, forcedSources):
diaSource.assign(forcedSource, mapper)

def calculateMetrics(self, science, difference, diaSources):
def calculateMetrics(self, science, difference, diaSources, kernelSources):
Copy link
Contributor

@fred3m fred3m Jun 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update docstring

"""Add difference image QA metrics to the Task metadata.

This may be used to produce corresponding metrics (see
lsst.analysis.tools.tasks.diffimTaskDetectorVisitMetricAnalysis).

Parameters
----------
science : `lsst.afw.image.ExposureF`
science : `lsst.afw.image.ExposureF`
Science exposure that was subtracted.
difference : `lsst.afw.image.Exposure`
The target difference image to calculate metrics for.
diaSources : `lsst.afw.table.SourceCatalog`
The catalog of detected sources.
kernelSources : `lsst.afw.table.SourceCatalog`
Final selection of sources that was used for psf matching.
"""
mask = difference.mask
badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
Expand All @@ -823,75 +835,20 @@ def calculateMetrics(self, science, difference, diaSources):
self.metadata["%s_mask_fraction"%maskPlane.lower()] = -1
self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)

residualInfo = self._computeDifferenceResidual(science, difference, diaSources)
self.metadata["residualFootprintRatioMean"] = residualInfo.ratio
self.metadata["residualFootprintRatioStdev"] = residualInfo.ratio_std
if self.config.doSkySources:
skySources = diaSources[diaSources["sky_source"]]
else:
skySources = None
metrics = computeDifferenceImageMetrics(science, difference, kernelSources, sky_sources=skySources)

self.metadata["residualFootprintRatioMean"] = metrics.differenceFootprintRatioMean
self.metadata["residualFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
self.log.info("Mean, stdev of ratio of difference to science "
"pixels in star footprints: %5.4f, %5.4f",
residualInfo.ratio, residualInfo.ratio_std)

def _computeDifferenceResidual(self, science, difference, diaSources):
"""Compute quality metrics (saved to the task metadata) on the
difference image, at the locations of detected stars on the science
image (determined by the DETECTED mask plane). This restricts the metric
to locations that should be well-subtracted.

Parameters
----------
science : `lsst.afw.image.ExposureF`
Science exposure that was subtracted.
difference : `lsst.afw.image.ExposureF`
Result of subtracting template and science.
diaSources : `lsst.afw.table.SourceCatalog`
The catalog of detected sources.

Returns
-------
results: `lsst.pipe.base.Struct`
``ratio`` : `np.ndarray`
Mean of the ratio of the absolute value of the difference image
(with the mean absolute value of the sky regions on the
difference image removed) to the science image, computed in the
footprints of stars detected on the science image.
``ratio_std`` : `np.ndarray`
Standard Deviation across footprints of the above ratio.
"""
good = science.mask.array & science.mask.getPlaneBitMask('DETECTED') > 0
badPix = (difference.mask.array
& difference.mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
good &= ~badPix
if not self.config.doSkySources:
# Place sky sources to sample the background
skySourceFootprints = self.skySources.run(mask=science.mask, seed=difference.info.id,
catalog=afwTable.SourceCatalog(self.schema))
else:
skySourceFootprints = [sky.getFootprint() for sky in diaSources[diaSources["sky_source"]]]
n = len(skySourceFootprints)
sky_science = np.zeros(n)
sky_difference = np.zeros(n)
for i, footprint in enumerate(skySourceFootprints):
heavy = lsst.afw.detection.makeHeavyFootprint(footprint, science.maskedImage)
heavy_diff = lsst.afw.detection.makeHeavyFootprint(footprint, difference.maskedImage)
sky_science[i] = abs(heavy.getImageArray()).mean()
sky_difference[i] = abs(heavy_diff.getImageArray()).mean()

# Calculate the difference between the absolute value of the science
# image within detected footprints and the sky sources.
# Then do the same for the difference image
if np.sum(good) > 1:
science_mean = np.abs(np.mean(abs(science.image.array[good])) - sky_science.mean())
science_std = np.std(abs(science.image.array[good]))
difference_mean = np.abs(np.mean(abs(difference.image.array[good])) - sky_difference.mean())
difference_std = np.std(abs(difference.image.array[good]))
# Calculate stddev of the science and difference image values and
# propagate to the stddev of the quotient.
# This avoids potential divide by zero errors.
ratio = difference_mean/science_mean
ratio_std = ratio*np.sqrt((difference_std/difference_mean)**2 + (science_std/science_mean)**2)
else:
ratio = np.nan
ratio_std = np.nan
return pipeBase.Struct(ratio=ratio, ratio_std=ratio_std)
self.metadata["residualFootprintRatioMean"],
self.metadata["residualFootprintRatioStdev"])

def _runStreakMasking(self, difference):
"""Do streak masking and optionally save the resulting streak
Expand Down Expand Up @@ -995,7 +952,7 @@ class DetectAndMeasureScoreTask(DetectAndMeasureTask):
_DefaultName = "detectAndMeasureScore"

@timeMethod
def run(self, science, matchedTemplate, difference, scoreExposure,
def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources,
Copy link
Contributor

@fred3m fred3m Jun 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update docstring

idFactory=None):
"""Detect and measure sources on a score image.

Expand All @@ -1010,6 +967,8 @@ def run(self, science, matchedTemplate, difference, scoreExposure,
Result of subtracting template from the science image.
scoreExposure : `lsst.afw.image.ExposureF`
Score or maximum likelihood difference image
kernelSources : `lsst.afw.table.SourceCatalog`
Final selection of sources that was used for psf matching.
idFactory : `lsst.afw.table.IdFactory`, optional
Generator object used to assign ids to detected sources in the
difference image. Ids from this generator are not set until after
Expand Down Expand Up @@ -1047,7 +1006,7 @@ def run(self, science, matchedTemplate, difference, scoreExposure,
results.negative)

return self.processResults(science, matchedTemplate, difference,
sources, idFactory,
sources, idFactory, kernelSources,
positives=positives,
negatives=negatives)

Expand All @@ -1057,6 +1016,6 @@ def run(self, science, matchedTemplate, difference, scoreExposure,
negatives = afwTable.SourceCatalog(self.schema)
results.negative.makeSources(negatives)
return self.processResults(science, matchedTemplate, difference,
results.sources, idFactory,
results.sources, idFactory, kernelSources,
positives=positives,
negatives=negatives)
85 changes: 8 additions & 77 deletions python/lsst/ip/diffim/subtractImages.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
import lsst.afw.image
import lsst.afw.math
import lsst.geom
from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm
from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics
from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
import lsst.pex.config
import lsst.pipe.base
Expand Down Expand Up @@ -433,88 +433,19 @@ def run(self, template, science, sources, visitSummary=None):
# checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
raise e

self.computeImageMetrics(science, subtractResults.difference, sources)
metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources)

return subtractResults

def computeImageMetrics(self, science, difference, stars):
r"""Compute quality metrics (saved to the task metadata) on the
difference image, at the locations of detected stars on the science
image. This restricts the metric to locations that should be
well-subtracted.

Parameters
----------
science : `lsst.afw.image.ExposureF`
Science exposure that was subtracted.
difference : `lsst.afw.image.ExposureF`
Result of subtracting template and science.
stars : `lsst.afw.table.SourceCatalog`
Good calibration sources detected on science image; these
footprints are what the metrics are computed on.

Notes
-----
The task metadata does not include docstrings, so descriptions of the
computed metrics are given here:

differenceFootprintRatioMean
Mean of the ratio of the absolute value of the difference image
(with the mean absolute value of the sky regions on the difference
image removed) to the science image, computed in the footprints
of stars detected on the science image (the sums below are of the
pixels in each star or sky footprint):
:math:`\mathrm{mean}_{footprints}((\sum |difference| -
\mathrm{mean}(\sum |difference_{sky}|)) / \sum science)`
differenceFootprintRatioStdev
Standard Deviation across footprints of the above ratio.
differenceFootprintSkyRatioMean
Mean of the ratio of the absolute value of sky source regions on
the difference image to the science image (the sum below is of the
pixels in each sky source footprint):
:math:`\mathrm{mean}_{footprints}(\sum |difference_{sky}| / \sum science_{sky})`
differenceFootprintSkyRatioStdev
Standard Deivation across footprints of the above sky ratio.
"""
def footprint_mean(sources, sky=0):
"""Compute ratio of the absolute value of the diffim to the science
image, within each source footprint, subtracting the sky from the
diffim values if provided.
"""
n = len(sources)
science_footprints = np.zeros(n)
difference_footprints = np.zeros(n)
ratio = np.zeros(n)
for i, record in enumerate(sources):
footprint = record.getFootprint()
heavy = lsst.afw.detection.makeHeavyFootprint(footprint, science.maskedImage)
heavy_diff = lsst.afw.detection.makeHeavyFootprint(footprint, difference.maskedImage)
science_footprints[i] = abs(heavy.getImageArray()).mean()
difference_footprints[i] = abs(heavy_diff.getImageArray()).mean()
ratio[i] = abs((difference_footprints[i] - sky)/science_footprints[i])
return science_footprints, difference_footprints, ratio

sky = stars["sky_source"]
if np.count_nonzero(sky) > 0:
sky_science, sky_difference, sky_ratio = footprint_mean(stars[sky])
sky_mean = sky_ratio.mean()
sky_std = sky_ratio.std()
sky_difference = sky_difference.mean()
else:
sky_mean = np.nan
sky_std = np.nan
sky_difference = 0
science_footprints, difference_footprints, ratio = footprint_mean(stars[~sky], sky_difference)

self.metadata["differenceFootprintRatioMean"] = ratio.mean()
self.metadata["differenceFootprintRatioStdev"] = ratio.std()
self.metadata["differenceFootprintSkyRatioMean"] = sky_mean
self.metadata["differenceFootprintSkyRatioStdev"] = sky_std
self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean
self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
self.log.info("Mean, stdev of ratio of difference to science "
"pixels in star footprints: %5.4f, %5.4f",
self.metadata["differenceFootprintRatioMean"],
self.metadata["differenceFootprintRatioStdev"])

return subtractResults

def runConvolveTemplate(self, template, science, selectSources):
"""Convolve the template image with a PSF-matching kernel and subtract
from the science image.
Expand Down
Loading