Merge pull request #93 from lsst/tickets/DM-43176

DM-43176: Flag sources with no flux after re-distribution.

python/lsst/meas/extensions/scarlet/io.py

@@ -26,7 +26,8 @@
 
 from lsst.afw.table import SourceCatalog
 from lsst.afw.image import MaskedImage, Exposure
-from lsst.afw.detection import Footprint as afwFootprint
+from lsst.afw.detection import Footprint as afwFootprint, HeavyFootprintF
+from lsst.afw.geom import SpanSet, Span
 from lsst.geom import Box2I, Extent2I, Point2I
 import lsst.scarlet.lite as scl
 from lsst.scarlet.lite import Blend, Source, Box, Component, FixedParameter, FactorizedComponent, Image
@@ -333,9 +334,24 @@ def updateBlendRecords(
             blend=blend,
             useFlux=useFlux,
         )
-        sourceRecord.setFootprint(heavy)
 
         if updateFluxColumns:
+            if heavy.getArea() == 0:
+                # The source has no flux after being weighted with the PSF
+                # in this particular band (it might have flux in others).
+                sourceRecord.set("deblend_zeroFlux", True)
+                # Create a Footprint with a single pixel, set to zero,
+                # to avoid breakage in measurement algorithms.
+                center = Point2I(heavy.peaks[0]["i_x"], heavy.peaks[0]["i_y"])
+                spanList = [Span(center.y, center.x, center.x)]
+                footprint = afwFootprint(SpanSet(spanList))
+                footprint.setPeakCatalog(heavy.peaks)
+                heavy = HeavyFootprintF(footprint)
+                heavy.getImageArray()[0] = 0.0
+            else:
+                sourceRecord.set("deblend_zeroFlux", False)
+            sourceRecord.setFootprint(heavy)
+
             if useFlux:
                 # Set the fraction of pixels with valid data.
                 coverage = calculateFootprintCoverage(heavy, imageForRedistribution.mask)

python/lsst/meas/extensions/scarlet/scarletDeblendTask.py

@@ -783,6 +783,8 @@ def _addSchemaKeys(self, schema):
         self.coverageKey = schema.addField('deblend_dataCoverage', type=np.float32,
                                            doc='Fraction of pixels with data. '
                                                'In other words, 1 - fraction of pixels with NO_DATA set.')
+        self.zeroFluxKey = schema.addField("deblend_zeroFlux", type="Flag",
+                                           doc="Source has zero flux.")
         # Blendedness/classification metrics
         self.maxOverlapKey = schema.addField("deblend_maxOverlap", type=np.float32,
                                              doc="Maximum overlap with all of the other neighbors flux "

tests/test_deblend.py

@@ -77,6 +77,55 @@ def setUp(self):
         for b, coadd in enumerate(self.coadds):
             coadd.setPsf(psfs[b])
 
+    def _insert_blank_source(self, modelData, catalog):
+        # Add parent
+        parent = catalog.addNew()
+        parent.setParent(0)
+        parent["deblend_nChild"] = 1
+        parent["deblend_nPeaks"] = 1
+        ss = SpanSet.fromShape(5, Stencil.CIRCLE, offset=(30, 70))
+        footprint = Footprint(ss)
+        peak = footprint.addPeak(30, 70, 0)
+        parent.setFootprint(footprint)
+
+        # Add the zero flux source
+        dtype = np.float32
+        center = (70, 30)
+        origin = (center[0]-5, center[1]-5)
+        psf = list(modelData.blends.values())[0].psf
+        src = catalog.addNew()
+        src.setParent(parent.getId())
+        src["deblend_peak_center_x"] = center[1]
+        src["deblend_peak_center_y"] = center[0]
+        src["deblend_nPeaks"] = 1
+
+        sources = {
+            src.getId(): {
+                "components": [],
+                "factorized": [{
+                    "origin": origin,
+                    "peak": center,
+                    "spectrum": np.zeros((len(self.bands),), dtype=dtype),
+                    "morph": np.zeros((11, 11), dtype=dtype),
+                    "shape": (11, 11),
+                }],
+                "peak_id": peak.getId(),
+            }
+        }
+
+        blendData = scl.io.ScarletBlendData.from_dict({
+            "origin": origin,
+            "shape": (11, 11),
+            "psf_center": center,
+            "psf_shape": psf.shape,
+            "psf": psf.flatten(),
+            "sources": sources,
+            "bands": self.bands,
+        })
+        pid = parent.getId()
+        modelData.blends[pid] = blendData
+        return pid, src.getId()
+
     def _deblend(self, version):
         schema = SourceCatalog.Table.makeMinimalSchema()
         # Adjust config options to test skipping parents
@@ -117,6 +166,8 @@ def _deblend(self, version):
     def test_deblend_task(self):
         catalog, modelData, config = self._deblend("lite")
 
+        bad_blend_id, bad_src_id = self._insert_blank_source(modelData, catalog)
+
         # Attach the footprints in each band and compare to the full
         # data model. This is done in each band, both with and without
         # flux re-distribution to test all of the different possible
@@ -153,6 +204,8 @@ def test_deblend_task(self):
                     self.assertEqual(len(children), parent.get("deblend_nChild"))
                     # Check that parent columns are propagated
                     # to their children
+                    if parent.getId() == bad_blend_id:
+                        continue
                     for parentCol, childCol in config.columnInheritance.items():
                         np.testing.assert_array_equal(parent.get(parentCol), children[childCol])
 
@@ -256,6 +309,10 @@ def test_deblend_task(self):
         skipped = largeFootprint | denseFootprint
         np.testing.assert_array_equal(skipped, catalog["deblend_skipped"])
 
+        # Check that the zero flux source was flagged
+        for src in catalog:
+            np.testing.assert_equal(src["deblend_zeroFlux"], src.getId() == bad_src_id)
+
     def test_continuity(self):
         """This test ensures that lsst.scarlet.lite gives roughly the same
         result as scarlet.lite