Support missing STREAK and INJECTED mask planes

python/lsst/ip/diffim/detectAndMeasure.py

@@ -420,6 +420,9 @@ def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
             Warped and PSF-matched template that was used produce the
             difference image.
         """
+        # Ensure that the required mask planes are present
+        for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
+            difference.mask.addMaskPlane(mp)
         # Note that this may not be correct if we convolved the science image.
         # In the future we may wish to persist the matchedScience image.
         self.measurement.run(diaSources, difference, science, matchedTemplate)

python/lsst/ip/diffim/utils.py

@@ -1018,7 +1018,7 @@ def computeAveragePsf(exposure: afwImage.Exposure,
     return psfImage
 
 
-def detectTestSources(exposure):
+def detectTestSources(exposure, addMaskPlanes=None):
     """Minimal source detection wrapper suitable for unit tests.
 
     Parameters
@@ -1032,6 +1032,10 @@ def detectTestSources(exposure):
     selectSources :
         Source catalog containing candidates
     """
+    if addMaskPlanes is None:
+        # add empty streak mask plane in lieu of maskStreaksTask
+        # And add empty INJECTED and INJECTED_TEMPLATE mask planes
+        addMaskPlanes = ["STREAK", "INJECTED", "INJECTED_TEMPLATE"]
 
     schema = afwTable.SourceTable.makeMinimalSchema()
     selectDetection = measAlg.SourceDetectionTask(schema=schema)
@@ -1044,9 +1048,8 @@ def detectTestSources(exposure):
         sigma=None,  # The appropriate sigma is calculated from the PSF
         doSmooth=True
     )
-    exposure.mask.addMaskPlane("STREAK")  # add empty streak mask plane in lieu of maskStreaksTask
-    exposure.mask.addMaskPlane("INJECTED")  # add empty injected mask plane
-    exposure.mask.addMaskPlane("INJECTED_TEMPLATE")  # add empty injected template mask plane
+    for mp in addMaskPlanes:
+        exposure.mask.addMaskPlane(mp)
 
     selectSources = detRet.sources
     selectMeasurement.run(measCat=selectSources, exposure=exposure)
@@ -1078,6 +1081,7 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
                   yLoc=None,
                   flux=None,
                   clearEdgeMask=False,
+                  addMaskPlanes=None,
                   ):
     """Make a reproduceable PSF-convolved exposure for testing.
 
@@ -1119,6 +1123,8 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
         If specified, must have length equal to ``nSrc``
     clearEdgeMask : `bool`, optional
         Clear the "EDGE" mask plane after source detection.
+    addMaskPlanes : `list` of `str`, optional
+        Mask plane names to add to the image.
 
     Returns
     -------
@@ -1172,7 +1178,7 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
     modelExposure.image.array += noise
 
     # Run source detection to set up the mask plane
-    sourceCat = detectTestSources(modelExposure)
+    sourceCat = detectTestSources(modelExposure, addMaskPlanes=addMaskPlanes)
     if clearEdgeMask:
         modelExposure.mask &= ~modelExposure.mask.getPlaneBitMask("EDGE")
     modelExposure.setPhotoCalib(afwImage.PhotoCalib(calibration, 0., bbox))

tests/test_detectAndMeasure.py

@@ -282,6 +282,24 @@ def _detection_wrapper(positive=True):
         _detection_wrapper(positive=True)
         _detection_wrapper(positive=False)
 
+    def test_missing_mask_planes(self):
+        """Check that detection runs with missing mask planes.
+        """
+        # Set up the simulated images
+        noiseLevel = 1.
+        fluxLevel = 500
+        kwargs = {"psfSize": 2.4, "fluxLevel": fluxLevel, "addMaskPlanes": []}
+        # Use different seeds for the science and template so every source is a diaSource
+        science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
+        matchedTemplate, _ = makeTestImage(seed=6, noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
+
+        difference = science.clone()
+        difference.maskedImage -= matchedTemplate.maskedImage
+        detectionTask = self._setup_detection()
+
+        # Verify that detection runs without errors
+        detectionTask.run(science, matchedTemplate, difference)
+
     def test_detect_dipoles(self):
         """Run detection on a difference image containing dipoles.
         """

tests/test_subtractTask.py

@@ -173,6 +173,32 @@ def test_equal_images(self):
                                                 makeStats(), statistic=afwMath.STDEV)
         self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
 
+    def test_equal_images_missing_mask_planes(self):
+        """Test that running with enough sources produces reasonable output,
+        with the same size psf in the template and science and with missing
+        mask planes.
+        """
+        noiseLevel = 1.
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6, addMaskPlanes=[])
+        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+                                    templateBorderSize=20, doApplyCalibration=True, addMaskPlanes=[])
+        config = subtractImages.AlardLuptonSubtractTask.ConfigClass()
+        config.doSubtractBackground = False
+        task = subtractImages.AlardLuptonSubtractTask(config=config)
+        output = task.run(template, science, sources)
+        # There shoud be no NaN values in the difference image
+        self.assertTrue(np.all(np.isfinite(output.difference.image.array)))
+        # Mean of difference image should be close to zero.
+        meanError = noiseLevel/np.sqrt(output.difference.image.array.size)
+        # Make sure to include pixels with the DETECTED mask bit set.
+        statsCtrl = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA", "DETECTED", "DETECTED_NEGATIVE"))
+        differenceMean = computeRobustStatistics(output.difference.image, output.difference.mask, statsCtrl)
+        self.assertFloatsAlmostEqual(differenceMean, 0, atol=5*meanError)
+        # stddev of difference image should be close to expected value.
+        differenceStd = computeRobustStatistics(output.difference.image, output.difference.mask,
+                                                makeStats(), statistic=afwMath.STDEV)
+        self.assertFloatsAlmostEqual(differenceStd, np.sqrt(2)*noiseLevel, rtol=0.1)
+
     def test_psf_size(self):
         """Test that the image subtract task runs without failing, if
         fwhmExposureBuffer and fwhmExposureGrid parameters are set.