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