Merge pull request #301 from lsst/tickets/DM-42774

DM-42774: Rename difference image mask planes to include _TEMPLATE if they are set from the template

python/lsst/ip/diffim/detectAndMeasure.py

@@ -304,21 +304,21 @@ def processResults(self, science, matchedTemplate, difference, sources, table,
 
         Parameters
         ----------
-        sources : `lsst.afw.table.SourceCatalog`
-            Detected sources on the difference exposure.
-        positiveFootprints : `lsst.afw.detection.FootprintSet`, optional
-            Positive polarity footprints.
-        negativeFootprints : `lsst.afw.detection.FootprintSet`, optional
-            Negative polarity footprints.
-        table : `lsst.afw.table.SourceTable`
-            Table object that will be used to create the SourceCatalog.
         science : `lsst.afw.image.ExposureF`
             Science exposure that the template was subtracted from.
         matchedTemplate : `lsst.afw.image.ExposureF`
             Warped and PSF-matched template that was used produce the
             difference image.
         difference : `lsst.afw.image.ExposureF`
             Result of subtracting template from the science image.
+        sources : `lsst.afw.table.SourceCatalog`
+            Detected sources on the difference exposure.
+        table : `lsst.afw.table.SourceTable`
+            Table object that will be used to create the SourceCatalog.
+        positiveFootprints : `lsst.afw.detection.FootprintSet`, optional
+            Positive polarity footprints.
+        negativeFootprints : `lsst.afw.detection.FootprintSet`, optional
+            Negative polarity footprints.
 
         Returns
         -------
@@ -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/subtractImages.py

@@ -233,9 +233,15 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
     )
     preserveTemplateMask = lsst.pex.config.ListField(
         dtype=str,
-        default=("NO_DATA", "BAD", "SAT", "FAKE", "INJECTED", "INJECTED_CORE"),
+        default=("NO_DATA", "BAD",),
         doc="Mask planes from the template to propagate to the image difference."
     )
