From 4ce7458ce26d0aaf866eb32eac788e5f1068e8f1 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Tue, 21 Nov 2023 15:43:35 -0800
Subject: [PATCH] Remove unphysical diaSources from the output of
 detectAndMeasure

---
 python/lsst/ip/diffim/detectAndMeasure.py | 35 +++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 0dc10a4d..2800fbaf 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -20,6 +20,7 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from deprecated.sphinx import deprecated
+import numpy as np
 
 import lsst.afw.table as afwTable
 import lsst.daf.base as dafBase
@@ -140,6 +141,18 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
         target=SkyObjectsTask,
         doc="Generate sky sources",
     )
+    badSourceFlags = lsst.pex.config.ListField(
+        dtype=str,
+        doc="Flags that, if set, the associated source should not "
+        "be included in the output diaSource catalog.",
+        default=("base_PixelFlags_flag_offimage",
+                 "base_PixelFlags_flag_interpolatedCenterAll",
+                 "base_PixelFlags_flag_saturatedCenterAll",
+                 "base_PixelFlags_flag_crCenterAll",
+                 "base_PixelFlags_flag_badCenterAll",
+                 "base_PixelFlags_flag_suspectCenterAll",
+                 ),
+    )
     idGenerator = DetectorVisitIdGeneratorConfig.make_field()
 
     def setDefaults(self):
@@ -361,6 +374,7 @@ def processResults(self, science, matchedTemplate, difference, sources, table,
             self.addSkySources(diaSources, difference.mask, difference.info.id)
 
         self.measureDiaSources(diaSources, science, difference, matchedTemplate)
+        self.purgeSources(diaSources)
 
         if self.config.doForcedMeasurement:
             self.measureForcedSources(diaSources, science, difference.getWcs())
@@ -372,6 +386,27 @@ def processResults(self, science, matchedTemplate, difference, sources, table,
 
         return measurementResults
 
+    def purgeSources(self, diaSources):
+        """Remove bad diaSources from the catalog.
+
+        Parameters
+        ----------
+        diaSources : `lsst.afw.table.SourceCatalog`
+            The catalog of detected sources. Modified in place if any sources
+            have a flag in `badSourceFlags` set.
+        """
+        flags = np.ones(len(diaSources), dtype=bool)
+        for flag in self.config.badSourceFlags:
+            try:
+                flags *= ~diaSources[flag]
+            except Exception as e:
+                self.log.warning("Could not apply source flag: %s", e)
+        nPurged = np.sum(~flags)
+        if nPurged > 0:
+            self.log.warning(f"Found and purged {nPurged} unphysical sources.")
+        diaSources = diaSources[flags].copy(deep=True)
+        self.metadata.add("nPurgedSources", nPurged)
+
     def addSkySources(self, diaSources, mask, seed):
         """Add sources in empty regions of the difference image
         for measuring the background.