Coverage for python / lsst / ip / isr / isrTaskLSST.py: 12%
904 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-14 23:58 +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 },
591 )
592 brighterFatterLevel = pexConfig.ChoiceField(
593 dtype=str,
594 doc="The level at which to correct for brighter-fatter.",
595 allowed={
596 "AMP": "Every amplifier treated separately.",
597 "DETECTOR": "One kernel per detector.",
598 },
599 default="DETECTOR",
600 )
601 brighterFatterMaxIter = pexConfig.Field(
602 dtype=int,
603 doc="Maximum number of iterations for the brighter-fatter correction for "
604 "the COULTON18* correction types",
605 default=10,
606 )
607 brighterFatterThreshold = pexConfig.Field(
608 dtype=float,
609 doc="Threshold used to stop iterating the brighter-fatter correction for "
610 "the COULTON18* correction types. It is the absolute value of the difference "
611 "between the current corrected image and the one from the previous iteration "
612 "summed over all the pixels.",
613 default=1000,
614 )
615 brighterFatterMaskListToInterpolate = pexConfig.ListField(
616 dtype=str,
617 doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
618 "correction.",
619 default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
620 )
621 brighterFatterMaskGrowSize = pexConfig.Field(
622 dtype=int,
623 doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
624 "when brighter-fatter correction is applied.",
625 default=2,
626 )
627 brighterFatterFwhmForInterpolation = pexConfig.Field(
628 dtype=float,
629 doc="FWHM of PSF in arcseconds used for interpolation in brighter-fatter correction "
630 "(currently unused).",
631 default=1.0,
632 )
633 growSaturationFootprintSize = pexConfig.Field(
634 dtype=int,
635 doc="Number of pixels by which to grow the saturation footprints.",
636 default=1,
637 )
639 # Dark subtraction.
640 doDark = pexConfig.Field(
641 dtype=bool,
642 doc="Apply dark frame correction.",
643 default=True,
644 )
646 # Flat correction.
647 doFlat = pexConfig.Field(
648 dtype=bool,
649 doc="Apply flat field correction.",
650 default=True,
651 )
652 flatScalingType = pexConfig.ChoiceField(
653 dtype=str,
654 doc="The method for scaling the flat on the fly.",
655 default='USER',
656 allowed={
657 "USER": "Scale by flatUserScale",
658 "MEAN": "Scale by the inverse of the mean",
659 "MEDIAN": "Scale by the inverse of the median",
660 },
661 )
662 flatUserScale = pexConfig.Field(
663 dtype=float,
664 doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise.",
665 default=1.0,
666 )
668 # Calculate image quality statistics?
669 doStandardStatistics = pexConfig.Field(
670 dtype=bool,
671 doc="Should standard image quality statistics be calculated?",
672 default=True,
673 )
674 # Calculate additional statistics?
675 doCalculateStatistics = pexConfig.Field(
676 dtype=bool,
677 doc="Should additional ISR statistics be calculated?",
678 default=True,
679 )
680 isrStats = pexConfig.ConfigurableField(
681 target=IsrStatisticsTask,
682 doc="Task to calculate additional statistics.",
683 )
685 # Make binned images?
686 doBinnedExposures = pexConfig.Field(
687 dtype=bool,
688 doc="Should binned exposures be calculated?",
689 default=False,
690 )
691 binning = pexConfig.ConfigurableField(
692 target=BinImageDataTask,
693 doc="Task to bin the exposure.",
694 )
695 binFactor1 = pexConfig.Field(
696 dtype=int,
697 doc="Binning factor for first binned exposure. This is intended for a finely binned output.",
698 default=8,
699 check=lambda x: x > 1,
700 )
701 binFactor2 = pexConfig.Field(
702 dtype=int,
703 doc="Binning factor for second binned exposure. This is intended for a coarsely binned output.",
704 default=64,
705 check=lambda x: x > 1,
706 )
708 def validate(self):
709 super().validate()
711 if self.doBootstrap:
712 # Additional checks in bootstrap (no gains) mode.
713 if self.doApplyGains:
714 raise ValueError("Cannot run task with doBootstrap=True and doApplyGains=True.")
715 if self.doCorrectGains:
716 raise ValueError("Cannot run task with doBootstrap=True and doCorrectGains=True.")
717 if self.doCrosstalk and self.crosstalk.doQuadraticCrosstalkCorrection:
718 raise ValueError("Cannot apply quadratic crosstalk correction with doBootstrap=True.")
719 if numpy.isfinite(self.serialOverscanMedianShiftSigmaThreshold):
720 raise ValueError("Cannot do amp overscan level checks with doBootstrap=True.")
721 if numpy.isfinite(self.ampNoiseThreshold):
722 raise ValueError("Cannot do amp noise thresholds with doBootstrap=True.")
724 if self.doITLEdgeBleedMask and not self.doSaturation:
725 raise ValueError("Cannot do ITL edge bleed masking when doSaturation=False.")
726 if self.doE2VEdgeBleedMask and not self.doSaturation:
727 raise ValueError("Cannot do e2v edge bleed masking when doSaturation=False.")
729 def setDefaults(self):
730 super().setDefaults()
733class IsrTaskLSST(pipeBase.PipelineTask):
734 ConfigClass = IsrTaskLSSTConfig
735 _DefaultName = "isrLSST"
737 def __init__(self, **kwargs):
738 super().__init__(**kwargs)
739 self.makeSubtask("assembleCcd")
740 self.makeSubtask("deferredChargeCorrection")
741 self.makeSubtask("crosstalk")
742 self.makeSubtask("masking")
743 self.makeSubtask("isrStats")
744 self.makeSubtask("ampOffset")
745 self.makeSubtask("binning")
747 def runQuantum(self, butlerQC, inputRefs, outputRefs):
749 inputs = butlerQC.get(inputRefs)
750 self.validateInput(inputs)
752 if self.config.doHeaderProvenance:
753 # Add calibration provenanace info to header.
754 exposureMetadata = inputs['ccdExposure'].metadata
755 for inputName in sorted(list(inputs.keys())):
756 reference = getattr(inputRefs, inputName, None)
757 if reference is not None and hasattr(reference, "run"):
758 runKey = f"LSST CALIB RUN {inputName.upper()}"
759 runValue = reference.run
760 idKey = f"LSST CALIB UUID {inputName.upper()}"
761 idValue = str(reference.id)
762 dateKey = f"LSST CALIB DATE {inputName.upper()}"
763 dateValue = self.extractCalibDate(inputs[inputName])
764 if dateValue != "Unknown Unknown":
765 butlerQC.add_additional_provenance(reference, {"calib date": dateValue})
767 exposureMetadata[runKey] = runValue
768 exposureMetadata[idKey] = idValue
769 exposureMetadata[dateKey] = dateValue
771 outputs = self.run(**inputs)
772 butlerQC.put(outputs, outputRefs)
774 def validateInput(self, inputs):
775 """
776 This is a check that all the inputs required by the config
777 are available.
778 """
780 # Inputs depend on where the gains are coming from
781 doApplyGains = self.config.doApplyGains
782 useLinearizerGains = self.config.useGainsFrom == "LINEARIZER"
783 usePtcGains = not useLinearizerGains
785 # Inputs depend on the BF method
786 method = self.config.brighterFatterCorrectionMethod
788 inputMap = {
789 'dnlLUT': self.config.doDiffNonLinearCorrection,
790 'bias': self.config.doBias,
791 'deferredChargeCalib': self.config.doDeferredCharge,
792 # Some tasks require gains in order to be
793 # supplied regardless of whether
794 # self.config.doApplyGains is True or False.
795 'linearizer': (self.config.doLinearize or (doApplyGains and useLinearizerGains)),
796 'ptc': self.config.doApplyGains and usePtcGains,
797 'crosstalk': self.config.doCrosstalk,
798 'defects': self.config.doDefect,
799 'bfKernel': self.config.doBrighterFatter & method.startswith("COULTON18"),
800 'electroBfDistortionMatrix': self.config.doBrighterFatter & (method == "ASTIER23"),
801 'dark': self.config.doDark,
802 }
804 for calibrationFile, configValue in inputMap.items():
805 if configValue and inputs[calibrationFile] is None:
806 raise RuntimeError("Must supply ", calibrationFile)
808 def diffNonLinearCorrection(self, ccdExposure, dnlLUT, **kwargs):
809 # TODO DM 36636
810 # isrFunctions.diffNonLinearCorrection
811 pass
813 def maskFullAmplifiers(self, ccdExposure, detector, defects, gains=None):
814 """
815 Check for fully masked bad amplifiers and mask them.
817 This includes defects which cover full amplifiers, as well
818 as amplifiers with nan gain values which should be used
819 if self.config.doApplyGains=True.
821 Full defect masking happens later to allow for defects which
822 cross amplifier boundaries.
824 Parameters
825 ----------
826 ccdExposure : `lsst.afw.image.Exposure`
827 Input exposure to be masked.
828 detector : `lsst.afw.cameraGeom.Detector`
829 Detector object.
830 defects : `lsst.ip.isr.Defects`
831 List of defects. Used to determine if an entire
832 amplifier is bad.
833 gains : `dict` [`str`, `float`], optional
834 Dictionary of gains to check if
835 self.config.doApplyGains=True.
837 Returns
838 -------
839 badAmpDict : `str`[`bool`]
840 Dictionary of amplifiers, keyed by name, value is True if
841 amplifier is fully masked.
842 """
843 badAmpDict = {}
845 maskedImage = ccdExposure.getMaskedImage()
847 for amp in detector:
848 ampName = amp.getName()
849 badAmpDict[ampName] = False
851 # Check if entire amp region is defined as a defect
852 # NB: need to use amp.getBBox() for correct comparison with current
853 # defects definition.
854 if defects is not None:
855 badAmpDict[ampName] = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
857 if badAmpDict[ampName]:
858 self.log.warning("Amplifier %s is bad (completely covered with defects)", ampName)
860 if gains is not None and self.config.doApplyGains:
861 if not math.isfinite(gains[ampName]):
862 badAmpDict[ampName] = True
864 self.log.warning("Amplifier %s is bad (non-finite gain)", ampName)
866 # In the case of a bad amp, we will set mask to "BAD"
867 # (here use amp.getRawBBox() for correct association with pixels in
868 # current ccdExposure).
869 if badAmpDict[ampName]:
870 dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
871 afwImage.PARENT)
872 maskView = dataView.getMask()
873 maskView |= maskView.getPlaneBitMask("BAD")
874 del maskView
876 return badAmpDict
878 def checkAllBadAmps(self, badAmpDict, detector):
879 """Check if all amps are marked as bad.
881 Parameters
882 ----------
883 badAmpDict : `str`[`bool`]
884 Dictionary of amplifiers, keyed by name, value is True if
885 amplifier is fully masked.
886 detector : `lsst.afw.cameraGeom.Detector`
887 Detector object.
889 Raises
890 ------
891 UnprocessableDataError if all amps are bad and doCheckUnprocessableData
892 configuration is True.
893 """
894 if not self.config.doCheckUnprocessableData:
895 return
897 for amp in detector:
898 if not badAmpDict.get(amp.getName(), False):
899 return
901 raise UnprocessableDataError(f"All amps in the exposure {detector.getName()} are bad; skipping ISR.")
903 def checkAmpOverscanLevel(self, badAmpDict, exposure, ptc):
904 """Check if the amplifier overscan levels have changed.
906 Any amplifier that has an overscan median level that has changed
907 significantly will be masked as BAD and added to toe badAmpDict.
909 Parameters
910 ----------
911 badAmpDict : `str` [`bool`]
912 Dictionary of amplifiers, keyed by name, value is True if
913 amplifier is fully masked.
914 exposure : `lsst.afw.image.Exposure`
915 Input exposure to be masked (untrimmed).
916 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
917 PTC dataset with gains/read noises.
919 Returns
920 -------
921 badAmpDict : `str`[`bool`]
922 Dictionary of amplifiers, keyed by name.
924 """
925 if isTrimmedExposure(exposure):
926 raise RuntimeError("checkAmpOverscanLevel must be run on an untrimmed exposure.")
928 # We want to consolidate all the amps into one warning if necessary.
929 # This config should not be set to a finite threshold if the necessary
930 # data is not in the PTC.
931 missingWarnString = "No PTC overscan information for amplifier "
932 missingWarnFlag = False
933 for amp in exposure.getDetector():
934 ampName = amp.getName()
936 if not numpy.isfinite(ptc.overscanMedian[ampName]) or \
937 not numpy.isfinite(ptc.overscanMedianSigma[ampName]):
938 missingWarnString += f"{ampName},"
939 missingWarnFlag = True
940 else:
941 key = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
942 # If it is missing, just return the PTC value and it
943 # will be skipped.
944 overscanLevel = exposure.metadata.get(key, ptc.overscanMedian[ampName])
945 pull = (overscanLevel - ptc.overscanMedian[ampName])/ptc.overscanMedianSigma[ampName]
946 if numpy.abs(pull) > self.config.serialOverscanMedianShiftSigmaThreshold:
947 self.log.warning(
948 "Amplifier %s has an overscan level that is %.2f sigma from the expected level; "
949 "masking it as BAD.",
950 ampName,
951 pull,
952 )
954 badAmpDict[ampName] = True
955 exposure.mask[amp.getRawBBox()] |= exposure.mask.getPlaneBitMask("BAD")
957 if missingWarnFlag:
958 self.log.warning(missingWarnString)
960 return badAmpDict
962 def checkAmpNoise(self, badAmpDict, exposure, ptc):
963 """Check if amplifier noise levels are above threshold.
965 Any amplifier that is above the noise level will be masked as BAD
966 and added to the badAmpDict.
968 Parameters
969 ----------
970 badAmpDict : `str` [`bool`]
971 Dictionary of amplifiers, keyed by name, value is True if
972 amplifier is fully masked.
973 exposure : `lsst.afw.image.Exposure`
974 Input exposure to be masked (untrimmed).
975 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`
976 PTC dataset with gains/read noises.
978 Returns
979 -------
980 badAmpDict : `str`[`bool`]
981 Dictionary of amplifiers, keyed by name.
982 """
984 if isTrimmedExposure(exposure):
985 raise RuntimeError("checkAmpNoise must be run on an untrimmed exposure.")
987 for amp in exposure.getDetector():
988 ampName = amp.getName()
990 doMask = False
991 if ptc.noise[ampName] > self.config.ampNoiseThreshold:
992 self.log.info(
993 "Amplifier %s has a PTC noise level of %.2f e-, above threshold.",
994 ampName,
995 ptc.noise[ampName],
996 )
997 doMask = True
998 else:
999 key = f"LSST ISR OVERSCAN RESIDUAL SERIAL STDEV {ampName}"
1000 overscanNoise = exposure.metadata.get(key, numpy.nan)
1001 if overscanNoise * ptc.gain[ampName] > self.config.ampNoiseThreshold:
1002 self.log.warning(
1003 "Amplifier %s has an overscan read noise level of %.2f e-, above threshold.",
1004 ampName,
1005 overscanNoise * ptc.gain[ampName],
1006 )
1007 doMask = True
1009 if doMask:
1010 badAmpDict[ampName] = True
1012 exposure.mask[amp.getRawBBox()] |= exposure.mask.getPlaneBitMask("BAD")
1014 return badAmpDict
1016 def maskSaturatedPixels(self, badAmpDict, ccdExposure, detector, detectorConfig, ptc=None):
1017 """
1018 Mask SATURATED and SUSPECT pixels and check if any amplifiers
1019 are fully masked.
1021 Parameters
1022 ----------
1023 badAmpDict : `str` [`bool`]
1024 Dictionary of amplifiers, keyed by name, value is True if
1025 amplifier is fully masked.
1026 ccdExposure : `lsst.afw.image.Exposure`
1027 Input exposure to be masked.
1028 detector : `lsst.afw.cameraGeom.Detector`
1029 Detector object.
1030 defects : `lsst.ip.isr.Defects`
1031 List of defects. Used to determine if an entire
1032 amplifier is bad.
1033 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1034 Per-amplifier configurations.
1035 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
1036 PTC dataset (used if configured to use PTCTURNOFF).
1038 Returns
1039 -------
1040 badAmpDict : `str`[`bool`]
1041 Dictionary of amplifiers, keyed by name.
1042 """
1043 maskedImage = ccdExposure.getMaskedImage()
1045 metadata = ccdExposure.metadata
1047 if self.config.doSaturation and self.config.defaultSaturationSource == "PTCTURNOFF" and ptc is None:
1048 raise RuntimeError("Must provide ptc if using PTCTURNOFF as saturation source.")
1049 if self.config.doSuspect and self.config.defaultSuspectSource == "PTCTURNOFF" and ptc is None:
1050 raise RuntimeError("Must provide ptc if using PTCTURNOFF as suspect source.")
1052 for amp in detector:
1053 ampName = amp.getName()
1055 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
1057 if badAmpDict[ampName]:
1058 # No need to check fully bad amplifiers.
1059 continue
1061 # Mask saturated and suspect pixels.
1062 limits = {}
1063 if self.config.doSaturation:
1064 if self.config.defaultSaturationSource == "PTCTURNOFF":
1065 limits.update({self.config.saturatedMaskName: ptc.ptcTurnoff[amp.getName()]})
1066 elif self.config.defaultSaturationSource == "CAMERAMODEL":
1067 # Set to the default from the camera model.
1068 limits.update({self.config.saturatedMaskName: amp.getSaturation()})
1069 elif self.config.defaultSaturationSource == "NONE":
1070 limits.update({self.config.saturatedMaskName: numpy.inf})
1072 # And update if it is set in the config.
1073 if math.isfinite(ampConfig.saturation):
1074 limits.update({self.config.saturatedMaskName: ampConfig.saturation})
1075 metadata[f"LSST ISR SATURATION LEVEL {ampName}"] = limits[self.config.saturatedMaskName]
1077 if self.config.doSuspect:
1078 if self.config.defaultSuspectSource == "PTCTURNOFF":
1079 limits.update({self.config.suspectMaskName: ptc.ptcTurnoff[amp.getName()]})
1080 elif self.config.defaultSuspectSource == "CAMERAMODEL":
1081 # Set to the default from the camera model.
1082 limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
1083 elif self.config.defaultSuspectSource == "NONE":
1084 limits.update({self.config.suspectMaskName: numpy.inf})
1086 # And update if it set in the config.
1087 if math.isfinite(ampConfig.suspectLevel):
1088 limits.update({self.config.suspectMaskName: ampConfig.suspectLevel})
1089 metadata[f"LSST ISR SUSPECT LEVEL {ampName}"] = limits[self.config.suspectMaskName]
1091 for maskName, maskThreshold in limits.items():
1092 if not math.isnan(maskThreshold):
1093 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
1094 toMask = (dataView.image.array >= maskThreshold)
1095 dataView.mask.array[toMask] |= dataView.mask.getPlaneBitMask(maskName)
1097 # Determine if we've fully masked this amplifier with SUSPECT and
1098 # SAT pixels.
1099 maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
1100 afwImage.PARENT)
1101 maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
1102 self.config.suspectMaskName])
1103 if numpy.all(maskView.getArray() & maskVal > 0):
1104 self.log.warning("Amplifier %s is bad (completely SATURATED or SUSPECT)", ampName)
1105 badAmpDict[ampName] = True
1106 maskView |= maskView.getPlaneBitMask("BAD")
1108 return badAmpDict
1110 def maskITLSatEdgesAndColumns(self, exposure, badAmpDict):
1112 # The following steps will rely on the footprint of saturated
1113 # cores with large areas.
1114 thresh = afwDetection.Threshold(exposure.mask.getPlaneBitMask("SAT"),
1115 afwDetection.Threshold.BITMASK
1116 )
1117 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
1119 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
1120 largeAreas, = numpy.where((satAreas >= self.config.itlEdgeBleedSatMinArea)
1121 & (satAreas < self.config.itlEdgeBleedSatMaxArea))
1123 for largeAreasIndex in largeAreas:
1125 fpCore = fpList[largeAreasIndex]
1127 # Edge bleed masking
1128 if self.config.doITLEdgeBleedMask:
1129 isrFunctions.maskITLEdgeBleed(ccdExposure=exposure,
1130 badAmpDict=badAmpDict,
1131 fpCore=fpCore,
1132 itlEdgeBleedSatMinArea=self.config.itlEdgeBleedSatMinArea,
1133 itlEdgeBleedSatMaxArea=self.config.itlEdgeBleedSatMaxArea,
1134 itlEdgeBleedThreshold=self.config.itlEdgeBleedThreshold,
1135 itlEdgeBleedModelConstant=self.config.itlEdgeBleedModelConstant,
1136 saturatedMaskName=self.config.saturatedMaskName,
1137 log=self.log
1138 )
1139 if self.config.doITLSatSagMask:
1140 isrFunctions.maskITLSatSag(ccdExposure=exposure, fpCore=fpCore,
1141 saturatedMaskName=self.config.saturatedMaskName)
1143 def overscanCorrection(self, mode, detectorConfig, detector, badAmpDict, ccdExposure):
1144 """Apply serial overscan correction in place to all amps.
1146 The actual overscan subtraction is performed by the
1147 `lsst.ip.isr.overscan.OverscanTask`, which is called here.
1149 Parameters
1150 ----------
1151 mode : `str`
1152 Must be `SERIAL` or `PARALLEL`.
1153 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1154 Per-amplifier configurations.
1155 detector : `lsst.afw.cameraGeom.Detector`
1156 Detector object.
1157 badAmpDict : `dict`
1158 Dictionary of amp name to whether it is a bad amp.
1159 ccdExposure : `lsst.afw.image.Exposure`
1160 Exposure to have overscan correction performed.
1162 Returns
1163 -------
1164 overscans : `list` [`lsst.pipe.base.Struct` or None]
1165 Overscan measurements (always in adu).
1166 Each result struct has components:
1168 ``imageFit``
1169 Value or fit subtracted from the amplifier image data.
1170 (scalar or `lsst.afw.image.Image`)
1171 ``overscanFit``
1172 Value or fit subtracted from the overscan image data.
1173 (scalar or `lsst.afw.image.Image`)
1174 ``overscanImage``
1175 Image of the overscan region with the overscan
1176 correction applied. This quantity is used to estimate
1177 the amplifier read noise empirically.
1178 (`lsst.afw.image.Image`)
1179 ``overscanMean``
1180 Mean overscan fit value. (`float`)
1181 ``overscanMedian``
1182 Median overscan fit value. (`float`)
1183 ``overscanSigma``
1184 Clipped standard deviation of the overscan fit. (`float`)
1185 ``residualMean``
1186 Mean of the overscan after fit subtraction. (`float`)
1187 ``residualMedian``
1188 Median of the overscan after fit subtraction. (`float`)
1189 ``residualSigma``
1190 Clipped standard deviation of the overscan after fit
1191 subtraction. (`float`)
1193 See Also
1194 --------
1195 lsst.ip.isr.overscan.OverscanTask
1196 """
1197 if mode not in ["SERIAL", "PARALLEL"]:
1198 raise ValueError("Mode must be SERIAL or PARALLEL")
1200 # This returns a list in amp order, with None for uncorrected amps.
1201 overscans = []
1203 for i, amp in enumerate(detector):
1204 ampName = amp.getName()
1206 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
1208 if mode == "SERIAL" and not ampConfig.doSerialOverscan:
1209 self.log.debug(
1210 "ISR_OSCAN: Amplifier %s/%s configured to skip serial overscan.",
1211 detector.getName(),
1212 ampName,
1213 )
1214 results = None
1215 elif mode == "PARALLEL" and not ampConfig.doParallelOverscan:
1216 self.log.debug(
1217 "ISR_OSCAN: Amplifier %s configured to skip parallel overscan.",
1218 detector.getName(),
1219 ampName,
1220 )
1221 results = None
1222 elif badAmpDict[ampName] or not ccdExposure.getBBox().contains(amp.getBBox()):
1223 results = None
1224 else:
1225 # This check is to confirm that we are not trying to run
1226 # overscan on an already trimmed image.
1227 if isTrimmedExposure(ccdExposure):
1228 self.log.warning(
1229 "ISR_OSCAN: No overscan region for amp %s. Not performing overscan correction.",
1230 ampName,
1231 )
1232 results = None
1233 else:
1234 if mode == "SERIAL":
1235 # We need to set up the subtask here with a custom
1236 # configuration.
1237 serialOverscan = SerialOverscanCorrectionTask(config=ampConfig.serialOverscanConfig)
1238 results = serialOverscan.run(ccdExposure, amp)
1239 else:
1240 config = ampConfig.parallelOverscanConfig
1241 parallelOverscan = ParallelOverscanCorrectionTask(
1242 config=config,
1243 )
1245 metadata = ccdExposure.metadata
1247 # We need to know the saturation level that was used
1248 # for the parallel overscan masking. If it isn't set
1249 # then the configured parallelOverscanSaturationLevel
1250 # will be used instead (assuming
1251 # doParallelOverscanSaturation is True). Note that
1252 # this will have the correct units (adu or electron)
1253 # depending on whether the gain has been applied.
1254 if self.config.doSaturation:
1255 saturationLevel = metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
1256 saturationLevel *= config.parallelOverscanSaturationLevelAdjustmentFactor
1257 else:
1258 saturationLevel = config.parallelOverscanSaturationLevel
1259 if ccdExposure.metadata["LSST ISR UNITS"] == "electron":
1260 # Need to convert to electron from adu.
1261 saturationLevel *= metadata[f"LSST ISR GAIN {amp.getName()}"]
1263 self.log.debug(
1264 "Using saturation level of %.2f for parallel overscan amp %s",
1265 saturationLevel,
1266 amp.getName(),
1267 )
1269 parallelOverscan.maskParallelOverscanAmp(
1270 ccdExposure,
1271 amp,
1272 saturationLevel=saturationLevel,
1273 )
1275 results = parallelOverscan.run(ccdExposure, amp)
1277 metadata = ccdExposure.metadata
1278 keyBase = "LSST ISR OVERSCAN"
1280 # The overscan is always in adu for the serial mode,
1281 # but, it may be electron in the parallel mode if
1282 # doApplyGains==True. If doApplyGains==True, then the
1283 # gains are applied to the untrimmed image, so the
1284 # overscan statistics units here will always match the
1285 # units of the image at this point.
1286 metadata[f"{keyBase} {mode} UNITS"] = ccdExposure.metadata["LSST ISR UNITS"]
1287 metadata[f"{keyBase} {mode} MEAN {ampName}"] = results.overscanMean
1288 metadata[f"{keyBase} {mode} MEDIAN {ampName}"] = results.overscanMedian
1289 metadata[f"{keyBase} {mode} STDEV {ampName}"] = results.overscanSigma
1291 metadata[f"{keyBase} RESIDUAL {mode} MEAN {ampName}"] = results.residualMean
1292 metadata[f"{keyBase} RESIDUAL {mode} MEDIAN {ampName}"] = results.residualMedian
1293 metadata[f"{keyBase} RESIDUAL {mode} STDEV {ampName}"] = results.residualSigma
1295 overscans.append(results)
1297 # Question: should this be finer grained?
1298 ccdExposure.metadata.set("OVERSCAN", "Overscan corrected")
1300 return overscans
1302 def maskNegativeVariance(self, exposure):
1303 """Identify and mask pixels with negative variance values.
1305 Parameters
1306 ----------
1307 exposure : `lsst.afw.image.Exposure`
1308 Exposure to process.
1310 See Also
1311 --------
1312 lsst.ip.isr.isrFunctions.updateVariance
1313 """
1314 maskPlane = exposure.getMask().getPlaneBitMask(self.config.negativeVarianceMaskName)
1315 bad = numpy.where(exposure.getVariance().getArray() <= 0.0)
1316 exposure.mask.array[bad] |= maskPlane
1318 def addVariancePlane(self, exposure, detector):
1319 """Add the variance plane to the image.
1321 The gain and read noise per amp must have been set in the
1322 exposure metadata as ``LSST ISR GAIN ampName`` and
1323 ``LSST ISR READNOISE ampName`` with the units of the image.
1324 Unit conversions for the variance plane will be done as
1325 necessary based on the exposure units.
1327 The units of the variance plane will always be of the same
1328 type as the units of the input image itself
1329 (``LSST ISR UNITS``^2).
1331 Parameters
1332 ----------
1333 exposure : `lsst.afw.image.Exposure`
1334 The exposure to add the variance plane.
1335 detector : `lsst.afw.cameraGeom.Detector`
1336 Detector with geometry info.
1337 """
1338 # NOTE: this will fail if the exposure is not trimmed.
1339 if not isTrimmedExposure(exposure):
1340 raise RuntimeError("Exposure must be trimmed to add variance plane.")
1342 isElectrons = (exposure.metadata["LSST ISR UNITS"] == "electron")
1344 for amp in detector:
1345 if exposure.getBBox().contains(amp.getBBox()):
1346 self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1347 ampExposure = exposure.Factory(exposure, amp.getBBox())
1349 # The effective gain is 1.0 if we are in electron units.
1350 # The metadata read noise is in the same units as the image.
1351 gain = exposure.metadata[f"LSST ISR GAIN {amp.getName()}"] if not isElectrons else 1.0
1352 readNoise = exposure.metadata[f"LSST ISR READNOISE {amp.getName()}"]
1354 isrFunctions.updateVariance(
1355 maskedImage=ampExposure.maskedImage,
1356 gain=gain,
1357 readNoise=readNoise,
1358 replace=False,
1359 )
1361 if self.config.qa is not None and self.config.qa.saveStats is True:
1362 qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1363 afwMath.MEDIAN | afwMath.STDEVCLIP)
1364 self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1365 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1366 qaStats.getValue(afwMath.STDEVCLIP))
1368 if self.config.maskNegativeVariance:
1369 self.maskNegativeVariance(exposure)
1371 def maskDefects(self, exposure, defectBaseList):
1372 """Mask defects using mask plane "BAD", in place.
1374 Parameters
1375 ----------
1376 exposure : `lsst.afw.image.Exposure`
1377 Exposure to process.
1379 defectBaseList : defect-type
1380 List of defects to mask. Can be of type `lsst.ip.isr.Defects`
1381 or `list` of `lsst.afw.image.DefectBase`.
1382 """
1383 maskedImage = exposure.getMaskedImage()
1384 if not isinstance(defectBaseList, Defects):
1385 # Promotes DefectBase to Defect
1386 defectList = Defects(defectBaseList)
1387 else:
1388 defectList = defectBaseList
1389 defectList.maskPixels(maskedImage, maskName="BAD")
1391 if len(self.config.badAmps) == 0:
1392 return
1394 detector = exposure.getDetector()
1395 mask = maskedImage.mask
1396 for badAmp in self.config.badAmps:
1397 if badAmp.startswith(detector.getName()):
1398 # Split on the full detector name plus _, which
1399 # gives us an empty string and the amp name.
1400 ampName = badAmp.split(detector.getName() + "_")[-1]
1401 self.log.info("Masking amplifier %s as bad via config.", ampName)
1402 mask[detector[ampName].getBBox()].array[:, :] |= mask.getPlaneBitMask("BAD")
1404 def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
1405 """Mask edge pixels with applicable mask plane.
1407 Parameters
1408 ----------
1409 exposure : `lsst.afw.image.Exposure`
1410 Exposure to process.
1411 numEdgePixels : `int`, optional
1412 Number of edge pixels to mask.
1413 maskPlane : `str`, optional
1414 Mask plane name to use.
1415 level : `str`, optional
1416 Level at which to mask edges.
1417 """
1418 maskedImage = exposure.getMaskedImage()
1419 maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
1421 if numEdgePixels > 0:
1422 if level == 'DETECTOR':
1423 boxes = [maskedImage.getBBox()]
1424 elif level == 'AMP':
1425 boxes = [amp.getBBox() for amp in exposure.getDetector()]
1427 for box in boxes:
1428 # This makes a bbox numEdgeSuspect pixels smaller than the
1429 # image on each side
1430 subImage = maskedImage[box]
1431 box.grow(-numEdgePixels)
1432 # Mask pixels outside box
1433 SourceDetectionTask.setEdgeBits(
1434 subImage,
1435 box,
1436 maskBitMask)
1438 def maskNan(self, exposure):
1439 """Mask NaNs using mask plane "UNMASKEDNAN", in place.
1441 Parameters
1442 ----------
1443 exposure : `lsst.afw.image.Exposure`
1444 Exposure to process.
1446 Notes
1447 -----
1448 We mask over all non-finite values (NaN, inf), including those
1449 that are masked with other bits (because those may or may not be
1450 interpolated over later, and we want to remove all NaN/infs).
1451 Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
1452 preserve the historical name.
1453 """
1454 maskedImage = exposure.getMaskedImage()
1456 # Find and mask NaNs
1457 maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
1458 numNans = maskNans(maskedImage, maskVal)
1459 self.metadata["NUMNANS"] = numNans
1460 if numNans > 0:
1461 self.log.warning("There were %d unmasked NaNs.", numNans)
1463 def setBadRegions(self, exposure):
1464 """Set bad regions from large contiguous regions.
1466 Parameters
1467 ----------
1468 exposure : `lsst.afw.Exposure`
1469 Exposure to set bad regions.
1471 Notes
1472 -----
1473 Reset and interpolate bad pixels.
1475 Large contiguous bad regions (which should have the BAD mask
1476 bit set) should have their values set to the image median.
1477 This group should include defects and bad amplifiers. As the
1478 area covered by these defects are large, there's little
1479 reason to expect that interpolation would provide a more
1480 useful value.
1482 Smaller defects can be safely interpolated after the larger
1483 regions have had their pixel values reset. This ensures
1484 that the remaining defects adjacent to bad amplifiers (as an
1485 example) do not attempt to interpolate extreme values.
1486 """
1487 badPixelCount, badPixelValue = isrFunctions.setBadRegions(exposure)
1488 if badPixelCount > 0:
1489 self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1491 @contextmanager
1492 def flatContext(self, exp, flat, dark=None):
1493 """Context manager that applies and removes flats and darks,
1494 if the task is configured to apply them.
1496 Parameters
1497 ----------
1498 exp : `lsst.afw.image.Exposure`
1499 Exposure to process.
1500 flat : `lsst.afw.image.Exposure`
1501 Flat exposure the same size as ``exp``.
1502 dark : `lsst.afw.image.Exposure`, optional
1503 Dark exposure the same size as ``exp``.
1505 Yields
1506 ------
1507 exp : `lsst.afw.image.Exposure`
1508 The flat and dark corrected exposure.
1509 """
1510 if self.config.doDark and dark is not None:
1511 self.darkCorrection(exp, dark)
1512 if self.config.doFlat and flat is not None:
1513 self.flatCorrection(exp, flat)
1514 try:
1515 yield exp
1516 finally:
1517 if self.config.doFlat and flat is not None:
1518 self.flatCorrection(exp, flat, invert=True)
1519 if self.config.doDark and dark is not None:
1520 self.darkCorrection(exp, dark, invert=True)
1522 def getBrighterFatterKernel(self, detector, bfKernel):
1523 detName = detector.getName()
1525 # This is expected to be a dictionary of amp-wise gains.
1526 bfGains = bfKernel.gain
1527 if bfKernel.level == 'DETECTOR':
1528 if detName in bfKernel.detKernels:
1529 bfKernelOut = bfKernel.detKernels[detName]
1530 return bfKernelOut, bfGains
1531 else:
1532 raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1533 elif bfKernel.level == 'AMP':
1534 self.log.info("Making DETECTOR level kernel from AMP based brighter "
1535 "fatter kernels.")
1536 bfKernel.makeDetectorKernelFromAmpwiseKernels(detName)
1537 bfKernelOut = bfKernel.detKernels[detName]
1538 return bfKernelOut, bfGains
1540 def applyElectrostaticBrighterFatterCorrection(self, ccdExposure, flat, dark, electroBfDistortionMatrix,
1541 brighterFatterApplyGain, bfGains):
1542 """Apply an electrostatic brighter fatter correction to the image
1543 using the method defined in Astier et al. 2023.
1545 Note that this correction requires that the image is in units
1546 electrons.
1548 Parameters
1549 ----------
1550 ccdExposure : `lsst.afw.image.Exposure`
1551 Exposure to process.
1552 flat : `lsst.afw.image.Exposure`
1553 Flat exposure the same size as ``exp``.
1554 dark : `lsst.afw.image.Exposure`, optional
1555 Dark exposure the same size as ``exp``.
1556 electroBfDistortionMatrix : `lsst.ip.isr.ElectrostaticBrighterFatter`
1557 The brighter-fatter kernel.
1558 brighterFatterApplyGain : `bool`
1559 Apply the gain to convert the image to electrons?
1560 bfGains : `dict`
1561 The gains to use if brighterFatterApplyGain = True.
1563 Yields
1564 ------
1565 exp : `lsst.afw.image.Exposure`
1566 The flat and dark corrected exposure.
1567 """
1568 interpExp = ccdExposure.clone()
1570 # We need to interpolate before we do B-F. Note that
1571 # brighterFatterFwhmForInterpolation is currently unused.
1572 isrFunctions.interpolateFromMask(
1573 maskedImage=interpExp.getMaskedImage(),
1574 fwhm=self.config.brighterFatterFwhmForInterpolation,
1575 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1576 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1577 useLegacyInterp=self.config.useLegacyInterp,
1578 )
1579 bfExp = interpExp.clone()
1581 ccdExposure = electrostaticBrighterFatterCorrection(
1582 bfExp,
1583 electroBfDistortionMatrix,
1584 brighterFatterApplyGain,
1585 bfGains,
1586 )
1588 # Applying the brighter-fatter correction applies a
1589 # convolution to the science image. At the edges this
1590 # convolution may not have sufficient valid pixels to
1591 # produce a valid correction. Mark pixels within the size
1592 # of the brighter-fatter kernel as EDGE to warn of this
1593 # fact.
1594 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1595 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(electroBfDistortionMatrix.aMatrix.shape) // 2,
1596 maskPlane="EDGE")
1598 if self.config.brighterFatterMaskGrowSize > 0:
1599 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1600 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1601 isrFunctions.growMasks(ccdExposure.getMask(),
1602 radius=self.config.brighterFatterMaskGrowSize,
1603 maskNameList=maskPlane,
1604 maskValue=maskPlane)
1606 return ccdExposure
1608 def applyFluxConservingBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel,
1609 brighterFatterApplyGain, bfGains):
1610 """Apply a brighter fatter correction to the image using the
1611 method defined in Coulton et al. 2019 with flux-conserving
1612 corrections.
1614 Note that this correction requires that the image is in units
1615 electrons.
1617 Parameters
1618 ----------
1619 ccdExposure : `lsst.afw.image.Exposure`
1620 Exposure to process.
1621 flat : `lsst.afw.image.Exposure`
1622 Flat exposure the same size as ``exp``.
1623 dark : `lsst.afw.image.Exposure`, optional
1624 Dark exposure the same size as ``exp``.
1625 bfKernel : `lsst.ip.isr.BrighterFatterKernel`
1626 The brighter-fatter kernel.
1627 brighterFatterApplyGain : `bool`
1628 Apply the gain to convert the image to electrons?
1629 bfGains : `dict`
1630 The gains to use if brighterFatterApplyGain = True.
1632 Yields
1633 ------
1634 exp : `lsst.afw.image.Exposure`
1635 The flat and dark corrected exposure.
1636 """
1637 interpExp = ccdExposure.clone()
1639 # We need to interpolate before we do B-F. Note that
1640 # brighterFatterFwhmForInterpolation is currently unused.
1641 isrFunctions.interpolateFromMask(
1642 maskedImage=interpExp.getMaskedImage(),
1643 fwhm=self.config.brighterFatterFwhmForInterpolation,
1644 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1645 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1646 useLegacyInterp=self.config.useLegacyInterp,
1647 )
1649 bfExp = interpExp.clone()
1651 bfResults = fluxConservingBrighterFatterCorrection(
1652 bfExp, bfKernel,
1653 self.config.brighterFatterMaxIter,
1654 self.config.brighterFatterThreshold,
1655 brighterFatterApplyGain,
1656 bfGains,
1657 )
1659 bfCorrIters = bfResults[1]
1660 if bfCorrIters == self.config.brighterFatterMaxIter:
1661 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1662 bfResults[0])
1663 else:
1664 self.log.info("Finished brighter-fatter correction in %d iterations.",
1665 bfResults[1])
1667 image = ccdExposure.getMaskedImage().getImage()
1668 bfCorr = bfExp.getMaskedImage().getImage()
1669 bfCorr -= interpExp.getMaskedImage().getImage()
1670 image += bfCorr
1672 # Applying the brighter-fatter correction applies a
1673 # convolution to the science image. At the edges this
1674 # convolution may not have sufficient valid pixels to
1675 # produce a valid correction. Mark pixels within the size
1676 # of the brighter-fatter kernel as EDGE to warn of this
1677 # fact.
1678 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1679 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1680 maskPlane="EDGE")
1682 if self.config.brighterFatterMaskGrowSize > 0:
1683 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1684 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1685 isrFunctions.growMasks(ccdExposure.getMask(),
1686 radius=self.config.brighterFatterMaskGrowSize,
1687 maskNameList=maskPlane,
1688 maskValue=maskPlane)
1690 # Set the metadata here.
1691 ccdExposure.metadata["LSST ISR BF ITERS"] = bfCorrIters
1693 return ccdExposure
1695 def applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, brighterFatterApplyGain,
1696 bfGains):
1697 """Apply a brighter fatter correction to the image using the
1698 method defined in Coulton et al. 2019.
1700 Note that this correction requires that the image is in units
1701 electrons.
1703 Parameters
1704 ----------
1705 ccdExposure : `lsst.afw.image.Exposure`
1706 Exposure to process.
1707 flat : `lsst.afw.image.Exposure`
1708 Flat exposure the same size as ``exp``.
1709 dark : `lsst.afw.image.Exposure`, optional
1710 Dark exposure the same size as ``exp``.
1711 bfKernel : `lsst.ip.isr.BrighterFatterKernel`
1712 The brighter-fatter kernel.
1713 brighterFatterApplyGain : `bool`
1714 Apply the gain to convert the image to electrons?
1715 bfGains : `dict`
1716 The gains to use if brighterFatterApplyGain = True.
1718 Yields
1719 ------
1720 exp : `lsst.afw.image.Exposure`
1721 The flat and dark corrected exposure.
1722 """
1723 interpExp = ccdExposure.clone()
1725 # We need to interpolate before we do B-F. Note that
1726 # brighterFatterFwhmForInterpolation is currently unused.
1727 isrFunctions.interpolateFromMask(
1728 maskedImage=interpExp.getMaskedImage(),
1729 fwhm=self.config.brighterFatterFwhmForInterpolation,
1730 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1731 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1732 useLegacyInterp=self.config.useLegacyInterp,
1733 )
1734 bfExp = interpExp.clone()
1735 bfResults = brighterFatterCorrection(bfExp, bfKernel,
1736 self.config.brighterFatterMaxIter,
1737 self.config.brighterFatterThreshold,
1738 brighterFatterApplyGain,
1739 bfGains)
1740 bfCorrIters = bfResults[1]
1741 if bfCorrIters == self.config.brighterFatterMaxIter:
1742 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1743 bfResults[0])
1744 else:
1745 self.log.info("Finished brighter-fatter correction in %d iterations.",
1746 bfResults[1])
1748 image = ccdExposure.getMaskedImage().getImage()
1749 bfCorr = bfExp.getMaskedImage().getImage()
1750 bfCorr -= interpExp.getMaskedImage().getImage()
1751 image += bfCorr
1753 # Applying the brighter-fatter correction applies a
1754 # convolution to the science image. At the edges this
1755 # convolution may not have sufficient valid pixels to
1756 # produce a valid correction. Mark pixels within the size
1757 # of the brighter-fatter kernel as EDGE to warn of this
1758 # fact.
1759 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1760 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1761 maskPlane="EDGE")
1763 if self.config.brighterFatterMaskGrowSize > 0:
1764 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1765 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1766 isrFunctions.growMasks(ccdExposure.getMask(),
1767 radius=self.config.brighterFatterMaskGrowSize,
1768 maskNameList=maskPlane,
1769 maskValue=maskPlane)
1771 # Set the metadata here.
1772 ccdExposure.metadata["LSST ISR BF ITERS"] = bfCorrIters
1774 return ccdExposure
1776 def darkCorrection(self, exposure, darkExposure, invert=False):
1777 """Apply dark correction in place.
1779 Parameters
1780 ----------
1781 exposure : `lsst.afw.image.Exposure`
1782 Exposure to process.
1783 darkExposure : `lsst.afw.image.Exposure`
1784 Dark exposure of the same size as ``exposure``.
1785 invert : `Bool`, optional
1786 If True, re-add the dark to an already corrected image.
1788 Raises
1789 ------
1790 RuntimeError
1791 Raised if either ``exposure`` or ``darkExposure`` do not
1792 have their dark time defined.
1794 See Also
1795 --------
1796 lsst.ip.isr.isrFunctions.darkCorrection
1797 """
1798 expScale = exposure.visitInfo.darkTime
1799 if math.isnan(expScale):
1800 raise RuntimeError("Exposure darktime is NAN.")
1801 if darkExposure.visitInfo is not None \
1802 and not math.isnan(darkExposure.visitInfo.darkTime):
1803 darkScale = darkExposure.visitInfo.darkTime
1804 else:
1805 # DM-17444: darkExposure.visitInfo is None
1806 # so darkTime does not exist.
1807 self.log.warning("darkExposure.visitInfo does not exist. Using darkScale = 1.0.")
1808 darkScale = 1.0
1810 isrFunctions.darkCorrection(
1811 maskedImage=exposure.maskedImage,
1812 darkMaskedImage=darkExposure.maskedImage,
1813 expScale=expScale,
1814 darkScale=darkScale,
1815 invert=invert,
1816 )
1818 @staticmethod
1819 def extractCalibDate(calib):
1820 """Extract common calibration metadata values that will be written to
1821 output header.
1823 Parameters
1824 ----------
1825 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1826 Calibration to pull date information from.
1828 Returns
1829 -------
1830 dateString : `str`
1831 Calibration creation date string to add to header.
1832 """
1833 if hasattr(calib, "getMetadata"):
1834 if 'CALIB_CREATION_DATE' in calib.metadata:
1835 return " ".join((calib.metadata.get("CALIB_CREATION_DATE", "Unknown"),
1836 calib.metadata.get("CALIB_CREATION_TIME", "Unknown")))
1837 else:
1838 return " ".join((calib.metadata.get("CALIB_CREATE_DATE", "Unknown"),
1839 calib.metadata.get("CALIB_CREATE_TIME", "Unknown")))
1840 else:
1841 return "Unknown Unknown"
1843 def compareUnits(self, calibMetadata, calibName):
1844 """Compare units from calibration to ISR units.
1846 This compares calibration units (adu or electron) to whether
1847 doApplyGain is set.
1849 Parameters
1850 ----------
1851 calibMetadata : `lsst.daf.base.PropertyList`
1852 Calibration metadata from header.
1853 calibName : `str`
1854 Calibration name for log message.
1855 """
1856 calibUnits = calibMetadata.get("LSST ISR UNITS", "adu")
1857 isrUnits = "electron" if self.config.doApplyGains else "adu"
1858 if calibUnits != isrUnits:
1859 if self.config.doRaiseOnCalibMismatch:
1860 raise RuntimeError(
1861 "Unit mismatch: isr has %s units but %s has %s units",
1862 isrUnits,
1863 calibName,
1864 calibUnits,
1865 )
1866 else:
1867 self.log.warning(
1868 "Unit mismatch: isr has %s units but %s has %s units",
1869 isrUnits,
1870 calibName,
1871 calibUnits,
1872 )
1874 def convertIntToFloat(self, exposure):
1875 """Convert exposure image from uint16 to float.
1877 If the exposure does not need to be converted, the input is
1878 immediately returned. For exposures that are converted to use
1879 floating point pixels, the variance is set to unity and the
1880 mask to zero.
1882 Parameters
1883 ----------
1884 exposure : `lsst.afw.image.Exposure`
1885 The raw exposure to be converted.
1887 Returns
1888 -------
1889 newexposure : `lsst.afw.image.Exposure`
1890 The input ``exposure``, converted to floating point pixels.
1892 Raises
1893 ------
1894 RuntimeError
1895 Raised if the exposure type cannot be converted to float.
1897 """
1898 if isinstance(exposure, afwImage.ExposureF):
1899 # Nothing to be done
1900 self.log.debug("Exposure already of type float.")
1901 return exposure
1902 if not hasattr(exposure, "convertF"):
1903 raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
1905 newexposure = exposure.convertF()
1906 newexposure.variance[:] = 1
1907 newexposure.mask[:] = 0x0
1909 return newexposure
1911 def ditherCounts(self, exposure, detectorConfig, fallbackSeed=12345):
1912 """Dither the counts in the exposure.
1914 Parameters
1915 ----------
1916 exposure : `lsst.afw.image.Exposure`
1917 The raw exposure to be dithered.
1918 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1919 Configuration for overscan/etc for this detector.
1920 fallbackSeed : `int`, optional
1921 Random seed to fall back to if exposure.getInfo().getId() is
1922 not set.
1923 """
1924 if detectorConfig.integerDitherMode == "NONE":
1925 # Nothing to do here.
1926 return
1928 # This ID is a unique combination of {exposure, detector} for a raw
1929 # image as we have here. We additionally need to take the lower
1930 # 32 bits to be used as a random seed.
1931 if exposure.info.id is not None:
1932 seed = exposure.info.id & 0xFFFFFFFF
1933 else:
1934 seed = fallbackSeed
1935 self.log.warning("No exposure ID found; using fallback random seed.")
1937 self.log.info("Seeding dithering random number generator with %d.", seed)
1938 rng = numpy.random.RandomState(seed=seed)
1940 if detectorConfig.integerDitherMode == "POSITIVE":
1941 low = 0.0
1942 high = 1.0
1943 elif detectorConfig.integerDitherMode == "NEGATIVE":
1944 low = -1.0
1945 high = 0.0
1946 elif detectorConfig.integerDitherMode == "SYMMETRIC":
1947 low = -0.5
1948 high = 0.5
1949 else:
1950 raise RuntimeError("Invalid config")
1952 exposure.image.array[:, :] += rng.uniform(low=low, high=high, size=exposure.image.array.shape)
1954 def checkBssVoltage(self, exposure):
1955 """Check the back-side bias voltage to see if the detector is on.
1957 Parameters
1958 ----------
1959 exposure : `lsst.afw.image.ExposureF`
1960 Input exposure.
1962 Raises
1963 ------
1964 `UnprocessableDataError` if voltage is off.
1965 """
1966 voltage = exposure.metadata.get(self.config.bssVoltageKeyword, None)
1967 if voltage is None or not numpy.isfinite(voltage):
1968 self.log.warning(
1969 "Back-side bias voltage %s not found in metadata.",
1970 self.config.bssVoltageKeyword,
1971 )
1972 return
1974 hv = exposure.metadata.get(self.config.hvBiasKeyword, None)
1975 if hv is None:
1976 self.log.warning(
1977 "HV bias on %s not found in metadata.",
1978 self.config.hvBiasKeyword,
1979 )
1980 return
1982 if voltage < self.config.bssVoltageMinimum or hv == "OFF":
1983 detector = exposure.getDetector()
1984 raise UnprocessableDataError(
1985 f"Back-side bias voltage is turned off for {detector.getName()}; skipping ISR.",
1986 )
1988 @deprecated(
1989 reason=(
1990 "makeBinnedImages is no longer used. "
1991 "Please subtask lsst.ip.isr.BinImageDataTask instead."
1992 ),
1993 version="v28", category=FutureWarning
1994 )
1995 def makeBinnedImages(self, exposure):
1996 """Make visualizeVisit style binned exposures.
1998 Parameters
1999 ----------
2000 exposure : `lsst.afw.image.Exposure`
2001 Exposure to bin.
2003 Returns
2004 -------
2005 bin1 : `lsst.afw.image.Exposure`
2006 Binned exposure using binFactor1.
2007 bin2 : `lsst.afw.image.Exposure`
2008 Binned exposure using binFactor2.
2009 """
2010 mi = exposure.getMaskedImage()
2012 bin1 = afwMath.binImage(mi, self.config.binFactor1)
2013 bin2 = afwMath.binImage(mi, self.config.binFactor2)
2015 bin1 = afwImage.makeExposure(bin1)
2016 bin2 = afwImage.makeExposure(bin2)
2018 bin1.setInfo(exposure.getInfo())
2019 bin2.setInfo(exposure.getInfo())
2021 return bin1, bin2
2023 def run(
2024 self,
2025 ccdExposure,
2026 *,
2027 dnlLUT=None,
2028 bias=None,
2029 deferredChargeCalib=None,
2030 linearizer=None,
2031 ptc=None,
2032 gainCorrection=None,
2033 crosstalk=None,
2034 defects=None,
2035 bfKernel=None,
2036 electroBfDistortionMatrix=None,
2037 dark=None,
2038 flat=None,
2039 camera=None,
2040 ):
2041 """Run the IsrTaskLSST task.
2043 Parameters
2044 ----------
2045 ccdExposure : `lsst.afw.image.Exposure`
2046 Exposure to run ISR.
2047 dnlLUT : `None`, optional
2048 DNL lookup table; placeholder, unused.
2049 bias : `lsst.afw.image.Exposure`, optional
2050 Bias frame.
2051 deferredChargeCalib : `lsst.ip.isr.DeferredChargeCalib`, optional
2052 Deferred charge calibration.
2053 linearizer : `lsst.ip.isr.Linearizer`, optional
2054 Linearizer calibration.
2055 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
2056 PTC dataset.
2057 gainCorrection : `lsst.ip.isr.GainCorrection`, optional
2058 Gain correction dataset.
2059 crosstalk : `lsst.ip.isr.CrosstalkCalib`, optional
2060 Crosstalk calibration dataset.
2061 defects : `lsst.ip.isr.Defects`, optional
2062 Defects dataset.
2063 bfKernel : `lsst.ip.isr.BrighterFatterKernel`, optional
2064 Brighter-fatter kernel dataset.
2065 dark : `lsst.afw.image.Exposure`, optional
2066 Dark frame.
2067 flat : `lsst.afw.image.Exposure`, optional
2068 Flat-field frame.
2069 camera : `lsst.afw.cameraGeom.Camera`, optional
2070 Camera object.
2072 Returns
2073 -------
2074 result : `lsst.pipe.base.Struct`
2075 Struct with fields:
2076 ``exposure``: `lsst.afw.image.Exposure`
2077 Calibrated exposure.
2078 ``outputBin1Exposure``: `lsst.afw.image.Exposure`
2079 Binned exposure (bin1 config).
2080 ``outputBin2Exposure``: `lsst.afw.image.Exposure`
2081 Binned exposure (bin2 config).
2082 ``outputExposure``: `lsst.afw.image.Exposure`
2083 Calibrated exposure (same as ``exposure``).
2084 ``outputStatistics``: `lsst.ip.isr.isrStatistics`
2085 Calibrated exposure statistics.
2086 """
2087 detector = ccdExposure.getDetector()
2089 overscanDetectorConfig = self.config.overscanCamera.getOverscanDetectorConfig(detector)
2091 if self.config.doBootstrap:
2092 if ptc is not None:
2093 self.log.warning("Task configured with doBootstrap=True. Ignoring provided PTC.")
2094 ptc = None
2095 else:
2096 if self.config.useGainsFrom == "LINEARIZER":
2097 if linearizer is None:
2098 raise RuntimeError("doBootstrap==False and useGainsFrom == 'LINEARIZER' but "
2099 "no linearizer provided.")
2100 elif self.config.useGainsFrom == "PTC":
2101 if ptc is None:
2102 raise RuntimeError("doBootstrap==False and useGainsFrom == 'PTC' but no PTC provided.")
2104 # Validation step: check inputs match exposure configuration.
2105 exposureMetadata = ccdExposure.metadata
2106 doRaise = self.config.doRaiseOnCalibMismatch
2107 keywords = self.config.cameraKeywordsToCompare
2108 if not self.config.doBootstrap:
2109 if self.config.useGainsFrom == "LINEARIZER":
2110 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer,
2111 "LINEARIZER", log=self.log)
2112 elif self.config.useGainsFrom == "PTC":
2113 compareCameraKeywords(doRaise, keywords, exposureMetadata, ptc, "PTC",
2114 log=self.log)
2115 # Note that doCorrectGains can be True without a gainCorrection.
2116 if self.config.doCorrectGains and gainCorrection is not None:
2117 compareCameraKeywords(
2118 doRaise,
2119 keywords,
2120 exposureMetadata,
2121 gainCorrection,
2122 "gain_correction",
2123 log=self.log,
2124 )
2125 else:
2126 if self.config.doCorrectGains:
2127 raise RuntimeError("doCorrectGains is True but no ptc provided.")
2128 if self.config.doDiffNonLinearCorrection:
2129 if dnlLUT is None:
2130 raise RuntimeError("doDiffNonLinearCorrection is True but no dnlLUT provided.")
2131 compareCameraKeywords(doRaise, keywords, exposureMetadata, dnlLUT, "dnlLUT", log=self.log)
2132 if self.config.doLinearize:
2133 if linearizer is None:
2134 raise RuntimeError("doLinearize is True but no linearizer provided.")
2135 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer, "linearizer", log=self.log)
2136 if self.config.doBias:
2137 if bias is None:
2138 raise RuntimeError("doBias is True but no bias provided.")
2139 compareCameraKeywords(doRaise, keywords, exposureMetadata, bias, "bias", log=self.log)
2140 self.compareUnits(bias.metadata, "bias")
2141 if self.config.doCrosstalk:
2142 if crosstalk is None:
2143 raise RuntimeError("doCrosstalk is True but no crosstalk provided.")
2144 compareCameraKeywords(doRaise, keywords, exposureMetadata, crosstalk, "crosstalk", log=self.log)
2145 if self.config.doDeferredCharge:
2146 if deferredChargeCalib is None:
2147 raise RuntimeError("doDeferredCharge is True but no deferredChargeCalib provided.")
2148 compareCameraKeywords(
2149 doRaise,
2150 keywords,
2151 exposureMetadata,
2152 deferredChargeCalib,
2153 "CTI",
2154 log=self.log,
2155 )
2156 if self.config.doDefect:
2157 if defects is None:
2158 raise RuntimeError("doDefect is True but no defects provided.")
2159 compareCameraKeywords(doRaise, keywords, exposureMetadata, defects, "defects", log=self.log)
2160 if self.config.doDark:
2161 if dark is None:
2162 raise RuntimeError("doDark is True but no dark frame provided.")
2163 compareCameraKeywords(doRaise, keywords, exposureMetadata, dark, "dark", log=self.log)
2164 self.compareUnits(bias.metadata, "dark")
2165 if self.config.doBrighterFatter:
2166 if self.config.brighterFatterCorrectionMethod == "ASTIER23":
2167 if electroBfDistortionMatrix is None:
2168 raise RuntimeError("Must supply an electroBfDistortionMatrix if BF "
2169 "correction method is ASTIER23.")
2170 compareCameraKeywords(doRaise, keywords, exposureMetadata,
2171 electroBfDistortionMatrix, "bf", log=self.log)
2172 elif self.config.brighterFatterCorrectionMethod in ["COULTON18", "COULTON18_FLUX_CONSERVING"]:
2173 if bfKernel is None:
2174 raise RuntimeError("Must supply an kernel if BF correction method is COULTON*.")
2175 compareCameraKeywords(doRaise, keywords, exposureMetadata, bfKernel, "bf", log=self.log)
2176 else:
2177 raise ValueError("%s is not a known brighter-fatter correction "
2178 "method." % self.config.brighterFatterCorrectionMethod)
2179 if self.config.doFlat:
2180 if flat is None:
2181 raise RuntimeError("doFlat is True but no flat provided.")
2182 compareCameraKeywords(doRaise, keywords, exposureMetadata, flat, "flat", log=self.log)
2184 if self.config.doSaturation:
2185 if self.config.defaultSaturationSource in ["PTCTURNOFF",]:
2186 if ptc is None:
2187 raise RuntimeError(
2188 "doSaturation is True and defaultSaturationSource is "
2189 f"{self.config.defaultSaturationSource}, but no ptc provided."
2190 )
2191 if self.config.doSuspect:
2192 if self.config.defaultSuspectSource in ["PTCTURNOFF",]:
2193 if ptc is None:
2194 raise RuntimeError(
2195 "doSuspect is True and defaultSuspectSource is "
2196 f"{self.config.defaultSuspectSource}, but no ptc provided."
2197 )
2199 if self.config.doCheckUnprocessableData and self.config.bssVoltageMinimum > 0.0:
2200 self.checkBssVoltage(ccdExposure)
2202 # FIXME: Make sure that if linearity is done then it is matched
2203 # with the right PTC.
2205 # We keep track of units: start in adu.
2206 exposureMetadata["LSST ISR UNITS"] = "adu"
2207 exposureMetadata["LSST ISR GAINCORRECTION APPLIED"] = False
2208 exposureMetadata["LSST ISR CROSSTALK APPLIED"] = False
2209 exposureMetadata["LSST ISR OVERSCANLEVEL CHECKED"] = False
2210 exposureMetadata["LSST ISR NOISE CHECKED"] = False
2211 exposureMetadata["LSST ISR LINEARIZER APPLIED"] = False
2212 exposureMetadata["LSST ISR CTI APPLIED"] = False
2213 exposureMetadata["LSST ISR BIAS APPLIED"] = False
2214 exposureMetadata["LSST ISR DARK APPLIED"] = False
2215 exposureMetadata["LSST ISR BF APPLIED"] = False
2216 exposureMetadata["LSST ISR FLAT APPLIED"] = False
2217 exposureMetadata["LSST ISR DEFECTS APPLIED"] = False
2219 if self.config.doBootstrap:
2220 self.log.info("Configured using doBootstrap=True; using gain of 1.0 (adu units)")
2221 ptc = PhotonTransferCurveDataset([amp.getName() for amp in detector], "NOMINAL_PTC", 1)
2222 for amp in detector:
2223 ptc.gain[amp.getName()] = 1.0
2224 ptc.noise[amp.getName()] = 0.0
2225 elif self.config.useGainsFrom == "LINEARIZER":
2226 self.log.info("Using gains from linearizer.")
2227 # Create a dummy ptc object to hold the gains from the linearizer.
2228 ptc = PhotonTransferCurveDataset([amp.getName() for amp in detector], "NOMINAL_PTC", 1)
2229 for amp in detector:
2230 ptc.gain[amp.getName()] = linearizer.inputGain[amp.getName()]
2231 ptc.noise[amp.getName()] = 0.0
2233 exposureMetadata["LSST ISR BOOTSTRAP"] = self.config.doBootstrap
2235 # Choose the gains to use
2236 gains = ptc.gain
2238 # And check if we have configured gains to override. This is
2239 # also a warning, since it should not be typical usage.
2240 for amp in detector:
2241 if not math.isnan(gain := overscanDetectorConfig.getOverscanAmpConfig(amp).gain):
2242 gains[amp.getName()] = gain
2243 self.log.warning(
2244 "Overriding gain for amp %s with configured value of %.3f.",
2245 amp.getName(),
2246 gain,
2247 )
2249 # First we convert the exposure to floating point values
2250 # (if necessary).
2251 self.log.debug("Converting exposure to floating point values.")
2252 ccdExposure = self.convertIntToFloat(ccdExposure)
2254 # Then we mark which amplifiers are completely bad from defects.
2255 badAmpDict = self.maskFullAmplifiers(ccdExposure, detector, defects, gains=gains)
2257 self.checkAllBadAmps(badAmpDict, detector)
2259 # Now we go through ISR steps.
2261 # Differential non-linearity correction.
2262 # Units: adu
2263 if self.config.doDiffNonLinearCorrection:
2264 self.diffNonLinearCorrection(ccdExposure, dnlLUT)
2266 # Dither the integer counts.
2267 # Input units: integerized adu
2268 # Output units: floating-point adu
2269 self.ditherCounts(ccdExposure, overscanDetectorConfig)
2271 # Serial overscan correction.
2272 # Input units: adu
2273 # Output units: adu
2274 if overscanDetectorConfig.doAnySerialOverscan:
2275 serialOverscans = self.overscanCorrection(
2276 "SERIAL",
2277 overscanDetectorConfig,
2278 detector,
2279 badAmpDict,
2280 ccdExposure,
2281 )
2283 if self.config.doBootstrap or self.config.useGainsFrom == "LINEARIZER":
2284 # Get the empirical read noise
2285 for amp, serialOverscan in zip(detector, serialOverscans):
2286 if serialOverscan is None:
2287 ptc.noise[amp.getName()] = 0.0
2288 else:
2289 # All PhotonTransferCurveDataset objects should contain
2290 # noise attributes in units of electrons. The read
2291 # noise measured from overscans is always in adu, so we
2292 # scale it by the gain.
2293 # Note that in bootstrap mode, these gains will always
2294 # be 1.0, but we put this conversion here for clarity.
2295 ptc.noise[amp.getName()] = serialOverscan.residualSigma * gains[amp.getName()]
2296 else:
2297 serialOverscans = [None]*len(detector)
2299 # After serial overscan correction, we can mask SATURATED and
2300 # SUSPECT pixels. This updates badAmpDict if any amplifier
2301 # is fully saturated after serial overscan correction.
2303 # The saturation is currently assumed to be recorded in
2304 # overscan-corrected adu.
2305 badAmpDict = self.maskSaturatedPixels(
2306 badAmpDict,
2307 ccdExposure,
2308 detector,
2309 overscanDetectorConfig,
2310 ptc=ptc,
2311 )
2313 self.checkAllBadAmps(badAmpDict, detector)
2315 if self.config.doCorrectGains and gainCorrection is not None:
2316 self.log.info("Correcting gains based on input GainCorrection.")
2317 gainCorrection.correctGains(gains, exposure=ccdExposure)
2318 exposureMetadata["LSST ISR GAINCORRECTION APPLIED"] = True
2319 elif self.config.doCorrectGains:
2320 self.log.info("Skipping gain correction because no GainCorrection available.")
2322 # Do gain normalization.
2323 # Input units: adu
2324 # Output units: electron
2325 if self.config.doApplyGains:
2326 self.log.info("Using gain values to convert from adu to electron units.")
2327 isrFunctions.applyGains(ccdExposure, normalizeGains=False, ptcGains=gains, isTrimmed=False)
2328 # The units are now electron.
2329 exposureMetadata["LSST ISR UNITS"] = "electron"
2331 # Update the saturation units in the metadata if there.
2332 # These will always have the same units as the image.
2333 for amp in detector:
2334 ampName = amp.getName()
2335 if (key := f"LSST ISR SATURATION LEVEL {ampName}") in exposureMetadata:
2336 exposureMetadata[key] *= gains[ampName]
2337 if (key := f"LSST ISR SUSPECT LEVEL {ampName}") in exposureMetadata:
2338 exposureMetadata[key] *= gains[ampName]
2340 # Record gain and read noise in header.
2341 metadata = ccdExposure.metadata
2342 metadata["LSST ISR READNOISE UNITS"] = "electron"
2343 metadata["LSST ISR GAIN SOURCE"] = self.config.useGainsFrom
2344 for amp in detector:
2345 # This includes any gain correction (if applied).
2346 metadata[f"LSST ISR GAIN {amp.getName()}"] = gains[amp.getName()]
2348 # At this stage, the read noise is always in electrons.
2349 noise = ptc.noise[amp.getName()]
2350 metadata[f"LSST ISR READNOISE {amp.getName()}"] = noise
2352 # Do crosstalk correction in the full region.
2353 # Output units: electron (adu if doBootstrap=True)
2354 if self.config.doCrosstalk:
2355 self.log.info("Applying crosstalk corrections to full amplifier region.")
2356 if self.config.doBootstrap and numpy.any(crosstalk.fitGains != 0):
2357 crosstalkGains = None
2358 else:
2359 crosstalkGains = gains
2360 self.crosstalk.run(
2361 ccdExposure,
2362 crosstalk=crosstalk,
2363 gains=crosstalkGains,
2364 fullAmplifier=True,
2365 badAmpDict=badAmpDict,
2366 ignoreVariance=True,
2367 )
2368 ccdExposure.metadata["LSST ISR CROSSTALK APPLIED"] = True
2370 # After crosstalk, we check for amplifier noise and state changes.
2371 if numpy.isfinite(self.config.serialOverscanMedianShiftSigmaThreshold):
2372 badAmpDict = self.checkAmpOverscanLevel(badAmpDict, ccdExposure, ptc)
2373 ccdExposure.metadata["LSST ISR OVERSCANLEVEL CHECKED"] = True
2375 if numpy.isfinite(self.config.ampNoiseThreshold):
2376 badAmpDict = self.checkAmpNoise(badAmpDict, ccdExposure, ptc)
2377 ccdExposure.metadata["LSST ISR NOISE CHECKED"] = True
2379 if numpy.isfinite(self.config.serialOverscanMedianShiftSigmaThreshold) or \
2380 numpy.isfinite(self.config.ampNoiseThreshold):
2381 self.checkAllBadAmps(badAmpDict, detector)
2383 # Parallel overscan correction.
2384 # Output units: electron (adu if doBootstrap=True)
2385 parallelOverscans = None
2386 if overscanDetectorConfig.doAnyParallelOverscan:
2387 # At the moment we do not use the return values from this task.
2388 parallelOverscans = self.overscanCorrection(
2389 "PARALLEL",
2390 overscanDetectorConfig,
2391 detector,
2392 badAmpDict,
2393 ccdExposure,
2394 )
2396 # Linearity correction
2397 # Output units: electron (adu if doBootstrap=True)
2398 if self.config.doLinearize:
2399 self.log.info("Applying linearizer.")
2400 # The linearizer is in units of adu.
2401 # If our units are electron, then pass in the gains
2402 # for conversion.
2403 if exposureMetadata["LSST ISR UNITS"] == "electron":
2404 linearityGains = gains
2405 else:
2406 linearityGains = None
2407 linearizer.applyLinearity(
2408 image=ccdExposure.image,
2409 detector=detector,
2410 log=self.log,
2411 gains=linearityGains,
2412 )
2413 ccdExposure.metadata["LSST ISR LINEARIZER APPLIED"] = True
2415 # Serial CTI (deferred charge) correction
2416 # This will be performed in electron units
2417 # Output units: same as input units
2418 if self.config.doDeferredCharge:
2419 if self.config.doBootstrap:
2420 self.log.info("Applying deferred charge correction with doBootstrap=True: "
2421 "will need to use deferredChargeCalib.inputGain to apply "
2422 "CTI correction in electron units.")
2423 deferredChargeGains = deferredChargeCalib.inputGain
2424 if numpy.all(numpy.isnan(list(deferredChargeGains.values()))):
2425 self.log.warning("All gains contained in the deferredChargeCalib are "
2426 "NaN, approximating with gain of 1.0.")
2427 deferredChargeGains = gains
2428 else:
2429 deferredChargeGains = gains
2430 self.deferredChargeCorrection.run(
2431 ccdExposure,
2432 deferredChargeCalib,
2433 gains=deferredChargeGains,
2434 )
2435 ccdExposure.metadata["LSST ISR CTI APPLIED"] = True
2437 # Save the untrimmed version for later statistics,
2438 # which still contains the overscan information
2439 untrimmedCcdExposure = None
2440 if self.config.isrStats.doCtiStatistics:
2441 untrimmedCcdExposure = ccdExposure.clone()
2443 # Assemble/trim
2444 # Output units: electron (adu if doBootstrap=True)
2445 if self.config.doAssembleCcd:
2446 self.log.info("Assembling CCD from amplifiers.")
2447 ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
2449 if self.config.expectWcs and not ccdExposure.getWcs():
2450 self.log.warning("No WCS found in input exposure.")
2452 # E2V edge bleed
2453 if self.config.doE2VEdgeBleedMask and detector.getPhysicalType() == "E2V":
2454 isrFunctions.maskE2VEdgeBleed(
2455 exposure=ccdExposure,
2456 e2vEdgeBleedSatMinArea=self.config.e2vEdgeBleedSatMinArea,
2457 e2vEdgeBleedSatMaxArea=self.config.e2vEdgeBleedSatMaxArea,
2458 e2vEdgeBleedYMax=self.config.e2vEdgeBleedYMax,
2459 saturatedMaskName=self.config.saturatedMaskName,
2460 log=self.log,
2461 )
2463 # ITL Dip Masking
2464 for maskPlane in self.config.itlDipMaskPlanes:
2465 if maskPlane not in ccdExposure.mask.getMaskPlaneDict():
2466 self.log.info("Adding %s mask plane to image.", maskPlane)
2467 ccdExposure.mask.addMaskPlane(maskPlane)
2469 if self.config.doITLDipMask:
2470 isrFunctions.maskITLDip(
2471 exposure=ccdExposure,
2472 detectorConfig=overscanDetectorConfig,
2473 log=self.log,
2474 maskPlaneNames=self.config.itlDipMaskPlanes,
2475 )
2477 if (self.config.doITLSatSagMask or self.config.doITLEdgeBleedMask) \
2478 and detector.getPhysicalType() == 'ITL':
2479 self.maskITLSatEdgesAndColumns(exposure=ccdExposure,
2480 badAmpDict=badAmpDict)
2482 # Bias subtraction
2483 # Output units: electron (adu if doBootstrap=True)
2484 if self.config.doBias:
2485 self.log.info("Applying bias correction.")
2486 # Bias frame and ISR unit consistency is checked at the top of
2487 # the run method.
2488 isrFunctions.biasCorrection(ccdExposure.maskedImage, bias.maskedImage)
2489 ccdExposure.metadata["LSST ISR BIAS APPLIED"] = True
2491 # Dark subtraction
2492 # Output units: electron (adu if doBootstrap=True)
2493 if self.config.doDark:
2494 self.log.info("Applying dark subtraction.")
2495 # Dark frame and ISR unit consistency is checked at the top of
2496 # the run method.
2497 self.darkCorrection(ccdExposure, dark)
2498 ccdExposure.metadata["LSST ISR DARK APPLIED"] = True
2500 # Defect masking
2501 # Masking block (defects, NAN pixels and trails).
2502 # Saturated and suspect pixels have already been masked.
2503 # Output units: electron (adu if doBootstrap=True)
2504 if self.config.doDefect:
2505 self.log.info("Applying defect masking.")
2506 self.maskDefects(ccdExposure, defects)
2507 ccdExposure.metadata["LSST ISR DEFECTS APPLIED"] = True
2509 self.log.info("Adding UNMASKEDNAN mask plane to image.")
2510 ccdExposure.mask.addMaskPlane("UNMASKEDNAN")
2511 if self.config.doNanMasking:
2512 self.log.info("Masking non-finite (NAN, inf) value pixels.")
2513 self.maskNan(ccdExposure)
2515 if self.config.doWidenSaturationTrails:
2516 self.log.info("Widening saturation trails.")
2517 isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
2519 # Brighter-Fatter
2520 # Output units: electron (adu if doBootstrap=True)
2521 if self.config.doBrighterFatter:
2522 if self.config.brighterFatterCorrectionMethod == "ASTIER23":
2523 # Do the Astier et al. 2023 Brighter-Fatter correction
2524 self.log.info("Applying electrostatic brighter-fatter "
2525 "correction.")
2526 bfCalib = electroBfDistortionMatrix
2527 bfGains = electroBfDistortionMatrix.gain
2528 bfCorrFunction = self.applyElectrostaticBrighterFatterCorrection
2529 elif self.config.brighterFatterCorrectionMethod == "COULTON18":
2530 # Do the Coulton et al. 2018 Brighter-Fatter correction
2531 self.log.info("Applying brighter-fatter correction.")
2533 # Use the original gains used to compute the BFE
2534 bfCalib, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
2535 bfCorrFunction = self.applyBrighterFatterCorrection
2536 elif self.config.brighterFatterCorrectionMethod == "COULTON18_FLUX_CONSERVING":
2537 # Do the Coulton et al. 2018 correction with flux conserving
2538 # corrections
2539 bfCalib, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
2540 bfCorrFunction = self.applyFluxConservingBrighterFatterCorrection
2541 else:
2542 # This is an unknown BF correction type
2543 raise RuntimeError("Cannot perform brighter-fatter correction with unknown"
2544 "correction type.")
2546 # Needs to be done in electrons; applyBrighterFatterCorrection
2547 # will convert the image if necessary.
2548 if exposureMetadata["LSST ISR UNITS"] == "electron":
2549 brighterFatterApplyGain = False
2550 else:
2551 brighterFatterApplyGain = True
2553 if brighterFatterApplyGain and (ptc is not None) and (bfGains != gains):
2554 # The supplied ptc should be the same as the ptc used to
2555 # generate the bfKernel, in which case they will have the
2556 # same stored amp-keyed dictionary of gains. If not, there
2557 # is a mismatch in the calibrations being used. This should
2558 # not be always be a fatal error, but ideally, everything
2559 # should to be consistent.
2560 self.log.warning("Need to apply gain for brighter-fatter, but the stored"
2561 "gains in the kernel are not the same as the gains used"
2562 f"by {self.config.useGainsFrom}. Using the gains stored"
2563 "in the kernel.")
2565 # Apply the BF correction
2566 ccdExposure = bfCorrFunction(
2567 ccdExposure,
2568 flat,
2569 dark,
2570 bfCalib,
2571 brighterFatterApplyGain,
2572 bfGains,
2573 )
2575 ccdExposure.metadata["LSST ISR BF APPLIED"] = True
2576 ccdExposure.metadata["LSST ISR BF CORR METHOD"] = self.config.brighterFatterCorrectionMethod
2578 # Variance plane creation
2579 # Output units: electron (adu if doBootstrap=True)
2580 if self.config.doVariance:
2581 self.addVariancePlane(ccdExposure, detector)
2583 # Flat-fielding
2584 # This may move elsewhere, but this is the most convenient
2585 # location for simple flat-fielding for attractive backgrounds.
2586 # Output units: electron (adu if doBootstrap=True)
2587 if self.config.doFlat:
2588 self.log.info("Applying flat correction.")
2589 isrFunctions.flatCorrection(
2590 maskedImage=ccdExposure.maskedImage,
2591 flatMaskedImage=flat.maskedImage,
2592 scalingType=self.config.flatScalingType,
2593 userScale=self.config.flatUserScale,
2594 )
2596 # Copy over valid polygon from flat if it is
2597 # available, and set NO_DATA to 0.0 (which may
2598 # be inherited from the flat in the fully
2599 # vignetted region).
2600 if (validPolygon := flat.info.getValidPolygon()) is not None:
2601 ccdExposure.info.setValidPolygon(validPolygon)
2603 noData = (ccdExposure.mask.array & ccdExposure.mask.getPlaneBitMask("NO_DATA")) > 0
2604 ccdExposure.image.array[noData] = 0.0
2605 ccdExposure.variance.array[noData] = 0.0
2607 ccdExposure.metadata["LSST ISR FLAT APPLIED"] = True
2608 ccdExposure.metadata["LSST ISR FLAT SOURCE"] = flat.metadata.get("FLATSRC", "UNKNOWN")
2610 # Pixel values for masked regions are set here
2611 preInterpExp = None
2612 if self.config.doSaveInterpPixels:
2613 preInterpExp = ccdExposure.clone()
2615 if self.config.doSetBadRegions:
2616 self.log.info('Setting values in large contiguous bad regions.')
2617 self.setBadRegions(ccdExposure)
2619 if self.config.doInterpolate:
2620 self.log.info("Interpolating masked pixels.")
2621 isrFunctions.interpolateFromMask(
2622 maskedImage=ccdExposure.getMaskedImage(),
2623 fwhm=self.config.brighterFatterFwhmForInterpolation,
2624 growSaturatedFootprints=self.config.growSaturationFootprintSize,
2625 maskNameList=list(self.config.maskListToInterpolate),
2626 useLegacyInterp=self.config.useLegacyInterp,
2627 )
2629 # Calculate amp offset corrections within the CCD.
2630 if self.config.doAmpOffset:
2631 if self.config.ampOffset.doApplyAmpOffset:
2632 self.log.info("Measuring and applying amp offset corrections.")
2633 else:
2634 self.log.info("Measuring amp offset corrections only, without applying them.")
2635 self.ampOffset.run(ccdExposure)
2637 # Calculate standard image quality statistics
2638 if self.config.doStandardStatistics:
2639 metadata = ccdExposure.metadata
2640 for amp in detector:
2641 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
2642 ampName = amp.getName()
2643 metadata[f"LSST ISR MASK SAT {ampName}"] = isrFunctions.countMaskedPixels(
2644 ampExposure.getMaskedImage(),
2645 [self.config.saturatedMaskName]
2646 )
2647 metadata[f"LSST ISR MASK BAD {ampName}"] = isrFunctions.countMaskedPixels(
2648 ampExposure.getMaskedImage(),
2649 ["BAD"]
2650 )
2651 metadata[f"LSST ISR MASK SUSPECT {ampName}"] = isrFunctions.countMaskedPixels(
2652 ampExposure.getMaskedImage(),
2653 ["SUSPECT"],
2654 )
2655 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
2656 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP)
2658 metadata[f"LSST ISR FINAL MEAN {ampName}"] = qaStats.getValue(afwMath.MEAN)
2659 metadata[f"LSST ISR FINAL MEDIAN {ampName}"] = qaStats.getValue(afwMath.MEDIAN)
2660 metadata[f"LSST ISR FINAL STDEV {ampName}"] = qaStats.getValue(afwMath.STDEVCLIP)
2662 k1 = f"LSST ISR FINAL MEDIAN {ampName}"
2663 k2 = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
2664 if overscanDetectorConfig.doAnySerialOverscan and k1 in metadata and k2 in metadata:
2665 metadata[f"LSST ISR LEVEL {ampName}"] = metadata[k1] - metadata[k2]
2666 else:
2667 metadata[f"LSST ISR LEVEL {ampName}"] = numpy.nan
2669 # calculate additional statistics.
2670 outputStatistics = None
2671 if self.config.doCalculateStatistics:
2672 outputStatistics = self.isrStats.run(ccdExposure,
2673 untrimmedInputExposure=untrimmedCcdExposure,
2674 serialOverscanResults=serialOverscans,
2675 parallelOverscanResults=parallelOverscans,
2676 bias=bias, dark=dark, flat=flat,
2677 ptc=ptc, defects=defects).results
2679 # do image binning.
2680 outputBin1Exposure = None
2681 outputBin2Exposure = None
2682 if self.config.doBinnedExposures:
2683 self.log.info("Creating binned exposures.")
2684 outputBin1Exposure = self.binning.run(
2685 ccdExposure,
2686 binFactor=self.config.binFactor1,
2687 ).outputData
2688 outputBin2Exposure = self.binning.run(
2689 ccdExposure,
2690 binFactor=self.config.binFactor2,
2691 ).outputData
2693 return pipeBase.Struct(
2694 exposure=ccdExposure,
2696 outputBin1Exposure=outputBin1Exposure,
2697 outputBin2Exposure=outputBin2Exposure,
2699 preInterpExposure=preInterpExp,
2700 outputExposure=ccdExposure,
2701 outputStatistics=outputStatistics,
2702 )