Merge pull request #352 from lsst/tickets/DM-46852

DM-46852: Make detectAndMeasure robust to invalid PSF sigma

python/lsst/ip/diffim/detectAndMeasure.py

@@ -189,6 +189,11 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
                  "base_PixelFlags_flag_saturatedCenterAll",
                  ),
     )
+    clearMaskPlanes = lsst.pex.config.ListField(
+        dtype=str,
+        doc="Mask planes to clear before running detection.",
+        default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"),
+    )
     idGenerator = DetectorVisitIdGeneratorConfig.make_field()
 
     def setDefaults(self):
@@ -339,13 +344,7 @@ def run(self, science, matchedTemplate, difference,
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
-        # Ensure that we start with an empty detection and deblended mask.
-        mask = difference.mask
-        clearMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"]
-        for mp in clearMaskPlanes:
-            if mp not in mask.getMaskPlaneDict():
-                mask.addMaskPlane(mp)
-        mask &= ~mask.getPlaneBitMask(clearMaskPlanes)
+        self._prepareInputs(difference)
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't
@@ -365,6 +364,31 @@ def run(self, science, matchedTemplate, difference,
                                    positiveFootprints=positives,
                                    negativeFootprints=negatives)
 
+    def _prepareInputs(self, difference):
+        """Ensure that we start with an empty detection and deblended mask.
+
+        Parameters
+        ----------
+        difference : `lsst.afw.image.ExposureF`
+            The difference image that will be used for detecting diaSources.
+            The mask plane will be modified in place.
+
+        Raises
+        ------
+        lsst.pipe.base.UpstreamFailureNoWorkFound
+            If the PSF is not usable for measurement.
+        """
+        # Check that we have a valid PSF now before we do more work
+        sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
+        if np.isnan(sigma):
+            raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.")
+        # Ensure that we start with an empty detection and deblended mask.
+        mask = difference.mask
+        for mp in self.config.clearMaskPlanes:
+            if mp not in mask.getMaskPlaneDict():
+                mask.addMaskPlane(mp)
+        mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)
+
     def processResults(self, science, matchedTemplate, difference, sources, idFactory,
                        positiveFootprints=None, negativeFootprints=None,):
         """Measure and process the results of source detection.
@@ -733,9 +757,7 @@ def run(self, science, matchedTemplate, difference, scoreExposure,
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
-        # Ensure that we start with an empty detection mask.
-        mask = scoreExposure.mask
-        mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
+        self._prepareInputs(scoreExposure)
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't

python/lsst/ip/diffim/imageDecorrelation.py

@@ -274,7 +274,10 @@ def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatching
             corr = self.computeDiffimCorrection(kArr, varianceMean, targetVarianceMean)
 
         correctedImage = self.computeCorrectedImage(corr.corrft, subtractedExposure.image.array)
-        correctedPsf = self.computeCorrectedDiffimPsf(corr.corrft, psfImg.array)
+        # TODO DM-47461: only decorrelate the PSF if it is calculated for the
+        # difference image. If it is taken directly from the science (or template)
+        # image the decorrelation calculation is incorrect.
+        # correctedPsf = self.computeCorrectedDiffimPsf(corr.corrft, psfImg.array)
 
         # The subtracted exposure variance plane is already correlated, we cannot propagate
         # it through another convolution; instead we need to use the uncorrelated originals
@@ -300,7 +303,8 @@ def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatching
             correctedVariance /= kSumSq
         subtractedExposure.image.array[...] = correctedImage  # Allow for numpy type casting
         subtractedExposure.variance.array[...] = correctedVariance
-        subtractedExposure.setPsf(correctedPsf)
+        # TODO DM-47461
+        # subtractedExposure.setPsf(correctedPsf)
 
         newVarMean = self.computeVarianceMean(subtractedExposure)
         self.log.info("Variance plane mean of corrected diffim: %.5e", newVarMean)

tests/test_detectAndMeasure.py

@@ -22,9 +22,12 @@
 import numpy as np
 import unittest
 
+import lsst.afw.image as afwImage
+import lsst.afw.math as afwMath
 import lsst.geom
 from lsst.ip.diffim import detectAndMeasure, subtractImages
-from lsst.pipe.base import InvalidQuantumError
+import lsst.meas.algorithms as measAlg
+from lsst.pipe.base import InvalidQuantumError, UpstreamFailureNoWorkFound
 import lsst.utils.tests
 
 from utils import makeTestImage
@@ -161,6 +164,34 @@ def test_detection_xy0(self):
         self.assertTrue(all(output.diaSources['id'] > 1000000000))
         self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)
 
+    def test_raise_bad_psf(self):
+        """Detection should raise if the PSF width is NaN
+        """
+        # Set up the simulated images
+        noiseLevel = 1.
+        staticSeed = 1
+        fluxLevel = 500
+        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
+        science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
+        matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
+        difference = science.clone()
+
+        psfShape = difference.getPsf().computeBBox(lsst.geom.Point2D(0, 0)).getDimensions()
+        psfcI = afwImage.ImageD(lsst.geom.Extent2I(psfShape[1], psfShape[0]))
+        psfNew = np.zeros(psfShape)
+        psfNew[0, :] = 1
+        psfcI.array = psfNew
+        psfcK = afwMath.FixedKernel(psfcI)
+        psf = measAlg.KernelPsf(psfcK)
+        difference.setPsf(psf)
+
+        # Configure the detection Task
+        detectionTask = self._setup_detection()
+
+        # Run detection and check the results
+        with self.assertRaises(UpstreamFailureNoWorkFound):
+            detectionTask.run(science, matchedTemplate, difference)
+
     def test_measurements_finite(self):
         """Measured fluxes and centroids should always be finite.
         """