Coverage for python / lsst / ip / isr / isrTaskLSST.py: 12%
910 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-28 08:55 +0000
1__all__ = ["IsrTaskLSST", "IsrTaskLSSTConfig"]
3import numpy
4import math
6from . import isrFunctions
7from . import isrQa
8from .defects import Defects
10from contextlib import contextmanager
11from deprecated.sphinx import deprecated
12import lsst.pex.config as pexConfig
13import lsst.afw.math as afwMath
14import lsst.pipe.base as pipeBase
15import lsst.afw.image as afwImage
16import lsst.pipe.base.connectionTypes as cT
17from lsst.pipe.base import UnprocessableDataError
18from lsst.meas.algorithms.detection import SourceDetectionTask
19import lsst.afw.detection as afwDetection
21from .ampOffset import AmpOffsetTask
22from .binImageDataTask import BinImageDataTask
23from .overscan import SerialOverscanCorrectionTask, ParallelOverscanCorrectionTask
24from .overscanAmpConfig import OverscanCameraConfig
25from .assembleCcdTask import AssembleCcdTask
26from .deferredCharge import DeferredChargeTask
27from .crosstalk import CrosstalkTask
28from .masking import MaskingTask
29from .isrStatistics import IsrStatisticsTask
30from .isr import maskNans
31from .ptcDataset import PhotonTransferCurveDataset
32from .isrFunctions import isTrimmedExposure, compareCameraKeywords
33from .brighterFatterKernel import (brighterFatterCorrection,
34 fluxConservingBrighterFatterCorrection)
35from .electrostaticBrighterFatter import electrostaticBrighterFatterCorrection
38class IsrTaskLSSTConnections(pipeBase.PipelineTaskConnections,
39 dimensions={"instrument", "exposure", "detector"},
40 defaultTemplates={}):
41 ccdExposure = cT.Input(
42 name="raw",
43 doc="Input exposure to process.",
44 storageClass="Exposure",
45 dimensions=["instrument", "exposure", "detector"],
46 )
47 camera = cT.PrerequisiteInput(
48 name="camera",
49 storageClass="Camera",
50 doc="Input camera to construct complete exposures.",
51 dimensions=["instrument"],
52 isCalibration=True,
53 )
54 dnlLUT = cT.PrerequisiteInput(
55 name="dnlLUT",
56 doc="Look-up table for differential non-linearity.",
57 storageClass="IsrCalib",
58 dimensions=["instrument", "exposure", "detector"],
59 isCalibration=True,
60 # TODO DM 36636
61 )
62 bias = cT.PrerequisiteInput(
63 name="bias",
64 doc="Input bias calibration.",
65 storageClass="ExposureF",
66 dimensions=["instrument", "detector"],
67 isCalibration=True,
68 )
69 deferredChargeCalib = cT.PrerequisiteInput(
70 name="cti",
71 doc="Deferred charge/CTI correction dataset.",
72 storageClass="IsrCalib",
73 dimensions=["instrument", "detector"],
74 isCalibration=True,
75 )
76 linearizer = cT.PrerequisiteInput(
77 name='linearizer',
78 storageClass="Linearizer",
79 doc="Linearity correction calibration.",
80 dimensions=["instrument", "detector"],
81 isCalibration=True,
82 )
83 ptc = cT.PrerequisiteInput(
84 name="ptc",
85 doc="Input Photon Transfer Curve dataset",
86 storageClass="PhotonTransferCurveDataset",
87 dimensions=["instrument", "detector"],
88 isCalibration=True,
89 )
90 gainCorrection = cT.PrerequisiteInput(
91 name="gain_correction",
92 doc="Gain correction dataset",
93 storageClass="IsrCalib",
94 dimensions=["instrument", "detector"],
95 isCalibration=True,
96 minimum=0,
97 )
98 crosstalk = cT.PrerequisiteInput(
99 name="crosstalk",
100 doc="Input crosstalk object",
101 storageClass="CrosstalkCalib",
102 dimensions=["instrument", "detector"],
103 isCalibration=True,
104 )
105 defects = cT.PrerequisiteInput(
106 name='defects',
107 doc="Input defect tables.",
108 storageClass="Defects",
109 dimensions=["instrument", "detector"],
110 isCalibration=True,
111 )
112 bfKernel = cT.PrerequisiteInput(
113 name="bfk",
114 doc="Complete kernel + gain solutions.",
115 storageClass="BrighterFatterKernel",
116 dimensions=["instrument", "detector"],
117 isCalibration=True,
118 )
119 electroBfDistortionMatrix = cT.PrerequisiteInput(
120 name="electroBfDistortionMatrix",
121 doc="Electrostatic BF solution",
122 storageClass="IsrCalib",
123 dimensions=["instrument", "detector"],
124 isCalibration=True,
125 )
126 dark = cT.PrerequisiteInput(
127 name='dark',
128 doc="Input dark calibration.",
129 storageClass="ExposureF",
130 dimensions=["instrument", "detector"],
131 isCalibration=True,
132 )
133 flat = cT.PrerequisiteInput(
134 name="flat",
135 doc="Input flat calibration.",
136 storageClass="ExposureF",
137 dimensions=["instrument", "detector", "physical_filter"],
138 isCalibration=True,
139 )
140 outputExposure = cT.Output(
141 name='postISRCCD',
142 doc="Output ISR processed exposure.",
143 storageClass="Exposure",
144 dimensions=["instrument", "exposure", "detector"],
145 )
146 preInterpExposure = cT.Output(
147 name='preInterpISRCCD',
148 doc="Output ISR processed exposure, with pixels left uninterpolated.",
149 storageClass="ExposureF",
150 dimensions=["instrument", "exposure", "detector"],
151 )
152 outputBin1Exposure = cT.Output(
153 name="postIsrBin1",
154 doc="First binned image.",
155 storageClass="ExposureF",
156 dimensions=["instrument", "exposure", "detector"],
157 )
158 outputBin2Exposure = cT.Output(
159 name="postIsrBin2",
160 doc="Second binned image.",
161 storageClass="ExposureF",
162 dimensions=["instrument", "exposure", "detector"],
163 )
165 outputStatistics = cT.Output(
166 name="isrStatistics",
167 doc="Output of additional statistics table.",
168 storageClass="StructuredDataDict",
169 dimensions=["instrument", "exposure", "detector"],
170 )
172 def __init__(self, *, config=None):
173 super().__init__(config=config)
175 doApplyGains = config.doApplyGains
176 useLinearizerGains = config.useGainsFrom == "LINEARIZER"
178 if config.doBootstrap or (doApplyGains and useLinearizerGains):
179 del self.ptc
180 if not config.doCorrectGains:
181 del self.gainCorrection
182 if config.doDiffNonLinearCorrection is not True:
183 del self.dnlLUT
184 if config.doBias is not True:
185 del self.bias
186 if config.doDeferredCharge is not True:
187 del self.deferredChargeCalib
188 if (config.doLinearize or (doApplyGains and useLinearizerGains)) is not True:
189 del self.linearizer
191 if not config.doCrosstalk:
192 del self.crosstalk
193 if config.doDefect is not True:
194 del self.defects
195 if config.doBrighterFatter is not True:
196 del self.bfKernel
197 del self.electroBfDistortionMatrix
198 else:
199 if config.brighterFatterCorrectionMethod.startswith("COULTON"):
200 del self.electroBfDistortionMatrix
201 else:
202 del self.bfKernel
203 if config.doDark is not True:
204 del self.dark
205 if config.doFlat is not True:
206 del self.flat
208 if config.doBinnedExposures is not True:
209 del self.outputBin1Exposure
210 del self.outputBin2Exposure
211 if config.doSaveInterpPixels is not True:
212 del self.preInterpExposure
214 if config.doCalculateStatistics is not True:
215 del self.outputStatistics
218class IsrTaskLSSTConfig(pipeBase.PipelineTaskConfig,
219 pipelineConnections=IsrTaskLSSTConnections):
220 """Configuration parameters for IsrTaskLSST.
222 Items are grouped in the order in which they are executed by the task.
223 """
224 expectWcs = pexConfig.Field(
225 dtype=bool,
226 default=True,
227 doc="Expect input science images to have a WCS (set False for e.g. spectrographs)."
228 )
229 qa = pexConfig.ConfigField(
230 dtype=isrQa.IsrQaConfig,
231 doc="QA related configuration options.",
232 )
233 doHeaderProvenance = pexConfig.Field(
234 dtype=bool,
235 default=True,
236 doc="Write calibration identifiers into output exposure header.",
237 )
239 # Calib checking configuration:
240 doRaiseOnCalibMismatch = pexConfig.Field(
241 dtype=bool,
242 default=False,
243 doc="Should IsrTaskLSST halt if exposure and calibration header values do not match?",
244 )
245 cameraKeywordsToCompare = pexConfig.ListField(
246 dtype=str,
247 doc="List of header keywords to compare between exposure and calibrations.",
248 default=[],
249 )
251 # Differential non-linearity correction.
252 doDiffNonLinearCorrection = pexConfig.Field(
253 dtype=bool,
254 doc="Do differential non-linearity correction?",
255 default=False,
256 )
258 doBootstrap = pexConfig.Field(
259 dtype=bool,
260 default=False,
261 doc="Is this task to be run in a ``bootstrap`` fashion that does not require "
262 "a PTC or full calibrations?",
263 )
265 doCheckUnprocessableData = pexConfig.Field(
266 dtype=bool,
267 default=True,
268 doc="Check if this image is completely unprocessable due to all bad amps.",
269 )
271 overscanCamera = pexConfig.ConfigField(
272 dtype=OverscanCameraConfig,
273 doc="Per-detector and per-amplifier overscan configurations.",
274 )
276 serialOverscanMedianShiftSigmaThreshold = pexConfig.Field(
277 dtype=float,
278 default=numpy.inf,
279 doc="Number of sigma difference from per-amp overscan median (as compared to PTC) to "
280 "check if an amp is in a different state than the baseline PTC calib and should "
281 "be marked BAD. Set to np.inf/np.nan to turn off overscan median checking.",
282 )
284 ampNoiseThreshold = pexConfig.Field(
285 dtype=float,
286 default=25.0,
287 doc="Maximum amplifier noise (e-) that is allowed before an amp is masked as bad. "
288 "Set to np.inf/np.nan to turn off noise checking.",
289 )
291 bssVoltageMinimum = pexConfig.Field(
292 dtype=float,
293 default=5.0,
294 doc="Minimum back-side bias voltage. Below this the detector is ``off`` and an "
295 "UnprocessableDataError will be logged. Check will be skipped if doCheckUnprocessableData "
296 "is False or if value is less than or equal to 0.",
297 )
298 bssVoltageKeyword = pexConfig.Field(
299 dtype=str,
300 default="BSSVBS",
301 doc="Back-side bias voltage header keyword. Only checked if doCheckUnprocessableData is True "
302 "and bssVoltageMinimum is greater than 0.",
303 )
304 hvBiasKeyword = pexConfig.Field(
305 dtype=str,
306 default="HVBIAS",
307 doc="Back-side bias voltage on/off header keyword. Only checked if doCheckUnprocessableData is True "
308 "and bssVoltageMinimum is greater than 0.",
309 )
311 # Amplifier to CCD assembly configuration.
312 doAssembleCcd = pexConfig.Field(
313 dtype=bool,
314 default=True,
315 doc="Assemble amp-level exposures into a ccd-level exposure?"
316 )
317 assembleCcd = pexConfig.ConfigurableField(
318 target=AssembleCcdTask,
319 doc="CCD assembly task.",
320 )
322 # Bias subtraction.
323 doBias = pexConfig.Field(
324 dtype=bool,
325 doc="Apply bias frame correction?",
326 default=True,
327 )
329 # Deferred charge correction.
330 doDeferredCharge = pexConfig.Field(
331 dtype=bool,
332 doc="Apply deferred charge correction?",
333 default=True,
334 )
335 deferredChargeCorrection = pexConfig.ConfigurableField(
336 target=DeferredChargeTask,
337 doc="Deferred charge correction task.",
338 )
340 # Linearization.
341 doLinearize = pexConfig.Field(
342 dtype=bool,
343 doc="Correct for nonlinearity of the detector's response?",
344 default=True,
345 )
347 # Gains.
348 doCorrectGains = pexConfig.Field(
349 dtype=bool,
350 doc="Apply gain corrections from detector restarts?",
351 default=True,
352 )
353 doApplyGains = pexConfig.Field(
354 dtype=bool,
355 doc="Apply gains to the image?",
356 default=True,
357 )
358 useGainsFrom = pexConfig.ChoiceField(
359 dtype=str,
360 doc="Where to retrieve the gains. Unused if doBootstrap is True.",
361 allowed={
362 "PTC": "Use the gains from the inputPtc calibration.",
363 "LINEARIZER": "Use the gains from the linearizer calibration.",
364 },
365 default="PTC",
366 )
368 # Variance construction.
369 doVariance = pexConfig.Field(
370 dtype=bool,
371 doc="Calculate variance?",
372 default=True
373 )
374 maskNegativeVariance = pexConfig.Field(
375 dtype=bool,
376 doc="Mask pixels that claim a negative variance. This likely indicates a failure "
377 "in the measurement of the overscan at an edge due to the data falling off faster "
378 "than the overscan model can account for it.",
379 default=True,
380 )
381 negativeVarianceMaskName = pexConfig.Field(
382 dtype=str,
383 doc="Mask plane to use to mark pixels with negative variance, if `maskNegativeVariance` is True.",
384 default="BAD",
385 )
386 doSaturation = pexConfig.Field(
387 dtype=bool,
388 doc="Mask saturated pixels? NB: this is totally independent of the"
389 " interpolation option - this is ONLY setting the bits in the mask."
390 " To have them interpolated make sure doInterpolate=True and"
391 " maskListToInterpolate includes SAT.",
392 default=True,
393 )
394 saturatedMaskName = pexConfig.Field(
395 dtype=str,
396 doc="Name of mask plane to use in saturation detection and interpolation.",
397 default="SAT",
398 )
399 defaultSaturationSource = pexConfig.ChoiceField(
400 dtype=str,
401 doc="Source to retrieve default amp-level saturation values.",
402 allowed={
403 "NONE": "No default saturation values; only config overrides will be used.",
404 "CAMERAMODEL": "Use the default from the camera model (old defaults).",
405 "PTCTURNOFF": "Use the ptcTurnoff value as the saturation level.",
406 },
407 default="PTCTURNOFF",
408 )
409 doSuspect = pexConfig.Field(
410 dtype=bool,
411 doc="Mask suspect pixels?",
412 default=True,
413 )
414 suspectMaskName = pexConfig.Field(
415 dtype=str,
416 doc="Name of mask plane to use for suspect pixels.",
417 default="SUSPECT",
418 )
419 defaultSuspectSource = pexConfig.ChoiceField(
420 dtype=str,
421 doc="Source to retrieve default amp-level suspect values.",
422 allowed={
423 "NONE": "No default suspect values; only config overrides will be used.",
424 "CAMERAMODEL": "Use the default from the camera model (old defaults).",
425 "PTCTURNOFF": "Use the ptcTurnoff value as the suspect level.",
426 },
427 default="PTCTURNOFF",
428 )
430 # Crosstalk.
431 doCrosstalk = pexConfig.Field(
432 dtype=bool,
433 doc="Apply intra-CCD crosstalk correction?",
434 default=True,
435 )
436 crosstalk = pexConfig.ConfigurableField(
437 target=CrosstalkTask,
438 doc="Intra-CCD crosstalk correction.",
439 )
441 # Masking options.
442 doITLDipMask = pexConfig.Field(
443 dtype=bool,
444 doc="Apply ``ITL dip`` masking. The ``itlDipMaskPlane`` mask plane "
445 "will be added even if this configuration is False.",
446 default=True,
447 )
448 itlDipMaskPlanes = pexConfig.ListField(
449 dtype=str,
450 doc="Mask plane to use for ITL dip pixels.",
451 default=["SUSPECT", "ITL_DIP"],
452 )
453 doDefect = pexConfig.Field(
454 dtype=bool,
455 doc="Apply correction for CCD defects, e.g. hot pixels?",
456 default=True,
457 )
458 badAmps = pexConfig.ListField(
459 dtype=str,
460 doc="List of bad amps that should be masked as BAD in the defect code. "
461 "Value should be of form {detector_name}_{amp_name}, e.g. ``R42_S21_C07``. "
462 "Only used if doDefect is True.",
463 default=[],
464 )
465 doNanMasking = pexConfig.Field(
466 dtype=bool,
467 doc="Mask non-finite (NAN, inf) pixels. The UNMASKEDNAN mask plane "
468 "will be added even if this configuration is False.",
469 default=True,
470 )
471 doWidenSaturationTrails = pexConfig.Field(
472 dtype=bool,
473 doc="Widen bleed trails based on their width.",
474 default=False,
475 )
476 masking = pexConfig.ConfigurableField(
477 target=MaskingTask,
478 doc="Masking task."
479 )
480 doE2VEdgeBleedMask = pexConfig.Field(
481 dtype=bool,
482 doc="Mask flag-like edge bleeds from saturated columns "
483 "in E2V amplifiers.",
484 default=True,
485 )
486 e2vEdgeBleedSatMinArea = pexConfig.Field(
487 dtype=int,
488 doc="Minimum limit of saturated cores footprint area to apply edge"
489 "bleed masking in E2V amplifiers.",
490 default=10000,
491 )
492 e2vEdgeBleedSatMaxArea = pexConfig.Field(
493 dtype=int,
494 doc="Maximum limit of saturated cores footprint area to apply edge"
495 "bleed masking in E2V amplifiers.",
496 default=100000,
497 )
498 e2vEdgeBleedYMax = pexConfig.Field(
499 dtype=int,
500 doc="Height in pixels of edge bleed masking in E2V amplifiers (width"
501 "is the width of the amplifier).",
502 default=350,
503 )
504 doITLEdgeBleedMask = pexConfig.Field(
505 dtype=bool,
506 doc="Mask edge bleeds from saturated columns in ITL amplifiers.",
507 default=True,
508 )
509 doITLSatSagMask = pexConfig.Field(
510 dtype=bool,
511 doc="Mask columns presenting saturation sag.",
512 default=True,
513 )
514 itlEdgeBleedSatMinArea = pexConfig.Field(
515 dtype=int,
516 doc="Minimum limit for saturated cores footprint area.",
517 default=10000,
518 )
519 itlEdgeBleedSatMaxArea = pexConfig.Field(
520 dtype=int,
521 doc="Maximum limit for saturated cores footprint area.",
522 default=100000,
523 )
524 itlEdgeBleedThreshold = pexConfig.Field(
525 dtype=float,
526 doc="Sky background threshold for edge bleed detection.",
527 default=5000.,
528 )
529 itlEdgeBleedModelConstant = pexConfig.Field(
530 dtype=float,
531 doc="Constant in the edge bleed exponential decay model.",
532 default=0.02,
533 )
535 # Interpolation options.
536 doInterpolate = pexConfig.Field(
537 dtype=bool,
538 doc="Interpolate masked pixels?",
539 default=True,
540 )
541 maskListToInterpolate = pexConfig.ListField(
542 dtype=str,
543 doc="List of mask planes that should be interpolated.",
544 default=["SAT", "BAD", "UNMASKEDNAN"],
545 )
546 doSaveInterpPixels = pexConfig.Field(
547 dtype=bool,
548 doc="Save a copy of the pre-interpolated pixel values?",
549 default=False,
550 )
551 useLegacyInterp = pexConfig.Field(
552 dtype=bool,
553 doc="Use the legacy interpolation algorithm. If False use Gaussian Process.",
554 default=True,
555 )
557 # Amp offset correction.
558 doAmpOffset = pexConfig.Field(
559 doc="Calculate amp offset corrections?",
560 dtype=bool,
561 default=False,
562 )
563 ampOffset = pexConfig.ConfigurableField(
564 doc="Amp offset correction task.",
565 target=AmpOffsetTask,
566 )
568 # Initial masking options.
569 doSetBadRegions = pexConfig.Field(
570 dtype=bool,
571 doc="Should we set the level of all BAD patches of the chip to the chip's average value?",
572 default=True,
573 )
575 # Brighter-Fatter correction.
576 doBrighterFatter = pexConfig.Field(
577 dtype=bool,
578 doc="Apply the brighter-fatter correction?",
579 default=True,
580 )
581 brighterFatterCorrectionMethod = pexConfig.ChoiceField(
582 dtype=str,
583 doc="The method for brighter-fatter correction.",
584 default="COULTON18",
585 allowed={
586 "COULTON18": "Coulton et al. 2018 BF correction with kernel",
587 "COULTON18_FLUX_CONSERVING": "Coulton et al. 2018 BF correction "
588 "with kernel + Flux conserving corrections",
589 "ASTIER23": "Astier & Regenault 2023 electrostatic BF correction",
590 "ASTIER23+FILTERCORRECTION": "Astier & Regenault 2023 electrostatic BF "
591 "correction + color correction",
592 },
593 )
594 brighterFatterLevel = pexConfig.ChoiceField(
595 dtype=str,
596 doc="The level at which to correct for brighter-fatter.",
597 allowed={
598 "AMP": "Every amplifier treated separately.",
599 "DETECTOR": "One kernel per detector.",
600 },
601 default="DETECTOR",
602 )
603 brighterFatterMaxIter = pexConfig.Field(
604 dtype=int,
605 doc="Maximum number of iterations for the brighter-fatter correction for "
606 "the COULTON18* correction types",
607 default=10,
608 )
609 brighterFatterThreshold = pexConfig.Field(
610 dtype=float,
611 doc="Threshold used to stop iterating the brighter-fatter correction for "
612 "the COULTON18* correction types. It is the absolute value of the difference "
613 "between the current corrected image and the one from the previous iteration "
614 "summed over all the pixels.",
615 default=1000,
616 )
617 brighterFatterMaskListToInterpolate = pexConfig.ListField(
618 dtype=str,
619 doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
620 "correction.",
621 default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
622 )
623 brighterFatterMaskGrowSize = pexConfig.Field(
624 dtype=int,
625 doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
626 "when brighter-fatter correction is applied.",
627 default=2,
628 )
629 brighterFatterFwhmForInterpolation = pexConfig.Field(
630 dtype=float,
631 doc="FWHM of PSF in arcseconds used for interpolation in brighter-fatter correction "
632 "(currently unused).",
633 default=1.0,
634 )
635 growSaturationFootprintSize = pexConfig.Field(
636 dtype=int,
637 doc="Number of pixels by which to grow the saturation footprints.",
638 default=1,
639 )
641 # Dark subtraction.
642 doDark = pexConfig.Field(
643 dtype=bool,
644 doc="Apply dark frame correction.",
645 default=True,
646 )
648 # Flat correction.
649 doFlat = pexConfig.Field(
650 dtype=bool,
651 doc="Apply flat field correction.",
652 default=True,
653 )
654 flatScalingType = pexConfig.ChoiceField(
655 dtype=str,
656 doc="The method for scaling the flat on the fly.",
657 default='USER',
658 allowed={
659 "USER": "Scale by flatUserScale",
660 "MEAN": "Scale by the inverse of the mean",
661 "MEDIAN": "Scale by the inverse of the median",
662 },
663 )
664 flatUserScale = pexConfig.Field(
665 dtype=float,
666 doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise.",
667 default=1.0,
668 )
669 noDataVariance = pexConfig.Field(
670 dtype=float,
671 doc="Value to fill for variance plane in NO_DATA regions (which are generally "
672 "inherited from the flat in fully vignetted regions.)",
673 default=1e10,
674 )
676 # Calculate image quality statistics?
677 doStandardStatistics = pexConfig.Field(
678 dtype=bool,
679 doc="Should standard image quality statistics be calculated?",
680 default=True,
681 )
682 # Calculate additional statistics?
683 doCalculateStatistics = pexConfig.Field(
684 dtype=bool,
685 doc="Should additional ISR statistics be calculated?",
686 default=True,
687 )
688 isrStats = pexConfig.ConfigurableField(
689 target=IsrStatisticsTask,
690 doc="Task to calculate additional statistics.",
691 )
693 # Make binned images?
694 doBinnedExposures = pexConfig.Field(
695 dtype=bool,
696 doc="Should binned exposures be calculated?",
697 default=False,
698 )
699 binning = pexConfig.ConfigurableField(
700 target=BinImageDataTask,
701 doc="Task to bin the exposure.",
702 )
703 binFactor1 = pexConfig.Field(
704 dtype=int,
705 doc="Binning factor for first binned exposure. This is intended for a finely binned output.",
706 default=8,
707 check=lambda x: x > 1,
708 )
709 binFactor2 = pexConfig.Field(
710 dtype=int,
711 doc="Binning factor for second binned exposure. This is intended for a coarsely binned output.",
712 default=64,
713 check=lambda x: x > 1,
714 )
716 def validate(self):
717 super().validate()
719 if self.doBootstrap:
720 # Additional checks in bootstrap (no gains) mode.
721 if self.doApplyGains:
722 raise ValueError("Cannot run task with doBootstrap=True and doApplyGains=True.")
723 if self.doCorrectGains:
724 raise ValueError("Cannot run task with doBootstrap=True and doCorrectGains=True.")
725 if self.doCrosstalk and self.crosstalk.doQuadraticCrosstalkCorrection:
726 raise ValueError("Cannot apply quadratic crosstalk correction with doBootstrap=True.")
727 if numpy.isfinite(self.serialOverscanMedianShiftSigmaThreshold):
728 raise ValueError("Cannot do amp overscan level checks with doBootstrap=True.")
729 if numpy.isfinite(self.ampNoiseThreshold):
730 raise ValueError("Cannot do amp noise thresholds with doBootstrap=True.")
732 if self.doITLEdgeBleedMask and not self.doSaturation:
733 raise ValueError("Cannot do ITL edge bleed masking when doSaturation=False.")
734 if self.doE2VEdgeBleedMask and not self.doSaturation:
735 raise ValueError("Cannot do e2v edge bleed masking when doSaturation=False.")
737 def setDefaults(self):
738 super().setDefaults()
741class IsrTaskLSST(pipeBase.PipelineTask):
742 ConfigClass = IsrTaskLSSTConfig
743 _DefaultName = "isrLSST"
745 def __init__(self, **kwargs):
746 super().__init__(**kwargs)
747 self.makeSubtask("assembleCcd")
748 self.makeSubtask("deferredChargeCorrection")
749 self.makeSubtask("crosstalk")
750 self.makeSubtask("masking")
751 self.makeSubtask("isrStats")
752 self.makeSubtask("ampOffset")
753 self.makeSubtask("binning")
755 def runQuantum(self, butlerQC, inputRefs, outputRefs):
757 inputs = butlerQC.get(inputRefs)
758 self.validateInput(inputs)
760 if self.config.doHeaderProvenance:
761 # Add calibration provenanace info to header.
762 exposureMetadata = inputs['ccdExposure'].metadata
763 for inputName in sorted(list(inputs.keys())):
764 reference = getattr(inputRefs, inputName, None)
765 if reference is not None and hasattr(reference, "run"):
766 runKey = f"LSST CALIB RUN {inputName.upper()}"
767 runValue = reference.run
768 idKey = f"LSST CALIB UUID {inputName.upper()}"
769 idValue = str(reference.id)
770 dateKey = f"LSST CALIB DATE {inputName.upper()}"
771 dateValue = self.extractCalibDate(inputs[inputName])
772 if dateValue != "Unknown Unknown":
773 butlerQC.add_additional_provenance(reference, {"calib date": dateValue})
775 exposureMetadata[runKey] = runValue
776 exposureMetadata[idKey] = idValue
777 exposureMetadata[dateKey] = dateValue
779 outputs = self.run(**inputs)
780 butlerQC.put(outputs, outputRefs)
782 def validateInput(self, inputs):
783 """
784 This is a check that all the inputs required by the config
785 are available.
786 """
788 # Inputs depend on where the gains are coming from
789 doApplyGains = self.config.doApplyGains
790 useLinearizerGains = self.config.useGainsFrom == "LINEARIZER"
791 usePtcGains = not useLinearizerGains
793 # Inputs depend on the BF method
794 method = self.config.brighterFatterCorrectionMethod
796 inputMap = {
797 'dnlLUT': self.config.doDiffNonLinearCorrection,
798 'bias': self.config.doBias,
799 'deferredChargeCalib': self.config.doDeferredCharge,
800 # Some tasks require gains in order to be
801 # supplied regardless of whether
802 # self.config.doApplyGains is True or False.
803 'linearizer': (self.config.doLinearize or (doApplyGains and useLinearizerGains)),
804 'ptc': self.config.doApplyGains and usePtcGains,
805 'crosstalk': self.config.doCrosstalk,
806 'defects': self.config.doDefect,
807 'bfKernel': self.config.doBrighterFatter & method.startswith("COULTON18"),
808 'electroBfDistortionMatrix': self.config.doBrighterFatter & (method == "ASTIER23"),
809 'dark': self.config.doDark,
810 }
812 for calibrationFile, configValue in inputMap.items():
813 if configValue and inputs[calibrationFile] is None:
814 raise RuntimeError("Must supply ", calibrationFile)
816 def diffNonLinearCorrection(self, ccdExposure, dnlLUT, **kwargs):
817 # TODO DM 36636
818 # isrFunctions.diffNonLinearCorrection
819 pass
821 def maskFullAmplifiers(self, ccdExposure, detector, defects, gains=None):
822 """
823 Check for fully masked bad amplifiers and mask them.
825 This includes defects which cover full amplifiers, as well
826 as amplifiers with nan gain values which should be used
827 if self.config.doApplyGains=True.
829 Full defect masking happens later to allow for defects which
830 cross amplifier boundaries.
832 Parameters
833 ----------
834 ccdExposure : `lsst.afw.image.Exposure`
835 Input exposure to be masked.
836 detector : `lsst.afw.cameraGeom.Detector`
837 Detector object.
838 defects : `lsst.ip.isr.Defects`
839 List of defects. Used to determine if an entire
840 amplifier is bad.
841 gains : `dict` [`str`, `float`], optional
842 Dictionary of gains to check if
843 self.config.doApplyGains=True.
845 Returns
846 -------
847 badAmpDict : `str`[`bool`]
848 Dictionary of amplifiers, keyed by name, value is True if
849 amplifier is fully masked.
850 """
851 badAmpDict = {}
853 maskedImage = ccdExposure.getMaskedImage()
855 for amp in detector:
856 ampName = amp.getName()
857 badAmpDict[ampName] = False
859 # Check if entire amp region is defined as a defect
860 # NB: need to use amp.getBBox() for correct comparison with current
861 # defects definition.
862 if defects is not None:
863 badAmpDict[ampName] = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
865 if badAmpDict[ampName]:
866 self.log.warning("Amplifier %s is bad (completely covered with defects)", ampName)
868 if gains is not None and self.config.doApplyGains:
869 if not math.isfinite(gains[ampName]):
870 badAmpDict[ampName] = True
872 self.log.warning("Amplifier %s is bad (non-finite gain)", ampName)
874 # In the case of a bad amp, we will set mask to "BAD"
875 # (here use amp.getRawBBox() for correct association with pixels in
876 # current ccdExposure).
877 if badAmpDict[ampName]:
878 dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
879 afwImage.PARENT)
880 maskView = dataView.getMask()
881 maskView |= maskView.getPlaneBitMask("BAD")
882 del maskView
884 return badAmpDict
886 def checkAllBadAmps(self, badAmpDict, detector):
887 """Check if all amps are marked as bad.
889 Parameters
890 ----------
891 badAmpDict : `str`[`bool`]
892 Dictionary of amplifiers, keyed by name, value is True if
893 amplifier is fully masked.
894 detector : `lsst.afw.cameraGeom.Detector`
895 Detector object.
897 Raises
898 ------
899 UnprocessableDataError if all amps are bad and doCheckUnprocessableData
900 configuration is True.
901 """
902 if not self.config.doCheckUnprocessableData:
903 return
905 for amp in detector:
906 if not badAmpDict.get(amp.getName(), False):
907 return
909 raise UnprocessableDataError(f"All amps in the exposure {detector.getName()} are bad; skipping ISR.")
911 def checkAmpOverscanLevel(self, badAmpDict, exposure, ptc):
912 """Check if the amplifier overscan levels have changed.
914 Any amplifier that has an overscan median level that has changed
915 significantly will be masked as BAD and added to toe badAmpDict.
917 Parameters
918 ----------
919 badAmpDict : `str` [`bool`]
920 Dictionary of amplifiers, keyed by name, value is True if
921 amplifier is fully masked.
922 exposure : `lsst.afw.image.Exposure`
923 Input exposure to be masked (untrimmed).
924 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
925 PTC dataset with gains/read noises.
927 Returns
928 -------
929 badAmpDict : `str`[`bool`]
930 Dictionary of amplifiers, keyed by name.
932 """
933 if isTrimmedExposure(exposure):
934 raise RuntimeError("checkAmpOverscanLevel must be run on an untrimmed exposure.")
936 # We want to consolidate all the amps into one warning if necessary.
937 # This config should not be set to a finite threshold if the necessary
938 # data is not in the PTC.
939 missingWarnString = "No PTC overscan information for amplifier "
940 missingWarnFlag = False
941 for amp in exposure.getDetector():
942 ampName = amp.getName()
944 if not numpy.isfinite(ptc.overscanMedian[ampName]) or \
945 not numpy.isfinite(ptc.overscanMedianSigma[ampName]):
946 missingWarnString += f"{ampName},"
947 missingWarnFlag = True
948 else:
949 key = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
950 # If it is missing, just return the PTC value and it
951 # will be skipped.
952 overscanLevel = exposure.metadata.get(key, ptc.overscanMedian[ampName])
953 pull = (overscanLevel - ptc.overscanMedian[ampName])/ptc.overscanMedianSigma[ampName]
954 if numpy.abs(pull) > self.config.serialOverscanMedianShiftSigmaThreshold:
955 self.log.warning(
956 "Amplifier %s has an overscan level that is %.2f sigma from the expected level; "
957 "masking it as BAD.",
958 ampName,
959 pull,
960 )
962 badAmpDict[ampName] = True
963 exposure.mask[amp.getRawBBox()] |= exposure.mask.getPlaneBitMask("BAD")
965 if missingWarnFlag:
966 self.log.warning(missingWarnString)
968 return badAmpDict
970 def checkAmpNoise(self, badAmpDict, exposure, ptc):
971 """Check if amplifier noise levels are above threshold.
973 Any amplifier that is above the noise level will be masked as BAD
974 and added to the badAmpDict.
976 Parameters
977 ----------
978 badAmpDict : `str` [`bool`]
979 Dictionary of amplifiers, keyed by name, value is True if
980 amplifier is fully masked.
981 exposure : `lsst.afw.image.Exposure`
982 Input exposure to be masked (untrimmed).
983 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
984 PTC dataset with gains/read noises.
986 Returns
987 -------
988 badAmpDict : `str`[`bool`]
989 Dictionary of amplifiers, keyed by name.
990 """
992 if isTrimmedExposure(exposure):
993 raise RuntimeError("checkAmpNoise must be run on an untrimmed exposure.")
995 for amp in exposure.getDetector():
996 ampName = amp.getName()
998 doMask = False
999 if ptc.noise[ampName] > self.config.ampNoiseThreshold:
1000 self.log.info(
1001 "Amplifier %s has a PTC noise level of %.2f e-, above threshold.",
1002 ampName,
1003 ptc.noise[ampName],
1004 )
1005 doMask = True
1006 else:
1007 key = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"
1008 overscanNoise = exposure.metadata.get(key, numpy.nan)
1009 if overscanNoise * ptc.gain[ampName] > self.config.ampNoiseThreshold:
1010 self.log.warning(
1011 "Amplifier %s has an overscan read noise level of %.2f e-, above threshold.",
1012 ampName,
1013 overscanNoise * ptc.gain[ampName],
1014 )
1015 doMask = True
1017 if doMask:
1018 badAmpDict[ampName] = True
1020 exposure.mask[amp.getRawBBox()] |= exposure.mask.getPlaneBitMask("BAD")
1022 return badAmpDict
1024 def maskSaturatedPixels(self, badAmpDict, ccdExposure, detector, detectorConfig, ptc=None):
1025 """
1026 Mask SATURATED and SUSPECT pixels and check if any amplifiers
1027 are fully masked.
1029 Parameters
1030 ----------
1031 badAmpDict : `str` [`bool`]
1032 Dictionary of amplifiers, keyed by name, value is True if
1033 amplifier is fully masked.
1034 ccdExposure : `lsst.afw.image.Exposure`
1035 Input exposure to be masked.
1036 detector : `lsst.afw.cameraGeom.Detector`
1037 Detector object.
1038 defects : `lsst.ip.isr.Defects`
1039 List of defects. Used to determine if an entire
1040 amplifier is bad.
1041 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1042 Per-amplifier configurations.
1043 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
1044 PTC dataset (used if configured to use PTCTURNOFF).
1046 Returns
1047 -------
1048 badAmpDict : `str`[`bool`]
1049 Dictionary of amplifiers, keyed by name.
1050 """
1051 maskedImage = ccdExposure.getMaskedImage()
1053 metadata = ccdExposure.metadata
1055 if self.config.doSaturation and self.config.defaultSaturationSource == "PTCTURNOFF" and ptc is None:
1056 raise RuntimeError("Must provide ptc if using PTCTURNOFF as saturation source.")
1057 if self.config.doSuspect and self.config.defaultSuspectSource == "PTCTURNOFF" and ptc is None:
1058 raise RuntimeError("Must provide ptc if using PTCTURNOFF as suspect source.")
1060 for amp in detector:
1061 ampName = amp.getName()
1063 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
1065 if badAmpDict[ampName]:
1066 # No need to check fully bad amplifiers.
1067 continue
1069 # Mask saturated and suspect pixels.
1070 limits = {}
1071 if self.config.doSaturation:
1072 if self.config.defaultSaturationSource == "PTCTURNOFF":
1073 limits.update({self.config.saturatedMaskName: ptc.ptcTurnoff[amp.getName()]})
1074 elif self.config.defaultSaturationSource == "CAMERAMODEL":
1075 # Set to the default from the camera model.
1076 limits.update({self.config.saturatedMaskName: amp.getSaturation()})
1077 elif self.config.defaultSaturationSource == "NONE":
1078 limits.update({self.config.saturatedMaskName: numpy.inf})
1080 # And update if it is set in the config.
1081 if math.isfinite(ampConfig.saturation):
1082 limits.update({self.config.saturatedMaskName: ampConfig.saturation})
1083 metadata[f"LSST ISR SATURATION LEVEL {ampName}"] = limits[self.config.saturatedMaskName]
1085 if self.config.doSuspect:
1086 if self.config.defaultSuspectSource == "PTCTURNOFF":
1087 limits.update({self.config.suspectMaskName: ptc.ptcTurnoff[amp.getName()]})
1088 elif self.config.defaultSuspectSource == "CAMERAMODEL":
1089 # Set to the default from the camera model.
1090 limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
1091 elif self.config.defaultSuspectSource == "NONE":
1092 limits.update({self.config.suspectMaskName: numpy.inf})
1094 # And update if it set in the config.
1095 if math.isfinite(ampConfig.suspectLevel):
1096 limits.update({self.config.suspectMaskName: ampConfig.suspectLevel})
1097 metadata[f"LSST ISR SUSPECT LEVEL {ampName}"] = limits[self.config.suspectMaskName]
1099 for maskName, maskThreshold in limits.items():
1100 if not math.isnan(maskThreshold):
1101 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
1102 toMask = (dataView.image.array >= maskThreshold)
1103 dataView.mask.array[toMask] |= dataView.mask.getPlaneBitMask(maskName)
1105 # Determine if we've fully masked this amplifier with SUSPECT and
1106 # SAT pixels.
1107 maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
1108 afwImage.PARENT)
1109 maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
1110 self.config.suspectMaskName])
1111 if numpy.all(maskView.getArray() & maskVal > 0):
1112 self.log.warning("Amplifier %s is bad (completely SATURATED or SUSPECT)", ampName)
1113 badAmpDict[ampName] = True
1114 maskView |= maskView.getPlaneBitMask("BAD")
1116 return badAmpDict
1118 def maskITLSatEdgesAndColumns(self, exposure, badAmpDict):
1120 # The following steps will rely on the footprint of saturated
1121 # cores with large areas.
1122 thresh = afwDetection.Threshold(exposure.mask.getPlaneBitMask("SAT"),
1123 afwDetection.Threshold.BITMASK
1124 )
1125 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
1127 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
1128 largeAreas, = numpy.where((satAreas >= self.config.itlEdgeBleedSatMinArea)
1129 & (satAreas < self.config.itlEdgeBleedSatMaxArea))
1131 for largeAreasIndex in largeAreas:
1133 fpCore = fpList[largeAreasIndex]
1135 # Edge bleed masking
1136 if self.config.doITLEdgeBleedMask:
1137 isrFunctions.maskITLEdgeBleed(ccdExposure=exposure,
1138 badAmpDict=badAmpDict,
1139 fpCore=fpCore,
1140 itlEdgeBleedSatMinArea=self.config.itlEdgeBleedSatMinArea,
1141 itlEdgeBleedSatMaxArea=self.config.itlEdgeBleedSatMaxArea,
1142 itlEdgeBleedThreshold=self.config.itlEdgeBleedThreshold,
1143 itlEdgeBleedModelConstant=self.config.itlEdgeBleedModelConstant,
1144 saturatedMaskName=self.config.saturatedMaskName,
1145 log=self.log
1146 )
1147 if self.config.doITLSatSagMask:
1148 isrFunctions.maskITLSatSag(ccdExposure=exposure, fpCore=fpCore,
1149 saturatedMaskName=self.config.saturatedMaskName)
1151 def overscanCorrection(self, mode, detectorConfig, detector, badAmpDict, ccdExposure):
1152 """Apply serial overscan correction in place to all amps.
1154 The actual overscan subtraction is performed by the
1155 `lsst.ip.isr.overscan.OverscanTask`, which is called here.
1157 Parameters
1158 ----------
1159 mode : `str`
1160 Must be `SERIAL` or `PARALLEL`.
1161 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1162 Per-amplifier configurations.
1163 detector : `lsst.afw.cameraGeom.Detector`
1164 Detector object.
1165 badAmpDict : `dict`
1166 Dictionary of amp name to whether it is a bad amp.
1167 ccdExposure : `lsst.afw.image.Exposure`
1168 Exposure to have overscan correction performed.
1170 Returns
1171 -------
1172 overscans : `list` [`lsst.pipe.base.Struct` or None]
1173 Overscan measurements (always in adu).
1174 Each result struct has components:
1176 ``imageFit``
1177 Value or fit subtracted from the amplifier image data.
1178 (scalar or `lsst.afw.image.Image`)
1179 ``overscanFit``
1180 Value or fit subtracted from the overscan image data.
1181 (scalar or `lsst.afw.image.Image`)
1182 ``overscanImage``
1183 Image of the overscan region with the overscan
1184 correction applied. This quantity is used to estimate
1185 the amplifier read noise empirically.
1186 (`lsst.afw.image.Image`)
1187 ``overscanMean``
1188 Mean overscan fit value. (`float`)
1189 ``overscanMedian``
1190 Median overscan fit value. (`float`)
1191 ``overscanSigma``
1192 Clipped standard deviation of the overscan fit. (`float`)
1193 ``residualMean``
1194 Mean of the overscan after fit subtraction. (`float`)
1195 ``residualMedian``
1196 Median of the overscan after fit subtraction. (`float`)
1197 ``residualSigma``
1198 Clipped standard deviation of the overscan after fit
1199 subtraction. (`float`)
1201 See Also
1202 --------
1203 lsst.ip.isr.overscan.OverscanTask
1204 """
1205 if mode not in ["SERIAL", "PARALLEL"]:
1206 raise ValueError("Mode must be SERIAL or PARALLEL")
1208 # This returns a list in amp order, with None for uncorrected amps.
1209 overscans = []
1211 for i, amp in enumerate(detector):
1212 ampName = amp.getName()
1214 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
1216 if mode == "SERIAL" and not ampConfig.doSerialOverscan:
1217 self.log.debug(
1218 "ISR_OSCAN: Amplifier %s/%s configured to skip serial overscan.",
1219 detector.getName(),
1220 ampName,
1221 )
1222 results = None
1223 elif mode == "PARALLEL" and not ampConfig.doParallelOverscan:
1224 self.log.debug(
1225 "ISR_OSCAN: Amplifier %s configured to skip parallel overscan.",
1226 detector.getName(),
1227 ampName,
1228 )
1229 results = None
1230 elif badAmpDict[ampName] or not ccdExposure.getBBox().contains(amp.getBBox()):
1231 results = None
1232 else:
1233 # This check is to confirm that we are not trying to run
1234 # overscan on an already trimmed image.
1235 if isTrimmedExposure(ccdExposure):
1236 self.log.warning(
1237 "ISR_OSCAN: No overscan region for amp %s. Not performing overscan correction.",
1238 ampName,
1239 )
1240 results = None
1241 else:
1242 if mode == "SERIAL":
1243 # We need to set up the subtask here with a custom
1244 # configuration.
1245 serialOverscan = SerialOverscanCorrectionTask(config=ampConfig.serialOverscanConfig)
1246 results = serialOverscan.run(ccdExposure, amp)
1247 else:
1248 config = ampConfig.parallelOverscanConfig
1249 parallelOverscan = ParallelOverscanCorrectionTask(
1250 config=config,
1251 )
1253 metadata = ccdExposure.metadata
1255 # We need to know the saturation level that was used
1256 # for the parallel overscan masking. If it isn't set
1257 # then the configured parallelOverscanSaturationLevel
1258 # will be used instead (assuming
1259 # doParallelOverscanSaturation is True). Note that
1260 # this will have the correct units (adu or electron)
1261 # depending on whether the gain has been applied.
1262 if self.config.doSaturation:
1263 saturationLevel = metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
1264 saturationLevel *= config.parallelOverscanSaturationLevelAdjustmentFactor
1265 else:
1266 saturationLevel = config.parallelOverscanSaturationLevel
1267 if ccdExposure.metadata["LSST ISR UNITS"] == "electron":
1268 # Need to convert to electron from adu.
1269 saturationLevel *= metadata[f"LSST ISR GAIN {amp.getName()}"]
1271 self.log.debug(
1272 "Using saturation level of %.2f for parallel overscan amp %s",
1273 saturationLevel,
1274 amp.getName(),
1275 )
1277 parallelOverscan.maskParallelOverscanAmp(
1278 ccdExposure,
1279 amp,
1280 saturationLevel=saturationLevel,
1281 )
1283 results = parallelOverscan.run(ccdExposure, amp)
1285 metadata = ccdExposure.metadata
1286 keyBase = "LSST ISR OVERSCAN"
1288 # The overscan is always in adu for the serial mode,
1289 # but, it may be electron in the parallel mode if
1290 # doApplyGains==True. If doApplyGains==True, then the
1291 # gains are applied to the untrimmed image, so the
1292 # overscan statistics units here will always match the
1293 # units of the image at this point.
1294 metadata[f"{keyBase} {mode} UNITS"] = ccdExposure.metadata["LSST ISR UNITS"]
1295 metadata[f"{keyBase} {mode} MEAN {ampName}"] = results.overscanMean
1296 metadata[f"{keyBase} {mode} MEDIAN {ampName}"] = results.overscanMedian
1297 metadata[f"{keyBase} {mode} STDEV {ampName}"] = results.overscanSigma
1299 metadata[f"{keyBase} RESIDUAL {mode} MEAN {ampName}"] = results.residualMean
1300 metadata[f"{keyBase} RESIDUAL {mode} MEDIAN {ampName}"] = results.residualMedian
1301 metadata[f"{keyBase} RESIDUAL {mode} STDEV {ampName}"] = results.residualSigma
1303 overscans.append(results)
1305 # Question: should this be finer grained?
1306 ccdExposure.metadata.set("OVERSCAN", "Overscan corrected")
1308 return overscans
1310 def maskNegativeVariance(self, exposure):
1311 """Identify and mask pixels with negative variance values.
1313 Parameters
1314 ----------
1315 exposure : `lsst.afw.image.Exposure`
1316 Exposure to process.
1318 See Also
1319 --------
1320 lsst.ip.isr.isrFunctions.updateVariance
1321 """
1322 maskPlane = exposure.getMask().getPlaneBitMask(self.config.negativeVarianceMaskName)
1323 bad = numpy.where(exposure.getVariance().getArray() <= 0.0)
1324 exposure.mask.array[bad] |= maskPlane
1326 def addVariancePlane(self, exposure, detector):
1327 """Add the variance plane to the image.
1329 The gain and read noise per amp must have been set in the
1330 exposure metadata as ``LSST ISR GAIN ampName`` and
1331 ``LSST ISR READNOISE ampName`` with the units of the image.
1332 Unit conversions for the variance plane will be done as
1333 necessary based on the exposure units.
1335 The units of the variance plane will always be of the same
1336 type as the units of the input image itself
1337 (``LSST ISR UNITS``^2).
1339 Parameters
1340 ----------
1341 exposure : `lsst.afw.image.Exposure`
1342 The exposure to add the variance plane.
1343 detector : `lsst.afw.cameraGeom.Detector`
1344 Detector with geometry info.
1345 """
1346 # NOTE: this will fail if the exposure is not trimmed.
1347 if not isTrimmedExposure(exposure):
1348 raise RuntimeError("Exposure must be trimmed to add variance plane.")
1350 isElectrons = (exposure.metadata["LSST ISR UNITS"] == "electron")
1352 for amp in detector:
1353 if exposure.getBBox().contains(amp.getBBox()):
1354 self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1355 ampExposure = exposure.Factory(exposure, amp.getBBox())
1357 # The effective gain is 1.0 if we are in electron units.
1358 # The metadata read noise is in the same units as the image.
1359 gain = exposure.metadata[f"LSST ISR GAIN {amp.getName()}"] if not isElectrons else 1.0
1360 readNoise = exposure.metadata[f"LSST ISR READNOISE {amp.getName()}"]
1362 isrFunctions.updateVariance(
1363 maskedImage=ampExposure.maskedImage,
1364 gain=gain,
1365 readNoise=readNoise,
1366 replace=False,
1367 )
1369 if self.config.qa is not None and self.config.qa.saveStats is True:
1370 qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1371 afwMath.MEDIAN | afwMath.STDEVCLIP)
1372 self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1373 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1374 qaStats.getValue(afwMath.STDEVCLIP))
1376 if self.config.maskNegativeVariance:
1377 self.maskNegativeVariance(exposure)
1379 def maskDefects(self, exposure, defectBaseList):
1380 """Mask defects using mask plane "BAD", in place.
1382 Parameters
1383 ----------
1384 exposure : `lsst.afw.image.Exposure`
1385 Exposure to process.
1387 defectBaseList : defect-type
1388 List of defects to mask. Can be of type `lsst.ip.isr.Defects`
1389 or `list` of `lsst.afw.image.DefectBase`.
1390 """
1391 maskedImage = exposure.getMaskedImage()
1392 if not isinstance(defectBaseList, Defects):
1393 # Promotes DefectBase to Defect
1394 defectList = Defects(defectBaseList)
1395 else:
1396 defectList = defectBaseList
1397 defectList.maskPixels(maskedImage, maskName="BAD")
1399 if len(self.config.badAmps) == 0:
1400 return
1402 detector = exposure.getDetector()
1403 mask = maskedImage.mask
1404 for badAmp in self.config.badAmps:
1405 if badAmp.startswith(detector.getName()):
1406 # Split on the full detector name plus _, which
1407 # gives us an empty string and the amp name.
1408 ampName = badAmp.split(detector.getName() + "_")[-1]
1409 self.log.info("Masking amplifier %s as bad via config.", ampName)
1410 mask[detector[ampName].getBBox()].array[:, :] |= mask.getPlaneBitMask("BAD")
1412 def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
1413 """Mask edge pixels with applicable mask plane.
1415 Parameters
1416 ----------
1417 exposure : `lsst.afw.image.Exposure`
1418 Exposure to process.
1419 numEdgePixels : `int`, optional
1420 Number of edge pixels to mask.
1421 maskPlane : `str`, optional
1422 Mask plane name to use.
1423 level : `str`, optional
1424 Level at which to mask edges.
1425 """
1426 maskedImage = exposure.getMaskedImage()
1427 maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
1429 if numEdgePixels > 0:
1430 if level == 'DETECTOR':
1431 boxes = [maskedImage.getBBox()]
1432 elif level == 'AMP':
1433 boxes = [amp.getBBox() for amp in exposure.getDetector()]
1435 for box in boxes:
1436 # This makes a bbox numEdgeSuspect pixels smaller than the
1437 # image on each side
1438 subImage = maskedImage[box]
1439 box.grow(-numEdgePixels)
1440 # Mask pixels outside box
1441 SourceDetectionTask.setEdgeBits(
1442 subImage,
1443 box,
1444 maskBitMask)
1446 def maskNan(self, exposure):
1447 """Mask NaNs using mask plane "UNMASKEDNAN", in place.
1449 Parameters
1450 ----------
1451 exposure : `lsst.afw.image.Exposure`
1452 Exposure to process.
1454 Notes
1455 -----
1456 We mask over all non-finite values (NaN, inf), including those
1457 that are masked with other bits (because those may or may not be
1458 interpolated over later, and we want to remove all NaN/infs).
1459 Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
1460 preserve the historical name.
1461 """
1462 maskedImage = exposure.getMaskedImage()
1464 # Find and mask NaNs
1465 maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
1466 numNans = maskNans(maskedImage, maskVal)
1467 self.metadata["NUMNANS"] = numNans
1468 if numNans > 0:
1469 self.log.warning("There were %d unmasked NaNs.", numNans)
1471 def setBadRegions(self, exposure):
1472 """Set bad regions from large contiguous regions.
1474 Parameters
1475 ----------
1476 exposure : `lsst.afw.Exposure`
1477 Exposure to set bad regions.
1479 Notes
1480 -----
1481 Reset and interpolate bad pixels.
1483 Large contiguous bad regions (which should have the BAD mask
1484 bit set) should have their values set to the image median.
1485 This group should include defects and bad amplifiers. As the
1486 area covered by these defects are large, there's little
1487 reason to expect that interpolation would provide a more
1488 useful value.
1490 Smaller defects can be safely interpolated after the larger
1491 regions have had their pixel values reset. This ensures
1492 that the remaining defects adjacent to bad amplifiers (as an
1493 example) do not attempt to interpolate extreme values.
1494 """
1495 badPixelCount, badPixelValue = isrFunctions.setBadRegions(exposure)
1496 if badPixelCount > 0:
1497 self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1499 @contextmanager
1500 def flatContext(self, exp, flat, dark=None):
1501 """Context manager that applies and removes flats and darks,
1502 if the task is configured to apply them.
1504 Parameters
1505 ----------
1506 exp : `lsst.afw.image.Exposure`
1507 Exposure to process.
1508 flat : `lsst.afw.image.Exposure`
1509 Flat exposure the same size as ``exp``.
1510 dark : `lsst.afw.image.Exposure`, optional
1511 Dark exposure the same size as ``exp``.
1513 Yields
1514 ------
1515 exp : `lsst.afw.image.Exposure`
1516 The flat and dark corrected exposure.
1517 """
1518 if self.config.doDark and dark is not None:
1519 self.darkCorrection(exp, dark)
1520 if self.config.doFlat and flat is not None:
1521 self.flatCorrection(exp, flat)
1522 try:
1523 yield exp
1524 finally:
1525 if self.config.doFlat and flat is not None:
1526 self.flatCorrection(exp, flat, invert=True)
1527 if self.config.doDark and dark is not None:
1528 self.darkCorrection(exp, dark, invert=True)
1530 def getBrighterFatterKernel(self, detector, bfKernel):
1531 detName = detector.getName()
1533 # This is expected to be a dictionary of amp-wise gains.
1534 bfGains = bfKernel.gain
1535 if bfKernel.level == 'DETECTOR':
1536 if detName in bfKernel.detKernels:
1537 bfKernelOut = bfKernel.detKernels[detName]
1538 return bfKernelOut, bfGains
1539 else:
1540 raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1541 elif bfKernel.level == 'AMP':
1542 self.log.info("Making DETECTOR level kernel from AMP based brighter "
1543 "fatter kernels.")
1544 bfKernel.makeDetectorKernelFromAmpwiseKernels(detName)
1545 bfKernelOut = bfKernel.detKernels[detName]
1546 return bfKernelOut, bfGains
1548 def applyElectrostaticBrighterFatterCorrection(self, ccdExposure, flat, dark, electroBfDistortionMatrix,
1549 brighterFatterApplyGain, bfGains):
1550 """Apply an electrostatic brighter fatter correction to the image
1551 using the method defined in Astier et al. 2023.
1553 Note that this correction requires that the image is in units
1554 electrons.
1556 Parameters
1557 ----------
1558 ccdExposure : `lsst.afw.image.Exposure`
1559 Exposure to process.
1560 flat : `lsst.afw.image.Exposure`
1561 Flat exposure the same size as ``exp``.
1562 dark : `lsst.afw.image.Exposure`, optional
1563 Dark exposure the same size as ``exp``.
1564 electroBfDistortionMatrix : `lsst.ip.isr.ElectrostaticBrighterFatter`
1565 The brighter-fatter kernel.
1566 brighterFatterApplyGain : `bool`
1567 Apply the gain to convert the image to electrons?
1568 bfGains : `dict`
1569 The gains to use if brighterFatterApplyGain = True.
1571 Yields
1572 ------
1573 exp : `lsst.afw.image.Exposure`
1574 The flat and dark corrected exposure.
1575 """
1576 interpExp = ccdExposure.clone()
1577 method = self.config.brighterFatterCorrectionMethod
1579 # We need to interpolate before we do B-F. Note that
1580 # brighterFatterFwhmForInterpolation is currently unused.
1581 isrFunctions.interpolateFromMask(
1582 maskedImage=interpExp.getMaskedImage(),
1583 fwhm=self.config.brighterFatterFwhmForInterpolation,
1584 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1585 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1586 useLegacyInterp=self.config.useLegacyInterp,
1587 )
1588 bfExp = interpExp.clone()
1590 ccdExposure = electrostaticBrighterFatterCorrection(
1591 self.log,
1592 bfExp,
1593 electroBfDistortionMatrix,
1594 brighterFatterApplyGain,
1595 bfGains,
1596 applyFilterCorrection=(method == "ASTIER23+FILTERCORRECTION"),
1597 )
1599 # Applying the brighter-fatter correction applies a
1600 # convolution to the science image. At the edges this
1601 # convolution may not have sufficient valid pixels to
1602 # produce a valid correction. Mark pixels within the size
1603 # of the brighter-fatter kernel as EDGE to warn of this
1604 # fact.
1605 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1606 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(electroBfDistortionMatrix.aMatrix.shape) // 2,
1607 maskPlane="EDGE")
1609 if self.config.brighterFatterMaskGrowSize > 0:
1610 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1611 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1612 isrFunctions.growMasks(ccdExposure.getMask(),
1613 radius=self.config.brighterFatterMaskGrowSize,
1614 maskNameList=maskPlane,
1615 maskValue=maskPlane)
1617 return ccdExposure
1619 def applyFluxConservingBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel,
1620 brighterFatterApplyGain, bfGains):
1621 """Apply a brighter fatter correction to the image using the
1622 method defined in Coulton et al. 2019 with flux-conserving
1623 corrections.
1625 Note that this correction requires that the image is in units
1626 electrons.
1628 Parameters
1629 ----------
1630 ccdExposure : `lsst.afw.image.Exposure`
1631 Exposure to process.
1632 flat : `lsst.afw.image.Exposure`
1633 Flat exposure the same size as ``exp``.
1634 dark : `lsst.afw.image.Exposure`, optional
1635 Dark exposure the same size as ``exp``.
1636 bfKernel : `lsst.ip.isr.BrighterFatterKernel`
1637 The brighter-fatter kernel.
1638 brighterFatterApplyGain : `bool`
1639 Apply the gain to convert the image to electrons?
1640 bfGains : `dict`
1641 The gains to use if brighterFatterApplyGain = True.
1643 Yields
1644 ------
1645 exp : `lsst.afw.image.Exposure`
1646 The flat and dark corrected exposure.
1647 """
1648 interpExp = ccdExposure.clone()
1650 # We need to interpolate before we do B-F. Note that
1651 # brighterFatterFwhmForInterpolation is currently unused.
1652 isrFunctions.interpolateFromMask(
1653 maskedImage=interpExp.getMaskedImage(),
1654 fwhm=self.config.brighterFatterFwhmForInterpolation,
1655 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1656 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1657 useLegacyInterp=self.config.useLegacyInterp,
1658 )
1660 bfExp = interpExp.clone()
1662 bfResults = fluxConservingBrighterFatterCorrection(
1663 bfExp, bfKernel,
1664 self.config.brighterFatterMaxIter,
1665 self.config.brighterFatterThreshold,
1666 brighterFatterApplyGain,
1667 bfGains,
1668 )
1670 bfCorrIters = bfResults[1]
1671 if bfCorrIters == self.config.brighterFatterMaxIter:
1672 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1673 bfResults[0])
1674 else:
1675 self.log.info("Finished brighter-fatter correction in %d iterations.",
1676 bfResults[1])
1678 image = ccdExposure.getMaskedImage().getImage()
1679 bfCorr = bfExp.getMaskedImage().getImage()
1680 bfCorr -= interpExp.getMaskedImage().getImage()
1681 image += bfCorr
1683 # Applying the brighter-fatter correction applies a
1684 # convolution to the science image. At the edges this
1685 # convolution may not have sufficient valid pixels to
1686 # produce a valid correction. Mark pixels within the size
1687 # of the brighter-fatter kernel as EDGE to warn of this
1688 # fact.
1689 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1690 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1691 maskPlane="EDGE")
1693 if self.config.brighterFatterMaskGrowSize > 0:
1694 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1695 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1696 isrFunctions.growMasks(ccdExposure.getMask(),
1697 radius=self.config.brighterFatterMaskGrowSize,
1698 maskNameList=maskPlane,
1699 maskValue=maskPlane)
1701 # Set the metadata here.
1702 ccdExposure.metadata["LSST ISR BF ITERS"] = bfCorrIters
1704 return ccdExposure
1706 def applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, brighterFatterApplyGain,
1707 bfGains):
1708 """Apply a brighter fatter correction to the image using the
1709 method defined in Coulton et al. 2019.
1711 Note that this correction requires that the image is in units
1712 electrons.
1714 Parameters
1715 ----------
1716 ccdExposure : `lsst.afw.image.Exposure`
1717 Exposure to process.
1718 flat : `lsst.afw.image.Exposure`
1719 Flat exposure the same size as ``exp``.
1720 dark : `lsst.afw.image.Exposure`, optional
1721 Dark exposure the same size as ``exp``.
1722 bfKernel : `lsst.ip.isr.BrighterFatterKernel`
1723 The brighter-fatter kernel.
1724 brighterFatterApplyGain : `bool`
1725 Apply the gain to convert the image to electrons?
1726 bfGains : `dict`
1727 The gains to use if brighterFatterApplyGain = True.
1729 Yields
1730 ------
1731 exp : `lsst.afw.image.Exposure`
1732 The flat and dark corrected exposure.
1733 """
1734 interpExp = ccdExposure.clone()
1736 # We need to interpolate before we do B-F. Note that
1737 # brighterFatterFwhmForInterpolation is currently unused.
1738 isrFunctions.interpolateFromMask(
1739 maskedImage=interpExp.getMaskedImage(),
1740 fwhm=self.config.brighterFatterFwhmForInterpolation,
1741 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1742 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1743 useLegacyInterp=self.config.useLegacyInterp,
1744 )
1745 bfExp = interpExp.clone()
1746 bfResults = brighterFatterCorrection(bfExp, bfKernel,
1747 self.config.brighterFatterMaxIter,
1748 self.config.brighterFatterThreshold,
1749 brighterFatterApplyGain,
1750 bfGains)
1751 bfCorrIters = bfResults[1]
1752 if bfCorrIters == self.config.brighterFatterMaxIter:
1753 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1754 bfResults[0])
1755 else:
1756 self.log.info("Finished brighter-fatter correction in %d iterations.",
1757 bfResults[1])
1759 image = ccdExposure.getMaskedImage().getImage()
1760 bfCorr = bfExp.getMaskedImage().getImage()
1761 bfCorr -= interpExp.getMaskedImage().getImage()
1762 image += bfCorr
1764 # Applying the brighter-fatter correction applies a
1765 # convolution to the science image. At the edges this
1766 # convolution may not have sufficient valid pixels to
1767 # produce a valid correction. Mark pixels within the size
1768 # of the brighter-fatter kernel as EDGE to warn of this
1769 # fact.
1770 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1771 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1772 maskPlane="EDGE")
1774 if self.config.brighterFatterMaskGrowSize > 0:
1775 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1776 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1777 isrFunctions.growMasks(ccdExposure.getMask(),
1778 radius=self.config.brighterFatterMaskGrowSize,
1779 maskNameList=maskPlane,
1780 maskValue=maskPlane)
1782 # Set the metadata here.
1783 ccdExposure.metadata["LSST ISR BF ITERS"] = bfCorrIters
1785 return ccdExposure
1787 def darkCorrection(self, exposure, darkExposure, invert=False):
1788 """Apply dark correction in place.
1790 Parameters
1791 ----------
1792 exposure : `lsst.afw.image.Exposure`
1793 Exposure to process.
1794 darkExposure : `lsst.afw.image.Exposure`
1795 Dark exposure of the same size as ``exposure``.
1796 invert : `Bool`, optional
1797 If True, re-add the dark to an already corrected image.
1799 Raises
1800 ------
1801 RuntimeError
1802 Raised if either ``exposure`` or ``darkExposure`` do not
1803 have their dark time defined.
1805 See Also
1806 --------
1807 lsst.ip.isr.isrFunctions.darkCorrection
1808 """
1809 expScale = exposure.visitInfo.darkTime
1810 if math.isnan(expScale):
1811 raise RuntimeError("Exposure darktime is NAN.")
1812 if darkExposure.visitInfo is not None \
1813 and not math.isnan(darkExposure.visitInfo.darkTime):
1814 darkScale = darkExposure.visitInfo.darkTime
1815 else:
1816 # DM-17444: darkExposure.visitInfo is None
1817 # so darkTime does not exist.
1818 self.log.warning("darkExposure.visitInfo does not exist. Using darkScale = 1.0.")
1819 darkScale = 1.0
1821 isrFunctions.darkCorrection(
1822 maskedImage=exposure.maskedImage,
1823 darkMaskedImage=darkExposure.maskedImage,
1824 expScale=expScale,
1825 darkScale=darkScale,
1826 invert=invert,
1827 )
1829 @staticmethod
1830 def extractCalibDate(calib):
1831 """Extract common calibration metadata values that will be written to
1832 output header.
1834 Parameters
1835 ----------
1836 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1837 Calibration to pull date information from.
1839 Returns
1840 -------
1841 dateString : `str`
1842 Calibration creation date string to add to header.
1843 """
1844 if hasattr(calib, "getMetadata"):
1845 if 'CALIB_CREATION_DATE' in calib.metadata:
1846 return " ".join((calib.metadata.get("CALIB_CREATION_DATE", "Unknown"),
1847 calib.metadata.get("CALIB_CREATION_TIME", "Unknown")))
1848 else:
1849 return " ".join((calib.metadata.get("CALIB_CREATE_DATE", "Unknown"),
1850 calib.metadata.get("CALIB_CREATE_TIME", "Unknown")))
1851 else:
1852 return "Unknown Unknown"
1854 def compareUnits(self, calibMetadata, calibName):
1855 """Compare units from calibration to ISR units.
1857 This compares calibration units (adu or electron) to whether
1858 doApplyGain is set.
1860 Parameters
1861 ----------
1862 calibMetadata : `lsst.daf.base.PropertyList`
1863 Calibration metadata from header.
1864 calibName : `str`
1865 Calibration name for log message.
1866 """
1867 calibUnits = calibMetadata.get("LSST ISR UNITS", "adu")
1868 isrUnits = "electron" if self.config.doApplyGains else "adu"
1869 if calibUnits != isrUnits:
1870 if self.config.doRaiseOnCalibMismatch:
1871 raise RuntimeError(
1872 "Unit mismatch: isr has %s units but %s has %s units",
1873 isrUnits,
1874 calibName,
1875 calibUnits,
1876 )
1877 else:
1878 self.log.warning(
1879 "Unit mismatch: isr has %s units but %s has %s units",
1880 isrUnits,
1881 calibName,
1882 calibUnits,
1883 )
1885 def convertIntToFloat(self, exposure):
1886 """Convert exposure image from uint16 to float.
1888 If the exposure does not need to be converted, the input is
1889 immediately returned. For exposures that are converted to use
1890 floating point pixels, the variance is set to unity and the
1891 mask to zero.
1893 Parameters
1894 ----------
1895 exposure : `lsst.afw.image.Exposure`
1896 The raw exposure to be converted.
1898 Returns
1899 -------
1900 newexposure : `lsst.afw.image.Exposure`
1901 The input ``exposure``, converted to floating point pixels.
1903 Raises
1904 ------
1905 RuntimeError
1906 Raised if the exposure type cannot be converted to float.
1908 """
1909 if isinstance(exposure, afwImage.ExposureF):
1910 # Nothing to be done
1911 self.log.debug("Exposure already of type float.")
1912 return exposure
1913 if not hasattr(exposure, "convertF"):
1914 raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
1916 newexposure = exposure.convertF()
1917 newexposure.variance[:] = 1
1918 newexposure.mask[:] = 0x0
1920 return newexposure
1922 def ditherCounts(self, exposure, detectorConfig, fallbackSeed=12345):
1923 """Dither the counts in the exposure.
1925 Parameters
1926 ----------
1927 exposure : `lsst.afw.image.Exposure`
1928 The raw exposure to be dithered.
1929 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1930 Configuration for overscan/etc for this detector.
1931 fallbackSeed : `int`, optional
1932 Random seed to fall back to if exposure.getInfo().getId() is
1933 not set.
1934 """
1935 if detectorConfig.integerDitherMode == "NONE":
1936 # Nothing to do here.
1937 return
1939 # This ID is a unique combination of {exposure, detector} for a raw
1940 # image as we have here. We additionally need to take the lower
1941 # 32 bits to be used as a random seed.
1942 if exposure.info.id is not None:
1943 seed = exposure.info.id & 0xFFFFFFFF
1944 else:
1945 seed = fallbackSeed
1946 self.log.warning("No exposure ID found; using fallback random seed.")
1948 self.log.info("Seeding dithering random number generator with %d.", seed)
1949 rng = numpy.random.RandomState(seed=seed)
1951 if detectorConfig.integerDitherMode == "POSITIVE":
1952 low = 0.0
1953 high = 1.0
1954 elif detectorConfig.integerDitherMode == "NEGATIVE":
1955 low = -1.0
1956 high = 0.0
1957 elif detectorConfig.integerDitherMode == "SYMMETRIC":
1958 low = -0.5
1959 high = 0.5
1960 else:
1961 raise RuntimeError("Invalid config")
1963 exposure.image.array[:, :] += rng.uniform(low=low, high=high, size=exposure.image.array.shape)
1965 def checkBssVoltage(self, exposure):
1966 """Check the back-side bias voltage to see if the detector is on.
1968 Parameters
1969 ----------
1970 exposure : `lsst.afw.image.ExposureF`
1971 Input exposure.
1973 Raises
1974 ------
1975 `UnprocessableDataError` if voltage is off.
1976 """
1977 voltage = exposure.metadata.get(self.config.bssVoltageKeyword, None)
1978 if voltage is None or not numpy.isfinite(voltage):
1979 self.log.warning(
1980 "Back-side bias voltage %s not found in metadata.",
1981 self.config.bssVoltageKeyword,
1982 )
1983 return
1985 hv = exposure.metadata.get(self.config.hvBiasKeyword, None)
1986 if hv is None:
1987 self.log.warning(
1988 "HV bias on %s not found in metadata.",
1989 self.config.hvBiasKeyword,
1990 )
1991 return
1993 if voltage < self.config.bssVoltageMinimum or hv == "OFF":
1994 detector = exposure.getDetector()
1995 raise UnprocessableDataError(
1996 f"Back-side bias voltage is turned off for {detector.getName()}; skipping ISR.",
1997 )
1999 @deprecated(
2000 reason=(
2001 "makeBinnedImages is no longer used. "
2002 "Please subtask lsst.ip.isr.BinImageDataTask instead."
2003 ),
2004 version="v28", category=FutureWarning
2005 )
2006 def makeBinnedImages(self, exposure):
2007 """Make visualizeVisit style binned exposures.
2009 Parameters
2010 ----------
2011 exposure : `lsst.afw.image.Exposure`
2012 Exposure to bin.
2014 Returns
2015 -------
2016 bin1 : `lsst.afw.image.Exposure`
2017 Binned exposure using binFactor1.
2018 bin2 : `lsst.afw.image.Exposure`
2019 Binned exposure using binFactor2.
2020 """
2021 mi = exposure.getMaskedImage()
2023 bin1 = afwMath.binImage(mi, self.config.binFactor1)
2024 bin2 = afwMath.binImage(mi, self.config.binFactor2)
2026 bin1 = afwImage.makeExposure(bin1)
2027 bin2 = afwImage.makeExposure(bin2)
2029 bin1.setInfo(exposure.getInfo())
2030 bin2.setInfo(exposure.getInfo())
2032 return bin1, bin2
2034 def run(
2035 self,
2036 ccdExposure,
2037 *,
2038 dnlLUT=None,
2039 bias=None,
2040 deferredChargeCalib=None,
2041 linearizer=None,
2042 ptc=None,
2043 gainCorrection=None,
2044 crosstalk=None,
2045 defects=None,
2046 bfKernel=None,
2047 electroBfDistortionMatrix=None,
2048 dark=None,
2049 flat=None,
2050 camera=None,
2051 ):
2052 """Run the IsrTaskLSST task.
2054 Parameters
2055 ----------
2056 ccdExposure : `lsst.afw.image.Exposure`
2057 Exposure to run ISR.
2058 dnlLUT : `None`, optional
2059 DNL lookup table; placeholder, unused.
2060 bias : `lsst.afw.image.Exposure`, optional
2061 Bias frame.
2062 deferredChargeCalib : `lsst.ip.isr.DeferredChargeCalib`, optional
2063 Deferred charge calibration.
2064 linearizer : `lsst.ip.isr.Linearizer`, optional
2065 Linearizer calibration.
2066 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
2067 PTC dataset.
2068 gainCorrection : `lsst.ip.isr.GainCorrection`, optional
2069 Gain correction dataset.
2070 crosstalk : `lsst.ip.isr.CrosstalkCalib`, optional
2071 Crosstalk calibration dataset.
2072 defects : `lsst.ip.isr.Defects`, optional
2073 Defects dataset.
2074 bfKernel : `lsst.ip.isr.BrighterFatterKernel`, optional
2075 Brighter-fatter kernel dataset.
2076 dark : `lsst.afw.image.Exposure`, optional
2077 Dark frame.
2078 flat : `lsst.afw.image.Exposure`, optional
2079 Flat-field frame.
2080 camera : `lsst.afw.cameraGeom.Camera`, optional
2081 Camera object.
2083 Returns
2084 -------
2085 result : `lsst.pipe.base.Struct`
2086 Struct with fields:
2087 ``exposure``: `lsst.afw.image.Exposure`
2088 Calibrated exposure.
2089 ``outputBin1Exposure``: `lsst.afw.image.Exposure`
2090 Binned exposure (bin1 config).
2091 ``outputBin2Exposure``: `lsst.afw.image.Exposure`
2092 Binned exposure (bin2 config).
2093 ``outputExposure``: `lsst.afw.image.Exposure`
2094 Calibrated exposure (same as ``exposure``).
2095 ``outputStatistics``: `lsst.ip.isr.isrStatistics`
2096 Calibrated exposure statistics.
2097 """
2098 detector = ccdExposure.getDetector()
2100 overscanDetectorConfig = self.config.overscanCamera.getOverscanDetectorConfig(detector)
2102 if self.config.doBootstrap:
2103 if ptc is not None:
2104 self.log.warning("Task configured with doBootstrap=True. Ignoring provided PTC.")
2105 ptc = None
2106 else:
2107 if self.config.useGainsFrom == "LINEARIZER":
2108 if linearizer is None:
2109 raise RuntimeError("doBootstrap==False and useGainsFrom == 'LINEARIZER' but "
2110 "no linearizer provided.")
2111 elif self.config.useGainsFrom == "PTC":
2112 if ptc is None:
2113 raise RuntimeError("doBootstrap==False and useGainsFrom == 'PTC' but no PTC provided.")
2115 # Validation step: check inputs match exposure configuration.
2116 exposureMetadata = ccdExposure.metadata
2117 doRaise = self.config.doRaiseOnCalibMismatch
2118 keywords = self.config.cameraKeywordsToCompare
2119 if not self.config.doBootstrap:
2120 if self.config.useGainsFrom == "LINEARIZER":
2121 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer,
2122 "LINEARIZER", log=self.log)
2123 elif self.config.useGainsFrom == "PTC":
2124 compareCameraKeywords(doRaise, keywords, exposureMetadata, ptc, "PTC",
2125 log=self.log)
2126 # Note that doCorrectGains can be True without a gainCorrection.
2127 if self.config.doCorrectGains and gainCorrection is not None:
2128 compareCameraKeywords(
2129 doRaise,
2130 keywords,
2131 exposureMetadata,
2132 gainCorrection,
2133 "gain_correction",
2134 log=self.log,
2135 )
2136 else:
2137 if self.config.doCorrectGains:
2138 raise RuntimeError("doCorrectGains is True but no ptc provided.")
2139 if self.config.doDiffNonLinearCorrection:
2140 if dnlLUT is None:
2141 raise RuntimeError("doDiffNonLinearCorrection is True but no dnlLUT provided.")
2142 compareCameraKeywords(doRaise, keywords, exposureMetadata, dnlLUT, "dnlLUT", log=self.log)
2143 if self.config.doLinearize:
2144 if linearizer is None:
2145 raise RuntimeError("doLinearize is True but no linearizer provided.")
2146 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer, "linearizer", log=self.log)
2147 if self.config.doBias:
2148 if bias is None:
2149 raise RuntimeError("doBias is True but no bias provided.")
2150 compareCameraKeywords(doRaise, keywords, exposureMetadata, bias, "bias", log=self.log)
2151 self.compareUnits(bias.metadata, "bias")
2152 if self.config.doCrosstalk:
2153 if crosstalk is None:
2154 raise RuntimeError("doCrosstalk is True but no crosstalk provided.")
2155 compareCameraKeywords(doRaise, keywords, exposureMetadata, crosstalk, "crosstalk", log=self.log)
2156 if self.config.doDeferredCharge:
2157 if deferredChargeCalib is None:
2158 raise RuntimeError("doDeferredCharge is True but no deferredChargeCalib provided.")
2159 compareCameraKeywords(
2160 doRaise,
2161 keywords,
2162 exposureMetadata,
2163 deferredChargeCalib,
2164 "CTI",
2165 log=self.log,
2166 )
2167 if self.config.doDefect:
2168 if defects is None:
2169 raise RuntimeError("doDefect is True but no defects provided.")
2170 compareCameraKeywords(doRaise, keywords, exposureMetadata, defects, "defects", log=self.log)
2171 if self.config.doDark:
2172 if dark is None:
2173 raise RuntimeError("doDark is True but no dark frame provided.")
2174 compareCameraKeywords(doRaise, keywords, exposureMetadata, dark, "dark", log=self.log)
2175 self.compareUnits(bias.metadata, "dark")
2176 if self.config.doBrighterFatter:
2177 if self.config.brighterFatterCorrectionMethod in ["ASTIER23", "ASTIER23+FILTERCORRECTION"]:
2178 if electroBfDistortionMatrix is None:
2179 raise RuntimeError("Must supply an electroBfDistortionMatrix if BF "
2180 "correction method is ASTIER23*.")
2181 compareCameraKeywords(doRaise, keywords, exposureMetadata,
2182 electroBfDistortionMatrix, "bf", log=self.log)
2183 elif self.config.brighterFatterCorrectionMethod in ["COULTON18", "COULTON18_FLUX_CONSERVING"]:
2184 if bfKernel is None:
2185 raise RuntimeError("Must supply an kernel if BF correction method is COULTON*.")
2186 compareCameraKeywords(doRaise, keywords, exposureMetadata, bfKernel, "bf", log=self.log)
2187 else:
2188 raise ValueError("%s is not a known brighter-fatter correction "
2189 "method." % self.config.brighterFatterCorrectionMethod)
2190 if self.config.doFlat:
2191 if flat is None:
2192 raise RuntimeError("doFlat is True but no flat provided.")
2193 compareCameraKeywords(doRaise, keywords, exposureMetadata, flat, "flat", log=self.log)
2195 if self.config.doSaturation:
2196 if self.config.defaultSaturationSource in ["PTCTURNOFF",]:
2197 if ptc is None:
2198 raise RuntimeError(
2199 "doSaturation is True and defaultSaturationSource is "
2200 f"{self.config.defaultSaturationSource}, but no ptc provided."
2201 )
2202 if self.config.doSuspect:
2203 if self.config.defaultSuspectSource in ["PTCTURNOFF",]:
2204 if ptc is None:
2205 raise RuntimeError(
2206 "doSuspect is True and defaultSuspectSource is "
2207 f"{self.config.defaultSuspectSource}, but no ptc provided."
2208 )
2210 if self.config.doCheckUnprocessableData and self.config.bssVoltageMinimum > 0.0:
2211 self.checkBssVoltage(ccdExposure)
2213 # FIXME: Make sure that if linearity is done then it is matched
2214 # with the right PTC.
2216 # We keep track of units: start in adu.
2217 exposureMetadata["LSST ISR UNITS"] = "adu"
2218 exposureMetadata["LSST ISR GAINCORRECTION APPLIED"] = False
2219 exposureMetadata["LSST ISR CROSSTALK APPLIED"] = False
2220 exposureMetadata["LSST ISR OVERSCANLEVEL CHECKED"] = False
2221 exposureMetadata["LSST ISR NOISE CHECKED"] = False
2222 exposureMetadata["LSST ISR LINEARIZER APPLIED"] = False
2223 exposureMetadata["LSST ISR CTI APPLIED"] = False
2224 exposureMetadata["LSST ISR BIAS APPLIED"] = False
2225 exposureMetadata["LSST ISR DARK APPLIED"] = False
2226 exposureMetadata["LSST ISR BF APPLIED"] = False
2227 exposureMetadata["LSST ISR FLAT APPLIED"] = False
2228 exposureMetadata["LSST ISR DEFECTS APPLIED"] = False
2230 if self.config.doBootstrap:
2231 self.log.info("Configured using doBootstrap=True; using gain of 1.0 (adu units)")
2232 ptc = PhotonTransferCurveDataset([amp.getName() for amp in detector], "NOMINAL_PTC", 1)
2233 for amp in detector:
2234 ptc.gain[amp.getName()] = 1.0
2235 ptc.noise[amp.getName()] = 0.0
2236 ptc.ptcTurnoff[amp.getName()] = numpy.inf
2237 elif self.config.useGainsFrom == "LINEARIZER":
2238 self.log.info("Using gains from linearizer.")
2239 # Create a dummy ptc object to hold the gains from the linearizer.
2240 ptc = PhotonTransferCurveDataset([amp.getName() for amp in detector], "NOMINAL_PTC", 1)
2241 for amp in detector:
2242 ptc.gain[amp.getName()] = linearizer.inputGain[amp.getName()]
2243 ptc.noise[amp.getName()] = 0.0
2244 ptc.ptcTurnoff[amp.getName()] = linearizer.inputTurnoff[amp.getName()]
2246 exposureMetadata["LSST ISR BOOTSTRAP"] = self.config.doBootstrap
2248 # Choose the gains to use
2249 gains = ptc.gain
2251 # And check if we have configured gains to override. This is
2252 # also a warning, since it should not be typical usage.
2253 for amp in detector:
2254 if not math.isnan(gain := overscanDetectorConfig.getOverscanAmpConfig(amp).gain):
2255 gains[amp.getName()] = gain
2256 self.log.warning(
2257 "Overriding gain for amp %s with configured value of %.3f.",
2258 amp.getName(),
2259 gain,
2260 )
2262 # First we convert the exposure to floating point values
2263 # (if necessary).
2264 self.log.debug("Converting exposure to floating point values.")
2265 ccdExposure = self.convertIntToFloat(ccdExposure)
2267 # Then we mark which amplifiers are completely bad from defects.
2268 badAmpDict = self.maskFullAmplifiers(ccdExposure, detector, defects, gains=gains)
2270 self.checkAllBadAmps(badAmpDict, detector)
2272 # Now we go through ISR steps.
2274 # Differential non-linearity correction.
2275 # Units: adu
2276 if self.config.doDiffNonLinearCorrection:
2277 self.diffNonLinearCorrection(ccdExposure, dnlLUT)
2279 # Dither the integer counts.
2280 # Input units: integerized adu
2281 # Output units: floating-point adu
2282 self.ditherCounts(ccdExposure, overscanDetectorConfig)
2284 # Serial overscan correction.
2285 # Input units: adu
2286 # Output units: adu
2287 if overscanDetectorConfig.doAnySerialOverscan:
2288 serialOverscans = self.overscanCorrection(
2289 "SERIAL",
2290 overscanDetectorConfig,
2291 detector,
2292 badAmpDict,
2293 ccdExposure,
2294 )
2296 if self.config.doBootstrap or self.config.useGainsFrom == "LINEARIZER":
2297 # Get the empirical read noise
2298 for amp, serialOverscan in zip(detector, serialOverscans):
2299 if serialOverscan is None:
2300 ptc.noise[amp.getName()] = 0.0
2301 else:
2302 # All PhotonTransferCurveDataset objects should contain
2303 # noise attributes in units of electrons. The read
2304 # noise measured from overscans is always in adu, so we
2305 # scale it by the gain.
2306 # Note that in bootstrap mode, these gains will always
2307 # be 1.0, but we put this conversion here for clarity.
2308 ptc.noise[amp.getName()] = serialOverscan.residualSigma * gains[amp.getName()]
2309 else:
2310 serialOverscans = [None]*len(detector)
2312 # After serial overscan correction, we can mask SATURATED and
2313 # SUSPECT pixels. This updates badAmpDict if any amplifier
2314 # is fully saturated after serial overscan correction.
2316 # The saturation is currently assumed to be recorded in
2317 # overscan-corrected adu.
2318 badAmpDict = self.maskSaturatedPixels(
2319 badAmpDict,
2320 ccdExposure,
2321 detector,
2322 overscanDetectorConfig,
2323 ptc=ptc,
2324 )
2326 self.checkAllBadAmps(badAmpDict, detector)
2328 if self.config.doCorrectGains and gainCorrection is not None:
2329 self.log.info("Correcting gains based on input GainCorrection.")
2330 gainCorrection.correctGains(gains, exposure=ccdExposure)
2331 exposureMetadata["LSST ISR GAINCORRECTION APPLIED"] = True
2332 elif self.config.doCorrectGains:
2333 self.log.info("Skipping gain correction because no GainCorrection available.")
2335 # Do gain normalization.
2336 # Input units: adu
2337 # Output units: electron
2338 if self.config.doApplyGains:
2339 self.log.info("Using gain values to convert from adu to electron units.")
2340 isrFunctions.applyGains(ccdExposure, normalizeGains=False, ptcGains=gains, isTrimmed=False)
2341 # The units are now electron.
2342 exposureMetadata["LSST ISR UNITS"] = "electron"
2344 # Update the saturation units in the metadata if there.
2345 # These will always have the same units as the image.
2346 for amp in detector:
2347 ampName = amp.getName()
2348 if (key := f"LSST ISR SATURATION LEVEL {ampName}") in exposureMetadata:
2349 exposureMetadata[key] *= gains[ampName]
2350 if (key := f"LSST ISR SUSPECT LEVEL {ampName}") in exposureMetadata:
2351 exposureMetadata[key] *= gains[ampName]
2353 # Record gain and read noise in header.
2354 metadata = ccdExposure.metadata
2355 metadata["LSST ISR READNOISE UNITS"] = "electron"
2356 metadata["LSST ISR PTCTURNOFF UNITS"] = "electron"
2357 metadata["LSST ISR GAIN SOURCE"] = self.config.useGainsFrom
2358 for amp in detector:
2359 # This includes any gain correction (if applied).
2360 metadata[f"LSST ISR GAIN {amp.getName()}"] = gains[amp.getName()]
2362 # At this stage, the read noise is always in electrons.
2363 metadata[f"LSST ISR READNOISE {amp.getName()}"] = ptc.noise[amp.getName()]
2365 # And the ptc turnoff will be in electrons.
2366 turnoff = ptc.ptcTurnoff[amp.getName()] * gains[amp.getName()]
2367 metadata[f"LSST ISR PTCTURNOFF {amp.getName()}"] = turnoff
2369 # Do crosstalk correction in the full region.
2370 # Output units: electron (adu if doBootstrap=True)
2371 if self.config.doCrosstalk:
2372 self.log.info("Applying crosstalk corrections to full amplifier region.")
2373 if self.config.doBootstrap and numpy.any(crosstalk.fitGains != 0):
2374 crosstalkGains = None
2375 else:
2376 crosstalkGains = gains
2377 self.crosstalk.run(
2378 ccdExposure,
2379 crosstalk=crosstalk,
2380 gains=crosstalkGains,
2381 fullAmplifier=True,
2382 badAmpDict=badAmpDict,
2383 ignoreVariance=True,
2384 )
2385 ccdExposure.metadata["LSST ISR CROSSTALK APPLIED"] = True
2387 # After crosstalk, we check for amplifier noise and state changes.
2388 if numpy.isfinite(self.config.serialOverscanMedianShiftSigmaThreshold):
2389 badAmpDict = self.checkAmpOverscanLevel(badAmpDict, ccdExposure, ptc)
2390 ccdExposure.metadata["LSST ISR OVERSCANLEVEL CHECKED"] = True
2392 if numpy.isfinite(self.config.ampNoiseThreshold):
2393 badAmpDict = self.checkAmpNoise(badAmpDict, ccdExposure, ptc)
2394 ccdExposure.metadata["LSST ISR NOISE CHECKED"] = True
2396 if numpy.isfinite(self.config.serialOverscanMedianShiftSigmaThreshold) or \
2397 numpy.isfinite(self.config.ampNoiseThreshold):
2398 self.checkAllBadAmps(badAmpDict, detector)
2400 # Parallel overscan correction.
2401 # Output units: electron (adu if doBootstrap=True)
2402 parallelOverscans = None
2403 if overscanDetectorConfig.doAnyParallelOverscan:
2404 # At the moment we do not use the return values from this task.
2405 parallelOverscans = self.overscanCorrection(
2406 "PARALLEL",
2407 overscanDetectorConfig,
2408 detector,
2409 badAmpDict,
2410 ccdExposure,
2411 )
2413 # Linearity correction
2414 # Output units: electron (adu if doBootstrap=True)
2415 if self.config.doLinearize:
2416 self.log.info("Applying linearizer.")
2417 # The linearizer is in units of adu.
2418 # If our units are electron, then pass in the gains
2419 # for conversion.
2420 if exposureMetadata["LSST ISR UNITS"] == "electron":
2421 linearityGains = gains
2422 else:
2423 linearityGains = None
2424 linearizer.applyLinearity(
2425 image=ccdExposure.image,
2426 detector=detector,
2427 log=self.log,
2428 gains=linearityGains,
2429 )
2430 ccdExposure.metadata["LSST ISR LINEARIZER APPLIED"] = True
2432 # Serial CTI (deferred charge) correction
2433 # This will be performed in electron units
2434 # Output units: same as input units
2435 if self.config.doDeferredCharge:
2436 if self.config.doBootstrap:
2437 self.log.info("Applying deferred charge correction with doBootstrap=True: "
2438 "will need to use deferredChargeCalib.inputGain to apply "
2439 "CTI correction in electron units.")
2440 deferredChargeGains = deferredChargeCalib.inputGain
2441 if numpy.all(numpy.isnan(list(deferredChargeGains.values()))):
2442 self.log.warning("All gains contained in the deferredChargeCalib are "
2443 "NaN, approximating with gain of 1.0.")
2444 deferredChargeGains = gains
2445 else:
2446 deferredChargeGains = gains
2447 self.deferredChargeCorrection.run(
2448 ccdExposure,
2449 deferredChargeCalib,
2450 gains=deferredChargeGains,
2451 )
2452 ccdExposure.metadata["LSST ISR CTI APPLIED"] = True
2454 # Save the untrimmed version for later statistics,
2455 # which still contains the overscan information
2456 untrimmedCcdExposure = None
2457 if self.config.isrStats.doCtiStatistics:
2458 untrimmedCcdExposure = ccdExposure.clone()
2460 # Assemble/trim
2461 # Output units: electron (adu if doBootstrap=True)
2462 if self.config.doAssembleCcd:
2463 self.log.info("Assembling CCD from amplifiers.")
2464 ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
2466 if self.config.expectWcs and not ccdExposure.getWcs():
2467 self.log.warning("No WCS found in input exposure.")
2469 # E2V edge bleed
2470 if self.config.doE2VEdgeBleedMask and detector.getPhysicalType() == "E2V":
2471 isrFunctions.maskE2VEdgeBleed(
2472 exposure=ccdExposure,
2473 e2vEdgeBleedSatMinArea=self.config.e2vEdgeBleedSatMinArea,
2474 e2vEdgeBleedSatMaxArea=self.config.e2vEdgeBleedSatMaxArea,
2475 e2vEdgeBleedYMax=self.config.e2vEdgeBleedYMax,
2476 saturatedMaskName=self.config.saturatedMaskName,
2477 log=self.log,
2478 )
2480 # ITL Dip Masking
2481 for maskPlane in self.config.itlDipMaskPlanes:
2482 if maskPlane not in ccdExposure.mask.getMaskPlaneDict():
2483 self.log.info("Adding %s mask plane to image.", maskPlane)
2484 ccdExposure.mask.addMaskPlane(maskPlane)
2486 if self.config.doITLDipMask:
2487 isrFunctions.maskITLDip(
2488 exposure=ccdExposure,
2489 detectorConfig=overscanDetectorConfig,
2490 log=self.log,
2491 maskPlaneNames=self.config.itlDipMaskPlanes,
2492 )
2494 if (self.config.doITLSatSagMask or self.config.doITLEdgeBleedMask) \
2495 and detector.getPhysicalType() == 'ITL':
2496 self.maskITLSatEdgesAndColumns(exposure=ccdExposure,
2497 badAmpDict=badAmpDict)
2499 # Bias subtraction
2500 # Output units: electron (adu if doBootstrap=True)
2501 if self.config.doBias:
2502 self.log.info("Applying bias correction.")
2503 # Bias frame and ISR unit consistency is checked at the top of
2504 # the run method.
2505 isrFunctions.biasCorrection(ccdExposure.maskedImage, bias.maskedImage)
2506 ccdExposure.metadata["LSST ISR BIAS APPLIED"] = True
2508 # Dark subtraction
2509 # Output units: electron (adu if doBootstrap=True)
2510 if self.config.doDark:
2511 self.log.info("Applying dark subtraction.")
2512 # Dark frame and ISR unit consistency is checked at the top of
2513 # the run method.
2514 self.darkCorrection(ccdExposure, dark)
2515 ccdExposure.metadata["LSST ISR DARK APPLIED"] = True
2517 # Defect masking
2518 # Masking block (defects, NAN pixels and trails).
2519 # Saturated and suspect pixels have already been masked.
2520 # Output units: electron (adu if doBootstrap=True)
2521 if self.config.doDefect:
2522 self.log.info("Applying defect masking.")
2523 self.maskDefects(ccdExposure, defects)
2524 ccdExposure.metadata["LSST ISR DEFECTS APPLIED"] = True
2526 self.log.info("Adding UNMASKEDNAN mask plane to image.")
2527 ccdExposure.mask.addMaskPlane("UNMASKEDNAN")
2528 if self.config.doNanMasking:
2529 self.log.info("Masking non-finite (NAN, inf) value pixels.")
2530 self.maskNan(ccdExposure)
2532 if self.config.doWidenSaturationTrails:
2533 self.log.info("Widening saturation trails.")
2534 isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
2536 # Brighter-Fatter
2537 # Output units: electron (adu if doBootstrap=True)
2538 if self.config.doBrighterFatter:
2539 if "ASTIER23" in self.config.brighterFatterCorrectionMethod:
2540 # Do the Astier et al. 2023 Brighter-Fatter correction
2541 self.log.info("Applying electrostatic brighter-fatter "
2542 "correction.")
2543 bfCalib = electroBfDistortionMatrix
2544 bfGains = electroBfDistortionMatrix.gain
2545 bfCorrFunction = self.applyElectrostaticBrighterFatterCorrection
2546 elif self.config.brighterFatterCorrectionMethod == "COULTON18":
2547 # Do the Coulton et al. 2018 Brighter-Fatter correction
2548 self.log.info("Applying brighter-fatter correction.")
2550 # Use the original gains used to compute the BFE
2551 bfCalib, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
2552 bfCorrFunction = self.applyBrighterFatterCorrection
2553 elif self.config.brighterFatterCorrectionMethod == "COULTON18_FLUX_CONSERVING":
2554 # Do the Coulton et al. 2018 correction with flux conserving
2555 # corrections
2556 bfCalib, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
2557 bfCorrFunction = self.applyFluxConservingBrighterFatterCorrection
2558 else:
2559 # This is an unknown BF correction type
2560 raise RuntimeError("Cannot perform brighter-fatter correction with unknown"
2561 "correction type.")
2563 # Needs to be done in electrons; applyBrighterFatterCorrection
2564 # will convert the image if necessary.
2565 if exposureMetadata["LSST ISR UNITS"] == "electron":
2566 brighterFatterApplyGain = False
2567 else:
2568 brighterFatterApplyGain = True
2570 if brighterFatterApplyGain and (ptc is not None) and (bfGains != gains):
2571 # The supplied ptc should be the same as the ptc used to
2572 # generate the bfKernel, in which case they will have the
2573 # same stored amp-keyed dictionary of gains. If not, there
2574 # is a mismatch in the calibrations being used. This should
2575 # not be always be a fatal error, but ideally, everything
2576 # should to be consistent.
2577 self.log.warning("Need to apply gain for brighter-fatter, but the stored"
2578 "gains in the kernel are not the same as the gains used"
2579 f"by {self.config.useGainsFrom}. Using the gains stored"
2580 "in the kernel.")
2582 # Apply the BF correction
2583 ccdExposure = bfCorrFunction(
2584 ccdExposure,
2585 flat,
2586 dark,
2587 bfCalib,
2588 brighterFatterApplyGain,
2589 bfGains,
2590 )
2592 ccdExposure.metadata["LSST ISR BF APPLIED"] = True
2593 ccdExposure.metadata["LSST ISR BF CORR METHOD"] = self.config.brighterFatterCorrectionMethod
2595 # Variance plane creation
2596 # Output units: electron (adu if doBootstrap=True)
2597 if self.config.doVariance:
2598 self.addVariancePlane(ccdExposure, detector)
2600 # Flat-fielding
2601 # This may move elsewhere, but this is the most convenient
2602 # location for simple flat-fielding for attractive backgrounds.
2603 # Output units: electron (adu if doBootstrap=True)
2604 if self.config.doFlat:
2605 self.log.info("Applying flat correction.")
2606 isrFunctions.flatCorrection(
2607 maskedImage=ccdExposure.maskedImage,
2608 flatMaskedImage=flat.maskedImage,
2609 scalingType=self.config.flatScalingType,
2610 userScale=self.config.flatUserScale,
2611 )
2613 # Copy over valid polygon from flat if it is
2614 # available, and set NO_DATA to 0.0 (which may
2615 # be inherited from the flat in the fully
2616 # vignetted region).
2617 if (validPolygon := flat.info.getValidPolygon()) is not None:
2618 ccdExposure.info.setValidPolygon(validPolygon)
2620 noData = (ccdExposure.mask.array & ccdExposure.mask.getPlaneBitMask("NO_DATA")) > 0
2621 ccdExposure.image.array[noData] = 0.0
2622 ccdExposure.variance.array[noData] = self.config.noDataVariance
2624 ccdExposure.metadata["LSST ISR FLAT APPLIED"] = True
2625 ccdExposure.metadata["LSST ISR FLAT SOURCE"] = flat.metadata.get("FLATSRC", "UNKNOWN")
2627 # Pixel values for masked regions are set here
2628 preInterpExp = None
2629 if self.config.doSaveInterpPixels:
2630 preInterpExp = ccdExposure.clone()
2632 if self.config.doSetBadRegions:
2633 self.log.info('Setting values in large contiguous bad regions.')
2634 self.setBadRegions(ccdExposure)
2636 if self.config.doInterpolate:
2637 self.log.info("Interpolating masked pixels.")
2638 isrFunctions.interpolateFromMask(
2639 maskedImage=ccdExposure.getMaskedImage(),
2640 fwhm=self.config.brighterFatterFwhmForInterpolation,
2641 growSaturatedFootprints=self.config.growSaturationFootprintSize,
2642 maskNameList=list(self.config.maskListToInterpolate),
2643 useLegacyInterp=self.config.useLegacyInterp,
2644 )
2646 # Calculate amp offset corrections within the CCD.
2647 if self.config.doAmpOffset:
2648 if self.config.ampOffset.doApplyAmpOffset:
2649 self.log.info("Measuring and applying amp offset corrections.")
2650 else:
2651 self.log.info("Measuring amp offset corrections only, without applying them.")
2652 self.ampOffset.run(ccdExposure)
2654 # Calculate standard image quality statistics
2655 if self.config.doStandardStatistics:
2656 metadata = ccdExposure.metadata
2657 for amp in detector:
2658 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
2659 ampName = amp.getName()
2660 metadata[f"LSST ISR MASK SAT {ampName}"] = isrFunctions.countMaskedPixels(
2661 ampExposure.getMaskedImage(),
2662 [self.config.saturatedMaskName]
2663 )
2664 metadata[f"LSST ISR MASK BAD {ampName}"] = isrFunctions.countMaskedPixels(
2665 ampExposure.getMaskedImage(),
2666 ["BAD"]
2667 )
2668 metadata[f"LSST ISR MASK SUSPECT {ampName}"] = isrFunctions.countMaskedPixels(
2669 ampExposure.getMaskedImage(),
2670 ["SUSPECT"],
2671 )
2672 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
2673 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP)
2675 metadata[f"LSST ISR FINAL MEAN {ampName}"] = qaStats.getValue(afwMath.MEAN)
2676 metadata[f"LSST ISR FINAL MEDIAN {ampName}"] = qaStats.getValue(afwMath.MEDIAN)
2677 metadata[f"LSST ISR FINAL STDEV {ampName}"] = qaStats.getValue(afwMath.STDEVCLIP)
2679 k1 = f"LSST ISR FINAL MEDIAN {ampName}"
2680 k2 = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
2681 if overscanDetectorConfig.doAnySerialOverscan and k1 in metadata and k2 in metadata:
2682 metadata[f"LSST ISR LEVEL {ampName}"] = metadata[k1] - metadata[k2]
2683 else:
2684 metadata[f"LSST ISR LEVEL {ampName}"] = numpy.nan
2686 # calculate additional statistics.
2687 outputStatistics = None
2688 if self.config.doCalculateStatistics:
2689 outputStatistics = self.isrStats.run(ccdExposure,
2690 untrimmedInputExposure=untrimmedCcdExposure,
2691 serialOverscanResults=serialOverscans,
2692 parallelOverscanResults=parallelOverscans,
2693 bias=bias, dark=dark, flat=flat,
2694 ptc=ptc, defects=defects).results
2696 # do image binning.
2697 outputBin1Exposure = None
2698 outputBin2Exposure = None
2699 if self.config.doBinnedExposures:
2700 self.log.info("Creating binned exposures.")
2701 outputBin1Exposure = self.binning.run(
2702 ccdExposure,
2703 binFactor=self.config.binFactor1,
2704 ).outputData
2705 outputBin2Exposure = self.binning.run(
2706 ccdExposure,
2707 binFactor=self.config.binFactor2,
2708 ).outputData
2710 return pipeBase.Struct(
2711 exposure=ccdExposure,
2713 outputBin1Exposure=outputBin1Exposure,
2714 outputBin2Exposure=outputBin2Exposure,
2716 preInterpExposure=preInterpExp,
2717 outputExposure=ccdExposure,
2718 outputStatistics=outputStatistics,
2719 )