Coverage for python / lsst / pipe / tasks / characterizeImage.py: 28%
193 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 08:53 +0000
1# This file is part of pipe_tasks.
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__ = ["CharacterizeImageConfig", "CharacterizeImageTask"]
24import numpy as np
26from lsstDebug import getDebugFrame
27import lsst.afw.table as afwTable
28import lsst.pex.config as pexConfig
29import lsst.pipe.base as pipeBase
30import lsst.daf.base as dafBase
31import lsst.pipe.base.connectionTypes as cT
32from lsst.afw.math import BackgroundList
33from lsst.afw.table import SourceTable
34from lsst.meas.algorithms import (
35 SubtractBackgroundTask,
36 SourceDetectionTask,
37 MeasureApCorrTask,
38 MeasureApCorrError,
39 MaskStreaksTask,
40 NormalizedCalibrationFluxTask,
41)
42from lsst.meas.algorithms.installGaussianPsf import InstallGaussianPsfTask
43from lsst.meas.astrom import displayAstrometry
44from lsst.meas.base import (
45 SingleFrameMeasurementTask,
46 ApplyApCorrTask,
47 CatalogCalculationTask,
48 IdGenerator,
49 DetectorVisitIdGeneratorConfig,
50)
51from lsst.meas.deblender import SourceDeblendTask
52import lsst.meas.extensions.shapeHSM # noqa: F401 needed for default shape plugin
53from .measurePsf import MeasurePsfTask
54from .repair import RepairTask
55from .computeExposureSummaryStats import ComputeExposureSummaryStatsTask
56from lsst.pex.exceptions import LengthError
57from lsst.utils.timer import timeMethod
60class CharacterizeImageConnections(pipeBase.PipelineTaskConnections,
61 dimensions=("instrument", "visit", "detector")):
62 exposure = cT.Input(
63 doc="Input exposure data",
64 name="postISRCCD",
65 storageClass="Exposure",
66 dimensions=["instrument", "exposure", "detector"],
67 )
68 characterized = cT.Output(
69 doc="Output characterized data.",
70 name="icExp",
71 storageClass="ExposureF",
72 dimensions=["instrument", "visit", "detector"],
73 )
74 sourceCat = cT.Output(
75 doc="Output source catalog.",
76 name="icSrc",
77 storageClass="SourceCatalog",
78 dimensions=["instrument", "visit", "detector"],
79 )
80 backgroundModel = cT.Output(
81 doc="Output background model.",
82 name="icExpBackground",
83 storageClass="Background",
84 dimensions=["instrument", "visit", "detector"],
85 )
86 outputSchema = cT.InitOutput(
87 doc="Schema of the catalog produced by CharacterizeImage",
88 name="icSrc_schema",
89 storageClass="SourceCatalog",
90 )
92 def adjustQuantum(self, inputs, outputs, label, dataId):
93 # Docstring inherited from PipelineTaskConnections
94 try:
95 return super().adjustQuantum(inputs, outputs, label, dataId)
96 except pipeBase.ScalarError as err:
97 raise pipeBase.ScalarError(
98 "CharacterizeImageTask can at present only be run on visits that are associated with "
99 "exactly one exposure. Either this is not a valid exposure for this pipeline, or the "
100 "snap-combination step you probably want hasn't been configured to run between ISR and "
101 "this task (as of this writing, that would be because it hasn't been implemented yet)."
102 ) from err
105class CharacterizeImageConfig(pipeBase.PipelineTaskConfig,
106 pipelineConnections=CharacterizeImageConnections):
107 """Config for CharacterizeImageTask."""
109 doMeasurePsf = pexConfig.Field(
110 dtype=bool,
111 default=True,
112 doc="Measure PSF? If False then for all subsequent operations use either existing PSF "
113 "model when present, or install simple PSF model when not (see installSimplePsf "
114 "config options)"
115 )
116 doWrite = pexConfig.Field(
117 dtype=bool,
118 default=True,
119 doc="Persist results?",
120 )
121 doWriteExposure = pexConfig.Field(
122 dtype=bool,
123 default=True,
124 doc="Write icExp and icExpBackground in addition to icSrc? Ignored if doWrite False.",
125 )
126 psfIterations = pexConfig.RangeField(
127 dtype=int,
128 default=2,
129 min=1,
130 doc="Number of iterations of detect sources, measure sources, "
131 "estimate PSF. If useSimplePsf is True then 2 should be plenty; "
132 "otherwise more may be wanted.",
133 )
134 background = pexConfig.ConfigurableField(
135 target=SubtractBackgroundTask,
136 doc="Configuration for initial background estimation",
137 )
138 detection = pexConfig.ConfigurableField(
139 target=SourceDetectionTask,
140 doc="Detect sources"
141 )
142 doDeblend = pexConfig.Field(
143 dtype=bool,
144 default=True,
145 doc="Run deblender input exposure"
146 )
147 deblend = pexConfig.ConfigurableField(
148 target=SourceDeblendTask,
149 doc="Split blended source into their components"
150 )
151 measurement = pexConfig.ConfigurableField(
152 target=SingleFrameMeasurementTask,
153 doc="Measure sources"
154 )
155 doNormalizedCalibration = pexConfig.Field(
156 dtype=bool,
157 default=True,
158 doc="Use normalized calibration flux (e.g. compensated tophats)?",
159 )
160 normalizedCalibrationFlux = pexConfig.ConfigurableField(
161 target=NormalizedCalibrationFluxTask,
162 doc="Task to normalize the calibration flux (e.g. compensated tophats).",
163 )
164 doApCorr = pexConfig.Field(
165 dtype=bool,
166 default=True,
167 doc="Run subtasks to measure and apply aperture corrections"
168 )
169 measureApCorr = pexConfig.ConfigurableField(
170 target=MeasureApCorrTask,
171 doc="Subtask to measure aperture corrections"
172 )
173 applyApCorr = pexConfig.ConfigurableField(
174 target=ApplyApCorrTask,
175 doc="Subtask to apply aperture corrections"
176 )
177 # If doApCorr is False, and the exposure does not have apcorrections already applied, the
178 # active plugins in catalogCalculation almost certainly should not contain the characterization plugin
179 catalogCalculation = pexConfig.ConfigurableField(
180 target=CatalogCalculationTask,
181 doc="Subtask to run catalogCalculation plugins on catalog"
182 )
183 doComputeSummaryStats = pexConfig.Field(
184 dtype=bool,
185 default=True,
186 doc="Run subtask to measure exposure summary statistics",
187 deprecated=("This subtask has been moved to CalibrateTask "
188 "with DM-30701.")
189 )
190 computeSummaryStats = pexConfig.ConfigurableField(
191 target=ComputeExposureSummaryStatsTask,
192 doc="Subtask to run computeSummaryStats on exposure",
193 deprecated=("This subtask has been moved to CalibrateTask "
194 "with DM-30701.")
195 )
196 useSimplePsf = pexConfig.Field(
197 dtype=bool,
198 default=True,
199 doc="Replace the existing PSF model with a simplified version that has the same sigma "
200 "at the start of each PSF determination iteration? Doing so makes PSF determination "
201 "converge more robustly and quickly.",
202 )
203 installSimplePsf = pexConfig.ConfigurableField(
204 target=InstallGaussianPsfTask,
205 doc="Install a simple PSF model",
206 )
207 measurePsf = pexConfig.ConfigurableField(
208 target=MeasurePsfTask,
209 doc="Measure PSF",
210 )
211 repair = pexConfig.ConfigurableField(
212 target=RepairTask,
213 doc="Remove cosmic rays",
214 )
215 requireCrForPsf = pexConfig.Field(
216 dtype=bool,
217 default=True,
218 doc="Require cosmic ray detection and masking to run successfully before measuring the PSF."
219 )
220 checkUnitsParseStrict = pexConfig.Field(
221 doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'",
222 dtype=str,
223 default="raise",
224 )
225 doMaskStreaks = pexConfig.Field(
226 doc="Mask streaks",
227 default=False,
228 dtype=bool,
229 deprecated=("This subtask has been moved to detectAndMeasureTask in "
230 "ip_diffim with DM-43370 and will be removed in DM-44658.")
231 )
232 maskStreaks = pexConfig.ConfigurableField(
233 target=MaskStreaksTask,
234 doc="Subtask for masking streaks. Only used if doMaskStreaks is True. "
235 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.",
236 deprecated=("This subtask has been moved to detectAndMeasureTask in "
237 "ip_diffim with DM-43370 and will be removed in DM-44658.")
238 )
239 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
241 def setDefaults(self):
242 super().setDefaults()
243 # Just detect bright stars.
244 # The thresholdValue sets the minimum flux in a pixel to be included in the
245 # footprint, while peaks are only detected when they are above
246 # thresholdValue * includeThresholdMultiplier. The low thresholdValue
247 # ensures that the footprints are large enough for the noise replacer
248 # to mask out faint undetected neighbors that are not to be measured.
249 self.detection.thresholdValue = 5.0
250 self.detection.includeThresholdMultiplier = 10.0
251 # do not deblend, as it makes a mess
252 self.doDeblend = False
253 # measure and apply aperture correction; note: measuring and applying aperture
254 # correction are disabled until the final measurement, after PSF is measured
255 self.doApCorr = True
256 # During characterization, we don't have full source measurement information,
257 # so must do the aperture correction with only psf stars, combined with the
258 # default signal-to-noise cuts in MeasureApCorrTask.
259 selector = self.measureApCorr.sourceSelector["science"]
260 selector.doUnresolved = False
261 selector.flags.good = ["calib_psf_used"]
262 selector.flags.bad = []
264 # minimal set of measurements needed to determine PSF
265 self.measurement.plugins.names = [
266 "base_PixelFlags",
267 "base_SdssCentroid",
268 "ext_shapeHSM_HsmSourceMoments",
269 "base_GaussianFlux",
270 "base_PsfFlux",
271 "base_CircularApertureFlux",
272 "base_CompensatedTophatFlux",
273 ]
274 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
275 self.measurement.algorithms["base_CompensatedTophatFlux"].apertures = [12]
277 def validate(self):
278 if self.doApCorr and not self.measurePsf:
279 raise RuntimeError("Must measure PSF to measure aperture correction, "
280 "because flags determined by PSF measurement are used to identify "
281 "sources used to measure aperture correction")
284class CharacterizeImageTask(pipeBase.PipelineTask):
285 """Measure bright sources and use this to estimate background and PSF of
286 an exposure.
288 Given an exposure with defects repaired (masked and interpolated over,
289 e.g. as output by `~lsst.ip.isr.IsrTask`):
290 - detect and measure bright sources
291 - repair cosmic rays
292 - measure and subtract background
293 - measure PSF
295 Parameters
296 ----------
297 schema : `lsst.afw.table.Schema`, optional
298 Initial schema for icSrc catalog.
299 **kwargs
300 Additional keyword arguments.
302 Notes
303 -----
304 Debugging:
305 CharacterizeImageTask has a debug dictionary with the following keys:
307 frame
308 int: if specified, the frame of first debug image displayed (defaults to 1)
309 repair_iter
310 bool; if True display image after each repair in the measure PSF loop
311 background_iter
312 bool; if True display image after each background subtraction in the measure PSF loop
313 measure_iter
314 bool; if True display image and sources at the end of each iteration of the measure PSF loop
315 See `~lsst.meas.astrom.displayAstrometry` for the meaning of the various symbols.
316 psf
317 bool; if True display image and sources after PSF is measured;
318 this will be identical to the final image displayed by measure_iter if measure_iter is true
319 repair
320 bool; if True display image and sources after final repair
321 measure
322 bool; if True display image and sources after final measurement
323 """
325 ConfigClass = CharacterizeImageConfig
326 _DefaultName = "characterizeImage"
328 def __init__(self, schema=None, **kwargs):
329 super().__init__(**kwargs)
331 if schema is None:
332 schema = SourceTable.makeMinimalSchema()
333 self.schema = schema
334 self.makeSubtask("background")
335 self.makeSubtask("installSimplePsf")
336 self.makeSubtask("repair")
337 # TODO: DM-44658, streak masking to happen only in ip_diffim
338 if self.config.doMaskStreaks:
339 self.makeSubtask("maskStreaks")
340 self.makeSubtask("measurePsf", schema=self.schema)
341 self.algMetadata = dafBase.PropertyList()
342 self.makeSubtask('detection', schema=self.schema)
343 if self.config.doDeblend:
344 self.makeSubtask("deblend", schema=self.schema)
345 self.makeSubtask('measurement', schema=self.schema, algMetadata=self.algMetadata)
346 if self.config.doNormalizedCalibration:
347 self.makeSubtask('normalizedCalibrationFlux', schema=self.schema)
348 if self.config.doApCorr:
349 self.makeSubtask('measureApCorr', schema=self.schema)
350 self.makeSubtask('applyApCorr', schema=self.schema)
351 self.makeSubtask('catalogCalculation', schema=self.schema)
352 self._initialFrame = getDebugFrame(self._display, "frame") or 1
353 self._frame = self._initialFrame
354 self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict)
355 afwTable.CoordKey.addErrorFields(self.schema)
356 self.outputSchema = afwTable.SourceCatalog(self.schema)
358 def runQuantum(self, butlerQC, inputRefs, outputRefs):
359 inputs = butlerQC.get(inputRefs)
360 if 'idGenerator' not in inputs.keys():
361 inputs['idGenerator'] = self.config.idGenerator.apply(butlerQC.quantum.dataId)
362 outputs = self.run(**inputs)
363 butlerQC.put(outputs, outputRefs)
365 @timeMethod
366 def run(self, exposure, background=None, idGenerator=None):
367 """Characterize a science image.
369 Peforms the following operations:
370 - Iterate the following config.psfIterations times, or once if config.doMeasurePsf false:
371 - detect and measure sources and estimate PSF (see detectMeasureAndEstimatePsf for details)
372 - interpolate over cosmic rays
373 - perform final measurement
375 Parameters
376 ----------
377 exposure : `lsst.afw.image.ExposureF`
378 Exposure to characterize.
379 background : `lsst.afw.math.BackgroundList`, optional
380 Initial model of background already subtracted from exposure.
381 idGenerator : `lsst.meas.base.IdGenerator`, optional
382 Object that generates source IDs and provides RNG seeds.
384 Returns
385 -------
386 result : `lsst.pipe.base.Struct`
387 Results as a struct with attributes:
389 ``exposure``
390 Characterized exposure (`lsst.afw.image.ExposureF`).
391 ``sourceCat``
392 Detected sources (`lsst.afw.table.SourceCatalog`).
393 ``background``
394 Model of subtracted background (`lsst.afw.math.BackgroundList`).
395 ``psfCellSet``
396 Spatial cells of PSF candidates (`lsst.afw.math.SpatialCellSet`).
397 ``characterized``
398 Another reference to ``exposure`` for compatibility.
399 ``backgroundModel``
400 Another reference to ``background`` for compatibility.
402 Raises
403 ------
404 RuntimeError
405 Raised if PSF sigma is NaN.
406 """
407 self._frame = self._initialFrame # reset debug display frame
409 if not self.config.doMeasurePsf and not exposure.hasPsf():
410 self.log.info("CharacterizeImageTask initialized with 'simple' PSF.")
411 self.installSimplePsf.run(exposure=exposure)
413 if idGenerator is None:
414 idGenerator = IdGenerator()
416 # subtract an initial estimate of background level
417 background = self.background.run(exposure).background
419 psfIterations = self.config.psfIterations if self.config.doMeasurePsf else 1
420 for i in range(psfIterations):
421 dmeRes = self.detectMeasureAndEstimatePsf(
422 exposure=exposure,
423 idGenerator=idGenerator,
424 background=background,
425 )
427 psf = dmeRes.exposure.getPsf()
428 # Just need a rough estimate; average positions are fine
429 psfAvgPos = psf.getAveragePosition()
430 psfSigma = psf.computeShape(psfAvgPos).getDeterminantRadius()
431 psfDimensions = psf.computeImage(psfAvgPos).getDimensions()
432 medBackground = np.median(dmeRes.background.getImage().getArray())
433 self.log.info("iter %s; PSF sigma=%0.4f, dimensions=%s; median background=%0.2f",
434 i + 1, psfSigma, psfDimensions, medBackground)
435 if np.isnan(psfSigma):
436 raise RuntimeError("PSF sigma is NaN, cannot continue PSF determination.")
438 self.display("psf", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
440 # perform final repair with final PSF
441 self.repair.run(exposure=dmeRes.exposure)
442 self.display("repair", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
444 # mask streaks
445 # TODO: Remove in DM-44658, streak masking to happen only in ip_diffim
446 if self.config.doMaskStreaks:
447 _ = self.maskStreaks.run(dmeRes.exposure)
449 # perform final measurement with final PSF, including measuring and applying aperture correction,
450 # if wanted
451 self.measurement.run(measCat=dmeRes.sourceCat, exposure=dmeRes.exposure,
452 exposureId=idGenerator.catalog_id)
454 if self.config.doNormalizedCalibration:
455 normApCorrMap = self.normalizedCalibrationFlux.run(
456 exposure=dmeRes.exposure,
457 catalog=dmeRes.sourceCat,
458 ).ap_corr_map
459 dmeRes.exposure.info.setApCorrMap(normApCorrMap)
460 else:
461 normApCorrMap = None
463 if self.config.doApCorr:
464 # This aperture correction is relative to slot_CalibFlux_instFlux
465 # which is now set to the normalized calibration flux if that
466 # has been run.
467 try:
468 apCorrMap = self.measureApCorr.run(
469 exposure=dmeRes.exposure,
470 catalog=dmeRes.sourceCat,
471 ).apCorrMap
472 except MeasureApCorrError:
473 # We have failed to get a valid aperture correction map.
474 # Proceed with processing, and image will be filtered
475 # downstream.
476 dmeRes.exposure.info.setApCorrMap(None)
477 else:
478 # Need to merge the aperture correction map from the normalization.
479 if normApCorrMap:
480 for key in normApCorrMap:
481 apCorrMap[key] = normApCorrMap[key]
482 dmeRes.exposure.info.setApCorrMap(apCorrMap)
483 self.applyApCorr.run(catalog=dmeRes.sourceCat, apCorrMap=exposure.getInfo().getApCorrMap())
485 self.catalogCalculation.run(dmeRes.sourceCat)
487 self.display("measure", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
489 return pipeBase.Struct(
490 exposure=dmeRes.exposure,
491 sourceCat=dmeRes.sourceCat,
492 background=dmeRes.background,
493 psfCellSet=dmeRes.psfCellSet,
495 characterized=dmeRes.exposure,
496 backgroundModel=dmeRes.background
497 )
499 @timeMethod
500 def detectMeasureAndEstimatePsf(self, exposure, idGenerator, background):
501 """Perform one iteration of detect, measure, and estimate PSF.
503 Performs the following operations:
505 - if config.doMeasurePsf or not exposure.hasPsf():
507 - install a simple PSF model (replacing the existing one, if need be)
509 - interpolate over cosmic rays with keepCRs=True
510 - estimate background and subtract it from the exposure
511 - detect, deblend and measure sources, and subtract a refined background model;
512 - if config.doMeasurePsf:
513 - measure PSF
515 Parameters
516 ----------
517 exposure : `lsst.afw.image.ExposureF`
518 Exposure to characterize.
519 idGenerator : `lsst.meas.base.IdGenerator`
520 Object that generates source IDs and provides RNG seeds.
521 background : `lsst.afw.math.BackgroundList`, optional
522 Initial model of background already subtracted from exposure.
524 Returns
525 -------
526 result : `lsst.pipe.base.Struct`
527 Results as a struct with attributes:
529 ``exposure``
530 Characterized exposure (`lsst.afw.image.ExposureF`).
531 ``sourceCat``
532 Detected sources (`lsst.afw.table.SourceCatalog`).
533 ``background``
534 Model of subtracted background (`lsst.afw.math.BackgroundList`).
535 ``psfCellSet``
536 Spatial cells of PSF candidates (`lsst.afw.math.SpatialCellSet`).
538 Raises
539 ------
540 LengthError
541 Raised if there are too many CR pixels.
542 """
543 # install a simple PSF model, if needed or wanted
544 if not exposure.hasPsf() or (self.config.doMeasurePsf and self.config.useSimplePsf):
545 self.log.info("PSF estimation initialized with 'simple' PSF")
546 self.installSimplePsf.run(exposure=exposure)
548 # run repair, but do not interpolate over cosmic rays (do that elsewhere, with the final PSF model)
549 if self.config.requireCrForPsf:
550 self.repair.run(exposure=exposure, keepCRs=True)
551 else:
552 try:
553 self.repair.run(exposure=exposure, keepCRs=True)
554 except LengthError:
555 self.log.warning("Skipping cosmic ray detection: Too many CR pixels (max %0.f)",
556 self.config.repair.cosmicray.nCrPixelMax)
558 self.display("repair_iter", exposure=exposure)
560 if background is None:
561 background = BackgroundList()
562 self.schema.addField(
563 'psf_max_value',
564 type=np.float32,
565 doc="PSF max value.",
566 doReplace=True,
567 )
568 sourceIdFactory = idGenerator.make_table_id_factory()
569 table = SourceTable.make(self.schema, sourceIdFactory)
570 table.setMetadata(self.algMetadata)
572 detRes = self.detection.run(table=table, exposure=exposure, doSmooth=True)
573 sourceCat = detRes.sources
574 if detRes.background:
575 for bg in detRes.background:
576 background.append(bg)
578 if self.config.doDeblend:
579 self.deblend.run(exposure=exposure, sources=sourceCat)
580 # We need the output catalog to be contiguous for further processing.
581 if not sourceCat.isContiguous():
582 sourceCat = sourceCat.copy(deep=True)
584 self.measurement.run(measCat=sourceCat, exposure=exposure, exposureId=idGenerator.catalog_id)
586 measPsfRes = pipeBase.Struct(cellSet=None)
587 if self.config.doMeasurePsf:
588 measPsfRes = self.measurePsf.run(exposure=exposure, sources=sourceCat,
589 expId=idGenerator.catalog_id)
590 self.display("measure_iter", exposure=exposure, sourceCat=sourceCat)
592 return pipeBase.Struct(
593 exposure=exposure,
594 sourceCat=sourceCat,
595 background=background,
596 psfCellSet=measPsfRes.cellSet,
597 )
599 def display(self, itemName, exposure, sourceCat=None):
600 """Display exposure and sources on next frame (for debugging).
602 Parameters
603 ----------
604 itemName : `str`
605 Name of item in ``debugInfo``.
606 exposure : `lsst.afw.image.ExposureF`
607 Exposure to display.
608 sourceCat : `lsst.afw.table.SourceCatalog`, optional
609 Catalog of sources detected on the exposure.
610 """
611 val = getDebugFrame(self._display, itemName)
612 if not val:
613 return
615 displayAstrometry(exposure=exposure, sourceCat=sourceCat, frame=self._frame, pause=False)
616 self._frame += 1