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