Skip to content

DM-48701: Fix incorrect footprint merging in diffim #379

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

Merged
parejkoj merged 2 commits into from
Feb 14, 2025
Merged

DM-48701: Fix incorrect footprint merging in diffim #379

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
159911c
Add basic tests that purely positive sources are deblended
parejkoj Feb 12, 2025
2cf1bc2
Use FootprintMergeList instead of FootprintSet.merge
parejkoj Feb 5, 2025
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
63 changes: 39 additions & 24 deletions python/lsst/ip/diffim/detectAndMeasure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -399,14 +399,18 @@ def run(self, science, matchedTemplate, difference,

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

else:
positives = afwTable.SourceCatalog(self.schema)
results.positive.makeSources(positives)
negatives = afwTable.SourceCatalog(self.schema)
results.negative.makeSources(negatives)
return self.processResults(science, matchedTemplate, difference,
results.sources, idFactory,
positiveFootprints=results.positive,
negativeFootprints=results.negative)
positives=positives,
negatives=negatives)

def _prepareInputs(self, difference):
"""Ensure that we start with an empty detection and deblended mask.
Expand Down Expand Up @@ -434,7 +438,7 @@ def _prepareInputs(self, difference):
mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)

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

Parameters
Expand All @@ -451,9 +455,9 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
idFactory : `lsst.afw.table.IdFactory`
Generator object used to assign ids to detected sources in the
difference image.
positiveFootprints : `lsst.afw.detection.FootprintSet`, optional
positives : `lsst.afw.table.SourceCatalog`, optional
Positive polarity footprints.
negativeFootprints : `lsst.afw.detection.FootprintSet`, optional
negatives : `lsst.afw.table.SourceCatalog`, optional
Negative polarity footprints.

Returns
Expand All @@ -467,11 +471,22 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
"""
self.metadata["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)
# preserve peak schema, if there are any footprints
if len(positives) > 0:
schema = positives[0].getFootprint().peaks.schema
elif len(negatives) > 0:
schema = negatives[0].getFootprint().peaks.schema
else:
schema = afwDetection.PeakTable.makeMinimalSchema()
mergeList = lsst.afw.detection.FootprintMergeList(afwTable.SourceTable.makeMinimalSchema(),
["positive", "negative"], schema)
initialDiaSources = lsst.afw.table.SourceCatalog(positives.schema)
# Start with positive, as FootprintMergeList will self-merge the
# subsequent added catalogs, and we want to try to preserve
# deblended positive sources.
mergeList.addCatalog(initialDiaSources.table, positives, "positive", minNewPeakDist=0)
mergeList.addCatalog(initialDiaSources.table, negatives, "negative", minNewPeakDist=0)
mergeList.getFinalSources(initialDiaSources)
self.log.info("Merging detections into %d sources", len(initialDiaSources))
else:
initialDiaSources = sources
Expand Down Expand Up @@ -528,14 +543,10 @@ def _deblend(self, difference, positiveFootprints, negativeFootprints):
-------
sources : `lsst.afw.table.SourceCatalog`
Positive and negative deblended children.
positives, negatives : `lsst.afw.detection.FootprintSet`
Deblended positive and negative polarity footprints measured on
``difference``.
positives, negatives : `lsst.afw.table.SourceCatalog`
Deblended positive and negative polarity sources with footprints
detected on ``difference``.
"""
def makeFootprints(sources):
footprints = afwDetection.FootprintSet(difference.getBBox())
footprints.setFootprints([src.getFootprint() for src in sources])
return footprints

def deblend(footprints, negative=False):
"""Deblend a positive or negative footprint set,
Expand Down Expand Up @@ -580,7 +591,7 @@ def deblend(footprints, negative=False):
sources.reserve(len(positives) + len(negatives))
sources.extend(positives, deep=True)
sources.extend(negatives, deep=True)
return sources, makeFootprints(positives), makeFootprints(negatives)
return sources, positives, negatives

def _removeBadSources(self, diaSources):
"""Remove unphysical diaSources from the catalog.
Expand Down Expand Up @@ -877,11 +888,15 @@ def run(self, science, matchedTemplate, difference, scoreExposure,

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

else:
positives = afwTable.SourceCatalog(self.schema)
results.positive.makeSources(positives)
negatives = afwTable.SourceCatalog(self.schema)
results.negative.makeSources(negatives)
return self.processResults(science, matchedTemplate, difference,
results.sources, idFactory,
positiveFootprints=results.positive,
negativeFootprints=results.negative)
positives=positives,
negatives=negatives)
102 changes: 93 additions & 9 deletions tests/test_detectAndMeasure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ def test_detection_xy0(self):
difference = science.clone()

# Configure the detection Task
detectionTask = self._setup_detection()
detectionTask = self._setup_detection(doDeblend=True)

# Run detection and check the results
output = detectionTask.run(science, matchedTemplate, difference,
Expand All @@ -166,6 +166,12 @@ def test_detection_xy0(self):
self.assertTrue(all(output.diaSources['id'] > 1000000000))
self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)

# all of the sources should have been detected
self.assertEqual(len(output.diaSources), len(sources))
refIds = []
for source in sources:
self._check_diaSource(output.diaSources, source, refIds=refIds)

def test_raise_bad_psf(self):
"""Detection should raise if the PSF width is NaN
"""
Expand Down Expand Up @@ -647,7 +653,7 @@ def test_detection_xy0(self):
score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)

# Configure the detection Task
detectionTask = self._setup_detection()
detectionTask = self._setup_detection(doDeblend=True)

# Run detection and check the results
output = detectionTask.run(science, matchedTemplate, difference, score,
Expand All @@ -659,6 +665,16 @@ def test_detection_xy0(self):

self.assertImagesEqual(subtractedMeasuredExposure.image, difference.image)

# Not all of the sources will be detected: preconvolution results in
# a larger edge mask, so we miss an edge source.
self.assertEqual(len(output.diaSources), len(sources)-1)
# TODO DM-41496: restore this block once we handle detections on edge
# pixels better; at least one of these sources currently has a bad
# centroid because most of the source is rejected as EDGE.
# refIds = []
# for source in sources:
# self._check_diaSource(output.diaSources, source, refIds=refIds)

def test_measurements_finite(self):
"""Measured fluxes and centroids should always be finite.
"""
Expand Down Expand Up @@ -944,11 +960,10 @@ def testDeblendNegatives(self):

# DM-48704 fixed a problem where the peaks were in the footprints, but
# the spans were empty.
footprints = negatives.getFootprints()
self.assertEqual(len(negatives.getFootprints()), 2)
self.assertEqual(len(positives.getFootprints()), 0)
self.assertGreater(footprints[0].getSpans().getArea(), 0)
self.assertGreater(footprints[1].getSpans().getArea(), 0)
self.assertEqual(len(negatives), 2)
self.assertEqual(len(positives), 0)
self.assertGreater(negatives[0].getFootprint().getSpans().getArea(), 0)
self.assertGreater(negatives[1].getFootprint().getSpans().getArea(), 0)
# Deblended children are HeavyFootprints; we have to make sure the
# pixel values in those are correct (though DetectAndMeasureTask
# doesn't use the fact that they're Heavy).
Expand All @@ -957,11 +972,80 @@ def testDeblendNegatives(self):
# the maximum value in the science catalog footprint (ignoring the
# noise in the science and template, hence rtol).
# (catalog[0] is the parent; we want the children)
self.assertFloatsAlmostEqual(footprints[0].getImageArray().min(),
self.assertFloatsAlmostEqual(negatives[0].getFootprint().getImageArray().min(),
-catalog[1].getFootprint().getImageArray().max(), rtol=1e-4)
self.assertFloatsAlmostEqual(footprints[1].getImageArray().min(),
self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(),
-catalog[2].getFootprint().getImageArray().max(), rtol=1e-4)

def testMergeFootprints(self):
"""Test that merging footprints does not cause negative ones to
disappear (e.g. get merged into non-connected footprints).
"""
# Make a science image with multiple blends, designed to trigger the
# negative-footprint-related bug in `FootprintSet.merge`.
bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(400, 200))
dataset = lsst.meas.base.tests.TestDataset(bbox)
dataset.addSource(instFlux=.5E5, centroid=lsst.geom.Point2D(25, 26))
dataset.addSource(instFlux=.7E5, centroid=lsst.geom.Point2D(75, 24),
shape=lsst.afw.geom.Quadrupole(12, 7, 2))
delta = 10
with dataset.addBlend() as family:
family.addChild(instFlux=1E5, centroid=lsst.geom.Point2D(150, 72))
family.addChild(instFlux=1.5E5, centroid=lsst.geom.Point2D(150+delta, 74))
with dataset.addBlend() as family:
family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(250, 72))
family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(250+delta, 74))
with dataset.addBlend() as family:
family.addChild(instFlux=3E5, centroid=lsst.geom.Point2D(350, 72))
family.addChild(instFlux=3.5E5, centroid=lsst.geom.Point2D(350+delta, 74))
dataset.addSource(instFlux=4E5, centroid=lsst.geom.Point2D(75, 124))
dataset.addSource(instFlux=5E5, centroid=lsst.geom.Point2D(175, 124))
template, catalog = dataset.realize(noise=100.0,
schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
# These will be positive sources on the diffim.
template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(15, 15),
lsst.geom.Point2I(35, 35))).array *= -1
template.image.subset(lsst.geom.Box2I(lsst.geom.Point2I(230, 60),
lsst.geom.Point2I(270, 85))).array *= -1
dataset = lsst.meas.base.tests.TestDataset(bbox)
science, _ = dataset.realize(noise=100.0,
schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema())
difference = science.clone()
difference.image -= template.image

config = detectAndMeasure.DetectAndMeasureTask.ConfigClass()
config.doDeblend = True
task = detectAndMeasure.DetectAndMeasureTask(config=config)
result = task.run(science, template, difference)

# The original bug manifested as the (175,124) source disappaering in
# the merge step, and being included in the peaks of a duplicated
# (75,124) source, even though their footprints are not connected.
mask = np.isclose(result.diaSources["slot_Centroid_x"], 75) & \
np.isclose(result.diaSources["slot_Centroid_y"], 124)
self.assertEqual(mask.sum(), 1)
peaks = result.diaSources[mask][0].getFootprint().peaks
self.assertEqual(len(peaks), 1)
self.assertEqual(peaks[0].getIx(), 75)
self.assertEqual(peaks[0].getIy(), 124)

mask = np.isclose(result.diaSources["slot_Centroid_x"], 175) & \
np.isclose(result.diaSources["slot_Centroid_y"], 124)
self.assertEqual(mask.sum(), 1)
peaks = result.diaSources[mask][0].getFootprint().peaks
self.assertEqual(len(peaks), 1)
self.assertEqual(peaks[0].getIx(), 175)
self.assertEqual(peaks[0].getIy(), 124)
Comment on lines +1021 to +1038
Copy link
Contributor

@fred3m fred3m Feb 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure that this is the most useful test anymore. While it did catch the one specific bug in mergeFootprintSets it doesn't check for the general problem of having duplicate peaks. I would suggest something like

peak_x = np.array([peak['f_x'] for src in catalog[catalog["parent"]==0] for peak in src.getFootprint().peaks])
peak_y = np.array([peak['f_y'] for src in catalog[catalog["parent"]==0] for peak in src.getFootprint().peaks])
for i in range(len(peak_x)-1):
    for j in range(i, len(peak_x)):
        self.assertFalse((peak_x[i] == peak_x[j]) & (peak_y[i] == peak_y[j]))

Note: this code is untested and there is probably a way to vectorize it using np.unique. Something like (also untested)

coords = np.column_stack((peak_x, peak_y))
unique_coords, counts = np.unique(coords, axis=0, return_counts=True)
self.assertEqual(np.sum(counts > 1), 0)

Copy link
Contributor Author

@parejkoj parejkoj Feb 14, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Fred, that was a nearly correct bit of code you produced there! I've added that and checked that it does fail (with one count==2) for the original bug condition, but left the original tests with a comment about why.


# This checks that all the returned footprints have exactly one peak;
# it's a more general test than the above, but I'm keeping that for
# the more explicit test of the original bugged sources.
peak_x = np.array([peak['f_x'] for src in result.diaSources for peak in src.getFootprint().peaks])
peak_y = np.array([peak['f_y'] for src in result.diaSources for peak in src.getFootprint().peaks])
peaks = np.column_stack((peak_x, peak_y))
unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True)
self.assertEqual(np.sum(counts > 1), 0)


def setup_module(module):
lsst.utils.tests.init()
Expand Down