DM-52420: Add flag for diaSources in high variance regions

python/lsst/ip/diffim/detectAndMeasure.py

@@ -416,9 +416,9 @@ def setDefaults(self):
 
         # Keep track of which footprints contain streaks
         self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [
-            "STREAK", "INJECTED", "INJECTED_TEMPLATE"]
+            "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE"]
         self.measurement.plugins["base_PixelFlags"].masksFpCenter = [
-            "STREAK", "INJECTED", "INJECTED_TEMPLATE"]
+            "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE"]
         self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"]
 
     def validate(self):
@@ -808,7 +808,10 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
             # This option allows returning sky sources,
             # even if there are no diaSources
             measurementResults.diaSources = diaSources
-
+        self.log.info("Measured %d diaSources and %d sky sources",
+                      np.count_nonzero(~diaSources["sky_source"]),
+                      np.count_nonzero(diaSources["sky_source"])
+                      )
         return measurementResults
 
     def _deblend(self, difference, positiveFootprints, negativeFootprints):

python/lsst/ip/diffim/getTemplate.py

@@ -34,7 +34,7 @@
 
 from lsst.skymap import BaseSkyMap
 from lsst.ip.diffim.dcrModel import DcrModel
-from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig
+from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask
 from lsst.utils.timer import timeMethod
 
 __all__ = [
@@ -103,6 +103,17 @@ class GetTemplateConfig(
         doc="Configuration for CoaddPsf",
         dtype=CoaddPsfConfig,
     )
+    varianceBackground = pexConfig.ConfigurableField(
+        target=SubtractBackgroundTask,
+        doc="Task to estimate the background variance.",
+    )
+    highVarianceThreshold = pexConfig.RangeField(
+        dtype=float,
+        default=4,
+        min=1,
+        doc="Set the HIGH_VARIANCE mask plane for regions with variance"
+        " greater than the median by this factor.",
+    )
 
     def setDefaults(self):
         # Use a smaller cache: per SeparableKernel.computeCache, this should
@@ -115,6 +126,19 @@ def setDefaults(self):
         self.warp.warpingKernelName = "lanczos3"
         self.coaddPsf.warpingKernelName = self.warp.warpingKernelName
 
+        # Background subtraction of the variance plane
+        self.varianceBackground.algorithm = "LINEAR"
+        self.varianceBackground.binSize = 32
+        self.varianceBackground.useApprox = False
+        self.varianceBackground.statisticsProperty = "MEDIAN"
+        self.varianceBackground.doFilterSuperPixels = True
+        self.varianceBackground.ignoredPixelMask = ["BAD",
+                                                    "EDGE",
+                                                    "DETECTED",
+                                                    "DETECTED_NEGATIVE",
+                                                    "NO_DATA",
+                                                    ]
+
 
 class GetTemplateTask(pipeBase.PipelineTask):
     ConfigClass = GetTemplateConfig
@@ -137,6 +161,7 @@ def __init__(self, *args, **kwargs):
             type=float,
             doc="Weight for each exposure, used to make the CoaddPsf; should always be 1.",
         )
+        self.makeSubtask("varianceBackground")
 
     def runQuantum(self, butlerQC, inputRefs, outputRefs):
         inputs = butlerQC.get(inputRefs)
@@ -357,11 +382,30 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter):
         for c in catalogs:
             catalog.extend(c)
 
+        # Set a mask plane for any regions with exceptionally high variance.
+        self.checkHighVariance(template)
         template.setPsf(self._makePsf(template, catalog, wcs))
         template.setFilter(afwImage.FilterLabel(band, physical_filter))
         template.setPhotoCalib(photoCalib)
         return pipeBase.Struct(template=template)
 
+    def checkHighVariance(self, template):
+        """Set a mask plane for regions with unusually high variance.
+
+        Parameters
+        ----------
+        template : `lsst.afw.image.Exposure`
+            The warped template exposure, which will be modified in place.
+        """
+        highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE")
+        template.mask.getPlaneBitMask("HIGH_VARIANCE")
+        varianceExposure = template.clone()
+        varianceExposure.image.array = varianceExposure.variance.array
+        varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array
+        threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground)
+        highVariancePix = varianceBackground > threshold
+        template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit
+
     @staticmethod
     def _checkInputs(dataIds, coaddExposures):
         """Check that the all the dataIds are from the same band and that

python/lsst/ip/diffim/subtractImages.py

@@ -26,7 +26,7 @@
 import lsst.afw.image
 import lsst.afw.math
 import lsst.geom
-from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics
+from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics, checkMask
 from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
 import lsst.pex.config
 import lsst.pipe.base
@@ -268,9 +268,8 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
     )
     excludeMaskPlanes = lsst.pex.config.ListField(
         dtype=str,
-        default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE"),
-        doc="Mask planes to exclude when selecting sources for PSF matching.",
-        deprecated="No longer used. Will be removed after v30"
+        default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"),
+        doc="Template mask planes to exclude when selecting sources for PSF matching.",
     )
     badMaskPlanes = lsst.pex.config.ListField(
         dtype=str,
@@ -279,7 +278,7 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
     )
     preserveTemplateMask = lsst.pex.config.ListField(
         dtype=str,
-        default=("NO_DATA", "BAD",),
+        default=("NO_DATA", "BAD", "HIGH_VARIANCE"),
         doc="Mask planes from the template to propagate to the image difference."
     )
     renameTemplateMask = lsst.pex.config.ListField(
@@ -443,7 +442,7 @@ def run(self, template, science, sources, visitSummary=None):
 
         convolveTemplate = self.chooseConvolutionMethod(template, science)
 
-        selectSources = self._sourceSelector(sources, science.getBBox())
+        selectSources = self._sourceSelector(sources, science.getBBox(), template.mask)
 
         kernelResult = self.runMakeKernel(template, science, selectSources,
                                           convolveTemplate=convolveTemplate)
@@ -627,7 +626,7 @@ def runKernelSourceDetection(self, template, science):
                                                       scienceFwhmPix=self.sciencePsfSize)
 
         # return sources
-        return self._sourceSelector(sources, science.getBBox(), fallback=True)
+        return self._sourceSelector(sources, science.getBBox(), template.mask, fallback=True)
 
     def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None):
         """Convolve the template image with a PSF-matching kernel and subtract
