From 33d8c395dc18ca2fb2414299f1ee74cb34503d39 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 21 May 2025 15:03:42 -0700 Subject: [PATCH 1/2] Fix warnings in the diffim unit tests --- python/lsst/ip/diffim/subtractImages.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 05d1e644a..bad11777c 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -528,13 +528,21 @@ def footprint_mean(sources, sky=0): return science_footprints, difference_footprints, ratio sky = stars["sky_source"] - sky_science, sky_difference, sky_ratio = footprint_mean(stars[sky]) - science_footprints, difference_footprints, ratio = footprint_mean(stars[~sky], sky_difference.mean()) + 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_ratio.mean() - self.metadata["differenceFootprintSkyRatioStdev"] = sky_ratio.std() + self.metadata["differenceFootprintSkyRatioMean"] = sky_mean + self.metadata["differenceFootprintSkyRatioStdev"] = sky_std self.log.info("Mean, stdev of ratio of difference to science " "pixels in star footprints: %5.4f, %5.4f", self.metadata["differenceFootprintRatioMean"], From 1325264e6d96843d840fc8fcc70726026f398189 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 21 May 2025 15:03:54 -0700 Subject: [PATCH 2/2] Add residual QA metrics to detectAndMeasure --- python/lsst/ip/diffim/detectAndMeasure.py | 86 +++++++++++++++++++++-- 1 file changed, 81 insertions(+), 5 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 9d059a4b5..1bc909ed9 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -119,7 +119,8 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, doForcedMeasurement = pexConfig.Field( dtype=bool, default=True, - doc="Force photometer diaSource locations on PVI?") + doc="Force photometer diaSource locations on PVI?" + ) doAddMetrics = pexConfig.Field( dtype=bool, default=False, @@ -361,8 +362,9 @@ def __init__(self, **kwargs): self.schema.addField("refMatchId", "L", "unique id of reference catalog match") self.schema.addField("srcMatchId", "L", "unique id of source match") - if self.config.doSkySources: - self.makeSubtask("skySources", schema=self.schema) + # Create the sky source task for use by metrics, + # even if sky sources are not added to the diaSource catalog + self.makeSubtask("skySources", schema=self.schema) if self.config.doMaskStreaks: self.makeSubtask("maskStreaks") self.makeSubtask("streakDetection") @@ -597,7 +599,7 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor if self.config.doForcedMeasurement: self.measureForcedSources(diaSources, science, difference.getWcs()) - self.calculateMetrics(difference) + self.calculateMetrics(science, difference, diaSources) measurementResults = pipeBase.Struct( subtractedMeasuredExposure=difference, @@ -787,7 +789,7 @@ def measureForcedSources(self, diaSources, science, wcs): for diaSource, forcedSource in zip(diaSources, forcedSources): diaSource.assign(forcedSource, mapper) - def calculateMetrics(self, difference): + def calculateMetrics(self, science, difference, diaSources): """Add difference image QA metrics to the Task metadata. This may be used to produce corresponding metrics (see @@ -795,8 +797,12 @@ def calculateMetrics(self, difference): Parameters ---------- + 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. """ mask = difference.mask badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0 @@ -819,6 +825,76 @@ def calculateMetrics(self, difference): 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 + 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) + def _runStreakMasking(self, difference): """Do streak masking and optionally save the resulting streak fit parameters in a catalog.