+    renameTemplateMask = lsst.pex.config.ListField(
+        dtype=str,
+        default=("SAT", "INJECTED", "INJECTED_CORE",),
+        doc="Mask planes from the template to propagate to the image difference"
+        "with '_TEMPLATE' appended to the name."
+    )
     allowKernelSourceDetection = lsst.pex.config.Field(
         dtype=bool,
         default=False,
@@ -624,11 +630,6 @@ def finalize(self, template, science, difference, kernel,
         correctedExposure : `lsst.afw.image.ExposureF`
             The decorrelated image difference.
         """
-        # Erase existing detection mask planes.
-        #  We don't want the detection mask from the science image
-
-        self.updateMasks(template, science, difference)
-
         if self.config.doDecorrelation:
             self.log.info("Decorrelating image difference.")
             # We have cleared the template mask plane, so copy the mask plane of
@@ -644,35 +645,6 @@ def finalize(self, template, science, difference, kernel,
             correctedExposure = difference
         return correctedExposure
 
-    def updateMasks(self, template, science, difference):
-        """Update the mask planes on images for finalizing."""
-
-        bbox = science.getBBox()
-        mask = difference.mask
-        mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
-
-        self.log.info("Adding injected mask planes")
-        mask.addMaskPlane("INJECTED")
-        mask.addMaskPlane("INJECTED_TEMPLATE")
-
-        if "FAKE" in science.mask.getMaskPlaneDict().keys():
-            # propagate the mask plane related to Fake source injection
-            # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED
-            # NOTE: This can be removed in DM-40796
-            diffInjectedBitMask = mask.getPlaneBitMask("INJECTED")
-            diffInjTmpltBitMask = mask.getPlaneBitMask("INJECTED_TEMPLATE")
-
-            scienceFakeBitMask = science.mask.getPlaneBitMask('FAKE')
-            tmpltFakeBitMask = template[bbox].mask.getPlaneBitMask('FAKE')
-
-            injScienceMaskArray = ((science.mask.array & scienceFakeBitMask) > 0) * diffInjectedBitMask
-            injTemplateMaskArray = ((template[bbox].mask.array & tmpltFakeBitMask) > 0) * diffInjTmpltBitMask
-
-            mask.array |= injScienceMaskArray
-            mask.array |= injTemplateMaskArray
-
-        template[bbox].mask.array[...] = difference.mask.array[...]
-
     @staticmethod
     def _validateExposures(template, science):
         """Check that the WCS of the two Exposures match, and the template bbox
@@ -868,20 +840,87 @@ def _prepareInputs(self, template, science, visitSummary=None):
             self.metadata.add("scaleTemplateVarianceFactor", templateVarFactor)
             self.log.info("Science variance scaling factor: %.2f", sciVarFactor)
             self.metadata.add("scaleScienceVarianceFactor", sciVarFactor)
-        self._clearMask(template)
 
-    def _clearMask(self, template):
-        """Clear the mask plane of the template.
+        # Erase existing detection mask planes.
+        #  We don't want the detection mask from the science image
+        self.updateMasks(template, science)
+
+    def updateMasks(self, template, science):
+        """Update the science and template mask planes before differencing.
 
         Parameters
         ----------
-        template : `lsst.afw.image.ExposureF`
+        template : `lsst.afw.image.Exposure`
             Template exposure, warped to match the science exposure.
-            The mask plane will be modified in place.
+            The template mask planes will be erased, except for a few specified
+            in the task config.
+        science : `lsst.afw.image.Exposure`
+            Science exposure to subtract from the template.
+            The DETECTED and DETECTED_NEGATIVE mask planes of the science image
+            will be erased.
+        """
+        self._clearMask(science.mask, clearMaskPlanes=["DETECTED", "DETECTED_NEGATIVE"])
+
+        # We will clear ALL template mask planes, except for those specified
+        # via the `preserveTemplateMask` config. Mask planes specified via
+        # the `renameTemplateMask` config will be copied to new planes with
+        # "_TEMPLATE" appended to their names, and the original mask plane will
+        # be cleared.
+        clearMaskPlanes = [mp for mp in template.mask.getMaskPlaneDict().keys()
+                           if mp not in self.config.preserveTemplateMask]
+        renameMaskPlanes = [mp for mp in self.config.renameTemplateMask
+                            if mp in template.mask.getMaskPlaneDict().keys()]
+
+        # propagate the mask plane related to Fake source injection
+        # NOTE: the fake source injection sets FAKE plane, but it should be INJECTED
+        # NOTE: This can be removed in DM-40796
+        if "FAKE" in science.mask.getMaskPlaneDict().keys():
+            self.log.info("Adding injected mask plane to science image")
+            self._renameMaskPlanes(science.mask, "FAKE", "INJECTED")
+        if "FAKE" in template.mask.getMaskPlaneDict().keys():
+            self.log.info("Adding injected mask plane to template image")
+            self._renameMaskPlanes(template.mask, "FAKE", "INJECTED_TEMPLATE")
+            if "INJECTED" in renameMaskPlanes:
+                renameMaskPlanes.remove("INJECTED")
+            if "INJECTED_TEMPLATE" in clearMaskPlanes:
+                clearMaskPlanes.remove("INJECTED_TEMPLATE")
+
+        for maskPlane in renameMaskPlanes:
+            self._renameMaskPlanes(template.mask, maskPlane, maskPlane + "_TEMPLATE")
+        self._clearMask(template.mask, clearMaskPlanes=clearMaskPlanes)
+
+    @staticmethod
+    def _renameMaskPlanes(mask, maskPlane, newMaskPlane):
+        """Rename a mask plane by adding the new name and copying the data.
+
+        Parameters
+        ----------
+        mask : `lsst.afw.image.Mask`
+            The mask image to update in place.
+        maskPlane : `str`
+            The name of the existing mask plane to copy.
+        newMaskPlane : `str`
+            The new name of the mask plane that will be added.
+            If the mask plane already exists, it will be updated in place.
+        """
+        mask.addMaskPlane(newMaskPlane)
+        originBitMask = mask.getPlaneBitMask(maskPlane)
+        destinationBitMask = mask.getPlaneBitMask(newMaskPlane)
+        mask.array |= ((mask.array & originBitMask) > 0)*destinationBitMask
+
+    def _clearMask(self, mask, clearMaskPlanes=None):
+        """Clear the mask plane of the template.
+
+        Parameters
+        ----------
+        mask : `lsst.afw.image.Mask`
+            The mask plane to erase, which will be modified in place.
+        clearMaskPlanes : `list` of `str`, optional
+            Erase the specified mask planes.
+            If not supplied, the entire mask will be erased.
         """
-        mask = template.mask
-        clearMaskPlanes = [maskplane for maskplane in mask.getMaskPlaneDict().keys()
-                           if maskplane not in self.config.preserveTemplateMask]
+        if clearMaskPlanes is None:
+            clearMaskPlanes = list(mask.getMaskPlaneDict().keys())
 
         bitMaskToClear = mask.getPlaneBitMask(clearMaskPlanes)
         mask &= ~bitMaskToClear

python/lsst/ip/diffim/utils.py

@@ -1018,20 +1018,26 @@ def computeAveragePsf(exposure: afwImage.Exposure,
     return psfImage
 
 
-def detectTestSources(exposure):
+def detectTestSources(exposure, addMaskPlanes=None):
     """Minimal source detection wrapper suitable for unit tests.
 
     Parameters
     ----------
     exposure : `lsst.afw.image.Exposure`
         Exposure on which to run detection/measurement
         The exposure is modified in place to set the 'DETECTED' mask plane.
+    addMaskPlanes : `list` of `str`, optional
+        Additional mask planes to add to the maskedImage of the exposure.
 
     Returns
     -------
-    selectSources :
+    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 +1050,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 +1083,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 +1125,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 +1180,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.
         """
@@ -481,6 +499,9 @@ def test_fake_mask_plane_propagation(self):
         injTmplt_masked = (diff_mask.array & diff_mask.getPlaneBitMask("INJECTED_TEMPLATE")) > 0
 
         self.assertFloatsEqual(inj_masked.astype(int), science_fake_masked.astype(int))
+        # The template is convolved, so the INJECTED_TEMPLATE mask plane may
+        # include more pixels than the FAKE mask plane
+        injTmplt_masked &= template_fake_masked
         self.assertFloatsEqual(injTmplt_masked.astype(int), template_fake_masked.astype(int))
 
         # Now check that detection of fakes have the correct flag for injections

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.