Merge pull request #457 from lsst/tickets/DM-28542

DM-28542: Implement RFC-750

python/lsst/pipe/tasks/calibrate.py

@@ -703,7 +703,7 @@ def run(self, exposure, exposureIdInfo=None, background=None,
             )
         self.catalogCalculation.run(sourceCat)
 
-        self.setPrimaryFlags.run(sourceCat, includeDeblend=self.config.doDeblend)
+        self.setPrimaryFlags.run(sourceCat)
 
         if icSourceCat is not None and \
            len(self.config.icSourceFieldsToCopy) > 0:

python/lsst/pipe/tasks/multiBand.py

@@ -358,7 +358,7 @@ class DeblendCoaddSourcesConfig(Config):
     simultaneous = Field(dtype=bool,
                          default=True,
                          doc="Simultaneously deblend all bands? "
-                             "True uses 'singleBandDeblend' while False uses 'multibandDeblend'")
+                             "True uses `multibandDeblend` while False uses `singleBandDeblend`")
     coaddName = Field(dtype=str, default="deep", doc="Name of coadd")
     hasFakes = Field(dtype=bool,
                      default=False,
@@ -1089,7 +1089,7 @@ def run(self, exposure, sources, skyInfo, exposureId, ccdInputs=None, visitCatal
             self.catalogCalculation.run(sources)
 
         self.setPrimaryFlags.run(sources, skyMap=skyInfo.skyMap, tractInfo=skyInfo.tractInfo,
-                                 patchInfo=skyInfo.patchInfo, includeDeblend=self.deblended)
+                                 patchInfo=skyInfo.patchInfo)
         if self.config.doPropagateFlags:
             self.propagateFlags.run(butler, sources, ccdInputs, exposure.getWcs(), visitCatalogs, wcsUpdates)
 

python/lsst/pipe/tasks/setPrimaryFlags.py

@@ -20,15 +20,177 @@
 # the GNU General Public License along with this program.  If not,
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
-import numpy
+import numpy as np
 from lsst.pex.config import Config, Field, ListField
 from lsst.pipe.base import Task
 from lsst.geom import Box2D
 
 
+def getPatchInner(sources, patchInfo):
+    """Set a flag for each source if it is in the innerBBox of a patch.
+
+    Parameters
+    ----------
+    sources : `lsst.afw.table.SourceCatalog`
+        A sourceCatalog with pre-calculated centroids.
+    patchInfo : `lsst.skymap.PatchInfo`
+        Information about a `SkyMap` `Patch`.
+
+    Returns
+    --------
+    isPatchInner : array-like of `bool`
+        `True` for each source that has a centroid
+        in the inner region of a patch.
+    """
+    # Extract the centroid position for all the sources
+    x = sources["slot_Centroid_x"]
+    y = sources["slot_Centroid_y"]
+    centroidFlag = sources["slot_Centroid_flag"]
+
+    # set inner flags for each source and set primary flags for
+    # sources with no children (or all sources if deblend info not available)
+    innerFloatBBox = Box2D(patchInfo.getInnerBBox())
+    inInner = innerFloatBBox.contains(x, y)
+
+    # When the centroider fails, we can still fall back to the peak,
+    # but we don't trust that quite as much -
+    # so we use a slightly smaller box for the patch comparison.
+    shrunkInnerFloatBBox = Box2D(innerFloatBBox)
+    shrunkInnerFloatBBox.grow(-1)
+    inShrunkInner = shrunkInnerFloatBBox.contains(x, y)
+
+    # Flag sources contained in the inner region of a patch
+    isPatchInner = (centroidFlag & inShrunkInner) | (~centroidFlag & inInner)
+    return isPatchInner
+
+
+def getTractInner(sources, tractInfo, skyMap):
+    """Set a flag for each source that the skyMap includes in tractInfo.
+
+    Parameters
+    ----------
+    sources : `lsst.afw.table.SourceCatalog`
+        A sourceCatalog with pre-calculated centroids.
+    tractInfo : `lsst.skymap.TractInfo`
+        Tract object
+    skyMap : `lsst.skymap.BaseSkyMap`
+        Sky tessellation object
+
+    Returns
+    -------
+    isTractInner : array-like of `bool`
+        True if the skyMap.findTract method returns
+        the same tract as tractInfo.
+    """
+    tractId = tractInfo.getId()
+    isTractInner = np.array([skyMap.findTract(s.getCoord()).getId() == tractId for s in sources])
+    return isTractInner
+
+
+def getPseudoSources(sources, pseudoFilterList, schema, log):
+    """Get a flag that marks pseudo sources.
+
+    Some categories of sources, for example sky objects,
+    are not really detected sources and should not be considered primary
+    sources.
+
+    Parameters
+    ----------
+    sources : `lsst.afw.table.SourceCatalog`
+        The catalog of sources for which to identify "pseudo"
+        (e.g. sky) objects.
+    pseudoFilterList : `list` of `str`
+        Names of filters which should never be primary
+
+    Returns
+    -------
+    isPseudo : array-like of `bool`
+        True for each source that is a pseudo source.
+        Note: to remove pseudo sources use `~isPseudo`.
+    """
+    # Filter out sources that should never be primary
+    isPseudo = np.zeros(len(sources), dtype=bool)
+    for filt in pseudoFilterList:
+        try:
+            pseudoFilterKey = schema.find("merge_peak_%s" % filt).getKey()
+            isPseudo |= sources[pseudoFilterKey]
+        except KeyError:
+            log.warn("merge_peak is not set for pseudo-filter %s" % filt)
+    return isPseudo
+
+
+def getDeblendPrimaryFlags(sources):
+    """Get flags generated by the deblender
+
+    scarlet is different than meas_deblender in that it is not
+    (necessarily) flux conserving. For consistency in scarlet,
+    all of the parents with only a single child (isolated sources)
+    need to be deblended. This creates a question: which type
+    of isolated source should we make measurements on, the
+    undeblended "parent" or the deblended child?
+    For that reason we distinguish between a DeblendedSource,
+    which is a source that has no children and uses the
+    isolated parents, and a DeblendedModelSource, which uses
+    the scarlet models for both isolated and blended sources.
+    In the case of meas_deblender, DeblendedModelSource is
+    `None` because it is not contained in the output catalog.
+
+    Parameters
+    ----------
+    sources : `lsst.afw.table.SourceCatalog`
+        A sourceCatalog that has already been deblended using
+        either meas_extensions_scarlet or meas_deblender.
+
+    Returns
+    -------
+    fromBlend : array-like of `bool`
+        True for each source modeled by the deblender from a `Peak`
+        in a parent footprint that contained at least one other `Peak`.
+        While these models can be approximated as isolated,
+        and measurements are made on them as if that's the case,
+        we know deblending to introduce biases in the shape and centroid
+        of objects and it is important to know that the sources that these
+        models are based on are all bleneded in the true image.
+    isIsolated : array-like of `bool`
+        True for isolated sources, regardless of whether or not they
+        were modeled by the deblender.
+    isDeblendedSource : array-like of `bool`
+        True for each source that is a "DeblendedSource" as defined above.
+    isDeblendedModelSource : array-like of `bool`
+        True for each source that is a "DeblendedSourceModel"
+        as defined above.
+    """
+    nChildKey = "deblend_nChild"
+    nChild = sources[nChildKey]
+    parent = sources["parent"]
+
+    if "deblend_scarletFlux" in sources.schema:
+        # The number of peaks in the sources footprint.
+        # This (may be) different than nChild,
+        # the number of deblended sources in the catalog,
+        # because some peaks might have been culled during deblending.
+        nPeaks = sources["deblend_nPeaks"]
+        parentNChild = sources["deblend_parentNChild"]
+        # It is possible for a catalog to contain a hierarchy of sources,
+        # so we mark the leaves (end nodes of the hierarchy tree with no
+        # children).
+        isLeaf = nPeaks == 1
+        fromBlend = parentNChild > 1
+        isIsolated = isLeaf & ((parent == 0) | parentNChild == 1)
+        isDeblendedSource = (fromBlend & isLeaf) | (isIsolated & (parent == 0))
+        isDeblendedModelSource = (parent != 0) & isLeaf
+    else:
+        # Set the flags for meas_deblender
+        fromBlend = parent != 0
+        isIsolated = (nChild == 0) & (parent == 0)
+        isDeblendedSource = nChild == 0
+        isDeblendedModelSource = None
+    return fromBlend, isIsolated, isDeblendedSource, isDeblendedModelSource
+
+
 class SetPrimaryFlagsConfig(Config):
-    nChildKeyName = Field(dtype=str, default="deblend_nChild",
-                          doc="Name of field in schema with number of deblended children")
+    nChildKeyName = Field(dtype=str, default="deprecated",
+                          doc="Deprecated. This parameter is not used.")
     pseudoFilterList = ListField(dtype=str, default=['sky'],
                                  doc="Names of filters which should never be primary")
 
@@ -42,6 +204,9 @@ class SetPrimaryFlagsTask(Task):
         The input schema.
     isSingleFrame : `bool`
         Flag specifying if task is operating with single frame imaging.
+    includeDeblend : `bool`
+        Include deblend information in isPrimary and
+        add isDeblendedSource field?
     kwargs :
         Keyword arguments passed to the task.
     """
@@ -52,6 +217,7 @@ def __init__(self, schema, isSingleFrame=False, **kwargs):
         Task.__init__(self, **kwargs)
         self.schema = schema
         self.isSingleFrame = isSingleFrame
+        self.includeDeblend = False
         if not self.isSingleFrame:
             primaryDoc = ("true if source has no children and is in the inner region of a coadd patch "
                           "and is in the inner region of a coadd tract "
@@ -71,14 +237,35 @@ def __init__(self, schema, isSingleFrame=False, **kwargs):
             doc=primaryDoc,
         )
 
-    def run(self, sources, skyMap=None, tractInfo=None, patchInfo=None,
-            includeDeblend=True):
-        """Set is-patch-inner, is-tract-inner and is-primary flags on sources.
-        For coadded imaging, the is-primary flag returns True when an object
+        if "deblend_nChild" in schema.getNames():
+            self.includeDeblend = True
+            self.isDeblendedSourceKey = self.schema.addField(
+                "detect_isDeblendedSource", type="Flag",
+                doc=primaryDoc + " and is either an unblended isolated source or a "
+                                 "deblended child from a parent with 'deblend_nChild' > 1")
+            self.fromBlendKey = self.schema.addField(
+                "detect_fromBlend", type="Flag",
+                doc="This source is deblended from a parent with more than one child."
+            )
+            self.isIsolatedKey = self.schema.addField(
+                "detect_isIsolated", type="Flag",
+                doc="This source is not a part of a blend."
+            )
+            if "deblend_scarletFlux" in schema.getNames():
+                self.isDeblendedModelKey = self.schema.addField(
+                    "detect_isDeblendedModelSource", type="Flag",
+                    doc=primaryDoc + " and is a deblended child")
+            else:
+                self.isDeblendedModelKey = None
+
+    def run(self, sources, skyMap=None, tractInfo=None, patchInfo=None):
+        """Set isPrimary and related flags on sources.
+
+        For coadded imaging, the `isPrimary` flag returns True when an object
         has no children, is in the inner region of a coadd patch, is in the
         inner region of a coadd trach, and is not detected in a pseudo-filter
         (e.g., a sky_object).
-        For single frame imaging, the is-primary flag returns True when a
+        For single frame imaging, the isPrimary flag returns True when a
         source has no children and is not a sky source.
 
         Parameters
@@ -92,67 +279,34 @@ def run(self, sources, skyMap=None, tractInfo=None, patchInfo=None,
             Tract object
         patchInfo : `lsst.skymap.PatchInfo`
             Patch object
-        includeDeblend : `bool`
-            Include deblend information in isPrimary?
         """
-        nChildKey = None
-        if includeDeblend:
-            nChildKey = self.schema.find(self.config.nChildKeyName).key
-
-        # coadd case
+        # Mark whether sources are contained within the inner regions of the
+        # given tract/patch and are not "pseudo" (e.g. sky) sources.
         if not self.isSingleFrame:
-            # set inner flags for each source and set primary flags for sources with no children
-            # (or all sources if deblend info not available)
-            innerFloatBBox = Box2D(patchInfo.getInnerBBox())
-
-            # When the centroider fails, we can still fall back to the peak, but we don't trust
-            # that quite as much - so we use a slightly smaller box for the patch comparison.
-            # That's trickier for the tract comparison, so we just use the peak without extra
-            # care there.
-            shrunkInnerFloatBBox = Box2D(innerFloatBBox)
-            shrunkInnerFloatBBox.grow(-1)
-
-            pseudoFilterKeys = []
-            for filt in self.config.pseudoFilterList:
-                try:
-                    pseudoFilterKeys.append(self.schema.find("merge_peak_%s" % filt).getKey())
-                except Exception:
-                    self.log.warn("merge_peak is not set for pseudo-filter %s" % filt)
-
-            tractId = tractInfo.getId()
-            for source in sources:
-                centroidPos = source.getCentroid()
-                if numpy.any(numpy.isnan(centroidPos)):
-                    continue
-                if source.getCentroidFlag():
-                    # Use a slightly smaller box to guard against bad centroids (see above)
-                    isPatchInner = shrunkInnerFloatBBox.contains(centroidPos)
-                else:
-                    isPatchInner = innerFloatBBox.contains(centroidPos)
-                source.setFlag(self.isPatchInnerKey, isPatchInner)
-
-                skyPos = source.getCoord()
-                sourceInnerTractId = skyMap.findTract(skyPos).getId()
-                isTractInner = sourceInnerTractId == tractId
-                source.setFlag(self.isTractInnerKey, isTractInner)
-
-                if nChildKey is None or source.get(nChildKey) == 0:
-                    for pseudoFilterKey in pseudoFilterKeys:
-                        if source.get(pseudoFilterKey):
-                            isPseudo = True
-                            break
-                    else:
-                        isPseudo = False
-
-                    source.setFlag(self.isPrimaryKey, isPatchInner and isTractInner and not isPseudo)
-
-        # single frame case
+            isPatchInner = getPatchInner(sources, patchInfo)
+            isTractInner = getTractInner(sources, tractInfo, skyMap)
+            isPseudo = getPseudoSources(sources, self.config.pseudoFilterList, self.schema, self.log)
+            isPrimary = isTractInner & isPatchInner & ~isPseudo
+
+            sources[self.isPatchInnerKey] = isPatchInner
+            sources[self.isTractInnerKey] = isTractInner
         else:
-            hasSkySources = True if "sky_source" in sources.schema else False
-            for source in sources:
-                hasNoChildren = True if nChildKey is None or source.get(nChildKey) == 0 else False
-                isSkySource = False
-                if hasSkySources:
-                    if source["sky_source"]:
-                        isSkySource = True
-                source.setFlag(self.isPrimaryKey, hasNoChildren and not isSkySource)
+            # Mark all of the sky sources in SingleFrame images
+            # (if they were added)
+            if "sky_source" in sources.schema:
+                isSky = sources["sky_source"]
+            else:
+                isSky = np.zeros(len(sources), dtype=bool)
+            isPrimary = ~isSky
+
+        if self.includeDeblend:
+            result = getDeblendPrimaryFlags(sources)
+            fromBlend, isIsolated, isDeblendedSource, isDeblendedModelSource = result
+            sources[self.fromBlendKey] = fromBlend
+            sources[self.isIsolatedKey] = isIsolated
+            sources[self.isDeblendedSourceKey] = isDeblendedSource
+            if self.isDeblendedModelKey is not None:
+                sources[self.isDeblendedModelKey] = isDeblendedModelSource
+            isPrimary = isPrimary & isDeblendedSource
+
+        sources[self.isPrimaryKey] = isPrimary