Coverage for python / lsst / pipe / tasks / calibrateImage.py: 12%
784 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-20 01:17 -0800
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__ = ["CalibrateImageTask", "CalibrateImageConfig", "NoPsfStarsToStarsMatchError",
23 "AllCentroidsFlaggedError"]
25import math
26import numpy as np
27import requests
28import os
30from lsst.afw.geom import SpanSet
31import lsst.afw.table as afwTable
32import lsst.afw.image as afwImage
33import lsst.afw.math as afwMath
34from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache
35import lsst.meas.algorithms
36import lsst.meas.algorithms.installGaussianPsf
37import lsst.meas.algorithms.measureApCorr
38import lsst.meas.algorithms.setPrimaryFlags
39from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask
40import lsst.meas.base
41import lsst.meas.astrom
42import lsst.meas.deblender
43import lsst.meas.extensions.shapeHSM
44from lsst.obs.base import createInitialSkyWcs
45import lsst.pex.config as pexConfig
46import lsst.pipe.base as pipeBase
47from lsst.pipe.base import connectionTypes
48from lsst.utils.timer import timeMethod
50from . import measurePsf, repair, photoCal, computeExposureSummaryStats, snapCombine, diffractionSpikeMask
53class AllCentroidsFlaggedError(pipeBase.AlgorithmError):
54 """Raised when there are no valid centroids during psf fitting.
55 """
56 def __init__(self, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size):
57 msg = (f"All source centroids (out of {n_sources}) flagged during PSF fitting. "
58 "Original image PSF is likely unusable; best-fit PSF shape parameters: "
59 f"Ixx={psf_shape_ixx}, Iyy={psf_shape_iyy}, Ixy={psf_shape_ixy}, size={psf_size}"
60 )
61 super().__init__(msg)
62 self.n_sources = n_sources
63 self.psf_shape_ixx = psf_shape_ixx
64 self.psf_shape_iyy = psf_shape_iyy
65 self.psf_shape_ixy = psf_shape_ixy
66 self.psf_size = psf_size
68 @property
69 def metadata(self):
70 return {"n_sources": self.n_sources,
71 "psf_shape_ixx": self.psf_shape_ixx,
72 "psf_shape_iyy": self.psf_shape_iyy,
73 "psf_shape_ixy": self.psf_shape_ixy,
74 "psf_size": self.psf_size
75 }
78class NoPsfStarsToStarsMatchError(pipeBase.AlgorithmError):
79 """Raised when there are no matches between the psf_stars and stars
80 catalogs.
81 """
82 def __init__(self, *, n_psf_stars, n_stars):
83 msg = (f"No psf stars out of {n_psf_stars} matched {n_stars} calib stars."
84 " Downstream processes probably won't have useful {calib_type} stars in this case."
85 " Is `star_selector` too strict or is this a bad image?")
86 super().__init__(msg)
87 self.n_psf_stars = n_psf_stars
88 self.n_stars = n_stars
90 @property
91 def metadata(self):
92 return {"n_psf_stars": self.n_psf_stars,
93 "n_stars": self.n_stars
94 }
97class CalibrateImageConnections(pipeBase.PipelineTaskConnections,
98 dimensions=("instrument", "visit", "detector")):
100 astrometry_ref_cat = connectionTypes.PrerequisiteInput(
101 doc="Reference catalog to use for astrometric calibration.",
102 name="gaia_dr3_20230707",
103 storageClass="SimpleCatalog",
104 dimensions=("skypix",),
105 deferLoad=True,
106 multiple=True,
107 )
108 photometry_ref_cat = connectionTypes.PrerequisiteInput(
109 doc="Reference catalog to use for photometric calibration.",
110 name="ps1_pv3_3pi_20170110",
111 storageClass="SimpleCatalog",
112 dimensions=("skypix",),
113 deferLoad=True,
114 multiple=True
115 )
116 camera_model = connectionTypes.PrerequisiteInput(
117 doc="Camera distortion model to use for astrometric calibration.",
118 name="astrometry_camera",
119 storageClass="Camera",
120 dimensions=("instrument", "physical_filter"),
121 isCalibration=True,
122 minimum=0, # Can fall back on WCS already attached to exposure.
123 )
124 exposures = connectionTypes.Input(
125 doc="Exposure (or two snaps) to be calibrated, and detected and measured on.",
126 name="postISRCCD",
127 storageClass="Exposure",
128 multiple=True, # to handle 1 exposure or 2 snaps
129 dimensions=["instrument", "exposure", "detector"],
130 )
132 background_flat = connectionTypes.PrerequisiteInput(
133 name="flat",
134 doc="Flat calibration frame used for background correction.",
135 storageClass="ExposureF",
136 dimensions=["instrument", "detector", "physical_filter"],
137 isCalibration=True,
138 )
139 illumination_correction = connectionTypes.PrerequisiteInput(
140 name="illuminationCorrection",
141 doc="Illumination correction frame.",
142 storageClass="Exposure",
143 dimensions=["instrument", "detector", "physical_filter"],
144 isCalibration=True,
145 )
147 # outputs
148 initial_stars_schema = connectionTypes.InitOutput(
149 doc="Schema of the output initial stars catalog.",
150 name="initial_stars_schema",
151 storageClass="SourceCatalog",
152 )
154 # TODO DM-38732: We want some kind of flag on Exposures/Catalogs to make
155 # it obvious which components had failed to be computed/persisted.
156 exposure = connectionTypes.Output(
157 doc="Photometrically calibrated, background-subtracted exposure with fitted calibrations and "
158 "summary statistics. To recover the original exposure, first add the background "
159 "(`initial_pvi_background`), and then uncalibrate (divide by `initial_photoCalib_detector`).",
160 name="initial_pvi",
161 storageClass="ExposureF",
162 dimensions=("instrument", "visit", "detector"),
163 )
164 stars = connectionTypes.Output(
165 doc="Catalog of unresolved sources detected on the calibrated exposure.",
166 name="initial_stars_detector",
167 storageClass="ArrowAstropy",
168 dimensions=["instrument", "visit", "detector"],
169 )
170 stars_footprints = connectionTypes.Output(
171 doc="Catalog of unresolved sources detected on the calibrated exposure; "
172 "includes source footprints.",
173 name="initial_stars_footprints_detector",
174 storageClass="SourceCatalog",
175 dimensions=["instrument", "visit", "detector"],
176 )
177 applied_photo_calib = connectionTypes.Output(
178 doc=(
179 "Photometric calibration that was applied to exposure's pixels. "
180 "This connection is disabled when do_calibrate_pixels=False."
181 ),
182 name="initial_photoCalib_detector",
183 storageClass="PhotoCalib",
184 dimensions=("instrument", "visit", "detector"),
185 )
186 background = connectionTypes.Output(
187 doc="Background models estimated during calibration task; calibrated to be in nJy units.",
188 name="initial_pvi_background",
189 storageClass="Background",
190 dimensions=("instrument", "visit", "detector"),
191 )
192 background_to_photometric_ratio = connectionTypes.Output(
193 doc="Ratio of a background-flattened image to a photometric-flattened image. Only persisted "
194 "if do_illumination_correction is True.",
195 name="background_to_photometric_ratio",
196 storageClass="Image",
197 dimensions=("instrument", "visit", "detector"),
198 )
200 # Optional outputs
201 mask = connectionTypes.Output(
202 doc="The mask plane of the calibrated exposure, written as a separate "
203 "output so that it can be passed to downstream tasks even if the "
204 "exposure is removed to save storage.",
205 name="preliminary_visit_mask",
206 storageClass="Mask",
207 dimensions=("instrument", "visit", "detector"),
208 )
209 psf_stars_footprints = connectionTypes.Output(
210 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination; "
211 "includes source footprints.",
212 name="initial_psf_stars_footprints_detector",
213 storageClass="SourceCatalog",
214 dimensions=["instrument", "visit", "detector"],
215 )
216 psf_stars = connectionTypes.Output(
217 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination.",
218 name="initial_psf_stars_detector",
219 storageClass="ArrowAstropy",
220 dimensions=["instrument", "visit", "detector"],
221 )
222 astrometry_matches = connectionTypes.Output(
223 doc="Source to reference catalog matches from the astrometry solver.",
224 name="initial_astrometry_match_detector",
225 storageClass="Catalog",
226 dimensions=("instrument", "visit", "detector"),
227 )
228 photometry_matches = connectionTypes.Output(
229 doc="Source to reference catalog matches from the photometry solver.",
230 name="initial_photometry_match_detector",
231 storageClass="Catalog",
232 dimensions=("instrument", "visit", "detector"),
233 )
235 def __init__(self, *, config=None):
236 super().__init__(config=config)
237 if "psf_stars" not in config.optional_outputs:
238 del self.psf_stars
239 if "psf_stars_footprints" not in config.optional_outputs:
240 del self.psf_stars_footprints
241 if "astrometry_matches" not in config.optional_outputs:
242 del self.astrometry_matches
243 if "photometry_matches" not in config.optional_outputs:
244 del self.photometry_matches
245 if "mask" not in config.optional_outputs:
246 del self.mask
247 if not config.do_calibrate_pixels:
248 del self.applied_photo_calib
249 if not config.do_illumination_correction:
250 del self.background_flat
251 del self.illumination_correction
252 del self.background_to_photometric_ratio
253 if not config.useButlerCamera:
254 del self.camera_model
257class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CalibrateImageConnections):
258 optional_outputs = pexConfig.ListField(
259 doc="Which optional outputs to save (as their connection name)?",
260 dtype=str,
261 # TODO: note somewhere to disable this for benchmarking, but should
262 # we always have it on for production runs?
263 default=["psf_stars", "psf_stars_footprints", "astrometry_matches", "photometry_matches", "mask"],
264 optional=False
265 )
267 # To generate catalog ids consistently across subtasks.
268 id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()
270 snap_combine = pexConfig.ConfigurableField(
271 target=snapCombine.SnapCombineTask,
272 doc="Task to combine two snaps to make one exposure.",
273 )
275 # subtasks used during psf characterization
276 install_simple_psf = pexConfig.ConfigurableField(
277 target=lsst.meas.algorithms.installGaussianPsf.InstallGaussianPsfTask,
278 doc="Task to install a simple PSF model into the input exposure to use "
279 "when detecting bright sources for PSF estimation.",
280 )
281 psf_repair = pexConfig.ConfigurableField(
282 target=repair.RepairTask,
283 doc="Task to repair cosmic rays on the exposure before PSF determination.",
284 )
285 psf_subtract_background = pexConfig.ConfigurableField(
286 target=lsst.meas.algorithms.SubtractBackgroundTask,
287 doc="Task to perform initial background subtraction, before first detection pass.",
288 )
289 psf_detection = pexConfig.ConfigurableField(
290 target=lsst.meas.algorithms.SourceDetectionTask,
291 doc="Task to detect sources for PSF determination."
292 )
293 do_adaptive_threshold_detection = pexConfig.Field(
294 dtype=bool,
295 default=True,
296 doc="Implement the adaptive detection thresholding approach?",
297 )
298 psf_adaptive_threshold_detection = pexConfig.ConfigurableField(
299 target=AdaptiveThresholdDetectionTask,
300 doc="Task to adaptively detect sources for PSF determination.",
301 )
302 psf_source_measurement = pexConfig.ConfigurableField(
303 target=lsst.meas.base.SingleFrameMeasurementTask,
304 doc="Task to measure sources to be used for psf estimation."
305 )
306 psf_measure_psf = pexConfig.ConfigurableField(
307 target=measurePsf.MeasurePsfTask,
308 doc="Task to measure the psf on bright sources."
309 )
310 psf_normalized_calibration_flux = pexConfig.ConfigurableField(
311 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
312 doc="Task to normalize the calibration flux (e.g. compensated tophats) "
313 "for the bright stars used for psf estimation.",
314 )
316 # TODO DM-39203: we can remove aperture correction from this task once we
317 # are using the shape-based star/galaxy code.
318 measure_aperture_correction = pexConfig.ConfigurableField(
319 target=lsst.meas.algorithms.measureApCorr.MeasureApCorrTask,
320 doc="Task to compute the aperture correction from the bright stars."
321 )
323 # subtasks used during star measurement
324 star_background = pexConfig.ConfigurableField(
325 target=lsst.meas.algorithms.SubtractBackgroundTask,
326 doc="Task to perform final background subtraction, just before photoCal.",
327 )
328 star_background_min_footprints = pexConfig.Field(
329 dtype=int,
330 default=3,
331 doc="Minimum number of footprints in the detection mask for star_background measurement. "
332 "This number will get adjusted to the fraction config.star_background_peak_fraction "
333 "of the detected peaks if that number is larger. If the number of footprints is less "
334 "than the minimum, the detection threshold is iteratively increased until the "
335 "threshold is met.",
336 )
337 star_background_peak_fraction = pexConfig.Field(
338 dtype=float,
339 default=0.01,
340 doc="The minimum number of footprints in the detection mask for star_background measurement. "
341 "gets set to the maximum of this fraction of the detected peaks and the value set in "
342 "config.star_background_min_footprints. If the number of footprints is less than the "
343 "current minimum set, the detection threshold is iteratively increased until the "
344 "threshold is met.",
345 )
346 star_detection = pexConfig.ConfigurableField(
347 target=lsst.meas.algorithms.SourceDetectionTask,
348 doc="Task to detect stars to return in the output catalog."
349 )
350 star_sky_sources = pexConfig.ConfigurableField(
351 target=lsst.meas.algorithms.SkyObjectsTask,
352 doc="Task to generate sky sources ('empty' regions where there are no detections).",
353 )
354 do_downsample_footprints = pexConfig.Field(
355 dtype=bool,
356 default=False,
357 doc="Downsample footprints prior to deblending to optimize speed?",
358 )
359 downsample_max_footprints = pexConfig.Field(
360 dtype=int,
361 default=1000,
362 doc="Maximum number of non-sky-source footprints to use if do_downsample_footprints is True,",
363 )
364 star_deblend = pexConfig.ConfigurableField(
365 target=lsst.meas.deblender.SourceDeblendTask,
366 doc="Split blended sources into their components."
367 )
368 star_measurement = pexConfig.ConfigurableField(
369 target=lsst.meas.base.SingleFrameMeasurementTask,
370 doc="Task to measure stars to return in the output catalog."
371 )
372 star_normalized_calibration_flux = pexConfig.ConfigurableField(
373 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
374 doc="Task to apply the normalization for calibration fluxes (e.g. compensated tophats) "
375 "for the final output star catalog.",
376 )
377 star_apply_aperture_correction = pexConfig.ConfigurableField(
378 target=lsst.meas.base.ApplyApCorrTask,
379 doc="Task to apply aperture corrections to the selected stars."
380 )
381 star_catalog_calculation = pexConfig.ConfigurableField(
382 target=lsst.meas.base.CatalogCalculationTask,
383 doc="Task to compute extendedness values on the star catalog, "
384 "for the star selector to remove extended sources."
385 )
386 star_set_primary_flags = pexConfig.ConfigurableField(
387 target=lsst.meas.algorithms.setPrimaryFlags.SetPrimaryFlagsTask,
388 doc="Task to add isPrimary to the catalog."
389 )
390 star_selector = lsst.meas.algorithms.sourceSelectorRegistry.makeField(
391 default="science",
392 doc="Task to select reliable stars to use for calibration."
393 )
395 # final calibrations and statistics
396 astrometry = pexConfig.ConfigurableField(
397 target=lsst.meas.astrom.AstrometryTask,
398 doc="Task to perform astrometric calibration to fit a WCS.",
399 )
400 astrometry_ref_loader = pexConfig.ConfigField(
401 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
402 doc="Configuration of reference object loader for astrometric fit.",
403 )
404 photometry = pexConfig.ConfigurableField(
405 target=photoCal.PhotoCalTask,
406 doc="Task to perform photometric calibration to fit a PhotoCalib.",
407 )
408 photometry_ref_loader = pexConfig.ConfigField(
409 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
410 doc="Configuration of reference object loader for photometric fit.",
411 )
412 doMaskDiffractionSpikes = pexConfig.Field(
413 dtype=bool,
414 default=False,
415 doc="If True, load the photometric reference catalog again but select "
416 "only bright stars. Use the bright star catalog to set the SPIKE "
417 "mask for regions likely contaminated by diffraction spikes.",
418 )
419 diffractionSpikeMask = pexConfig.ConfigurableField(
420 target=diffractionSpikeMask.DiffractionSpikeMaskTask,
421 doc="Task to identify and mask the diffraction spikes of bright stars.",
422 )
424 compute_summary_stats = pexConfig.ConfigurableField(
425 target=computeExposureSummaryStats.ComputeExposureSummaryStatsTask,
426 doc="Task to to compute summary statistics on the calibrated exposure."
427 )
429 do_illumination_correction = pexConfig.Field(
430 dtype=bool,
431 default=False,
432 doc="If True, apply the illumination correction. This assumes that the "
433 "input image has already been flat-fielded such that it is suitable "
434 "for background subtraction.",
435 )
437 do_calibrate_pixels = pexConfig.Field(
438 dtype=bool,
439 default=True,
440 doc=(
441 "If True, apply the photometric calibration to the image pixels "
442 "and background model, and attach an identity PhotoCalib to "
443 "the output image to reflect this. If False`, leave the image "
444 "and background uncalibrated and attach the PhotoCalib that maps "
445 "them to physical units."
446 )
447 )
449 do_include_astrometric_errors = pexConfig.Field(
450 dtype=bool,
451 default=True,
452 doc="If True, include astrometric errors in the output catalog.",
453 )
455 background_stats_ignored_pixel_masks = pexConfig.ListField(
456 dtype=str,
457 default=["SAT", "SUSPECT", "SPIKE"],
458 doc="Pixel mask flags to ignore when calculating post-sky-subtraction "
459 "background statistics. These are added to those ignored by the "
460 "meas.algorithms.SubtractBackgroundConfig algorithm."
461 )
463 run_sattle = pexConfig.Field(
464 dtype=bool,
465 default=False,
466 doc="If True, the sattle service will populate a cache for later use "
467 "in ip_diffim.detectAndMeasure alert verification."
468 )
470 sattle_historical = pexConfig.Field(
471 dtype=bool,
472 default=False,
473 doc="If re-running a pipeline that requires sattle, this should be set "
474 "to True. This will populate sattle's cache with the historic data "
475 "closest in time to the exposure.",
476 )
478 useButlerCamera = pexConfig.Field(
479 dtype=bool,
480 default=False,
481 doc="If True, use a camera distortion model generated elsewhere in "
482 "the pipeline combined with the telescope boresight as a starting "
483 "point for fitting the WCS, instead of using the WCS attached to "
484 "the exposure, which is generated from the boresight and the "
485 "camera model from the obs_* package."
486 )
488 background_stats_flux_column = pexConfig.Field(
489 dtype=str,
490 default="base_CircularApertureFlux_12_0_flux",
491 doc="Column used to generate post-subtracted background stats."
492 )
494 def setDefaults(self):
495 super().setDefaults()
497 # Use a very broad PSF here, to thoroughly reject CRs.
498 # TODO investigation: a large initial psf guess may make stars look
499 # like CRs for very good seeing images.
500 self.install_simple_psf.fwhm = 4
502 # S/N>=50 sources for PSF determination, but detection to S/N=10.
503 # The thresholdValue sets the minimum flux in a pixel to be included
504 # in the footprint, while peaks are only detected when they are above
505 # thresholdValue * includeThresholdMultiplier. The low thresholdValue
506 # ensures that the footprints are large enough for the noise replacer
507 # to mask out faint undetected neighbors that are not to be measured.
508 self.psf_detection.thresholdValue = 10.0
509 self.psf_detection.includeThresholdMultiplier = 5.0
510 # TODO investigation: Probably want False here, but that may require
511 # tweaking the background spatial scale, to make it small enough to
512 # prevent extra peaks in the wings of bright objects.
513 self.psf_detection.doTempLocalBackground = False
514 self.psf_detection.reEstimateBackground = False
516 # Minimal measurement plugins for PSF determination.
517 # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use
518 # shapeHSM/moments for star/galaxy separation.
519 # TODO DM-39203: we can remove aperture correction from this task
520 # once we are using the shape-based star/galaxy code.
521 self.psf_source_measurement.plugins = ["base_PixelFlags",
522 "base_SdssCentroid",
523 "ext_shapeHSM_HsmSourceMoments",
524 "base_CircularApertureFlux",
525 "base_GaussianFlux",
526 "base_PsfFlux",
527 "base_CompensatedTophatFlux",
528 ]
529 self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
530 # Only measure apertures we need for PSF measurement.
531 self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
532 self.psf_source_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
533 # TODO DM-40843: Remove this line once this is the psfex default.
534 self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \
535 "base_CircularApertureFlux_12_0_instFlux"
537 # No extendedness information available: we need the aperture
538 # corrections to determine that.
539 self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False
540 self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]
541 self.measure_aperture_correction.sourceSelector["science"].flags.bad = []
543 # Detection for good S/N for astrometry/photometry and other
544 # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks.
545 self.star_detection.reEstimateBackground = False
546 self.star_detection.thresholdValue = 5.0
547 self.star_detection.includeThresholdMultiplier = 2.0
548 self.star_measurement.plugins = ["base_PixelFlags",
549 "base_SdssCentroid",
550 "ext_shapeHSM_HsmSourceMoments",
551 "ext_shapeHSM_HsmPsfMoments",
552 "base_GaussianFlux",
553 "base_PsfFlux",
554 "base_CircularApertureFlux",
555 "base_ClassificationSizeExtendedness",
556 "base_CompensatedTophatFlux",
557 ]
558 self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
559 self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
560 # Only measure the apertures we need for star selection.
561 self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
562 self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
564 # We measure and apply the normalization aperture correction with
565 # the psf_normalized_calibration_flux task, and we only apply the
566 # normalization aperture correction for the full list of stars.
567 self.star_normalized_calibration_flux.do_measure_ap_corr = False
569 # Select stars with reliable measurements and no bad flags.
570 self.star_selector["science"].doFlags = True
571 self.star_selector["science"].doUnresolved = True
572 self.star_selector["science"].doSignalToNoise = True
573 self.star_selector["science"].signalToNoise.minimum = 10.0
574 # Keep sky sources in the output catalog, even though they aren't
575 # wanted for calibration.
576 self.star_selector["science"].doSkySources = True
577 # Set the flux and error fields
578 self.star_selector["science"].signalToNoise.fluxField = "slot_CalibFlux_instFlux"
579 self.star_selector["science"].signalToNoise.errField = "slot_CalibFlux_instFluxErr"
581 # Use the affine WCS fitter (assumes we have a good camera geometry).
582 self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)
583 # phot_g_mean is the primary Gaia band for all input bands.
584 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"
586 # Only reject sky sources; we already selected good stars.
587 self.astrometry.sourceSelector["science"].doFlags = True
588 self.astrometry.sourceSelector["science"].flags.good = [] # ["calib_psf_candidate"]
589 # self.astrometry.sourceSelector["science"].flags.bad = []
590 self.astrometry.sourceSelector["science"].doUnresolved = False
591 self.astrometry.sourceSelector["science"].doIsolated = False
592 self.astrometry.sourceSelector["science"].doRequirePrimary = False
593 self.photometry.match.sourceSelection.doFlags = True
594 self.photometry.match.sourceSelection.flags.bad = ["sky_source"]
595 # Unset the (otherwise reasonable, but we've already made the
596 # selections we want above) selection settings in PhotoCalTask.
597 self.photometry.match.sourceSelection.doRequirePrimary = False
598 self.photometry.match.sourceSelection.doUnresolved = False
600 def validate(self):
601 super().validate()
603 # Ensure that the normalization calibration flux tasks
604 # are configured correctly.
605 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:
606 msg = ("psf_normalized_calibration_flux task must be configured with do_measure_ap_corr=True "
607 "or else the normalization and calibration flux will not be properly measured.")
608 raise pexConfig.FieldValidationError(
609 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,
610 )
611 if self.star_normalized_calibration_flux.do_measure_ap_corr:
612 msg = ("star_normalized_calibration_flux task must be configured with do_measure_ap_corr=False "
613 "to apply the previously measured normalization to the full catalog of calibration "
614 "fluxes.")
615 raise pexConfig.FieldValidationError(
616 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,
617 )
619 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,
620 # because they'd be running too early to pick up the fitted PhotoCalib
621 # and WCS.
622 if "base_LocalWcs" in self.psf_source_measurement.plugins.names:
623 raise pexConfig.FieldValidationError(
624 CalibrateImageConfig.psf_source_measurement,
625 self,
626 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
627 )
628 if "base_LocalWcs" in self.star_measurement.plugins.names:
629 raise pexConfig.FieldValidationError(
630 CalibrateImageConfig.star_measurement,
631 self,
632 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
633 )
634 if "base_LocalPhotoCalib" in self.psf_source_measurement.plugins.names:
635 raise pexConfig.FieldValidationError(
636 CalibrateImageConfig.psf_source_measurement,
637 self,
638 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
639 "as it would be run before the photometry fit."
640 )
641 if "base_LocalPhotoCalib" in self.star_measurement.plugins.names:
642 raise pexConfig.FieldValidationError(
643 CalibrateImageConfig.star_measurement,
644 self,
645 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
646 "as it would be run before the photometry fit."
647 )
649 # Check for illumination correction and background consistency.
650 if self.do_illumination_correction:
651 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:
652 raise pexConfig.FieldValidationError(
653 CalibrateImageConfig.psf_subtract_background,
654 self,
655 "CalibrateImageTask.psf_subtract_background must be configured with "
656 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
657 )
658 if self.psf_detection.reEstimateBackground:
659 if not self.psf_detection.doApplyFlatBackgroundRatio:
660 raise pexConfig.FieldValidationError(
661 CalibrateImageConfig.psf_detection,
662 self,
663 "CalibrateImageTask.psf_detection background must be configured with "
664 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
665 )
666 if self.star_detection.reEstimateBackground:
667 if not self.star_detection.doApplyFlatBackgroundRatio:
668 raise pexConfig.FieldValidationError(
669 CalibrateImageConfig.star_detection,
670 self,
671 "CalibrateImageTask.star_detection background must be configured with "
672 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
673 )
674 if not self.star_background.doApplyFlatBackgroundRatio:
675 raise pexConfig.FieldValidationError(
676 CalibrateImageConfig.star_background,
677 self,
678 "CalibrateImageTask.star_background background must be configured with "
679 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
680 )
682 if self.run_sattle:
683 if not os.getenv("SATTLE_URI_BASE"):
684 raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self,
685 "Sattle requested but URI environment variable not set.")
687 if not self.do_adaptive_threshold_detection:
688 if not self.psf_detection.reEstimateBackground:
689 raise pexConfig.FieldValidationError(
690 CalibrateImageConfig.psf_detection,
691 self,
692 "If not using the adaptive threshold detection implementation (i.e. "
693 "config.do_adaptive_threshold_detection = False), CalibrateImageTask.psf_detection "
694 "background must be configured with "
695 "reEstimateBackground = True to maintain original behavior."
696 )
697 if not self.star_detection.reEstimateBackground:
698 raise pexConfig.FieldValidationError(
699 CalibrateImageConfig.psf_detection,
700 self,
701 "If not using the adaptive threshold detection implementation "
702 "(i.e. config.do_adaptive_threshold_detection = False), "
703 "CalibrateImageTask.star_detection background must be configured "
704 "with reEstimateBackground = True to maintain original behavior."
705 )
708class CalibrateImageTask(pipeBase.PipelineTask):
709 """Compute the PSF, aperture corrections, astrometric and photometric
710 calibrations, and summary statistics for a single science exposure, and
711 produce a catalog of brighter stars that were used to calibrate it.
713 Parameters
714 ----------
715 initial_stars_schema : `lsst.afw.table.Schema`
716 Schema of the initial_stars output catalog.
717 """
718 _DefaultName = "calibrateImage"
719 ConfigClass = CalibrateImageConfig
721 def __init__(self, initial_stars_schema=None, **kwargs):
722 super().__init__(**kwargs)
723 self.makeSubtask("snap_combine")
725 # PSF determination subtasks
726 self.makeSubtask("install_simple_psf")
727 self.makeSubtask("psf_repair")
728 self.makeSubtask("psf_subtract_background")
729 self.psf_schema = afwTable.SourceTable.makeMinimalSchema()
730 self.psf_schema.addField(
731 'psf_max_value',
732 type=np.float32,
733 doc="PSF max value.",
734 )
735 afwTable.CoordKey.addErrorFields(self.psf_schema)
736 self.makeSubtask("psf_detection", schema=self.psf_schema)
737 self.makeSubtask("psf_source_measurement", schema=self.psf_schema)
738 self.makeSubtask("psf_adaptive_threshold_detection")
739 self.makeSubtask("psf_measure_psf", schema=self.psf_schema)
740 self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema)
742 self.makeSubtask("measure_aperture_correction", schema=self.psf_schema)
743 self.makeSubtask("astrometry", schema=self.psf_schema)
745 # star measurement subtasks
746 if initial_stars_schema is None:
747 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()
749 # These fields let us track which sources were used for psf modeling,
750 # astrometric fitting, and aperture correction calculations.
751 self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved",
752 "calib_astrometry_used",
753 # TODO DM-39203:
754 # these can be removed once apcorr is gone.
755 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",
756 "apcorr_base_PsfFlux_used",)
757 for field in self.psf_fields:
758 item = self.psf_schema.find(field)
759 initial_stars_schema.addField(item.getField())
760 id_type = self.psf_schema["id"].asField().getTypeString()
761 psf_max_value_type = self.psf_schema['psf_max_value'].asField().getTypeString()
762 initial_stars_schema.addField("psf_id",
763 type=id_type,
764 doc="id of this source in psf_stars; 0 if there is no match.")
765 initial_stars_schema.addField("psf_max_value",
766 type=psf_max_value_type,
767 doc="Maximum value in the star image used to train PSF.")
769 afwTable.CoordKey.addErrorFields(initial_stars_schema)
770 self.makeSubtask("star_background")
771 self.makeSubtask("star_detection", schema=initial_stars_schema)
772 self.makeSubtask("star_sky_sources", schema=initial_stars_schema)
773 self.makeSubtask("star_deblend", schema=initial_stars_schema)
774 self.makeSubtask("star_measurement", schema=initial_stars_schema)
775 self.makeSubtask("star_normalized_calibration_flux", schema=initial_stars_schema)
777 self.makeSubtask("star_apply_aperture_correction", schema=initial_stars_schema)
778 self.makeSubtask("star_catalog_calculation", schema=initial_stars_schema)
779 self.makeSubtask("star_set_primary_flags", schema=initial_stars_schema, isSingleFrame=True)
780 self.makeSubtask("star_selector")
781 self.makeSubtask("photometry", schema=initial_stars_schema)
782 if self.config.doMaskDiffractionSpikes:
783 self.makeSubtask("diffractionSpikeMask")
784 self.makeSubtask("compute_summary_stats")
786 # The final catalog will have calibrated flux columns, which we add to
787 # the init-output schema by calibrating our zero-length catalog with an
788 # arbitrary dummy PhotoCalib. We also use this schema to initialize
789 # the stars catalog in order to ensure it's the same even when we hit
790 # an error (and write partial outputs) before calibrating the catalog
791 # - note that calibrateCatalog will happily reuse existing output
792 # columns.
793 dummy_photo_calib = afwImage.PhotoCalib(1.0, 0, bbox=lsst.geom.Box2I())
794 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(
795 afwTable.SourceCatalog(initial_stars_schema)
796 )
798 def runQuantum(self, butlerQC, inputRefs, outputRefs):
799 inputs = butlerQC.get(inputRefs)
800 exposures = inputs.pop("exposures")
802 id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)
804 astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
805 dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],
806 refCats=inputs.pop("astrometry_ref_cat"),
807 name=self.config.connections.astrometry_ref_cat,
808 config=self.config.astrometry_ref_loader, log=self.log)
809 self.astrometry.setRefObjLoader(astrometry_loader)
811 photometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
812 dataIds=[ref.datasetRef.dataId for ref in inputRefs.photometry_ref_cat],
813 refCats=inputs.pop("photometry_ref_cat"),
814 name=self.config.connections.photometry_ref_cat,
815 config=self.config.photometry_ref_loader, log=self.log)
816 self.photometry.match.setRefObjLoader(photometry_loader)
818 if self.config.doMaskDiffractionSpikes:
819 # Use the same photometry reference catalog for the bright star
820 # mask.
821 self.diffractionSpikeMask.setRefObjLoader(photometry_loader)
823 if self.config.do_illumination_correction:
824 background_flat = inputs.pop("background_flat")
825 illumination_correction = inputs.pop("illumination_correction")
826 else:
827 background_flat = None
828 illumination_correction = None
830 camera_model = None
831 if self.config.useButlerCamera:
832 if "camera_model" in inputs:
833 camera_model = inputs.pop("camera_model")
834 else:
835 self.log.warning("useButlerCamera=True, but camera is not available for filter %s. The "
836 "astrometry fit will use the WCS already attached to the exposure.",
837 exposures[0].filter.bandLabel)
839 # This should not happen with a properly configured execution context.
840 assert not inputs, "runQuantum got more inputs than expected"
842 # Specify the fields that `annotate` needs below, to ensure they
843 # exist, even as None.
844 result = pipeBase.Struct(
845 exposure=None,
846 stars_footprints=None,
847 psf_stars_footprints=None,
848 background_to_photometric_ratio=None,
849 )
850 try:
851 self.run(
852 exposures=exposures,
853 result=result,
854 id_generator=id_generator,
855 background_flat=background_flat,
856 illumination_correction=illumination_correction,
857 camera_model=camera_model,
858 )
859 except pipeBase.AlgorithmError as e:
860 error = pipeBase.AnnotatedPartialOutputsError.annotate(
861 e,
862 self,
863 result.exposure,
864 result.psf_stars_footprints,
865 result.stars_footprints,
866 log=self.log
867 )
868 butlerQC.put(result, outputRefs)
869 raise error from e
871 butlerQC.put(result, outputRefs)
873 @timeMethod
874 def run(
875 self,
876 *,
877 exposures,
878 id_generator=None,
879 result=None,
880 background_flat=None,
881 illumination_correction=None,
882 camera_model=None,
883 ):
884 """Find stars and perform psf measurement, then do a deeper detection
885 and measurement and calibrate astrometry and photometry from that.
887 Parameters
888 ----------
889 exposures : `lsst.afw.image.Exposure` or \
890 `list` [`lsst.afw.image.Exposure`]
891 Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter.
892 Modified in-place during processing if only one is passed.
893 If two exposures are passed, treat them as snaps and combine
894 before doing further processing.
895 id_generator : `lsst.meas.base.IdGenerator`, optional
896 Object that generates source IDs and provides random seeds.
897 result : `lsst.pipe.base.Struct`, optional
898 Result struct that is modified to allow saving of partial outputs
899 for some failure conditions. If the task completes successfully,
900 this is also returned.
901 background_flat : `lsst.afw.image.Exposure`, optional
902 Background flat-field image.
903 illumination_correction : `lsst.afw.image.Exposure`, optional
904 Illumination correction image.
905 camera_model : `lsst.afw.cameraGeom.Camera`, optional
906 Camera to be used if constructing updated WCS.
908 Returns
909 -------
910 result : `lsst.pipe.base.Struct`
911 Results as a struct with attributes:
913 ``exposure``
914 Calibrated exposure, with pixels in nJy units.
915 (`lsst.afw.image.Exposure`)
916 ``stars``
917 Stars that were used to calibrate the exposure, with
918 calibrated fluxes and magnitudes.
919 (`astropy.table.Table`)
920 ``stars_footprints``
921 Footprints of stars that were used to calibrate the exposure.
922 (`lsst.afw.table.SourceCatalog`)
923 ``psf_stars``
924 Stars that were used to determine the image PSF.
925 (`astropy.table.Table`)
926 ``psf_stars_footprints``
927 Footprints of stars that were used to determine the image PSF.
928 (`lsst.afw.table.SourceCatalog`)
929 ``background``
930 Background that was fit to the exposure when detecting
931 ``stars``. (`lsst.afw.math.BackgroundList`)
932 ``applied_photo_calib``
933 Photometric calibration that was fit to the star catalog and
934 applied to the exposure. (`lsst.afw.image.PhotoCalib`)
935 This is `None` if ``config.do_calibrate_pixels`` is `False`.
936 ``astrometry_matches``
937 Reference catalog stars matches used in the astrometric fit.
938 (`list` [`lsst.afw.table.ReferenceMatch`] or
939 `lsst.afw.table.BaseCatalog`).
940 ``photometry_matches``
941 Reference catalog stars matches used in the photometric fit.
942 (`list` [`lsst.afw.table.ReferenceMatch`] or
943 `lsst.afw.table.BaseCatalog`).
944 ``mask``
945 Copy of the mask plane of `exposure`.
946 (`lsst.afw.image.Mask`)
947 """
948 if result is None:
949 result = pipeBase.Struct()
950 if id_generator is None:
951 id_generator = lsst.meas.base.IdGenerator()
953 result.exposure = self.snap_combine.run(exposures).exposure
954 self.log.info("Initial PhotoCalib: %s", result.exposure.getPhotoCalib())
956 result.exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = False
958 # Check input image processing.
959 if self.config.do_illumination_correction:
960 if not result.exposure.metadata.get("LSST ISR FLAT APPLIED", False):
961 raise pipeBase.InvalidQuantumError(
962 "Cannot use do_illumination_correction with an image that has not had a flat applied",
963 )
965 if camera_model:
966 if camera_model.get(result.exposure.detector.getId()):
967 self.log.info("Updating WCS with the provided camera model.")
968 self._update_wcs_with_camera_model(result.exposure, camera_model)
969 else:
970 self.log.warning(
971 "useButlerCamera=True, but detector %s is not available in the provided camera. The "
972 "astrometry fit will use the WCS already attached to the exposure.",
973 result.exposure.detector.getId())
975 result.background = None
976 result.background_to_photometric_ratio = None
977 summary_stat_catalog = None
978 # Some exposure components are set to initial placeholder objects
979 # while we try to bootstrap them. If we fail before we fit for them,
980 # we want to reset those components to None so the placeholders don't
981 # masquerade as the real thing.
982 have_fit_psf = False
983 have_fit_astrometry = False
984 have_fit_photometry = False
985 try:
986 result.background_to_photometric_ratio = self._apply_illumination_correction(
987 result.exposure,
988 background_flat,
989 illumination_correction,
990 )
992 result.psf_stars_footprints, _, adaptive_det_res_struct = self._compute_psf(
993 result,
994 id_generator,
995 )
996 have_fit_psf = True
998 # Check if all centroids have been flagged. This should happen
999 # after the PSF is fit and recorded because even a junky PSF
1000 # model is informative.
1001 if result.psf_stars_footprints["slot_Centroid_flag"].all():
1002 psf_shape = result.exposure.psf.computeShape(result.exposure.psf.getAveragePosition())
1003 raise AllCentroidsFlaggedError(
1004 n_sources=len(result.psf_stars_footprints),
1005 psf_shape_ixx=psf_shape.getIxx(),
1006 psf_shape_iyy=psf_shape.getIyy(),
1007 psf_shape_ixy=psf_shape.getIxy(),
1008 psf_size=psf_shape.getDeterminantRadius(),
1009 )
1011 if result.background is None:
1012 result.background = afwMath.BackgroundList()
1014 self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)
1015 result.psf_stars = result.psf_stars_footprints.asAstropy()
1016 # Run astrometry using PSF candidate stars.
1017 # Update "the psf_stars" source coordinates with the current wcs.
1018 afwTable.updateSourceCoords(
1019 result.exposure.wcs,
1020 sourceList=result.psf_stars_footprints,
1021 include_covariance=self.config.do_include_astrometric_errors,
1022 )
1023 astrometry_matches, astrometry_meta = self._fit_astrometry(
1024 result.exposure, result.psf_stars_footprints
1025 )
1026 num_astrometry_matches = len(astrometry_matches)
1027 if "astrometry_matches" in self.config.optional_outputs:
1028 result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches,
1029 astrometry_meta)
1030 result.psf_stars = result.psf_stars_footprints.asAstropy()
1032 # Add diffraction spike mask here so it can be used (i.e. avoided)
1033 # in the adaptive background estimation.
1034 if self.config.doMaskDiffractionSpikes:
1035 self.diffractionSpikeMask.run(result.exposure)
1037 self.metadata['adaptive_threshold_value'] = float("nan")
1038 if self.config.do_adaptive_threshold_detection:
1039 self._remeasure_star_background(
1040 result,
1041 background_to_photometric_ratio=result.background_to_photometric_ratio,
1042 )
1044 # Run the stars_detection subtask for the photometric calibration.
1045 result.stars_footprints = self._find_stars(
1046 result.exposure,
1047 result.background,
1048 id_generator,
1049 background_to_photometric_ratio=result.background_to_photometric_ratio,
1050 adaptive_det_res_struct=adaptive_det_res_struct,
1051 num_astrometry_matches=num_astrometry_matches,
1052 )
1053 psf = result.exposure.getPsf()
1054 psfSigma = psf.computeShape(result.exposure.getBBox().getCenter()).getDeterminantRadius()
1055 self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints,
1056 psfSigma=psfSigma)
1058 # Update the "stars" source coordinates with the current wcs.
1059 afwTable.updateSourceCoords(
1060 result.exposure.wcs,
1061 sourceList=result.stars_footprints,
1062 include_covariance=self.config.do_include_astrometric_errors,
1063 )
1065 summary_stat_catalog = result.stars_footprints
1066 result.stars = result.stars_footprints.asAstropy()
1067 self.metadata["star_count"] = np.sum(~result.stars["sky_source"])
1069 # Validate the astrometric fit. Send in the stars_footprints
1070 # catalog so that its coords get set to NaN if the fit is deemed
1071 # a failure.
1072 self.astrometry.check(result.exposure, result.stars_footprints, len(astrometry_matches))
1073 result.stars = result.stars_footprints.asAstropy()
1074 have_fit_astrometry = True
1076 result.stars_footprints, photometry_matches, \
1077 photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints)
1078 have_fit_photometry = True
1079 self.metadata["photometry_matches_count"] = len(photometry_matches)
1080 # fit_photometry returns a new catalog, so we need a new astropy
1081 # table view.
1082 result.stars = result.stars_footprints.asAstropy()
1083 # summary stats don't make use of the calibrated fluxes, but we
1084 # might as well use the best catalog we've got in case that
1085 # changes, and help the old one get garbage-collected.
1086 summary_stat_catalog = result.stars_footprints
1087 if "photometry_matches" in self.config.optional_outputs:
1088 result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,
1089 photometry_meta)
1090 if "mask" in self.config.optional_outputs:
1091 result.mask = result.exposure.mask.clone()
1092 except pipeBase.AlgorithmError:
1093 if not have_fit_psf:
1094 result.exposure.setPsf(None)
1095 if not have_fit_astrometry:
1096 result.exposure.setWcs(None)
1097 if not have_fit_photometry:
1098 result.exposure.setPhotoCalib(None)
1099 # Summary stat calculations can handle missing components
1100 # gracefully, but we want to run them as late as possible (but
1101 # still before we calibrate pixels, if we do that at all).
1102 # So we run them after we succeed or if we get an AlgorithmError.
1103 # We intentionally don't use 'finally' here because we don't
1104 # want to run them if we get some other kind of error.
1105 self._summarize(result.exposure, summary_stat_catalog, result.background)
1106 raise
1107 else:
1108 self._summarize(result.exposure, summary_stat_catalog, result.background)
1110 # Output post-subtraction background stats to task metadata:
1111 # Specify the pixels flags to ignore, starting with those ignored
1112 # by the subtraction.
1113 bgIgnoredPixelMasks = lsst.meas.algorithms.SubtractBackgroundConfig().ignoredPixelMask.list()
1114 for maskName in self.config.background_stats_ignored_pixel_masks:
1115 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
1116 and maskName not in bgIgnoredPixelMasks):
1117 bgIgnoredPixelMasks += [maskName]
1119 statsCtrl = afwMath.StatisticsControl()
1120 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(bgIgnoredPixelMasks))
1122 statObj = afwMath.makeStatistics(result.exposure.getMaskedImage(), afwMath.MEDIAN, statsCtrl)
1123 median_bg, _ = statObj.getResult(afwMath.MEDIAN)
1124 self.metadata["bg_subtracted_skyPixel_instFlux_median"] = median_bg
1126 statObj = afwMath.makeStatistics(result.exposure.getMaskedImage(), afwMath.STDEV, statsCtrl)
1127 stdev_bg, _ = statObj.getResult(afwMath.STDEV)
1128 self.metadata["bg_subtracted_skyPixel_instFlux_stdev"] = stdev_bg
1130 self.metadata["bg_subtracted_skySource_flux_median"] = (
1131 np.median(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column])
1132 )
1133 self.metadata["bg_subtracted_skySource_flux_stdev"] = (
1134 np.std(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column])
1135 )
1137 if self.config.do_calibrate_pixels:
1138 self._apply_photometry(
1139 result.exposure,
1140 result.background,
1141 background_to_photometric_ratio=result.background_to_photometric_ratio,
1142 )
1143 result.applied_photo_calib = photo_calib
1144 else:
1145 result.applied_photo_calib = None
1147 self._recordMaskedPixelFractions(result.exposure)
1149 if self.config.run_sattle:
1150 # send boresight and timing information to sattle so the cache
1151 # is populated by the time we reach ip_diffim detectAndMeasure.
1152 try:
1153 populate_sattle_visit_cache(result.exposure.getInfo().getVisitInfo(),
1154 historical=self.config.sattle_historical)
1155 self.log.info('Successfully triggered load of sattle visit cache')
1156 except requests.exceptions.HTTPError:
1157 self.log.exception("Sattle visit cache update failed; continuing with image processing")
1159 return result
1161 def _apply_illumination_correction(self, exposure, background_flat, illumination_correction):
1162 """Apply the illumination correction to a background-flattened image.
1164 Parameters
1165 ----------
1166 exposure : `lsst.afw.image.Exposure`
1167 Exposure to convert to a photometric-flattened image.
1168 background_flat : `lsst.afw.image.Exposure`
1169 Flat image that had previously been applied to exposure.
1170 illumination_correction : `lsst.afw.image.Exposure`
1171 Illumination correction image to convert to photometric-flattened
1172 image.
1174 Returns
1175 -------
1176 background_to_photometric_ratio : `lsst.afw.image.Image`
1177 Ratio image to convert a photometric-flattened image to/from
1178 a background-flattened image. Will be None if task not
1179 configured to use the illumination correction.
1180 """
1181 if not self.config.do_illumination_correction:
1182 return None
1184 # From a raw image to a background-flattened image, we have:
1185 # bfi = image / background_flat
1186 # From a raw image to a photometric-flattened image, we have:
1187 # pfi = image / reference_flux_flat
1188 # pfi = image / (dome_flat * illumination_correction),
1189 # where the illumination correction contains the jacobian
1190 # of the wcs, converting to fluence units.
1191 # Currently background_flat == dome_flat, so we have for the
1192 # "background_to_photometric_ratio", the ratio of the background-
1193 # flattened image to the photometric-flattened image:
1194 # bfi / pfi = illumination_correction.
1196 background_to_photometric_ratio = illumination_correction.image.clone()
1198 # Dividing the ratio will convert a background-flattened image to
1199 # a photometric-flattened image.
1200 exposure.maskedImage /= background_to_photometric_ratio
1202 exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = True
1204 return background_to_photometric_ratio
1206 def _compute_psf(self, result, id_generator):
1207 """Find bright sources detected on an exposure and fit a PSF model to
1208 them, repairing likely cosmic rays before detection.
1210 Repair, detect, measure, and compute PSF twice, to ensure the PSF
1211 model does not include contributions from cosmic rays.
1213 Parameters
1214 ----------
1215 result : `lsst.pipe.base.Struct`
1216 Result struct that is modified to allow saving of partial outputs
1217 for some failure conditions. Should contain at least the following
1218 attributes:
1220 - exposure : `lsst.afw.image.Exposure`
1221 Exposure to detect and measure bright stars on.
1222 - background : `lsst.afw.math.BackgroundList` | `None`
1223 Background that was fit to the exposure during detection.
1224 - background_to_photometric_ratio : `lsst.afw.image.Image` | `None`
1225 Image to convert photometric-flattened image to
1226 background-flattened image.
1227 id_generator : `lsst.meas.base.IdGenerator`
1228 Object that generates source IDs and provides random seeds.
1230 Returns
1231 -------
1232 sources : `lsst.afw.table.SourceCatalog`
1233 Catalog of detected bright sources.
1234 cell_set : `lsst.afw.math.SpatialCellSet`
1235 PSF candidates returned by the psf determiner.
1236 adaptive_det_res_struct : `lsst.pipe.base.Struct`
1237 Result struct from the adaptive threshold detection.
1239 Notes
1240 -----
1241 This method modifies the exposure, background and
1242 background_to_photometric_ratio attributes of the result struct
1243 in-place.
1244 """
1245 def log_psf(msg, addToMetadata=False):
1246 """Log the parameters of the psf and background, with a prepended
1247 message. There is also the option to add the PSF sigma to the task
1248 metadata.
1250 Parameters
1251 ----------
1252 msg : `str`
1253 Message to prepend the log info with.
1254 addToMetadata : `bool`, optional
1255 Whether to add the final psf sigma value to the task
1256 metadata (the default is False).
1257 """
1258 position = result.exposure.psf.getAveragePosition()
1259 sigma = result.exposure.psf.computeShape(position).getDeterminantRadius()
1260 dimensions = result.exposure.psf.computeImage(position).getDimensions()
1261 if result.background is not None:
1262 median_background = np.median(result.background.getImage().array)
1263 else:
1264 median_background = 0.0
1265 self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",
1266 msg, sigma, dimensions, median_background)
1267 if addToMetadata:
1268 self.metadata["final_psf_sigma"] = sigma
1270 self.log.info("First pass detection with Gaussian PSF FWHM=%s pixels",
1271 self.config.install_simple_psf.fwhm)
1272 self.install_simple_psf.run(exposure=result.exposure)
1274 result.background = self.psf_subtract_background.run(
1275 exposure=result.exposure,
1276 backgroundToPhotometricRatio=result.background_to_photometric_ratio,
1277 ).background
1278 log_psf("Initial PSF:")
1279 self.psf_repair.run(exposure=result.exposure, keepCRs=True)
1281 table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())
1282 if not self.config.do_adaptive_threshold_detection:
1283 adaptive_det_res_struct = None
1284 # Re-estimate the background during this detection step, so that
1285 # measurement uses the most accurate background-subtraction.
1286 detections = self.psf_detection.run(
1287 table=table,
1288 exposure=result.exposure,
1289 background=result.background,
1290 backgroundToPhotometricRatio=result.background_to_photometric_ratio,
1291 )
1292 else:
1293 initialThreshold = self.config.psf_detection.thresholdValue
1294 initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier
1295 adaptive_det_res_struct = self.psf_adaptive_threshold_detection.run(
1296 table, result.exposure,
1297 initialThreshold=initialThreshold,
1298 initialThresholdMultiplier=initialThresholdMultiplier,
1299 )
1300 detections = adaptive_det_res_struct.detections
1302 self.metadata["initial_psf_positive_footprint_count"] = detections.numPos
1303 self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg
1304 self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks
1305 self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks
1306 self.psf_source_measurement.run(detections.sources, result.exposure)
1307 psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)
1309 # This 2nd round of PSF fitting has been deemed unnecessary (and
1310 # sometimes even causing harm), so it is being skipped for the
1311 # adaptive threshold logic, but is left here for now to maintain
1312 # consistency with the previous logic for any users wanting to
1313 # revert to it.
1314 if not self.config.do_adaptive_threshold_detection:
1315 # Replace the initial PSF with something simpler for the second
1316 # repair/detect/measure/measure_psf step: this can help it
1317 # converge.
1318 self.install_simple_psf.run(exposure=result.exposure)
1320 log_psf("Rerunning with simple PSF:")
1321 # TODO investigation: Should we only re-run repair here, to use the
1322 # new PSF? Maybe we *do* need to re-run measurement with PsfFlux,
1323 # to use the fitted PSF?
1324 # TODO investigation: do we need a separate measurement task here
1325 # for the post-psf_measure_psf step, since we only want to do
1326 # PsfFlux and GaussianFlux *after* we have a PSF? Maybe that's not
1327 # relevant once DM-39203 is merged?
1328 self.psf_repair.run(exposure=result.exposure, keepCRs=True)
1329 # Re-estimate the background during this detection step, so that
1330 # measurement uses the most accurate background-subtraction.
1331 detections = self.psf_detection.run(
1332 table=table,
1333 exposure=result.exposure,
1334 background=result.background,
1335 backgroundToPhotometricRatio=result.background_to_photometric_ratio,
1336 )
1337 self.psf_source_measurement.run(detections.sources, result.exposure)
1338 psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)
1339 self.metadata["simple_psf_positive_footprint_count"] = detections.numPos
1340 self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg
1341 self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks
1342 self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks
1343 log_psf("Final PSF:", addToMetadata=True)
1345 # Final repair with final PSF, removing cosmic rays this time.
1346 self.psf_repair.run(exposure=result.exposure)
1347 # Final measurement with the CRs removed.
1348 self.psf_source_measurement.run(detections.sources, result.exposure)
1350 # PSF is set on exposure; candidates are returned to use for
1351 # calibration flux normalization and aperture corrections.
1352 return detections.sources, psf_result.cellSet, adaptive_det_res_struct
1354 def _measure_aperture_correction(self, exposure, bright_sources):
1355 """Measure and set the ApCorrMap on the Exposure, using
1356 previously-measured bright sources.
1358 This function first normalizes the calibration flux and then
1359 the full set of aperture corrections are measured relative
1360 to this normalized calibration flux.
1362 Parameters
1363 ----------
1364 exposure : `lsst.afw.image.Exposure`
1365 Exposure to set the ApCorrMap on.
1366 bright_sources : `lsst.afw.table.SourceCatalog`
1367 Catalog of detected bright sources; modified to include columns
1368 necessary for point source determination for the aperture
1369 correction calculation.
1370 """
1371 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(
1372 exposure=exposure,
1373 catalog=bright_sources,
1374 ).ap_corr_map
1376 ap_corr_map = self.measure_aperture_correction.run(exposure, bright_sources).apCorrMap
1378 # Need to merge the aperture correction map from the normalization.
1379 for key in norm_ap_corr_map:
1380 ap_corr_map[key] = norm_ap_corr_map[key]
1382 exposure.info.setApCorrMap(ap_corr_map)
1384 def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None,
1385 adaptive_det_res_struct=None, num_astrometry_matches=None):
1386 """Detect stars on an exposure that has a PSF model, and measure their
1387 PSF, circular aperture, compensated gaussian fluxes.
1389 Parameters
1390 ----------
1391 exposure : `lsst.afw.image.Exposure`
1392 Exposure to detect and measure stars on.
1393 background : `lsst.afw.math.BackgroundList`
1394 Background that was fit to the exposure during detection;
1395 modified in-place during subsequent detection.
1396 id_generator : `lsst.meas.base.IdGenerator`
1397 Object that generates source IDs and provides random seeds.
1398 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1399 Image to convert photometric-flattened image to
1400 background-flattened image.
1402 Returns
1403 -------
1404 stars : `SourceCatalog`
1405 Sources that are very likely to be stars, with a limited set of
1406 measurements performed on them.
1407 """
1408 table = afwTable.SourceTable.make(self.initial_stars_schema.schema,
1409 id_generator.make_table_id_factory())
1411 maxAdaptiveDetIter = 8
1412 threshFactor = 0.2
1413 if adaptive_det_res_struct is not None:
1414 for adaptiveDetIter in range(maxAdaptiveDetIter):
1415 adaptiveDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
1416 adaptiveDetectionConfig.reEstimateBackground = False
1417 adaptiveDetectionConfig.includeThresholdMultiplier = 2.0
1418 psfThreshold = (
1419 adaptive_det_res_struct.thresholdValue*adaptive_det_res_struct.includeThresholdMultiplier
1420 )
1421 adaptiveDetectionConfig.thresholdValue = max(
1422 self.config.star_detection.thresholdValue,
1423 threshFactor*psfThreshold/adaptiveDetectionConfig.includeThresholdMultiplier
1424 )
1425 self.log.info("Using adaptive threshold detection (nIter = %d) with "
1426 "thresholdValue = %.2f and multiplier = %.1f",
1427 adaptiveDetIter, adaptiveDetectionConfig.thresholdValue,
1428 adaptiveDetectionConfig.includeThresholdMultiplier)
1429 adaptiveDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
1430 config=adaptiveDetectionConfig
1431 )
1432 detections = adaptiveDetectionTask.run(
1433 table=table,
1434 exposure=exposure,
1435 background=background,
1436 backgroundToPhotometricRatio=background_to_photometric_ratio,
1437 )
1438 nFootprint = len(detections.sources)
1439 nPeak = 0
1440 nIsolated = 0
1441 for src in detections.sources:
1442 nPeakSrc = len(src.getFootprint().getPeaks())
1443 if nPeakSrc == 1:
1444 nIsolated += 1
1445 nPeak += nPeakSrc
1446 minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches))
1447 if nFootprint > 0:
1448 self.log.info("Adaptive threshold detection _find_stars nIter: %d; "
1449 "nPeak/nFootprint = %.2f (max is 800); nIsolated = %d (min is %.1f).",
1450 adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated)
1451 if nPeak/nFootprint > 800 or nIsolated < minIsolated:
1452 threshFactor = max(0.01, 1.5*threshFactor)
1453 self.log.warning("nPeak/nFootprint = %.2f (max is 800); nIsolated = %d "
1454 "(min is %.1f).", nPeak/nFootprint, nIsolated, minIsolated)
1455 else:
1456 break
1457 else:
1458 threshFactor *= 0.75
1459 self.log.warning("No footprints detected on image. Decreasing threshold "
1460 "factor to %.2f. and rerunning.", threshFactor)
1461 else:
1462 # Re-estimate the background during this detection step, so that
1463 # measurement uses the most accurate background-subtraction.
1464 detections = self.star_detection.run(
1465 table=table,
1466 exposure=exposure,
1467 background=background,
1468 backgroundToPhotometricRatio=background_to_photometric_ratio,
1469 )
1470 sources = detections.sources
1471 self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources)
1473 n_sky_sources = np.sum(sources["sky_source"])
1474 if (self.config.do_downsample_footprints
1475 and (len(sources) - n_sky_sources) > self.config.downsample_max_footprints):
1476 if exposure.info.id is None:
1477 self.log.warning("Exposure does not have a proper id; using 0 seed for downsample.")
1478 seed = 0
1479 else:
1480 seed = exposure.info.id & 0xFFFFFFFF
1482 gen = np.random.RandomState(seed)
1484 # Isolate the sky sources before downsampling.
1485 indices = np.arange(len(sources))[~sources["sky_source"]]
1486 indices = gen.choice(
1487 indices,
1488 size=self.config.downsample_max_footprints,
1489 replace=False,
1490 )
1491 skyIndices, = np.where(sources["sky_source"])
1492 indices = np.concatenate((indices, skyIndices))
1494 self.log.info("Downsampling from %d to %d non-sky-source footprints.", len(sources), len(indices))
1496 sel = np.zeros(len(sources), dtype=bool)
1497 sel[indices] = True
1498 sources = sources[sel]
1500 # TODO investigation: Could this deblender throw away blends of
1501 # non-PSF sources?
1502 self.star_deblend.run(exposure=exposure, sources=sources)
1503 # The deblender may not produce a contiguous catalog; ensure
1504 # contiguity for subsequent tasks.
1505 if not sources.isContiguous():
1506 sources = sources.copy(deep=True)
1508 # Measure everything, and use those results to select only stars.
1509 self.star_measurement.run(sources, exposure)
1510 self.metadata["post_deblend_source_count"] = np.sum(~sources["sky_source"])
1511 self.metadata["failed_deblend_source_count"] = np.sum(
1512 ~sources["sky_source"] & sources["deblend_failed"]
1513 )
1514 self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"])
1515 self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"])
1517 # Run the normalization calibration flux task to apply the
1518 # normalization correction to create normalized
1519 # calibration fluxes.
1520 self.star_normalized_calibration_flux.run(exposure=exposure, catalog=sources)
1521 self.star_apply_aperture_correction.run(sources, exposure.apCorrMap)
1522 self.star_catalog_calculation.run(sources)
1523 self.star_set_primary_flags.run(sources)
1525 result = self.star_selector.run(sources)
1526 # The star selector may not produce a contiguous catalog.
1527 if not result.sourceCat.isContiguous():
1528 return result.sourceCat.copy(deep=True)
1529 else:
1530 return result.sourceCat
1532 def _match_psf_stars(self, psf_stars, stars, psfSigma=None):
1533 """Match calibration stars to psf stars, to identify which were psf
1534 candidates, and which were used or reserved during psf measurement
1535 and the astrometric fit.
1537 Parameters
1538 ----------
1539 psf_stars : `lsst.afw.table.SourceCatalog`
1540 PSF candidate stars that were sent to the psf determiner and
1541 used in the astrometric fit. Used to populate psf and astrometry
1542 related flag fields.
1543 stars : `lsst.afw.table.SourceCatalog`
1544 Stars that will be used for calibration; psf-related fields will
1545 be updated in-place.
1547 Notes
1548 -----
1549 This code was adapted from CalibrateTask.copyIcSourceFields().
1550 """
1551 control = afwTable.MatchControl()
1552 matchRadius = 3.0 if psfSigma is None else max(3.0, psfSigma) # in pixels
1553 # Return all matched objects, to separate blends.
1554 control.findOnlyClosest = False
1555 matches = afwTable.matchXy(psf_stars, stars, matchRadius, control)
1556 deblend_key = stars.schema["deblend_nChild"].asKey()
1557 matches = [m for m in matches if m[1].get(deblend_key) == 0]
1559 # Because we had to allow multiple matches to handle parents, we now
1560 # need to prune to the best (closest) matches.
1561 # Closest matches is a dict of psf_stars source ID to Match record
1562 # (psf_stars source, sourceCat source, distance in pixels).
1563 best = {}
1564 for match_psf, match_stars, d in matches:
1565 match = best.get(match_psf.getId())
1566 if match is None or d <= match[2]:
1567 best[match_psf.getId()] = (match_psf, match_stars, d)
1568 matches = list(best.values())
1569 # We'll use this to construct index arrays into each catalog.
1570 ids = np.array([(match_psf.getId(), match_stars.getId()) for match_psf, match_stars, d in matches]).T
1572 if (n_matches := len(matches)) == 0:
1573 raise NoPsfStarsToStarsMatchError(n_psf_stars=len(psf_stars), n_stars=len(stars))
1575 self.log.info("%d psf/astrometry stars out of %d matched %d calib stars",
1576 n_matches, len(psf_stars), len(stars))
1577 self.metadata["matched_psf_star_count"] = n_matches
1579 # Check that no stars sources are listed twice; we already know
1580 # that each match has a unique psf_stars id, due to using as the key
1581 # in best above.
1582 n_unique = len(set(m[1].getId() for m in matches))
1583 if n_unique != n_matches:
1584 self.log.warning("%d psf_stars matched only %d stars", n_matches, n_unique)
1586 # The indices of the IDs, so we can update the flag fields as arrays.
1587 idx_psf_stars = np.searchsorted(psf_stars["id"], ids[0])
1588 idx_stars = np.searchsorted(stars["id"], ids[1])
1589 for field in self.psf_fields:
1590 result = np.zeros(len(stars), dtype=bool)
1591 result[idx_stars] = psf_stars[field][idx_psf_stars]
1592 stars[field] = result
1593 stars['psf_id'][idx_stars] = psf_stars['id'][idx_psf_stars]
1594 stars['psf_max_value'][idx_stars] = psf_stars['psf_max_value'][idx_psf_stars]
1596 def _fit_astrometry(self, exposure, stars):
1597 """Fit an astrometric model to the data and return the reference
1598 matches used in the fit, and the fitted WCS.
1600 Parameters
1601 ----------
1602 exposure : `lsst.afw.image.Exposure`
1603 Exposure that is being fit, to get PSF and other metadata from.
1604 Modified to add the fitted skyWcs.
1605 stars : `SourceCatalog`
1606 Good stars selected for use in calibration, with RA/Dec coordinates
1607 computed from the pixel positions and fitted WCS.
1609 Returns
1610 -------
1611 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1612 Reference/stars matches used in the fit.
1613 """
1614 initialWcs = exposure.wcs
1615 result = self.astrometry.run(stars, exposure)
1617 # Record astrometry stats to metadata
1618 self.metadata["astrometry_matches_count"] = len(result.matches)
1619 if exposure.wcs is not None:
1620 if initialWcs is not None:
1621 center = exposure.getBBox().getCenter()
1622 self.metadata['initial_to_final_wcs'] = (
1623 initialWcs.pixelToSky(center).separation(
1624 exposure.wcs.pixelToSky(center)
1625 ).asArcseconds()
1626 )
1627 else:
1628 self.metadata['initial_to_final_wcs'] = float("nan")
1629 self.metadata['astrom_offset_mean'] = exposure.metadata['SFM_ASTROM_OFFSET_MEAN']
1630 self.metadata['astrom_offset_std'] = exposure.metadata['SFM_ASTROM_OFFSET_STD']
1631 self.metadata['astrom_offset_median'] = exposure.metadata['SFM_ASTROM_OFFSET_MEDIAN']
1632 else:
1633 self.metadata['initial_to_final_wcs'] = float("nan")
1634 self.metadata['astrom_offset_mean'] = float("nan")
1635 self.metadata['astrom_offset_std'] = float("nan")
1636 self.metadata['astrom_offset_median'] = float("nan")
1638 return result.matches, result.matchMeta
1640 def _fit_photometry(self, exposure, stars):
1641 """Fit a photometric model to the data and return the reference
1642 matches used in the fit, and the fitted PhotoCalib.
1644 Parameters
1645 ----------
1646 exposure : `lsst.afw.image.Exposure`
1647 Exposure that is being fit, to get PSF and other metadata from.
1648 Has the fit `lsst.afw.image.PhotoCalib` attached, with pixel values
1649 unchanged.
1650 stars : `lsst.afw.table.SourceCatalog`
1651 Good stars selected for use in calibration.
1652 background : `lsst.afw.math.BackgroundList`
1653 Background that was fit to the exposure during detection of the
1654 above stars.
1656 Returns
1657 -------
1658 calibrated_stars : `lsst.afw.table.SourceCatalog`
1659 Star catalog with flux/magnitude columns computed from the fitted
1660 photoCalib (instFlux columns are retained as well).
1661 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1662 Reference/stars matches used in the fit.
1663 matchMeta : `lsst.daf.base.PropertyList`
1664 Metadata needed to unpersist matches, as returned by the matcher.
1665 photo_calib : `lsst.afw.image.PhotoCalib`
1666 Photometric calibration that was fit to the star catalog.
1667 """
1668 result = self.photometry.run(exposure, stars)
1669 calibrated_stars = result.photoCalib.calibrateCatalog(stars)
1670 exposure.setPhotoCalib(result.photoCalib)
1671 return calibrated_stars, result.matches, result.matchMeta, result.photoCalib
1673 def _apply_photometry(self, exposure, background, background_to_photometric_ratio=None):
1674 """Apply the photometric model attached to the exposure to the
1675 exposure's pixels and an associated background model.
1677 Parameters
1678 ----------
1679 exposure : `lsst.afw.image.Exposure`
1680 Exposure with the target `lsst.afw.image.PhotoCalib` attached.
1681 On return, pixel values will be calibrated and an identity
1682 photometric transform will be attached.
1683 background : `lsst.afw.math.BackgroundList`
1684 Background model to convert to nanojansky units in place.
1685 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1686 Image to convert photometric-flattened image to
1687 background-flattened image.
1688 """
1689 photo_calib = exposure.getPhotoCalib()
1690 exposure.maskedImage = photo_calib.calibrateImage(exposure.maskedImage)
1691 identity = afwImage.PhotoCalib(1.0,
1692 photo_calib.getCalibrationErr(),
1693 bbox=exposure.getBBox())
1694 exposure.setPhotoCalib(identity)
1695 exposure.metadata["BUNIT"] = "nJy"
1697 assert photo_calib._isConstant, \
1698 "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that."
1700 for bg in background:
1701 # The statsImage is a view, but we can't assign to a function call
1702 # in python.
1703 binned_image = bg[0].getStatsImage()
1704 binned_image *= photo_calib.getCalibrationMean()
1706 def _summarize(self, exposure, stars, background):
1707 """Compute summary statistics on the exposure and update in-place the
1708 calibrations attached to it.
1710 Parameters
1711 ----------
1712 exposure : `lsst.afw.image.Exposure`
1713 Exposure that was calibrated, to get PSF and other metadata from.
1714 Should be in instrumental units with the photometric calibration
1715 attached.
1716 Modified to contain the computed summary statistics.
1717 stars : `SourceCatalog`
1718 Good stars selected used in calibration.
1719 background : `lsst.afw.math.BackgroundList`
1720 Background that was fit to the exposure during detection of the
1721 above stars. Should be in instrumental units.
1722 """
1723 summary = self.compute_summary_stats.run(exposure, stars, background)
1724 exposure.info.setSummaryStats(summary)
1726 def _recordMaskedPixelFractions(self, exposure):
1727 """Record the fraction of all the pixels in an exposure
1728 that are masked with a given flag. Each fraction is
1729 recorded in the task metadata. One record per flag type.
1731 Parameters
1732 ----------
1733 exposure : `lsst.afw.image.ExposureF`
1734 The target exposure to calculate masked pixel fractions for.
1735 """
1737 mask = exposure.mask
1738 maskPlanes = list(mask.getMaskPlaneDict().keys())
1739 for maskPlane in maskPlanes:
1740 self.metadata[f"{maskPlane.lower()}_mask_fraction"] = (
1741 evaluateMaskFraction(mask, maskPlane)
1742 )
1744 def _update_wcs_with_camera_model(self, exposure, cameraModel):
1745 """Replace the existing WCS with one generated using the pointing
1746 and the input camera model.
1748 Parameters
1749 ----------
1750 exposure : `lsst.afw.image.ExposureF`
1751 The exposure whose WCS will be updated.
1752 cameraModel : `lsst.afw.cameraGeom.Camera`
1753 Camera to be used when constructing updated WCS.
1754 """
1755 detector = cameraModel[exposure.detector.getId()]
1756 newWcs = createInitialSkyWcs(exposure.visitInfo, detector)
1757 exposure.setWcs(newWcs)
1759 def _compute_mask_fraction(self, mask, mask_planes, bad_mask_planes):
1760 """Evaluate the fraction of masked pixels in a (set of) mask plane(s).
1762 Parameters
1763 ----------
1764 mask : `lsst.afw.image.Mask`
1765 The mask on which to evaluate the fraction.
1766 mask_planes : `list`, `str`
1767 The mask planes for which to evaluate the fraction.
1768 bad_mask_planes : `list`, `str`
1769 The mask planes to exclude from the computation.
1771 Returns
1772 -------
1773 detected_fraction : `float`
1774 The calculated fraction of masked pixels
1775 """
1776 bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes)
1777 n_good_pix = np.sum(mask.array & bad_pixel_mask == 0)
1778 if n_good_pix == 0:
1779 detected_fraction = float("nan")
1780 return detected_fraction
1781 detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes)
1782 n_detected_pix = np.sum((mask.array & detected_pixel_mask != 0)
1783 & (mask.array & bad_pixel_mask == 0))
1784 detected_fraction = n_detected_pix/n_good_pix
1785 return detected_fraction
1787 def _compute_per_amp_fraction(self, exposure, detected_fraction, mask_planes, bad_mask_planes):
1788 """Evaluate the maximum per-amplifier fraction of masked pixels.
1790 Parameters
1791 ----------
1792 exposure : `lsst.afw.image.ExposureF`
1793 The exposure on which to compute the per-amp masked fraction.
1794 detected_fraction : `float`
1795 The current detected_fraction of the ``mask_planes`` for the
1796 full detector.
1797 mask_planes : `list`, `str`
1798 The mask planes for which to evaluate the fraction.
1799 bad_mask_planes : `list`, `str`
1800 The mask planes to exclude from the computation.
1802 Returns
1803 -------
1804 n_above_max_per_amp : `int`
1805 The number of amplifiers with masked fractions above a maximum
1806 value (set by the current full-detector ``detected_fraction``).
1807 highest_detected_fraction_per_amp : `float`
1808 The highest value of the per-amplifier fraction of masked pixels.
1809 no_zero_det_amps : `bool`
1810 A boolean representing whether any of the amplifiers has zero
1811 masked pixels.
1812 """
1813 highest_detected_fraction_per_amp = -9.99
1814 n_above_max_per_amp = 0
1815 n_no_zero_det_amps = 0
1816 no_zero_det_amps = True
1817 amps = exposure.detector.getAmplifiers()
1818 if amps is not None:
1819 for ia, amp in enumerate(amps):
1820 amp_bbox = amp.getBBox()
1821 exp_bbox = exposure.getBBox()
1822 if not exp_bbox.contains(amp_bbox):
1823 self.log.info("Bounding box of amplifier (%s) does not fit in exposure's "
1824 "bounding box (%s). Skipping...", amp_bbox, exp_bbox)
1825 continue
1826 sub_image = exposure.subset(amp.getBBox())
1827 detected_fraction_amp = self._compute_mask_fraction(sub_image.mask,
1828 mask_planes,
1829 bad_mask_planes)
1830 self.log.debug("Current detected fraction for amplifier %s = %.3f",
1831 amp.getName(), detected_fraction_amp)
1832 if detected_fraction_amp < 0.002:
1833 n_no_zero_det_amps += 1
1834 if n_no_zero_det_amps > 2:
1835 no_zero_det_amps = False
1836 break
1837 highest_detected_fraction_per_amp = max(detected_fraction_amp,
1838 highest_detected_fraction_per_amp)
1839 if highest_detected_fraction_per_amp > min(0.998, max(0.8, 3.0*detected_fraction)):
1840 n_above_max_per_amp += 1
1841 if n_above_max_per_amp > 2:
1842 break
1843 else:
1844 self.log.info("No amplifier object for detector %d, so skipping per-amp "
1845 "detection fraction checks.", exposure.detector.getId())
1846 return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps
1848 def _remeasure_star_background(self, result, background_to_photometric_ratio=None):
1849 """Remeasure the exposure's background with iterative adaptive
1850 threshold detection.
1852 Parameters
1853 ----------
1854 result : `lsst.pipe.base.Struct`
1855 The current state of the result Strut from the run method of
1856 CalibrateImageTask. Will be modified in place.
1857 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1858 Image to convert photometric-flattened image to
1859 background-flattened image.
1861 Returns
1862 -------
1863 result : `lsst.pipe.base.Struct`
1864 The modified result Struct with the new background subtracted.
1865 """
1866 # Restore the previously measured background and remeasure it
1867 # using an adaptive threshold detection iteration to ensure a
1868 # "Goldilocks Zone" for the fraction of detected pixels.
1869 backgroundOrig = result.background.clone()
1870 median_background_orig = np.median(backgroundOrig.getImage().array)
1871 self.log.info("Original median_background = %.2f", median_background_orig)
1872 result.exposure.image.array += result.background.getImage().array
1873 result.background = afwMath.BackgroundList()
1875 origMask = result.exposure.mask.clone()
1876 bad_mask_planes = ["BAD", "EDGE", "NO_DATA"]
1877 detected_mask_planes = ["DETECTED", "DETECTED_NEGATIVE"]
1878 detected_fraction_orig = self._compute_mask_fraction(result.exposure.mask,
1879 detected_mask_planes,
1880 bad_mask_planes)
1881 minDetFracForFinalBg = 0.02
1882 maxDetFracForFinalBg = 0.93
1883 # Dilate the current detected mask planes and don't clear
1884 # them in the detection step.
1885 nPixToDilate = 10
1886 maxAdaptiveDetIter = 10
1887 for dilateIter in range(maxAdaptiveDetIter):
1888 dilatedMask = origMask.clone()
1889 for maskName in detected_mask_planes:
1890 # Compute the grown detection mask plane using SpanSet
1891 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)
1892 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)
1893 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)
1894 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())
1895 # Clear the detected mask plane
1896 detectedMask = dilatedMask.getMaskPlane(maskName)
1897 dilatedMask.clearMaskPlane(detectedMask)
1898 # Set the mask plane to the dilated one
1899 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)
1901 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,
1902 detected_mask_planes,
1903 bad_mask_planes)
1904 if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1:
1905 break
1906 else:
1907 nPixToDilate -= 1
1909 result.exposure.mask = dilatedMask
1910 self.log.warning("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f",
1911 detected_fraction_orig, detected_fraction_dilated)
1912 n_above_max_per_amp = -99
1913 highest_detected_fraction_per_amp = float("nan")
1915 # Check the per-amplifier detection fractions.
1916 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \
1917 self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated,
1918 detected_mask_planes, bad_mask_planes)
1919 self.log.warning("Dilated mask: n_above_max_per_amp = %d, "
1920 "highest_detected_fraction_per_amp = %.3f",
1921 n_above_max_per_amp, highest_detected_fraction_per_amp)
1923 bgIgnoreMasksToAdd = ["SAT", "SUSPECT", "SPIKE"]
1924 detected_fraction = 1.0
1925 nFootprintTemp = 1e12
1926 starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
1927 for maskName in bgIgnoreMasksToAdd:
1928 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
1929 and maskName not in starBackgroundDetectionConfig.background.ignoredPixelMask):
1930 starBackgroundDetectionConfig.background.ignoredPixelMask += [maskName]
1931 starBackgroundDetectionConfig.doTempLocalBackground = False
1932 starBackgroundDetectionConfig.nSigmaToGrow = 70.0
1933 starBackgroundDetectionConfig.reEstimateBackground = False
1934 starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0
1935 starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background_orig)
1936 starBackgroundDetectionConfig.thresholdType = "pixel_stdev"
1938 n_above_max_per_amp = -99
1939 highest_detected_fraction_per_amp = float("nan")
1940 doCheckPerAmpDetFraction = True
1942 minFootprints = self.config.star_background_min_footprints
1943 maxIter = 40
1944 for nIter in range(maxIter):
1945 currentThresh = starBackgroundDetectionConfig.thresholdValue
1946 nZeroEncountered = 0
1947 if nFootprintTemp == 0:
1948 zeroFactor = min(0.98, 0.9 + 0.01*nZeroEncountered)
1949 starBackgroundDetectionConfig.thresholdValue = zeroFactor*currentThresh
1950 self.log.warning("No footprints detected. Decreasing threshold to %.2f and rerunning.",
1951 starBackgroundDetectionConfig.thresholdValue)
1952 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
1953 config=starBackgroundDetectionConfig)
1954 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)
1955 tempDetections = starBackgroundDetectionTask.run(
1956 table=table, exposure=result.exposure, clearMask=True)
1957 nFootprintTemp = len(tempDetections.sources)
1958 minFootprints = max(self.config.star_background_min_footprints,
1959 int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))
1960 minFootprints = min(200, minFootprints)
1961 nZeroEncountered += 1
1962 if nFootprintTemp >= minFootprints:
1963 detected_fraction = self._compute_mask_fraction(result.exposure.mask,
1964 detected_mask_planes,
1965 bad_mask_planes)
1966 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
1967 "DETECTED_NEGATIVE in star_background_detection = %.3f "
1968 "(max is %.3f; min is %.3f) nFootprint = %d (current min is %d)",
1969 nIter, starBackgroundDetectionConfig.thresholdValue,
1970 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
1971 nFootprintTemp, minFootprints)
1972 self.metadata['adaptive_threshold_value'] = starBackgroundDetectionConfig.thresholdValue
1974 break
1975 else:
1976 # Still not enough footprints, so make sure this loop is
1977 # entered again.
1978 if nFootprintTemp > 0 and nFootprintTemp < minFootprints:
1979 nFootprintTemp = 0
1980 continue
1981 if detected_fraction > maxDetFracForFinalBg or nFootprintTemp <= minFootprints:
1982 starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh
1983 if nFootprintTemp < minFootprints and detected_fraction > 0.9*maxDetFracForFinalBg:
1984 if nFootprintTemp == 1:
1985 starBackgroundDetectionConfig.thresholdValue = 1.4*currentThresh
1986 else:
1987 starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh
1989 if n_above_max_per_amp > 1:
1990 starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh
1991 if detected_fraction < minDetFracForFinalBg:
1992 starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh
1993 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
1994 config=starBackgroundDetectionConfig)
1995 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)
1996 tempDetections = starBackgroundDetectionTask.run(
1997 table=table, exposure=result.exposure, clearMask=True)
1998 result.exposure.mask |= dilatedMask
1999 nFootprintTemp = len(tempDetections.sources)
2000 minFootprints = max(self.config.star_background_min_footprints,
2001 int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))
2002 minFootprints = min(200, minFootprints)
2003 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
2004 bad_mask_planes)
2005 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
2006 "DETECTED_NEGATIVE in star_background_detection = %.3f "
2007 "(max is %.3f; min is %.3f) nFootprint = %d (current min is %d)",
2008 nIter, starBackgroundDetectionConfig.thresholdValue,
2009 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
2010 nFootprintTemp, minFootprints)
2011 self.metadata['adaptive_threshold_value'] = starBackgroundDetectionConfig.thresholdValue
2013 n_amp = len(result.exposure.detector.getAmplifiers())
2014 if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg:
2015 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \
2016 self._compute_per_amp_fraction(result.exposure, detected_fraction,
2017 detected_mask_planes, bad_mask_planes)
2019 if not no_zero_det_amps:
2020 starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh
2022 if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg
2023 and n_above_max_per_amp < int(0.75*n_amp)
2024 and no_zero_det_amps
2025 and nFootprintTemp >= minFootprints):
2026 if (n_above_max_per_amp < max(1, int(0.15*n_amp))
2027 or detected_fraction < 0.85*maxDetFracForFinalBg):
2028 break
2029 else:
2030 self.log.warning("Making small tweak....")
2031 starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh
2032 self.log.warning("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp))
2034 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
2035 bad_mask_planes)
2036 self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f "
2037 "(highest per amp section = %.5f)",
2038 detected_fraction, highest_detected_fraction_per_amp)
2040 if detected_fraction > maxDetFracForFinalBg:
2041 result.exposure.mask = dilatedMask
2042 self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE "
2043 "was too large in star_background_detection = %.3f (max = %.3f). "
2044 "Reverting to dilated mask from PSF detection...",
2045 detected_fraction, maxDetFracForFinalBg)
2046 self.metadata['adaptive_threshold_value'] = float("nan")
2047 star_background = self.star_background.run(
2048 exposure=result.exposure,
2049 backgroundToPhotometricRatio=background_to_photometric_ratio,
2050 ).background
2051 result.background.append(star_background[0])
2053 median_background = np.median(result.background.getImage().array)
2054 self.log.info("New initial median_background = %.2f", median_background)
2056 # Perform one more round of background subtraction that is just an
2057 # overall pedestal (order = 0). This is intended to account for
2058 # any potential gross oversubtraction imposed by the higher-order
2059 # subtraction.
2060 # Dilate DETECTED mask a bit more if it's below 50% detected.
2061 nPixToDilate = 2
2062 if detected_fraction < 0.5:
2063 dilatedMask = result.exposure.mask.clone()
2064 for maskName in detected_mask_planes:
2065 # Compute the grown detection mask plane using SpanSet
2066 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)
2067 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)
2068 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)
2069 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())
2070 # Clear the detected mask plane
2071 detectedMask = dilatedMask.getMaskPlane(maskName)
2072 dilatedMask.clearMaskPlane(detectedMask)
2073 # Set the mask plane to the dilated one
2074 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)
2076 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,
2077 detected_mask_planes,
2078 bad_mask_planes)
2079 result.exposure.mask = dilatedMask
2080 self.log.debug("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f",
2081 detected_fraction_orig, detected_fraction_dilated)
2083 bbox = result.exposure.getBBox()
2085 # Now measure a final background pedestal level (largely to account
2086 # for the over-subtraction often seen in the first pass, especially
2087 # in high fill-factor fields). Since the pedestal measurement can be
2088 # sensitive to the bin size, start an iteration with a small bin size,
2089 # then double it on each iteration until the relative or absolute
2090 # change criterion is met. If those are never achieved, the iteration
2091 # stops when the bin size gets bigger than the exposure's bounding box.
2093 # Initialize the parameters for the pedestal iteration.
2094 pedestalBackgroundConfig = lsst.meas.algorithms.SubtractBackgroundConfig()
2095 for maskName in bgIgnoreMasksToAdd:
2096 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
2097 and maskName not in pedestalBackgroundConfig.ignoredPixelMask):
2098 pedestalBackgroundConfig.ignoredPixelMask += [maskName]
2099 pedestalBackgroundConfig.statisticsProperty = "MEDIAN"
2100 pedestalBackgroundConfig.approxOrderX = 0
2101 pedestalBackgroundConfig.binSize = min(32, int(math.ceil(max(bbox.width, bbox.height)/4.0)))
2102 self.log.info("Initial pedestal binSize = %d pixels", pedestalBackgroundConfig.binSize)
2103 cumulativePedestalLevel = 0.0
2104 # When the cumulative pedestal value changes by less than 5% from one
2105 # bin size to the next we assume we are converged.
2106 relativeStoppingCriterion = 0.05
2107 # When the cumulative pedestal value changes by less than 0.5 counts
2108 # (electrons or adu) from one bin size to the next, we assume we are
2109 # converged.
2110 absoluteStoppingCriterion = 0.5
2111 inPedestalIteration = True
2112 while inPedestalIteration:
2113 cumulativePedestalPrev = cumulativePedestalLevel
2114 pedestalBackgroundTask = lsst.meas.algorithms.SubtractBackgroundTask(
2115 config=pedestalBackgroundConfig)
2116 pedestalBackgroundList = pedestalBackgroundTask.run(
2117 exposure=result.exposure,
2118 background=result.background,
2119 backgroundToPhotometricRatio=background_to_photometric_ratio,
2120 ).background
2121 # Isolate the most recent pedestal background measured to log the
2122 # current computed value and add it to the cumulative value.
2123 pedestalBackground = afwMath.BackgroundList()
2124 pedestalBackground.append(pedestalBackgroundList[-1])
2125 pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0]
2126 cumulativePedestalLevel += pedestalBackgroundLevel
2127 absoluteDelta = abs(cumulativePedestalLevel - cumulativePedestalPrev)
2128 if cumulativePedestalPrev != 0.0:
2129 relativeDelta = abs(absoluteDelta/cumulativePedestalPrev)
2130 else:
2131 relativeDelta = 1.0
2132 self.log.info("Subtracted pedestalBackgroundLevel = %.4f (cumulative = %.4f; "
2133 "relativeDelta = %.4f; absoluteDelta = %.3f)",
2134 pedestalBackgroundLevel, cumulativePedestalLevel,
2135 relativeDelta, absoluteDelta)
2137 if (relativeDelta < relativeStoppingCriterion
2138 or absoluteDelta < absoluteStoppingCriterion):
2139 # We are converged.
2140 inPedestalIteration = False
2141 else:
2142 # We have not yet converged; grow the bin size if possible.
2143 if pedestalBackgroundConfig.binSize*2 < 2*max(bbox.width, bbox.height):
2144 pedestalBackgroundConfig.binSize = int(pedestalBackgroundConfig.binSize*2)
2145 self.log.info("Increasing pedestal binSize to %d pixels and remeasuring.",
2146 pedestalBackgroundConfig.binSize)
2147 else:
2148 # We have reached the maximum bin size.
2149 self.log.warning("Reached maximum bin size. Exiting pedestal loop without meeting "
2150 "the convergence criteria (relativeDelta <= %.2f or "
2151 "absoluteDelta <= %.2f).",
2152 relativeStoppingCriterion, absoluteStoppingCriterion)
2153 inPedestalIteration = False
2154 self.log.info("Final subtracted cumulativePedestalBackgroundLevel = %.4f", cumulativePedestalLevel)
2156 # Clear detected mask plane before final round of detection
2157 mask = result.exposure.mask
2158 for mp in detected_mask_planes:
2159 if mp not in mask.getMaskPlaneDict():
2160 mask.addMaskPlane(mp)
2161 mask &= ~mask.getPlaneBitMask(detected_mask_planes)
2163 return result