Merge pull request #282 from lsst/tickets/DM-41118

DM-41118: add diffim QA metrics

python/lsst/ip/diffim/detectAndMeasure.py

@@ -20,6 +20,7 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from deprecated.sphinx import deprecated
+import numpy as np
 
 import lsst.afw.table as afwTable
 import lsst.daf.base as dafBase
@@ -148,7 +149,7 @@ def setDefaults(self):
         self.detection.thresholdValue = 5.0
         self.detection.reEstimateBackground = False
         self.detection.thresholdType = "pixel_stdev"
-        self.detection.excludeMaskPlanes = ["EDGE"]
+        self.detection.excludeMaskPlanes = ["EDGE", "SAT", "BAD", "INTRP"]
 
         # Add filtered flux measurement, the correct measurement for pre-convolved images.
         self.measurement.algorithms.names.add('base_PeakLikelihoodFlux')
@@ -347,6 +348,7 @@ def processResults(self, science, matchedTemplate, difference, sources, table,
             ``diaSources``  : `lsst.afw.table.SourceCatalog`
                 The catalog of detected sources.
         """
+        self.metadata.add("nUnmergedDiaSources", len(sources))
         if self.config.doMerge:
             fpSet = positiveFootprints
             fpSet.merge(negativeFootprints, self.config.growFootprint,
@@ -356,6 +358,7 @@ def processResults(self, science, matchedTemplate, difference, sources, table,
             self.log.info("Merging detections into %d sources", len(diaSources))
         else:
             diaSources = sources
+        self.metadata.add("nMergedDiaSources", len(diaSources))
 
         if self.config.doSkySources:
             self.addSkySources(diaSources, difference.mask, difference.info.id)
@@ -369,6 +372,7 @@ def processResults(self, science, matchedTemplate, difference, sources, table,
             subtractedMeasuredExposure=difference,
             diaSources=diaSources,
         )
+        self.calculateMetrics(difference)
 
         return measurementResults
 
@@ -453,6 +457,27 @@ def measureForcedSources(self, diaSources, science, wcs):
         for diaSource, forcedSource in zip(diaSources, forcedSources):
             diaSource.assign(forcedSource, mapper)
 
+    def calculateMetrics(self, difference):
+        """Add image QA metrics to the Task metadata.
+
+        Parameters
+        ----------
+        difference : `lsst.afw.image.Exposure`
+            The target image to calculate metrics for.
+        """
+        mask = difference.mask
+        goodPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) == 0
+        self.metadata.add("nGoodPixels", np.sum(goodPix))
+        self.metadata.add("nBadPixels", np.sum(~goodPix))
+        detPosPix = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0
+        detNegPix = (mask.array & mask.getPlaneBitMask("DETECTED_NEGATIVE")) > 0
+        self.metadata.add("nPixelsDetectedPositive", np.sum(detPosPix))
+        self.metadata.add("nPixelsDetectedNegative", np.sum(detNegPix))
+        detPosPix *= ~goodPix
+        detNegPix *= ~goodPix
+        self.metadata.add("nBadPixelsDetectedPositive", np.sum(detPosPix))
+        self.metadata.add("nBadPixelsDetectedNegative", np.sum(detNegPix))
+
 
 class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
     scoreExposure = pipeBase.connectionTypes.Input(

python/lsst/ip/diffim/subtractImages.py

@@ -422,10 +422,11 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None,
         except (RuntimeError, lsst.pex.exceptions.Exception) as e:
             self.log.warn("Failed to match template. Checking coverage")
             #  Raise NoWorkFound if template fraction is insufficient
-            checkTemplateIsSufficient(template, self.log,
-                                      self.config.minTemplateFractionForExpectedSuccess,
-                                      exceptionMessage="Template coverage lower than expected to succeed."
-                                      f" Failure is tolerable: {e}")
+            self.checkTemplateIsSufficient(
+                template,
+                exceptionMessage="Template coverage lower than expected to succeed."
+                f" Failure is tolerable: {e}"
+            )
             #  checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
             raise e
 
@@ -815,10 +816,9 @@ def _prepareInputs(self, template, science, visitSummary=None):
         self._validateExposures(template, science)
         if visitSummary is not None:
             self._applyExternalCalibrations(science, visitSummary=visitSummary)
-        checkTemplateIsSufficient(template, self.log,
-                                  requiredTemplateFraction=self.config.requiredTemplateFraction,
-                                  exceptionMessage="Not attempting subtraction. To force subtraction,"
-                                  " set config requiredTemplateFraction=0")
+        self.checkTemplateIsSufficient(template,
+                                       exceptionMessage="Not attempting subtraction. To force subtraction,"
+                                       " set config requiredTemplateFraction=0")
 
         if self.config.doScaleVariance:
             # Scale the variance of the template and science images before
@@ -848,6 +848,39 @@ def _clearMask(self, template):
         bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
         mask &= ~bitMaskToClear
 
+    def checkTemplateIsSufficient(self, templateExposure, exceptionMessage=""):
+        """Raise NoWorkFound if template coverage < requiredTemplateFraction
+
+        Parameters
+        ----------
+        templateExposure : `lsst.afw.image.Exposure`
+            The template exposure to check
+        exceptionMessage : `str`, optional
+            Message to include in the exception raised if the template coverage
+            is insufficient.
+
+        Raises
+        ------
+        lsst.pipe.base.NoWorkFound
+            Raised if fraction of good pixels, defined as not having NO_DATA
+            set, is less than the configured requiredTemplateFraction
+        """
+        # Count the number of pixels with the NO_DATA mask bit set
+        # counting NaN pixels is insufficient because pixels without data are often intepolated over)
+        pixNoData = np.count_nonzero(templateExposure.mask.array
+                                     & templateExposure.mask.getPlaneBitMask('NO_DATA'))
+        pixGood = templateExposure.getBBox().getArea() - pixNoData
+        coverageFraction = pixGood/templateExposure.getBBox().getArea()
+        self.log.info("template has %d good pixels (%.1f%%)", pixGood,
+                      100*coverageFraction)
+        self.metadata.add("templateCoveragePercent", 100*coverageFraction)
+
+        if pixGood/templateExposure.getBBox().getArea() < self.config.requiredTemplateFraction:
+            message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
+                       100*coverageFraction,
+                       100*self.config.requiredTemplateFraction))
+            raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
+
 
 class AlardLuptonPreconvolveSubtractConnections(SubtractInputConnections,
                                                 SubtractScoreOutputConnections):
@@ -1007,44 +1040,6 @@ def runPreconvolve(self, template, science, matchedScience, selectSources, preCo
                                      psfMatchingKernel=kernelResult.psfMatchingKernel)
 
 
-def checkTemplateIsSufficient(templateExposure, logger, requiredTemplateFraction=0.,
-                              exceptionMessage=""):
-    """Raise NoWorkFound if template coverage < requiredTemplateFraction
-
-    Parameters
-    ----------
-    templateExposure : `lsst.afw.image.ExposureF`
-        The template exposure to check
-    logger : `lsst.log.Log`
-        Logger for printing output.
-    requiredTemplateFraction : `float`, optional
-        Fraction of pixels of the science image required to have coverage
-        in the template.
-    exceptionMessage : `str`, optional
-        Message to include in the exception raised if the template coverage
-        is insufficient.
-
-    Raises
-    ------
-    lsst.pipe.base.NoWorkFound
-        Raised if fraction of good pixels, defined as not having NO_DATA
-        set, is less than the requiredTemplateFraction
-    """
-    # Count the number of pixels with the NO_DATA mask bit set
-    # counting NaN pixels is insufficient because pixels without data are often intepolated over)
-    pixNoData = np.count_nonzero(templateExposure.mask.array
-                                 & templateExposure.mask.getPlaneBitMask('NO_DATA'))
-    pixGood = templateExposure.getBBox().getArea() - pixNoData
-    logger.info("template has %d good pixels (%.1f%%)", pixGood,
-                100*pixGood/templateExposure.getBBox().getArea())
-
-    if pixGood/templateExposure.getBBox().getArea() < requiredTemplateFraction:
-        message = ("Insufficient Template Coverage. (%.1f%% < %.1f%%)" % (
-                   100*pixGood/templateExposure.getBBox().getArea(),
-                   100*requiredTemplateFraction))
-        raise lsst.pipe.base.NoWorkFound(message + " " + exceptionMessage)
-
-
 def _subtractImages(science, template, backgroundModel=None):
     """Subtract template from science, propagating relevant metadata.