diff --git a/python/lsst/ip/diffim/__init__.py b/python/lsst/ip/diffim/__init__.py index 58f3879c..85f3cb6a 100755 --- a/python/lsst/ip/diffim/__init__.py +++ b/python/lsst/ip/diffim/__init__.py @@ -34,6 +34,7 @@ from .dipoleMeasurement import * from .getTemplate import * from .dipoleFitTask import * +from .deblendDifferenceFootprints import * from .detectAndMeasure import * from .imageDecorrelation import * from .imageMapReduce import * diff --git a/python/lsst/ip/diffim/deblendDifferenceFootprints.py b/python/lsst/ip/diffim/deblendDifferenceFootprints.py new file mode 100644 index 00000000..e9c5f344 --- /dev/null +++ b/python/lsst/ip/diffim/deblendDifferenceFootprints.py @@ -0,0 +1,747 @@ +# 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 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 + +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 + +# 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 _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. + 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): + """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") + # Fill the sparse matrix in "coordinate" (COO) format + 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 _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 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 + ---------- + clusters : `list` [`list` [`int`]] + Initial clusters, each a list of peak indices. + 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 + ------- + 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): + # 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 + 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, 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. A cluster that absorbs a trailed source is itself a + trailed source. + + 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. + 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 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() + + positions = np.array([(peak.getIx(), peak.getIy()) for peak in peaks], dtype=int) + amplitudes = np.array([peak["significance"] for peak in peaks], dtype=float) + + # 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) + # 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_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. + 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) + 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() + 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:], trails=trails) + + +class DeblendDifferenceFootprintsConfig(pexConfig.Config): + """Configuration for `DeblendDifferenceFootprintsTask`.""" + + clusterRadius = pexConfig.Field( + dtype=float, + 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, + default=100000, + doc="Skip deblending footprints larger than this, in pixels." + ) + maxFootprintBBoxArea = pexConfig.Field( + dtype=int, + 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=200, + 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." + ) + 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): + """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. + + 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. + + 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. 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, and for those left as detected because all their peaks + belong to one trailed source. + + 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() + 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") + notDeblendedBitmask = difference.mask.getPlaneBitMask("NOT_DEBLENDED") + + 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): + 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.getArea() > self.config.maxFootprintArea + 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 + 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 + 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.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. + + 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 ed6b9eca..0313f6f2 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,7 @@ from lsst.utils.timer import timeMethod from . import DipoleFitTask +from .deblendDifferenceFootprints import DeblendDifferenceFootprintsTask __all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] @@ -263,12 +263,20 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, ) doDeblend = pexConfig.Field( dtype=bool, - default=False, - doc="Deblend DIASources after detection?" + 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" ) deblend = pexConfig.ConfigurableField( - target=lsst.meas.deblender.SourceDeblendTask, - doc="Task to split blended sources into their components." + target=DeblendDifferenceFootprintsTask, + doc="Task to deblend merged footprints by clustering their peaks. " + "Only used if doDeblend=True." ) measurement = pexConfig.ConfigurableField( target=DipoleFitTask, @@ -287,11 +295,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, @@ -493,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"} @@ -522,6 +531,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 +542,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 +603,11 @@ 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.") + 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") @@ -786,17 +800,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 +831,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 +976,9 @@ 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) + + if self.config.doDeblend: + initialDiaSources = self.deblend.run(initialDiaSources, difference).diaSources else: initialDiaSources = sources @@ -1057,74 +1070,6 @@ 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. - - 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. - - 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``. - """ - - 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 - def _rejectBadMaskedDetections(self, initialDiaSources, mask): """Remove detections whose footprint peak lies on a bad-masked pixel. @@ -1704,25 +1649,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_deblendDifferenceFootprints.py b/tests/test_deblendDifferenceFootprints.py new file mode 100644 index 00000000..2b467c73 --- /dev/null +++ b/tests/test_deblendDifferenceFootprints.py @@ -0,0 +1,778 @@ +# 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, + _merge_trail_clusters, +) + + +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 + + +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.""" + + 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, max_contamination=0.05): + footprint = makeFootprint(self.peakSchema, peaks) + + clusters = _cluster_peak_indices(footprint.getPeaks(), linkingRadius) + result = _deblend_footprint(footprint, clusters, sigma=sigma, + max_contamination=max_contamination) + 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 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( + [(12, 20, False, 5.0), (16, 20, False, 100.0), (18, 20, True, 100.0)], + 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) + + 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.""" + + def setUp(self): + self.schema, self.peakSchema = makeSchemas() + + def makeExposure(self, psfSigma=2.0): + # 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(21, 21, psfSigma)) + exposure.variance.array[:] = 1.0 + 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() + difference = 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.""" + 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 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, difference) + + 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)) + + def test_run_skips_too_large_area(self): + """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 64**2 + difference = 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 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 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 + + +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 84b5d5aa..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, **kwargs): + badSourceFlags=None, clusterRadius=None, + doReplaceWithNoise=None, **kwargs): """Setup and configure the detection and measurement PipelineTask. Parameters @@ -148,6 +149,10 @@ 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 + if doReplaceWithNoise is not None: + config.measurement.doReplaceWithNoise = doReplaceWithNoise config.run_sattle = run_sattle @@ -193,7 +198,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, clusterRadius=1.5) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, sources, @@ -212,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 """ @@ -825,7 +888,8 @@ def _setup_sattle_tests(self): detectionTask = self._setup_detection(doDeblend=True, badSubtractionRatioThreshold=1., - doSkySources=False, run_sattle=True) + doSkySources=False, run_sattle=True, + clusterRadius=1.5) return science, matchedTemplate, difference, sources, detectionTask @@ -959,7 +1023,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, clusterRadius=1.5) # Run detection and check the results output = detectionTask.run(science, matchedTemplate, difference, score, sources, @@ -978,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. """ @@ -1406,57 +1532,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). @@ -1508,6 +1583,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) @@ -1541,10 +1619,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 ``clusterRadius``. + self.assertEqual((~result.diaSources["is_negative"]).sum(), 2) def testReorderPeaksBySignificance(self): """A merged footprint whose most significant peak is negative should @@ -1581,6 +1659,104 @@ 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( + 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.doDeblend = True + config.deblend.maxPeaks = 2 + task = detectAndMeasure.DetectAndMeasureTask(config=config) + # 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 maxPeaks. + diaSources = self._makeDeblendableCatalog(task.schema, bbox, [(8, 20), (20, 20), (32, 20)]) + result = task.deblend.run(diaSources, difference).diaSources + self.assertEqual(len(result), 1) + self.assertFalse(result[0]["is_deblended"]) + 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.doDeblend = True + config.deblend.maxFootprintArea = 100 + task = detectAndMeasure.DetectAndMeasureTask(config=config) + 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.deblend.run(diaSources, difference).diaSources + self.assertEqual(len(result), 1) + self.assertFalse(result[0]["is_deblended"]) + self.assertEqual(task.deblend.metadata["nSkipped"], 1) + def makeVisitInfo(id=1): """Return a non-NaN visitInfo."""