From 8cab480f6f478964d66060db09f4a7edc085b519 Mon Sep 17 00:00:00 2001 From: fred3m Date: Wed, 23 Sep 2026 07:58:31 -0700 Subject: [PATCH 1/6] Add starlet-based MultiResolutionDetectionTask Detect sources with the scarlet lite starlet transform as a compensated filter, on single-band Exposures or MultibandExposures, as a drop-in replacement for SourceDetectionTask that needs no background subtraction. The detection intermediate products (peaks, candidates, positions) are returned as astropy tables. Generated with AI Co-Authored-By: SLAC AI --- python/lsst/meas/algorithms/__init__.py | 1 + .../algorithms/multiResolutionDetection.py | 569 ++++++++++++++++++ tests/test_multiResolutionDetection.py | 138 +++++ ups/meas_algorithms.table | 1 + 4 files changed, 709 insertions(+) create mode 100644 python/lsst/meas/algorithms/multiResolutionDetection.py create mode 100644 tests/test_multiResolutionDetection.py diff --git a/python/lsst/meas/algorithms/__init__.py b/python/lsst/meas/algorithms/__init__.py index 59e9742b..2f45fa90 100644 --- a/python/lsst/meas/algorithms/__init__.py +++ b/python/lsst/meas/algorithms/__init__.py @@ -54,6 +54,7 @@ from .reserveSourcesTask import * from .skyObjects import * from .dynamicDetection import * +from .multiResolutionDetection import * from .makePsfCandidates import * from .stamps import * from .accumulator_mean_stack import * diff --git a/python/lsst/meas/algorithms/multiResolutionDetection.py b/python/lsst/meas/algorithms/multiResolutionDetection.py new file mode 100644 index 00000000..498e91e7 --- /dev/null +++ b/python/lsst/meas/algorithms/multiResolutionDetection.py @@ -0,0 +1,569 @@ +# This file is part of meas_algorithms. +# +# 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 . + +from __future__ import annotations + +__all__ = ("MultiResolutionDetectionConfig", "MultiResolutionDetectionTask") + +from contextlib import contextmanager +from typing import Any + +import numpy as np +from astropy.table import Table +from numpy.lib import recfunctions as rfn + +import lsst.afw.detection as afwDet +import lsst.afw.geom as afwGeom +import lsst.afw.image as afwImage +import lsst.afw.math as afwMath +import lsst.afw.table as afwTable +import lsst.geom as geom +import lsst.pex.config as pexConfig +import lsst.pipe.base as pipeBase +import lsst.scarlet.lite as scl +from lsst.utils.timer import timeMethod + +from .detection import SourceDetectionConfig, SourceDetectionTask + + +# Conversion from Gaussian sigma to full width at half maximum. +# This is also the default minimum FWHM when the exposure has no usable PSF. +SIGMA_TO_FWHM = 2.0*np.sqrt(2.0*np.log(2.0)) + + +def scarletFootprintToAfw(footprint: scl.detect.Footprint) -> afwDet.Footprint: + """Convert a scarlet lite footprint into an afw footprint. + + Parameters + ---------- + footprint: + The scarlet footprint to convert. + + Returns + ------- + newFootprint: + The converted afw footprint, with the scarlet peaks copied over. + + Notes + ----- + This duplicates ``scarletFootprintToAfw`` from meas_extensions_scarlet. + That package depends on meas_algorithms, so importing it here would create + a circular dependency. + """ + xy0 = geom.Point2I(int(footprint.bbox.origin[1]), int(footprint.bbox.origin[0])) + data = afwImage.Mask(footprint.data.astype(np.int32), xy0=xy0) + spans = afwGeom.SpanSet.fromMask(data) + newFootprint = afwDet.Footprint(spans) + for peak in footprint.peaks: + newFootprint.addPeak(peak.x, peak.y, peak.flux) + return newFootprint + + +class MultiResolutionDetectionConfig(SourceDetectionConfig): + """Configuration parameters for the MultiResolutionDetectionTask. + """ + starletScales = pexConfig.Field( + doc="The maximum number of starlet scales to use in the multi-resolution detection.", + dtype=int, default=3, + ) + firstScale = pexConfig.Field( + doc="The first starlet scale to use in the multi-resolution detection.", + dtype=int, default=1, + ) + generation = pexConfig.Field( + doc="The generation of the starlet transform to use (1 or 2).", + dtype=int, default=2, + ) + minPixels = pexConfig.Field( + doc="Detected footprints with fewer than the specified number of pixels will be ignored.", + dtype=int, default=8, + ) + peakThreshold = pexConfig.Field( + doc="Threshold for detecting peaks within footprints, in sigma.", + dtype=float, default=5.0, + ) + thresholdValue = pexConfig.Field( + doc="Threshold for detecting footprints, in sigma.", + dtype=float, default=2.0, + ) + saddleThreshold = pexConfig.Field( + doc="Minimum prominence of a peak above the saddle to a brighter peak, in sigma.", + dtype=float, default=3.0, + ) + + def setDefaults(self) -> None: + super().setDefaults() + # Starlets act as a compensated filter, so none of the background + # subtraction machinery of the base task is used or wanted. + self.reEstimateBackground = False + self.doTempLocalBackground = False + self.doTempWideBackground = False + + def validate(self) -> None: + super().validate() + if self.reEstimateBackground: + raise pexConfig.FieldValidationError( + self.__class__.reEstimateBackground, self, + "reEstimateBackground must be False for MultiResolutionDetectionTask" + ) + + +class MultiResolutionDetectionTask(SourceDetectionTask): + """Detect peaks and footprints of sources in an image. + + Parameters + ---------- + schema: + Schema object used to create the output `lsst.afw.table.SourceCatalog`. + **kwds: + Keyword arguments passed to `lsst.pipe.base.Task.__init__`. + + Notes + ----- + This task uses the starlet transform as a compensated filter to detect + sources at multiple scales, and (optionally) multiple bands. + Once the starlet coefficients are built, a significance map is made that + converts images into units of sigma and creates a chi^2 image + (if multiple bands are used). Candidate peaks are then detected in each + significance map using a watershed algorithm. Next the peaks are merged + across bands and scales to produce a final catalog of detected sources. + Finally the peaks are added to footprints and are optionally grown + according to the configuration. + + Because the starlet detection accounts for the noise in the exposure + itself, the exposure must not have its background re-estimated or its + variance rescaled before detection. + """ + ConfigClass = MultiResolutionDetectionConfig + _DefaultName = "multiResolutionDetection" + + def __init__(self, schema: afwTable.Schema | None, **kwds: Any) -> None: + # Skip SourceDetectionTask.__init__ so that none of the background + # subtasks are constructed; only the negative-detection flag is needed. + pipeBase.Task.__init__(self, **kwds) + if schema is not None and self.config.thresholdPolarity == "both": + self.negativeFlagKey = schema.addField( + "is_negative", type="Flag", + doc="Set if source peak was detected as negative." + ) + else: + if self.config.thresholdPolarity == "both": + self.log.warning("Detection polarity set to 'both', but no flag will be " + "set to distinguish between positive and negative detections") + self.negativeFlagKey = None + + @timeMethod + def run( + self, + table: afwTable.SourceTable, + exposure: afwImage.Exposure | afwImage.MultibandExposure, + doSmooth: bool = True, + sigma: float | None = None, + clearMask: bool = True, + expId: int | None = None, + background: afwMath.BackgroundList | None = None, + backgroundToPhotometricRatio: afwImage.Image | None = None, + ) -> pipeBase.Struct: + r"""Detect sources and return a catalog of detections. + + Parameters + ---------- + table: + Table object that will be used to create the SourceCatalog. + exposure: + Exposure to process; the DETECTED mask plane will be set in-place. + May be a single-band ``Exposure`` or a multi-band + ``MultibandExposure``. + doSmooth: + Unused; retained for compatibility with the base + ``SourceDetectionTask``. + sigma: + Gaussian sigma of the PSF (pixels), used to grow detections; if + `None` it is measured from the PSF of the ``exposure``. + clearMask: + Clear DETECTED{,_NEGATIVE} planes before running detection. + expId: + Unused; retained for compatibility with the base task. + background: + Unused; retained for compatibility with the base task. + backgroundToPhotometricRatio: + Unused; retained for compatibility with the base task. + + Returns + ------- + result : `lsst.pipe.base.Struct` + See `detectFootprints` for the contents, plus: + + ``sources`` + Detected sources on the exposure. + (`lsst.afw.table.SourceCatalog`) + + Raises + ------ + ValueError + Raised if flags.negative is needed, but isn't in table's schema. + + Notes + ----- + This override exists only to accept a `MultibandExposure` and to + document the parameters that this task actually uses; the base + implementation already delegates the detection to the overridden + `detectFootprints`. + """ + return super().run( + table=table, + exposure=exposure, + doSmooth=doSmooth, + sigma=sigma, + clearMask=clearMask, + expId=expId, + background=background, + backgroundToPhotometricRatio=backgroundToPhotometricRatio, + ) + + @timeMethod + def detectFootprints( + self, + exposure: afwImage.Exposure | afwImage.MultibandExposure, + doSmooth: bool = False, + sigma: float | None = None, + clearMask: bool = True, + expId: int | None = None, + background: afwMath.BackgroundList | None = None, + backgroundToPhotometricRatio: afwImage.Image | None = None, + ) -> pipeBase.Struct: + """Detect footprints on an exposure. + + Parameters + ---------- + exposure: + Exposure to process; the DETECTED{,_NEGATIVE} mask plane will be + set in-place. May be a single-band ``Exposure`` or a multi-band + ``MultibandExposure``. + doSmooth: + Unused; retained for compatibility with the base task. + sigma: + Gaussian sigma of the PSF (pixels), used to grow detections; if + `None` it is measured from the PSF of the ``exposure``. + clearMask: + Clear both DETECTED and DETECTED_NEGATIVE planes before running + detection. + expId: + Unused; retained for compatibility with the base task. + background: + Unused; retained for compatibility with the base task. + backgroundToPhotometricRatio: + Unused; retained for compatibility with the base task. + + Returns + ------- + results : `lsst.pipe.base.Struct` + See `scarletToAfwResults` for the contents. In addition + `finalizeFootprints` adds the ``numPos``, ``numPosPeaks``, + ``numNeg`` and ``numNegPeaks`` counts. + """ + multiband = hasattr(exposure, "bands") + x0, y0 = exposure.getBBox().getMin() + + if multiband: + singles = list(exposure.singles) + if clearMask: + for single in singles: + self.clearMask(single.mask) + self.removeBadPixels(single.maskedImage) + fwhm = self._getFwhm(singles, sigma) + images = exposure.image.array + variance = exposure.variance.array + mask = singles[0].mask + else: + maskedImage = exposure.maskedImage + if clearMask: + self.clearMask(maskedImage.getMask()) + self.removeBadPixels(maskedImage) + fwhm = self._getFwhm([exposure], sigma) + images = exposure.image.array[None] + variance = exposure.variance.array[None] + mask = maskedImage.getMask() + + positive = None + negative = None + if self.config.thresholdPolarity in ("positive", "both"): + positive = self._detectPeaks(images, variance, (y0, x0), fwhm) + if self.config.thresholdPolarity in ("negative", "both"): + negative = self._detectPeaks(-images, variance, (y0, x0), fwhm) + + results = self.scarletToAfwResults(exposure, positive, negative) + + # Grow the footprints and set the DETECTED mask planes. The footprints + # live in a single coordinate system, so growing once is enough; the + # mask plane is then set on every band. + sigmaToGrow = fwhm/SIGMA_TO_FWHM + self.finalizeFootprints(mask, results, sigmaToGrow, factor=1.0, factorNeg=1.0) + if multiband: + for single in singles[1:]: + for polarity, maskName in (("positive", "DETECTED"), ("negative", "DETECTED_NEGATIVE")): + fpSet = getattr(results, polarity) + if fpSet is not None: + fpSet.setMask(single.mask, maskName) + self.clearUnwantedResults(single.mask, results) + + self.clearUnwantedResults(mask, results) + + return results + + def _getFwhm(self, singles: list[afwImage.Exposure], sigma: float | None) -> float: + """Measure the PSF FWHM to use for growing and linking detections. + + Parameters + ---------- + singles: + The single-band exposures to measure the PSF from. For a multi-band + exposure the widest PSF is used, since it sets the most conservative + linking radius across bands. + sigma: + Gaussian sigma of the PSF (pixels); if `None` it is measured from + each exposure's PSF. + + Returns + ------- + fwhm: + The PSF full width at half maximum, in pixels. + """ + if sigma is not None: + return sigma*SIGMA_TO_FWHM + sigmas = [] + for single in singles: + try: + psf = self.getPsf(single, sigma=None) + sigmas.append(psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()) + except Exception: + continue + if not sigmas: + return SIGMA_TO_FWHM + return max(max(sigmas)*SIGMA_TO_FWHM, SIGMA_TO_FWHM) + + def _detectPeaks( + self, + images: np.ndarray, + variance: np.ndarray, + origin: tuple[int, int], + fwhm: float, + ) -> scl.detect.PeakDetectionResult: + """Run the scarlet lite peak detection on a set of images. + + Parameters + ---------- + images: + The images to detect on, with shape ``(n_bands, Ny, Nx)``. + variance: + The per-pixel variance, with the same shape as ``images``. + origin: + The ``(y, x)`` location of the lower corner of the images. + fwhm: + The PSF full width at half maximum, in pixels. + + Returns + ------- + result: + The detected peaks and the intermediate detection products. + """ + return scl.detect.detect_peaks( + images=images, + variance=variance, + scales=self.config.starletScales, + generation=self.config.generation, + first_scale=self.config.firstScale, + origin=origin, + min_separation=0, + min_area=self.config.minPixels, + peak_thresh=self.config.peakThreshold, + footprint_thresh=self.config.thresholdValue, + psf_fwhm=fwhm, + kappa=self.config.saddleThreshold, + ) + + def scarletToAfwResults( + self, + exposure: afwImage.Exposure | afwImage.MultibandExposure, + positive: scl.detect.PeakDetectionResult | None = None, + negative: scl.detect.PeakDetectionResult | None = None, + ) -> pipeBase.Struct: + """Convert scarlet detection results into an afw detection Struct. + + Parameters + ---------- + exposure: + The exposure detection was run on, used for the footprint region. + positive: + The positive polarity detection result. + negative: + The negative polarity detection result. + + Returns + ------- + results : `lsst.pipe.base.Struct` + A `~lsst.pipe.base.Struct` containing: + + ``positive`` + Positive polarity footprints. + (`lsst.afw.detection.FootprintSet` or `None`) + ``negative`` + Negative polarity footprints. + (`lsst.afw.detection.FootprintSet` or `None`) + ``factor`` + Multiplication factor applied to the positive threshold; always + 1.0. (`float`) + ``factorNeg`` + Multiplication factor applied to the negative threshold; always + 1.0. (`float`) + ``background`` + Always `None`; this task does not estimate a background. + ``peaks`` + The final detections, one row per source, with a ``polarity`` + column. (`astropy.table.Table`) + ``candidates`` + The peak candidates, before grouping, with a ``polarity`` + column. (`astropy.table.Table`) + ``positions`` + The unique candidate positions, with a ``polarity`` column. + (`astropy.table.Table`) + ``significanceMap`` + The per-scale significance map of the positive detection (or the + negative detection if only negative was run). (`numpy.ndarray` + or `None`) + ``starlets`` + The multiband starlet coefficients, as for ``significanceMap``. + (`numpy.ndarray` or `None`) + ``sigma`` + The per-scale, per-band coefficient noise, as for + ``significanceMap``. (`numpy.ndarray` or `None`) + """ + bbox = exposure.getBBox() + results = pipeBase.Struct( + positive=self._makeFootprintSet(positive, bbox), + negative=self._makeFootprintSet(negative, bbox), + factor=1.0, + factorNeg=1.0, + background=None, + peaks=self._joinPolarities(positive, negative, "peaks"), + candidates=self._joinPolarities(positive, negative, "candidates"), + positions=self._joinPolarities(positive, negative, "positions"), + ) + primary = positive if positive is not None else negative + if primary is not None: + results.significanceMap = primary.significance_map + results.starlets = primary.starlets + results.sigma = primary.sigma + else: + results.significanceMap = None + results.starlets = None + results.sigma = None + return results + + @staticmethod + def _makeFootprintSet( + result: scl.detect.PeakDetectionResult | None, + bbox: geom.Box2I, + ) -> afwDet.FootprintSet | None: + """Build an afw FootprintSet from a scarlet detection result. + + Parameters + ---------- + result: + The detection result to convert. + bbox: + The region the footprints were detected in. + + Returns + ------- + footprintSet: + The footprints, or `None` if ``result`` is `None`. + """ + if result is None: + return None + footprintSet = afwDet.FootprintSet(bbox) + footprintSet.setFootprints([scarletFootprintToAfw(fp) for fp in result.footprints]) + return footprintSet + + @staticmethod + def _joinPolarities( + positive: scl.detect.PeakDetectionResult | None, + negative: scl.detect.PeakDetectionResult | None, + attr: str, + ) -> Table: + """Join a structured-array product across detection polarities. + + Parameters + ---------- + positive: + The positive polarity detection result. + negative: + The negative polarity detection result. + attr: + The name of the structured-array attribute to join (``peaks``, + ``candidates`` or ``positions``). + + Returns + ------- + joined: + The concatenated rows as an `astropy.table.Table` with an added + ``polarity`` column (``1`` for positive, ``-1`` for negative). The + table is empty if neither result is set. + """ + arrays = [] + for result, polarity in ((positive, 1), (negative, -1)): + if result is None: + continue + array = getattr(result, attr) + array = rfn.append_fields( + array, "polarity", np.full(len(array), polarity, dtype=np.int32), usemask=False + ) + arrays.append(array) + if not arrays: + return Table() + if len(arrays) == 1: + return Table(arrays[0]) + return Table(np.concatenate(arrays)) + + def applyTempLocalBackground(self, exposure, middle, results): + raise NotImplementedError("applyTempLocalBackground is not used in MultiResolutionDetectionTask") + + def convolveImage(self, maskedImage, psf, doSmooth=True): + raise NotImplementedError("convolveImage is not used in MultiResolutionDetectionTask") + + def applyThreshold(self, middle, bbox, factor=1.0, factorNeg=None): + raise NotImplementedError("applyThreshold is not used in MultiResolutionDetectionTask") + + def reEstimateBackground(self, maskedImage, backgrounds, backgroundToPhotometricRatio=None): + raise NotImplementedError("reEstimateBackground is not used in MultiResolutionDetectionTask") + + def setPeakSignificance(self, exposure, footprints, threshold, negative=False): + raise NotImplementedError("setPeakSignificance is not used in MultiResolutionDetectionTask") + + def makeThreshold(self, image, thresholdParity, factor=1.0): + raise NotImplementedError("makeThreshold is not used in MultiResolutionDetectionTask") + + def updatePeaks(self, fpSet, image, threshold): + raise NotImplementedError("updatePeaks is not used in MultiResolutionDetectionTask") + + @contextmanager + def tempWideBackgroundContext(self, exposure): + raise NotImplementedError("tempWideBackgroundContext is not used in MultiResolutionDetectionTask") diff --git a/tests/test_multiResolutionDetection.py b/tests/test_multiResolutionDetection.py new file mode 100644 index 00000000..d90be86a --- /dev/null +++ b/tests/test_multiResolutionDetection.py @@ -0,0 +1,138 @@ +# This file is part of meas_algorithms. +# +# 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 + +from astropy.table import Table + +import lsst.afw.detection as afwDet +import lsst.afw.image as afwImage +import lsst.afw.table as afwTable +import lsst.geom +import lsst.pex.config +import lsst.utils.tests +from lsst.meas.algorithms import MultiResolutionDetectionConfig, MultiResolutionDetectionTask +from lsst.meas.algorithms.testUtils import plantSources + + +class MultiResolutionDetectionTaskTestCase(lsst.utils.tests.TestCase): + + def setUp(self) -> None: + self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(128, 128)) + self.kwid = 11 + self.sky = 100.0 + self.starSigma = 2.0 + # (x, y, counts, sigma) for each planted star. + self.coords = [(30, 30, 5000, self.starSigma), + (80, 90, 8000, self.starSigma), + (60, 20, 3000, self.starSigma)] + + def _makeExposure(self, offset: int = 0) -> afwImage.Exposure: + coords = [(x, y, counts + offset, sigma) for x, y, counts, sigma in self.coords] + return plantSources(self.bbox, self.kwid, self.sky, coords, addPoissonNoise=True) + + def testConfigValidation(self) -> None: + """The config must default to no background re-estimation and reject + turning it back on.""" + config = MultiResolutionDetectionConfig() + self.assertFalse(config.reEstimateBackground) + self.assertFalse(config.doTempLocalBackground) + self.assertFalse(config.doTempWideBackground) + config.reEstimateBackground = True + with self.assertRaises(lsst.pex.config.FieldValidationError): + config.validate() + + def testSingleBand(self) -> None: + """Detection on a single-band exposure recovers the planted sources.""" + exposure = self._makeExposure() + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema) + table = afwTable.SourceTable.make(schema) + results = task.run(table, exposure) + + self.assertGreaterEqual(len(results.sources), len(self.coords)) + self.assertEqual(results.numPos, len(results.positive.getFootprints())) + self.assertIsInstance(results.positive, afwDet.FootprintSet) + self.assertIsNone(results.negative) + self.assertIsNone(results.background) + + # The scarlet intermediate products are attached to the result as + # astropy tables. + for attr in ("peaks", "candidates", "positions"): + table = getattr(results, attr) + self.assertIsInstance(table, Table) + self.assertIn("polarity", table.colnames) + self.assertGreater(len(table), 0) + self.assertIsNotNone(results.significanceMap) + self.assertIsNotNone(results.starlets) + + # The DETECTED plane is set on the exposure. + mask = exposure.mask + detected = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0 + self.assertGreater(detected.sum(), 0) + + def testMultiBand(self) -> None: + """Detection on a MultibandExposure produces a single merged catalog + and sets the DETECTED plane in every band.""" + bands = ["g", "r", "i"] + exposures = [self._makeExposure(offset=1000 * i) for i in range(len(bands))] + mExposure = afwImage.MultibandExposure.fromExposures(bands, exposures) + + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema) + table = afwTable.SourceTable.make(schema) + results = task.run(table, mExposure) + + self.assertGreaterEqual(len(results.sources), len(self.coords)) + self.assertEqual(results.significanceMap.shape[1], len(bands) + 1) + self.assertEqual(results.starlets.shape[1], len(bands)) + for single in mExposure.singles: + mask = single.mask + detected = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0 + self.assertGreater(detected.sum(), 0) + + def testBothPolarities(self) -> None: + """With thresholdPolarity='both', negative footprints are returned and + the peak table records both polarities.""" + exposure = self._makeExposure() + config = MultiResolutionDetectionConfig() + config.thresholdPolarity = "both" + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + table = afwTable.SourceTable.make(schema) + results = task.run(table, exposure) + + self.assertIsInstance(results.positive, afwDet.FootprintSet) + self.assertIsInstance(results.negative, afwDet.FootprintSet) + self.assertIsNotNone(task.negativeFlagKey) + + +class TestMemory(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/ups/meas_algorithms.table b/ups/meas_algorithms.table index 57cd6961..d0a6c1d1 100644 --- a/ups/meas_algorithms.table +++ b/ups/meas_algorithms.table @@ -13,6 +13,7 @@ setupRequired(cpputils) setupRequired(kht) setupRequired(obs_base) setupRequired(shapelet) +setupRequired(scarlet_lite) envPrepend(LD_LIBRARY_PATH, ${PRODUCT_DIR}/lib) envPrepend(DYLD_LIBRARY_PATH, ${PRODUCT_DIR}/lib) From 13860ff05cf4516d98297296c25ad9ba8caccde1 Mon Sep 17 00:00:00 2001 From: Fred Moolekamp Date: Thu, 24 Sep 2026 16:33:07 -0400 Subject: [PATCH 2/6] Detect footprints on the image, not the wavelets Starlets are a compensated filter with no response to smooth flux, so the coefficient footprints of an extended source stop well inside it. Footprints now come from a chi^2 coadd of the PSF-correlated bands, with the starlet peaks inserted into whichever footprint contains them, and peakValue read off the image the footprint came from. That leaves an image in counts to re-estimate a background from, so the background handling of the base task is used rather than suppressed. Co-Authored-By: Claude Opus 5 (1M context) --- .../algorithms/multiResolutionDetection.py | 710 ++++++++++++++---- tests/test_multiResolutionDetection.py | 273 ++++++- 2 files changed, 827 insertions(+), 156 deletions(-) diff --git a/python/lsst/meas/algorithms/multiResolutionDetection.py b/python/lsst/meas/algorithms/multiResolutionDetection.py index 498e91e7..84817551 100644 --- a/python/lsst/meas/algorithms/multiResolutionDetection.py +++ b/python/lsst/meas/algorithms/multiResolutionDetection.py @@ -24,14 +24,13 @@ __all__ = ("MultiResolutionDetectionConfig", "MultiResolutionDetectionTask") from contextlib import contextmanager -from typing import Any +from typing import cast import numpy as np from astropy.table import Table from numpy.lib import recfunctions as rfn import lsst.afw.detection as afwDet -import lsst.afw.geom as afwGeom import lsst.afw.image as afwImage import lsst.afw.math as afwMath import lsst.afw.table as afwTable @@ -48,33 +47,9 @@ # This is also the default minimum FWHM when the exposure has no usable PSF. SIGMA_TO_FWHM = 2.0*np.sqrt(2.0*np.log(2.0)) - -def scarletFootprintToAfw(footprint: scl.detect.Footprint) -> afwDet.Footprint: - """Convert a scarlet lite footprint into an afw footprint. - - Parameters - ---------- - footprint: - The scarlet footprint to convert. - - Returns - ------- - newFootprint: - The converted afw footprint, with the scarlet peaks copied over. - - Notes - ----- - This duplicates ``scarletFootprintToAfw`` from meas_extensions_scarlet. - That package depends on meas_algorithms, so importing it here would create - a circular dependency. - """ - xy0 = geom.Point2I(int(footprint.bbox.origin[1]), int(footprint.bbox.origin[0])) - data = afwImage.Mask(footprint.data.astype(np.int32), xy0=xy0) - spans = afwGeom.SpanSet.fromMask(data) - newFootprint = afwDet.Footprint(spans) - for peak in footprint.peaks: - newFootprint.addPeak(peak.x, peak.y, peak.flux) - return newFootprint +# The scarlet lite name for each threshold type, so that the peaks and the +# footprints are standardized by the same noise. +VARIANCE_MODE = {"pixel_stdev": "pixel", "stdev": "median"} class MultiResolutionDetectionConfig(SourceDetectionConfig): @@ -101,7 +76,12 @@ class MultiResolutionDetectionConfig(SourceDetectionConfig): dtype=float, default=5.0, ) thresholdValue = pexConfig.Field( - doc="Threshold for detecting footprints, in sigma.", + doc="Threshold for detecting footprints on the chi^2 detection image, in sigma. Whether that " + "is the per-pixel or the per-band sigma is set by thresholdType.", + dtype=float, default=2.0, + ) + filterThreshold = pexConfig.Field( + doc="Threshold for the starlet coefficient footprints that peaks are detected on, in sigma.", dtype=float, default=2.0, ) saddleThreshold = pexConfig.Field( @@ -111,18 +91,34 @@ class MultiResolutionDetectionConfig(SourceDetectionConfig): def setDefaults(self) -> None: super().setDefaults() - # Starlets act as a compensated filter, so none of the background - # subtraction machinery of the base task is used or wanted. - self.reEstimateBackground = False + # Peaks come from the starlet coefficients, which are unaffected by a + # low frequency background, so neither local nor wide temporary + # background subtraction is implemented here. However, the background + # re-estimation of the base task is used, since footprints are + # detected on the image itself. self.doTempLocalBackground = False self.doTempWideBackground = False + # Unlike SourceDetectionTask, we want the default thresholdType to + # be 'stdev', since the detection gains on starlet coefficients by + # using `pixel_stdev` is minimal at the cost of significantly more + # memory usage. + self.thresholdType = "stdev" + def validate(self) -> None: super().validate() - if self.reEstimateBackground: + for name in ("doTempLocalBackground", "doTempWideBackground"): + if getattr(self, name): + raise pexConfig.FieldValidationError( + getattr(self.__class__, name), self, + f"{name} must be False for MultiResolutionDetectionTask" + ) + if self.thresholdType not in ("pixel_stdev", "stdev"): raise pexConfig.FieldValidationError( - self.__class__.reEstimateBackground, self, - "reEstimateBackground must be False for MultiResolutionDetectionTask" + self.__class__.thresholdType, self, + "thresholdType must be 'pixel_stdev' or 'stdev' for MultiResolutionDetectionTask: " + "the bands are divided by their noise before they are coadded, so the chi^2 detection " + "image is in sigma and there is no image in counts to threshold" ) @@ -145,31 +141,28 @@ class MultiResolutionDetectionTask(SourceDetectionTask): (if multiple bands are used). Candidate peaks are then detected in each significance map using a watershed algorithm. Next the peaks are merged across bands and scales to produce a final catalog of detected sources. - Finally the peaks are added to footprints and are optionally grown - according to the configuration. - Because the starlet detection accounts for the noise in the exposure - itself, the exposure must not have its background re-estimated or its - variance rescaled before detection. + The footprints come from the image itself rather than from the starlet + coefficients. Footprints built off of the compensated filter intentionally + remove signal from lower scales, so the footprints do not capture the + full extent of the source. Instead every band is correlated with its own + PSF, the bands are coadded into a chi^2 detection image in units of sigma, + and that image is thresholded as in the base task. + Each peak is then inserted into the footprint that contains it, and + footprints with no peak are dropped. + + References + ---------- + .. [1] Lupton, R. H., "chi^2 coadds", 2021, section 1.2.2. + + .. [2] Szalay, A. S., Connolly, A. J., and Szokoly, G. P., "Simultaneous + Multicolor Detection of Faint Galaxies in the Hubble Deep Field", + The Astronomical Journal, vol. 117, no. 1, pp. 68-74, 1999. + doi:10.1086/300689. """ ConfigClass = MultiResolutionDetectionConfig _DefaultName = "multiResolutionDetection" - def __init__(self, schema: afwTable.Schema | None, **kwds: Any) -> None: - # Skip SourceDetectionTask.__init__ so that none of the background - # subtasks are constructed; only the negative-detection flag is needed. - pipeBase.Task.__init__(self, **kwds) - if schema is not None and self.config.thresholdPolarity == "both": - self.negativeFlagKey = schema.addField( - "is_negative", type="Flag", - doc="Set if source peak was detected as negative." - ) - else: - if self.config.thresholdPolarity == "both": - self.log.warning("Detection polarity set to 'both', but no flag will be " - "set to distinguish between positive and negative detections") - self.negativeFlagKey = None - @timeMethod def run( self, @@ -179,7 +172,7 @@ def run( sigma: float | None = None, clearMask: bool = True, expId: int | None = None, - background: afwMath.BackgroundList | None = None, + background: afwMath.BackgroundList | list[afwMath.BackgroundList] | None = None, backgroundToPhotometricRatio: afwImage.Image | None = None, ) -> pipeBase.Struct: r"""Detect sources and return a catalog of detections. @@ -193,19 +186,25 @@ def run( May be a single-band ``Exposure`` or a multi-band ``MultibandExposure``. doSmooth: - Unused; retained for compatibility with the base - ``SourceDetectionTask``. + If True, convolve each band with a Gaussian of width ``sigma``, or + of the measured PSF width of ``exposure``, before the footprints + are detected. The peaks are unaffected either way. sigma: - Gaussian sigma of the PSF (pixels), used to grow detections; if - `None` it is measured from the PSF of the ``exposure``. + Gaussian sigma of the PSF (pixels), used for smoothing and to grow + footprint detections; if `None` it is measured from the PSF of the + ``exposure``. clearMask: Clear DETECTED{,_NEGATIVE} planes before running detection. expId: Unused; retained for compatibility with the base task. background: - Unused; retained for compatibility with the base task. + Background that was already subtracted from the exposure; will be + modified in-place if ``reEstimateBackground=True``. One per band + for a ``MultibandExposure``. backgroundToPhotometricRatio: - Unused; retained for compatibility with the base task. + Image to convert photometric-flattened image to + background-flattened image if ``reEstimateBackground=True`` and + exposure has been photometric-flattened. Returns ------- @@ -243,12 +242,14 @@ def run( def detectFootprints( self, exposure: afwImage.Exposure | afwImage.MultibandExposure, - doSmooth: bool = False, + doSmooth: bool = True, sigma: float | None = None, clearMask: bool = True, expId: int | None = None, - background: afwMath.BackgroundList | None = None, + background: afwMath.BackgroundList | list[afwMath.BackgroundList] | None = None, backgroundToPhotometricRatio: afwImage.Image | None = None, + factor: float = 1.0, + factorNeg: float | None = None, ) -> pipeBase.Struct: """Detect footprints on an exposure. @@ -259,7 +260,10 @@ def detectFootprints( set in-place. May be a single-band ``Exposure`` or a multi-band ``MultibandExposure``. doSmooth: - Unused; retained for compatibility with the base task. + If True, smooth the image before detection using a Gaussian + of width ``sigma``, or the measured PSF width of ``exposure``. + Set to False when running on e.g. a pre-convolved image, or a mask + plane. sigma: Gaussian sigma of the PSF (pixels), used to grow detections; if `None` it is measured from the PSF of the ``exposure``. @@ -269,40 +273,66 @@ def detectFootprints( expId: Unused; retained for compatibility with the base task. background: - Unused; retained for compatibility with the base task. + Background that was already subtracted from the exposure; will be + modified in-place if ``reEstimateBackground=True``. backgroundToPhotometricRatio: - Unused; retained for compatibility with the base task. + Image to convert photometric-flattened image to + background-flattened image if ``reEstimateBackground=True`` and + exposure has been photometric-flattened. + factor: + Multiplier for the configured threshold for positive detection + polarity. + factorNeg: + Multiplier for the configured threshold for negative detection + polarity. If `None`, will be set equal to ``factor`` (i.e. equal + to the factor used for positive detection polarity). Returns ------- results : `lsst.pipe.base.Struct` - See `scarletToAfwResults` for the contents. In addition + See `_mergePeaksAndFootprints` for the contents. In addition `finalizeFootprints` adds the ``numPos``, ``numPosPeaks``, ``numNeg`` and ``numNegPeaks`` counts. """ - multiband = hasattr(exposure, "bands") + isMultiband = hasattr(exposure, "bands") x0, y0 = exposure.getBBox().getMin() - if multiband: - singles = list(exposure.singles) - if clearMask: - for single in singles: - self.clearMask(single.mask) - self.removeBadPixels(single.maskedImage) - fwhm = self._getFwhm(singles, sigma) - images = exposure.image.array - variance = exposure.variance.array - mask = singles[0].mask + # The masked images are views, so setting their mask planes sets them + # on the caller's exposure. Everything that modifies pixels below works + # on a copy instead. + if isMultiband: + singles = list(cast(afwImage.MultibandExposure, exposure).singles) else: - maskedImage = exposure.maskedImage - if clearMask: + singles = [cast(afwImage.Exposure, exposure)] + maskedImages = [single.maskedImage for single in singles] + + if clearMask: + for maskedImage in maskedImages: self.clearMask(maskedImage.getMask()) - self.removeBadPixels(maskedImage) - fwhm = self._getFwhm([exposure], sigma) - images = exposure.image.array[None] + + psfs = self._getDetectionPsfs(singles, sigma) + # The widest band sets the growth radius and the radius that links + # peaks across scales, which is the most conservative choice. + sigma = max(psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() for psf in psfs) + fwhm = sigma*SIGMA_TO_FWHM + + # Instead of calling ``SourceDetectionTask.removeBadPixels`` + # individually for each band, here it is applied to every band at + # once. Only the image is copied, since it is the only plane the peak + # detection modifies. + if isMultiband: + images = exposure.image.array.copy() + variance = exposure.variance.array + masks = exposure.mask.array + else: + images = exposure.image.array[None].copy() variance = exposure.variance.array[None] - mask = maskedImage.getMask() + masks = exposure.mask.array[None] + badPixelMask = maskedImages[0].mask.getPlaneBitMask(self.config.excludeMaskPlanes) + if badPixelMask: + images[(masks & badPixelMask) > 0] = 0 + # Find peaks using starlets as compensated filters positive = None negative = None if self.config.thresholdPolarity in ("positive", "both"): @@ -310,55 +340,82 @@ def detectFootprints( if self.config.thresholdPolarity in ("negative", "both"): negative = self._detectPeaks(-images, variance, (y0, x0), fwhm) - results = self.scarletToAfwResults(exposure, positive, negative) + # Find footprints on the image itself, where extended sources keep the + # flux that the starlet coefficients suppress. + footprintResults = self._detectChi2Footprints( + maskedImages, psfs, doSmooth, factor, factorNeg, background + ) + + results = self._mergePeaksAndFootprints(exposure, footprintResults, positive, negative) # Grow the footprints and set the DETECTED mask planes. The footprints # live in a single coordinate system, so growing once is enough; the # mask plane is then set on every band. - sigmaToGrow = fwhm/SIGMA_TO_FWHM - self.finalizeFootprints(mask, results, sigmaToGrow, factor=1.0, factorNeg=1.0) - if multiband: - for single in singles[1:]: - for polarity, maskName in (("positive", "DETECTED"), ("negative", "DETECTED_NEGATIVE")): - fpSet = getattr(results, polarity) - if fpSet is not None: - fpSet.setMask(single.mask, maskName) - self.clearUnwantedResults(single.mask, results) - - self.clearUnwantedResults(mask, results) + self.finalizeFootprints(maskedImages[0].mask, results, sigma, factor=factor, factorNeg=factorNeg) + for maskedImage in maskedImages[1:]: + for polarity, maskName in (("positive", "DETECTED"), ("negative", "DETECTED_NEGATIVE")): + fpSet = getattr(results, polarity) + if fpSet is not None: + fpSet.setMask(maskedImage.mask, maskName) + + if self.config.reEstimateBackground: + for maskedImage, backgrounds in zip(maskedImages, results.background): + self.reEstimateBackground( + maskedImage, + backgrounds, + backgroundToPhotometricRatio=backgroundToPhotometricRatio, + ) + + # Drop the unwanted polarity only once every band has had its mask + # plane set, since the first call clears that polarity from ``results``. + for maskedImage in maskedImages: + self.clearUnwantedResults(maskedImage.getMask(), results) + + if not isMultiband: + results.background = results.background[0] return results - def _getFwhm(self, singles: list[afwImage.Exposure], sigma: float | None) -> float: - """Measure the PSF FWHM to use for growing and linking detections. + def _getDetectionPsfs( + self, + singles: list[afwImage.Exposure], + sigma: float | None, + ) -> list[afwDet.GaussianPsf]: + """Build the Gaussian PSF to convolve each band with. Parameters ---------- singles: - The single-band exposures to measure the PSF from. For a multi-band - exposure the widest PSF is used, since it sets the most conservative - linking radius across bands. + The single-band exposures to measure the PSF from. sigma: Gaussian sigma of the PSF (pixels); if `None` it is measured from - each exposure's PSF. + each exposure's own PSF. Returns ------- - fwhm: - The PSF full width at half maximum, in pixels. + psfs: + One PSF per band, in band order. + + Notes + ----- + The chi^2 coadd weights each band by the variance of its own matched + filter, so each band is correlated with its own PSF. + See Lupton, chi^2 coadds, equation 1.16. """ - if sigma is not None: - return sigma*SIGMA_TO_FWHM - sigmas = [] + psfs = [] for single in singles: - try: - psf = self.getPsf(single, sigma=None) - sigmas.append(psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()) - except Exception: - continue - if not sigmas: - return SIGMA_TO_FWHM - return max(max(sigmas)*SIGMA_TO_FWHM, SIGMA_TO_FWHM) + bandSigma = sigma + if bandSigma is None: + psf = single.getPsf() + if psf is not None: + bandSigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() + if bandSigma is None or not np.isfinite(bandSigma): + bandSigma = 1.0 + # A Gaussian narrower than a pixel is not representable on the + # pixel grid. + bandSigma = max(bandSigma, 1.0) + psfs.append(self.getPsf(single, sigma=bandSigma)) + return psfs def _detectPeaks( self, @@ -395,14 +452,291 @@ def _detectPeaks( min_separation=0, min_area=self.config.minPixels, peak_thresh=self.config.peakThreshold, - footprint_thresh=self.config.thresholdValue, + footprint_thresh=self.config.filterThreshold, psf_fwhm=fwhm, kappa=self.config.saddleThreshold, + variance_mode=VARIANCE_MODE[self.config.thresholdType], ) - def scarletToAfwResults( + def _bandNoise(self, convolvedImage: afwImage.MaskedImage) -> np.ndarray | float: + """Measure the noise of one convolved band. + + Parameters + ---------- + convolvedImage: + The convolved masked image of a single band. + + Returns + ------- + noise: + The standard deviation to divide the band by, per pixel for + ``thresholdType="pixel_stdev"`` and a single value for + ``thresholdType="stdev"``. + """ + if self.config.thresholdType == "pixel_stdev": + return np.sqrt(convolvedImage.variance.array) + # The median has to skip the pixels with no usable variance, which + # would otherwise drag it to zero over a region with no data. + statsControl = afwMath.StatisticsControl() + statsControl.setAndMask(convolvedImage.mask.getPlaneBitMask(self.config.statsMask)) + statistics = afwMath.makeStatistics(convolvedImage.variance, convolvedImage.mask, + afwMath.MEDIAN, statsControl) + return np.sqrt(statistics.getValue(afwMath.MEDIAN)) + + def _standardizeBands(self, convolvedImages: list[afwImage.MaskedImage]) -> np.ndarray: + """Divide each convolved band by the noise of its own matched filter. + + Parameters + ---------- + convolvedImages: + The convolved masked images, one per band, all sharing a bounding + box. + + Returns + ------- + standardized: + The per-band planes, with shape ``(nBands, Ny, Nx)``, each in + units of its own noise. + """ + standardized = np.zeros((len(convolvedImages),) + convolvedImages[0].image.array.shape, + dtype=np.float32) + for index, convolvedImage in enumerate(convolvedImages): + # A pixel with no variance carries no information, and dividing + # by it would spread NaN across the whole coadd. + noise = self._bandNoise(convolvedImage) + np.divide(convolvedImage.image.array, noise, out=standardized[index], where=noise > 0) + return standardized + + @staticmethod + def _assembleChi2Image( + convolvedImages: list[afwImage.MaskedImage], + standardized: np.ndarray, + ) -> afwImage.MaskedImage: + """Coadd the standardized bands into a chi^2 detection image. + + Parameters + ---------- + convolvedImages: + The convolved masked images the standardized bands came from, which + supply the bounding box and the mask planes. + standardized: + The per-band planes from `_standardizeBands`, negated to detect + negative sources. + + Returns + ------- + chi2Image: + The coadd, in units of Gaussian sigma, with the union of the + per-band masks. + + Notes + ----- + The variance plane of the coadd is left empty. Clipping the per-band + planes at zero means the chi^2 coadd has no meaningful per-pixel + variance; only the upper tail is calibrated, which is what + `build_chi2_significance` maps to sigma and is all that a positive + threshold uses. + """ + chi2Image = afwImage.MaskedImageF(convolvedImages[0].getBBox()) + for convolvedImage in convolvedImages: + chi2Image.mask.array |= convolvedImage.mask.array + chi2Image.image.array[:] = scl.detect.build_chi2_significance(standardized) + return chi2Image + + def makeThreshold( + self, + image: afwImage.MaskedImage, + thresholdParity: str, + factor: float = 1.0, + ) -> afwDet.Threshold: + """Make the threshold to apply to a chi^2 detection image. + + Parameters + ---------- + image: + Unused; the threshold does not depend on the image. + thresholdParity: + One of "positive" or "negative", to set the kind of fluctuations + the threshold will detect. + factor: + Factor by which to multiply the configured detection threshold. + + Returns + ------- + threshold: + Detection threshold. + + Notes + ----- + The bands are divided by their noise in `_standardizeBands`, so by the + time the chi^2 coadd exists it is already in the units that + ``config.thresholdType`` names and the threshold is applied to it + directly. The base task instead divides here, because it thresholds an + image in counts. + """ + threshold = afwDet.createThreshold(self.config.thresholdValue*factor, "value", + thresholdParity != "negative") + threshold.setIncludeMultiplier(self.config.includeThresholdMultiplier) + self.log.debug("Detection threshold: %s", threshold) + return threshold + + def applyThreshold( + self, + middle: afwImage.MaskedImage | list[afwImage.MaskedImage], + bbox: geom.Box2I, + factor: float = 1.0, + factorNeg: float | None = None, + ) -> pipeBase.Struct: + r"""Threshold the chi^2 coadd of the convolved bands. + + Parameters + ---------- + middle: + The convolved masked images to coadd, one per band, all sharing a + bounding box. A single masked image is treated as one band. + bbox: + Bounding box of the unconvolved image, which the returned + `~lsst.afw.detection.FootprintSet`\ s take as their region. + factor: + Multiplier for the configured threshold. + factorNeg: + Multiplier for the configured threshold for negative detection + polarity. If `None`, will be set equal to ``factor``. + + Returns + ------- + results : `lsst.pipe.base.Struct` + A `~lsst.pipe.base.Struct` containing: + + ``positive`` + Positive polarity footprints. + (`lsst.afw.detection.FootprintSet` or `None`) + ``negative`` + Negative polarity footprints. + (`lsst.afw.detection.FootprintSet` or `None`) + ``factor`` + Multiplication factor applied to the configured threshold + for positive detection polarity. (`float`) + ``factorNeg`` + Multiplication factor applied to the configured threshold + for negative detection polarity. (`float`) + ``negativeThreshold`` + Threshold used for negative detection polarity. (`float`) + ``positiveImage`` + Coadded image used for positive detection polarity. + (`lsst.afw.image.MaskedImage` or `None`) + ``negativeImage`` + Coadded image used for negative detection polarity. + (`lsst.afw.image.MaskedImage` or `None`) + + Notes + ----- + A chi^2 coadd is clipped at zero, so it has no negative side and each + polarity needs a coadd of its own: the negative detections are the + positive detections of the coadd built from the negated bands. Both + are therefore thresholded on their positive side. + """ + if factorNeg is None: + factorNeg = factor + self.log.info("Threshold scaling factor for positive detections is: %.3f. For negative " + "detections it is: %.3f", factor, factorNeg) + convolvedImages = list(middle) if isinstance(middle, (list, tuple)) else [middle] + results = pipeBase.Struct( + positive=None, + negative=None, + factor=factor, + factorNeg=factorNeg, + positiveThreshold=None, + negativeThreshold=None, + positiveImage=None, + negativeImage=None, + ) + + # Both polarities coadd the same standardized bands, so the negative + # pass negates them in place rather than building them again. + standardized = self._standardizeBands(convolvedImages) + for polarity, maskName, polarityFactor in (("positive", "DETECTED", factor), + ("negative", "DETECTED_NEGATIVE", factorNeg)): + if (not self.config.reEstimateBackground + and self.config.thresholdPolarity not in (polarity, "both")): + continue + if polarity == "negative": + np.negative(standardized, out=standardized) + image = self._assembleChi2Image(convolvedImages, standardized) + threshold = self.makeThreshold(image, "positive", factor=polarityFactor) + footprints = afwDet.FootprintSet(image, threshold, maskName, self.config.minPixels) + footprints.setRegion(bbox) + setattr(results, polarity, footprints) + setattr(results, polarity + "Threshold", threshold) + setattr(results, polarity + "Image", image) + return results + + def _detectChi2Footprints( + self, + maskedImages: list[afwImage.MaskedImage], + psfs: list[afwDet.GaussianPsf], + doSmooth: bool = True, + factor: float = 1.0, + factorNeg: float | None = None, + background: afwMath.BackgroundList | list[afwMath.BackgroundList] | None = None, + ) -> pipeBase.Struct: + r"""Convolve each band and threshold their chi^2 coadd. + + Parameters + ---------- + maskedImages: + The masked images to detect footprints on, one per band. + psfs: + The PSF to convolve each band with, in band order. + doSmooth: + Whether to convolve before detection. + factor: + The multiplication factor applied to the positive threshold. + factorNeg: + The multiplication factor applied to the negative threshold. + background: + Background that was already subtracted from the exposure; will be + modified in-place if ``reEstimateBackground=True``. One per band, + or a single one for a single-band exposure. + + Returns + ------- + results: `lsst.pipe.base.Struct` + What `applyThreshold` returns, plus: + + ``background`` + One `lsst.afw.math.BackgroundList` per band, to be filled by + `reEstimateBackground`. (`list`) + """ + bbox = maskedImages[0].getBBox() + convolvedImages = [] + for maskedImage, psf in zip(maskedImages, psfs): + middle = self.convolveImage(maskedImage, psf, doSmooth=doSmooth).middle + self.removeBadPixels(middle) + convolvedImages.append(middle) + + # A wider PSF loses a wider border to the convolution, so the bands can + # come back different sizes. The coadd is defined where all of them are. + common = geom.Box2I(convolvedImages[0].getBBox()) + for middle in convolvedImages[1:]: + common.clip(middle.getBBox()) + if any(middle.getBBox() != common for middle in convolvedImages): + convolvedImages = [middle.Factory(middle, common, afwImage.PARENT, False) + for middle in convolvedImages] + + results = self.applyThreshold(convolvedImages, bbox, factor=factor, factorNeg=factorNeg) + + if background is None: + background = [afwMath.BackgroundList() for _ in maskedImages] + elif isinstance(background, afwMath.BackgroundList): + background = [background] + results.background = background + return results + + def _mergePeaksAndFootprints( self, exposure: afwImage.Exposure | afwImage.MultibandExposure, + footprintResults: pipeBase.Struct, positive: scl.detect.PeakDetectionResult | None = None, negative: scl.detect.PeakDetectionResult | None = None, ) -> pipeBase.Struct: @@ -412,6 +746,9 @@ def scarletToAfwResults( ---------- exposure: The exposure detection was run on, used for the footprint region. + footprintResults: + The footprints detected on the chi^2 image, from + `_detectChi2Footprints`. positive: The positive polarity detection result. negative: @@ -429,16 +766,19 @@ def scarletToAfwResults( Negative polarity footprints. (`lsst.afw.detection.FootprintSet` or `None`) ``factor`` - Multiplication factor applied to the positive threshold; always - 1.0. (`float`) + Multiplication factor applied to the positive threshold. + (`float`) ``factorNeg`` - Multiplication factor applied to the negative threshold; always - 1.0. (`float`) + Multiplication factor applied to the negative threshold. + (`float`) ``background`` - Always `None`; this task does not estimate a background. + One `lsst.afw.math.BackgroundList` per band, to be filled by + `reEstimateBackground`. (`list`) ``peaks`` The final detections, one row per source, with a ``polarity`` - column. (`astropy.table.Table`) + column and a ``footprint`` column indexing the footprint the + peak was placed in, ``-1`` for a peak that fell outside every + footprint. (`astropy.table.Table`) ``candidates`` The peak candidates, before grouping, with a ``polarity`` column. (`astropy.table.Table`) @@ -457,12 +797,18 @@ def scarletToAfwResults( ``significanceMap``. (`numpy.ndarray` or `None`) """ bbox = exposure.getBBox() + # The peak tables record which footprint each peak landed in, so the + # footprint sets have to be built before the tables are joined. + positiveFootprints = self._insertPeaks(positive, footprintResults.positive, + footprintResults.positiveImage, bbox) + negativeFootprints = self._insertPeaks(negative, footprintResults.negative, + footprintResults.negativeImage, bbox) results = pipeBase.Struct( - positive=self._makeFootprintSet(positive, bbox), - negative=self._makeFootprintSet(negative, bbox), - factor=1.0, - factorNeg=1.0, - background=None, + positive=positiveFootprints, + negative=negativeFootprints, + factor=footprintResults.factor, + factorNeg=footprintResults.factorNeg, + background=footprintResults.background, peaks=self._joinPolarities(positive, negative, "peaks"), candidates=self._joinPolarities(positive, negative, "candidates"), positions=self._joinPolarities(positive, negative, "positions"), @@ -478,30 +824,104 @@ def scarletToAfwResults( results.sigma = None return results - @staticmethod - def _makeFootprintSet( + def _insertPeaks( + self, result: scl.detect.PeakDetectionResult | None, + footprintSet: afwDet.FootprintSet | None, + detectionImage: afwImage.MaskedImage | None, bbox: geom.Box2I, ) -> afwDet.FootprintSet | None: - """Build an afw FootprintSet from a scarlet detection result. + """Place the starlet peaks into the footprints that contain them. Parameters ---------- result: - The detection result to convert. + The scarlet detection result holding the peaks, or `None` if this + polarity was not detected. + footprintSet: + The footprints detected on the chi^2 image for this polarity. Their + own peaks, found by `lsst.afw.detection.FootprintSet` on the way + in, are discarded in favor of the starlet peaks. + detectionImage: + The chi^2 image the footprints were detected on, which supplies + each peak's ``peakValue``. bbox: The region the footprints were detected in. Returns ------- footprintSet: - The footprints, or `None` if ``result`` is `None`. + The footprints that hold at least one peak, brightest peak first, + or `None` if any input is `None`. + + Notes + ----- + A peak's ``peakValue`` is read off the detection image, as + `lsst.afw.detection.FootprintSet` does, rather than carried over from + the starlet significance map the peak was found in. The two are + different images, and the peak ordering within a footprint has to + follow the image the footprint came from. The starlet significance + stays available as ``peak_sigma`` in the output peak table. + + ``result.peaks["footprint"]`` is overwritten with the index of the + footprint each peak was placed in, and set to ``-1`` for a peak that + fell outside every footprint. Those peaks reach the output table but no + footprint, so a detection is never silently dropped. + + A footprint with no peak in it is dropped. The chi^2 threshold sits + well below the peak threshold, so such footprints are mostly noise, and + a parent with no peak is of no use to the deblender in any case. """ - if result is None: + if result is None or footprintSet is None or detectionImage is None: return None - footprintSet = afwDet.FootprintSet(bbox) - footprintSet.setFootprints([scarletFootprintToAfw(fp) for fp in result.footprints]) - return footprintSet + + peaks = result.peaks + footprints = footprintSet.getFootprints() + merged = afwDet.FootprintSet(bbox) + if len(peaks) == 0 or len(footprints) == 0: + peaks["footprint"] = -1 + return merged + + # Label every footprint's pixels with its index, so that each peak can + # be looked up in one pass rather than tested against every footprint. + # The peaks come from the full exposure, but the detection image loses + # a border to the convolution, so they are not all inside it. + imageBBox = detectionImage.getBBox() + x0, y0 = imageBBox.getMin() + width, height = imageBBox.getDimensions() + column = peaks["x"] - x0 + row = peaks["y"] - y0 + inside = (column >= 0) & (column < width) & (row >= 0) & (row < height) + + labels = afwImage.ImageI(imageBBox) + for index, footprint in enumerate(footprints): + footprint.spans.setImage(labels, index + 1, doClip=True) + labelOf = np.full(len(peaks), -1, dtype=int) + labelOf[inside] = labels.array[row[inside], column[inside]] - 1 + peakValue = np.zeros(len(peaks), dtype=float) + peakValue[inside] = detectionImage.image.array[row[inside], column[inside]] + + # Brightest first within a footprint, matching the afw convention. + order = np.lexsort((-peakValue, labelOf)) + kept: list[afwDet.Footprint] = [] + footprintOf = np.full(len(peaks), -1, dtype=peaks["footprint"].dtype) + current = -1 + for index in order: + label = labelOf[index] + if label < 0: + continue + if label != current: + footprints[label].getPeaks().clear() + kept.append(footprints[label]) + current = label + footprintOf[index] = len(kept) - 1 + kept[-1].addPeak(int(peaks["x"][index]), int(peaks["y"][index]), float(peakValue[index])) + peaks["footprint"] = footprintOf + + merged.setFootprints(kept) + self.log.info("Merged %d of %d peaks into %d of %d chi^2 footprints", + int(np.sum(footprintOf >= 0)), len(peaks), len(kept), len(footprints)) + return merged @staticmethod def _joinPolarities( @@ -546,21 +966,9 @@ def _joinPolarities( def applyTempLocalBackground(self, exposure, middle, results): raise NotImplementedError("applyTempLocalBackground is not used in MultiResolutionDetectionTask") - def convolveImage(self, maskedImage, psf, doSmooth=True): - raise NotImplementedError("convolveImage is not used in MultiResolutionDetectionTask") - - def applyThreshold(self, middle, bbox, factor=1.0, factorNeg=None): - raise NotImplementedError("applyThreshold is not used in MultiResolutionDetectionTask") - - def reEstimateBackground(self, maskedImage, backgrounds, backgroundToPhotometricRatio=None): - raise NotImplementedError("reEstimateBackground is not used in MultiResolutionDetectionTask") - def setPeakSignificance(self, exposure, footprints, threshold, negative=False): raise NotImplementedError("setPeakSignificance is not used in MultiResolutionDetectionTask") - def makeThreshold(self, image, thresholdParity, factor=1.0): - raise NotImplementedError("makeThreshold is not used in MultiResolutionDetectionTask") - def updatePeaks(self, fpSet, image, threshold): raise NotImplementedError("updatePeaks is not used in MultiResolutionDetectionTask") diff --git a/tests/test_multiResolutionDetection.py b/tests/test_multiResolutionDetection.py index d90be86a..0eb16c30 100644 --- a/tests/test_multiResolutionDetection.py +++ b/tests/test_multiResolutionDetection.py @@ -20,16 +20,22 @@ # along with this program. If not, see . import unittest +import unittest.mock +import numpy as np from astropy.table import Table import lsst.afw.detection as afwDet import lsst.afw.image as afwImage +import lsst.afw.math as afwMath import lsst.afw.table as afwTable import lsst.geom import lsst.pex.config +import lsst.scarlet.lite as scl import lsst.utils.tests from lsst.meas.algorithms import MultiResolutionDetectionConfig, MultiResolutionDetectionTask +from lsst.meas.algorithms import multiResolutionDetection +from lsst.meas.algorithms.multiResolutionDetection import SIGMA_TO_FWHM from lsst.meas.algorithms.testUtils import plantSources @@ -49,16 +55,269 @@ def _makeExposure(self, offset: int = 0) -> afwImage.Exposure: coords = [(x, y, counts + offset, sigma) for x, y, counts, sigma in self.coords] return plantSources(self.bbox, self.kwid, self.sky, coords, addPoissonNoise=True) + def _makeExtendedExposure(self) -> afwImage.Exposure: + """Plant one compact and one broad source, at a non-zero origin.""" + bbox = lsst.geom.Box2I(lsst.geom.Point2I(1000, 2000), lsst.geom.Extent2I(180, 180)) + coords = [(1040, 2040, 6000, self.starSigma), + (1110, 2110, 60000, 9.0)] + return plantSources(bbox, 15, self.sky, coords, addPoissonNoise=True) + def testConfigValidation(self) -> None: - """The config must default to no background re-estimation and reject - turning it back on.""" + """The config must reject the background steps this task does not + implement, and allow the base task's background re-estimation.""" config = MultiResolutionDetectionConfig() - self.assertFalse(config.reEstimateBackground) self.assertFalse(config.doTempLocalBackground) self.assertFalse(config.doTempWideBackground) config.reEstimateBackground = True - with self.assertRaises(lsst.pex.config.FieldValidationError): + config.validate() + for name in ("doTempLocalBackground", "doTempWideBackground"): + config = MultiResolutionDetectionConfig() + setattr(config, name, True) + with self.assertRaises(lsst.pex.config.FieldValidationError): + config.validate() + + def testExtendedSourceFootprints(self) -> None: + """Footprints must follow the extended flux that the starlet + coefficients suppress, and stay in the exposure's coordinate system.""" + exposure = self._makeExtendedExposure() + config = MultiResolutionDetectionConfig() + config.reEstimateBackground = False + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + results = task.run(afwTable.SourceTable.make(schema), exposure) + + self.assertEqual(results.positive.getRegion(), exposure.getBBox()) + footprints = results.positive.getFootprints() + self.assertEqual(len(footprints), 2) + for footprint in footprints: + self.assertTrue(exposure.getBBox().contains(footprint.getBBox())) + + starletResult = scl.detect.detect_peaks( + images=exposure.image.array[None], + variance=exposure.variance.array[None], + scales=config.starletScales, + generation=config.generation, + first_scale=config.firstScale, + origin=(exposure.getY0(), exposure.getX0()), + min_separation=0, + min_area=config.minPixels, + peak_thresh=config.peakThreshold, + footprint_thresh=config.filterThreshold, + psf_fwhm=self.starSigma*SIGMA_TO_FWHM, + kappa=config.saddleThreshold, + ) + # The starlet coefficient footprints stop inside the source, so the + # broadest of them is still smaller than the smallest chi^2 footprint. + widestStarlet = max(int(fp.data.sum()) for fp in starletResult.footprints) + self.assertGreater(min(fp.getArea() for fp in footprints), widestStarlet) + + def testPeaksAssignedToFootprints(self) -> None: + """Every footprint must hold at least one peak, and the peak table must + index the footprint each peak was placed in.""" + exposure = self._makeExtendedExposure() + config = MultiResolutionDetectionConfig() + config.reEstimateBackground = False + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + results = task.run(afwTable.SourceTable.make(schema), exposure) + + footprints = results.positive.getFootprints() + for footprint in footprints: + self.assertGreater(len(footprint.getPeaks()), 0) + peakValues = [peak["peakValue"] for peak in footprint.getPeaks()] + self.assertEqual(peakValues, sorted(peakValues, reverse=True)) + + self.assertIn("footprint", results.peaks.colnames) + for row in results.peaks: + if row["footprint"] < 0: + continue + footprint = footprints[row["footprint"]] + self.assertTrue(footprint.contains(lsst.geom.Point2I(row["x"], row["y"]))) + self.assertEqual(sum(len(fp.getPeaks()) for fp in footprints), + int((results.peaks["footprint"] >= 0).sum())) + + def testPerBandPsfWidths(self) -> None: + """Bands with different PSF widths are each correlated with their own + PSF, so the convolution trims them by different amounts.""" + bands = ["g", "r", "i"] + psfSigmas = [1.5, 2.5, 4.0] + exposures = [] + for bandSigma in psfSigmas: + coords = [(x, y, counts, bandSigma) for x, y, counts, _ in self.coords] + exposure = plantSources(self.bbox, self.kwid, self.sky, coords, addPoissonNoise=True) + exposure.setPsf(afwDet.GaussianPsf(self.kwid, self.kwid, bandSigma)) + exposures.append(exposure) + mExposure = afwImage.MultibandExposure.fromExposures(bands, exposures) + + config = MultiResolutionDetectionConfig() + config.reEstimateBackground = False + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + + singles = list(mExposure.singles) + psfs = task._getDetectionPsfs(singles, None) + radii = [psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() for psf in psfs] + self.assertFloatsAlmostEqual(np.array(radii), np.array(psfSigmas), atol=1e-6) + self.assertGreater(len(set(psf.getDimensions().getX() for psf in psfs)), 1) + + results = task.run(afwTable.SourceTable.make(schema), mExposure) + self.assertGreater(len(results.positive.getFootprints()), 0) + for footprint in results.positive.getFootprints(): + self.assertGreater(len(footprint.getPeaks()), 0) + + def testPeakValueFromDetectionImage(self) -> None: + """A peak's value must be read off the image its footprint came from, + not carried over from the starlet significance map.""" + exposure = self._makeExtendedExposure() + config = MultiResolutionDetectionConfig() + config.reEstimateBackground = False + config.nSigmaToGrow = 0 + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + results = task.run(afwTable.SourceTable.make(schema), exposure) + + peakRows = {(row["x"], row["y"]): row for row in results.peaks if row["footprint"] >= 0} + self.assertGreater(len(peakRows), 0) + differs = False + for footprint in results.positive.getFootprints(): + for peak in footprint.getPeaks(): + row = peakRows[(peak.getIx(), peak.getIy())] + # Every peak is above the footprint threshold on that image. + self.assertGreater(peak["peakValue"], config.thresholdValue) + if abs(peak["peakValue"] - row["peak_sigma"]) > 1e-4: + differs = True + self.assertTrue(differs, "peakValue should not simply mirror peak_sigma") + + def testNoRedundantThresholding(self) -> None: + """Each polarity is thresholded once, on the coadd built for it, even + when background re-estimation asks for both.""" + exposure = self._makeExposure() + cases = (("positive", True, 2), ("positive", False, 1), + ("both", True, 2), ("both", False, 2), + ("negative", True, 2), ("negative", False, 1)) + for polarity, reEstimateBackground, expected in cases: + config = MultiResolutionDetectionConfig() + config.thresholdPolarity = polarity + config.reEstimateBackground = reEstimateBackground + # Growing also builds a FootprintSet, which would be counted below. + config.nSigmaToGrow = 0 + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + + original = multiResolutionDetection.afwDet.FootprintSet + scans = [] + + def counting(*args, **kwargs): + # The thresholding constructor takes an image and a threshold; + # FootprintSet(bbox) takes neither. + if len(args) > 1: + scans.append(args[1]) + return original(*args, **kwargs) + + with unittest.mock.patch.object(multiResolutionDetection.afwDet, "FootprintSet", counting): + task.run(afwTable.SourceTable.make(schema), exposure.clone()) + self.assertEqual(len(scans), expected, + f"thresholdPolarity={polarity}, " + f"reEstimateBackground={reEstimateBackground}") + + def testThresholdTypeValidation(self) -> None: + """thresholdType selects the noise the bands are divided by, so only + the two standard deviation choices mean anything here.""" + self.assertEqual(MultiResolutionDetectionConfig().thresholdType, "stdev") + for thresholdType in ("pixel_stdev", "stdev"): + config = MultiResolutionDetectionConfig() + config.thresholdType = thresholdType config.validate() + for thresholdType in ("value", "variance"): + config = MultiResolutionDetectionConfig() + config.thresholdType = thresholdType + with self.assertRaises(lsst.pex.config.FieldValidationError): + config.validate() + + def testBandNoise(self) -> None: + """The per-pixel noise follows the variance plane, and the single + value skips the pixels with no usable variance.""" + convolved = afwImage.MaskedImageF(self.bbox) + convolved.variance.array[:] = 4.0 + convolved.variance.array[:40] = 0.0 + convolved.mask.array[:40] |= convolved.mask.getPlaneBitMask("NO_DATA") + schema = afwTable.SourceTable.makeMinimalSchema() + + config = MultiResolutionDetectionConfig() + config.thresholdType = "pixel_stdev" + task = MultiResolutionDetectionTask(schema=schema, config=config) + perPixel = task._bandNoise(convolved) + self.assertEqual(np.shape(perPixel), convolved.image.array.shape) + self.assertFloatsAlmostEqual(perPixel[40:], 2.0) + + config = MultiResolutionDetectionConfig() + config.thresholdType = "stdev" + task = MultiResolutionDetectionTask(schema=schema, config=config) + single = task._bandNoise(convolved) + self.assertEqual(np.shape(single), ()) + # A raw median over the whole plane would be dragged to zero. + self.assertFloatsAlmostEqual(float(single), 2.0) + + def testThresholdTypeReachesPeakDetection(self) -> None: + """The starlet coefficients are standardized by the same noise as the + footprints, so the choice reaches scarlet lite too.""" + exposure = self._makeExtendedExposure() + expected = {"pixel_stdev": "pixel", "stdev": "median"} + for thresholdType, varianceMode in expected.items(): + config = MultiResolutionDetectionConfig() + config.thresholdType = thresholdType + config.reEstimateBackground = False + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + seen = [] + original = multiResolutionDetection.scl.detect.detect_peaks + + def recording(*args, **kwargs): + seen.append(kwargs["variance_mode"]) + return original(*args, **kwargs) + + with unittest.mock.patch.object(multiResolutionDetection.scl.detect, + "detect_peaks", recording): + results = task.run(afwTable.SourceTable.make(schema), exposure.clone()) + self.assertEqual(seen, [varianceMode]) + # The noise the coefficients were divided by follows the choice. + self.assertEqual(results.sigma.ndim, 4) + perPixel = varianceMode == "pixel" + self.assertEqual(results.sigma.shape[-2:] != (1, 1), perPixel) + + def testThresholdTypeChangesFootprints(self) -> None: + """A depth step across the image is treated differently by the two + noise choices, on the same pixels.""" + base = self._makeExtendedExposure() + base.variance.array[:, :90] *= 0.25 + detected = {} + for thresholdType in ("pixel_stdev", "stdev"): + config = MultiResolutionDetectionConfig() + config.thresholdType = thresholdType + config.reEstimateBackground = False + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + exposure = base.clone() + results = task.run(afwTable.SourceTable.make(schema), exposure) + self.assertGreater(len(results.positive.getFootprints()), 0) + mask = exposure.mask + detected[thresholdType] = int(((mask.array & mask.getPlaneBitMask("DETECTED")) > 0).sum()) + self.assertNotEqual(detected["pixel_stdev"], detected["stdev"]) + + def testImageUnchanged(self) -> None: + """Removing the bad pixels must act on a copy, leaving the exposure's + image plane alone.""" + exposure = self._makeExposure() + exposure.mask.array[10:20, 10:20] |= exposure.mask.getPlaneBitMask("BAD") + original = exposure.image.array.copy() + config = MultiResolutionDetectionConfig() + config.reEstimateBackground = False + config.excludeMaskPlanes = ["BAD"] + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + task.run(afwTable.SourceTable.make(schema), exposure) + + self.assertFloatsEqual(exposure.image.array, original) def testSingleBand(self) -> None: """Detection on a single-band exposure recovers the planted sources.""" @@ -72,7 +331,7 @@ def testSingleBand(self) -> None: self.assertEqual(results.numPos, len(results.positive.getFootprints())) self.assertIsInstance(results.positive, afwDet.FootprintSet) self.assertIsNone(results.negative) - self.assertIsNone(results.background) + self.assertIsInstance(results.background, afwMath.BackgroundList) # The scarlet intermediate products are attached to the result as # astropy tables. @@ -104,6 +363,10 @@ def testMultiBand(self) -> None: self.assertGreaterEqual(len(results.sources), len(self.coords)) self.assertEqual(results.significanceMap.shape[1], len(bands) + 1) self.assertEqual(results.starlets.shape[1], len(bands)) + # One background per band, so that each can be added back to its own. + self.assertEqual(len(results.background), len(bands)) + for background in results.background: + self.assertIsInstance(background, afwMath.BackgroundList) for single in mExposure.singles: mask = single.mask detected = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0 From b4e38b6d723d3a7c46c56ace5073862da864dc5d Mon Sep 17 00:00:00 2001 From: fred3m Date: Mon, 28 Sep 2026 11:34:26 -0700 Subject: [PATCH 3/6] Update config options for MultiResolutionDetectionTask --- .../algorithms/multiResolutionDetection.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/python/lsst/meas/algorithms/multiResolutionDetection.py b/python/lsst/meas/algorithms/multiResolutionDetection.py index 84817551..f6f96fed 100644 --- a/python/lsst/meas/algorithms/multiResolutionDetection.py +++ b/python/lsst/meas/algorithms/multiResolutionDetection.py @@ -67,19 +67,10 @@ class MultiResolutionDetectionConfig(SourceDetectionConfig): doc="The generation of the starlet transform to use (1 or 2).", dtype=int, default=2, ) - minPixels = pexConfig.Field( - doc="Detected footprints with fewer than the specified number of pixels will be ignored.", - dtype=int, default=8, - ) peakThreshold = pexConfig.Field( doc="Threshold for detecting peaks within footprints, in sigma.", dtype=float, default=5.0, ) - thresholdValue = pexConfig.Field( - doc="Threshold for detecting footprints on the chi^2 detection image, in sigma. Whether that " - "is the per-pixel or the per-band sigma is set by thresholdType.", - dtype=float, default=2.0, - ) filterThreshold = pexConfig.Field( doc="Threshold for the starlet coefficient footprints that peaks are detected on, in sigma.", dtype=float, default=2.0, @@ -99,6 +90,15 @@ def setDefaults(self) -> None: self.doTempLocalBackground = False self.doTempWideBackground = False + # Since we're detecting on a chi^2 image, we don't need to grow the + # pixels by as much to get the same effective detection area. + self.nSigmaToGrow = 1.0 + + # Since peak detection is in starlet space, merged across many bands + # and scales, setting a higher minimum number of pixels for a detected + # footprint helps remove spurious detections due to noise fluctuations. + self.minPixels = 8 + # Unlike SourceDetectionTask, we want the default thresholdType to # be 'stdev', since the detection gains on starlet coefficients by # using `pixel_stdev` is minimal at the cost of significantly more From 145d96e106fcb69e776bcabf4ef62ddb3caeedd0 Mon Sep 17 00:00:00 2001 From: fred3m Date: Tue, 29 Sep 2026 12:11:33 -0700 Subject: [PATCH 4/6] Suppress peaks near bad pixels in detection Feed the excludeMaskPlanes mask to the starlet peak finder so that detections ringing around saturated stars are dropped in the bands that flag them, while clean bands still detect the source. The image is still zeroed at those pixels to keep the transform finite; badPixelDilation sets how far the suppression reaches. Also default thresholdType to pixel_stdev. The median stdev underestimates the noise on images with a long tail, where a plain stdev threshold yields more false detections than real ones. Generated with AI Co-Authored-By: SLAC AI --- .../algorithms/multiResolutionDetection.py | 36 +++++++++--- tests/test_multiResolutionDetection.py | 55 ++++++++++++++++++- 2 files changed, 82 insertions(+), 9 deletions(-) diff --git a/python/lsst/meas/algorithms/multiResolutionDetection.py b/python/lsst/meas/algorithms/multiResolutionDetection.py index f6f96fed..67d0016c 100644 --- a/python/lsst/meas/algorithms/multiResolutionDetection.py +++ b/python/lsst/meas/algorithms/multiResolutionDetection.py @@ -79,6 +79,14 @@ class MultiResolutionDetectionConfig(SourceDetectionConfig): doc="Minimum prominence of a peak above the saddle to a brighter peak, in sigma.", dtype=float, default=3.0, ) + badPixelDilation = pexConfig.Field( + doc="Base radius in pixels by which the excludeMaskPlanes mask is grown before it suppresses " + "peaks. Starlet ringing seeds false peaks in a halo around a bad region, so the " + "suppression reaches beyond the flagged pixels. The kernel doubles in width with each " + "starlet scale, so the radius applied at a scale grows with it. The current value is " + "calibrated to DP2. Grow this only if false peaks persist around bright sources.", + dtype=float, default=3.0, + ) def setDefaults(self) -> None: super().setDefaults() @@ -99,11 +107,11 @@ def setDefaults(self) -> None: # footprint helps remove spurious detections due to noise fluctuations. self.minPixels = 8 - # Unlike SourceDetectionTask, we want the default thresholdType to - # be 'stdev', since the detection gains on starlet coefficients by - # using `pixel_stdev` is minimal at the cost of significantly more - # memory usage. - self.thresholdType = "stdev" + # We want to remove peak detections from bands where they have + # a bad pixel mask. This removes a lot of junk around saturated + # stars but allows through detections from bands without masked + # pixels. + self.excludeMaskPlanes = ["BAD", "SAT", "EDGE", "NO_DATA"] def validate(self) -> None: super().validate() @@ -329,16 +337,21 @@ def detectFootprints( variance = exposure.variance.array[None] masks = exposure.mask.array[None] badPixelMask = maskedImages[0].mask.getPlaneBitMask(self.config.excludeMaskPlanes) + badPixels = None if badPixelMask: - images[(masks & badPixelMask) > 0] = 0 + badPixels = (masks & badPixelMask) > 0 + # Zeroing the flagged pixels keeps the starlet transform finite and + # bounds its ringing; badPixels then suppresses the peaks that + # ringing seeds around a bad region. + images[badPixels] = 0 # Find peaks using starlets as compensated filters positive = None negative = None if self.config.thresholdPolarity in ("positive", "both"): - positive = self._detectPeaks(images, variance, (y0, x0), fwhm) + positive = self._detectPeaks(images, variance, (y0, x0), fwhm, badPixels) if self.config.thresholdPolarity in ("negative", "both"): - negative = self._detectPeaks(-images, variance, (y0, x0), fwhm) + negative = self._detectPeaks(-images, variance, (y0, x0), fwhm, badPixels) # Find footprints on the image itself, where extended sources keep the # flux that the starlet coefficients suppress. @@ -423,6 +436,7 @@ def _detectPeaks( variance: np.ndarray, origin: tuple[int, int], fwhm: float, + badPixels: np.ndarray | None = None, ) -> scl.detect.PeakDetectionResult: """Run the scarlet lite peak detection on a set of images. @@ -436,6 +450,10 @@ def _detectPeaks( The ``(y, x)`` location of the lower corner of the images. fwhm: The PSF full width at half maximum, in pixels. + badPixels: + Per-band boolean mask with the same shape as ``images``, `True` + where a band is bad. Peaks near a bad region are suppressed in the + bands that flag it. `None` disables the suppression. Returns ------- @@ -456,6 +474,8 @@ def _detectPeaks( psf_fwhm=fwhm, kappa=self.config.saddleThreshold, variance_mode=VARIANCE_MODE[self.config.thresholdType], + bad_pixel_mask=badPixels, + dilation_radius=self.config.badPixelDilation, ) def _bandNoise(self, convolvedImage: afwImage.MaskedImage) -> np.ndarray | float: diff --git a/tests/test_multiResolutionDetection.py b/tests/test_multiResolutionDetection.py index 0eb16c30..9961aa4c 100644 --- a/tests/test_multiResolutionDetection.py +++ b/tests/test_multiResolutionDetection.py @@ -223,7 +223,7 @@ def counting(*args, **kwargs): def testThresholdTypeValidation(self) -> None: """thresholdType selects the noise the bands are divided by, so only the two standard deviation choices mean anything here.""" - self.assertEqual(MultiResolutionDetectionConfig().thresholdType, "stdev") + self.assertEqual(MultiResolutionDetectionConfig().thresholdType, "pixel_stdev") for thresholdType in ("pixel_stdev", "stdev"): config = MultiResolutionDetectionConfig() config.thresholdType = thresholdType @@ -319,6 +319,59 @@ def testImageUnchanged(self) -> None: self.assertFloatsEqual(exposure.image.array, original) + def testBadPixelSuppressionKeepsCleanBand(self) -> None: + """A source flagged bad in every band but one is still detected, seeded + by the clean band, and the other sources are untouched.""" + bands = ["g", "r", "i"] + exposures = [self._makeExposure(offset=1000 * i) for i in range(len(bands))] + mExposure = afwImage.MultibandExposure.fromExposures(bands, exposures) + + # Flag the first source's core as saturated in every band but g. + x, y = self.coords[0][0], self.coords[0][1] + for single in mExposure.singles[1:]: + single.mask.array[y - 3:y + 4, x - 3:x + 4] |= single.mask.getPlaneBitMask("SAT") + + config = MultiResolutionDetectionConfig() + config.excludeMaskPlanes = ["SAT"] + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + results = task.run(afwTable.SourceTable.make(schema), mExposure) + + peaks = results.peaks + matched = peaks[(np.abs(peaks["y"] - y) <= 2) & (np.abs(peaks["x"] - x) <= 2)] + self.assertEqual(len(matched), 1) + # Only the clean g band seeds it: no other band and no chi coadd plane. + self.assertEqual(int(matched["band_flags"][0]), 1 << 0) + # The other planted sources are still recovered. + for sx, sy, _, _ in self.coords[1:]: + other = peaks[(np.abs(peaks["y"] - sy) <= 2) & (np.abs(peaks["x"] - sx) <= 2)] + self.assertGreaterEqual(len(other), 1) + + def testBadPixelDilationSuppressesPeak(self) -> None: + """Growing badPixelDilation suppresses a peak that a smaller radius, + which does not reach the source, leaves in place.""" + exposure = self._makeExposure() + x, y = self.coords[0][0], self.coords[0][1] + # Flag a pixel a few pixels off the source center in the only band. + exposure.mask.array[y, x - 4] |= exposure.mask.getPlaneBitMask("SAT") + + def countNearFirst(dilation: float) -> int: + single = exposure.clone() + config = MultiResolutionDetectionConfig() + config.excludeMaskPlanes = ["SAT"] + config.badPixelDilation = dilation + schema = afwTable.SourceTable.makeMinimalSchema() + task = MultiResolutionDetectionTask(schema=schema, config=config) + results = task.run(afwTable.SourceTable.make(schema), single) + peaks = results.peaks + near = peaks[(np.abs(peaks["y"] - y) <= 2) & (np.abs(peaks["x"] - x) <= 2)] + return len(near) + + # A radius short of the source keeps the peak; one that reaches its + # center drops it. + self.assertEqual(countNearFirst(1.0), 1) + self.assertEqual(countNearFirst(6.0), 0) + def testSingleBand(self) -> None: """Detection on a single-band exposure recovers the planted sources.""" exposure = self._makeExposure() From a477176535706a931b2afc3a2914b47a146aef07 Mon Sep 17 00:00:00 2001 From: Fred Moolekamp Date: Thu, 8 Oct 2026 18:09:08 -0400 Subject: [PATCH 5/6] Adopt the per-band scarlet peak detection API scarlet_lite's detect_peaks now takes a PSF FWHM per band, derives the minimum separation from the narrowest one, and returns the pairs and components tables that describe how positions were grouped. The tables are joined across polarities like the others, and every table carries the detection settings in its meta. Co-Authored-By: Claude Fable 5.1 --- .../algorithms/multiResolutionDetection.py | 37 +++++++++++++------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/python/lsst/meas/algorithms/multiResolutionDetection.py b/python/lsst/meas/algorithms/multiResolutionDetection.py index 67d0016c..64e8aa29 100644 --- a/python/lsst/meas/algorithms/multiResolutionDetection.py +++ b/python/lsst/meas/algorithms/multiResolutionDetection.py @@ -319,10 +319,11 @@ def detectFootprints( self.clearMask(maskedImage.getMask()) psfs = self._getDetectionPsfs(singles, sigma) - # The widest band sets the growth radius and the radius that links - # peaks across scales, which is the most conservative choice. - sigma = max(psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() for psf in psfs) - fwhm = sigma*SIGMA_TO_FWHM + bandSigma = [psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() for psf in psfs] + fwhm = [s*SIGMA_TO_FWHM for s in bandSigma] + # The widest band sets the growth radius, which is the most + # conservative choice. + sigma = max(bandSigma) # Instead of calling ``SourceDetectionTask.removeBadPixels`` # individually for each band, here it is applied to every band at @@ -435,7 +436,7 @@ def _detectPeaks( images: np.ndarray, variance: np.ndarray, origin: tuple[int, int], - fwhm: float, + fwhm: list[float], badPixels: np.ndarray | None = None, ) -> scl.detect.PeakDetectionResult: """Run the scarlet lite peak detection on a set of images. @@ -449,7 +450,7 @@ def _detectPeaks( origin: The ``(y, x)`` location of the lower corner of the images. fwhm: - The PSF full width at half maximum, in pixels. + The PSF full width at half maximum of each band, in pixels. badPixels: Per-band boolean mask with the same shape as ``images``, `True` where a band is bad. Peaks near a bad region are suppressed in the @@ -467,7 +468,7 @@ def _detectPeaks( generation=self.config.generation, first_scale=self.config.firstScale, origin=origin, - min_separation=0, + min_separation=None, min_area=self.config.minPixels, peak_thresh=self.config.peakThreshold, footprint_thresh=self.config.filterThreshold, @@ -805,6 +806,12 @@ def _mergePeaksAndFootprints( ``positions`` The unique candidate positions, with a ``polarity`` column. (`astropy.table.Table`) + ``pairs`` + The pairs of positions within the largest link radius, with a + ``polarity`` column. (`astropy.table.Table`) + ``components`` + The connected components of the position graph, with a + ``polarity`` column. (`astropy.table.Table`) ``significanceMap`` The per-scale significance map of the positive detection (or the negative detection if only negative was run). (`numpy.ndarray` @@ -832,6 +839,8 @@ def _mergePeaksAndFootprints( peaks=self._joinPolarities(positive, negative, "peaks"), candidates=self._joinPolarities(positive, negative, "candidates"), positions=self._joinPolarities(positive, negative, "positions"), + pairs=self._joinPolarities(positive, negative, "pairs"), + components=self._joinPolarities(positive, negative, "components"), ) primary = positive if positive is not None else negative if primary is not None: @@ -959,14 +968,16 @@ def _joinPolarities( The negative polarity detection result. attr: The name of the structured-array attribute to join (``peaks``, - ``candidates`` or ``positions``). + ``candidates``, ``positions``, ``pairs`` or ``components``). Returns ------- joined: The concatenated rows as an `astropy.table.Table` with an added - ``polarity`` column (``1`` for positive, ``-1`` for negative). The - table is empty if neither result is set. + ``polarity`` column (``1`` for positive, ``-1`` for negative), and + the detection settings in its ``meta``. Row indices such as + ``candidates["position"]`` count within a polarity. The table is + empty if neither result is set. """ arrays = [] for result, polarity in ((positive, 1), (negative, -1)): @@ -979,9 +990,11 @@ def _joinPolarities( arrays.append(array) if not arrays: return Table() + # Both polarities are detected with the same settings. + primary = positive if positive is not None else negative if len(arrays) == 1: - return Table(arrays[0]) - return Table(np.concatenate(arrays)) + return Table(arrays[0], meta=primary.metadata) + return Table(np.concatenate(arrays), meta=primary.metadata) def applyTempLocalBackground(self, exposure, middle, results): raise NotImplementedError("applyTempLocalBackground is not used in MultiResolutionDetectionTask") From 951132a2023aa7c48dd0d2eebb67e4fa8166962d Mon Sep 17 00:00:00 2001 From: Fred Moolekamp Date: Mon, 19 Oct 2026 15:55:55 -0400 Subject: [PATCH 6/6] Expose the peak structure threshold in detection config Co-Authored-By: Claude Fable 5.1 --- .../algorithms/multiResolutionDetection.py | 14 ++++++++++- tests/test_multiResolutionDetection.py | 23 +++++++++++++++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/python/lsst/meas/algorithms/multiResolutionDetection.py b/python/lsst/meas/algorithms/multiResolutionDetection.py index 64e8aa29..bca7186b 100644 --- a/python/lsst/meas/algorithms/multiResolutionDetection.py +++ b/python/lsst/meas/algorithms/multiResolutionDetection.py @@ -87,6 +87,14 @@ class MultiResolutionDetectionConfig(SourceDetectionConfig): "calibrated to DP2. Grow this only if false peaks persist around bright sources.", dtype=float, default=3.0, ) + structureThreshold = pexConfig.Field( + doc="Threshold in sigma that connects the chi^2 coadd of the first starlet scale into " + "structures. The peaks that share a structure, e.g. the arms and star forming regions " + "of a galaxy, get a common structure_id in the peak table and its number of peaks in " + "structure_peaks. Structures of two peaks are mostly blended stars, structures of three " + "or more mostly shredded galaxies.", + dtype=float, default=1.0, + ) def setDefaults(self) -> None: super().setDefaults() @@ -477,6 +485,7 @@ def _detectPeaks( variance_mode=VARIANCE_MODE[self.config.thresholdType], bad_pixel_mask=badPixels, dilation_radius=self.config.badPixelDilation, + structure_thresh=self.config.structureThreshold, ) def _bandNoise(self, convolvedImage: afwImage.MaskedImage) -> np.ndarray | float: @@ -799,7 +808,10 @@ def _mergePeaksAndFootprints( The final detections, one row per source, with a ``polarity`` column and a ``footprint`` column indexing the footprint the peak was placed in, ``-1`` for a peak that fell outside every - footprint. (`astropy.table.Table`) + footprint. The ``structure_id`` and ``structure_peaks`` + columns label peaks that are likely pieces of a larger + structure; see `lsst.scarlet.lite.detect.find_structures`. + (`astropy.table.Table`) ``candidates`` The peak candidates, before grouping, with a ``polarity`` column. (`astropy.table.Table`) diff --git a/tests/test_multiResolutionDetection.py b/tests/test_multiResolutionDetection.py index 9961aa4c..e8dd3167 100644 --- a/tests/test_multiResolutionDetection.py +++ b/tests/test_multiResolutionDetection.py @@ -425,6 +425,29 @@ def testMultiBand(self) -> None: detected = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0 self.assertGreater(detected.sum(), 0) + def testStructureColumns(self) -> None: + """The structure threshold must reach the peak table, and the peak + count must agree with the ids.""" + exposure = self._makeExposure() + schema = afwTable.SourceTable.makeMinimalSchema() + for threshold, expectStructures in ((1.0, True), (1e6, False)): + config = MultiResolutionDetectionConfig() + config.reEstimateBackground = False + config.structureThreshold = threshold + task = MultiResolutionDetectionTask(schema=schema, config=config) + results = task.run(afwTable.SourceTable.make(schema), exposure) + + peaks = results.peaks + self.assertEqual(peaks.meta["structure_thresh"], threshold) + structureId = np.asarray(peaks["structure_id"]) + structurePeaks = np.asarray(peaks["structure_peaks"]) + self.assertEqual(bool((structureId > 0).any()), expectStructures) + self.assertEqual(bool((structurePeaks > 0).any()), expectStructures) + for structure in np.unique(structureId[structureId > 0]): + members = structureId == structure + self.assertTrue((structurePeaks[members] == members.sum()).all()) + self.assertTrue((structurePeaks[structureId == 0] == 0).all()) + def testBothPolarities(self) -> None: """With thresholdPolarity='both', negative footprints are returned and the peak table records both polarities."""