@@ -908,7 +907,7 @@ def _convolveExposure(self, exposure, kernel, convolutionControl,
         else:
             return convolvedExposure[bbox]
 
-    def _sourceSelector(self, sources, bbox, fallback=False):
+    def _sourceSelector(self, sources, bbox, mask, fallback=False):
         """Select sources from a catalog that meet the selection criteria.
         The selection criteria include any configured parameters of the
         `sourceSelector` subtask, as well as distance from the edge if
@@ -920,6 +919,13 @@ def _sourceSelector(self, sources, bbox, fallback=False):
             Input source catalog to select sources from.
         bbox : `lsst.geom.Box2I`
             Bounding box of the science image.
+        mask : `lsst.afw.image.Mask`
+            The mask plane of the template to use to reject kernel candidates.
+        fallback : `bool`, optional
+            Switch indicating the source selector is being called after
+            running the fallback source detection subtask, which does not run a
+            full set of measurement plugins and can't use the same settings for
+            the source selector.
 
         Returns
         -------
@@ -929,6 +935,9 @@ def _sourceSelector(self, sources, bbox, fallback=False):
 
         Raises
         ------
+        InsufficientKernelSourcesError
+            An AlgorithmError that is raised if there are not enough PSF
+            candidates to construct the PSF matching kernel.
         RuntimeError
             If there are too few sources to compute the PSF matching kernel
             remaining after source selection.
@@ -945,6 +954,9 @@ def _sourceSelector(self, sources, bbox, fallback=False):
                           np.count_nonzero(~bboxSelected), rejectRadius)
             selected &= bboxSelected
         selectSources = sources[selected].copy(deep=True)
+        # Optionally remove sources that land on masked pixels
+        maskSelected = checkMask(mask, selectSources, self.config.excludeMaskPlanes, checkAdjacent=True)
+        selectSources = selectSources[maskSelected].copy(deep=True)
         # Trim selectSources if they exceed ``maxKernelSources``.
         # Keep the highest signal-to-noise sources of those selected.
         if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
@@ -1211,7 +1223,7 @@ def run(self, template, science, sources, visitSummary=None):
                                                 interpolateBadMaskPlanes=True)
         self.metadata["convolvedExposure"] = "Preconvolution"
         try:
-            selectSources = self._sourceSelector(sources, science.getBBox())
+            selectSources = self._sourceSelector(sources, science.getBBox(), template.mask)
             subtractResults = self.runPreconvolve(template, science, matchedScience,
                                                   selectSources, scienceKernel)
 

python/lsst/ip/diffim/utils.py

@@ -23,7 +23,7 @@
 
 
 __all__ = ["evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement",
-           "computeDifferenceImageMetrics",
+           "computeDifferenceImageMetrics", "checkMask"
            ]
 
 import itertools
@@ -454,3 +454,44 @@ def populate_sattle_visit_cache(visit_info, historical=False):
               "historical": historical})
 
     r.raise_for_status()
+
+
+def checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
+    """Exclude sources that are located on or adjacent to masked pixels.
+
+    Parameters
+    ----------
+    mask : `lsst.afw.image.Mask`
+        The image mask plane to use to reject sources
+        based on the location of their centroid on the ccd.
+    sources : `lsst.afw.table.SourceCatalog`
+        The source catalog to evaluate.
+    excludeMaskPlanes : `list` of `str`
+        List of the names of the mask planes to exclude.
+
+    Returns
+    -------
+    flags : `numpy.ndarray` of `bool`
+        Array indicating whether each source in the catalog should be
+        kept (True) or rejected (False) based on the value of the
+        mask plane at its location.
+    """
+    setExcludeMaskPlanes = [
+        maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict()
+    ]
+
+    excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
+
+    xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
+    yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
+
+    flags = np.ones(len(sources), dtype=bool)
+    if checkAdjacent:
+        pixRange = (0, -1, 1)
+    else:
+        pixRange = (0,)
+    for j in pixRange:
+        for i in pixRange:
+            mv = mask.array[yv + j, xv + i]
+            flags *= np.bitwise_and(mv, excludePixelMask) == 0
+    return flags

tests/test_subtractTask.py

@@ -117,9 +117,9 @@ def test_incomplete_template_coverage(self):
         border = 20
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6, nSrc=50,
+        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6, nSrc=100,
                                          xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, nSrc=50,
+        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, nSrc=100,
                                     templateBorderSize=border, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
 
@@ -485,7 +485,7 @@ def _run_and_check_sources(sourcesIn, maxKernelSources=1000, minKernelSources=3,
             signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr()
             signalToNoise = signalToNoise[~sources[badSourceFlag]]
             signalToNoise.sort()
-            selectSources = task._sourceSelector(sources, science.getBBox())
+            selectSources = task._sourceSelector(sources, science.getBBox(), template.mask)
             self.assertEqual(nGoodSources, len(selectSources))
             signalToNoiseOut = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
             signalToNoiseOut.sort()