Coverage for python / lsst / meas / algorithms / adaptive_thresholds.py: 11%
123 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 08:51 +0000
1# This file is part of meas_algorithms.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = [
23 "AdaptiveThresholdDetectionConfig",
24 "AdaptiveThresholdDetectionTask",
25]
27import numpy as np
29from lsst.pex.config import Field, Config, DictField, FieldValidationError
30from lsst.pipe.base import Struct, Task
32from .detection import SourceDetectionConfig, SourceDetectionTask
35class AdaptiveThresholdDetectionConfig(Config):
36 """Configuration for AdaptiveThresholdDetectionTask
37 """
38 maxAdaptiveDetIter = Field(dtype=int, default=20,
39 doc="Maximum number of adaptive threshold detection iterations.")
40 maxNumPeakPerBand = DictField(
41 keytype=str,
42 itemtype=int,
43 default={
44 "u": 3000,
45 "g": 3000,
46 "r": 3000,
47 "i": 5000,
48 "z": 5000,
49 "y": 3000,
50 "fallback": 5000,
51 },
52 doc="Maximum number of peaks per band. If the current band for the exposure "
53 "is not included as a key in this dict, the value associated with the "
54 "'fallback' key will be used.",
55 )
56 minFootprint = Field(dtype=int, default=15,
57 doc="Minimum number of footprints considered sufficient to exit the "
58 "iteration loop. This should be larger than or equal to ``minIsolated`` as "
59 "it does not take into account whether any given footprint is multi-peaked "
60 "(i.e. blended thus not providing any isolated sources for PSF modeling.)")
61 maxPeakToFootRatio = Field(dtype=float, default=800.0,
62 doc="Maximum ratio of the number of peaks per footprint considered "
63 "sufficient to exit the iteration loop. This helps guard against "
64 "large contiguous footprints leaving no isolated sources for "
65 "inclusion in the PSF modeling.")
66 minIsolated = Field(dtype=int, default=6,
67 doc="Minimum number of single-peaked (i.e. isolated) footprints "
68 "considered sufficient to exit the iteration loop. This should be "
69 "larger than the minimum number of sources desired for PSF modeling "
70 "since some may be rejected by the source selector.")
71 sufficientIsolated = Field(dtype=int, default=100,
72 doc="Number of single-peaked (isolated) footprints considered "
73 "sufficient to exit the iteration loop. Must be larger than "
74 "``minIsolated``.")
75 sufficientFractionIsolated = Field(dtype=float, default=0.45,
76 doc="Fraction of single-peaked (isolated) footprints considered "
77 "sufficient to exit the iteration loop.")
79 def validate(self):
80 super().validate()
81 if "fallback" not in self.maxNumPeakPerBand:
82 msg = ("Must include a \"fallback\" key in the config.maxNumPeakPerBand config dict. "
83 f"It is currently: {self.maxNumPeakPerBand}.")
84 raise FieldValidationError(self.__class__.maxNumPeakPerBand, self, msg)
85 if self.minFootprint < self.minIsolated:
86 msg = (f"The config.minFootprint (= {self.minFootprint}) must be >= that of "
87 f"config.minIsolated (= {self.minIsolated}).")
88 raise FieldValidationError(self.__class__.minFootprint, self, msg)
89 if self.sufficientIsolated < self.minIsolated:
90 msg = (f"The config.sufficientIsolated (= {self.sufficientIsolated}) must be >= that of "
91 f"config.minIsolated (= {self.minIsolated}).")
92 raise FieldValidationError(self.__class__.sufficientIsolated, self, msg)
95class AdaptiveThresholdDetectionTask(Task):
96 """Detection of sources on an image using an adaptive scheme for
97 the detection threshold.
98 """
99 ConfigClass = AdaptiveThresholdDetectionConfig
100 _DefaultName = "adaptiveThresholdDetection"
102 def run(self, table, exposure, initialThreshold=None, initialThresholdMultiplier=2.0):
103 """Perform detection with an adaptive threshold detection scheme
104 conditioned to maximize the likelihood of a successful PSF model fit
105 for any given "scene".
107 In particular, we'd like to be able to handle different scenes, from
108 sparsely populated ones through very crowded ones, and possibly high
109 fill-factor nebulosity along the way, with a single pipeline (and set
110 of configs). This requires some flexibility in setting the detection
111 thresholds in order to detect enough sources suitable for PSF modeling
112 (e.g. crowded fields require higher thresholds to ensure the detections
113 don't end up overlapping into a single or very small number of blended
114 footprints.
116 We first detect sources using the default threshold and multiplier.
117 Then, cycling through a series of criteria based on the DETECTED mask
118 planes (number of footprints, number of peaks, number of isolated
119 footprints, number of peaks-per-footrpint, etc.) conditioned to identify
120 a "Goldilocks Zone" where we have enough isolated peaks from which to
121 measure the PSF, we iterate while adjusting the detection thresholds
122 in the appropriate direction until all criteria are met (or the maximum
123 number of iterations is reached).
125 Parameters
126 ----------
127 table : `lsst.afw.table.SourceTable`
128 Table object that will be used to create the SourceCatalog.
129 exposure : `lsst.afw.image.Exposure`
130 Exposure to process; DETECTED mask plane will be set in-place.
131 initialThreshold : `float`, optional
132 Initial threshold for detection of PSF sources.
133 initialThresholdMultiplier : `float`, optional
134 Initial threshold for detection of PSF sources.
136 Returns
137 -------
138 results : `lsst.pipe.base.Struct`
139 The adaptive threshold detection results as a struct with
140 attributes:
142 ``detections``
143 Results of the final round of detection as a struch with
144 attributes:
146 ``sources``
147 Detected sources on the exposure
148 (`lsst.afw.table.SourceCatalog`).
149 ``positive``
150 Positive polarity footprints
151 (`lsst.afw.detection.FootprintSet` or `None`).
152 ``negative``
153 Negative polarity footprints
154 (`lsst.afw.detection.FootprintSet` or `None`).
155 ``numPos``
156 Number of footprints in positive or 0 if detection polarity was
157 negative (`int`).
158 ``numNeg``
159 Number of footprints in negative or 0 if detection polarity was
160 positive (`int`).
161 ``background``
162 Always `None`; provided for compatibility with
163 `SourceDetectionTask`.
164 ``factor``
165 Multiplication factor applied to the configured detection
166 threshold. (`float`).
167 ``thresholdValue``
168 The final threshold value used to the configure the final round
169 of detection (`float`).
170 ``includeThresholdMultiplier``
171 The final multiplication factor applied to the configured detection
172 threshold. (`float`).
173 """
174 band = "fallback"
175 if exposure.filter is not None:
176 if exposure.filter.hasBandLabel():
177 band = exposure.filter.bandLabel
178 if band in self.config.maxNumPeakPerBand:
179 maxNumPeak = self.config.maxNumPeakPerBand[band]
180 else:
181 maxNumPeak = self.config.maxNumPeakPerBand["fallback"]
183 # Set up and configure the adaptive detection task on first iteration.
184 thresholdFactor = 1.0
185 if initialThreshold is None:
186 maxSn = float(np.nanmax(exposure.image.array/np.sqrt(exposure.variance.array)))
187 adaptiveDetThreshold = min(maxSn, 5.0)
188 else:
189 adaptiveDetThreshold = initialThreshold
190 adaptiveDetectionConfig = SourceDetectionConfig()
191 adaptiveDetectionConfig.thresholdValue = adaptiveDetThreshold
192 adaptiveDetectionConfig.includeThresholdMultiplier = initialThresholdMultiplier
193 adaptiveDetectionConfig.reEstimateBackground = False
194 adaptiveDetectionConfig.doTempWideBackground = True
195 adaptiveDetectionConfig.tempWideBackground.binSize = 512
196 adaptiveDetectionConfig.thresholdPolarity = "both"
197 self.log.info("Using adaptive detection with thresholdValue = %.2f and multiplier = %.1f",
198 adaptiveDetectionConfig.thresholdValue,
199 adaptiveDetectionConfig.includeThresholdMultiplier)
200 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
202 maxNumNegFactor = 1.0
203 # We use a 1-indexed iteration variable just to make logs intuitive.
204 for nAdaptiveDetIter in range(1, self.config.maxAdaptiveDetIter + 1):
205 detRes = adaptiveDetectionTask.run(table=table, exposure=exposure, doSmooth=True)
206 sourceCat = detRes.sources
207 nFootprint = len(sourceCat)
208 nPeak = 0
209 nPosPeak = detRes.numPosPeaks
210 nNegPeak = detRes.numNegPeaks # Often detected in high nebulosity scenes.
211 maxNumNegPeak = max(15, int(0.025*nPosPeak))*maxNumNegFactor
212 nIsolated = 0
213 nPeakPerSrcMax = 0
214 maxNumPeakPerSrcMax = 0.2*maxNumPeak
215 for src in sourceCat:
216 nPeakSrc = len(src.getFootprint().getPeaks())
217 if nPeakSrc == 1:
218 nIsolated += 1
219 nPeak += nPeakSrc
220 nPeakPerSrcMax = max(nPeakPerSrcMax, nPeakSrc)
221 if nFootprint > 0.0:
222 fractionIsolated = nIsolated/nFootprint
223 avgPeakPerFoot = nPeak/nFootprint
224 else:
225 fractionIsolated = float("nan")
226 avgPeakPerFoot = float("nan")
227 self.log.info("In adaptive detection iter %d: nFootprints = %d, nPosPeak = %d (max is %d), "
228 "nNegPeak = %d (max is %d), nPeak/nFoot = %.1f (max is %.1f), "
229 "nPeakPerkSrcMax = %d (max is %d), nIsolated = %d, fractionIsolated = %.2f",
230 nAdaptiveDetIter, nFootprint, nPeak, maxNumPeak, nNegPeak, maxNumNegPeak,
231 avgPeakPerFoot, self.config.maxPeakToFootRatio, nPeakPerSrcMax,
232 maxNumPeakPerSrcMax, nIsolated, fractionIsolated)
233 if (nIsolated > self.config.sufficientIsolated
234 and fractionIsolated > self.config.sufficientFractionIsolated
235 and (nAdaptiveDetIter > 1 or self.config.maxAdaptiveDetIter == 1)):
236 if ((nIsolated > 5.0*self.config.sufficientIsolated and nPeak < 2.5*maxNumPeak
237 and nNegPeak < 100.0*maxNumNegPeak)
238 or (nNegPeak < 5.0*maxNumNegPeak and nPeak < 1.2*maxNumPeak)):
239 self.log.info("Sufficient isolated footprints (%d > %d) and fraction of isolated "
240 "footprints (%.2f > %.2f) for PSF modeling. Exiting adaptive detection "
241 "at iter: %d.",
242 nIsolated, self.config.sufficientIsolated, fractionIsolated,
243 self.config.sufficientFractionIsolated, nAdaptiveDetIter)
244 break
246 if nFootprint == 0 or nPosPeak == 0:
247 thresholdFactor = 0.25
248 maxNumNegFactor *= 10
249 self.log.warning("Adaptive threshold increase went too far (nFootprint = 0). "
250 "Decreasing threshold to %.2f and rerunning.",
251 thresholdFactor*adaptiveDetectionConfig.thresholdValue)
252 adaptiveDetectionConfig.thresholdValue = (
253 thresholdFactor*adaptiveDetectionConfig.thresholdValue)
254 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
255 continue
257 if ((nPeak/nFootprint > self.config.maxPeakToFootRatio and nIsolated < self.config.minIsolated)
258 or nNegPeak > maxNumNegPeak):
259 if nNegPeak > 2*maxNumNegPeak:
260 thresholdFactor = 1.5
261 else:
262 thresholdFactor = 1.25
263 thresholdFactor *= adaptiveDetectionConfig.includeThresholdMultiplier
264 adaptiveDetectionConfig.includeThresholdMultiplier = 1.0
265 self.log.warning("Adaptive detection iter %d: catalog had nPeak/nFootprint = "
266 "%.1f (max is %.1f) and %d negative peaks (max is %d). "
267 "Increasing threshold to %.2f and setting multiplier "
268 "to %.1f and rerunning.",
269 nAdaptiveDetIter, nPeak/nFootprint, self.config.maxPeakToFootRatio,
270 nNegPeak, maxNumNegPeak,
271 thresholdFactor*adaptiveDetectionConfig.thresholdValue,
272 adaptiveDetectionConfig.includeThresholdMultiplier)
273 adaptiveDetectionConfig.thresholdValue = (
274 thresholdFactor*adaptiveDetectionConfig.thresholdValue)
275 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
276 continue
278 if (nPeak > maxNumPeak or nPeakPerSrcMax > maxNumPeakPerSrcMax
279 or nFootprint <= self.config.minFootprint):
280 if nPeak > maxNumPeak or nPeakPerSrcMax > 0.25*maxNumPeak:
281 if nAdaptiveDetIter < 0.5*self.config.maxAdaptiveDetIter:
282 if nPeak > 3*maxNumPeak or nPeakPerSrcMax > maxNumPeak:
283 thresholdFactor = 1.7
284 elif nPeak > 2*maxNumPeak or nPeakPerSrcMax > 0.5*maxNumPeak:
285 thresholdFactor = 1.4
286 else:
287 thresholdFactor = 1.2
288 else:
289 thresholdFactor = 1.2
290 thresholdFactor *= adaptiveDetectionConfig.includeThresholdMultiplier
291 newThresholdMultiplier = max(1.0, 0.5*adaptiveDetectionConfig.includeThresholdMultiplier)
292 adaptiveDetectionConfig.includeThresholdMultiplier = newThresholdMultiplier
293 adaptiveDetectionConfig.thresholdValue = (
294 thresholdFactor*adaptiveDetectionConfig.thresholdValue)
295 self.log.warning("Adaptive detection iter %d catalog had nPeak = %d (max = %d) "
296 "and nPeakPerSrcMax = %d (max = %d). Increasing threshold to %.2f "
297 "and setting multiplier to %.1f and rerunning.",
298 nAdaptiveDetIter, nPeak, maxNumPeak, nPeakPerSrcMax, maxNumPeakPerSrcMax,
299 adaptiveDetectionConfig.thresholdValue,
300 adaptiveDetectionConfig.includeThresholdMultiplier)
301 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
302 continue
304 if nFootprint <= self.config.minFootprint:
305 maxNumNegFactor *= 10 # Allow more -ve detections at this point.
306 thresholdFactor = min(0.85, 0.4*np.log10(10*(nFootprint + 1)))
307 self.log.warning("Adaptive detection iter %d catalog had only %d footprints. "
308 "Lowering threshold to %.2f and increasing the allowance for "
309 "negative detections and rerunning.",
310 nAdaptiveDetIter, nFootprint,
311 thresholdFactor*adaptiveDetectionConfig.thresholdValue)
312 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
313 adaptiveDetectionConfig.thresholdValue = (
314 thresholdFactor*adaptiveDetectionConfig.thresholdValue
315 )
316 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
317 else:
318 break
319 # Final round of detection with positive polarity
320 adaptiveDetectionConfig.thresholdPolarity = "positive"
321 adaptiveDetectionTask = SourceDetectionTask(config=adaptiveDetectionConfig)
322 self.log.info("Perfomring final round of detection with threshold %.2f and multiplier %.1f",
323 adaptiveDetectionConfig.thresholdValue,
324 adaptiveDetectionConfig.includeThresholdMultiplier)
325 detRes = adaptiveDetectionTask.run(table=table, exposure=exposure, doSmooth=True)
326 return Struct(
327 detections=detRes,
328 thresholdValue=adaptiveDetectionConfig.thresholdValue,
329 includeThresholdMultiplier=adaptiveDetectionConfig.includeThresholdMultiplier,
330 )