From 4761f75892b203d602a6d02d0a4bf673f4d1b113 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Fri, 13 Dec 2024 15:32:43 -0800 Subject: [PATCH] Add image-based diffim metric, with test --- python/lsst/ip/diffim/subtractImages.py | 73 ++++++++++++++++++++++++- tests/test_subtractTask.py | 16 ++++++ 2 files changed, 88 insertions(+), 1 deletion(-) diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py index 8cf0f76b..e02c38b6 100644 --- a/python/lsst/ip/diffim/subtractImages.py +++ b/python/lsst/ip/diffim/subtractImages.py @@ -370,7 +370,6 @@ def run(self, template, science, sources, visitSummary=None): Raised if fraction of good pixels, defined as not having NO_DATA set, is less then the configured requiredTemplateFraction """ - self._prepareInputs(template, science, visitSummary=visitSummary) # In the event that getPsfFwhm fails, evaluate the PSF on a grid. @@ -465,8 +464,80 @@ 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) + 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] = heavy.getImageArray().sum() + difference_footprints[i] = abs(heavy_diff.getImageArray()).sum() + ratio[i] = (difference_footprints[i] - sky) / science_footprints[i] + 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()) + + self.metadata["differenceFootprintRatioMean"] = ratio.mean() + self.metadata["differenceFootprintRatioStdev"] = ratio.std() + self.metadata["differenceFootprintSkyRatioMean"] = sky_ratio.mean() + self.metadata["differenceFootprintSkyRatioStdev"] = sky_ratio.std() + 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"]) + def runConvolveTemplate(self, template, science, selectSources): """Convolve the template image with a PSF-matching kernel and subtract from the science image. diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py index 930d5d52..17382f11 100644 --- a/tests/test_subtractTask.py +++ b/tests/test_subtractTask.py @@ -895,6 +895,18 @@ def test_metadata_metrics(self): template_good, _ = makeTestImage(psfSize=2.4, doApplyCalibration=True) template_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True) + # Add a few sky objects; sky footprints are needed for some metrics. + config = measAlg.SkyObjectsTask.ConfigClass() + config.nSources = 3 + skyTask = measAlg.SkyObjectsTask(config=config, name="skySources") + skyTask.skySourceKey = sources.schema["sky_source"].asKey() + skyTask.run(science.mask, 10, catalog=sources) + sources = sources.copy(deep=True) + # Add centroids, since these sources were added post-measurement. + for record in sources[sources["sky_source"]]: + record["truth_x"] = record.getFootprint().getPeaks()[0].getFx() + record["truth_y"] = record.getFootprint().getPeaks()[0].getFy() + # The metadata fields are attached to the subtractTask, so we do # need to run that; run it for both "good" and "bad" seeing templates @@ -955,6 +967,10 @@ def test_metadata_metrics(self): self.assertIn('scienceLimitingMagnitude', subtractTask_good.metadata) self.assertIn('templateLimitingMagnitude', subtractTask_good.metadata) + # The mean ratio metric should be much worse on the "bad" subtraction. + self.assertLess(subtractTask_good.metadata['differenceFootprintRatioMean'], 0.2) + self.assertGreater(subtractTask_bad.metadata['differenceFootprintRatioMean'], 1.0) + class AlardLuptonPreconvolveSubtractTest(AlardLuptonSubtractTestBase, lsst.utils.tests.TestCase): subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask