From 6790bf12432e631fdcbc60b083cccec1bd0a4b64 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 13 Jul 2026 22:00:18 +0200 Subject: [PATCH 01/12] Add peak clustering within difference image footprints Cluster detected peaks by single-linkage on their positions, ignoring polarity so dipole lobes group together. Isolated positive peaks can then be extracted as their own diaSources. --- python/lsst/ip/diffim/peakClustering.py | 260 ++++++++++++++++++++++++ tests/test_peakClustering.py | 244 ++++++++++++++++++++++ 2 files changed, 504 insertions(+) create mode 100644 python/lsst/ip/diffim/peakClustering.py create mode 100644 tests/test_peakClustering.py diff --git a/python/lsst/ip/diffim/peakClustering.py b/python/lsst/ip/diffim/peakClustering.py new file mode 100644 index 00000000..583a6fca --- /dev/null +++ b/python/lsst/ip/diffim/peakClustering.py @@ -0,0 +1,260 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Identify spatial clusters of peaks within difference image footprints. +""" + +__all__ = ["PeakCluster", "FootprintPeakClusters", "cluster_peaks", + "find_peak_clusters"] + +import dataclasses + +import numpy as np +from scipy.sparse import coo_matrix +from scipy.sparse.csgraph import connected_components +from scipy.spatial import cKDTree + +import lsst.afw.detection +import lsst.afw.table + + +# Per-peak flag set by `FootprintMergeList` marking an individual peak as +# originating from the negative-polarity detection footprint set. +_NEGATIVE_PEAK_FLAG = "merge_peak_negative" + + +@dataclasses.dataclass +class PeakCluster: + """A group of peaks within a single footprint that lie close together. + + Peaks are grouped by single-linkage clustering on their positions, + ignoring polarity, so a dipole's positive and negative lobes should fall + into the same cluster. A cluster with a single member is "isolated": no other + peak in the footprint lies within the linking radius. + + Parameters + ---------- + peak_indices : `list` [`int`] + Rows into the parent footprint's peak catalog + (``footprint.getPeaks()``) that belong to this cluster. + positions : `numpy.ndarray` + ``(n, 2)`` array of integer peak-pixel positions ``(i_x, i_y)`` in + the parent (image) coordinate frame. + peak_values : `numpy.ndarray` + ``(n,)`` array of the signed peak values. + is_negative : `numpy.ndarray` + ``(n,)`` boolean array, `True` for peaks from the negative detection + set (from the ``merge_peak_negative`` flag). + """ + + peak_indices: list + positions: np.ndarray + peak_values: np.ndarray + is_negative: np.ndarray + + @property + def is_singleton(self): + """Whether this cluster holds exactly one peak (`bool`).""" + return len(self.peak_indices) == 1 + + @property + def is_isolated_positive(self): + """Whether this cluster is a single positive peak (`bool`). + + These are the peaks the deblender extracts as new footprints. + """ + return self.is_singleton and not bool(self.is_negative[0]) + + +@dataclasses.dataclass +class FootprintPeakClusters: + """The peak clusters found within one multi-peak footprint. + + Parameters + ---------- + source : `lsst.afw.table.SourceRecord` + The source record whose footprint was clustered. This is the live + record from the catalog, so a downstream deblender can modify it (and + its footprint) in place without a separate lookup. + clusters : `list` [`PeakCluster`] + The clusters of peaks within the footprint. + """ + + source: lsst.afw.table.SourceRecord + clusters: list + + @property + def footprint(self): + """The footprint whose peaks were clustered + (`lsst.afw.detection.Footprint`). + """ + return self.source.getFootprint() + + @property + def source_id(self): + """Id of the source record (`int`). + + Only meaningful once ids have been assigned to the catalog; while + clustering runs inside the merge step, ids are not yet set. + """ + return self.source.getId() + + @property + def n_peaks(self): + """Total number of peaks in the footprint (`int`).""" + return sum(len(cluster.peak_indices) for cluster in self.clusters) + + @property + def isolated_positive(self): + """Single-peak positive clusters (`list` [`PeakCluster`]). + + The candidates for extracting into their own footprints. + """ + return [cluster for cluster in self.clusters if cluster.is_isolated_positive] + + @property + def retained(self): + """Clusters whose pixels stay in the footprint (`list` [`PeakCluster`]). + + Everything that is not an isolated positive peak: multi-peak blends, + dipoles, and isolated negative peaks. + """ + return [cluster for cluster in self.clusters if not cluster.is_isolated_positive] + + @property + def is_deblendable(self): + """Whether the footprint can be split by cluster deblending. + """ + return bool(self.isolated_positive) and len(self.clusters) >= 2 + + +def _single_linkage_labels(positions, linking_radius): + """Assign a cluster label to each position by single-linkage clustering. + + Two positions are linked when they lie within ``linking_radius`` of each + other; clusters are the connected components of the resulting graph. + + Parameters + ---------- + positions : `numpy.ndarray` + ``(n, 2)`` array of positions. + linking_radius : `float` + Maximum distance, in pixels, between two peaks for them to be linked + into the same cluster. + + Returns + ------- + labels : `numpy.ndarray` + ``(n,)`` integer array of cluster labels. + """ + n = len(positions) + if n <= 1: + return np.zeros(n, dtype=int) + pairs = cKDTree(positions).query_pairs(r=linking_radius, output_type="ndarray") + if len(pairs) == 0: + graph = coo_matrix((n, n)) + else: + graph = coo_matrix((np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n)) + _, labels = connected_components(graph, directed=False) + return labels + + +def cluster_peaks(peaks, linking_radius): + """Group the peaks of one footprint into spatial clusters. + + Peaks are clustered by single-linkage on their positions, ignoring + polarity: two peaks within ``linking_radius`` of each other join the same + cluster, and clustering is transitive (chaining). A peak with no + neighbour within the radius forms its own singleton cluster. + + Parameters + ---------- + peaks : `lsst.afw.detection.PeakCatalog` + Peaks of a single (merged) footprint. The schema must carry the + ``merge_peak_negative`` flag. + linking_radius : `float` + Maximum distance, in pixels, between two peaks for them to be + clustered together. + + Returns + ------- + clusters : `list` [`PeakCluster`] + One entry per cluster. Empty if ``peaks`` is empty. + """ + n = len(peaks) + positions = np.empty((n, 2), dtype=float) + peak_values = np.empty(n, dtype=float) + is_negative = np.empty(n, dtype=bool) + negative_key = peaks.getSchema().find(_NEGATIVE_PEAK_FLAG).key + for i, peak in enumerate(peaks): + # Detection peaks store the integer local-maximum pixel; f_x/f_y equal + # i_x/i_y here, so use the integer accessors to be explicit. + positions[i] = (peak.getIx(), peak.getIy()) + peak_values[i] = peak.getPeakValue() + is_negative[i] = peak.get(negative_key) + + labels = _single_linkage_labels(positions, linking_radius) + clusters = [] + for label in np.unique(labels): + (members,) = np.nonzero(labels == label) + clusters.append(PeakCluster( + peak_indices=members.tolist(), + positions=positions[members], + peak_values=peak_values[members], + is_negative=is_negative[members], + )) + return clusters + + +def find_peak_clusters(catalog, linking_radius, *, min_peaks=2): + """Find peak clusters in every multi-peak footprint of a catalog. + + Parameters + ---------- + catalog : `lsst.afw.table.SourceCatalog` + Sources whose footprints are to be searched, e.g. the merged + ``initialDiaSources``. Footprint peak catalogs must carry the + ``merge_peak_negative`` flag. + linking_radius : `float` + Maximum distance, in pixels, between two peaks for them to be + clustered together. + min_peaks : `int`, optional + Only footprints with at least this many peaks are considered; those + with fewer are not deblending candidates and are skipped. + + Returns + ------- + results : `list` [`FootprintPeakClusters`] + One entry per footprint with at least ``min_peaks`` peaks. + """ + results = [] + for source in catalog: + footprint = source.getFootprint() + if footprint is None: + continue + peaks = footprint.getPeaks() + if len(peaks) < min_peaks: + continue + results.append(FootprintPeakClusters( + source=source, + clusters=cluster_peaks(peaks, linking_radius), + )) + return results diff --git a/tests/test_peakClustering.py b/tests/test_peakClustering.py new file mode 100644 index 00000000..c88261f1 --- /dev/null +++ b/tests/test_peakClustering.py @@ -0,0 +1,244 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest + +import lsst.afw.detection as afwDetection +import lsst.afw.geom as afwGeom +import lsst.afw.table as afwTable +import lsst.geom as geom +import lsst.utils.tests +from lsst.ip.diffim.peakClustering import cluster_peaks, find_peak_clusters + + +def makePeakSchema(): + """Build a source schema and the matching merged-peak schema. + + The peak schema carries the ``merge_peak_positive``/``merge_peak_negative`` + flags, mimicking the footprints produced by ``FootprintMergeList`` in + ``DetectAndMeasureTask``. + """ + schema = afwTable.SourceTable.makeMinimalSchema() + mergeList = afwDetection.FootprintMergeList(schema, ["positive", "negative"]) + return schema, mergeList.getPeakSchema() + + +def makeFootprint(peakSchema, peaks, size=100): + """Make a footprint spanning a ``size`` x ``size`` box with the given peaks. + + Parameters + ---------- + peakSchema : `lsst.afw.table.Schema` + Schema carrying the merge-peak flags. + peaks : `list` [`tuple`] + Each entry is ``(fx, fy, value, isNegative)``. + """ + spans = afwGeom.SpanSet(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(size, size))) + footprint = afwDetection.Footprint(spans, peakSchema) + negativeKey = peakSchema.find("merge_peak_negative").key + positiveKey = peakSchema.find("merge_peak_positive").key + for fx, fy, value, isNegative in peaks: + peak = footprint.addPeak(fx, fy, value) + peak.set(negativeKey, isNegative) + peak.set(positiveKey, not isNegative) + return footprint + + +def clusterContaining(clusters, peakIndex): + """Return the single cluster whose members include ``peakIndex``.""" + matches = [c for c in clusters if peakIndex in c.peak_indices] + assert len(matches) == 1, f"peak {peakIndex} found in {len(matches)} clusters" + return matches[0] + + +class PeakClusteringTestCase(lsst.utils.tests.TestCase): + + def setUp(self): + self.schema, self.peakSchema = makePeakSchema() + self.radius = 5.0 + + def test_two_peaks_far_apart(self): + """Peaks beyond the linking radius form separate singleton clusters.""" + footprint = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (30, 30, 120, False)]) + clusters = cluster_peaks(footprint.getPeaks(), self.radius) + self.assertEqual(len(clusters), 2) + for cluster in clusters: + self.assertTrue(cluster.is_singleton) + self.assertTrue(cluster.is_isolated_positive) + + def test_two_peaks_close(self): + """Peaks within the linking radius merge into one cluster.""" + footprint = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (12, 10, 90, False)]) + clusters = cluster_peaks(footprint.getPeaks(), self.radius) + self.assertEqual(len(clusters), 1) + self.assertFalse(clusters[0].is_singleton) + self.assertFalse(clusters[0].is_isolated_positive) + self.assertEqual(sorted(clusters[0].peak_indices), [0, 1]) + + def test_transitive_chain(self): + """Single-linkage chains A-B-C even when A and C exceed the radius.""" + # 10->14->18: adjacent gaps are 4 (< 5), but 10->18 is 8 (> 5). + footprint = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (14, 10, 100, False), + (18, 10, 100, False)]) + clusters = cluster_peaks(footprint.getPeaks(), self.radius) + self.assertEqual(len(clusters), 1) + self.assertEqual(sorted(clusters[0].peak_indices), [0, 1, 2]) + + def test_dipole_single_cluster(self): + """A close positive/negative pair clusters together, not isolated.""" + footprint = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (12, 10, -90, True)]) + clusters = cluster_peaks(footprint.getPeaks(), self.radius) + self.assertEqual(len(clusters), 1) + cluster = clusters[0] + self.assertFalse(cluster.is_isolated_positive) + self.assertEqual(cluster.is_negative.tolist(), [False, True]) + self.assertEqual(cluster.peak_values.tolist(), [100, -90]) + + def test_isolated_positive_beside_dipole(self): + """An isolated positive peak is separated from a nearby dipole.""" + footprint = makeFootprint(self.peakSchema, + [(50, 50, 200, False), # isolated positive + (10, 10, 100, False), (12, 10, -90, True)]) # dipole + clusters = cluster_peaks(footprint.getPeaks(), self.radius) + self.assertEqual(len(clusters), 2) + isolated = clusterContaining(clusters, 0) + dipole = clusterContaining(clusters, 1) + self.assertTrue(isolated.is_isolated_positive) + self.assertFalse(dipole.is_isolated_positive) + self.assertEqual(sorted(dipole.peak_indices), [1, 2]) + + def test_isolated_negative_not_extracted(self): + """A lone negative peak is a singleton but not an isolated positive.""" + footprint = makeFootprint(self.peakSchema, + [(50, 50, 200, False), (10, 10, -80, True)]) + clusters = cluster_peaks(footprint.getPeaks(), self.radius) + negative = clusterContaining(clusters, 1) + self.assertTrue(negative.is_singleton) + self.assertFalse(negative.is_isolated_positive) + + def test_positions_and_values(self): + """Cluster arrays carry the peaks' positions and signed values.""" + footprint = makeFootprint(self.peakSchema, [(10, 20, 100, False)]) + (cluster,) = cluster_peaks(footprint.getPeaks(), self.radius) + self.assertEqual(cluster.positions.tolist(), [[10, 20]]) + self.assertEqual(cluster.peak_values.tolist(), [100]) + self.assertEqual(cluster.is_negative.tolist(), [False]) + + +class FindPeakClustersTestCase(lsst.utils.tests.TestCase): + + def setUp(self): + self.schema, self.peakSchema = makePeakSchema() + self.radius = 5.0 + + def makeCatalog(self, footprints): + catalog = afwTable.SourceCatalog(self.schema) + for i, footprint in enumerate(footprints): + record = catalog.addNew() + record.setId(i + 1) + record.setFootprint(footprint) + return catalog + + def test_skips_single_peak_footprints(self): + """Footprints with fewer than min_peaks peaks are not returned.""" + single = makeFootprint(self.peakSchema, [(10, 10, 100, False)]) + multi = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (30, 30, 120, False)]) + catalog = self.makeCatalog([single, multi]) + results = find_peak_clusters(catalog, self.radius) + self.assertEqual(len(results), 1) + # Only the second (id=2) footprint qualifies. + self.assertEqual(results[0].source_id, 2) + self.assertEqual(results[0].n_peaks, 2) + + def test_handles_missing_footprint(self): + """Sources with no footprint are skipped rather than raising.""" + multi = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (30, 30, 120, False)]) + catalog = afwTable.SourceCatalog(self.schema) + catalog.addNew().setId(1) # no footprint attached + record = catalog.addNew() + record.setId(2) + record.setFootprint(multi) + results = find_peak_clusters(catalog, self.radius) + self.assertEqual(len(results), 1) + self.assertEqual(results[0].source_id, 2) + + def test_deblendable_classification(self): + """is_deblendable requires both an extraction candidate and a + remainder. + """ + # Isolated positive + dipole: deblendable. + mixed = makeFootprint(self.peakSchema, + [(50, 50, 200, False), + (10, 10, 100, False), (12, 10, -90, True)]) + # Two isolated positives beyond the radius: deblendable into two + # sources even though there is no remainder cluster. + allIsolated = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (50, 50, 120, False)]) + # One tight blend: a remainder but nothing isolated to extract. + blend = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (12, 10, 90, False)]) + catalog = self.makeCatalog([mixed, allIsolated, blend]) + results = find_peak_clusters(catalog, self.radius) + byId = {r.source_id: r for r in results} + + self.assertTrue(byId[1].is_deblendable) + self.assertEqual(len(byId[1].isolated_positive), 1) + self.assertEqual(len(byId[1].retained), 1) + + self.assertTrue(byId[2].is_deblendable) + self.assertEqual(len(byId[2].isolated_positive), 2) + self.assertEqual(len(byId[2].retained), 0) + + self.assertFalse(byId[3].is_deblendable) + self.assertEqual(len(byId[3].isolated_positive), 0) + self.assertEqual(len(byId[3].retained), 1) + + def test_source_handle_is_live(self): + """The stored source is the live catalog record, so deblender edits + made through it land on the catalog. + """ + multi = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (30, 30, 120, False)]) + catalog = self.makeCatalog([multi]) + results = find_peak_clusters(catalog, self.radius) + # Mutating through the stored handle must be visible in the catalog. + results[0].source.setId(4242) + self.assertEqual(catalog[0].getId(), 4242) + self.assertEqual(results[0].footprint.getBBox(), catalog[0].getFootprint().getBBox()) + + +class MemoryTester(lsst.utils.tests.MemoryTestCase): + pass + + +def setup_module(module): + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main() From a1c55329df7c50a7538a344fe51ec2b72e5879cf Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 13 Jul 2026 22:00:25 +0200 Subject: [PATCH 02/12] Add cluster-based footprint deblender for difference images Extract isolated positive peaks out of multi-peak footprints into new footprints, leaving the remainder on the original source. --- python/lsst/ip/diffim/clusterDeblend.py | 251 ++++++++++++++++++++++++ tests/test_clusterDeblend.py | 195 ++++++++++++++++++ 2 files changed, 446 insertions(+) create mode 100644 python/lsst/ip/diffim/clusterDeblend.py create mode 100644 tests/test_clusterDeblend.py diff --git a/python/lsst/ip/diffim/clusterDeblend.py b/python/lsst/ip/diffim/clusterDeblend.py new file mode 100644 index 00000000..2fa73e82 --- /dev/null +++ b/python/lsst/ip/diffim/clusterDeblend.py @@ -0,0 +1,251 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Simple footprint deblender for difference images. +""" + +__all__ = ["DeblendedPeaks", "deblend_clustered_peaks"] + +import dataclasses + +import numpy as np + +import lsst.afw.detection +import lsst.afw.geom +import lsst.afw.image + + +# Per-peak flag marking an individual peak as originating from the negative +# detection set (see peakClustering._NEGATIVE_PEAK_FLAG). +_NEGATIVE_PEAK_FLAG = "merge_peak_negative" + + +@dataclasses.dataclass +class DeblendedPeaks: + """The footprints resulting from deblending one multi-peak footprint. + + Parameters + ---------- + parent : `lsst.afw.detection.Footprint` + Footprint to set back on the original source record. This is the + remainder (all non-carved clusters) if there is one; otherwise, when + the footprint is entirely single positive peaks, it is one of those + carved peaks. + children : `list` [`lsst.afw.detection.Footprint`] + One new footprint per remaining carved single positive peak, each + carrying that single peak. + """ + + parent: lsst.afw.detection.Footprint + children: list + + +def _peak_amplitude(image, mask, col, row, bad_bitmask): + """Summed absolute image flux in the 3x3 box around a peak. + + Parameters + ---------- + image, mask : `numpy.ndarray` + Difference-image and mask arrays. + col, row : `int` + Peak location in array coordinates (i.e. image coordinates minus the + image origin). + bad_bitmask : `int` + Bitmask of mask planes to exclude from the sum. + + Returns + ------- + amplitude : `float` + """ + height, width = image.shape + c0, c1 = max(col - 1, 0), min(col + 2, width) + r0, r1 = max(row - 1, 0), min(row + 2, height) + sub_image = image[r0:r1, c0:c1] + good = (mask[r0:r1, c0:c1] & bad_bitmask) == 0 + return float(np.sum(np.abs(sub_image[good]))) + + +def _merge_peak_clusters(peak_clusters, weighted, n_peaks): + """Merge clusters of peaks until every peak owns (dominates) its own peak + pixel. + + Parameters + ---------- + peaks : `list` [`list` [`int`]] + Initial peaks, each a list of peak indices. + weighted : `numpy.ndarray` + ``(n_peaks, n_peaks)`` array where ``weighted[p, j]`` is peak ``p``'s + influence evaluated at peak ``j``'s pixel. + n_peaks : `int` + Number of peaks. + + Returns + ------- + peaks : `list` [`list` [`int`]] + The merged peaks. + """ + peak_clusters = [list(peak_cluster) for peak_cluster in peak_clusters] + while True: + owner = [0]*n_peaks + for index, peak_cluster in enumerate(peak_clusters): + for peak in peak_cluster: + owner[peak] = index + merged = False + for j in range(n_peaks): + # Check if the current peak is the most significant influence at + # its location. If not, add it to the more significant cluster. + influence = np.array([weighted[peak_cluster, j].sum() for peak_cluster in peak_clusters]) + best = int(np.argmax(influence)) + # Merge only if another cluster of peaks *strictly* dominates the + # owner at this peak pixel. + if influence[best] > influence[owner[j]]: + low, high = sorted((owner[j], best)) + peak_clusters[low] = peak_clusters[low] + peak_clusters[high] + del peak_clusters[high] + merged = True + break + if not merged: + return peak_clusters + + +def _spans_to_bool(span_set, bbox): + """Rasterize a SpanSet to a boolean array over ``bbox``.""" + array = np.zeros((bbox.getHeight(), bbox.getWidth()), dtype=bool) + min_x, min_y = bbox.getMinX(), bbox.getMinY() + for span in span_set: + array[span.getY() - min_y, span.getX0() - min_x:span.getX1() - min_x + 1] = True + return array + + +def _footprint_from_bool(pixels, bbox, peak_schema, region): + """Build a peak-schema-preserving footprint from a boolean pixel array.""" + mask = lsst.afw.image.Mask(bbox) + bit = mask.getPlaneBitMask("DETECTED") + mask.array[:] = 0 + mask.array[pixels] = bit + footprint = lsst.afw.detection.Footprint(lsst.afw.geom.SpanSet.fromMask(mask, bit), peak_schema) + footprint.setRegion(region) + return footprint + + +def deblend_clustered_peaks(footprint_clusters, image, mask, image_xy0, sigma, bad_bitmask): + """Extract the isolated positive peaks out of one footprint. + + Parameters + ---------- + footprint_clusters : `~lsst.ip.diffim.peakClustering.FootprintPeakClusters` + Clustering result for one multi-peak footprint. + image, mask : `numpy.ndarray` + Difference-image and mask arrays used to weight the peaks. + image_xy0 : `tuple` [`int`, `int`] + Image origin ``(x0, y0)``, so peak image coordinates can index the + arrays. + sigma : `float` + Width (in pixels) of the 2D Gaussian to use to determine the relative + significance of each peak to the pixels in the footprint. + bad_bitmask : `int` + Bitmask of mask planes to exclude from the peak amplitude sums. + + Returns + ------- + result : `DeblendedPeaks` or `None` + The parent footprint (kept on the original record) and the new + footprints, or `None` if the footprint should be left unmodified. + """ + footprint = footprint_clusters.footprint + peaks = footprint.getPeaks() + n_peaks = len(peaks) + peak_schema = peaks.getSchema() + negative_key = peak_schema.find(_NEGATIVE_PEAK_FLAG).key + x0, y0 = image_xy0 + + positions = np.zeros((n_peaks, 2), dtype=int) + amplitudes = np.zeros(n_peaks, dtype=float) + is_negative = np.zeros(n_peaks, dtype=bool) + for i, peak in enumerate(peaks): + ix, iy = peak.getIx(), peak.getIy() + positions[i] = (ix, iy) + is_negative[i] = peak.get(negative_key) + amplitudes[i] = _peak_amplitude(image, mask, ix - x0, iy - y0, bad_bitmask) + + # weighted[p, j]: influence of peak p evaluated at peak j's pixel. + inv_two_sigma_sq = 1.0/(2.0*sigma*sigma) + delta = positions[:, np.newaxis, :] - positions[np.newaxis, :, :] + dist_sq = np.sum(delta*delta, axis=2) + weighted = amplitudes[:, np.newaxis]*np.exp(-dist_sq*inv_two_sigma_sq) + + peak_clusters = _merge_peak_clusters( + [cluster.peak_indices for cluster in footprint_clusters.clusters], weighted, n_peaks) + + def is_isolated_positive(peak_cluster): + return len(peak_cluster) == 1 and not is_negative[peak_cluster[0]] + + # Extract each single positive peak into its own source. + # Every other cluster pools into one remainder. + extracted_peaks = [peak_cluster for peak_cluster in peak_clusters if is_isolated_positive(peak_cluster)] + remainder_peak_clusters = [peak_cluster for peak_cluster in peak_clusters + if not is_isolated_positive(peak_cluster)] + n_sources = len(extracted_peaks) + (1 if remainder_peak_clusters else 0) + if not extracted_peaks or n_sources < 2: + return None + + # Assign every footprint pixel to the cluster of highest influence. + bbox = footprint.getBBox() + min_x, min_y = bbox.getMinX(), bbox.getMinY() + footprint_mask = _spans_to_bool(footprint.getSpans(), bbox) + xx, yy = np.meshgrid(np.arange(min_x, bbox.getMaxX() + 1), + np.arange(min_y, bbox.getMaxY() + 1)) + ordered = extracted_peaks + remainder_peak_clusters + influence = np.empty((len(ordered), bbox.getHeight(), bbox.getWidth()), dtype=float) + for index, peak_cluster in enumerate(ordered): + accumulator = np.zeros(footprint_mask.shape, dtype=float) + for peak in peak_cluster: + r_sq = (xx - positions[peak, 0])**2 + (yy - positions[peak, 1])**2 + accumulator += amplitudes[peak]*np.exp(-r_sq*inv_two_sigma_sq) + influence[index] = accumulator + labels = np.argmax(influence, axis=0) + # Guarantee each cluster owns its own peak pixels (breaks any influence + # ties in the owner's favour, so no cluster ends up empty). + for index, peak_cluster in enumerate(ordered): + for peak in peak_cluster: + labels[positions[peak, 1] - min_y, positions[peak, 0] - min_x] = index + + region = footprint.getRegion() + # One footprint per extracted single positive peak. + extracted_footprints = [] + for index, peak_cluster in enumerate(extracted_peaks): + fp = _footprint_from_bool(footprint_mask & (labels == index), bbox, peak_schema, region) + fp.getPeaks().append(peaks[peak_cluster[0]]) + extracted_footprints.append(fp) + + if remainder_peak_clusters: + remainder_pixels = np.zeros(footprint_mask.shape, dtype=bool) + for index in range(len(extracted_peaks), len(ordered)): + remainder_pixels |= footprint_mask & (labels == index) + remainder = _footprint_from_bool(remainder_pixels, bbox, peak_schema, region) + for peak_cluster in remainder_peak_clusters: + for peak in peak_cluster: + remainder.getPeaks().append(peaks[peak]) + return DeblendedPeaks(parent=remainder, children=extracted_footprints) + + # No remainder cluster: the footprint is entirely single positive peaks, + # so the first stays on the original record and the rest are new records. + return DeblendedPeaks(parent=extracted_footprints[0], children=extracted_footprints[1:]) diff --git a/tests/test_clusterDeblend.py b/tests/test_clusterDeblend.py new file mode 100644 index 00000000..1be97221 --- /dev/null +++ b/tests/test_clusterDeblend.py @@ -0,0 +1,195 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest + +import numpy as np + +import lsst.afw.detection as afwDetection +import lsst.afw.geom as afwGeom +import lsst.afw.table as afwTable +import lsst.geom as geom +import lsst.utils.tests +from lsst.ip.diffim.clusterDeblend import deblend_clustered_peaks +from lsst.ip.diffim.peakClustering import find_peak_clusters + +SIZE = 40 + + +def makeSchemas(): + schema = afwTable.SourceTable.makeMinimalSchema() + mergeList = afwDetection.FootprintMergeList(schema, ["positive", "negative"]) + return schema, mergeList.getPeakSchema() + + +def makeFootprintAndImage(peakSchema, peaks, ampValue=100.0, badPeaks=()): + """Build a full-frame footprint with the given peaks, and a matching image. + + Parameters + ---------- + peaks : `list` of `tuple` + Each ``(ix, iy, isNegative)``. + ampValue : `float` + Magnitude placed in the 3x3 block around each peak (negated for + negative peaks) in the difference image. + badPeaks : `set` of `int` + Indices of peaks whose 3x3 block is flagged bad in the mask. + """ + box = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(SIZE, SIZE)) + footprint = afwDetection.Footprint(afwGeom.SpanSet(box), peakSchema) + negativeKey = peakSchema.find("merge_peak_negative").key + positiveKey = peakSchema.find("merge_peak_positive").key + + image = np.zeros((SIZE, SIZE), dtype=np.float32) + mask = np.zeros((SIZE, SIZE), dtype=np.int32) + badBit = 1 + for i, (ix, iy, isNegative) in enumerate(peaks): + peak = footprint.addPeak(ix, iy, -ampValue if isNegative else ampValue) + peak.set(negativeKey, isNegative) + peak.set(positiveKey, not isNegative) + block = (slice(iy - 1, iy + 2), slice(ix - 1, ix + 2)) + image[block] = -ampValue if isNegative else ampValue + if i in badPeaks: + mask[block] = badBit + return footprint, image, mask, badBit + + +def makeClusters(schema, footprint, linkingRadius): + catalog = afwTable.SourceCatalog(schema) + record = catalog.addNew() + record.setId(1) + record.setFootprint(footprint) + results = find_peak_clusters(catalog, linkingRadius) + return results[0] if results else None + + +class ClusterDeblendTestCase(lsst.utils.tests.TestCase): + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + + def test_extract_isolated_positive_from_dipole(self): + """An isolated positive peak is split off; the dipole stays behind.""" + # Isolated positive at (8, 20); dipole at (30, 20)/(32, 20). + footprint, image, mask, _ = makeFootprintAndImage( + self.peakSchema, [(8, 20, False), (30, 20, False), (32, 20, True)]) + fc = makeClusters(self.schema, footprint, linkingRadius=5.0) + self.assertTrue(fc.is_deblendable) + + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) + self.assertIsNotNone(result) + self.assertEqual(len(result.children), 1) + + child = result.children[0] + self.assertEqual(len(child.getPeaks()), 1) + self.assertEqual((child.getPeaks()[0].getIx(), child.getPeaks()[0].getIy()), (8, 20)) + self.assertTrue(child.getSpans().contains(geom.Point2I(8, 20))) + + # The dipole is the remainder (parent), keeping both its peaks. + self.assertEqual(len(result.parent.getPeaks()), 2) + self.assertTrue(result.parent.getSpans().contains(geom.Point2I(30, 20))) + self.assertTrue(result.parent.getSpans().contains(geom.Point2I(32, 20))) + + # The extracted and remainder pixels exactly tile the original footprint. + self.assertEqual(child.getSpans().getArea() + result.parent.getSpans().getArea(), + footprint.getSpans().getArea()) + # Peak schema (with merge flags) is preserved on the outputs. + self.assertIn("merge_peak_negative", child.getPeaks().schema.getNames()) + + def test_splits_beyond_radius_positives(self): + """Two positive peaks beyond the linking radius split into two sources + even with no remainder cluster. + """ + footprint, image, mask, _ = makeFootprintAndImage( + self.peakSchema, [(8, 20, False), (32, 20, False)]) + fc = makeClusters(self.schema, footprint, linkingRadius=5.0) + self.assertTrue(fc.is_deblendable) + + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) + self.assertIsNotNone(result) + # One peak stays on the parent, the other becomes a new source. + self.assertEqual(len(result.children), 1) + self.assertEqual(len(result.parent.getPeaks()), 1) + self.assertEqual(len(result.children[0].getPeaks()), 1) + self.assertEqual(result.parent.getSpans().getArea() + + result.children[0].getSpans().getArea(), + footprint.getSpans().getArea()) + + def test_within_radius_positive_blend_stays(self): + """Two positive peaks within the linking radius are one cluster and + are left as a single source. + """ + footprint, image, mask, _ = makeFootprintAndImage( + self.peakSchema, [(20, 20, False), (22, 20, False)]) + fc = makeClusters(self.schema, footprint, linkingRadius=5.0) + self.assertFalse(fc.is_deblendable) + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) + self.assertIsNone(result) + + def test_no_isolated_positive_returns_none(self): + """A footprint with only a dipole has nothing to extract.""" + footprint, image, mask, _ = makeFootprintAndImage( + self.peakSchema, [(20, 20, False), (22, 20, True)]) + fc = makeClusters(self.schema, footprint, linkingRadius=5.0) + self.assertFalse(fc.is_deblendable) + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) + self.assertIsNone(result) + + def test_dominated_isolated_peak_is_merged(self): + """A faint 'isolated' peak dominated by a bright neighbour is merged + back in, leaving nothing to extract. + """ + # Faint positive at (12, 20); bright dipole at (16, 20)/(18, 20). + peaks = [(12, 20, False), (16, 20, False), (18, 20, True)] + footprint, image, mask, _ = makeFootprintAndImage(self.peakSchema, peaks, ampValue=100.0) + # Make the first peak faint: overwrite its 3x3 block with a small value. + image[19:22, 11:14] = 5.0 + fc = makeClusters(self.schema, footprint, linkingRadius=3.0) + self.assertTrue(fc.is_deblendable) + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=2.0, bad_bitmask=0) + self.assertIsNone(result) + + def test_bad_masked_amplitude_excluded(self): + """A peak whose 3x3 block is entirely bad-masked gets zero amplitude + and is merged away (proving bad pixels are excluded from the sum). + """ + footprint, image, mask, badBit = makeFootprintAndImage( + self.peakSchema, [(8, 20, False), (30, 20, False), (32, 20, True)], + ampValue=100.0, badPeaks={0}) + fc = makeClusters(self.schema, footprint, linkingRadius=5.0) + self.assertTrue(fc.is_deblendable) + # With the bright isolated peak's flux excluded, it has no influence + # and cannot be extracted. + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=badBit) + self.assertIsNone(result) + + +class MemoryTester(lsst.utils.tests.MemoryTestCase): + pass + + +def setup_module(module): + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main() From f016f654ba9e81ba4e783a2fa67eb0f414ec59b7 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 13 Jul 2026 22:00:32 +0200 Subject: [PATCH 03/12] Use cluster deblending in DetectAndMeasure instead of meas.deblender Drop the SourceDeblendTask subtask and deblend merged footprints by extracting isolated positive peaks, flagging affected records is_deblended. --- python/lsst/ip/diffim/__init__.py | 2 + python/lsst/ip/diffim/detectAndMeasure.py | 206 +++++++++++----------- tests/test_detectAndMeasure.py | 66 +------ 3 files changed, 112 insertions(+), 162 deletions(-) diff --git a/python/lsst/ip/diffim/__init__.py b/python/lsst/ip/diffim/__init__.py index 58f3879c..cee6f8eb 100755 --- a/python/lsst/ip/diffim/__init__.py +++ b/python/lsst/ip/diffim/__init__.py @@ -34,6 +34,8 @@ from .dipoleMeasurement import * from .getTemplate import * from .dipoleFitTask import * +from .peakClustering import * +from .clusterDeblend import * from .detectAndMeasure import * from .imageDecorrelation import * from .imageMapReduce import * diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index ed6b9eca..8a526a55 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -39,7 +39,6 @@ from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask from lsst.meas.algorithms import FindGlintTrailsTask, FindCosmicRaysConfig, findCosmicRays from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig -import lsst.meas.deblender import lsst.meas.extensions.trailedSources # noqa: F401 import lsst.meas.extensions.shapeHSM import lsst.pex.config as pexConfig @@ -49,6 +48,8 @@ from lsst.utils.timer import timeMethod from . import DipoleFitTask +from .peakClustering import find_peak_clusters +from .clusterDeblend import deblend_clustered_peaks __all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] @@ -264,11 +265,20 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, doDeblend = pexConfig.Field( dtype=bool, default=False, - doc="Deblend DIASources after detection?" + doc="Deblend merged footprints by extracting isolated positive peaks " + "into their own diaSources? Only used if ``doMerge`` is set." ) - deblend = pexConfig.ConfigurableField( - target=lsst.meas.deblender.SourceDeblendTask, - doc="Task to split blended sources into their components." + growFootprint = pexConfig.Field( + dtype=int, + default=2, + doc="Grow positive and negative footprints by this many pixels before merging" + ) + peakClusterRadius = pexConfig.Field( + dtype=float, + default=3.0, + doc="Radius (in PSF FWHM) to cluster detected peaks within footprints. " + "Positive peaks that are isolated after clustering can be deblended " + "if doDeblend=True." ) measurement = pexConfig.ConfigurableField( target=DipoleFitTask, @@ -287,11 +297,6 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, target=ForcedMeasurementTask, doc="Task to force photometer science image at diaSource locations.", ) - growFootprint = pexConfig.Field( - dtype=int, - default=2, - doc="Grow positive and negative footprints by this many pixels before merging" - ) diaSourceMatchRadius = pexConfig.Field( dtype=float, default=0.5, @@ -522,6 +527,8 @@ class DetectAndMeasureTask(lsst.pipe.base.PipelineTask): ConfigClass = DetectAndMeasureConfig _DefaultName = "detectAndMeasure" + sigma = None + def __init__(self, **kwargs): super().__init__(**kwargs) self.schema = afwTable.SourceTable.makeMinimalSchema() @@ -531,8 +538,6 @@ def __init__(self, **kwargs): self.makeSubtask("subtractInitialBackground") self.makeSubtask("subtractFinalBackground") self.makeSubtask("detection", schema=self.schema) - if self.config.doDeblend: - self.makeSubtask("deblend", schema=self.schema) self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True) self.makeSubtask("measurement", schema=self.schema, algMetadata=self.algMetadata) @@ -594,6 +599,9 @@ def __init__(self, **kwargs): self.makeSubtask("streakDetection") self.makeSubtask("findGlints") self.schema.addField("glint_trail", "Flag", "DiaSource is part of a glint trail.") + self.schema.addField("is_deblended", "Flag", + "The original detected footprint was deblended, producing multiple diaSources." + "The deblended diaSources split the original footprint.") self.schema.addField("reliability", type="F", doc="Reliability score of the DiaSource") self.schema.addField("reliabilityVersion", type=str, size=7, doc="Version of the reliability model") @@ -786,17 +794,11 @@ def run(self, science, matchedTemplate, difference, kernelSources=None, ) measurementResults.differenceBackground = background - if self.config.doDeblend: - sources, positives, negatives = self._deblend(difference, - results.positive, - results.negative) - - else: - positives = afwTable.SourceCatalog(self.schema) - results.positive.makeSources(positives) - negatives = afwTable.SourceCatalog(self.schema) - results.negative.makeSources(negatives) - sources = results.sources + positives = afwTable.SourceCatalog(self.schema) + results.positive.makeSources(positives) + negatives = afwTable.SourceCatalog(self.schema) + results.negative.makeSources(negatives) + sources = results.sources self.processResults(science, matchedTemplate, difference, sources, idFactory, kernelSources, @@ -823,6 +825,8 @@ def _prepareInputs(self, difference): sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius() if np.isnan(sigma): raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.") + if self.sigma is None: + self.sigma = sigma # Ensure that we start with an empty detection and deblended mask. mask = difference.mask for mp in self.config.clearMaskPlanes: @@ -966,6 +970,10 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor # All positive peaks were added before any negative peaks, so # re-order the peaks in each footprint by significance. self._reorderPeaksBySignificance(initialDiaSources) + + # Extract isolated and purely positive peaks from large footprints + if self.config.doDeblend: + initialDiaSources = self._deblendClusteredPeaks(initialDiaSources, difference) else: initialDiaSources = sources @@ -1057,73 +1065,75 @@ def _reorderPeaksBySignificance(self, diaSources): peaks.clear() source.getFootprint().setPeakCatalog(newPeaks) - 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 _deblendClusteredPeaks(self, diaSources, difference): + """Extract isolated and purely positive peaks from large footprints. + + For each footprint that contains an isolated positive peak alongside + at least one other peak, the isolated peak's pixels are assigned to + a new footprint (a new diaSource), and the original source's footprint + is reduced to the remaining pixels. Both the original and the new + records are flagged ``is_deblended``. Footprints without an isolated + positive peak are left untouched. Parameters ---------- - difference : `lsst.afw.image.Exposure` - Result of subtracting template from the science image. - positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet` - Positive and negative polarity footprints measured on - ``difference`` to be deblended separately. + diaSources : `lsst.afw.table.SourceCatalog` + The merged diaSources; modified in place and possibly extended. + difference : `lsst.afw.image.ExposureF` + Difference image the footprints were detected on, used to weight + the peaks. Returns ------- - sources : `lsst.afw.table.SourceCatalog` - Positive and negative deblended children. - positives, negatives : `lsst.afw.table.SourceCatalog` - Deblended positive and negative polarity sources with footprints - detected on ``difference``. + diaSources : `lsst.afw.table.SourceCatalog` + The updated, contiguous catalog. """ - - def deblend(footprints, negative=False): - """Deblend a positive or negative footprint set, - and return the deblended children. - - Parameters - ---------- - footprints : `lsst.afw.detection.FootprintSet` - negative : `bool` - Set True if the footprints contain negative fluxes - - Returns - ------- - sources : `lsst.afw.table.SourceCatalog` - """ - sources = afwTable.SourceCatalog(self.schema) - footprints.makeSources(sources) - if negative: - # Invert the image so the deblender can run on positive peaks - difference_inverted = difference.clone() - difference_inverted.image *= -1 - self.deblend.run(exposure=difference_inverted, sources=sources) - children = sources[sources["parent"] != 0] - # Set the heavy footprint pixel values back to reality - for child in children: - footprint = child.getFootprint() - array = footprint.getImageArray() - array *= -1 - else: - 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, negative=True) - - sources = afwTable.SourceCatalog(self.schema) - sources.reserve(len(positives) + len(negatives)) - sources.extend(positives, deep=True) - sources.extend(negatives, deep=True) - if len(negatives) > 0: - sources[-len(negatives):]["is_negative"] = True - return sources, positives, negatives + sigmaToFwhm = np.sqrt(8.0*np.log(2.0)) + linkingRadius = self.config.peakClusterRadius*self.sigma*sigmaToFwhm + # Identify clusters of peaks within each footprint. + # Positive peaks that share a footprint with other peaks may be + # isolated enough that we can extract them as separate diaSources. + peakClusters = find_peak_clusters(diaSources, linkingRadius) + + # Precalculate some mask information up front, so the pseudo-deblending + # function can be as lightweight as possible + presentBadPlanes = [mp for mp in self.config.badMaskPlanes + if mp in difference.mask.getMaskPlaneDict()] + badBitmask = difference.mask.getPlaneBitMask(presentBadPlanes) if presentBadPlanes else 0 + imageArray = difference.image.array + maskArray = difference.mask.array + imageXY0 = (difference.getBBox().getMinX(), difference.getBBox().getMinY()) + isDeblendedKey = self.schema.find("is_deblended").key + + nFootprints = 0 + nChildren = 0 + # Loop over each individual footprint, checking whether there are + # positive peaks that could be extracted and measured independently + for footprintClusters in peakClusters: + if not footprintClusters.is_deblendable: + continue + result = deblend_clustered_peaks(footprintClusters, imageArray, maskArray, + imageXY0, self.sigma, badBitmask) + # Exit early without modifying the footprint catalog if no changes + # are needed for this footprint + if result is None: + continue + footprintClusters.source.setFootprint(result.parent) + footprintClusters.source.set(isDeblendedKey, True) + for childFootprint in result.children: + child = diaSources.addNew() + child.setFootprint(childFootprint) + child.set(isDeblendedKey, True) + nFootprints += 1 + nChildren += len(result.children) + + self.metadata["nClusterDeblendedFootprints"] = nFootprints + self.metadata["nClusterDeblendedDiaSources"] = nChildren + if nFootprints: + self.log.info("Deblended %d footprints into %d diaSources", nFootprints, nChildren) + if not diaSources.isContiguous(): + diaSources = diaSources.copy(deep=True) + return diaSources def _rejectBadMaskedDetections(self, initialDiaSources, mask): """Remove detections whose footprint peak lies on a bad-masked pixel. @@ -1704,25 +1714,13 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources edgePix = scoreExposure.mask.array & edgeBit > 0 difference.mask.array[edgePix] |= edgeBit - if self.config.doDeblend: - sources, positives, negatives = self._deblend(difference, - detectResults.positive, - detectResults.negative) - - self.processResults(science, matchedTemplate, difference, - sources, idFactory, kernelSources, - positives=positives, - negatives=negatives, - measurementResults=measurementResults) - - else: - positives = afwTable.SourceCatalog(self.schema) - detectResults.positive.makeSources(positives) - negatives = afwTable.SourceCatalog(self.schema) - detectResults.negative.makeSources(negatives) - self.processResults(science, matchedTemplate, difference, - detectResults.sources, idFactory, kernelSources, - positives=positives, - negatives=negatives, - measurementResults=measurementResults) + positives = afwTable.SourceCatalog(self.schema) + detectResults.positive.makeSources(positives) + negatives = afwTable.SourceCatalog(self.schema) + detectResults.negative.makeSources(negatives) + self.processResults(science, matchedTemplate, difference, + detectResults.sources, idFactory, kernelSources, + positives=positives, + negatives=negatives, + measurementResults=measurementResults) return measurementResults diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 84b5d5aa..7e7cbe97 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -193,7 +193,7 @@ def test_detection_xy0(self): # Set the subtraction residual metric threshold high, since the template # is not actually subtracted from the science image. detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., - doSkySources=False) + doSkySources=False, peakClusterRadius=1.0) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, sources, @@ -825,7 +825,8 @@ def _setup_sattle_tests(self): detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., - doSkySources=False, run_sattle=True) + doSkySources=False, run_sattle=True, + peakClusterRadius=1.0) return science, matchedTemplate, difference, sources, detectionTask @@ -959,7 +960,7 @@ def test_detection_xy0(self): # Set the subtraction residual metric threshold high, since the template # is not actually subtracted from the science image. detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., - doSkySources=False) + doSkySources=False, peakClusterRadius=1.0) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, score, sources, @@ -1406,57 +1407,6 @@ class TestNegativePeaks(lsst.utils.tests.TestCase): various problems found on DM-48596. """ - def testDeblendNegatives(self): - """Test that negative peaks get deblended and not destroyed: DM-48704. - This is only a test of deblending, not of merging. - """ - # Make a science image with one blend of two positive sources, to - # subtract from an empty template, resulting in a negative diffim. - bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(100, 100)) - dataset = lsst.meas.base.tests.TestDataset(bbox) - delta = 10 - with dataset.addBlend() as family: - family.addChild(instFlux=2E5, centroid=lsst.geom.Point2D(50, 72)) - family.addChild(instFlux=2.5E5, centroid=lsst.geom.Point2D(50+delta, 74)) - science, catalog = dataset.realize(noise=1.0, - schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema()) - dataset = lsst.meas.base.tests.TestDataset(bbox) - template, _ = dataset.realize(noise=1.0, - schema=lsst.meas.base.tests.TestDataset.makeMinimalSchema()) - difference = template.clone() - difference.image -= science.image - - config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() - config.doDeblend = True - task = detectAndMeasure.DetectAndMeasureTask(config=config) - # prelude steps taken from `detectAndMeasure.run` - task._prepareInputs(difference) - table = lsst.afw.table.SourceTable.make(task.schema) - results = task.detection.run(table=table, exposure=difference, doSmooth=True) - # Just run the deblend step so we can check the footprints independently. - sources, positives, negatives = task._deblend(difference, results.positive, results.negative) - - # DM-48704 fixed a problem where the peaks were in the footprints, but - # the spans were empty. - 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). - # The sources are positive in the science image, and negative in the - # diffim, so the minimum value in the deblended negative footprint is - # 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(negatives[0].getFootprint().getImageArray().min(), - -catalog[1].getFootprint().getImageArray().max(), rtol=1e-4) - self.assertFloatsAlmostEqual(negatives[1].getFootprint().getImageArray().min(), - -catalog[2].getFootprint().getImageArray().max(), rtol=1e-4) - - self.assertEqual(sources["is_negative"].sum(), 2) - def testMergeFootprints(self): """Test that merging footprints does not cause negative ones to disappear (e.g. get merged into non-connected footprints). @@ -1541,10 +1491,10 @@ def testMergeFootprints(self): unique_peaks, counts = np.unique(peaks, axis=0, return_counts=True) self.assertEqual(np.sum(counts > 1), 0) - # There are three positive sources that should not have `is_negative` - # set, independent of how deblending/merging of negative sources is - # handled. - self.assertEqual((~result.diaSources["is_negative"]).sum(), 3) + # Two positive sources should not have `is_negative` set: the isolated + # source at (25, 26) and the positive blended pair, which stays merged + # as a single source because its peaks are within ``peakClusterRadius``. + self.assertEqual((~result.diaSources["is_negative"]).sum(), 2) def testReorderPeaksBySignificance(self): """A merged footprint whose most significant peak is negative should From 63224e713dfcd8047b2c1bf9f6630ae791a048ff Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Tue, 14 Jul 2026 15:55:19 +0200 Subject: [PATCH 04/12] Remove isolated positive constraint --- python/lsst/ip/diffim/clusterDeblend.py | 68 +++++++---------------- python/lsst/ip/diffim/detectAndMeasure.py | 39 +++++++------ python/lsst/ip/diffim/peakClustering.py | 38 ++----------- tests/test_clusterDeblend.py | 49 +++++++++++----- tests/test_peakClustering.py | 63 ++++++++------------- 5 files changed, 102 insertions(+), 155 deletions(-) diff --git a/python/lsst/ip/diffim/clusterDeblend.py b/python/lsst/ip/diffim/clusterDeblend.py index 2fa73e82..830302f3 100644 --- a/python/lsst/ip/diffim/clusterDeblend.py +++ b/python/lsst/ip/diffim/clusterDeblend.py @@ -33,11 +33,6 @@ import lsst.afw.image -# Per-peak flag marking an individual peak as originating from the negative -# detection set (see peakClustering._NEGATIVE_PEAK_FLAG). -_NEGATIVE_PEAK_FLAG = "merge_peak_negative" - - @dataclasses.dataclass class DeblendedPeaks: """The footprints resulting from deblending one multi-peak footprint. @@ -45,13 +40,10 @@ class DeblendedPeaks: Parameters ---------- parent : `lsst.afw.detection.Footprint` - Footprint to set back on the original source record. This is the - remainder (all non-carved clusters) if there is one; otherwise, when - the footprint is entirely single positive peaks, it is one of those - carved peaks. + Footprint of the first cluster, set back on the original source record. children : `list` [`lsst.afw.detection.Footprint`] - One new footprint per remaining carved single positive peak, each - carrying that single peak. + One new footprint per remaining cluster, each carrying that cluster's + peaks. """ parent: lsst.afw.detection.Footprint @@ -147,7 +139,7 @@ def _footprint_from_bool(pixels, bbox, peak_schema, region): def deblend_clustered_peaks(footprint_clusters, image, mask, image_xy0, sigma, bad_bitmask): - """Extract the isolated positive peaks out of one footprint. + """Extract every cluster of peaks from a footprint into new footprints. Parameters ---------- @@ -168,22 +160,20 @@ def deblend_clustered_peaks(footprint_clusters, image, mask, image_xy0, sigma, b ------- result : `DeblendedPeaks` or `None` The parent footprint (kept on the original record) and the new - footprints, or `None` if the footprint should be left unmodified. + footprints, or `None` if the footprint did not split into more than + one cluster. """ footprint = footprint_clusters.footprint peaks = footprint.getPeaks() n_peaks = len(peaks) peak_schema = peaks.getSchema() - negative_key = peak_schema.find(_NEGATIVE_PEAK_FLAG).key x0, y0 = image_xy0 positions = np.zeros((n_peaks, 2), dtype=int) amplitudes = np.zeros(n_peaks, dtype=float) - is_negative = np.zeros(n_peaks, dtype=bool) for i, peak in enumerate(peaks): ix, iy = peak.getIx(), peak.getIy() positions[i] = (ix, iy) - is_negative[i] = peak.get(negative_key) amplitudes[i] = _peak_amplitude(image, mask, ix - x0, iy - y0, bad_bitmask) # weighted[p, j]: influence of peak p evaluated at peak j's pixel. @@ -195,16 +185,9 @@ def deblend_clustered_peaks(footprint_clusters, image, mask, image_xy0, sigma, b peak_clusters = _merge_peak_clusters( [cluster.peak_indices for cluster in footprint_clusters.clusters], weighted, n_peaks) - def is_isolated_positive(peak_cluster): - return len(peak_cluster) == 1 and not is_negative[peak_cluster[0]] - - # Extract each single positive peak into its own source. - # Every other cluster pools into one remainder. - extracted_peaks = [peak_cluster for peak_cluster in peak_clusters if is_isolated_positive(peak_cluster)] - remainder_peak_clusters = [peak_cluster for peak_cluster in peak_clusters - if not is_isolated_positive(peak_cluster)] - n_sources = len(extracted_peaks) + (1 if remainder_peak_clusters else 0) - if not extracted_peaks or n_sources < 2: + # Extract every cluster into its own source; skip footprints that did not + # split into more than one cluster. + if len(peak_clusters) < 2: return None # Assign every footprint pixel to the cluster of highest influence. @@ -213,9 +196,8 @@ def is_isolated_positive(peak_cluster): footprint_mask = _spans_to_bool(footprint.getSpans(), bbox) xx, yy = np.meshgrid(np.arange(min_x, bbox.getMaxX() + 1), np.arange(min_y, bbox.getMaxY() + 1)) - ordered = extracted_peaks + remainder_peak_clusters - influence = np.empty((len(ordered), bbox.getHeight(), bbox.getWidth()), dtype=float) - for index, peak_cluster in enumerate(ordered): + influence = np.empty((len(peak_clusters), bbox.getHeight(), bbox.getWidth()), dtype=float) + for index, peak_cluster in enumerate(peak_clusters): accumulator = np.zeros(footprint_mask.shape, dtype=float) for peak in peak_cluster: r_sq = (xx - positions[peak, 0])**2 + (yy - positions[peak, 1])**2 @@ -224,28 +206,18 @@ def is_isolated_positive(peak_cluster): labels = np.argmax(influence, axis=0) # Guarantee each cluster owns its own peak pixels (breaks any influence # ties in the owner's favour, so no cluster ends up empty). - for index, peak_cluster in enumerate(ordered): + for index, peak_cluster in enumerate(peak_clusters): for peak in peak_cluster: labels[positions[peak, 1] - min_y, positions[peak, 0] - min_x] = index + # One footprint per cluster, carrying that cluster's peaks. region = footprint.getRegion() - # One footprint per extracted single positive peak. - extracted_footprints = [] - for index, peak_cluster in enumerate(extracted_peaks): + footprints = [] + for index, peak_cluster in enumerate(peak_clusters): fp = _footprint_from_bool(footprint_mask & (labels == index), bbox, peak_schema, region) - fp.getPeaks().append(peaks[peak_cluster[0]]) - extracted_footprints.append(fp) - - if remainder_peak_clusters: - remainder_pixels = np.zeros(footprint_mask.shape, dtype=bool) - for index in range(len(extracted_peaks), len(ordered)): - remainder_pixels |= footprint_mask & (labels == index) - remainder = _footprint_from_bool(remainder_pixels, bbox, peak_schema, region) - for peak_cluster in remainder_peak_clusters: - for peak in peak_cluster: - remainder.getPeaks().append(peaks[peak]) - return DeblendedPeaks(parent=remainder, children=extracted_footprints) + for peak in peak_cluster: + fp.getPeaks().append(peaks[peak]) + footprints.append(fp) - # No remainder cluster: the footprint is entirely single positive peaks, - # so the first stays on the original record and the rest are new records. - return DeblendedPeaks(parent=extracted_footprints[0], children=extracted_footprints[1:]) + # The first cluster stays on the original record; the rest are new records. + return DeblendedPeaks(parent=footprints[0], children=footprints[1:]) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 8a526a55..5b9297ef 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -1066,14 +1066,14 @@ def _reorderPeaksBySignificance(self, diaSources): source.getFootprint().setPeakCatalog(newPeaks) def _deblendClusteredPeaks(self, diaSources, difference): - """Extract isolated and purely positive peaks from large footprints. + """Extract clusters of peaks from large footprints into new diaSources. - For each footprint that contains an isolated positive peak alongside - at least one other peak, the isolated peak's pixels are assigned to - a new footprint (a new diaSource), and the original source's footprint - is reduced to the remaining pixels. Both the original and the new - records are flagged ``is_deblended``. Footprints without an isolated - positive peak are left untouched. + For each footprint whose peaks split into more than one cluster, every + cluster is assigned its own footprint: the first stays on the original + record and the rest become new diaSources. All resulting records are + flagged ``is_deblended``, and ``is_negative`` is set from the polarity + of each new footprint's peaks. Footprints with a single cluster are + left untouched. Parameters ---------- @@ -1090,9 +1090,8 @@ def _deblendClusteredPeaks(self, diaSources, difference): """ sigmaToFwhm = np.sqrt(8.0*np.log(2.0)) linkingRadius = self.config.peakClusterRadius*self.sigma*sigmaToFwhm - # Identify clusters of peaks within each footprint. - # Positive peaks that share a footprint with other peaks may be - # isolated enough that we can extract them as separate diaSources. + # Identify clusters of peaks within each footprint. Footprints whose + # peaks split into more than one cluster can be deblended. peakClusters = find_peak_clusters(diaSources, linkingRadius) # Precalculate some mask information up front, so the pseudo-deblending @@ -1104,11 +1103,12 @@ def _deblendClusteredPeaks(self, diaSources, difference): maskArray = difference.mask.array imageXY0 = (difference.getBBox().getMinX(), difference.getBBox().getMinY()) isDeblendedKey = self.schema.find("is_deblended").key + isNegativeKey = self.schema.find("is_negative").key nFootprints = 0 nChildren = 0 - # Loop over each individual footprint, checking whether there are - # positive peaks that could be extracted and measured independently + # Loop over each individual footprint, checking whether its peaks split + # into clusters that can be extracted and measured independently for footprintClusters in peakClusters: if not footprintClusters.is_deblendable: continue @@ -1118,12 +1118,15 @@ def _deblendClusteredPeaks(self, diaSources, difference): # are needed for this footprint if result is None: continue - footprintClusters.source.setFootprint(result.parent) - footprintClusters.source.set(isDeblendedKey, True) - for childFootprint in result.children: - child = diaSources.addNew() - child.setFootprint(childFootprint) - child.set(isDeblendedKey, True) + # The first cluster stays on the original record; the rest are new. + records = [footprintClusters.source] + records += [diaSources.addNew() for _ in result.children] + negativePeakKey = result.parent.getPeaks().getSchema().find("merge_peak_negative").key + for record, newFootprint in zip(records, [result.parent] + result.children): + record.setFootprint(newFootprint) + record.set(isDeblendedKey, True) + # A source is negative if all of its peaks are negative. + record.set(isNegativeKey, all(peak.get(negativePeakKey) for peak in newFootprint.getPeaks())) nFootprints += 1 nChildren += len(result.children) diff --git a/python/lsst/ip/diffim/peakClustering.py b/python/lsst/ip/diffim/peakClustering.py index 583a6fca..27ef01e2 100644 --- a/python/lsst/ip/diffim/peakClustering.py +++ b/python/lsst/ip/diffim/peakClustering.py @@ -47,8 +47,7 @@ class PeakCluster: Peaks are grouped by single-linkage clustering on their positions, ignoring polarity, so a dipole's positive and negative lobes should fall - into the same cluster. A cluster with a single member is "isolated": no other - peak in the footprint lies within the linking radius. + into the same cluster. Parameters ---------- @@ -70,19 +69,6 @@ class PeakCluster: peak_values: np.ndarray is_negative: np.ndarray - @property - def is_singleton(self): - """Whether this cluster holds exactly one peak (`bool`).""" - return len(self.peak_indices) == 1 - - @property - def is_isolated_positive(self): - """Whether this cluster is a single positive peak (`bool`). - - These are the peaks the deblender extracts as new footprints. - """ - return self.is_singleton and not bool(self.is_negative[0]) - @dataclasses.dataclass class FootprintPeakClusters: @@ -122,28 +108,12 @@ def n_peaks(self): """Total number of peaks in the footprint (`int`).""" return sum(len(cluster.peak_indices) for cluster in self.clusters) - @property - def isolated_positive(self): - """Single-peak positive clusters (`list` [`PeakCluster`]). - - The candidates for extracting into their own footprints. - """ - return [cluster for cluster in self.clusters if cluster.is_isolated_positive] - - @property - def retained(self): - """Clusters whose pixels stay in the footprint (`list` [`PeakCluster`]). - - Everything that is not an isolated positive peak: multi-peak blends, - dipoles, and isolated negative peaks. - """ - return [cluster for cluster in self.clusters if not cluster.is_isolated_positive] - @property def is_deblendable(self): - """Whether the footprint can be split by cluster deblending. + """Whether the footprint splits into more than one cluster, so that + deblending would extract more than one source (`bool`). """ - return bool(self.isolated_positive) and len(self.clusters) >= 2 + return len(self.clusters) >= 2 def _single_linkage_labels(positions, linking_radius): diff --git a/tests/test_clusterDeblend.py b/tests/test_clusterDeblend.py index 1be97221..61c86760 100644 --- a/tests/test_clusterDeblend.py +++ b/tests/test_clusterDeblend.py @@ -86,9 +86,9 @@ class ClusterDeblendTestCase(lsst.utils.tests.TestCase): def setUp(self): self.schema, self.peakSchema = makeSchemas() - def test_extract_isolated_positive_from_dipole(self): - """An isolated positive peak is split off; the dipole stays behind.""" - # Isolated positive at (8, 20); dipole at (30, 20)/(32, 20). + def test_extract_clusters_from_footprint(self): + """A footprint that splits into two clusters yields two sources.""" + # Lone peak at (8, 20); dipole at (30, 20)/(32, 20). footprint, image, mask, _ = makeFootprintAndImage( self.peakSchema, [(8, 20, False), (30, 20, False), (32, 20, True)]) fc = makeClusters(self.schema, footprint, linkingRadius=5.0) @@ -98,18 +98,19 @@ def test_extract_isolated_positive_from_dipole(self): self.assertIsNotNone(result) self.assertEqual(len(result.children), 1) - child = result.children[0] - self.assertEqual(len(child.getPeaks()), 1) - self.assertEqual((child.getPeaks()[0].getIx(), child.getPeaks()[0].getIy()), (8, 20)) - self.assertTrue(child.getSpans().contains(geom.Point2I(8, 20))) + # The first cluster (the lone peak at (8, 20)) stays on the parent. + self.assertEqual(len(result.parent.getPeaks()), 1) + self.assertEqual((result.parent.getPeaks()[0].getIx(), + result.parent.getPeaks()[0].getIy()), (8, 20)) - # The dipole is the remainder (parent), keeping both its peaks. - self.assertEqual(len(result.parent.getPeaks()), 2) - self.assertTrue(result.parent.getSpans().contains(geom.Point2I(30, 20))) - self.assertTrue(result.parent.getSpans().contains(geom.Point2I(32, 20))) + # The dipole is the other cluster, keeping both its peaks. + child = result.children[0] + self.assertEqual(len(child.getPeaks()), 2) + self.assertTrue(child.getSpans().contains(geom.Point2I(30, 20))) + self.assertTrue(child.getSpans().contains(geom.Point2I(32, 20))) - # The extracted and remainder pixels exactly tile the original footprint. - self.assertEqual(child.getSpans().getArea() + result.parent.getSpans().getArea(), + # The extracted pixels exactly tile the original footprint. + self.assertEqual(result.parent.getSpans().getArea() + child.getSpans().getArea(), footprint.getSpans().getArea()) # Peak schema (with merge flags) is preserved on the outputs. self.assertIn("merge_peak_negative", child.getPeaks().schema.getNames()) @@ -144,8 +145,8 @@ def test_within_radius_positive_blend_stays(self): result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) self.assertIsNone(result) - def test_no_isolated_positive_returns_none(self): - """A footprint with only a dipole has nothing to extract.""" + def test_single_cluster_dipole_returns_none(self): + """A footprint that is a single cluster (here a dipole) is untouched.""" footprint, image, mask, _ = makeFootprintAndImage( self.peakSchema, [(20, 20, False), (22, 20, True)]) fc = makeClusters(self.schema, footprint, linkingRadius=5.0) @@ -153,6 +154,24 @@ def test_no_isolated_positive_returns_none(self): result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) self.assertIsNone(result) + def test_extract_negative_cluster(self): + """Negative and positive clusters are both extracted (polarity is not + a constraint in this version). + """ + # Lone negative at (8, 20); lone positive at (32, 20); beyond radius. + footprint, image, mask, _ = makeFootprintAndImage( + self.peakSchema, [(8, 20, True), (32, 20, False)]) + fc = makeClusters(self.schema, footprint, linkingRadius=5.0) + self.assertTrue(fc.is_deblendable) + + result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) + self.assertIsNotNone(result) + self.assertEqual(len(result.children), 1) + # Both peaks are extracted, one per output footprint. + outputs = [result.parent] + result.children + extracted = {(fp.getPeaks()[0].getIx(), fp.getPeaks()[0].getIy()) for fp in outputs} + self.assertEqual(extracted, {(8, 20), (32, 20)}) + def test_dominated_isolated_peak_is_merged(self): """A faint 'isolated' peak dominated by a bright neighbour is merged back in, leaving nothing to extract. diff --git a/tests/test_peakClustering.py b/tests/test_peakClustering.py index c88261f1..18e87688 100644 --- a/tests/test_peakClustering.py +++ b/tests/test_peakClustering.py @@ -82,8 +82,7 @@ def test_two_peaks_far_apart(self): clusters = cluster_peaks(footprint.getPeaks(), self.radius) self.assertEqual(len(clusters), 2) for cluster in clusters: - self.assertTrue(cluster.is_singleton) - self.assertTrue(cluster.is_isolated_positive) + self.assertEqual(len(cluster.peak_indices), 1) def test_two_peaks_close(self): """Peaks within the linking radius merge into one cluster.""" @@ -91,8 +90,6 @@ def test_two_peaks_close(self): [(10, 10, 100, False), (12, 10, 90, False)]) clusters = cluster_peaks(footprint.getPeaks(), self.radius) self.assertEqual(len(clusters), 1) - self.assertFalse(clusters[0].is_singleton) - self.assertFalse(clusters[0].is_isolated_positive) self.assertEqual(sorted(clusters[0].peak_indices), [0, 1]) def test_transitive_chain(self): @@ -106,37 +103,35 @@ def test_transitive_chain(self): self.assertEqual(sorted(clusters[0].peak_indices), [0, 1, 2]) def test_dipole_single_cluster(self): - """A close positive/negative pair clusters together, not isolated.""" + """A close positive/negative pair clusters together.""" footprint = makeFootprint(self.peakSchema, [(10, 10, 100, False), (12, 10, -90, True)]) clusters = cluster_peaks(footprint.getPeaks(), self.radius) self.assertEqual(len(clusters), 1) cluster = clusters[0] - self.assertFalse(cluster.is_isolated_positive) self.assertEqual(cluster.is_negative.tolist(), [False, True]) self.assertEqual(cluster.peak_values.tolist(), [100, -90]) - def test_isolated_positive_beside_dipole(self): - """An isolated positive peak is separated from a nearby dipole.""" + def test_lone_peak_beside_dipole(self): + """A lone peak is separated from a nearby dipole into two clusters.""" footprint = makeFootprint(self.peakSchema, - [(50, 50, 200, False), # isolated positive + [(50, 50, 200, False), # lone peak (10, 10, 100, False), (12, 10, -90, True)]) # dipole clusters = cluster_peaks(footprint.getPeaks(), self.radius) self.assertEqual(len(clusters), 2) - isolated = clusterContaining(clusters, 0) + lone = clusterContaining(clusters, 0) dipole = clusterContaining(clusters, 1) - self.assertTrue(isolated.is_isolated_positive) - self.assertFalse(dipole.is_isolated_positive) + self.assertEqual(lone.peak_indices, [0]) self.assertEqual(sorted(dipole.peak_indices), [1, 2]) - def test_isolated_negative_not_extracted(self): - """A lone negative peak is a singleton but not an isolated positive.""" + def test_lone_negative_is_own_cluster(self): + """A lone negative peak forms its own singleton cluster.""" footprint = makeFootprint(self.peakSchema, [(50, 50, 200, False), (10, 10, -80, True)]) clusters = cluster_peaks(footprint.getPeaks(), self.radius) negative = clusterContaining(clusters, 1) - self.assertTrue(negative.is_singleton) - self.assertFalse(negative.is_isolated_positive) + self.assertEqual(negative.peak_indices, [1]) + self.assertEqual(negative.is_negative.tolist(), [True]) def test_positions_and_values(self): """Cluster arrays carry the peaks' positions and signed values.""" @@ -187,35 +182,23 @@ def test_handles_missing_footprint(self): self.assertEqual(results[0].source_id, 2) def test_deblendable_classification(self): - """is_deblendable requires both an extraction candidate and a - remainder. - """ - # Isolated positive + dipole: deblendable. - mixed = makeFootprint(self.peakSchema, - [(50, 50, 200, False), - (10, 10, 100, False), (12, 10, -90, True)]) - # Two isolated positives beyond the radius: deblendable into two - # sources even though there is no remainder cluster. - allIsolated = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (50, 50, 120, False)]) - # One tight blend: a remainder but nothing isolated to extract. - blend = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (12, 10, 90, False)]) - catalog = self.makeCatalog([mixed, allIsolated, blend]) + """is_deblendable is true when the footprint splits into >1 cluster.""" + # A lone peak plus a dipole: two clusters, deblendable. + twoClusters = makeFootprint(self.peakSchema, + [(50, 50, 200, False), + (10, 10, 100, False), (12, 10, -90, True)]) + # One tight blend: a single cluster, not deblendable. + oneCluster = makeFootprint(self.peakSchema, + [(10, 10, 100, False), (12, 10, 90, False)]) + catalog = self.makeCatalog([twoClusters, oneCluster]) results = find_peak_clusters(catalog, self.radius) byId = {r.source_id: r for r in results} self.assertTrue(byId[1].is_deblendable) - self.assertEqual(len(byId[1].isolated_positive), 1) - self.assertEqual(len(byId[1].retained), 1) - - self.assertTrue(byId[2].is_deblendable) - self.assertEqual(len(byId[2].isolated_positive), 2) - self.assertEqual(len(byId[2].retained), 0) + self.assertEqual(len(byId[1].clusters), 2) - self.assertFalse(byId[3].is_deblendable) - self.assertEqual(len(byId[3].isolated_positive), 0) - self.assertEqual(len(byId[3].retained), 1) + self.assertFalse(byId[2].is_deblendable) + self.assertEqual(len(byId[2].clusters), 1) def test_source_handle_is_live(self): """The stored source is the live catalog record, so deblender edits From 91ffef6a99eae4c7eb27558eadbc374e9b940ab9 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 17 Jul 2026 19:40:44 +0200 Subject: [PATCH 05/12] Add maximum footprint size constraint for deblending --- python/lsst/ip/diffim/detectAndMeasure.py | 27 ++++++++++++++ tests/test_detectAndMeasure.py | 44 +++++++++++++++++++++++ 2 files changed, 71 insertions(+) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index 5b9297ef..c601d308 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -280,6 +280,21 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, "Positive peaks that are isolated after clustering can be deblended " "if doDeblend=True." ) + maxDeblendFootprintArea = pexConfig.Field( + dtype=int, + default=250000, + doc="Footprints whose bounding-box area (pixels) exceeds this are left " + "un-deblended, bounding the per-footprint deblending cost so a " + "single large footprint cannot blow the processing time budget. " + "Only used if doDeblend=True." + ) + maxDeblendPeaks = pexConfig.Field( + dtype=int, + default=200, + doc="Footprints with more than this many peaks are left un-deblended, " + "bounding the per-footprint deblending cost. Only used if " + "doDeblend=True." + ) measurement = pexConfig.ConfigurableField( target=DipoleFitTask, doc="Task to measure sources on the difference image.", @@ -1107,11 +1122,19 @@ def _deblendClusteredPeaks(self, diaSources, difference): nFootprints = 0 nChildren = 0 + nSkipped = 0 # Loop over each individual footprint, checking whether its peaks split # into clusters that can be extracted and measured independently for footprintClusters in peakClusters: if not footprintClusters.is_deblendable: continue + # Leave pathologically large footprints un-deblended: the deblending + # cost grows steeply with footprint area and peak count, so a single + # huge footprint could otherwise dominate the processing time. + if (footprintClusters.footprint.getBBox().getArea() > self.config.maxDeblendFootprintArea + or footprintClusters.n_peaks > self.config.maxDeblendPeaks): + nSkipped += 1 + continue result = deblend_clustered_peaks(footprintClusters, imageArray, maskArray, imageXY0, self.sigma, badBitmask) # Exit early without modifying the footprint catalog if no changes @@ -1132,8 +1155,12 @@ def _deblendClusteredPeaks(self, diaSources, difference): self.metadata["nClusterDeblendedFootprints"] = nFootprints self.metadata["nClusterDeblendedDiaSources"] = nChildren + self.metadata["nClusterDeblendSkipped"] = nSkipped if nFootprints: self.log.info("Deblended %d footprints into %d diaSources", nFootprints, nChildren) + if nSkipped: + self.log.info("Left %d oversized footprints un-deblended (area > %d or peaks > %d)", + nSkipped, self.config.maxDeblendFootprintArea, self.config.maxDeblendPeaks) if not diaSources.isContiguous(): diaSources = diaSources.copy(deep=True) return diaSources diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 7e7cbe97..f48f67b7 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -1531,6 +1531,50 @@ def testReorderPeaksBySignificance(self): self.assertEqual(peaks[1].getPeakValue(), 30.0) self.assertEqual(peaks[1]["significance"], 5.0) + def _makeDeblendableCatalog(self, schema, bbox, peakXY): + """A single-source catalog whose footprint has well-separated peaks.""" + peakSchema = afwDetection.FootprintMergeList( + afwTable.SourceTable.makeMinimalSchema(), ["positive", "negative"]).getPeakSchema() + negativeKey = peakSchema.find("merge_peak_negative").key + positiveKey = peakSchema.find("merge_peak_positive").key + catalog = afwTable.SourceCatalog(schema) + source = catalog.addNew() + footprint = afwDetection.Footprint(afwGeom.SpanSet(bbox), peakSchema) + for ix, iy in peakXY: + peak = footprint.addPeak(ix, iy, 100.0) + peak.set(negativeKey, False) + peak.set(positiveKey, True) + source.setFootprint(footprint) + return catalog + + def testDeblendGuardrailSkipsManyPeaks(self): + """A deblendable footprint with more peaks than the cap is left whole.""" + config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() + config.maxDeblendPeaks = 2 + task = detectAndMeasure.DetectAndMeasureTask(config=config) + task.sigma = 1.0 + bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(40, 40)) + # Three peaks, spaced beyond the linking radius, so the footprint is + # deblendable but exceeds maxDeblendPeaks. + diaSources = self._makeDeblendableCatalog(task.schema, bbox, [(8, 20), (20, 20), (32, 20)]) + result = task._deblendClusteredPeaks(diaSources, afwImage.ExposureF(bbox)) + self.assertEqual(len(result), 1) + self.assertFalse(result[0]["is_deblended"]) + self.assertEqual(task.metadata["nClusterDeblendSkipped"], 1) + + def testDeblendGuardrailSkipsLargeArea(self): + """A deblendable footprint bigger than the area cap is left whole.""" + config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() + config.maxDeblendFootprintArea = 100 + task = detectAndMeasure.DetectAndMeasureTask(config=config) + task.sigma = 1.0 + bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(40, 40)) # area 1600 > 100 + diaSources = self._makeDeblendableCatalog(task.schema, bbox, [(8, 20), (20, 20), (32, 20)]) + result = task._deblendClusteredPeaks(diaSources, afwImage.ExposureF(bbox)) + self.assertEqual(len(result), 1) + self.assertFalse(result[0]["is_deblended"]) + self.assertEqual(task.metadata["nClusterDeblendSkipped"], 1) + def makeVisitInfo(id=1): """Return a non-NaN visitInfo.""" From aecdbd7c77f56837bbbe11fe3883a73be9c7a0c9 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 17 Jul 2026 21:42:56 +0200 Subject: [PATCH 06/12] Refactor deblending as a subtask --- python/lsst/ip/diffim/__init__.py | 3 +- python/lsst/ip/diffim/clusterDeblend.py | 223 ----------- .../ip/diffim/deblendDifferenceFootprints.py | 346 ++++++++++++++++++ python/lsst/ip/diffim/detectAndMeasure.py | 123 +------ python/lsst/ip/diffim/peakClustering.py | 230 ------------ tests/test_clusterDeblend.py | 214 ----------- tests/test_deblendDifferenceFootprints.py | 308 ++++++++++++++++ tests/test_detectAndMeasure.py | 84 ++++- tests/test_peakClustering.py | 227 ------------ 9 files changed, 736 insertions(+), 1022 deletions(-) delete mode 100644 python/lsst/ip/diffim/clusterDeblend.py create mode 100644 python/lsst/ip/diffim/deblendDifferenceFootprints.py delete mode 100644 python/lsst/ip/diffim/peakClustering.py delete mode 100644 tests/test_clusterDeblend.py create mode 100644 tests/test_deblendDifferenceFootprints.py delete mode 100644 tests/test_peakClustering.py diff --git a/python/lsst/ip/diffim/__init__.py b/python/lsst/ip/diffim/__init__.py index cee6f8eb..85f3cb6a 100755 --- a/python/lsst/ip/diffim/__init__.py +++ b/python/lsst/ip/diffim/__init__.py @@ -34,8 +34,7 @@ from .dipoleMeasurement import * from .getTemplate import * from .dipoleFitTask import * -from .peakClustering import * -from .clusterDeblend import * +from .deblendDifferenceFootprints import * from .detectAndMeasure import * from .imageDecorrelation import * from .imageMapReduce import * diff --git a/python/lsst/ip/diffim/clusterDeblend.py b/python/lsst/ip/diffim/clusterDeblend.py deleted file mode 100644 index 830302f3..00000000 --- a/python/lsst/ip/diffim/clusterDeblend.py +++ /dev/null @@ -1,223 +0,0 @@ -# This file is part of ip_diffim. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -"""Simple footprint deblender for difference images. -""" - -__all__ = ["DeblendedPeaks", "deblend_clustered_peaks"] - -import dataclasses - -import numpy as np - -import lsst.afw.detection -import lsst.afw.geom -import lsst.afw.image - - -@dataclasses.dataclass -class DeblendedPeaks: - """The footprints resulting from deblending one multi-peak footprint. - - Parameters - ---------- - parent : `lsst.afw.detection.Footprint` - Footprint of the first cluster, set back on the original source record. - children : `list` [`lsst.afw.detection.Footprint`] - One new footprint per remaining cluster, each carrying that cluster's - peaks. - """ - - parent: lsst.afw.detection.Footprint - children: list - - -def _peak_amplitude(image, mask, col, row, bad_bitmask): - """Summed absolute image flux in the 3x3 box around a peak. - - Parameters - ---------- - image, mask : `numpy.ndarray` - Difference-image and mask arrays. - col, row : `int` - Peak location in array coordinates (i.e. image coordinates minus the - image origin). - bad_bitmask : `int` - Bitmask of mask planes to exclude from the sum. - - Returns - ------- - amplitude : `float` - """ - height, width = image.shape - c0, c1 = max(col - 1, 0), min(col + 2, width) - r0, r1 = max(row - 1, 0), min(row + 2, height) - sub_image = image[r0:r1, c0:c1] - good = (mask[r0:r1, c0:c1] & bad_bitmask) == 0 - return float(np.sum(np.abs(sub_image[good]))) - - -def _merge_peak_clusters(peak_clusters, weighted, n_peaks): - """Merge clusters of peaks until every peak owns (dominates) its own peak - pixel. - - Parameters - ---------- - peaks : `list` [`list` [`int`]] - Initial peaks, each a list of peak indices. - weighted : `numpy.ndarray` - ``(n_peaks, n_peaks)`` array where ``weighted[p, j]`` is peak ``p``'s - influence evaluated at peak ``j``'s pixel. - n_peaks : `int` - Number of peaks. - - Returns - ------- - peaks : `list` [`list` [`int`]] - The merged peaks. - """ - peak_clusters = [list(peak_cluster) for peak_cluster in peak_clusters] - while True: - owner = [0]*n_peaks - for index, peak_cluster in enumerate(peak_clusters): - for peak in peak_cluster: - owner[peak] = index - merged = False - for j in range(n_peaks): - # Check if the current peak is the most significant influence at - # its location. If not, add it to the more significant cluster. - influence = np.array([weighted[peak_cluster, j].sum() for peak_cluster in peak_clusters]) - best = int(np.argmax(influence)) - # Merge only if another cluster of peaks *strictly* dominates the - # owner at this peak pixel. - if influence[best] > influence[owner[j]]: - low, high = sorted((owner[j], best)) - peak_clusters[low] = peak_clusters[low] + peak_clusters[high] - del peak_clusters[high] - merged = True - break - if not merged: - return peak_clusters - - -def _spans_to_bool(span_set, bbox): - """Rasterize a SpanSet to a boolean array over ``bbox``.""" - array = np.zeros((bbox.getHeight(), bbox.getWidth()), dtype=bool) - min_x, min_y = bbox.getMinX(), bbox.getMinY() - for span in span_set: - array[span.getY() - min_y, span.getX0() - min_x:span.getX1() - min_x + 1] = True - return array - - -def _footprint_from_bool(pixels, bbox, peak_schema, region): - """Build a peak-schema-preserving footprint from a boolean pixel array.""" - mask = lsst.afw.image.Mask(bbox) - bit = mask.getPlaneBitMask("DETECTED") - mask.array[:] = 0 - mask.array[pixels] = bit - footprint = lsst.afw.detection.Footprint(lsst.afw.geom.SpanSet.fromMask(mask, bit), peak_schema) - footprint.setRegion(region) - return footprint - - -def deblend_clustered_peaks(footprint_clusters, image, mask, image_xy0, sigma, bad_bitmask): - """Extract every cluster of peaks from a footprint into new footprints. - - Parameters - ---------- - footprint_clusters : `~lsst.ip.diffim.peakClustering.FootprintPeakClusters` - Clustering result for one multi-peak footprint. - image, mask : `numpy.ndarray` - Difference-image and mask arrays used to weight the peaks. - image_xy0 : `tuple` [`int`, `int`] - Image origin ``(x0, y0)``, so peak image coordinates can index the - arrays. - sigma : `float` - Width (in pixels) of the 2D Gaussian to use to determine the relative - significance of each peak to the pixels in the footprint. - bad_bitmask : `int` - Bitmask of mask planes to exclude from the peak amplitude sums. - - Returns - ------- - result : `DeblendedPeaks` or `None` - The parent footprint (kept on the original record) and the new - footprints, or `None` if the footprint did not split into more than - one cluster. - """ - footprint = footprint_clusters.footprint - peaks = footprint.getPeaks() - n_peaks = len(peaks) - peak_schema = peaks.getSchema() - x0, y0 = image_xy0 - - positions = np.zeros((n_peaks, 2), dtype=int) - amplitudes = np.zeros(n_peaks, dtype=float) - for i, peak in enumerate(peaks): - ix, iy = peak.getIx(), peak.getIy() - positions[i] = (ix, iy) - amplitudes[i] = _peak_amplitude(image, mask, ix - x0, iy - y0, bad_bitmask) - - # weighted[p, j]: influence of peak p evaluated at peak j's pixel. - inv_two_sigma_sq = 1.0/(2.0*sigma*sigma) - delta = positions[:, np.newaxis, :] - positions[np.newaxis, :, :] - dist_sq = np.sum(delta*delta, axis=2) - weighted = amplitudes[:, np.newaxis]*np.exp(-dist_sq*inv_two_sigma_sq) - - peak_clusters = _merge_peak_clusters( - [cluster.peak_indices for cluster in footprint_clusters.clusters], weighted, n_peaks) - - # Extract every cluster into its own source; skip footprints that did not - # split into more than one cluster. - if len(peak_clusters) < 2: - return None - - # Assign every footprint pixel to the cluster of highest influence. - bbox = footprint.getBBox() - min_x, min_y = bbox.getMinX(), bbox.getMinY() - footprint_mask = _spans_to_bool(footprint.getSpans(), bbox) - xx, yy = np.meshgrid(np.arange(min_x, bbox.getMaxX() + 1), - np.arange(min_y, bbox.getMaxY() + 1)) - influence = np.empty((len(peak_clusters), bbox.getHeight(), bbox.getWidth()), dtype=float) - for index, peak_cluster in enumerate(peak_clusters): - accumulator = np.zeros(footprint_mask.shape, dtype=float) - for peak in peak_cluster: - r_sq = (xx - positions[peak, 0])**2 + (yy - positions[peak, 1])**2 - accumulator += amplitudes[peak]*np.exp(-r_sq*inv_two_sigma_sq) - influence[index] = accumulator - labels = np.argmax(influence, axis=0) - # Guarantee each cluster owns its own peak pixels (breaks any influence - # ties in the owner's favour, so no cluster ends up empty). - for index, peak_cluster in enumerate(peak_clusters): - for peak in peak_cluster: - labels[positions[peak, 1] - min_y, positions[peak, 0] - min_x] = index - - # One footprint per cluster, carrying that cluster's peaks. - region = footprint.getRegion() - footprints = [] - for index, peak_cluster in enumerate(peak_clusters): - fp = _footprint_from_bool(footprint_mask & (labels == index), bbox, peak_schema, region) - for peak in peak_cluster: - fp.getPeaks().append(peaks[peak]) - footprints.append(fp) - - # The first cluster stays on the original record; the rest are new records. - return DeblendedPeaks(parent=footprints[0], children=footprints[1:]) diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py new file mode 100644 index 00000000..6e0178a8 --- /dev/null +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -0,0 +1,346 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Deblend merged difference-image footprints by clustering their peaks. +""" + +__all__ = ["DeblendDifferenceFootprintsConfig", "DeblendDifferenceFootprintsTask"] + +import dataclasses + +import numpy as np +from scipy.sparse import coo_matrix +from scipy.sparse.csgraph import connected_components +from scipy.spatial import cKDTree + +import lsst.afw.detection +import lsst.afw.geom +import lsst.afw.image +import lsst.pex.config as pexConfig +import lsst.pipe.base as pipeBase + +# Per-peak flag set by `FootprintMergeList` marking an individual peak as +# originating from the negative-polarity detection footprint set. +_NEGATIVE_PEAK_FLAG = "merge_peak_negative" + +# Conversion from a Gaussian sigma to its full width at half maximum. +_SIGMA_TO_FWHM = np.sqrt(8.0*np.log(2.0)) + + +@dataclasses.dataclass +class _DeblendedFootprint: + """The footprints resulting from deblending one multi-peak footprint. + + Parameters + ---------- + parent : `lsst.afw.detection.Footprint` + Footprint of the first cluster, set back on the original source record. + children : `list` [`lsst.afw.detection.Footprint`] + One new footprint per remaining cluster, each carrying that cluster's + peaks. + """ + + parent: lsst.afw.detection.Footprint + children: list + + +def _cluster_peak_indices(peaks, linking_radius): + """Group a footprint's peaks into spatial clusters. + + Peaks are clustered with a KD tree with spacing equal to the linking radius. + + Parameters + ---------- + peaks : `lsst.afw.detection.PeakCatalog` + Peaks of a single (merged) footprint. + linking_radius : `float` + Maximum distance, in pixels, between two peaks for them to be + clustered together. + + Returns + ------- + clusters : `list` [`list` [`int`]] + One entry per cluster, each a list of rows into ``peaks``. + """ + positions = np.array([(peak.getIx(), peak.getIy()) for peak in peaks], dtype=float) + n = len(positions) + if n <= 1: + return [list(range(n))] + pairs = cKDTree(positions).query_pairs(r=linking_radius, output_type="ndarray") + if len(pairs) == 0: + graph = coo_matrix((n, n)) + else: + graph = coo_matrix((np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n)) + _, labels = connected_components(graph, directed=False) + return [np.nonzero(labels == label)[0].tolist() for label in np.unique(labels)] + + +def _merge_dominated_clusters(clusters, weighted, n_peaks): + """Merge clusters until every peak owns (dominates) its own peak pixel. + + A cluster is absorbed into another when the other has a strictly greater + summed influence at one of the first cluster's own peak pixels. + + Parameters + ---------- + clusters : `list` [`list` [`int`]] + Initial clusters, each a list of peak indices. + weighted : `numpy.ndarray` + ``(n_peaks, n_peaks)`` array where ``weighted[p, j]`` is peak ``p``'s + influence evaluated at peak ``j``'s pixel. + n_peaks : `int` + Number of peaks. + + Returns + ------- + clusters : `list` [`list` [`int`]] + The merged clusters. + """ + clusters = [list(cluster) for cluster in clusters] + while True: + owner = [0]*n_peaks + for index, cluster in enumerate(clusters): + for peak in cluster: + owner[peak] = index + merged = False + for j in range(n_peaks): + # Check whether the owning cluster has the most significant + # influence at peak j's pixel; if not, merge it into the one that + # does. + influence = np.array([weighted[cluster, j].sum() for cluster in clusters]) + best = int(np.argmax(influence)) + if influence[best] > influence[owner[j]]: + low, high = sorted((owner[j], best)) + clusters[low] = clusters[low] + clusters[high] + del clusters[high] + merged = True + break + if not merged: + return clusters + + +def _spans_to_bool(span_set, bbox): + """Rasterize a SpanSet to a boolean array over ``bbox``.""" + array = np.zeros((bbox.getHeight(), bbox.getWidth()), dtype=bool) + min_x, min_y = bbox.getMinX(), bbox.getMinY() + for span in span_set: + array[span.getY() - min_y, span.getX0() - min_x:span.getX1() - min_x + 1] = True + return array + + +def _footprint_from_bool(pixels, bbox, peak_schema, region): + """Build a peak-schema-preserving footprint from a boolean pixel array.""" + mask = lsst.afw.image.Mask(bbox) + bit = mask.getPlaneBitMask("DETECTED") + mask.array[:] = 0 + mask.array[pixels] = bit + footprint = lsst.afw.detection.Footprint(lsst.afw.geom.SpanSet.fromMask(mask, bit), peak_schema) + footprint.setRegion(region) + return footprint + + +def _deblend_footprint(footprint, clusters, sigma): + """Extract every cluster of peaks from a footprint into new footprints. + + Each peak's influence is calculated as a gaussian-approximation of the PSF + with amplitude equal to the absolute value of the signal to noise of the + peak on the detection image (```peak['significance']```). + + Parameters + ---------- + footprint : `lsst.afw.detection.Footprint` + The multi-peak footprint to deblend. + clusters : `list` [`list` [`int`]] + Clusters of peaks within the footprint, as rows into + ``footprint.getPeaks()``. + sigma : `float` + Width (in pixels) of the 2D Gaussian used to weight each peak's + influence on the footprint pixels. + + Returns + ------- + result : `_DeblendedFootprint` or `None` + The parent footprint (kept on the original record) and the new + footprints, or `None` if the footprint did not split into more than + one cluster. + """ + peaks = footprint.getPeaks() + n_peaks = len(peaks) + peak_schema = peaks.getSchema() + + positions = np.array([(peak.getIx(), peak.getIy()) for peak in peaks], dtype=int) + amplitudes = np.array([peak["significance"] for peak in peaks], dtype=float) + + # weighted[p, j]: influence of peak p evaluated at peak j's pixel. + inv_two_sigma_sq = 1.0/(2.0*sigma*sigma) + delta = positions[:, np.newaxis, :] - positions[np.newaxis, :, :] + dist_sq = np.sum(delta*delta, axis=2) + weighted = amplitudes[:, np.newaxis]*np.exp(-dist_sq*inv_two_sigma_sq) + + clusters = _merge_dominated_clusters(clusters, weighted, n_peaks) + + # Extract every cluster into its own source; skip footprints that did not + # split into more than one cluster. + if len(clusters) < 2: + return None + + # Assign every footprint pixel to the cluster of highest influence. + bbox = footprint.getBBox() + min_x, min_y = bbox.getMinX(), bbox.getMinY() + footprint_mask = _spans_to_bool(footprint.getSpans(), bbox) + xx, yy = np.meshgrid(np.arange(min_x, bbox.getMaxX() + 1), + np.arange(min_y, bbox.getMaxY() + 1)) + influence = np.empty((len(clusters), bbox.getHeight(), bbox.getWidth()), dtype=float) + for index, cluster in enumerate(clusters): + accumulator = np.zeros(footprint_mask.shape, dtype=float) + for peak in cluster: + r_sq = (xx - positions[peak, 0])**2 + (yy - positions[peak, 1])**2 + accumulator += amplitudes[peak]*np.exp(-r_sq*inv_two_sigma_sq) + influence[index] = accumulator + labels = np.argmax(influence, axis=0) + # Guarantee each cluster owns its own peak + for index, cluster in enumerate(clusters): + for peak in cluster: + labels[positions[peak, 1] - min_y, positions[peak, 0] - min_x] = index + + # One footprint per cluster, carrying that cluster's peaks. + region = footprint.getRegion() + footprints = [] + for index, cluster in enumerate(clusters): + fp = _footprint_from_bool(footprint_mask & (labels == index), bbox, peak_schema, region) + for peak in cluster: + fp.getPeaks().append(peaks[peak]) + footprints.append(fp) + + # The first cluster stays on the original record; the rest are new records. + return _DeblendedFootprint(parent=footprints[0], children=footprints[1:]) + + +class DeblendDifferenceFootprintsConfig(pexConfig.Config): + """Configuration for `DeblendDifferenceFootprintsTask`.""" + + clusterRadius = pexConfig.Field( + dtype=float, + default=3.0, + doc="Radius (in PSF FWHM) within which detected peaks are grouped into " + "one cluster. Each cluster is extracted into its own diaSource." + ) + maxFootprintArea = pexConfig.Field( + dtype=int, + default=50000, + doc="Footprints whose bounding-box area (pixels) exceeds this are left " + "un-deblended, bounding the per-footprint cost so a single large " + "footprint cannot dominate the processing time." + ) + maxPeaks = pexConfig.Field( + dtype=int, + default=40, + doc="Footprints with more than this many peaks are left un-deblended, " + "bounding the per-footprint deblending cost." + ) + + +class DeblendDifferenceFootprintsTask(pipeBase.Task): + """Deblend merged difference-image footprints by clustering their peaks. + + Each footprint whose peaks split into more than one spatial cluster is + replaced by one footprint per cluster: the first cluster stays on the + original record and the rest become new diaSources. Footprint pixels are + assigned to the cluster of highest Gaussian-weighted influence. Footprints + that are too large or have too many peaks are left un-deblended. + """ + + ConfigClass = DeblendDifferenceFootprintsConfig + _DefaultName = "deblendDifferenceFootprints" + + def run(self, diaSources, difference): + """Deblend the multi-cluster footprints of a diaSource catalog. + + Parameters + ---------- + diaSources : `lsst.afw.table.SourceCatalog` + The merged diaSources; modified in place and possibly extended. + The schema must carry the ``is_deblended`` and ``is_negative`` + flags, and footprint peak catalogs the ``merge_peak_negative`` and + ``significance`` fields. + difference : `lsst.afw.image.ExposureF` + Difference image the footprints were detected on. Only its PSF is + used, to set the clustering and influence scale. + + Returns + ------- + result : `lsst.pipe.base.Struct` + Results as a struct with attributes: + + ``diaSources`` + The updated catalog (`lsst.afw.table.SourceCatalog`). + """ + sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius() + linkingRadius = self.config.clusterRadius*sigma*_SIGMA_TO_FWHM + isDeblendedKey = diaSources.schema.find("is_deblended").key + isNegativeKey = diaSources.schema.find("is_negative").key + + nFootprints = 0 + nChildren = 0 + nSkipped = 0 + # Snapshot the records up front so appending children below does not + # disturb this iteration. + for source in list(diaSources): + footprint = source.getFootprint() + if footprint is None or len(footprint.getPeaks()) < 2: + continue + clusters = _cluster_peak_indices(footprint.getPeaks(), linkingRadius) + if len(clusters) < 2: + continue + # Leave very large footprints un-deblended + if (footprint.getBBox().getArea() > self.config.maxFootprintArea + or len(footprint.getPeaks()) > self.config.maxPeaks): + nSkipped += 1 + continue + result = _deblend_footprint(footprint, clusters, sigma) + if result is None: + continue + self._extractRecords(diaSources, source, result, isDeblendedKey, isNegativeKey) + nFootprints += 1 + nChildren += len(result.children) + + self.metadata["nDeblendedFootprints"] = nFootprints + self.metadata["nDeblendedDiaSources"] = nChildren + self.metadata["nSkipped"] = nSkipped + self.log.info("Deblended %d footprints into %d additional diaSources (%d skipped)", + nFootprints, nChildren, nSkipped) + return pipeBase.Struct(diaSources=diaSources) + + def _extractRecords(self, diaSources, source, result, isDeblendedKey, isNegativeKey): + """Write one deblended footprint per cluster back to the catalog. + + The first cluster stays on ``source``; the rest become new records + appended to ``diaSources``. Every resulting record is flagged + ``is_deblended``, and ``is_negative`` is set when all of that + footprint's peaks are negative. + """ + records = [source] + [diaSources.addNew() for _ in result.children] + negativePeakKey = result.parent.getPeaks().getSchema().find(_NEGATIVE_PEAK_FLAG).key + for record, newFootprint in zip(records, [result.parent] + result.children): + record.setFootprint(newFootprint) + record.set(isDeblendedKey, True) + record.set(isNegativeKey, + all(peak.get(negativePeakKey) for peak in newFootprint.getPeaks())) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index c601d308..a81ff790 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -48,8 +48,7 @@ from lsst.utils.timer import timeMethod from . import DipoleFitTask -from .peakClustering import find_peak_clusters -from .clusterDeblend import deblend_clustered_peaks +from .deblendDifferenceFootprints import DeblendDifferenceFootprintsTask __all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] @@ -264,37 +263,21 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, ) doDeblend = pexConfig.Field( dtype=bool, - default=False, - doc="Deblend merged footprints by extracting isolated positive peaks " - "into their own diaSources? Only used if ``doMerge`` is set." + default=True, + doc="Deblend merged footprints by clustering their peaks and " + "extracting each cluster into its own diaSource? Only used if " + "``doMerge`` is set." ) growFootprint = pexConfig.Field( dtype=int, default=2, doc="Grow positive and negative footprints by this many pixels before merging" ) - peakClusterRadius = pexConfig.Field( - dtype=float, - default=3.0, - doc="Radius (in PSF FWHM) to cluster detected peaks within footprints. " - "Positive peaks that are isolated after clustering can be deblended " - "if doDeblend=True." - ) - maxDeblendFootprintArea = pexConfig.Field( - dtype=int, - default=250000, - doc="Footprints whose bounding-box area (pixels) exceeds this are left " - "un-deblended, bounding the per-footprint deblending cost so a " - "single large footprint cannot blow the processing time budget. " + deblend = pexConfig.ConfigurableField( + target=DeblendDifferenceFootprintsTask, + doc="Task to deblend merged footprints by clustering their peaks. " "Only used if doDeblend=True." ) - maxDeblendPeaks = pexConfig.Field( - dtype=int, - default=200, - doc="Footprints with more than this many peaks are left un-deblended, " - "bounding the per-footprint deblending cost. Only used if " - "doDeblend=True." - ) measurement = pexConfig.ConfigurableField( target=DipoleFitTask, doc="Task to measure sources on the difference image.", @@ -617,6 +600,8 @@ def __init__(self, **kwargs): self.schema.addField("is_deblended", "Flag", "The original detected footprint was deblended, producing multiple diaSources." "The deblended diaSources split the original footprint.") + if self.config.doDeblend: + self.makeSubtask("deblend") self.schema.addField("reliability", type="F", doc="Reliability score of the DiaSource") self.schema.addField("reliabilityVersion", type=str, size=7, doc="Version of the reliability model") @@ -986,9 +971,8 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor # re-order the peaks in each footprint by significance. self._reorderPeaksBySignificance(initialDiaSources) - # Extract isolated and purely positive peaks from large footprints if self.config.doDeblend: - initialDiaSources = self._deblendClusteredPeaks(initialDiaSources, difference) + initialDiaSources = self.deblend.run(initialDiaSources, difference).diaSources else: initialDiaSources = sources @@ -1080,91 +1064,6 @@ def _reorderPeaksBySignificance(self, diaSources): peaks.clear() source.getFootprint().setPeakCatalog(newPeaks) - def _deblendClusteredPeaks(self, diaSources, difference): - """Extract clusters of peaks from large footprints into new diaSources. - - For each footprint whose peaks split into more than one cluster, every - cluster is assigned its own footprint: the first stays on the original - record and the rest become new diaSources. All resulting records are - flagged ``is_deblended``, and ``is_negative`` is set from the polarity - of each new footprint's peaks. Footprints with a single cluster are - left untouched. - - Parameters - ---------- - diaSources : `lsst.afw.table.SourceCatalog` - The merged diaSources; modified in place and possibly extended. - difference : `lsst.afw.image.ExposureF` - Difference image the footprints were detected on, used to weight - the peaks. - - Returns - ------- - diaSources : `lsst.afw.table.SourceCatalog` - The updated, contiguous catalog. - """ - sigmaToFwhm = np.sqrt(8.0*np.log(2.0)) - linkingRadius = self.config.peakClusterRadius*self.sigma*sigmaToFwhm - # Identify clusters of peaks within each footprint. Footprints whose - # peaks split into more than one cluster can be deblended. - peakClusters = find_peak_clusters(diaSources, linkingRadius) - - # Precalculate some mask information up front, so the pseudo-deblending - # function can be as lightweight as possible - presentBadPlanes = [mp for mp in self.config.badMaskPlanes - if mp in difference.mask.getMaskPlaneDict()] - badBitmask = difference.mask.getPlaneBitMask(presentBadPlanes) if presentBadPlanes else 0 - imageArray = difference.image.array - maskArray = difference.mask.array - imageXY0 = (difference.getBBox().getMinX(), difference.getBBox().getMinY()) - isDeblendedKey = self.schema.find("is_deblended").key - isNegativeKey = self.schema.find("is_negative").key - - nFootprints = 0 - nChildren = 0 - nSkipped = 0 - # Loop over each individual footprint, checking whether its peaks split - # into clusters that can be extracted and measured independently - for footprintClusters in peakClusters: - if not footprintClusters.is_deblendable: - continue - # Leave pathologically large footprints un-deblended: the deblending - # cost grows steeply with footprint area and peak count, so a single - # huge footprint could otherwise dominate the processing time. - if (footprintClusters.footprint.getBBox().getArea() > self.config.maxDeblendFootprintArea - or footprintClusters.n_peaks > self.config.maxDeblendPeaks): - nSkipped += 1 - continue - result = deblend_clustered_peaks(footprintClusters, imageArray, maskArray, - imageXY0, self.sigma, badBitmask) - # Exit early without modifying the footprint catalog if no changes - # are needed for this footprint - if result is None: - continue - # The first cluster stays on the original record; the rest are new. - records = [footprintClusters.source] - records += [diaSources.addNew() for _ in result.children] - negativePeakKey = result.parent.getPeaks().getSchema().find("merge_peak_negative").key - for record, newFootprint in zip(records, [result.parent] + result.children): - record.setFootprint(newFootprint) - record.set(isDeblendedKey, True) - # A source is negative if all of its peaks are negative. - record.set(isNegativeKey, all(peak.get(negativePeakKey) for peak in newFootprint.getPeaks())) - nFootprints += 1 - nChildren += len(result.children) - - self.metadata["nClusterDeblendedFootprints"] = nFootprints - self.metadata["nClusterDeblendedDiaSources"] = nChildren - self.metadata["nClusterDeblendSkipped"] = nSkipped - if nFootprints: - self.log.info("Deblended %d footprints into %d diaSources", nFootprints, nChildren) - if nSkipped: - self.log.info("Left %d oversized footprints un-deblended (area > %d or peaks > %d)", - nSkipped, self.config.maxDeblendFootprintArea, self.config.maxDeblendPeaks) - if not diaSources.isContiguous(): - diaSources = diaSources.copy(deep=True) - return diaSources - def _rejectBadMaskedDetections(self, initialDiaSources, mask): """Remove detections whose footprint peak lies on a bad-masked pixel. diff --git a/python/lsst/ip/diffim/peakClustering.py b/python/lsst/ip/diffim/peakClustering.py deleted file mode 100644 index 27ef01e2..00000000 --- a/python/lsst/ip/diffim/peakClustering.py +++ /dev/null @@ -1,230 +0,0 @@ -# This file is part of ip_diffim. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -"""Identify spatial clusters of peaks within difference image footprints. -""" - -__all__ = ["PeakCluster", "FootprintPeakClusters", "cluster_peaks", - "find_peak_clusters"] - -import dataclasses - -import numpy as np -from scipy.sparse import coo_matrix -from scipy.sparse.csgraph import connected_components -from scipy.spatial import cKDTree - -import lsst.afw.detection -import lsst.afw.table - - -# Per-peak flag set by `FootprintMergeList` marking an individual peak as -# originating from the negative-polarity detection footprint set. -_NEGATIVE_PEAK_FLAG = "merge_peak_negative" - - -@dataclasses.dataclass -class PeakCluster: - """A group of peaks within a single footprint that lie close together. - - Peaks are grouped by single-linkage clustering on their positions, - ignoring polarity, so a dipole's positive and negative lobes should fall - into the same cluster. - - Parameters - ---------- - peak_indices : `list` [`int`] - Rows into the parent footprint's peak catalog - (``footprint.getPeaks()``) that belong to this cluster. - positions : `numpy.ndarray` - ``(n, 2)`` array of integer peak-pixel positions ``(i_x, i_y)`` in - the parent (image) coordinate frame. - peak_values : `numpy.ndarray` - ``(n,)`` array of the signed peak values. - is_negative : `numpy.ndarray` - ``(n,)`` boolean array, `True` for peaks from the negative detection - set (from the ``merge_peak_negative`` flag). - """ - - peak_indices: list - positions: np.ndarray - peak_values: np.ndarray - is_negative: np.ndarray - - -@dataclasses.dataclass -class FootprintPeakClusters: - """The peak clusters found within one multi-peak footprint. - - Parameters - ---------- - source : `lsst.afw.table.SourceRecord` - The source record whose footprint was clustered. This is the live - record from the catalog, so a downstream deblender can modify it (and - its footprint) in place without a separate lookup. - clusters : `list` [`PeakCluster`] - The clusters of peaks within the footprint. - """ - - source: lsst.afw.table.SourceRecord - clusters: list - - @property - def footprint(self): - """The footprint whose peaks were clustered - (`lsst.afw.detection.Footprint`). - """ - return self.source.getFootprint() - - @property - def source_id(self): - """Id of the source record (`int`). - - Only meaningful once ids have been assigned to the catalog; while - clustering runs inside the merge step, ids are not yet set. - """ - return self.source.getId() - - @property - def n_peaks(self): - """Total number of peaks in the footprint (`int`).""" - return sum(len(cluster.peak_indices) for cluster in self.clusters) - - @property - def is_deblendable(self): - """Whether the footprint splits into more than one cluster, so that - deblending would extract more than one source (`bool`). - """ - return len(self.clusters) >= 2 - - -def _single_linkage_labels(positions, linking_radius): - """Assign a cluster label to each position by single-linkage clustering. - - Two positions are linked when they lie within ``linking_radius`` of each - other; clusters are the connected components of the resulting graph. - - Parameters - ---------- - positions : `numpy.ndarray` - ``(n, 2)`` array of positions. - linking_radius : `float` - Maximum distance, in pixels, between two peaks for them to be linked - into the same cluster. - - Returns - ------- - labels : `numpy.ndarray` - ``(n,)`` integer array of cluster labels. - """ - n = len(positions) - if n <= 1: - return np.zeros(n, dtype=int) - pairs = cKDTree(positions).query_pairs(r=linking_radius, output_type="ndarray") - if len(pairs) == 0: - graph = coo_matrix((n, n)) - else: - graph = coo_matrix((np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n)) - _, labels = connected_components(graph, directed=False) - return labels - - -def cluster_peaks(peaks, linking_radius): - """Group the peaks of one footprint into spatial clusters. - - Peaks are clustered by single-linkage on their positions, ignoring - polarity: two peaks within ``linking_radius`` of each other join the same - cluster, and clustering is transitive (chaining). A peak with no - neighbour within the radius forms its own singleton cluster. - - Parameters - ---------- - peaks : `lsst.afw.detection.PeakCatalog` - Peaks of a single (merged) footprint. The schema must carry the - ``merge_peak_negative`` flag. - linking_radius : `float` - Maximum distance, in pixels, between two peaks for them to be - clustered together. - - Returns - ------- - clusters : `list` [`PeakCluster`] - One entry per cluster. Empty if ``peaks`` is empty. - """ - n = len(peaks) - positions = np.empty((n, 2), dtype=float) - peak_values = np.empty(n, dtype=float) - is_negative = np.empty(n, dtype=bool) - negative_key = peaks.getSchema().find(_NEGATIVE_PEAK_FLAG).key - for i, peak in enumerate(peaks): - # Detection peaks store the integer local-maximum pixel; f_x/f_y equal - # i_x/i_y here, so use the integer accessors to be explicit. - positions[i] = (peak.getIx(), peak.getIy()) - peak_values[i] = peak.getPeakValue() - is_negative[i] = peak.get(negative_key) - - labels = _single_linkage_labels(positions, linking_radius) - clusters = [] - for label in np.unique(labels): - (members,) = np.nonzero(labels == label) - clusters.append(PeakCluster( - peak_indices=members.tolist(), - positions=positions[members], - peak_values=peak_values[members], - is_negative=is_negative[members], - )) - return clusters - - -def find_peak_clusters(catalog, linking_radius, *, min_peaks=2): - """Find peak clusters in every multi-peak footprint of a catalog. - - Parameters - ---------- - catalog : `lsst.afw.table.SourceCatalog` - Sources whose footprints are to be searched, e.g. the merged - ``initialDiaSources``. Footprint peak catalogs must carry the - ``merge_peak_negative`` flag. - linking_radius : `float` - Maximum distance, in pixels, between two peaks for them to be - clustered together. - min_peaks : `int`, optional - Only footprints with at least this many peaks are considered; those - with fewer are not deblending candidates and are skipped. - - Returns - ------- - results : `list` [`FootprintPeakClusters`] - One entry per footprint with at least ``min_peaks`` peaks. - """ - results = [] - for source in catalog: - footprint = source.getFootprint() - if footprint is None: - continue - peaks = footprint.getPeaks() - if len(peaks) < min_peaks: - continue - results.append(FootprintPeakClusters( - source=source, - clusters=cluster_peaks(peaks, linking_radius), - )) - return results diff --git a/tests/test_clusterDeblend.py b/tests/test_clusterDeblend.py deleted file mode 100644 index 61c86760..00000000 --- a/tests/test_clusterDeblend.py +++ /dev/null @@ -1,214 +0,0 @@ -# This file is part of ip_diffim. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -import unittest - -import numpy as np - -import lsst.afw.detection as afwDetection -import lsst.afw.geom as afwGeom -import lsst.afw.table as afwTable -import lsst.geom as geom -import lsst.utils.tests -from lsst.ip.diffim.clusterDeblend import deblend_clustered_peaks -from lsst.ip.diffim.peakClustering import find_peak_clusters - -SIZE = 40 - - -def makeSchemas(): - schema = afwTable.SourceTable.makeMinimalSchema() - mergeList = afwDetection.FootprintMergeList(schema, ["positive", "negative"]) - return schema, mergeList.getPeakSchema() - - -def makeFootprintAndImage(peakSchema, peaks, ampValue=100.0, badPeaks=()): - """Build a full-frame footprint with the given peaks, and a matching image. - - Parameters - ---------- - peaks : `list` of `tuple` - Each ``(ix, iy, isNegative)``. - ampValue : `float` - Magnitude placed in the 3x3 block around each peak (negated for - negative peaks) in the difference image. - badPeaks : `set` of `int` - Indices of peaks whose 3x3 block is flagged bad in the mask. - """ - box = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(SIZE, SIZE)) - footprint = afwDetection.Footprint(afwGeom.SpanSet(box), peakSchema) - negativeKey = peakSchema.find("merge_peak_negative").key - positiveKey = peakSchema.find("merge_peak_positive").key - - image = np.zeros((SIZE, SIZE), dtype=np.float32) - mask = np.zeros((SIZE, SIZE), dtype=np.int32) - badBit = 1 - for i, (ix, iy, isNegative) in enumerate(peaks): - peak = footprint.addPeak(ix, iy, -ampValue if isNegative else ampValue) - peak.set(negativeKey, isNegative) - peak.set(positiveKey, not isNegative) - block = (slice(iy - 1, iy + 2), slice(ix - 1, ix + 2)) - image[block] = -ampValue if isNegative else ampValue - if i in badPeaks: - mask[block] = badBit - return footprint, image, mask, badBit - - -def makeClusters(schema, footprint, linkingRadius): - catalog = afwTable.SourceCatalog(schema) - record = catalog.addNew() - record.setId(1) - record.setFootprint(footprint) - results = find_peak_clusters(catalog, linkingRadius) - return results[0] if results else None - - -class ClusterDeblendTestCase(lsst.utils.tests.TestCase): - - def setUp(self): - self.schema, self.peakSchema = makeSchemas() - - def test_extract_clusters_from_footprint(self): - """A footprint that splits into two clusters yields two sources.""" - # Lone peak at (8, 20); dipole at (30, 20)/(32, 20). - footprint, image, mask, _ = makeFootprintAndImage( - self.peakSchema, [(8, 20, False), (30, 20, False), (32, 20, True)]) - fc = makeClusters(self.schema, footprint, linkingRadius=5.0) - self.assertTrue(fc.is_deblendable) - - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) - self.assertIsNotNone(result) - self.assertEqual(len(result.children), 1) - - # The first cluster (the lone peak at (8, 20)) stays on the parent. - self.assertEqual(len(result.parent.getPeaks()), 1) - self.assertEqual((result.parent.getPeaks()[0].getIx(), - result.parent.getPeaks()[0].getIy()), (8, 20)) - - # The dipole is the other cluster, keeping both its peaks. - child = result.children[0] - self.assertEqual(len(child.getPeaks()), 2) - self.assertTrue(child.getSpans().contains(geom.Point2I(30, 20))) - self.assertTrue(child.getSpans().contains(geom.Point2I(32, 20))) - - # The extracted pixels exactly tile the original footprint. - self.assertEqual(result.parent.getSpans().getArea() + child.getSpans().getArea(), - footprint.getSpans().getArea()) - # Peak schema (with merge flags) is preserved on the outputs. - self.assertIn("merge_peak_negative", child.getPeaks().schema.getNames()) - - def test_splits_beyond_radius_positives(self): - """Two positive peaks beyond the linking radius split into two sources - even with no remainder cluster. - """ - footprint, image, mask, _ = makeFootprintAndImage( - self.peakSchema, [(8, 20, False), (32, 20, False)]) - fc = makeClusters(self.schema, footprint, linkingRadius=5.0) - self.assertTrue(fc.is_deblendable) - - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) - self.assertIsNotNone(result) - # One peak stays on the parent, the other becomes a new source. - self.assertEqual(len(result.children), 1) - self.assertEqual(len(result.parent.getPeaks()), 1) - self.assertEqual(len(result.children[0].getPeaks()), 1) - self.assertEqual(result.parent.getSpans().getArea() - + result.children[0].getSpans().getArea(), - footprint.getSpans().getArea()) - - def test_within_radius_positive_blend_stays(self): - """Two positive peaks within the linking radius are one cluster and - are left as a single source. - """ - footprint, image, mask, _ = makeFootprintAndImage( - self.peakSchema, [(20, 20, False), (22, 20, False)]) - fc = makeClusters(self.schema, footprint, linkingRadius=5.0) - self.assertFalse(fc.is_deblendable) - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) - self.assertIsNone(result) - - def test_single_cluster_dipole_returns_none(self): - """A footprint that is a single cluster (here a dipole) is untouched.""" - footprint, image, mask, _ = makeFootprintAndImage( - self.peakSchema, [(20, 20, False), (22, 20, True)]) - fc = makeClusters(self.schema, footprint, linkingRadius=5.0) - self.assertFalse(fc.is_deblendable) - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) - self.assertIsNone(result) - - def test_extract_negative_cluster(self): - """Negative and positive clusters are both extracted (polarity is not - a constraint in this version). - """ - # Lone negative at (8, 20); lone positive at (32, 20); beyond radius. - footprint, image, mask, _ = makeFootprintAndImage( - self.peakSchema, [(8, 20, True), (32, 20, False)]) - fc = makeClusters(self.schema, footprint, linkingRadius=5.0) - self.assertTrue(fc.is_deblendable) - - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=0) - self.assertIsNotNone(result) - self.assertEqual(len(result.children), 1) - # Both peaks are extracted, one per output footprint. - outputs = [result.parent] + result.children - extracted = {(fp.getPeaks()[0].getIx(), fp.getPeaks()[0].getIy()) for fp in outputs} - self.assertEqual(extracted, {(8, 20), (32, 20)}) - - def test_dominated_isolated_peak_is_merged(self): - """A faint 'isolated' peak dominated by a bright neighbour is merged - back in, leaving nothing to extract. - """ - # Faint positive at (12, 20); bright dipole at (16, 20)/(18, 20). - peaks = [(12, 20, False), (16, 20, False), (18, 20, True)] - footprint, image, mask, _ = makeFootprintAndImage(self.peakSchema, peaks, ampValue=100.0) - # Make the first peak faint: overwrite its 3x3 block with a small value. - image[19:22, 11:14] = 5.0 - fc = makeClusters(self.schema, footprint, linkingRadius=3.0) - self.assertTrue(fc.is_deblendable) - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=2.0, bad_bitmask=0) - self.assertIsNone(result) - - def test_bad_masked_amplitude_excluded(self): - """A peak whose 3x3 block is entirely bad-masked gets zero amplitude - and is merged away (proving bad pixels are excluded from the sum). - """ - footprint, image, mask, badBit = makeFootprintAndImage( - self.peakSchema, [(8, 20, False), (30, 20, False), (32, 20, True)], - ampValue=100.0, badPeaks={0}) - fc = makeClusters(self.schema, footprint, linkingRadius=5.0) - self.assertTrue(fc.is_deblendable) - # With the bright isolated peak's flux excluded, it has no influence - # and cannot be extracted. - result = deblend_clustered_peaks(fc, image, mask, (0, 0), sigma=4.0, bad_bitmask=badBit) - self.assertIsNone(result) - - -class MemoryTester(lsst.utils.tests.MemoryTestCase): - pass - - -def setup_module(module): - lsst.utils.tests.init() - - -if __name__ == "__main__": - lsst.utils.tests.init() - unittest.main() diff --git a/tests/test_deblendDifferenceFootprints.py b/tests/test_deblendDifferenceFootprints.py new file mode 100644 index 00000000..d84d34e3 --- /dev/null +++ b/tests/test_deblendDifferenceFootprints.py @@ -0,0 +1,308 @@ +# This file is part of ip_diffim. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import unittest + +import numpy as np + +import lsst.afw.detection as afwDetection +import lsst.afw.geom as afwGeom +import lsst.afw.image as afwImage +import lsst.afw.table as afwTable +import lsst.geom as geom +import lsst.utils.tests +from lsst.ip.diffim.deblendDifferenceFootprints import ( + DeblendDifferenceFootprintsTask, + _cluster_peak_indices, + _deblend_footprint, +) + + +def makeSchemas(): + """Build a source schema with the deblending flags and the matching merged + peak schema. + """ + schema = afwTable.SourceTable.makeMinimalSchema() + schema.addField("is_negative", type="Flag", doc="Source is negative.") + schema.addField("is_deblended", type="Flag", doc="Source was deblended.") + mergeList = afwDetection.FootprintMergeList(schema, ["positive", "negative"]) + peakSchema = mergeList.getPeakSchema() + peakSchema.addField("significance", type=np.float64, + doc="Ratio of peak value to the standard deviation.") + return schema, peakSchema + + +def makeFootprint(peakSchema, peaks): + """Build a full-frame footprint carrying the given peaks. + + Parameters + ---------- + peaks : `list` of `tuple` + Each ``(ix, iy, isNegative, significance)``. + """ + box = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64)) + footprint = afwDetection.Footprint(afwGeom.SpanSet(box), peakSchema) + negativeKey = peakSchema.find("merge_peak_negative").key + positiveKey = peakSchema.find("merge_peak_positive").key + significanceKey = peakSchema.find("significance").key + for ix, iy, isNegative, significance in peaks: + peak = footprint.addPeak(ix, iy, -significance if isNegative else significance) + peak.set(negativeKey, isNegative) + peak.set(positiveKey, not isNegative) + peak.set(significanceKey, significance) + return footprint + + +class ClusterPeaksTestCase(lsst.utils.tests.TestCase): + """Tests of the initial single-linkage peak clustering.""" + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + self.radius = 5.0 + + def _cluster(self, peaks): + # peaks: (ix, iy, isNegative); significance is irrelevant to clustering. + footprint = makeFootprint(self.peakSchema, [(ix, iy, neg, 100.0) for ix, iy, neg in peaks]) + return _cluster_peak_indices(footprint.getPeaks(), self.radius) + + def test_two_peaks_far_apart(self): + """Peaks beyond the linking radius form separate singleton clusters.""" + clusters = self._cluster([(10, 10, False), (30, 30, False)]) + self.assertEqual(sorted(clusters), [[0], [1]]) + + def test_two_peaks_close(self): + """Peaks within the linking radius merge into one cluster.""" + clusters = self._cluster([(10, 10, False), (12, 10, False)]) + self.assertEqual(len(clusters), 1) + self.assertEqual(sorted(clusters[0]), [0, 1]) + + def test_transitive_chain(self): + """Single-linkage chains A-B-C even when A and C exceed the radius.""" + # 10->14->18: adjacent gaps are 4 (< 5), but 10->18 is 8 (> 5). + A = (10, 10, False) + B = (14, 10, False) + C = (18, 10, False) + + # A and C alone should not merge + clustersAC = self._cluster([A, C]) + self.assertEqual(len(clustersAC), 2) + clustersABC = self._cluster([A, B, C]) + self.assertEqual(len(clustersABC), 1) + self.assertEqual(sorted(clustersABC[0]), [0, 1, 2]) + + def test_dipole_single_cluster(self): + """A close positive/negative pair clusters together.""" + clusters = self._cluster([(10, 10, False), (12, 10, True)]) + self.assertEqual(len(clusters), 1) + self.assertEqual(sorted(clusters[0]), [0, 1]) + + def test_lone_peak_beside_dipole(self): + """A lone peak is separated from a nearby dipole into two clusters.""" + clusters = self._cluster([(50, 50, False), # lone peak + (10, 10, False), (15, 10, True)]) # dipole + lone = [c for c in clusters if 0 in c] + dipole = [c for c in clusters if 1 in c] + self.assertEqual(len(clusters), 2) + self.assertEqual(lone[0], [0]) + self.assertEqual(sorted(dipole[0]), [1, 2]) + + +class DeblendFootprintTestCase(lsst.utils.tests.TestCase): + """Tests of the per-footprint deblending routine.""" + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + + def _deblend(self, peaks, linkingRadius=5.0, sigma=4.0): + footprint = makeFootprint(self.peakSchema, peaks) + + clusters = _cluster_peak_indices(footprint.getPeaks(), linkingRadius) + result = _deblend_footprint(footprint, clusters, sigma=sigma) + return footprint, clusters, result + + def test_extract_clusters_from_footprint(self): + """A footprint that splits into two clusters yields two footprints.""" + # Lone peak at (8, 20); dipole at (30, 20)/(35, 20). + footprint, clusters, result = self._deblend( + [(8, 20, False, 100.0), (30, 20, False, 100.0), (35, 20, True, 100.0)]) + self.assertEqual(len(clusters), 2) + self.assertIsNotNone(result) + self.assertEqual(len(result.children), 1) + + # The first cluster (the lone peak at (8, 20)) stays on the parent. + self.assertEqual(len(result.parent.getPeaks()), 1) + self.assertEqual((result.parent.getPeaks()[0].getIx(), + result.parent.getPeaks()[0].getIy()), (8, 20)) + + # The dipole is the other cluster, keeping both its peaks. + child = result.children[0] + self.assertEqual(len(child.getPeaks()), 2) + self.assertTrue(child.getSpans().contains(geom.Point2I(30, 20))) + self.assertTrue(child.getSpans().contains(geom.Point2I(32, 20))) + + # The extracted pixels exactly tile the original footprint. + self.assertEqual(result.parent.getSpans().getArea() + child.getSpans().getArea(), + footprint.getSpans().getArea()) + # Peak schema (with merge flags) is preserved on the outputs. + self.assertIn("merge_peak_negative", child.getPeaks().schema.getNames()) + + def test_splits_beyond_radius(self): + """Two peaks beyond the linking radius split into two footprints.""" + _, _, result = self._deblend([(8, 20, False, 100.0), (32, 20, False, 100.0)]) + self.assertIsNotNone(result) + self.assertEqual(len(result.children), 1) + self.assertEqual(len(result.parent.getPeaks()), 1) + self.assertEqual(len(result.children[0].getPeaks()), 1) + + def test_within_radius_blend_stays(self): + """Two peaks within the linking radius are one cluster: no split.""" + footprint, clusters, result = self._deblend([(20, 20, False, 100.0), (22, 20, False, 100.0)]) + self.assertEqual(len(clusters), 1) + self.assertIsNone(result) + + def test_single_cluster_dipole_returns_none(self): + """A footprint that is a single cluster (a dipole) is untouched.""" + footprint, clusters, result = self._deblend([(20, 20, False, 100.0), (22, 20, True, 100.0)]) + self.assertEqual(len(clusters), 1) + self.assertIsNone(result) + + def test_extract_negative_cluster(self): + """Negative and positive clusters are both extracted.""" + _, _, result = self._deblend([(8, 20, True, 100.0), (32, 20, False, 100.0)]) + self.assertIsNotNone(result) + self.assertEqual(len(result.children), 1) + outputs = [result.parent] + result.children + extracted = {(fp.getPeaks()[0].getIx(), fp.getPeaks()[0].getIy()) for fp in outputs} + self.assertEqual(extracted, {(8, 20), (32, 20)}) + + def test_low_significance_peak_is_merged(self): + """A faint (low-significance) peak dominated by a bright neighbour is + merged back in, leaving nothing to extract. + """ + # Faint positive at (12, 20); bright dipole at (16, 20)/(18, 20). + _, _, result = self._deblend( + [(12, 20, False, 5.0), (16, 20, False, 100.0), (18, 20, True, 100.0)], + linkingRadius=3.0, sigma=2.0) + self.assertIsNone(result) + + +class RunTestCase(lsst.utils.tests.TestCase): + """Tests of the task-level ``run`` over a diaSource catalog.""" + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + + def makeExposure(self, psfSigma=2.0): + # Only the PSF is used by the task, to set the clustering scale. + bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64)) + exposure = afwImage.ExposureF(bbox) + exposure.setPsf(afwDetection.GaussianPsf(11, 11, psfSigma)) + return exposure + + def makeCatalog(self, footprints): + catalog = afwTable.SourceCatalog(self.schema) + for footprint in footprints: + catalog.addNew().setFootprint(footprint) + return catalog + + def makeTask(self, **configKwargs): + config = DeblendDifferenceFootprintsTask.ConfigClass() + # sigma ~= 2 px, FWHM ~= 4.7 px, so clusterRadius=2 links within ~9 px. + config.clusterRadius = 2.0 + for key, value in configKwargs.items(): + setattr(config, key, value) + return DeblendDifferenceFootprintsTask(config=config) + + def test_run_splits_multicluster_footprint(self): + """A footprint with two well-separated peaks becomes two diaSources.""" + footprint = makeFootprint(self.peakSchema, [(10, 20, False, 100.0), (50, 20, False, 100.0)]) + catalog = self.makeCatalog([footprint]) + task = self.makeTask() + + result = task.run(catalog, self.makeExposure()) + + self.assertEqual(len(result.diaSources), 2) + self.assertTrue(all(record["is_deblended"] for record in result.diaSources)) + self.assertEqual(task.metadata["nDeblendedFootprints"], 1) + self.assertEqual(task.metadata["nDeblendedDiaSources"], 1) + self.assertEqual(task.metadata["nSkipped"], 0) + + def test_run_leaves_single_cluster_untouched(self): + """A footprint that is one cluster is not modified or flagged.""" + footprint = makeFootprint(self.peakSchema, [(20, 20, False, 100.0), (22, 20, False, 100.0)]) + catalog = self.makeCatalog([footprint]) + task = self.makeTask() + + result = task.run(catalog, self.makeExposure()) + + self.assertEqual(len(result.diaSources), 1) + self.assertFalse(result.diaSources[0]["is_deblended"]) + self.assertEqual(task.metadata["nDeblendedFootprints"], 0) + + def test_run_sets_is_negative(self): + """An extracted purely-negative cluster is flagged is_negative.""" + footprint = makeFootprint(self.peakSchema, [(10, 20, False, 100.0), (50, 20, True, 100.0)]) + catalog = self.makeCatalog([footprint]) + task = self.makeTask() + + result = task.run(catalog, self.makeExposure()) + + self.assertEqual(len(result.diaSources), 2) + byNegative = sorted(record["is_negative"] for record in result.diaSources) + self.assertEqual(byNegative, [False, True]) + + def test_run_skips_too_many_peaks(self): + """A deblendable footprint with more peaks than maxPeaks is skipped.""" + footprint = makeFootprint( + self.peakSchema, [(10, 20, False, 100.0), (30, 20, False, 100.0), (50, 20, False, 100.0)]) + catalog = self.makeCatalog([footprint]) + task = self.makeTask(maxPeaks=2) + + result = task.run(catalog, self.makeExposure()) + + self.assertEqual(len(result.diaSources), 1) + self.assertFalse(result.diaSources[0]["is_deblended"]) + self.assertEqual(task.metadata["nSkipped"], 1) + + def test_run_skips_too_large_area(self): + """A deblendable footprint bigger than maxFootprintArea is skipped.""" + footprint = makeFootprint(self.peakSchema, [(10, 20, False, 100.0), (50, 20, False, 100.0)]) + catalog = self.makeCatalog([footprint]) + task = self.makeTask(maxFootprintArea=100) # footprint area is SIZE**2 + + result = task.run(catalog, self.makeExposure()) + + self.assertEqual(len(result.diaSources), 1) + self.assertFalse(result.diaSources[0]["is_deblended"]) + self.assertEqual(task.metadata["nSkipped"], 1) + + +class MemoryTester(lsst.utils.tests.MemoryTestCase): + pass + + +def setup_module(module): + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main() diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index f48f67b7..04c074e8 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -126,7 +126,7 @@ def _check_values(self, values, minValue=None, maxValue=None): def _setup_detection(self, doSkySources=True, nSkySources=5, doSubtractBackground=False, run_sattle=False, - badSourceFlags=None, **kwargs): + badSourceFlags=None, clusterRadius=None, **kwargs): """Setup and configure the detection and measurement PipelineTask. Parameters @@ -148,6 +148,8 @@ def _setup_detection(self, doSkySources=True, nSkySources=5, if doSkySources: config.skySources.nSources = nSkySources config.update(**kwargs) + if clusterRadius is not None: + config.deblend.clusterRadius = clusterRadius config.run_sattle = run_sattle @@ -193,7 +195,7 @@ def test_detection_xy0(self): # Set the subtraction residual metric threshold high, since the template # is not actually subtracted from the science image. detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., - doSkySources=False, peakClusterRadius=1.0) + doSkySources=False, clusterRadius=1.5) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, sources, @@ -826,7 +828,7 @@ def _setup_sattle_tests(self): detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., doSkySources=False, run_sattle=True, - peakClusterRadius=1.0) + clusterRadius=1.5) return science, matchedTemplate, difference, sources, detectionTask @@ -960,7 +962,7 @@ def test_detection_xy0(self): # Set the subtraction residual metric threshold high, since the template # is not actually subtracted from the science image. detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., - doSkySources=False, peakClusterRadius=1.0) + doSkySources=False, clusterRadius=1.5) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, score, sources, @@ -1493,7 +1495,7 @@ def testMergeFootprints(self): # Two positive sources should not have `is_negative` set: the isolated # source at (25, 26) and the positive blended pair, which stays merged - # as a single source because its peaks are within ``peakClusterRadius``. + # as a single source because its peaks are within ``clusterRadius``. self.assertEqual((~result.diaSources["is_negative"]).sum(), 2) def testReorderPeaksBySignificance(self): @@ -1531,6 +1533,56 @@ def testReorderPeaksBySignificance(self): self.assertEqual(peaks[1].getPeakValue(), 30.0) self.assertEqual(peaks[1]["significance"], 5.0) + def testMergedPeakSignificance(self): + """Merged diaSource peaks carry the signal-to-noise of the smoothed + difference image at the peak, positive for negative peaks too. + """ + bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(200, 100)) + schema = lsst.meas.base.tests.TestDataset.makeMinimalSchema() + schema.addField("sky_source", "Flag", "Sky source.") + science, catalog = lsst.meas.base.tests.TestDataset(bbox).realize(noise=10.0, schema=schema) + template, _ = lsst.meas.base.tests.TestDataset(bbox).realize(noise=10.0, schema=schema) + # A noise-free difference image with a uniform variance plane, so the + # expected significance follows from the source fluxes alone. + difference = science.clone() + psf = difference.psf + sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() + variance = 25.0 + difference.image.array[:] = 0.0 + difference.variance.array[:] = variance + # (x, y, flux) of an isolated positive source, an isolated negative + # source, and a dipole whose lobes merge into one footprint. + sources = [(40, 30, 4000.0), (100, 30, -6000.0), (150, 60, 5000.0), (158, 60, -3000.0)] + for x, y, flux in sources: + stamp = psf.computeImage(lsst.geom.Point2D(x, y)) + difference.image[stamp.getBBox()].array[:] += flux*stamp.array + + config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() + config.doSkySources = False + config.raiseOnBadSubtractionRatio = False + config.raiseOnNoDiaSources = False + task = detectAndMeasure.DetectAndMeasureTask(config=config) + diaSources = task.run(science, template, difference, catalog).diaSources + + def expected(ix, iy): + # Detection smooths with a Gaussian of the PSF width, which turns + # each source into a Gaussian of variance 2 sigma^2 and divides the + # white noise by 2 sqrt(pi) sigma. + value = sum(flux/(4*np.pi*sigma**2)*np.exp(-((ix - x)**2 + (iy - y)**2)/(4*sigma**2)) + for x, y, flux in sources) + return abs(value)*2*np.sqrt(np.pi)*sigma/np.sqrt(variance) + + self.assertIn("significance", diaSources[0].getFootprint().peaks.schema) + # The dipole's lobes share one merged footprint. + self.assertEqual(sorted(len(source.getFootprint().peaks) for source in diaSources), [1, 1, 2]) + peaks = [peak for source in diaSources for peak in source.getFootprint().peaks] + self.assertEqual(sum(peak["merge_peak_negative"] for peak in peaks), 2) + for peak in peaks: + with self.subTest(x=peak.getIx(), y=peak.getIy()): + self.assertGreater(peak["significance"], 0) + self.assertFloatsAlmostEqual(peak["significance"], expected(peak.getIx(), peak.getIy()), + rtol=1e-3) + def _makeDeblendableCatalog(self, schema, bbox, peakXY): """A single-source catalog whose footprint has well-separated peaks.""" peakSchema = afwDetection.FootprintMergeList( @@ -1550,30 +1602,34 @@ def _makeDeblendableCatalog(self, schema, bbox, peakXY): def testDeblendGuardrailSkipsManyPeaks(self): """A deblendable footprint with more peaks than the cap is left whole.""" config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() - config.maxDeblendPeaks = 2 + config.doDeblend = True + config.deblend.maxPeaks = 2 task = detectAndMeasure.DetectAndMeasureTask(config=config) - task.sigma = 1.0 + # psfSize=1 keeps the linking radius small enough that the peaks below + # split into separate clusters, so the footprint is deblendable. + difference, _ = makeTestImage(psfSize=1., x0=0, y0=0) bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(40, 40)) # Three peaks, spaced beyond the linking radius, so the footprint is - # deblendable but exceeds maxDeblendPeaks. + # deblendable but exceeds maxPeaks. diaSources = self._makeDeblendableCatalog(task.schema, bbox, [(8, 20), (20, 20), (32, 20)]) - result = task._deblendClusteredPeaks(diaSources, afwImage.ExposureF(bbox)) + result = task.deblend.run(diaSources, difference).diaSources self.assertEqual(len(result), 1) self.assertFalse(result[0]["is_deblended"]) - self.assertEqual(task.metadata["nClusterDeblendSkipped"], 1) + self.assertEqual(task.deblend.metadata["nSkipped"], 1) def testDeblendGuardrailSkipsLargeArea(self): """A deblendable footprint bigger than the area cap is left whole.""" config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() - config.maxDeblendFootprintArea = 100 + config.doDeblend = True + config.deblend.maxFootprintArea = 100 task = detectAndMeasure.DetectAndMeasureTask(config=config) - task.sigma = 1.0 + difference, _ = makeTestImage(psfSize=1., x0=0, y0=0) bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(40, 40)) # area 1600 > 100 diaSources = self._makeDeblendableCatalog(task.schema, bbox, [(8, 20), (20, 20), (32, 20)]) - result = task._deblendClusteredPeaks(diaSources, afwImage.ExposureF(bbox)) + result = task.deblend.run(diaSources, difference).diaSources self.assertEqual(len(result), 1) self.assertFalse(result[0]["is_deblended"]) - self.assertEqual(task.metadata["nClusterDeblendSkipped"], 1) + self.assertEqual(task.deblend.metadata["nSkipped"], 1) def makeVisitInfo(id=1): diff --git a/tests/test_peakClustering.py b/tests/test_peakClustering.py deleted file mode 100644 index 18e87688..00000000 --- a/tests/test_peakClustering.py +++ /dev/null @@ -1,227 +0,0 @@ -# This file is part of ip_diffim. -# -# Developed for the LSST Data Management System. -# This product includes software developed by the LSST Project -# (https://www.lsst.org). -# See the COPYRIGHT file at the top-level directory of this distribution -# for details of code ownership. -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -import unittest - -import lsst.afw.detection as afwDetection -import lsst.afw.geom as afwGeom -import lsst.afw.table as afwTable -import lsst.geom as geom -import lsst.utils.tests -from lsst.ip.diffim.peakClustering import cluster_peaks, find_peak_clusters - - -def makePeakSchema(): - """Build a source schema and the matching merged-peak schema. - - The peak schema carries the ``merge_peak_positive``/``merge_peak_negative`` - flags, mimicking the footprints produced by ``FootprintMergeList`` in - ``DetectAndMeasureTask``. - """ - schema = afwTable.SourceTable.makeMinimalSchema() - mergeList = afwDetection.FootprintMergeList(schema, ["positive", "negative"]) - return schema, mergeList.getPeakSchema() - - -def makeFootprint(peakSchema, peaks, size=100): - """Make a footprint spanning a ``size`` x ``size`` box with the given peaks. - - Parameters - ---------- - peakSchema : `lsst.afw.table.Schema` - Schema carrying the merge-peak flags. - peaks : `list` [`tuple`] - Each entry is ``(fx, fy, value, isNegative)``. - """ - spans = afwGeom.SpanSet(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(size, size))) - footprint = afwDetection.Footprint(spans, peakSchema) - negativeKey = peakSchema.find("merge_peak_negative").key - positiveKey = peakSchema.find("merge_peak_positive").key - for fx, fy, value, isNegative in peaks: - peak = footprint.addPeak(fx, fy, value) - peak.set(negativeKey, isNegative) - peak.set(positiveKey, not isNegative) - return footprint - - -def clusterContaining(clusters, peakIndex): - """Return the single cluster whose members include ``peakIndex``.""" - matches = [c for c in clusters if peakIndex in c.peak_indices] - assert len(matches) == 1, f"peak {peakIndex} found in {len(matches)} clusters" - return matches[0] - - -class PeakClusteringTestCase(lsst.utils.tests.TestCase): - - def setUp(self): - self.schema, self.peakSchema = makePeakSchema() - self.radius = 5.0 - - def test_two_peaks_far_apart(self): - """Peaks beyond the linking radius form separate singleton clusters.""" - footprint = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (30, 30, 120, False)]) - clusters = cluster_peaks(footprint.getPeaks(), self.radius) - self.assertEqual(len(clusters), 2) - for cluster in clusters: - self.assertEqual(len(cluster.peak_indices), 1) - - def test_two_peaks_close(self): - """Peaks within the linking radius merge into one cluster.""" - footprint = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (12, 10, 90, False)]) - clusters = cluster_peaks(footprint.getPeaks(), self.radius) - self.assertEqual(len(clusters), 1) - self.assertEqual(sorted(clusters[0].peak_indices), [0, 1]) - - def test_transitive_chain(self): - """Single-linkage chains A-B-C even when A and C exceed the radius.""" - # 10->14->18: adjacent gaps are 4 (< 5), but 10->18 is 8 (> 5). - footprint = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (14, 10, 100, False), - (18, 10, 100, False)]) - clusters = cluster_peaks(footprint.getPeaks(), self.radius) - self.assertEqual(len(clusters), 1) - self.assertEqual(sorted(clusters[0].peak_indices), [0, 1, 2]) - - def test_dipole_single_cluster(self): - """A close positive/negative pair clusters together.""" - footprint = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (12, 10, -90, True)]) - clusters = cluster_peaks(footprint.getPeaks(), self.radius) - self.assertEqual(len(clusters), 1) - cluster = clusters[0] - self.assertEqual(cluster.is_negative.tolist(), [False, True]) - self.assertEqual(cluster.peak_values.tolist(), [100, -90]) - - def test_lone_peak_beside_dipole(self): - """A lone peak is separated from a nearby dipole into two clusters.""" - footprint = makeFootprint(self.peakSchema, - [(50, 50, 200, False), # lone peak - (10, 10, 100, False), (12, 10, -90, True)]) # dipole - clusters = cluster_peaks(footprint.getPeaks(), self.radius) - self.assertEqual(len(clusters), 2) - lone = clusterContaining(clusters, 0) - dipole = clusterContaining(clusters, 1) - self.assertEqual(lone.peak_indices, [0]) - self.assertEqual(sorted(dipole.peak_indices), [1, 2]) - - def test_lone_negative_is_own_cluster(self): - """A lone negative peak forms its own singleton cluster.""" - footprint = makeFootprint(self.peakSchema, - [(50, 50, 200, False), (10, 10, -80, True)]) - clusters = cluster_peaks(footprint.getPeaks(), self.radius) - negative = clusterContaining(clusters, 1) - self.assertEqual(negative.peak_indices, [1]) - self.assertEqual(negative.is_negative.tolist(), [True]) - - def test_positions_and_values(self): - """Cluster arrays carry the peaks' positions and signed values.""" - footprint = makeFootprint(self.peakSchema, [(10, 20, 100, False)]) - (cluster,) = cluster_peaks(footprint.getPeaks(), self.radius) - self.assertEqual(cluster.positions.tolist(), [[10, 20]]) - self.assertEqual(cluster.peak_values.tolist(), [100]) - self.assertEqual(cluster.is_negative.tolist(), [False]) - - -class FindPeakClustersTestCase(lsst.utils.tests.TestCase): - - def setUp(self): - self.schema, self.peakSchema = makePeakSchema() - self.radius = 5.0 - - def makeCatalog(self, footprints): - catalog = afwTable.SourceCatalog(self.schema) - for i, footprint in enumerate(footprints): - record = catalog.addNew() - record.setId(i + 1) - record.setFootprint(footprint) - return catalog - - def test_skips_single_peak_footprints(self): - """Footprints with fewer than min_peaks peaks are not returned.""" - single = makeFootprint(self.peakSchema, [(10, 10, 100, False)]) - multi = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (30, 30, 120, False)]) - catalog = self.makeCatalog([single, multi]) - results = find_peak_clusters(catalog, self.radius) - self.assertEqual(len(results), 1) - # Only the second (id=2) footprint qualifies. - self.assertEqual(results[0].source_id, 2) - self.assertEqual(results[0].n_peaks, 2) - - def test_handles_missing_footprint(self): - """Sources with no footprint are skipped rather than raising.""" - multi = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (30, 30, 120, False)]) - catalog = afwTable.SourceCatalog(self.schema) - catalog.addNew().setId(1) # no footprint attached - record = catalog.addNew() - record.setId(2) - record.setFootprint(multi) - results = find_peak_clusters(catalog, self.radius) - self.assertEqual(len(results), 1) - self.assertEqual(results[0].source_id, 2) - - def test_deblendable_classification(self): - """is_deblendable is true when the footprint splits into >1 cluster.""" - # A lone peak plus a dipole: two clusters, deblendable. - twoClusters = makeFootprint(self.peakSchema, - [(50, 50, 200, False), - (10, 10, 100, False), (12, 10, -90, True)]) - # One tight blend: a single cluster, not deblendable. - oneCluster = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (12, 10, 90, False)]) - catalog = self.makeCatalog([twoClusters, oneCluster]) - results = find_peak_clusters(catalog, self.radius) - byId = {r.source_id: r for r in results} - - self.assertTrue(byId[1].is_deblendable) - self.assertEqual(len(byId[1].clusters), 2) - - self.assertFalse(byId[2].is_deblendable) - self.assertEqual(len(byId[2].clusters), 1) - - def test_source_handle_is_live(self): - """The stored source is the live catalog record, so deblender edits - made through it land on the catalog. - """ - multi = makeFootprint(self.peakSchema, - [(10, 10, 100, False), (30, 30, 120, False)]) - catalog = self.makeCatalog([multi]) - results = find_peak_clusters(catalog, self.radius) - # Mutating through the stored handle must be visible in the catalog. - results[0].source.setId(4242) - self.assertEqual(catalog[0].getId(), 4242) - self.assertEqual(results[0].footprint.getBBox(), catalog[0].getFootprint().getBBox()) - - -class MemoryTester(lsst.utils.tests.MemoryTestCase): - pass - - -def setup_module(module): - lsst.utils.tests.init() - - -if __name__ == "__main__": - lsst.utils.tests.init() - unittest.main() From b1f7be26d4c602865ddb534c00aa2cffe72d5c10 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sat, 18 Jul 2026 00:23:17 +0200 Subject: [PATCH 07/12] Set NOT_DEBLENDED flag when clusters are too large --- .../ip/diffim/deblendDifferenceFootprints.py | 16 +++++----- tests/test_deblendDifferenceFootprints.py | 30 ++++++++++++++----- 2 files changed, 32 insertions(+), 14 deletions(-) diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py index 6e0178a8..4b8d33ac 100644 --- a/python/lsst/ip/diffim/deblendDifferenceFootprints.py +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -246,15 +246,12 @@ class DeblendDifferenceFootprintsConfig(pexConfig.Config): maxFootprintArea = pexConfig.Field( dtype=int, default=50000, - doc="Footprints whose bounding-box area (pixels) exceeds this are left " - "un-deblended, bounding the per-footprint cost so a single large " - "footprint cannot dominate the processing time." + doc="Skip deblending footprints whose area exceeds this threshold." ) maxPeaks = pexConfig.Field( dtype=int, default=40, - doc="Footprints with more than this many peaks are left un-deblended, " - "bounding the per-footprint deblending cost." + doc="Footprints with more than this many peaks are left un-deblended." ) @@ -282,8 +279,10 @@ def run(self, diaSources, difference): flags, and footprint peak catalogs the ``merge_peak_negative`` and ``significance`` fields. difference : `lsst.afw.image.ExposureF` - Difference image the footprints were detected on. Only its PSF is - used, to set the clustering and influence scale. + Difference image the footprints were detected on. Its PSF sets the + clustering and influence scale, and its mask is updated in place: + the ``NOT_DEBLENDED`` plane is set for footprints that can't be + deblended. Returns ------- @@ -297,6 +296,8 @@ def run(self, diaSources, difference): linkingRadius = self.config.clusterRadius*sigma*_SIGMA_TO_FWHM isDeblendedKey = diaSources.schema.find("is_deblended").key isNegativeKey = diaSources.schema.find("is_negative").key + difference.mask.addMaskPlane("NOT_DEBLENDED") + notDeblendedBitmask = difference.mask.getPlaneBitMask("NOT_DEBLENDED") nFootprints = 0 nChildren = 0 @@ -314,6 +315,7 @@ def run(self, diaSources, difference): if (footprint.getBBox().getArea() > self.config.maxFootprintArea or len(footprint.getPeaks()) > self.config.maxPeaks): nSkipped += 1 + footprint.getSpans().setMask(difference.mask, notDeblendedBitmask) continue result = _deblend_footprint(footprint, clusters, sigma) if result is None: diff --git a/tests/test_deblendDifferenceFootprints.py b/tests/test_deblendDifferenceFootprints.py index d84d34e3..28b1f883 100644 --- a/tests/test_deblendDifferenceFootprints.py +++ b/tests/test_deblendDifferenceFootprints.py @@ -211,7 +211,8 @@ def setUp(self): self.schema, self.peakSchema = makeSchemas() def makeExposure(self, psfSigma=2.0): - # Only the PSF is used by the task, to set the clustering scale. + # The PSF sets the clustering scale; the mask is updated in place with + # the NOT_DEBLENDED plane for any guardrail-skipped footprints. bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64)) exposure = afwImage.ExposureF(bbox) exposure.setPsf(afwDetection.GaussianPsf(11, 11, psfSigma)) @@ -236,14 +237,18 @@ def test_run_splits_multicluster_footprint(self): footprint = makeFootprint(self.peakSchema, [(10, 20, False, 100.0), (50, 20, False, 100.0)]) catalog = self.makeCatalog([footprint]) task = self.makeTask() + difference = self.makeExposure() - result = task.run(catalog, self.makeExposure()) + result = task.run(catalog, difference) self.assertEqual(len(result.diaSources), 2) self.assertTrue(all(record["is_deblended"] for record in result.diaSources)) self.assertEqual(task.metadata["nDeblendedFootprints"], 1) self.assertEqual(task.metadata["nDeblendedDiaSources"], 1) self.assertEqual(task.metadata["nSkipped"], 0) + # Nothing was skipped, so no pixels are flagged NOT_DEBLENDED. + notDeblended = difference.mask.getPlaneBitMask("NOT_DEBLENDED") + self.assertTrue(np.all((difference.mask.array & notDeblended) == 0)) def test_run_leaves_single_cluster_untouched(self): """A footprint that is one cluster is not modified or flagged.""" @@ -270,29 +275,40 @@ def test_run_sets_is_negative(self): self.assertEqual(byNegative, [False, True]) def test_run_skips_too_many_peaks(self): - """A deblendable footprint with more peaks than maxPeaks is skipped.""" + """A deblendable footprint with more peaks than maxPeaks is skipped and + its pixels are flagged NOT_DEBLENDED. + """ footprint = makeFootprint( self.peakSchema, [(10, 20, False, 100.0), (30, 20, False, 100.0), (50, 20, False, 100.0)]) catalog = self.makeCatalog([footprint]) task = self.makeTask(maxPeaks=2) + difference = self.makeExposure() - result = task.run(catalog, self.makeExposure()) + result = task.run(catalog, difference) self.assertEqual(len(result.diaSources), 1) self.assertFalse(result.diaSources[0]["is_deblended"]) self.assertEqual(task.metadata["nSkipped"], 1) + # The whole (full-frame) footprint is flagged NOT_DEBLENDED. + notDeblended = difference.mask.getPlaneBitMask("NOT_DEBLENDED") + self.assertTrue(np.all((difference.mask.array & notDeblended) != 0)) def test_run_skips_too_large_area(self): - """A deblendable footprint bigger than maxFootprintArea is skipped.""" + """A deblendable footprint bigger than maxFootprintArea is skipped and + its pixels are flagged NOT_DEBLENDED. + """ footprint = makeFootprint(self.peakSchema, [(10, 20, False, 100.0), (50, 20, False, 100.0)]) catalog = self.makeCatalog([footprint]) - task = self.makeTask(maxFootprintArea=100) # footprint area is SIZE**2 + task = self.makeTask(maxFootprintArea=100) # footprint area is 64**2 + difference = self.makeExposure() - result = task.run(catalog, self.makeExposure()) + result = task.run(catalog, difference) self.assertEqual(len(result.diaSources), 1) self.assertFalse(result.diaSources[0]["is_deblended"]) self.assertEqual(task.metadata["nSkipped"], 1) + notDeblended = difference.mask.getPlaneBitMask("NOT_DEBLENDED") + self.assertTrue(np.all((difference.mask.array & notDeblended) != 0)) class MemoryTester(lsst.utils.tests.MemoryTestCase): From f959b01f65a315bceb785b00a3a633fb24f0d6d6 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Wed, 22 Jul 2026 19:59:04 -0700 Subject: [PATCH 08/12] Test using overlap for nearby clusters --- .../ip/diffim/deblendDifferenceFootprints.py | 79 +++++++++++++------ tests/test_deblendDifferenceFootprints.py | 31 +++++++- tests/test_detectAndMeasure.py | 3 + 3 files changed, 83 insertions(+), 30 deletions(-) diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py index 4b8d33ac..28f0678d 100644 --- a/python/lsst/ip/diffim/deblendDifferenceFootprints.py +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -93,21 +93,26 @@ def _cluster_peak_indices(peaks, linking_radius): return [np.nonzero(labels == label)[0].tolist() for label in np.unique(labels)] -def _merge_dominated_clusters(clusters, weighted, n_peaks): - """Merge clusters until every peak owns (dominates) its own peak pixel. +def _merge_clusters(clusters, contamination, n_peaks, max_contamination): + """Merge clusters whose mutual PSF-flux contamination exceeds the limit. - A cluster is absorbed into another when the other has a strictly greater - summed influence at one of the first cluster's own peak pixels. + Two clusters are merged when, at one cluster's peak, the summed + contamination from the other cluster exceeds ``max_contamination`` (a + fraction of that peak's own PSF flux). Merging repeats until every + surviving cluster's peaks are contaminated by less than the limit. Parameters ---------- clusters : `list` [`list` [`int`]] Initial clusters, each a list of peak indices. - weighted : `numpy.ndarray` - ``(n_peaks, n_peaks)`` array where ``weighted[p, j]`` is peak ``p``'s - influence evaluated at peak ``j``'s pixel. + contamination : `numpy.ndarray` + ``(n_peaks, n_peaks)`` array where ``contamination[p, j]`` is the + fractional PSF-flux contamination peak ``p`` imposes on peak ``j``. n_peaks : `int` Number of peaks. + max_contamination : `float` + Maximum fractional contamination a peak may receive from another + cluster before the two clusters are merged. Returns ------- @@ -122,13 +127,13 @@ def _merge_dominated_clusters(clusters, weighted, n_peaks): owner[peak] = index merged = False for j in range(n_peaks): - # Check whether the owning cluster has the most significant - # influence at peak j's pixel; if not, merge it into the one that - # does. - influence = np.array([weighted[cluster, j].sum() for cluster in clusters]) - best = int(np.argmax(influence)) - if influence[best] > influence[owner[j]]: - low, high = sorted((owner[j], best)) + # Total contamination of peak j from each cluster; a peak is not + # contaminated by its own cluster. + received = np.array([contamination[cluster, j].sum() for cluster in clusters]) + received[owner[j]] = 0.0 + worst = int(np.argmax(received)) + if received[worst] > max_contamination: + low, high = sorted((owner[j], worst)) clusters[low] = clusters[low] + clusters[high] del clusters[high] merged = True @@ -157,12 +162,14 @@ def _footprint_from_bool(pixels, bbox, peak_schema, region): return footprint -def _deblend_footprint(footprint, clusters, sigma): +def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05): """Extract every cluster of peaks from a footprint into new footprints. - Each peak's influence is calculated as a gaussian-approximation of the PSF - with amplitude equal to the absolute value of the signal to noise of the - peak on the detection image (```peak['significance']```). + Nearby clusters are first merged when the PSF-flux contamination one would + impose on the other exceeds ``max_contamination``; each remaining cluster + is then extracted into its own footprint. Peak amplitudes are the + detection significance (``peak['significance']``), used as a proxy for + relative flux. Parameters ---------- @@ -174,6 +181,9 @@ def _deblend_footprint(footprint, clusters, sigma): sigma : `float` Width (in pixels) of the 2D Gaussian used to weight each peak's influence on the footprint pixels. + max_contamination : `float`, optional + Maximum fractional PSF-flux contamination a source may receive from a + neighboring cluster before the two are merged. Returns ------- @@ -189,13 +199,18 @@ def _deblend_footprint(footprint, clusters, sigma): positions = np.array([(peak.getIx(), peak.getIy()) for peak in peaks], dtype=int) amplitudes = np.array([peak["significance"] for peak in peaks], dtype=float) - # weighted[p, j]: influence of peak p evaluated at peak j's pixel. + # inv_two_sigma_sq weights the PSF *value* at a pixel (for pixel + # assignment); inv_four_sigma_sq is the PSF-flux *overlap* of two peaks + # (for the contamination merge). inv_two_sigma_sq = 1.0/(2.0*sigma*sigma) + inv_four_sigma_sq = 1.0/(4.0*sigma*sigma) delta = positions[:, np.newaxis, :] - positions[np.newaxis, :, :] dist_sq = np.sum(delta*delta, axis=2) - weighted = amplitudes[:, np.newaxis]*np.exp(-dist_sq*inv_two_sigma_sq) + # contamination[p, j]: fractional PSF-flux bias peak p imposes on peak j -- + # the Gaussian-PSF flux overlap scaled by the significance (flux) ratio. + contamination = (amplitudes[:, np.newaxis]/amplitudes[np.newaxis, :])*np.exp(-dist_sq*inv_four_sigma_sq) - clusters = _merge_dominated_clusters(clusters, weighted, n_peaks) + clusters = _merge_clusters(clusters, contamination, n_peaks, max_contamination) # Extract every cluster into its own source; skip footprints that did not # split into more than one cluster. @@ -239,9 +254,11 @@ class DeblendDifferenceFootprintsConfig(pexConfig.Config): clusterRadius = pexConfig.Field( dtype=float, - default=3.0, - doc="Radius (in PSF FWHM) within which detected peaks are grouped into " - "one cluster. Each cluster is extracted into its own diaSource." + default=2., + doc="Initial clustering radius (in PSF FWHM) to group tightly-packed " + "peaks. Clusters may ultimately merge with significantly " + "contaminating neighbors. Each final cluster is extracted into its " + "own diaSource." ) maxFootprintArea = pexConfig.Field( dtype=int, @@ -250,9 +267,19 @@ class DeblendDifferenceFootprintsConfig(pexConfig.Config): ) maxPeaks = pexConfig.Field( dtype=int, - default=40, + default=100, doc="Footprints with more than this many peaks are left un-deblended." ) + maxContamination = pexConfig.Field( + dtype=float, + default=0.05, + doc="Maximum fractional PSF-flux contamination a diaSource may receive " + "from a neighboring cluster before the two are merged into one " + "(un-deblended) source. For example, 0.05 keeps each deblended " + "source's PSF flux biased by less than 5% by its neighbors. " + "Contamination is the Gaussian-PSF flux overlap scaled by the " + "neighbor-to-source significance ratio." + ) class DeblendDifferenceFootprintsTask(pipeBase.Task): @@ -317,7 +344,7 @@ def run(self, diaSources, difference): nSkipped += 1 footprint.getSpans().setMask(difference.mask, notDeblendedBitmask) continue - result = _deblend_footprint(footprint, clusters, sigma) + result = _deblend_footprint(footprint, clusters, sigma, self.config.maxContamination) if result is None: continue self._extractRecords(diaSources, source, result, isDeblendedKey, isNegativeKey) diff --git a/tests/test_deblendDifferenceFootprints.py b/tests/test_deblendDifferenceFootprints.py index 28b1f883..9a10c3fc 100644 --- a/tests/test_deblendDifferenceFootprints.py +++ b/tests/test_deblendDifferenceFootprints.py @@ -131,11 +131,12 @@ class DeblendFootprintTestCase(lsst.utils.tests.TestCase): def setUp(self): self.schema, self.peakSchema = makeSchemas() - def _deblend(self, peaks, linkingRadius=5.0, sigma=4.0): + def _deblend(self, peaks, linkingRadius=5.0, sigma=4.0, max_contamination=0.05): footprint = makeFootprint(self.peakSchema, peaks) clusters = _cluster_peak_indices(footprint.getPeaks(), linkingRadius) - result = _deblend_footprint(footprint, clusters, sigma=sigma) + result = _deblend_footprint(footprint, clusters, sigma=sigma, + max_contamination=max_contamination) return footprint, clusters, result def test_extract_clusters_from_footprint(self): @@ -194,8 +195,8 @@ def test_extract_negative_cluster(self): self.assertEqual(extracted, {(8, 20), (32, 20)}) def test_low_significance_peak_is_merged(self): - """A faint (low-significance) peak dominated by a bright neighbour is - merged back in, leaving nothing to extract. + """A faint (low-significance) peak heavily contaminated by a bright + neighbour is merged back in, leaving nothing to extract. """ # Faint positive at (12, 20); bright dipole at (16, 20)/(18, 20). _, _, result = self._deblend( @@ -203,6 +204,28 @@ def test_low_significance_peak_is_merged(self): linkingRadius=3.0, sigma=2.0) self.assertIsNone(result) + def test_contamination_merges_beyond_cluster_radius(self): + """Peaks the KD step leaves separate are still merged when their PSF + flux contamination exceeds ``max_contamination``. + """ + # d=5 px, sigma=2 -> overlap exp(-25/16)=0.21 > 0.05, but the peaks are + # beyond linkingRadius=3 px, so the contamination test drives the merge. + _, clusters, result = self._deblend( + [(20, 20, False, 100.0), (25, 20, False, 100.0)], + linkingRadius=3.0, sigma=2.0, max_contamination=0.05) + self.assertEqual(len(clusters), 2) # KD kept them separate + self.assertIsNone(result) # contamination merged them + + def test_no_merge_when_contamination_below_limit(self): + """Peaks whose mutual contamination is below the limit stay separate.""" + # d=10 px, sigma=2 -> overlap exp(-100/16)=0.0019 < 0.05. + _, clusters, result = self._deblend( + [(20, 20, False, 100.0), (30, 20, False, 100.0)], + linkingRadius=3.0, sigma=2.0, max_contamination=0.05) + self.assertEqual(len(clusters), 2) + self.assertIsNotNone(result) + self.assertEqual(len(result.children), 1) + class RunTestCase(lsst.utils.tests.TestCase): """Tests of the task-level ``run`` over a diaSource catalog.""" diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index 04c074e8..a3c9169d 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -1460,6 +1460,9 @@ def testMergeFootprints(self): config = detectAndMeasure.DetectAndMeasureTask.ConfigClass() config.doDeblend = True + # Keep the positive blended pair (~10 px apart) in a single cluster, so + # this test exercises footprint merging rather than deblend splitting. + config.deblend.clusterRadius = 3.0 config.raiseOnBadSubtractionRatio = False config.raiseOnNoDiaSources = False task = detectAndMeasure.DetectAndMeasureTask(config=config) From 17401e6f9dbf9e13a542234959c47c49e30e244b Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 24 Jul 2026 11:23:13 -0700 Subject: [PATCH 09/12] Use the footprint noise replacer to make deblended diaSources less biased --- python/lsst/ip/diffim/detectAndMeasure.py | 6 ++ tests/test_detectAndMeasure.py | 125 +++++++++++++++++++++- 2 files changed, 130 insertions(+), 1 deletion(-) diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py index a81ff790..0313f6f2 100644 --- a/python/lsst/ip/diffim/detectAndMeasure.py +++ b/python/lsst/ip/diffim/detectAndMeasure.py @@ -496,6 +496,12 @@ def setDefaults(self): self.measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments" self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments" self.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 5.0 + # Replace neighboring sources with noise while measuring each diaSource. + # Important when deblending so that the deblended diaSource flux + # measurements are unbiased. Note that this deblending merely divides + # the diaSource footprints - it does not fit fractional flux from + # overlapping diaSources in contaminated pixels. + self.measurement.doReplaceWithNoise = True self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"] self.forcedMeasurement.copyColumns = { "id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"} diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py index a3c9169d..7f9866c7 100644 --- a/tests/test_detectAndMeasure.py +++ b/tests/test_detectAndMeasure.py @@ -126,7 +126,8 @@ def _check_values(self, values, minValue=None, maxValue=None): def _setup_detection(self, doSkySources=True, nSkySources=5, doSubtractBackground=False, run_sattle=False, - badSourceFlags=None, clusterRadius=None, **kwargs): + badSourceFlags=None, clusterRadius=None, + doReplaceWithNoise=None, **kwargs): """Setup and configure the detection and measurement PipelineTask. Parameters @@ -150,6 +151,8 @@ def _setup_detection(self, doSkySources=True, nSkySources=5, config.update(**kwargs) if clusterRadius is not None: config.deblend.clusterRadius = clusterRadius + if doReplaceWithNoise is not None: + config.measurement.doReplaceWithNoise = doReplaceWithNoise config.run_sattle = run_sattle @@ -214,6 +217,64 @@ def test_detection_xy0(self): for source in sources: self._check_diaSource(output.diaSources, source, refIds=refIds) + def test_deblended_measurement_limited_to_footprint(self): + """Deblended sources should be measured from only their own footprint + pixels: with neighbor pixels replaced by noise, a source's aperture + flux recovers its injected value, but should be biased without + replacement. + """ + # Set up the simulated images. + noiseLevel = 1. + staticSeed = 1 + transientSeed = 6 + psfSize = 2.4 + apColumn = "base_CircularApertureFlux_12_0_instFlux" + kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0, "xSize": 256, "ySize": 256} + science, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6, nSrc=1, **kwargs) + science.getInfo().setVisitInfo(makeVisitInfo()) + matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7, + nSrc=1, **kwargs) + + # Two equal-flux transients, close enough to share a footprint (so they + # are deblended) but far enough apart to be detected as separate peaks + # and split into two clusters. Their separation is well within the 12 + # pixel aperture, so each source's aperture overlaps the other. + fluxLevel = 3000. + separation = 10. + center = 128. + xLoc = [center - separation/2, center + separation/2] + yLoc = [center, center] + transients, _ = makeTestImage(seed=transientSeed, nSrc=2, xLoc=xLoc, yLoc=yLoc, + flux=[fluxLevel, fluxLevel], + noiseLevel=noiseLevel, noiseSeed=8, **kwargs) + difference = science.clone() + difference.maskedImage -= matchedTemplate.maskedImage + difference.maskedImage += transients.maskedImage + + def _run(doReplaceWithNoise): + detectionTask = self._setup_detection(doDeblend=True, doSkySources=False, + clusterRadius=1.0, badSubtractionRatioThreshold=1., + doReplaceWithNoise=doReplaceWithNoise) + diaSources = detectionTask.run( + science, matchedTemplate, difference.clone(), + idFactory=self.idGenerator.make_table_id_factory()).diaSources + # The blend should split into two separate deblended diaSources. + self.assertEqual(len(diaSources), 2) + self.assertTrue(all(diaSources["is_deblended"])) + return np.array(diaSources[apColumn]) + + apFluxReplaced = _run(doReplaceWithNoise=True) + apFluxContaminated = _run(doReplaceWithNoise=False) + + # Every source's aperture flux is inflated by its neighbor when the + # neighbor pixels are left in place. + self.assertTrue(np.all(apFluxContaminated > apFluxReplaced)) + # With the neighbor replaced by noise, aperture flux recovers the + # injected flux; without it, the aperture flux is inflated well beyond + # tolerance. + self.assertFloatsAlmostEqual(np.mean(apFluxReplaced), fluxLevel, rtol=0.1) + self.assertGreater(np.mean(apFluxContaminated), fluxLevel*1.3) + def test_raise_bad_psf(self): """Detection should raise if the PSF width is NaN """ @@ -981,6 +1042,68 @@ def test_detection_xy0(self): for source in sources: self._check_diaSource(output.diaSources, source, refIds=refIds) + def test_deblended_measurement_limited_to_footprint(self): + """Deblended sources should be measured from only their own footprint + pixels: with neighbor pixels replaced by noise, a source's aperture + flux recovers its injected value, but should be biased without + replacement. + """ + # Set up the simulated images. + noiseLevel = 1. + staticSeed = 1 + transientSeed = 6 + psfSize = 2.4 + apColumn = "base_CircularApertureFlux_12_0_instFlux" + kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0, "xSize": 256, "ySize": 256} + science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6, + nSrc=1, **kwargs) + science.getInfo().setVisitInfo(makeVisitInfo()) + matchedTemplate, _ = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel/4, noiseSeed=7, + nSrc=1, **kwargs) + + # Two equal-flux transients, close enough to share a footprint (so they + # are deblended) but far enough apart to be detected as separate peaks + # and split into two clusters. Their separation is well within the 12 + # pixel aperture, so each source's aperture overlaps the other. + fluxLevel = 3000. + separation = 10. + center = 128. + xLoc = [center - separation/2, center + separation/2] + yLoc = [center, center] + transients, _ = makeTestImage(seed=transientSeed, nSrc=2, xLoc=xLoc, yLoc=yLoc, + flux=[fluxLevel, fluxLevel], + noiseLevel=noiseLevel, noiseSeed=8, **kwargs) + difference = science.clone() + difference.maskedImage -= matchedTemplate.maskedImage + difference.maskedImage += transients.maskedImage + subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask() + scienceKernel = science.psf.getKernel() + score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl) + + def _run(doReplaceWithNoise): + detectionTask = self._setup_detection(doDeblend=True, doSkySources=False, + clusterRadius=1.0, badSubtractionRatioThreshold=1., + doReplaceWithNoise=doReplaceWithNoise) + diaSources = detectionTask.run( + science, matchedTemplate, difference.clone(), score.clone(), sources, + idFactory=self.idGenerator.make_table_id_factory()).diaSources + # The blend should split into two separate deblended diaSources. + self.assertEqual(len(diaSources), 2) + self.assertTrue(all(diaSources["is_deblended"])) + return np.array(diaSources[apColumn]) + + apFluxReplaced = _run(doReplaceWithNoise=True) + apFluxContaminated = _run(doReplaceWithNoise=False) + + # Every source's aperture flux is inflated by its neighbor when the + # neighbor pixels are left in place. + self.assertTrue(np.all(apFluxContaminated > apFluxReplaced)) + # With the neighbor replaced by noise, aperture flux recovers the + # injected flux; without it, the aperture flux is inflated well beyond + # tolerance. + self.assertFloatsAlmostEqual(np.mean(apFluxReplaced), fluxLevel, rtol=0.1) + self.assertGreater(np.mean(apFluxContaminated), fluxLevel*1.3) + def test_measurements_finite(self): """Measured fluxes and centroids should always be finite. """ From ec375af6c08b832ccafe4cc5aba71bf464124c5d Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Sun, 26 Jul 2026 12:48:00 -0700 Subject: [PATCH 10/12] Refine the deblending limits on footprint size. Use a smaller check on the actual size of the footprint spans, and a larger limit on the enclosing bbox. --- .../ip/diffim/deblendDifferenceFootprints.py | 13 ++- tests/test_deblendDifferenceFootprints.py | 82 +++++++++++++++++++ 2 files changed, 92 insertions(+), 3 deletions(-) diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py index 28f0678d..9bd01a5d 100644 --- a/python/lsst/ip/diffim/deblendDifferenceFootprints.py +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -262,8 +262,14 @@ class DeblendDifferenceFootprintsConfig(pexConfig.Config): ) maxFootprintArea = pexConfig.Field( dtype=int, - default=50000, - doc="Skip deblending footprints whose area exceeds this threshold." + default=10000, + doc="Skip deblending footprints larger than this, in pixels." + ) + maxFootprintBBoxArea = pexConfig.Field( + dtype=int, + default=100000, + doc="Additional check to skip deblending footprints with an enclosing " + "bounding box larger than this, in pixels." ) maxPeaks = pexConfig.Field( dtype=int, @@ -339,7 +345,8 @@ def run(self, diaSources, difference): if len(clusters) < 2: continue # Leave very large footprints un-deblended - if (footprint.getBBox().getArea() > self.config.maxFootprintArea + if (footprint.getArea() > self.config.maxFootprintArea + or footprint.getBBox().getArea() > self.config.maxFootprintBBoxArea or len(footprint.getPeaks()) > self.config.maxPeaks): nSkipped += 1 footprint.getSpans().setMask(difference.mask, notDeblendedBitmask) diff --git a/tests/test_deblendDifferenceFootprints.py b/tests/test_deblendDifferenceFootprints.py index 9a10c3fc..b3b582ff 100644 --- a/tests/test_deblendDifferenceFootprints.py +++ b/tests/test_deblendDifferenceFootprints.py @@ -71,6 +71,25 @@ def makeFootprint(peakSchema, peaks): return footprint +def makeDiagonalFootprint(peakSchema, peaks, width=3): + """Build a thin diagonal footprint carrying the given peaks. + """ + mask = afwImage.Mask(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64))) + bit = mask.getPlaneBitMask("DETECTED") + rows, columns = np.mgrid[0:64, 0:64] + mask.array[np.abs(rows - columns) <= width] = bit + footprint = afwDetection.Footprint(afwGeom.SpanSet.fromMask(mask, bit), peakSchema) + negativeKey = peakSchema.find("merge_peak_negative").key + positiveKey = peakSchema.find("merge_peak_positive").key + significanceKey = peakSchema.find("significance").key + for ix, iy, isNegative, significance in peaks: + peak = footprint.addPeak(ix, iy, -significance if isNegative else significance) + peak.set(negativeKey, isNegative) + peak.set(positiveKey, not isNegative) + peak.set(significanceKey, significance) + return footprint + + class ClusterPeaksTestCase(lsst.utils.tests.TestCase): """Tests of the initial single-linkage peak clustering.""" @@ -334,6 +353,69 @@ def test_run_skips_too_large_area(self): self.assertTrue(np.all((difference.mask.array & notDeblended) != 0)) +class FootprintAreaTestCase(lsst.utils.tests.TestCase): + """Tests for the limits on footprint size.""" + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + # Well beyond the KD tree linking radius, and both on the diagonal. + self.peaks = [(10, 10, False, 100.0), (50, 50, False, 100.0)] + + def makeExposure(self): + bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64)) + exposure = afwImage.ExposureF(bbox) + exposure.setPsf(afwDetection.GaussianPsf(21, 21, 2.0)) + exposure.variance.array[:] = 1.0 + return exposure + + def runTask(self, footprint, **configKwargs): + catalog = afwTable.SourceCatalog(self.schema) + catalog.addNew().setFootprint(footprint) + config = DeblendDifferenceFootprintsTask.ConfigClass() + config.clusterRadius = 2.0 + for key, value in configKwargs.items(): + setattr(config, key, value) + task = DeblendDifferenceFootprintsTask(config=config) + return task, task.run(catalog, self.makeExposure()).diaSources + + def test_diagonal_footprint_is_deblended(self): + """A narrow but long diagonal footprint should be deblended.""" + # The footprint should have ~400 pixels + footprint = makeDiagonalFootprint(self.peakSchema, self.peaks) + self.assertLess(footprint.getArea(), 1000) + self.assertGreater(footprint.getBBox().getArea(), 1000) + + task, diaSources = self.runTask(footprint, maxFootprintArea=1000) + + self.assertEqual(len(diaSources), 2) + self.assertEqual(task.metadata["nSkipped"], 0) + + def test_compact_footprint_of_same_bbox_is_skipped(self): + """An equivalent filled footprint will be skipped.""" + footprint = makeFootprint(self.peakSchema, self.peaks) + self.assertGreater(footprint.getArea(), 1000) + self.assertGreater(footprint.getBBox().getArea(), 1000) + + task, diaSources = self.runTask(footprint, maxFootprintArea=1000) + + self.assertEqual(len(diaSources), 1) + self.assertEqual(task.metadata["nSkipped"], 1) + + def test_bbox_backstop_skips_widely_spread_footprint(self): + """A narrow diagonal footprint will be skipped if the bbox is large + enough. + """ + footprint = makeDiagonalFootprint(self.peakSchema, self.peaks) + self.assertLess(footprint.getArea(), 1000) + self.assertGreater(footprint.getBBox().getArea(), 1000) + + task, diaSources = self.runTask(footprint, maxFootprintArea=1000, + maxFootprintBBoxArea=1000) + + self.assertEqual(len(diaSources), 1) + self.assertEqual(task.metadata["nSkipped"], 1) + + class MemoryTester(lsst.utils.tests.MemoryTestCase): pass From 43a5653c5c90d8f1c0c0fb8c243da2d5c73959ba Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Mon, 27 Jul 2026 10:18:39 -0700 Subject: [PATCH 11/12] Raise the deblending size and peak limits --- python/lsst/ip/diffim/deblendDifferenceFootprints.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py index 9bd01a5d..199158a7 100644 --- a/python/lsst/ip/diffim/deblendDifferenceFootprints.py +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -262,18 +262,18 @@ class DeblendDifferenceFootprintsConfig(pexConfig.Config): ) maxFootprintArea = pexConfig.Field( dtype=int, - default=10000, + default=100000, doc="Skip deblending footprints larger than this, in pixels." ) maxFootprintBBoxArea = pexConfig.Field( dtype=int, - default=100000, + default=600000, doc="Additional check to skip deblending footprints with an enclosing " "bounding box larger than this, in pixels." ) maxPeaks = pexConfig.Field( dtype=int, - default=100, + default=200, doc="Footprints with more than this many peaks are left un-deblended." ) maxContamination = pexConfig.Field( @@ -349,6 +349,9 @@ def run(self, diaSources, difference): or footprint.getBBox().getArea() > self.config.maxFootprintBBoxArea or len(footprint.getPeaks()) > self.config.maxPeaks): nSkipped += 1 + self.log.debug("Skipped footprint with area %d, bounding box area %d, and %d peaks.", + footprint.getArea(), footprint.getBBox().getArea(), + len(footprint.getPeaks())) footprint.getSpans().setMask(difference.mask, notDeblendedBitmask) continue result = _deblend_footprint(footprint, clusters, sigma, self.config.maxContamination) From 89aa701a93dd31d6455e9a03032e1a1f5237e1e7 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Tue, 27 Oct 2026 16:48:26 -0700 Subject: [PATCH 12/12] Keep trailed sources together when deblending Use FindGlintTrailsTask from meas_algorithms to identify collinear peaks that should be kept together as a long trailed source. --- .../ip/diffim/deblendDifferenceFootprints.py | 396 +++++++++++++++++- tests/test_deblendDifferenceFootprints.py | 351 +++++++++++++++- 2 files changed, 729 insertions(+), 18 deletions(-) diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py index 199158a7..e9c5f344 100644 --- a/python/lsst/ip/diffim/deblendDifferenceFootprints.py +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -25,8 +25,11 @@ __all__ = ["DeblendDifferenceFootprintsConfig", "DeblendDifferenceFootprintsTask"] import dataclasses +import functools +from astropy.stats import gaussian_sigma_to_fwhm import numpy as np +from scipy import ndimage from scipy.sparse import coo_matrix from scipy.sparse.csgraph import connected_components from scipy.spatial import cKDTree @@ -34,6 +37,9 @@ import lsst.afw.detection import lsst.afw.geom import lsst.afw.image +import lsst.afw.table +import lsst.geom +from lsst.meas.algorithms import FindGlintTrailsTask, fit_line import lsst.pex.config as pexConfig import lsst.pipe.base as pipeBase @@ -41,9 +47,6 @@ # originating from the negative-polarity detection footprint set. _NEGATIVE_PEAK_FLAG = "merge_peak_negative" -# Conversion from a Gaussian sigma to its full width at half maximum. -_SIGMA_TO_FWHM = np.sqrt(8.0*np.log(2.0)) - @dataclasses.dataclass class _DeblendedFootprint: @@ -56,10 +59,14 @@ class _DeblendedFootprint: children : `list` [`lsst.afw.detection.Footprint`] One new footprint per remaining cluster, each carrying that cluster's peaks. + trails : `list` [`bool`] + Whether each footprint of ``[parent] + children`` holds clusters that + were merged as one trailed source. """ parent: lsst.afw.detection.Footprint children: list + trails: list def _cluster_peak_indices(peaks, linking_radius): @@ -85,6 +92,7 @@ def _cluster_peak_indices(peaks, linking_radius): if n <= 1: return [list(range(n))] pairs = cKDTree(positions).query_pairs(r=linking_radius, output_type="ndarray") + # Fill the sparse matrix in "coordinate" (COO) format if len(pairs) == 0: graph = coo_matrix((n, n)) else: @@ -93,12 +101,251 @@ def _cluster_peak_indices(peaks, linking_radius): return [np.nonzero(labels == label)[0].tolist() for label in np.unique(labels)] +def _segment_on_mask(mask, origin, scale, start, end): + """Test whether the straight segment between two points lies entirely on + set pixels of a mask. + + Parameters + ---------- + mask : `numpy.ndarray` [`bool`] + Pixels to test against. + origin : `numpy.ndarray` + ``(x, y)`` pixel position of ``mask[0, 0]``. + scale : `float` + Pixels per unit of ``start`` and ``end``. + start, end : `numpy.ndarray` + ``(x, y)`` positions of the ends of the segment, as pixel coordinates + divided by ``scale``. + + Returns + ------- + on_mask : `bool` + Whether every pixel on an unbroken path of pixels along the segment + is set in ``mask``. + """ + start = start*scale - origin + end = end*scale - origin + # Samples at most half a pixel apart round to the same or adjacent + # (possibly diagonal) pixels, so the pixels tested form an unbroken path. + # A pixel whose corner the segment barely clips may not be tested. + n_samples = int(np.ceil(2.0*np.hypot(*(end - start)))) + 1 + x, y = np.linspace(start, end, n_samples).T + return bool(np.all(mask[np.rint(y).astype(int), np.rint(x).astype(int)])) + + +def _is_continuous(line, trail, candidate, connected): + """Test whether a candidate cluster joins a trail without crossing a gap + in its footprint. + + Parameters + ---------- + line : `lsst.meas.algorithms.FittedLine` + Line fit to the points of ``trail`` and ``candidate`` together. + trail : `numpy.ndarray` + ``(n, 2)`` array of the trail's points, as pixel coordinates divided + by the PSF FWHM. + candidate : `numpy.ndarray` + ``(m, 2)`` array of the candidate cluster's points, in the same units + as ``trail``. + connected : callable + Takes two positions and returns whether the segment between them lies + on the footprint. + + Returns + ------- + continuous : `bool` + Whether each candidate point is joined to its neighbors along the line. + """ + points = np.concatenate((trail, candidate)) + n_old = len(trail) + # Sort the points along the line. Each point is projected onto the line's + # direction, giving its distance along the trail from the line's center. + order = np.argsort((points - line.center) @ line.direction) + return all(connected(points[a], points[b]) for a, b in zip(order[:-1], order[1:]) + if a >= n_old or b >= n_old) + + +def _split_at_gaps(members, points, connected): + """Split a line of single-point clusters wherever the footprint has a gap + between neighbors along the line. + + Parameters + ---------- + members : `list` [`int`] + Clusters on the line, each standing as one point. + points : `list` [`numpy.ndarray`] + ``(1, 2)`` array of each cluster's point, as pixel coordinates divided + by the PSF FWHM. + connected : callable + Takes two positions and returns whether the segment between them lies + on the footprint. + + Returns + ------- + runs : `list` [`list` [`int`]] + The members of each continuous piece of the line, in order along it. + """ + xy = np.concatenate([points[index] for index in members]) + line = fit_line(*xy.T) + order = np.argsort((xy - line.center) @ line.direction) + runs = [[members[order[0]]]] + for a, b in zip(order[:-1], order[1:]): + if not connected(xy[a], xy[b]): + runs.append([]) + runs[-1].append(members[b]) + return runs + + +def _grow_trail(members, points, absorbed, positive, tolerance, connected): + """Find the positive clusters that extend a trail while all its points + stay collinear and joined by the footprint. + + A cluster is added when every point of the trail and the cluster lies + within ``tolerance`` of the line fit to them all, and the footprint has no + gaps between the cluster's points and its neighbors. Clusters are added one + at a time, best fit first, so that clusters on the line fix its direction + before a nearby peak off the line is tested; such a peak then no longer + fits and is left out of the trail. + + Parameters + ---------- + members : `list` [`int`] + Clusters already in the trail. + points : `list` [`numpy.ndarray`] + ``(n, 2)`` array of the points that stand for each cluster on a line, + as pixel coordinates divided by the PSF FWHM. + absorbed : `list` [`bool`] + Whether each cluster already belongs to another. + positive : `list` [`bool`] + Whether all of each cluster's peaks are positive. + tolerance : `float` + Maximum distance, in PSF FWHM, of any point from the trail's line. + connected : callable + Takes two positions and returns whether the segment between them lies + on the footprint. + + Returns + ------- + members : `list` [`int`] + Clusters of the grown trail, starting with the original ``members``. + """ + members = list(members) + trail = np.concatenate([points[index] for index in members]) + while True: + best = None + best_residual = tolerance + for index in range(len(points)): + if index in members or absorbed[index] or not positive[index]: + continue + # Fit a new line to the trail and this cluster. A line fit to only + # a few closely spaced points on a trail can miss a distant cluster + # on the same trail, even if they have a good simultaneous fit. + xy = np.concatenate((trail, points[index])) + line = fit_line(*xy.T) + residual = np.max(line.distance(*xy.T)) + if residual <= best_residual and _is_continuous(line, trail, points[index], connected): + best = index + best_residual = residual + if best is None: + return members + members.append(best) + trail = np.concatenate((trail, points[best])) + + +def _merge_trail_clusters(clusters, positions, negative, find_lines, connected, min_points, tolerance, + compact_span): + """Merge clusters of positive peaks that lie along one line. + + Find lines among the positive clusters of peaks, and grow each into a + trail that absorbs the positive clusters on its line (see `_grow_trail`). + A cluster whose peaks span less than ``compact_span`` is treated as a + single point at its mean position, while a wider cluster starts a line + along its own peaks. A trail is kept only if it holds at least + ``min_points`` clusters. Clusters with a negative peak are never part of a + trail, and a trail never spans a gap in the footprint. + + Parameters + ---------- + clusters : `list` [`list` [`int`]] + Clusters of peaks within one footprint, as rows into its peaks. + positions : `numpy.ndarray` + ``(n_peaks, 2)`` array of the peaks' positions, as pixel coordinates + divided by the PSF FWHM. + negative : `numpy.ndarray` [`bool`] + Whether each peak is negative. + find_lines : callable + Takes an ``(n, 2)`` array of positions and returns a list of lines, + each a list of rows into that array. + connected : callable + Takes two ``(x, y)`` positions, in the same units as ``positions``, + and returns whether the straight segment between them lies on the + footprint. + min_points : `int` + Minimum number of clusters in a trail. + tolerance : `float` + Maximum distance, in PSF FWHM, of any point from a trail's line. + compact_span : `float` + Clusters whose peaks span less than this distance (in PSF FWHM) will + be considered as a single point for line fitting. + + Returns + ------- + clusters : `list` [`list` [`int`]] + The merged clusters. + trails : `list` [`bool`] + Whether each returned cluster was merged as a trailed source. + """ + clusters = [list(cluster) for cluster in clusters] + positive = [not np.any(negative[cluster]) for cluster in clusters] + absorbed = [False]*len(clusters) + trails = [False]*len(clusters) + compact = [np.max(np.hypot(*(positions[cluster][:, np.newaxis] - positions[cluster]).T)) < compact_span + for cluster in clusters] + points = [positions[cluster].mean(axis=0, keepdims=True) if small else positions[cluster] + for cluster, small in zip(clusters, compact)] + + def join(members): + """Merge the clusters of a trail into its first.""" + seed = members[0] + for index in members[1:]: + clusters[seed] = clusters[seed] + clusters[index] + points[seed] = np.concatenate((points[seed], points[index])) + absorbed[index] = True + trails[seed] = True + + # Larger clusters fix their line better, so they seed first. + for index in np.argsort([-len(cluster) for cluster in clusters], kind="stable"): + if not compact[index] and positive[index] and not absorbed[index]: + members = _grow_trail([index], points, absorbed, positive, tolerance, connected) + if len(members) >= min_points: + join(members) + + candidates = [index for index in range(len(clusters)) + if compact[index] and positive[index] and not absorbed[index]] + if len(candidates) >= min_points: + for line in find_lines(np.concatenate([points[index] for index in candidates])): + # Lines may share points, and a trail grown from an earlier line + # may already hold some of this one. + members = [candidates[row] for row in line + if not absorbed[candidates[row]] and not trails[candidates[row]]] + if len(members) < min_points: + continue + for run in _split_at_gaps(members, points, connected): + # A trail grown from an earlier run may hold some of this one. + run = [index for index in run if not absorbed[index]] + if len(run) >= min_points: + join(_grow_trail(run, points, absorbed, positive, tolerance, connected)) + + keep = [index for index in range(len(clusters)) if not absorbed[index]] + return [clusters[index] for index in keep], [trails[index] for index in keep] + + def _merge_clusters(clusters, contamination, n_peaks, max_contamination): """Merge clusters whose mutual PSF-flux contamination exceeds the limit. Two clusters are merged when, at one cluster's peak, the summed - contamination from the other cluster exceeds ``max_contamination`` (a - fraction of that peak's own PSF flux). Merging repeats until every + contamination from the other cluster exceeds a ``max_contamination`` + fraction of that peak's own PSF flux. Merging repeats until every surviving cluster's peaks are contaminated by less than the limit. Parameters @@ -162,14 +409,15 @@ def _footprint_from_bool(pixels, bbox, peak_schema, region): return footprint -def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05): +def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05, trails=None): """Extract every cluster of peaks from a footprint into new footprints. Nearby clusters are first merged when the PSF-flux contamination one would impose on the other exceeds ``max_contamination``; each remaining cluster is then extracted into its own footprint. Peak amplitudes are the detection significance (``peak['significance']``), used as a proxy for - relative flux. + relative flux. A cluster that absorbs a trailed source is itself a + trailed source. Parameters ---------- @@ -184,14 +432,21 @@ def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05): max_contamination : `float`, optional Maximum fractional PSF-flux contamination a source may receive from a neighboring cluster before the two are merged. + trails : `list` [`bool`], optional + Whether each of ``clusters`` is a trailed source; none are if `None`. Returns ------- result : `_DeblendedFootprint` or `None` The parent footprint (kept on the original record) and the new footprints, or `None` if the footprint did not split into more than - one cluster. + one cluster and holds no trailed source. A footprint whose clusters + all merged into one trailed source is returned unchanged as the + parent, with no children. """ + if trails is None: + trails = [False]*len(clusters) + trail_peaks = {peak for cluster, trail in zip(clusters, trails) if trail for peak in cluster} peaks = footprint.getPeaks() n_peaks = len(peaks) peak_schema = peaks.getSchema() @@ -211,10 +466,13 @@ def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05): contamination = (amplitudes[:, np.newaxis]/amplitudes[np.newaxis, :])*np.exp(-dist_sq*inv_four_sigma_sq) clusters = _merge_clusters(clusters, contamination, n_peaks, max_contamination) + trails = [any(peak in trail_peaks for peak in cluster) for cluster in clusters] # Extract every cluster into its own source; skip footprints that did not # split into more than one cluster. if len(clusters) < 2: + if trails[0]: + return _DeblendedFootprint(parent=footprint, children=[], trails=[True]) return None # Assign every footprint pixel to the cluster of highest influence. @@ -229,12 +487,37 @@ def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05): for peak in cluster: r_sq = (xx - positions[peak, 0])**2 + (yy - positions[peak, 1])**2 accumulator += amplitudes[peak]*np.exp(-r_sq*inv_two_sigma_sq) + if trails[index]: + # Spread a trail's influence along the segment between its end + # peaks, so that a neighbor beside a gap between its crests does + # not cut its footprint in two. + rows = [peak for peak in cluster if peak in trail_peaks] + line = fit_line(positions[rows, 0], positions[rows, 1]) + along = (positions[rows] - line.center) @ line.direction + dx = xx - line.center[0] + dy = yy - line.center[1] + t = np.clip(dx*line.direction[0] + dy*line.direction[1], along.min(), along.max()) + d_sq = (dx - t*line.direction[0])**2 + (dy - t*line.direction[1])**2 + accumulator = np.maximum(accumulator, + np.median(amplitudes[rows])*np.exp(-d_sq*inv_two_sigma_sq)) influence[index] = accumulator labels = np.argmax(influence, axis=0) # Guarantee each cluster owns its own peak for index, cluster in enumerate(clusters): for peak in cluster: labels[positions[peak, 1] - min_y, positions[peak, 0] - min_x] = index + # Far from every peak the influences are tiny and can cross over, leaving + # a trail stray pixels cut off from its line. Give any piece of a trail + # holding none of its peaks to the next most influential cluster. + for index, cluster in enumerate(clusters): + if not trails[index]: + continue + pieces, _ = ndimage.label(footprint_mask & (labels == index)) + held = [pieces[positions[peak, 1] - min_y, positions[peak, 0] - min_x] for peak in cluster] + stray = (pieces > 0) & ~np.isin(pieces, held) + if np.any(stray): + influence[index][stray] = -np.inf + labels[stray] = np.argmax(influence[:, stray], axis=0) # One footprint per cluster, carrying that cluster's peaks. region = footprint.getRegion() @@ -246,7 +529,7 @@ def _deblend_footprint(footprint, clusters, sigma, max_contamination=0.05): footprints.append(fp) # The first cluster stays on the original record; the rest are new records. - return _DeblendedFootprint(parent=footprints[0], children=footprints[1:]) + return _DeblendedFootprint(parent=footprints[0], children=footprints[1:], trails=trails) class DeblendDifferenceFootprintsConfig(pexConfig.Config): @@ -286,6 +569,29 @@ class DeblendDifferenceFootprintsConfig(pexConfig.Config): "Contamination is the Gaussian-PSF flux overlap scaled by the " "neighbor-to-source significance ratio." ) + compactClusterSpan = pexConfig.Field( + dtype=float, + default=1.5, + doc="Clusters whose peaks span less than this distance (in PSF FWHM) " + "will be treated as one point at their mean position." + ) + findTrails = pexConfig.ConfigurableField( + target=FindGlintTrailsTask, + doc="Subtask to find lines of isolated positive peaks within a " + "footprint, which are kept together as one trailed source. Peak " + "positions are passed in units of the PSF FWHM, so its ``radius`` " + "and ``threshold`` are in FWHM rather than pixels; ``threshold`` " + "also limits how far any peak of a trail may lie from its line." + ) + + def setDefaults(self): + super().setDefaults() + self.findTrails.bad_flags = [] + # Larger than any footprint the size limits allow. + self.findTrails.radius = 1000.0 + self.findTrails.min_points = 3 + # Peaks of detected trails lie within ~0.25 FWHM of the trail's line. + self.findTrails.threshold = 0.25 class DeblendDifferenceFootprintsTask(pipeBase.Task): @@ -296,11 +602,23 @@ class DeblendDifferenceFootprintsTask(pipeBase.Task): original record and the rest become new diaSources. Footprint pixels are assigned to the cluster of highest Gaussian-weighted influence. Footprints that are too large or have too many peaks are left un-deblended. + + Clusters of positive peaks lying along one line are merged first, so that + a trailed source is not split into several diaSources: the trailed source + measurement handles them as long as their peaks share one footprint. """ ConfigClass = DeblendDifferenceFootprintsConfig _DefaultName = "deblendDifferenceFootprints" + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.makeSubtask("findTrails") + self._lineSchema = lsst.afw.table.SourceTable.makeMinimalSchema() + self._lineCentroidKey = lsst.afw.table.Point2DKey.addFields( + self._lineSchema, "centroid", "Peak position, as pixel coordinates divided by the PSF FWHM", "") + self._lineSchema.getAliasMap().set("slot_Centroid", "centroid") + def run(self, diaSources, difference): """Deblend the multi-cluster footprints of a diaSource catalog. @@ -315,7 +633,8 @@ def run(self, diaSources, difference): Difference image the footprints were detected on. Its PSF sets the clustering and influence scale, and its mask is updated in place: the ``NOT_DEBLENDED`` plane is set for footprints that can't be - deblended. + deblended, and for those left as detected because all their peaks + belong to one trailed source. Returns ------- @@ -326,7 +645,8 @@ def run(self, diaSources, difference): The updated catalog (`lsst.afw.table.SourceCatalog`). """ sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius() - linkingRadius = self.config.clusterRadius*sigma*_SIGMA_TO_FWHM + fwhm = sigma*gaussian_sigma_to_fwhm + linkingRadius = self.config.clusterRadius*fwhm isDeblendedKey = diaSources.schema.find("is_deblended").key isNegativeKey = diaSources.schema.find("is_negative").key difference.mask.addMaskPlane("NOT_DEBLENDED") @@ -335,6 +655,7 @@ def run(self, diaSources, difference): nFootprints = 0 nChildren = 0 nSkipped = 0 + nTrails = 0 # Snapshot the records up front so appending children below does not # disturb this iteration. for source in list(diaSources): @@ -354,20 +675,61 @@ def run(self, diaSources, difference): len(footprint.getPeaks())) footprint.getSpans().setMask(difference.mask, notDeblendedBitmask) continue - result = _deblend_footprint(footprint, clusters, sigma, self.config.maxContamination) + peaks = footprint.getPeaks() + negativeKey = peaks.getSchema().find(_NEGATIVE_PEAK_FLAG).key + positions = np.array([(peak.getIx(), peak.getIy()) for peak in peaks], dtype=float)/fwhm + negative = np.array([peak.get(negativeKey) for peak in peaks], dtype=bool) + bbox = footprint.getBBox() + connected = functools.partial(_segment_on_mask, _spans_to_bool(footprint.getSpans(), bbox), + np.array([bbox.getMinX(), bbox.getMinY()]), fwhm) + clusters, trails = _merge_trail_clusters(clusters, positions, negative, self._findLines, + connected, self.findTrails.config.min_points, + self.findTrails.config.threshold, + self.config.compactClusterSpan) + result = _deblend_footprint(footprint, clusters, sigma, self.config.maxContamination, trails) if result is None: continue - self._extractRecords(diaSources, source, result, isDeblendedKey, isNegativeKey) - nFootprints += 1 - nChildren += len(result.children) + nTrails += sum(result.trails) + if result.children: + self._extractRecords(diaSources, source, result, isDeblendedKey, isNegativeKey) + nFootprints += 1 + nChildren += len(result.children) + else: + # All of the peaks belong to one trailed source, left as + # detected for the trailed source measurement. + footprint.getSpans().setMask(difference.mask, notDeblendedBitmask) self.metadata["nDeblendedFootprints"] = nFootprints self.metadata["nDeblendedDiaSources"] = nChildren self.metadata["nSkipped"] = nSkipped - self.log.info("Deblended %d footprints into %d additional diaSources (%d skipped)", - nFootprints, nChildren, nSkipped) + self.metadata["nTrailedDiaSources"] = nTrails + self.log.info("Deblended %d footprints into %d additional diaSources " + "(%d skipped, %d trailed sources kept together)", + nFootprints, nChildren, nSkipped, nTrails) return pipeBase.Struct(diaSources=diaSources) + def _findLines(self, positions): + """Find lines of points with the ``findTrails`` subtask. + + Parameters + ---------- + positions : `numpy.ndarray` + ``(n, 2)`` array of the points' positions, as pixel coordinates + divided by the PSF FWHM. + + Returns + ------- + lines : `list` [`list` [`int`]] + Rows of ``positions`` on each line found. A row may appear on + more than one line. + """ + catalog = lsst.afw.table.SourceCatalog(self._lineSchema) + for x, y in positions: + catalog.addNew().set(self._lineCentroidKey, lsst.geom.Point2D(x, y)) + rows = {record.getId(): row for row, record in enumerate(catalog)} + trails = self.findTrails.run(catalog).trails + return [[rows[sourceId] for sourceId in trail["id"]] for trail in trails] + def _extractRecords(self, diaSources, source, result, isDeblendedKey, isNegativeKey): """Write one deblended footprint per cluster back to the catalog. diff --git a/tests/test_deblendDifferenceFootprints.py b/tests/test_deblendDifferenceFootprints.py index b3b582ff..2b467c73 100644 --- a/tests/test_deblendDifferenceFootprints.py +++ b/tests/test_deblendDifferenceFootprints.py @@ -33,6 +33,7 @@ DeblendDifferenceFootprintsTask, _cluster_peak_indices, _deblend_footprint, + _merge_trail_clusters, ) @@ -245,6 +246,19 @@ def test_no_merge_when_contamination_below_limit(self): self.assertIsNotNone(result) self.assertEqual(len(result.children), 1) + def test_trail_survives_contamination_merge(self): + """A cluster that absorbs a trailed source through the contamination + merge is still a trailed source, and is left as detected. + """ + footprint = makeFootprint(self.peakSchema, [(10, 20, False, 100.0), (22, 20, False, 100.0), + (34, 20, False, 100.0), (38, 20, False, 100.0)]) + # 4 px from the trail's last peak, with sigma=2: overlap exp(-1) > 0.05. + result = _deblend_footprint(footprint, [[3], [0, 1, 2]], sigma=2.0, trails=[False, True]) + self.assertIsNotNone(result) + self.assertEqual(result.children, []) + self.assertEqual(result.trails, [True]) + self.assertIs(result.parent, footprint) + class RunTestCase(lsst.utils.tests.TestCase): """Tests of the task-level ``run`` over a diaSource catalog.""" @@ -257,7 +271,8 @@ def makeExposure(self, psfSigma=2.0): # the NOT_DEBLENDED plane for any guardrail-skipped footprints. bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64)) exposure = afwImage.ExposureF(bbox) - exposure.setPsf(afwDetection.GaussianPsf(11, 11, psfSigma)) + exposure.setPsf(afwDetection.GaussianPsf(21, 21, psfSigma)) + exposure.variance.array[:] = 1.0 return exposure def makeCatalog(self, footprints): @@ -331,6 +346,8 @@ def test_run_skips_too_many_peaks(self): self.assertEqual(len(result.diaSources), 1) self.assertFalse(result.diaSources[0]["is_deblended"]) self.assertEqual(task.metadata["nSkipped"], 1) + # The size limits apply before the collinear peaks could form a trail. + self.assertEqual(task.metadata["nTrailedDiaSources"], 0) # The whole (full-frame) footprint is flagged NOT_DEBLENDED. notDeblended = difference.mask.getPlaneBitMask("NOT_DEBLENDED") self.assertTrue(np.all((difference.mask.array & notDeblended) != 0)) @@ -416,6 +433,338 @@ def test_bbox_backstop_skips_widely_spread_footprint(self): self.assertEqual(task.metadata["nSkipped"], 1) +class TrailTestCase(lsst.utils.tests.TestCase): + """Tests that clusters of positive peaks along one line are kept together + as a trailed source. + + With a 2 pixel Gaussian PSF the FWHM is ~4.7 pixels, so peaks 12 pixels + apart are beyond the ~9.4 pixel linking radius: each forms its own cluster, + and only the trail stage can keep them together. + """ + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + + def runTask(self, peaks, footprint=None, **configKwargs): + """Deblend one full-frame footprint carrying the given peaks. + + Parameters + ---------- + peaks : `list` of `tuple` + Each ``(ix, iy)`` for a positive peak of significance 20, or + ``(ix, iy, isNegative[, significance])``. + footprint : `lsst.afw.detection.Footprint`, optional + Footprint to deblend instead, carrying its own peaks. + """ + bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64)) + self.difference = afwImage.ExposureF(bbox) + self.difference.setPsf(afwDetection.GaussianPsf(21, 21, 2.0)) + config = DeblendDifferenceFootprintsTask.ConfigClass() + config.clusterRadius = 2.0 + for key, value in configKwargs.items(): + setattr(config, key, value) + catalog = afwTable.SourceCatalog(self.schema) + if footprint is None: + footprintPeaks = [(peak[0], peak[1], len(peak) > 2 and peak[2], + peak[3] if len(peak) > 3 else 20.0) for peak in peaks] + footprint = makeFootprint(self.peakSchema, footprintPeaks) + catalog.addNew().setFootprint(footprint) + self.task = DeblendDifferenceFootprintsTask(config=config) + return self.task.run(catalog, self.difference).diaSources + + def assertPeakCounts(self, diaSources, counts): + self.assertEqual(sorted(len(record.getFootprint().getPeaks()) for record in diaSources), counts) + + def notDeblendedPixels(self): + bit = self.difference.mask.getPlaneBitMask("NOT_DEBLENDED") + return (self.difference.mask.array & bit) != 0 + + def test_isolated_peaks_on_a_line(self): + """A line of isolated peaks stays one diaSource, left as detected.""" + diaSources = self.runTask([(14, 32), (26, 32), (38, 32), (50, 32)]) + + self.assertPeakCounts(diaSources, [4]) + self.assertFalse(diaSources[0]["is_deblended"]) + self.assertTrue(np.all(self.notDeblendedPixels())) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 1) + + def test_three_isolated_peaks_on_a_line(self): + """A line of exactly the minimum number of isolated peaks is found.""" + diaSources = self.runTask([(20, 32), (32, 32), (44, 32)]) + + self.assertPeakCounts(diaSources, [3]) + + def test_vertical_line(self): + """A vertical line of isolated peaks stays one diaSource.""" + diaSources = self.runTask([(32, 14), (32, 26), (32, 38), (32, 50)]) + + self.assertPeakCounts(diaSources, [4]) + + def test_two_isolated_peaks_split(self): + """Two isolated peaks are too few to make a line, so they deblend.""" + diaSources = self.runTask([(20, 32), (32, 32)]) + + self.assertPeakCounts(diaSources, [1, 1]) + self.assertTrue(all(record["is_deblended"] for record in diaSources)) + self.assertFalse(np.any(self.notDeblendedPixels())) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 0) + + def test_cluster_extends_along_its_line(self): + """A cluster of collinear positive peaks picks up isolated peaks on its + line beyond the linking radius. + + Two isolated peaks are too few for a line of isolated peaks, so only + the cluster's own line can pick them up. + """ + # The cluster's peaks are 8 pixels (1.7 FWHM) apart, enough to fix a + # direction. + diaSources = self.runTask([(10, 32), (18, 32), (30, 32), (42, 32)]) + + self.assertPeakCounts(diaSources, [4]) + + def test_cluster_and_one_peak_split(self): + """A cluster and one isolated peak on its line are too few clusters to + make a trail, so they deblend. + """ + diaSources = self.runTask([(10, 32), (18, 32), (30, 32)]) + + self.assertPeakCounts(diaSources, [1, 2]) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 0) + + def test_cluster_absorbs_collinear_clusters(self): + """A cluster absorbs whole clusters of positive peaks on its line.""" + diaSources = self.runTask([(4, 32), (12, 32), (24, 32), (32, 32), (44, 32), (52, 32)]) + + self.assertPeakCounts(diaSources, [6]) + + def test_negative_peak_blocks_absorption(self): + """A cluster on the line that holds a negative peak stays separate, + and the trail deblended from it is not flagged NOT_DEBLENDED. + """ + # The negative peak is on the line, within the linking radius of the + # last crest, so only its sign keeps the cluster out of the trail. + diaSources = self.runTask([(8, 32), (20, 32), (32, 32), (44, 32), (56, 32), (59, 32, True)]) + + self.assertPeakCounts(diaSources, [2, 4]) + self.assertTrue(all(record["is_deblended"] for record in diaSources)) + self.assertFalse(np.any(self.notDeblendedPixels())) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 1) + + def test_isolated_negative_peak_not_in_line(self): + """An isolated negative peak on the line is not part of the trail.""" + diaSources = self.runTask([(8, 32), (20, 32), (32, 32), (44, 32), (56, 32, True)]) + + self.assertPeakCounts(diaSources, [1, 4]) + + def test_two_peak_cluster_refit(self): + """A two-peak cluster picks up peaks that lie within tolerance of the + line refit to the trail, though not of the pair's own axis. + """ + # (40, 34) is 0.42 FWHM off the pair's axis, and fits best of the two. + diaSources = self.runTask([(16, 32), (24, 32), (40, 34), (56, 36)]) + + self.assertPeakCounts(diaSources, [4]) + + def test_two_peak_cluster_off_axis(self): + """A two-peak cluster does not pick up a peak far off its axis.""" + diaSources = self.runTask([(16, 32), (24, 32), (28, 46)]) + + self.assertPeakCounts(diaSources, [1, 2]) + + def test_compact_cluster_joins_line_as_a_point(self): + """A cluster spanning less than 1.5 FWHM joins a line of isolated peaks + at its mean position. + """ + diaSources = self.runTask([(8, 32), (20, 32), (31, 31), (33, 33), (44, 32)]) + + self.assertPeakCounts(diaSources, [5]) + + def test_compact_cluster_span(self): + """A cluster spread across a line joins it as a point only when its + span is below ``compactClusterSpan``. + """ + # The cluster's peaks are 6 pixels (1.3 FWHM) apart, each 0.64 FWHM + # off the line, but their mean is on it. + peaks = [(20, 32), (32, 29), (32, 35), (44, 32), (56, 32)] + self.assertPeakCounts(self.runTask(peaks), [5]) + self.assertPeakCounts(self.runTask(peaks, compactClusterSpan=1.0), [2, 3]) + + def test_compact_cluster_absorbed_at_its_mean(self): + """A trail grown from a wider cluster absorbs a compact cluster whose + mean lies on its line, though its peaks do not. + """ + # The middle pair's peaks are 6 pixels (1.3 FWHM) apart and 0.64 FWHM + # off the line; the first pair's are 8 pixels (1.7 FWHM) apart, so only + # the first fixes a direction. + peaks = [(8, 32), (16, 32), (28, 29), (28, 35), (40, 32), (52, 32)] + self.assertPeakCounts(self.runTask(peaks), [6]) + self.assertPeakCounts(self.runTask(peaks, compactClusterSpan=1.0), [2, 4]) + + def test_tolerance_in_fwhm(self): + """A peak 1 pixel (0.21 FWHM) off the line joins the trail, and one 2 + pixels (0.42 FWHM) off does not. + """ + line = [(12, 32), (32, 32), (42, 32), (52, 32)] + self.assertPeakCounts(self.runTask(line + [(22, 33)]), [5]) + self.assertPeakCounts(self.runTask(line + [(22, 34)]), [1, 4]) + + def test_off_line_peak_blocks_absorption(self): + """A cluster with a peak off the line stays separate.""" + # 8 pixels from the line: within the linking radius of the last crest, + # but well beyond the 0.25 FWHM tolerance. + diaSources = self.runTask([(8, 32), (20, 32), (32, 32), (44, 32), (56, 32), (56, 40)]) + + self.assertPeakCounts(diaSources, [2, 4]) + + def test_wide_gap_joins_line(self): + """A collinear peak far along the line still joins the trail when the + footprint covers the space between them. + """ + # 32 pixels, ~6.8 FWHM, from the nearest other peak. + diaSources = self.runTask([(4, 32), (16, 32), (28, 32), (60, 32)]) + + self.assertPeakCounts(diaSources, [4]) + + def notchedFootprint(self, columns, peaks): + """Build a footprint whose band along row 32 is cut by a notch, with + the two sides joined below it. + + Parameters + ---------- + columns : `tuple` [`int`, `int`] + First and last column of the notch. + peaks : `list` of `tuple` + Each ``(ix, iy)`` for a positive peak of significance 20. + """ + mask = afwImage.Mask(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64))) + bit = mask.getPlaneBitMask("DETECTED") + rows, cols = np.mgrid[0:64, 0:64] + notch = (rows <= 39) & (cols >= columns[0]) & (cols <= columns[1]) + mask.array[(rows >= 30) & (rows <= 44) & ~notch] = bit + footprint = afwDetection.Footprint(afwGeom.SpanSet.fromMask(mask, bit), self.peakSchema) + for ix, iy in peaks: + footprint.addPeak(ix, iy, 20.0)["significance"] = 20.0 + return footprint + + def test_footprint_gap_splits_line(self): + """A gap in the footprint splits a line of isolated peaks into a trail + on each side. + """ + footprint = self.notchedFootprint((30, 36), [(3, 32), (15, 32), (27, 32), + (39, 32), (51, 32), (63, 32)]) + + diaSources = self.runTask([], footprint=footprint) + + self.assertPeakCounts(diaSources, [3, 3]) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 2) + + def test_footprint_gap_stops_cluster_trail(self): + """A cluster does not extend its line across a gap in the footprint.""" + footprint = self.notchedFootprint((22, 26), [(10, 32), (18, 32), (30, 32), (42, 32)]) + + diaSources = self.runTask([], footprint=footprint) + + self.assertPeakCounts(diaSources, [1, 1, 2]) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 0) + + def test_trail_footprint_is_contiguous(self): + """A neighbor beside a gap between a trail's crests does not cut the + trail's footprint in two. + """ + # A narrow band holding the trail, with a bulge for the neighbor. + mask = afwImage.Mask(geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(64, 64))) + bit = mask.getPlaneBitMask("DETECTED") + rows, columns = np.mgrid[0:64, 0:64] + mask.array[(np.abs(rows - 32) <= 2) & (columns >= 2) & (columns <= 61)] = bit + mask.array[(rows - 37)**2 + (columns - 32)**2 <= 9] = bit + footprint = afwDetection.Footprint(afwGeom.SpanSet.fromMask(mask, bit), self.peakSchema) + for ix, iy in [(8, 32), (20, 32), (44, 32), (56, 32), (32, 37)]: + footprint.addPeak(ix, iy, 20.0)["significance"] = 20.0 + + diaSources = self.runTask([], footprint=footprint) + + self.assertPeakCounts(diaSources, [1, 4]) + trail = [record for record in diaSources if len(record.getFootprint().getPeaks()) == 4][0] + self.assertEqual(len(trail.getFootprint().getSpans().split()), 1) + + def test_split_line_keeps_every_peak(self): + """Every peak ends in exactly one cluster when the trail grown from a + line also absorbs a peak just off that line. + """ + # Positions in FWHM: the peak at 9.5 is off the line the finder + # returns, but within tolerance of the line refit to the trail. + positions = np.array([(0, 0), (3, 0), (6, 0), (9.5, 0.27), (13, 0), (16, 0), (19, 0), (22, -0.2)]) + + def findLines(points): + return [[0, 1, 2, 4, 5, 6, 7]] + + def connected(start, end): + return True + + clusters, trails = _merge_trail_clusters([[index] for index in range(8)], positions, + np.zeros(8, dtype=bool), findLines, connected, + min_points=3, tolerance=0.25, compact_span=1.5) + + self.assertEqual(sorted(peak for cluster in clusters for peak in cluster), list(range(8))) + self.assertEqual(sum(trails), 1) + + def test_short_trail_does_not_pivot(self): + """A short trail absorbs the far peaks of its own line before a nearer + peak off the line that a refit could also reach. + """ + # Positions in FWHM along and across a trail; the last peak is a + # bright neighbor about 1.2 FWHM off the line. + positions = np.array([(2.85, 0.15), (4.06, -0.15), (5.56, -0.15), + (-2.25, -0.15), (-4.5, 0.0), (0.0, 1.21)]) + + def findLines(points): + return [] + + def connected(start, end): + return True + + clusters, trails = _merge_trail_clusters([[0, 1, 2], [3], [4], [5]], positions, + np.zeros(6, dtype=bool), findLines, connected, + min_points=3, tolerance=0.25, compact_span=1.5) + + flags = {tuple(sorted(cluster)): trail for cluster, trail in zip(clusters, trails)} + self.assertEqual(flags, {(0, 1, 2, 3, 4): True, (5,): False}) + + def test_contamination_merge_into_trail(self): + """A faint peak merged into a trail by the contamination merge leaves + the footprint whole, as one trailed source. + """ + # 10 pixels off the line, beyond the linking radius, but 40 times + # fainter, so its PSF flux is contaminated by more than 5%. + diaSources = self.runTask([(14, 32), (26, 32), (38, 32), (50, 32), (26, 42, False, 0.5)]) + + self.assertPeakCounts(diaSources, [5]) + self.assertFalse(diaSources[0]["is_deblended"]) + self.assertTrue(np.all(self.notDeblendedPixels())) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 1) + + def test_two_trails_in_one_footprint(self): + """Two separate lines of peaks become two trailed sources.""" + diaSources = self.runTask([(8, 12), (20, 12), (32, 12), (44, 12), + (8, 52), (20, 52), (32, 52), (44, 52)]) + + self.assertPeakCounts(diaSources, [4, 4]) + self.assertFalse(np.any(self.notDeblendedPixels())) + self.assertEqual(self.task.metadata["nTrailedDiaSources"], 2) + self.assertEqual(self.task.metadata["nDeblendedFootprints"], 1) + self.assertEqual(self.task.metadata["nDeblendedDiaSources"], 1) + + def test_trail_and_separate_transient(self): + """A trail is kept together while an unrelated transient is deblended + from it, and neither is flagged NOT_DEBLENDED. + """ + diaSources = self.runTask([(14, 32), (26, 32), (38, 32), (50, 32), (32, 50)]) + + self.assertPeakCounts(diaSources, [1, 4]) + self.assertTrue(all(record["is_deblended"] for record in diaSources)) + self.assertFalse(np.any(self.notDeblendedPixels())) + + class MemoryTester(lsst.utils.tests.MemoryTestCase): pass