From 4c7ac346ae9c40ee154f92b61d930845b29a9c2f Mon Sep 17 00:00:00 2001
From: John Parejko <parejkoj@uw.edu>
Date: Fri, 13 Dec 2024 15:32:43 -0800
Subject: [PATCH] Add image-based diffim metric, with test

---
 python/lsst/ip/diffim/subtractImages.py | 50 ++++++++++++++++++++++++-
 tests/test_subtractTask.py              | 16 ++++++++
 2 files changed, 65 insertions(+), 1 deletion(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 8cf0f76b..89d3f318 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -370,7 +370,8 @@ 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
         """
-
+        # Save this mask for use in metrics calculations, before we clear it.
+        science_mask = science.mask.clone()
         self._prepareInputs(template, science, visitSummary=visitSummary)
 
         # In the event that getPsfFwhm fails, evaluate the PSF on a grid.
@@ -465,8 +466,55 @@ def run(self, template, science, sources, visitSummary=None):
             #  checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
             raise e
 
+        self._compute_image_metrics(science, science_mask, subtractResults.difference, sources)
+
         return subtractResults
 
+    def _compute_image_metrics(self, science, science_mask, difference, stars):
+        """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.
+        science_mask : `lsst.afw.image.Mask`
+            Original, unmodified science mask.
+        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.
+        """
+        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()
+
     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