Skip to content

tickets/DM-44456: Re-implement deblender refactor from DM-44319 #320

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fred3m wants to merge 4 commits into
base: main
Choose a base branch
from
Open

tickets/DM-44456: Re-implement deblender refactor from DM-44319 #320

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
61a09b2
Revert "Merge pull request #317 from lsst/tickets/DM-44319-revert"
fred3m May 21, 2024
e94c74d
Ensure that the deblended catalog is contiguous
fred3m May 21, 2024
746b71c
Re-add merged/unmerged metadata
fred3m May 21, 2024
0daec4f
Remove parents and children of bad sources
fred3m May 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
197 changes: 112 additions & 85 deletions python/lsst/ip/diffim/detectAndMeasure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,6 @@

import numpy as np

import lsst.afw.detection as afwDetection
import lsst.afw.table as afwTable
import lsst.daf.base as dafBase
import lsst.geom
Expand Down Expand Up @@ -77,6 +76,12 @@ class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections,
storageClass="SourceCatalog",
name="{fakesType}{coaddName}Diff_diaSrc",
)
removedSources = pipeBase.connectionTypes.Output(
doc="Sources and blends removed from the diaSources catalog.",
dimensions=("instrument", "visit", "detector"),
storageClass="SourceCatalog",
name="{fakesType}{coaddName}Diff_diaSrc_removed",
)
subtractedMeasuredExposure = pipeBase.connectionTypes.Output(
doc="Difference image with detection mask plane filled in.",
dimensions=("instrument", "visit", "detector"),
Expand Down Expand Up @@ -104,7 +109,8 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
dtype=bool,
default=True,
doc="Merge positive and negative diaSources with grow radius "
"set by growFootprint"
"set by growFootprint",
deprecated="This field is no longer used and will be removed after v28."
)
doForcedMeasurement = pexConfig.Field(
dtype=bool,
Expand Down Expand Up @@ -143,7 +149,8 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
growFootprint = pexConfig.Field(
dtype=int,
default=2,
doc="Grow positive and negative footprints by this many pixels before merging"
doc="Grow positive and negative footprints by this many pixels before merging",
deprecated="This field is no longer used and will be removed after v28."
)
diaSourceMatchRadius = pexConfig.Field(
dtype=float,
Expand Down Expand Up @@ -350,16 +357,11 @@ def run(self, science, matchedTemplate, difference,
doSmooth=True,
)

sources, positives, negatives = self._deblend(difference,
results.positive,
results.negative)
sources = self._buildCatalogAndDeblend(difference, results.positive, results.negative, idFactory)

return self.processResults(science, matchedTemplate, difference, sources, idFactory,
positiveFootprints=positives,
negativeFootprints=negatives)
return self._measureSources(science, matchedTemplate, difference, sources)

def processResults(self, science, matchedTemplate, difference, sources, idFactory,
positiveFootprints=None, negativeFootprints=None,):
def _measureSources(self, science, matchedTemplate, difference, initialDiaSources):
"""Measure and process the results of source detection.

Parameters
Expand All @@ -371,15 +373,8 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
difference image.
difference : `lsst.afw.image.ExposureF`
Result of subtracting template from the science image.
sources : `lsst.afw.table.SourceCatalog`
initialDiaSources : `lsst.afw.table.SourceCatalog`
Detected sources on the difference exposure.
idFactory : `lsst.afw.table.IdFactory`
Generator object used to assign ids to detected sources in the
difference image.
positiveFootprints : `lsst.afw.detection.FootprintSet`, optional
Positive polarity footprints.
negativeFootprints : `lsst.afw.detection.FootprintSet`, optional
Negative polarity footprints.

Returns
-------
Expand All @@ -390,38 +385,14 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
``diaSources`` : `lsst.afw.table.SourceCatalog`
The catalog of detected sources.
"""
self.metadata.add("nUnmergedDiaSources", len(sources))
if self.config.doMerge:
fpSet = positiveFootprints
fpSet.merge(negativeFootprints, self.config.growFootprint,
self.config.growFootprint, False)
initialDiaSources = afwTable.SourceCatalog(self.schema)
fpSet.makeSources(initialDiaSources)
self.log.info("Merging detections into %d sources", len(initialDiaSources))
else:
initialDiaSources = sources

# Assign source ids at the end: deblend/merge mean that we don't keep
# track of parents and children, we only care about the final ids.
for source in initialDiaSources:
source.setId(idFactory())
# Ensure sources added after this get correct ids.
initialDiaSources.getTable().setIdFactory(idFactory)
initialDiaSources.setMetadata(self.algMetadata)

self.metadata.add("nMergedDiaSources", len(initialDiaSources))
initialDiaSources.setMetadata(self.algMetadata)

if self.config.doMaskStreaks:
streakInfo = self._runStreakMasking(difference.maskedImage)

if self.config.doSkySources:
self.addSkySources(initialDiaSources, difference.mask, difference.info.id)

if not initialDiaSources.isContiguous():
initialDiaSources = initialDiaSources.copy(deep=True)

self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
diaSources = self._removeBadSources(initialDiaSources)
diaSources, removedSources = self._removeBadSources(initialDiaSources)

if self.config.doForcedMeasurement:
self.measureForcedSources(diaSources, science, difference.getWcs())
Expand All @@ -431,15 +402,21 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
measurementResults = pipeBase.Struct(
subtractedMeasuredExposure=difference,
diaSources=diaSources,
removedSources=removedSources,
)
if self.config.doMaskStreaks and self.config.writeStreakInfo:
measurementResults.mergeItems(streakInfo, 'maskedStreaks')

return measurementResults

def _deblend(self, difference, positiveFootprints, negativeFootprints):
"""Deblend the positive and negative footprints and return a catalog
containing just the children, and the deblended footprints.
def _buildCatalogAndDeblend(self, difference, positiveFootprints, negativeFootprints, idFactory):
"""Create a source catalog and deblend when possible.

This method creates a source catalog from the positive and negative
footprints, and then deblends the footprints that have only positive
peaks. This is because `lsst.meas.deblender.SourceDeblendTask` is not
designed to handle footprints that have negative flux, resulting in
garbage output.

Parameters
----------
Expand All @@ -448,6 +425,9 @@ def _deblend(self, difference, positiveFootprints, negativeFootprints):
positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet`
Positive and negative polarity footprints measured on
``difference`` to be deblended separately.
idFactory : `lsst.afw.table.IdFactory`
Generator object used to assign ids to detected sources in the
difference image.

Returns
-------
Expand All @@ -457,33 +437,67 @@ def _deblend(self, difference, positiveFootprints, negativeFootprints):
Deblended positive and negative polarity footprints measured on
``difference``.
"""
def makeFootprints(sources):
footprints = afwDetection.FootprintSet(difference.getBBox())
footprints.setFootprints([src.getFootprint() for src in sources])
return footprints

def deblend(footprints):
"""Deblend a positive or negative footprint set,
and return the deblended children.
"""
sources = afwTable.SourceCatalog(self.schema)
footprints.makeSources(sources)
self.deblend.run(exposure=difference, sources=sources)
self.setPrimaryFlags.run(sources)
children = sources["detect_isDeblendedSource"] == 1
sources = sources[children].copy(deep=True)
# Clear parents, so that measurement plugins behave correctly.
sources['parent'] = 0
return sources.copy(deep=True)

positives = deblend(positiveFootprints)
negatives = deblend(negativeFootprints)

sources = afwTable.SourceCatalog(self.schema)
sources.reserve(len(positives) + len(negatives))
sources.extend(positives, deep=True)
sources.extend(negatives, deep=True)
return sources, makeFootprints(positives), makeFootprints(negatives)
self.metadata.add(
"nUnmergedDiaSources",
len(positiveFootprints.getFootprints()) + len(negativeFootprints.getFootprints())
)
# Merge the positive and negative footprints.
# The original detection FootprintSets already grew each detection,
# so there is no need to grow them again.
merged_footprints = positiveFootprints
merged_footprints.merge(negativeFootprints, 0, 0, False)

# Create a source catalog from the footprints.
table = afwTable.SourceTable.make(self.schema, idFactory)
sources = afwTable.SourceCatalog(table)
merged_footprints.makeSources(sources)

# Sky sources must be added before deblending, otherwise the
# sources with parent == 0 will be out of order and
# downstream measurement tasks cannot run.
if self.config.doSkySources:
self.addSkySources(sources, difference.mask, difference.info.id)

# Find the footprints with only positive peaks and no negative peaks.
footprints = [src.getFootprint() for src in sources]
nPeaks = np.array([len(fp.peaks) for fp in footprints])
blend_ids = []
skipped_ids = []
for src in sources[nPeaks > 1]:
peaks = src.getFootprint().peaks
peak_flux = np.array([peak.getPeakValue() for peak in peaks])
if np.all(peak_flux > 0) or np.all(peak_flux < 0):
blend_ids.append(src.getId())
else:
skipped_ids.append(src.getId())

# Deblend the footprints that can be deblended with SourceDeblendTask.
temp_cat = afwTable.SourceCatalog(sources.table)
for sid in blend_ids:
temp_cat.append(sources.find(sid))
first_child_index = len(temp_cat)
self.deblend.run(exposure=difference, sources=temp_cat)

# Update the deblended parents and add their children
# to the sources catalog.
sources.extend(temp_cat[first_child_index:], deep=False)

# Set deblending flags.
# Since SourceDeblendTask already has a method for setting
# flags for skipped parents, we just call that method here
# instead of setting the flags manually.
for sid in skipped_ids:
src = sources.find(sid)
self.deblend.skipParent(src, difference.mask)

# Set detection and primary flags
self.setPrimaryFlags.run(sources)

table = afwTable.SourceTable.make(self.schema, idFactory)
_sources = afwTable.SourceCatalog(table)
_sources.extend(sources, deep=True)

return _sources

def _removeBadSources(self, diaSources):
"""Remove unphysical diaSources from the catalog.
Expand All @@ -499,17 +513,33 @@ def _removeBadSources(self, diaSources):
The updated catalog of detected sources, with any source that has a
flag in ``config.badSourceFlags`` set removed.
"""
selector = np.ones(len(diaSources), dtype=bool)
bad_sources = np.zeros(len(diaSources), dtype=bool)
for flag in self.config.badSourceFlags:
flags = diaSources[flag]
nBad = np.count_nonzero(flags)
if nBad > 0:
self.log.debug("Found %d unphysical sources with flag %s.", nBad, flag)
selector &= ~flags
nBadTotal = np.count_nonzero(~selector)
bad_sources |= flags
# Remove parents of bad children, and children of bad parents.
# This is due to the measurement plugins, where it is assumed that
# parent blends have all of their children in the catalog, and vice versa.
parents = (
(diaSources["parent"] == 0)
& (diaSources["deblend_nChild"] > 1)
& (~diaSources["deblend_skipped"])
)
children = (diaSources["parent"] != 0)
parents_to_remove = np.unique(diaSources["parent"][children & bad_sources])
bad_parents = parents & bad_sources
bad_parents |= np.in1d(diaSources["id"], parents_to_remove)
children_to_remove = np.in1d(diaSources["parent"], diaSources["id"][bad_parents])
bad_children = (bad_sources & children) | children_to_remove
bad_sources |= bad_parents | bad_children

nBadTotal = np.count_nonzero(bad_sources)
self.metadata.add("nRemovedBadFlaggedSources", nBadTotal)
self.log.info("Removed %d unphysical sources.", nBadTotal)
return diaSources[selector].copy(deep=True)
self.log.info("Removed %d unphysical sources and blends.", nBadTotal)
return diaSources[~bad_sources].copy(deep=True), diaSources[bad_sources]

def addSkySources(self, diaSources, mask, seed,
subtask=None):
Expand Down Expand Up @@ -731,9 +761,6 @@ def run(self, science, matchedTemplate, difference, scoreExposure,
# Copy the detection mask from the Score image to the difference image
difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())

sources, positives, negatives = self._deblend(difference,
results.positive,
results.negative)
sources = self._buildCatalogAndDeblend(difference, results.positive, results.negative, idFactory)

return self.processResults(science, matchedTemplate, difference, sources, idFactory,
positiveFootprints=positives, negativeFootprints=negatives)
return self._measureSources(science, matchedTemplate, difference, sources)
14 changes: 10 additions & 4 deletions tests/test_detectAndMeasure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -247,8 +247,9 @@ def test_remove_unphysical(self):
badDiaSrc = ~bbox.contains(diaSources.getX(), diaSources.getY())
nBad = np.count_nonzero(badDiaSrc)
self.assertEqual(nBad, nSetBad)
diaSourcesNoBad = detectionTask._removeBadSources(diaSources)
diaSourcesNoBad, removedSources = detectionTask._removeBadSources(diaSources)
badDiaSrcNoBad = ~bbox.contains(diaSourcesNoBad.getX(), diaSourcesNoBad.getY())
self.assertEqual(len(removedSources), nSetBad)

# Verify that no sources outside the image bounding box remain
self.assertEqual(np.count_nonzero(badDiaSrcNoBad), 0)
Expand Down Expand Up @@ -288,7 +289,9 @@ def _detection_wrapper(positive=True):
output = detectionTask.run(science.clone(), matchedTemplate, difference)
refIds = []
scale = 1. if positive else -1.
for diaSource in output.diaSources:
# Deblending provides multiple copies of the main parent source,
# so we only check the parents. Deblending is checked in another test.
for diaSource in output.diaSources[output.diaSources["parent"] == 0]:
self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
_detection_wrapper(positive=True)
_detection_wrapper(positive=False)
Expand Down Expand Up @@ -517,7 +520,7 @@ def test_fake_mask_plane_propagation(self):

sci_refIds = []
tmpl_refIds = []
for diaSrc in output.diaSources:
for diaSrc in output.diaSources[output.diaSources["parent"] == 0]:
if diaSrc['base_PsfFlux_instFlux'] > 0:
self._check_diaSource(science_fake_sources, diaSrc, scale=1, refIds=sci_refIds)
self.assertTrue(diaSrc['base_PixelFlags_flag_injected'])
Expand Down Expand Up @@ -687,7 +690,10 @@ def _detection_wrapper(positive=True):
scale = 1. if positive else -1.
# sources near the edge may have untrustworthy centroids
goodSrcFlags = ~output.diaSources['base_PixelFlags_flag_edge']
for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
# Deblending provides multiple copies of the main parent source,
# so we only check the parents. Deblending is checked in another test.
for diaSource, goodSrcFlag in zip(output.diaSources[output.diaSources["parent"] == 0],
goodSrcFlags):
if goodSrcFlag:
self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
_detection_wrapper(positive=True)
Expand Down