Add a image-based diffim quality metric

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.

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