DM-39976: Replace masked image pixels with median values before convolution

python/lsst/ip/diffim/subtractImages.py

@@ -372,6 +372,8 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None,
                                                  )
         self.log.info("Science PSF FWHM: %f pixels", sciencePsfSize)
         self.log.info("Template PSF FWHM: %f pixels", templatePsfSize)
+        self.metadata.add("sciencePsfSize", sciencePsfSize)
+        self.metadata.add("templatePsfSize", templatePsfSize)
         selectSources = self._sourceSelector(sources, science.mask)
 
         if self.config.mode == "auto":
@@ -397,8 +399,10 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None,
             raise RuntimeError("Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
 
         if convolveTemplate:
+            self.metadata.add("convolvedExposure", "Template")
             subtractResults = self.runConvolveTemplate(template, science, selectSources)
         else:
+            self.metadata.add("convolvedExposure", "Science")
             subtractResults = self.runConvolveScience(template, science, selectSources)
 
         return subtractResults
@@ -608,11 +612,12 @@ def _validateExposures(template, science):
         assert templateBBox.contains(scienceBBox),\
             "Template bbox does not contain all of the science image."
 
-    @staticmethod
-    def _convolveExposure(exposure, kernel, convolutionControl,
+    def _convolveExposure(self, exposure, kernel, convolutionControl,
                           bbox=None,
                           psf=None,
-                          photoCalib=None):
+                          photoCalib=None,
+                          interpolateBadMaskPlanes=False,
+                          ):
         """Convolve an exposure with the given kernel.
 
         Parameters
@@ -640,8 +645,12 @@ def _convolveExposure(exposure, kernel, convolutionControl,
             convolvedExposure.setPsf(psf)
         if photoCalib is not None:
             convolvedExposure.setPhotoCalib(photoCalib)
-        convolvedImage = lsst.afw.image.MaskedImageF(exposure.getBBox())
-        lsst.afw.math.convolve(convolvedImage, exposure.maskedImage, kernel, convolutionControl)
+        if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None:
+            nInterp = _interpolateImage(convolvedExposure.maskedImage,
+                                        self.config.badMaskPlanes)
+            self.metadata.add("nInterpolated", nInterp)
+        convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox())
+        lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)
         convolvedExposure.setMaskedImage(convolvedImage)
         if bbox is None:
             return convolvedExposure
@@ -688,6 +697,7 @@ def _sourceSelector(self, sources, mask):
                            "%i selected but %i needed for the calculation.",
                            len(selectSources), self.config.makeKernel.nStarPerCell)
             raise RuntimeError("Cannot compute PSF matching kernel: too few sources selected.")
+        self.metadata.add("nPsfSources", len(selectSources))
 
         return selectSources.copy(deep=True)
 
@@ -849,8 +859,10 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None, visitS
 
         # TODO: DM-37212 we need to mirror the kernel in order to get correct cross correlation
         scienceKernel = science.psf.getKernel()
-        matchedScience = self._convolveExposure(science, scienceKernel, self.convolutionControl)
+        matchedScience = self._convolveExposure(science, scienceKernel, self.convolutionControl,
+                                                interpolateBadMaskPlanes=True)
         selectSources = self._sourceSelector(sources, matchedScience.mask)
+        self.metadata.add("convolvedExposure", "Preconvolution")
 
         subtractResults = self.runPreconvolve(template, science, matchedScience, selectSources, scienceKernel)
 
@@ -910,6 +922,7 @@ def runPreconvolve(self, template, science, matchedScience, selectSources, preCo
                                                  self.convolutionControl,
                                                  bbox=bbox,
                                                  psf=science.psf,
+                                                 interpolateBadMaskPlanes=True,
                                                  photoCalib=science.photoCalib)
         score = _subtractImages(matchedScience, matchedTemplate,
                                 backgroundModel=(kernelResult.backgroundModel
@@ -1025,3 +1038,31 @@ def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
     xTest = shape1[0] <= shape2[0]
     yTest = shape1[1] <= shape2[1]
     return xTest | yTest
+
+
+def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
+    """Replace masked image pixels with interpolated values.
+
+    Parameters
+    ----------
+    maskedImage : `lsst.afw.image.MaskedImage`
+        Image on which to perform interpolation.
+    badMaskPlanes : `list` of `str`
+        List of mask planes to interpolate over.
+    fallbackValue : `float`, optional
+        Value to set when interpolation fails.
+
+    Returns
+    -------
+    result: `float`
+        The number of masked pixels that were replaced.
+    """
+    image = maskedImage.image.array
+    badPixels = (maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(badMaskPlanes)) > 0
+    image[badPixels] = np.nan
+    if fallbackValue is None:
+        fallbackValue = np.nanmedian(image)
+    # For this initial implementation, skip the interpolation and just fill with
+    # the median value.
+    image[badPixels] = fallbackValue
+    return np.sum(badPixels)