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