Merge pull request #340 from lsst/tickets/DM-46141

DM-46141: Use ScienceSourceSelector to select kernel candidates

python/lsst/ip/diffim/subtractImages.py

@@ -27,7 +27,7 @@
 import lsst.afw.math
 import lsst.geom
 from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm
-from lsst.meas.algorithms import ScaleVarianceTask
+from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
 import lsst.pex.config
 import lsst.pipe.base
 import lsst.pex.exceptions
@@ -173,6 +173,10 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
         dtype=bool,
         default=False,
     )
+    sourceSelector = lsst.pex.config.ConfigurableField(
+        target=ScienceSourceSelectorTask,
+        doc="Task to select sources to be used for PSF matching.",
+    )
     detectionThreshold = lsst.pex.config.Field(
         dtype=float,
         default=10,
@@ -206,6 +210,8 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
                  "base_PixelFlags_flag_saturated",
                  "base_PixelFlags_flag_bad",
                  ),
+        deprecated="This field is no longer used and will be removed after v28."
+        "Set the equivalent field for the sourceSelector subtask instead.",
     )
     excludeMaskPlanes = lsst.pex.config.ListField(
         dtype=str,
@@ -240,6 +246,13 @@ def setDefaults(self):
         self.makeKernel.kernel.active.fitForBackground = self.doSubtractBackground
         self.makeKernel.kernel.active.spatialKernelOrder = 1
         self.makeKernel.kernel.active.spatialBgOrder = 2
+        self.sourceSelector.doUnresolved = True  # apply star-galaxy separation
+        self.sourceSelector.doIsolated = True  # apply isolated star selection
+        self.sourceSelector.doRequirePrimary = True  # apply primary flag selection
+        self.sourceSelector.doSkySources = False  # Do not include sky sources
+        self.sourceSelector.doSignalToNoise = True  # apply signal to noise filter
+        self.sourceSelector.signalToNoise.minimum = 10
+        self.sourceSelector.signalToNoise.maximum = 500
 
 
 class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
@@ -265,6 +278,7 @@ def __init__(self, **kwargs):
         super().__init__(**kwargs)
         self.makeSubtask("decorrelate")
         self.makeSubtask("makeKernel")
+        self.makeSubtask("sourceSelector")
         if self.config.doScaleVariance:
             self.makeSubtask("scaleVariance")
 
@@ -772,26 +786,23 @@ def _sourceSelector(self, sources, mask):
             If there are too few sources to compute the PSF matching kernel
             remaining after source selection.
         """
-        flags = np.ones(len(sources), dtype=bool)
-        for flag in self.config.badSourceFlags:
-            try:
-                flags *= ~sources[flag]
-            except Exception as e:
-                self.log.warning("Could not apply source flag: %s", e)
-        signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr()
-        sToNFlag = signalToNoise > self.config.detectionThreshold
-        flags *= sToNFlag
-        sToNFlagMax = signalToNoise < self.config.detectionThresholdMax
-        flags *= sToNFlagMax
-        flags *= self._checkMask(mask, sources, self.config.excludeMaskPlanes)
-        selectSources = sources[flags].copy(deep=True)
+
+        selected = self.sourceSelector.selectSources(sources).selected
+        nInitialSelected = np.count_nonzero(selected)
+        selected *= self._checkMask(mask, sources, self.config.excludeMaskPlanes)
+        nSelected = np.count_nonzero(selected)
+        self.log.info("Rejecting %i candidate sources: an excluded template mask plane is set.",
+                      nInitialSelected - nSelected)
+        selectSources = sources[selected].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):
             signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
             indices = np.argsort(signalToNoise)
             indices = indices[-self.config.maxKernelSources:]
-            flags = np.zeros(len(selectSources), dtype=bool)
-            flags[indices] = True
-            selectSources = selectSources[flags].copy(deep=True)
+            selected = np.zeros(len(selectSources), dtype=bool)
+            selected[indices] = True
+            selectSources = selectSources[selected].copy(deep=True)
 
         self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
                       len(selectSources), len(sources), 100*len(selectSources)/len(sources))
@@ -805,8 +816,8 @@ def _sourceSelector(self, sources, mask):
         return selectSources
 
     @staticmethod
-    def _checkMask(mask, sources, excludeMaskPlanes):
-        """Exclude sources that are located on masked pixels.
+    def _checkMask(mask, sources, excludeMaskPlanes, checkAdjacent=True):
+        """Exclude sources that are located on or adjacent to masked pixels.
 
         Parameters
         ----------
@@ -831,11 +842,18 @@ def _checkMask(mask, sources, excludeMaskPlanes):
 
         excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)
 
-        xv = np.rint(sources.getX() - mask.getX0())
-        yv = np.rint(sources.getY() - mask.getY0())
+        xv = (np.rint(sources.getX() - mask.getX0())).astype(int)
+        yv = (np.rint(sources.getY() - mask.getY0())).astype(int)
 
-        mv = mask.array[yv.astype(int), xv.astype(int)]
-        flags = np.bitwise_and(mv, excludePixelMask) == 0
+        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
 
     def _prepareInputs(self, template, science, visitSummary=None):

python/lsst/ip/diffim/utils.py

@@ -1095,6 +1095,8 @@ def detectTestSources(exposure, addMaskPlanes=None):
         addMaskPlanes = ["STREAK", "INJECTED", "INJECTED_TEMPLATE"]
 
     schema = afwTable.SourceTable.makeMinimalSchema()
+    # Add coordinate error fields:
+    afwTable.CoordKey.addErrorFields(schema)
     selectDetection = measAlg.SourceDetectionTask(schema=schema)
     selectMeasurement = measBase.SingleFrameMeasurementTask(schema=schema)
     table = afwTable.SourceTable.make(schema)
@@ -1267,6 +1269,8 @@ def _makeTruthSchema():
         Calib, Ap, and Psf flux slots all are set to ``truth_instFlux``.
     """
     schema = afwTable.SourceTable.makeMinimalSchema()
+    # Add coordinate error fields:
+    afwTable.CoordKey.addErrorFields(schema)
     keys = {}
     # Don't use a FluxResultKey so we can manage the flux and err separately.
     keys["instFlux"] = schema.addField("truth_instFlux", type=np.float64,
@@ -1280,6 +1284,11 @@ def _makeTruthSchema():
     schema.addField("base_PixelFlags_flag_interpolated", "Flag", "testing flag.")
     schema.addField("base_PixelFlags_flag_saturated", "Flag", "testing flag.")
     schema.addField("base_PixelFlags_flag_bad", "Flag", "testing flag.")
+    schema.addField("base_PixelFlags_flag_edge", "Flag", "testing flag.")
+    schema.addField("base_PsfFlux_flag", "Flag", "testing flag.")
+    schema.addField("base_ClassificationSizeExtendedness_value", "Flag", "testing flag.")
+    schema.addField("deblend_nChild", "Flag", "testing flag.")
+    schema.addField("detect_isPrimary", "Flag", "testing flag.")
     schema.getAliasMap().set("slot_Centroid", "truth")
     schema.getAliasMap().set("slot_CalibFlux", "truth")
     schema.getAliasMap().set("slot_ApFlux", "truth")
@@ -1316,6 +1325,11 @@ def _fillTruthCatalog(injectList):
         footprint.addPeak(x, y, flux)
         record.setFootprint(footprint)
 
+        # Set source records for isolated stars
+        record["base_ClassificationSizeExtendedness_value"] = 0
+        record["deblend_nChild"] = 0
+        record["detect_isPrimary"] = True
+
     return catalog
 
 

tests/test_detectAndMeasure.py

@@ -488,6 +488,8 @@ def test_fake_mask_plane_propagation(self):
         template_fake_masked = (template.mask.array & template_fake_bitmask) > 0
 
         subtractConfig = subtractImages.AlardLuptonSubtractTask.ConfigClass()
+        subtractConfig.sourceSelector.signalToNoise.fluxField = "truth_instFlux"
+        subtractConfig.sourceSelector.signalToNoise.errField = "truth_instFluxErr"
         subtractTask = subtractImages.AlardLuptonSubtractTask(config=subtractConfig)
         subtraction = subtractTask.run(template, science, sources)