From 159911c65d6188843776ab66973006e63c9ceb9f Mon Sep 17 00:00:00 2001 From: John Parejko Date: Wed, 12 Feb 2025 03:00:27 -0800 Subject: [PATCH 1/2] Add basic tests that purely positive sources are deblended This is a change I somehow didn't make on DM-42696, but should have, to demonstrate that we really are deblending in the "easy" case. That ticket made it so that we could test against a proper truth catalog, but never implemented said test. --- tests/test_detectAndMeasure.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 31a7351ad..4b8c58e37 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -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, @@ -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 """ @@ -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, @@ -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. """ From 2cf1bc213441d069bb76498d31045a0707304257 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Wed, 5 Feb 2025 15:43:26 -0800 Subject: [PATCH 2/2] Use FootprintMergeList instead of FootprintSet.merge This avoids the (unresolved) bug in merge(), and leaves the positive footprints unmerged in the output catalog. It does require an API change to processResults, since we now have to use SourceCatalogs instead of FootprintSets, but I think that's an overall improvement (though it does make the `if doDeblend` branches rather messier). --- python/lsst/ip/diffim/detectAndMeasure.py | 63 ++++++++++------- tests/test_detectAndMeasure.py | 82 +++++++++++++++++++++-- 2 files changed, 114 insertions(+), 31 deletions(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 2c2261c49..7dd227ed9 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -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. @@ -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 @@ -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 @@ -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 @@ -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, @@ -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. @@ -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) diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 4b8c58e37..8812a0e9b 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -960,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). @@ -973,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) + + # 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()