Loading [MathJax]/extensions/tex2jax.js
LSST Applications g0fba68d861+05816baf74,g1ec0fe41b4+f536777771,g1fd858c14a+a9301854fb,g35bb328faa+fcb1d3bbc8,g4af146b050+a5c07d5b1d,g4d2262a081+6e5fcc2a4e,g53246c7159+fcb1d3bbc8,g56a49b3a55+9c12191793,g5a012ec0e7+3632fc3ff3,g60b5630c4e+ded28b650d,g67b6fd64d1+ed4b5058f4,g78460c75b0+2f9a1b4bcd,g786e29fd12+cf7ec2a62a,g8352419a5c+fcb1d3bbc8,g87b7deb4dc+7b42cf88bf,g8852436030+e5453db6e6,g89139ef638+ed4b5058f4,g8e3bb8577d+d38d73bdbd,g9125e01d80+fcb1d3bbc8,g94187f82dc+ded28b650d,g989de1cb63+ed4b5058f4,g9d31334357+ded28b650d,g9f33ca652e+50a8019d8c,gabe3b4be73+1e0a283bba,gabf8522325+fa80ff7197,gb1101e3267+d9fb1f8026,gb58c049af0+f03b321e39,gb665e3612d+2a0c9e9e84,gb89ab40317+ed4b5058f4,gcf25f946ba+e5453db6e6,gd6cbbdb0b4+bb83cc51f8,gdd1046aedd+ded28b650d,gde0f65d7ad+941d412827,ge278dab8ac+d65b3c2b70,ge410e46f29+ed4b5058f4,gf23fb2af72+b7cae620c0,gf5e32f922b+fcb1d3bbc8,gf67bdafdda+ed4b5058f4,w.2025.16
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
isrTaskLSST.py
Go to the documentation of this file.
1__all__ = ["IsrTaskLSST", "IsrTaskLSSTConfig"]
2
3import numpy
4import math
5
6from . import isrFunctions
7from . import isrQa
8from .defects import Defects
9
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.meas.algorithms.detection import SourceDetectionTask
18import lsst.afw.detection as afwDetection
19
20from .ampOffset import AmpOffsetTask
21from .binExposureTask import BinExposureTask
22from .overscan import SerialOverscanCorrectionTask, ParallelOverscanCorrectionTask
23from .overscanAmpConfig import OverscanCameraConfig
24from .assembleCcdTask import AssembleCcdTask
25from .deferredCharge import DeferredChargeTask
26from .crosstalk import CrosstalkTask
27from .masking import MaskingTask
28from .isrStatistics import IsrStatisticsTask
29from .isr import maskNans
30from .ptcDataset import PhotonTransferCurveDataset
31from .isrFunctions import isTrimmedExposure, compareCameraKeywords
32
33
34class IsrTaskLSSTConnections(pipeBase.PipelineTaskConnections,
35 dimensions={"instrument", "exposure", "detector"},
36 defaultTemplates={}):
37 ccdExposure = cT.Input(
38 name="raw",
39 doc="Input exposure to process.",
40 storageClass="Exposure",
41 dimensions=["instrument", "exposure", "detector"],
42 )
43 camera = cT.PrerequisiteInput(
44 name="camera",
45 storageClass="Camera",
46 doc="Input camera to construct complete exposures.",
47 dimensions=["instrument"],
48 isCalibration=True,
49 )
50 dnlLUT = cT.PrerequisiteInput(
51 name="dnlLUT",
52 doc="Look-up table for differential non-linearity.",
53 storageClass="IsrCalib",
54 dimensions=["instrument", "exposure", "detector"],
55 isCalibration=True,
56 # TODO DM 36636
57 )
58 bias = cT.PrerequisiteInput(
59 name="bias",
60 doc="Input bias calibration.",
61 storageClass="ExposureF",
62 dimensions=["instrument", "detector"],
63 isCalibration=True,
64 )
65 deferredChargeCalib = cT.PrerequisiteInput(
66 name="cti",
67 doc="Deferred charge/CTI correction dataset.",
68 storageClass="IsrCalib",
69 dimensions=["instrument", "detector"],
70 isCalibration=True,
71 )
72 linearizer = cT.PrerequisiteInput(
73 name='linearizer',
74 storageClass="Linearizer",
75 doc="Linearity correction calibration.",
76 dimensions=["instrument", "detector"],
77 isCalibration=True,
78 )
79 ptc = cT.PrerequisiteInput(
80 name="ptc",
81 doc="Input Photon Transfer Curve dataset",
82 storageClass="PhotonTransferCurveDataset",
83 dimensions=["instrument", "detector"],
84 isCalibration=True,
85 )
86 crosstalk = cT.PrerequisiteInput(
87 name="crosstalk",
88 doc="Input crosstalk object",
89 storageClass="CrosstalkCalib",
90 dimensions=["instrument", "detector"],
91 isCalibration=True,
92 )
93 defects = cT.PrerequisiteInput(
94 name='defects',
95 doc="Input defect tables.",
96 storageClass="Defects",
97 dimensions=["instrument", "detector"],
98 isCalibration=True,
99 )
100 bfKernel = cT.PrerequisiteInput(
101 name="bfk",
102 doc="Complete kernel + gain solutions.",
103 storageClass="BrighterFatterKernel",
104 dimensions=["instrument", "detector"],
105 isCalibration=True,
106 )
107 dark = cT.PrerequisiteInput(
108 name='dark',
109 doc="Input dark calibration.",
110 storageClass="ExposureF",
111 dimensions=["instrument", "detector"],
112 isCalibration=True,
113 )
114 flat = cT.PrerequisiteInput(
115 name="flat",
116 doc="Input flat calibration.",
117 storageClass="ExposureF",
118 dimensions=["instrument", "detector", "physical_filter"],
119 isCalibration=True,
120 )
121 outputExposure = cT.Output(
122 name='postISRCCD',
123 doc="Output ISR processed exposure.",
124 storageClass="Exposure",
125 dimensions=["instrument", "exposure", "detector"],
126 )
127 preInterpExposure = cT.Output(
128 name='preInterpISRCCD',
129 doc="Output ISR processed exposure, with pixels left uninterpolated.",
130 storageClass="ExposureF",
131 dimensions=["instrument", "exposure", "detector"],
132 )
133 outputBin1Exposure = cT.Output(
134 name="postIsrBin1",
135 doc="First binned image.",
136 storageClass="ExposureF",
137 dimensions=["instrument", "exposure", "detector"],
138 )
139 outputBin2Exposure = cT.Output(
140 name="postIsrBin2",
141 doc="Second binned image.",
142 storageClass="ExposureF",
143 dimensions=["instrument", "exposure", "detector"],
144 )
145
146 outputStatistics = cT.Output(
147 name="isrStatistics",
148 doc="Output of additional statistics table.",
149 storageClass="StructuredDataDict",
150 dimensions=["instrument", "exposure", "detector"],
151 )
152
153 def __init__(self, *, config=None):
154 super().__init__(config=config)
155
156 if config.doBootstrap:
157 del self.ptc
158 if config.doDiffNonLinearCorrection is not True:
159 del self.dnlLUT
160 if config.doBias is not True:
161 del self.bias
162 if config.doDeferredCharge is not True:
163 del self.deferredChargeCalib
164 if config.doLinearize is not True:
165 del self.linearizer
166 if not config.doCrosstalk:
167 del self.crosstalk
168 if config.doDefect is not True:
169 del self.defects
170 if config.doBrighterFatter is not True:
171 del self.bfKernel
172 if config.doDark is not True:
173 del self.dark
174 if config.doFlat is not True:
175 del self.flat
176
177 if config.doBinnedExposures is not True:
178 del self.outputBin1Exposure
179 del self.outputBin2Exposure
180 if config.doSaveInterpPixels is not True:
181 del self.preInterpExposure
182
183 if config.doCalculateStatistics is not True:
184 del self.outputStatistics
185
186
187class IsrTaskLSSTConfig(pipeBase.PipelineTaskConfig,
188 pipelineConnections=IsrTaskLSSTConnections):
189 """Configuration parameters for IsrTaskLSST.
190
191 Items are grouped in the order in which they are executed by the task.
192 """
193 expectWcs = pexConfig.Field(
194 dtype=bool,
195 default=True,
196 doc="Expect input science images to have a WCS (set False for e.g. spectrographs)."
197 )
198 qa = pexConfig.ConfigField(
199 dtype=isrQa.IsrQaConfig,
200 doc="QA related configuration options.",
201 )
202 doHeaderProvenance = pexConfig.Field(
203 dtype=bool,
204 default=True,
205 doc="Write calibration identifiers into output exposure header.",
206 )
207
208 # Calib checking configuration:
209 doRaiseOnCalibMismatch = pexConfig.Field(
210 dtype=bool,
211 default=False,
212 doc="Should IsrTaskLSST halt if exposure and calibration header values do not match?",
213 )
214 cameraKeywordsToCompare = pexConfig.ListField(
215 dtype=str,
216 doc="List of header keywords to compare between exposure and calibrations.",
217 default=[],
218 )
219
220 # Differential non-linearity correction.
221 doDiffNonLinearCorrection = pexConfig.Field(
222 dtype=bool,
223 doc="Do differential non-linearity correction?",
224 default=False,
225 )
226
227 doBootstrap = pexConfig.Field(
228 dtype=bool,
229 default=False,
230 doc="Is this task to be run in a ``bootstrap`` fashion that does not require "
231 "a PTC or full calibrations?",
232 )
233
234 overscanCamera = pexConfig.ConfigField(
235 dtype=OverscanCameraConfig,
236 doc="Per-detector and per-amplifier overscan configurations.",
237 )
238
239 # Amplifier to CCD assembly configuration.
240 doAssembleCcd = pexConfig.Field(
241 dtype=bool,
242 default=True,
243 doc="Assemble amp-level exposures into a ccd-level exposure?"
244 )
245 assembleCcd = pexConfig.ConfigurableField(
246 target=AssembleCcdTask,
247 doc="CCD assembly task.",
248 )
249
250 # Bias subtraction.
251 doBias = pexConfig.Field(
252 dtype=bool,
253 doc="Apply bias frame correction?",
254 default=True,
255 )
256
257 # Deferred charge correction.
258 doDeferredCharge = pexConfig.Field(
259 dtype=bool,
260 doc="Apply deferred charge correction?",
261 default=True,
262 )
263 deferredChargeCorrection = pexConfig.ConfigurableField(
264 target=DeferredChargeTask,
265 doc="Deferred charge correction task.",
266 )
267
268 # Linearization.
269 doLinearize = pexConfig.Field(
270 dtype=bool,
271 doc="Correct for nonlinearity of the detector's response?",
272 default=True,
273 )
274
275 # Gains.
276 doCorrectGains = pexConfig.Field(
277 dtype=bool,
278 doc="Apply temperature correction to the gains?",
279 default=False,
280 )
281 doApplyGains = pexConfig.Field(
282 dtype=bool,
283 doc="Apply gains to the image?",
284 default=True,
285 )
286
287 # Variance construction.
288 doVariance = pexConfig.Field(
289 dtype=bool,
290 doc="Calculate variance?",
291 default=True
292 )
293 maskNegativeVariance = pexConfig.Field(
294 dtype=bool,
295 doc="Mask pixels that claim a negative variance. This likely indicates a failure "
296 "in the measurement of the overscan at an edge due to the data falling off faster "
297 "than the overscan model can account for it.",
298 default=True,
299 )
300 negativeVarianceMaskName = pexConfig.Field(
301 dtype=str,
302 doc="Mask plane to use to mark pixels with negative variance, if `maskNegativeVariance` is True.",
303 default="BAD",
304 )
305 doSaturation = pexConfig.Field(
306 dtype=bool,
307 doc="Mask saturated pixels? NB: this is totally independent of the"
308 " interpolation option - this is ONLY setting the bits in the mask."
309 " To have them interpolated make sure doInterpolate=True and"
310 " maskListToInterpolate includes SAT.",
311 default=True,
312 )
313 saturatedMaskName = pexConfig.Field(
314 dtype=str,
315 doc="Name of mask plane to use in saturation detection and interpolation.",
316 default="SAT",
317 )
318 defaultSaturationSource = pexConfig.ChoiceField(
319 dtype=str,
320 doc="Source to retrieve default amp-level saturation values.",
321 allowed={
322 "NONE": "No default saturation values; only config overrides will be used.",
323 "CAMERAMODEL": "Use the default from the camera model (old defaults).",
324 "PTCTURNOFF": "Use the ptcTurnoff value as the saturation level.",
325 },
326 default="PTCTURNOFF",
327 )
328 doSuspect = pexConfig.Field(
329 dtype=bool,
330 doc="Mask suspect pixels?",
331 default=True,
332 )
333 suspectMaskName = pexConfig.Field(
334 dtype=str,
335 doc="Name of mask plane to use for suspect pixels.",
336 default="SUSPECT",
337 )
338 defaultSuspectSource = pexConfig.ChoiceField(
339 dtype=str,
340 doc="Source to retrieve default amp-level suspect values.",
341 allowed={
342 "NONE": "No default suspect values; only config overrides will be used.",
343 "CAMERAMODEL": "Use the default from the camera model (old defaults).",
344 "PTCTURNOFF": "Use the ptcTurnoff value as the suspect level.",
345 },
346 default="PTCTURNOFF",
347 )
348
349 # Crosstalk.
350 doCrosstalk = pexConfig.Field(
351 dtype=bool,
352 doc="Apply intra-CCD crosstalk correction?",
353 default=True,
354 )
355 crosstalk = pexConfig.ConfigurableField(
356 target=CrosstalkTask,
357 doc="Intra-CCD crosstalk correction.",
358 )
359
360 # Masking options.
361 doITLDipMask = pexConfig.Field(
362 dtype=bool,
363 doc="Apply ``ITL dip`` masking. The ``itlDipMaskPlane`` mask plane "
364 "will be added even if this configuration is False.",
365 default=True,
366 )
367 itlDipMaskPlanes = pexConfig.ListField(
368 dtype=str,
369 doc="Mask plane to use for ITL dip pixels.",
370 default=["SUSPECT", "ITL_DIP"],
371 )
372 doDefect = pexConfig.Field(
373 dtype=bool,
374 doc="Apply correction for CCD defects, e.g. hot pixels?",
375 default=True,
376 )
377 doNanMasking = pexConfig.Field(
378 dtype=bool,
379 doc="Mask non-finite (NAN, inf) pixels. The UNMASKEDNAN mask plane "
380 "will be added even if this configuration is False.",
381 default=True,
382 )
383 doWidenSaturationTrails = pexConfig.Field(
384 dtype=bool,
385 doc="Widen bleed trails based on their width.",
386 default=False,
387 )
388 masking = pexConfig.ConfigurableField(
389 target=MaskingTask,
390 doc="Masking task."
391 )
392 doITLEdgeBleedMask = pexConfig.Field(
393 dtype=bool,
394 doc="Mask edge bleeds from saturated columns in ITL amplifiers.",
395 default=True,
396 )
397 doITLSatSagMask = pexConfig.Field(
398 dtype=bool,
399 doc="Mask columns presenting saturation sag.",
400 default=True,
401 )
402 itlEdgeBleedSatMinArea = pexConfig.Field(
403 dtype=int,
404 doc="Minimum limit for saturated cores footprint area.",
405 default=10000,
406 )
407 itlEdgeBleedSatMaxArea = pexConfig.Field(
408 dtype=int,
409 doc="Maximum limit for saturated cores footprint area.",
410 default=100000,
411 )
412 itlEdgeBleedThreshold = pexConfig.Field(
413 dtype=float,
414 doc="Sky background threshold for edge bleed detection.",
415 default=5000.,
416 )
417 itlEdgeBleedModelConstant = pexConfig.Field(
418 dtype=float,
419 doc="Constant in the edge bleed exponential decay model.",
420 default=0.02,
421 )
422
423 # Interpolation options.
424 doInterpolate = pexConfig.Field(
425 dtype=bool,
426 doc="Interpolate masked pixels?",
427 default=True,
428 )
429 maskListToInterpolate = pexConfig.ListField(
430 dtype=str,
431 doc="List of mask planes that should be interpolated.",
432 default=["SAT", "BAD", "UNMASKEDNAN"],
433 )
434 doSaveInterpPixels = pexConfig.Field(
435 dtype=bool,
436 doc="Save a copy of the pre-interpolated pixel values?",
437 default=False,
438 )
439 useLegacyInterp = pexConfig.Field(
440 dtype=bool,
441 doc="Use the legacy interpolation algorithm. If False use Gaussian Process.",
442 default=True,
443 )
444
445 # Amp offset correction.
446 doAmpOffset = pexConfig.Field(
447 doc="Calculate amp offset corrections?",
448 dtype=bool,
449 default=False,
450 )
451 ampOffset = pexConfig.ConfigurableField(
452 doc="Amp offset correction task.",
453 target=AmpOffsetTask,
454 )
455
456 # Initial masking options.
457 doSetBadRegions = pexConfig.Field(
458 dtype=bool,
459 doc="Should we set the level of all BAD patches of the chip to the chip's average value?",
460 default=True,
461 )
462
463 # Brighter-Fatter correction.
464 doBrighterFatter = pexConfig.Field(
465 dtype=bool,
466 doc="Apply the brighter-fatter correction?",
467 default=True,
468 )
469 brighterFatterLevel = pexConfig.ChoiceField(
470 dtype=str,
471 doc="The level at which to correct for brighter-fatter.",
472 allowed={
473 "AMP": "Every amplifier treated separately.",
474 "DETECTOR": "One kernel per detector.",
475 },
476 default="DETECTOR",
477 )
478 brighterFatterMaxIter = pexConfig.Field(
479 dtype=int,
480 doc="Maximum number of iterations for the brighter-fatter correction.",
481 default=10,
482 )
483 brighterFatterThreshold = pexConfig.Field(
484 dtype=float,
485 doc="Threshold used to stop iterating the brighter-fatter correction. It is the "
486 "absolute value of the difference between the current corrected image and the one "
487 "from the previous iteration summed over all the pixels.",
488 default=1000,
489 )
490 brighterFatterMaskListToInterpolate = pexConfig.ListField(
491 dtype=str,
492 doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
493 "correction.",
494 default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
495 )
496 brighterFatterMaskGrowSize = pexConfig.Field(
497 dtype=int,
498 doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
499 "when brighter-fatter correction is applied.",
500 default=2,
501 )
502 brighterFatterFwhmForInterpolation = pexConfig.Field(
503 dtype=float,
504 doc="FWHM of PSF in arcseconds used for interpolation in brighter-fatter correction "
505 "(currently unused).",
506 default=1.0,
507 )
508 growSaturationFootprintSize = pexConfig.Field(
509 dtype=int,
510 doc="Number of pixels by which to grow the saturation footprints.",
511 default=1,
512 )
513
514 # Dark subtraction.
515 doDark = pexConfig.Field(
516 dtype=bool,
517 doc="Apply dark frame correction.",
518 default=True,
519 )
520
521 # Flat correction.
522 doFlat = pexConfig.Field(
523 dtype=bool,
524 doc="Apply flat field correction.",
525 default=True,
526 )
527 flatScalingType = pexConfig.ChoiceField(
528 dtype=str,
529 doc="The method for scaling the flat on the fly.",
530 default='USER',
531 allowed={
532 "USER": "Scale by flatUserScale",
533 "MEAN": "Scale by the inverse of the mean",
534 "MEDIAN": "Scale by the inverse of the median",
535 },
536 )
537 flatUserScale = pexConfig.Field(
538 dtype=float,
539 doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise.",
540 default=1.0,
541 )
542
543 # Calculate image quality statistics?
544 doStandardStatistics = pexConfig.Field(
545 dtype=bool,
546 doc="Should standard image quality statistics be calculated?",
547 default=True,
548 )
549 # Calculate additional statistics?
550 doCalculateStatistics = pexConfig.Field(
551 dtype=bool,
552 doc="Should additional ISR statistics be calculated?",
553 default=True,
554 )
555 isrStats = pexConfig.ConfigurableField(
556 target=IsrStatisticsTask,
557 doc="Task to calculate additional statistics.",
558 )
559
560 # Make binned images?
561 doBinnedExposures = pexConfig.Field(
562 dtype=bool,
563 doc="Should binned exposures be calculated?",
564 default=False,
565 )
566 binning = pexConfig.ConfigurableField(
567 target=BinExposureTask,
568 doc="Task to bin the exposure.",
569 )
570 binFactor1 = pexConfig.Field(
571 dtype=int,
572 doc="Binning factor for first binned exposure. This is intended for a finely binned output.",
573 default=8,
574 check=lambda x: x > 1,
575 )
576 binFactor2 = pexConfig.Field(
577 dtype=int,
578 doc="Binning factor for second binned exposure. This is intended for a coarsely binned output.",
579 default=64,
580 check=lambda x: x > 1,
581 )
582
583 def validate(self):
584 super().validate()
585
586 if self.doBootstrap:
587 # Additional checks in bootstrap (no PTC/gains) mode.
588 if self.doApplyGains:
589 raise ValueError("Cannot run task with doBootstrap=True and doApplyGains=True.")
590 if self.doCorrectGains:
591 raise ValueError("Cannot run task with doBootstrap=True and doCorrectGains=True.")
592 if self.doCrosstalk and self.crosstalk.doQuadraticCrosstalkCorrection:
593 raise ValueError("Cannot apply quadratic crosstalk correction with doBootstrap=True.")
594 if self.doITLEdgeBleedMask and not self.doSaturation:
595 raise ValueError("Cannot do edge bleed masking when doSaturation=False.")
596
597 def setDefaults(self):
598 super().setDefaults()
599
600
601class IsrTaskLSST(pipeBase.PipelineTask):
602 ConfigClass = IsrTaskLSSTConfig
603 _DefaultName = "isrLSST"
604
605 def __init__(self, **kwargs):
606 super().__init__(**kwargs)
607 self.makeSubtask("assembleCcd")
608 self.makeSubtask("deferredChargeCorrection")
609 self.makeSubtask("crosstalk")
610 self.makeSubtask("masking")
611 self.makeSubtask("isrStats")
612 self.makeSubtask("ampOffset")
613 self.makeSubtask("binning")
614
615 def runQuantum(self, butlerQC, inputRefs, outputRefs):
616
617 inputs = butlerQC.get(inputRefs)
618 self.validateInput(inputs)
619
620 if self.config.doHeaderProvenance:
621 # Add calibration provenanace info to header.
622 exposureMetadata = inputs['ccdExposure'].metadata
623 for inputName in sorted(list(inputs.keys())):
624 reference = getattr(inputRefs, inputName, None)
625 if reference is not None and hasattr(reference, "run"):
626 runKey = f"LSST CALIB RUN {inputName.upper()}"
627 runValue = reference.run
628 idKey = f"LSST CALIB UUID {inputName.upper()}"
629 idValue = str(reference.id)
630 dateKey = f"LSST CALIB DATE {inputName.upper()}"
631 dateValue = self.extractCalibDate(inputs[inputName])
632 if dateValue != "Unknown Unknown":
633 butlerQC.add_additional_provenance(reference, {"calib date": dateValue})
634
635 exposureMetadata[runKey] = runValue
636 exposureMetadata[idKey] = idValue
637 exposureMetadata[dateKey] = dateValue
638
639 outputs = self.run(**inputs)
640 butlerQC.put(outputs, outputRefs)
641
642 def validateInput(self, inputs):
643 """
644 This is a check that all the inputs required by the config
645 are available.
646 """
647
648 inputMap = {'dnlLUT': self.config.doDiffNonLinearCorrection,
649 'bias': self.config.doBias,
650 'deferredChargeCalib': self.config.doDeferredCharge,
651 'linearizer': self.config.doLinearize,
652 'ptc': self.config.doApplyGains,
653 'crosstalk': self.config.doCrosstalk,
654 'defects': self.config.doDefect,
655 'bfKernel': self.config.doBrighterFatter,
656 'dark': self.config.doDark,
657 }
658
659 for calibrationFile, configValue in inputMap.items():
660 if configValue and inputs[calibrationFile] is None:
661 raise RuntimeError("Must supply ", calibrationFile)
662
663 def diffNonLinearCorrection(self, ccdExposure, dnlLUT, **kwargs):
664 # TODO DM 36636
665 # isrFunctions.diffNonLinearCorrection
666 pass
667
668 def maskFullAmplifiers(self, ccdExposure, detector, defects, gains=None):
669 """
670 Check for fully masked bad amplifiers and mask them.
671
672 This includes defects which cover full amplifiers, as well
673 as amplifiers with nan gain values which should be used
674 if self.config.doApplyGains=True.
675
676 Full defect masking happens later to allow for defects which
677 cross amplifier boundaries.
678
679 Parameters
680 ----------
681 ccdExposure : `lsst.afw.image.Exposure`
682 Input exposure to be masked.
683 detector : `lsst.afw.cameraGeom.Detector`
684 Detector object.
685 defects : `lsst.ip.isr.Defects`
686 List of defects. Used to determine if an entire
687 amplifier is bad.
688 gains : `dict` [`str`, `float`], optional
689 Dictionary of gains to check if
690 self.config.doApplyGains=True.
691
692 Returns
693 -------
694 badAmpDict : `str`[`bool`]
695 Dictionary of amplifiers, keyed by name, value is True if
696 amplifier is fully masked.
697 """
698 badAmpDict = {}
699
700 maskedImage = ccdExposure.getMaskedImage()
701
702 for amp in detector:
703 ampName = amp.getName()
704 badAmpDict[ampName] = False
705
706 # Check if entire amp region is defined as a defect
707 # NB: need to use amp.getBBox() for correct comparison with current
708 # defects definition.
709 if defects is not None:
710 badAmpDict[ampName] = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
711
712 if badAmpDict[ampName]:
713 self.log.warning("Amplifier %s is bad (completely covered with defects)", ampName)
714
715 if gains is not None and self.config.doApplyGains:
716 if not math.isfinite(gains[ampName]):
717 badAmpDict[ampName] = True
718
719 self.log.warning("Amplifier %s is bad (non-finite gain)", ampName)
720
721 # In the case of a bad amp, we will set mask to "BAD"
722 # (here use amp.getRawBBox() for correct association with pixels in
723 # current ccdExposure).
724 if badAmpDict[ampName]:
725 dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
726 afwImage.PARENT)
727 maskView = dataView.getMask()
728 maskView |= maskView.getPlaneBitMask("BAD")
729 del maskView
730
731 return badAmpDict
732
733 def maskSaturatedPixels(self, badAmpDict, ccdExposure, detector, detectorConfig, ptc=None):
734 """
735 Mask SATURATED and SUSPECT pixels and check if any amplifiers
736 are fully masked.
737
738 Parameters
739 ----------
740 badAmpDict : `str` [`bool`]
741 Dictionary of amplifiers, keyed by name, value is True if
742 amplifier is fully masked.
743 ccdExposure : `lsst.afw.image.Exposure`
744 Input exposure to be masked.
745 detector : `lsst.afw.cameraGeom.Detector`
746 Detector object.
747 defects : `lsst.ip.isr.Defects`
748 List of defects. Used to determine if an entire
749 amplifier is bad.
750 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
751 Per-amplifier configurations.
752 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
753 PTC dataset (used if configured to use PTCTURNOFF).
754
755 Returns
756 -------
757 badAmpDict : `str`[`bool`]
758 Dictionary of amplifiers, keyed by name.
759 """
760 maskedImage = ccdExposure.getMaskedImage()
761
762 metadata = ccdExposure.metadata
763
764 if self.config.doSaturation and self.config.defaultSaturationSource == "PTCTURNOFF" and ptc is None:
765 raise RuntimeError("Must provide ptc if using PTCTURNOFF as saturation source.")
766 if self.config.doSuspect and self.config.defaultSuspectSource == "PTCTURNOFF" and ptc is None:
767 raise RuntimeError("Must provide ptc if using PTCTURNOFF as suspect source.")
768
769 for amp in detector:
770 ampName = amp.getName()
771
772 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
773
774 if badAmpDict[ampName]:
775 # No need to check fully bad amplifiers.
776 continue
777
778 # Mask saturated and suspect pixels.
779 limits = {}
780 if self.config.doSaturation:
781 if self.config.defaultSaturationSource == "PTCTURNOFF":
782 limits.update({self.config.saturatedMaskName: ptc.ptcTurnoff[amp.getName()]})
783 elif self.config.defaultSaturationSource == "CAMERAMODEL":
784 # Set to the default from the camera model.
785 limits.update({self.config.saturatedMaskName: amp.getSaturation()})
786 elif self.config.defaultSaturationSource == "NONE":
787 limits.update({self.config.saturatedMaskName: numpy.inf})
788
789 # And update if it is set in the config.
790 if math.isfinite(ampConfig.saturation):
791 limits.update({self.config.saturatedMaskName: ampConfig.saturation})
792 metadata[f"LSST ISR SATURATION LEVEL {ampName}"] = limits[self.config.saturatedMaskName]
793
794 if self.config.doSuspect:
795 if self.config.defaultSuspectSource == "PTCTURNOFF":
796 limits.update({self.config.suspectMaskName: ptc.ptcTurnoff[amp.getName()]})
797 elif self.config.defaultSuspectSource == "CAMERAMODEL":
798 # Set to the default from the camera model.
799 limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
800 elif self.config.defaultSuspectSource == "NONE":
801 limits.update({self.config.suspectMaskName: numpy.inf})
802
803 # And update if it set in the config.
804 if math.isfinite(ampConfig.suspectLevel):
805 limits.update({self.config.suspectMaskName: ampConfig.suspectLevel})
806 metadata[f"LSST ISR SUSPECT LEVEL {ampName}"] = limits[self.config.suspectMaskName]
807
808 for maskName, maskThreshold in limits.items():
809 if not math.isnan(maskThreshold):
810 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
811 toMask = (dataView.image.array >= maskThreshold)
812 dataView.mask.array[toMask] |= dataView.mask.getPlaneBitMask(maskName)
813
814 # Determine if we've fully masked this amplifier with SUSPECT and
815 # SAT pixels.
816 maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
817 afwImage.PARENT)
818 maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
819 self.config.suspectMaskName])
820 if numpy.all(maskView.getArray() & maskVal > 0):
821 self.log.warning("Amplifier %s is bad (completely SATURATED or SUSPECT)", ampName)
822 badAmpDict[ampName] = True
823 maskView |= maskView.getPlaneBitMask("BAD")
824
825 return badAmpDict
826
827 def maskITLSatEdgesAndColumns(self, exposure, badAmpDict):
828
829 # The following steps will rely on the footprint of saturated
830 # cores with large areas.
831 thresh = afwDetection.Threshold(exposure.mask.getPlaneBitMask("SAT"),
832 afwDetection.Threshold.BITMASK
833 )
834 fpList = afwDetection.FootprintSet(exposure.mask, thresh).getFootprints()
835
836 satAreas = numpy.asarray([fp.getArea() for fp in fpList])
837 largeAreas, = numpy.where((satAreas >= self.config.itlEdgeBleedSatMinArea)
838 & (satAreas < self.config.itlEdgeBleedSatMaxArea))
839
840 for largeAreasIndex in largeAreas:
841
842 fpCore = fpList[largeAreasIndex]
843
844 # Edge bleed masking
845 if self.config.doITLEdgeBleedMask:
846 isrFunctions.maskITLEdgeBleed(ccdExposure=exposure,
847 badAmpDict=badAmpDict,
848 fpCore=fpCore,
849 itlEdgeBleedSatMinArea=self.config.itlEdgeBleedSatMinArea,
850 itlEdgeBleedSatMaxArea=self.config.itlEdgeBleedSatMaxArea,
851 itlEdgeBleedThreshold=self.config.itlEdgeBleedThreshold,
852 itlEdgeBleedModelConstant=self.config.itlEdgeBleedModelConstant,
853 saturatedMaskName=self.config.saturatedMaskName,
854 log=self.log
855 )
856 if self.config.doITLSatSagMask:
857 isrFunctions.maskITLSatSag(ccdExposure=exposure, fpCore=fpCore,
858 saturatedMaskName=self.config.saturatedMaskName)
859
860 def overscanCorrection(self, mode, detectorConfig, detector, badAmpDict, ccdExposure):
861 """Apply serial overscan correction in place to all amps.
862
863 The actual overscan subtraction is performed by the
864 `lsst.ip.isr.overscan.OverscanTask`, which is called here.
865
866 Parameters
867 ----------
868 mode : `str`
869 Must be `SERIAL` or `PARALLEL`.
870 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
871 Per-amplifier configurations.
872 detector : `lsst.afw.cameraGeom.Detector`
873 Detector object.
874 badAmpDict : `dict`
875 Dictionary of amp name to whether it is a bad amp.
876 ccdExposure : `lsst.afw.image.Exposure`
877 Exposure to have overscan correction performed.
878
879 Returns
880 -------
881 overscans : `list` [`lsst.pipe.base.Struct` or None]
882 Overscan measurements (always in adu).
883 Each result struct has components:
884
885 ``imageFit``
886 Value or fit subtracted from the amplifier image data.
887 (scalar or `lsst.afw.image.Image`)
888 ``overscanFit``
889 Value or fit subtracted from the overscan image data.
890 (scalar or `lsst.afw.image.Image`)
891 ``overscanImage``
892 Image of the overscan region with the overscan
893 correction applied. This quantity is used to estimate
894 the amplifier read noise empirically.
895 (`lsst.afw.image.Image`)
896 ``overscanMean``
897 Mean overscan fit value. (`float`)
898 ``overscanMedian``
899 Median overscan fit value. (`float`)
900 ``overscanSigma``
901 Clipped standard deviation of the overscan fit. (`float`)
902 ``residualMean``
903 Mean of the overscan after fit subtraction. (`float`)
904 ``residualMedian``
905 Median of the overscan after fit subtraction. (`float`)
906 ``residualSigma``
907 Clipped standard deviation of the overscan after fit
908 subtraction. (`float`)
909
910 See Also
911 --------
912 lsst.ip.isr.overscan.OverscanTask
913 """
914 if mode not in ["SERIAL", "PARALLEL"]:
915 raise ValueError("Mode must be SERIAL or PARALLEL")
916
917 # This returns a list in amp order, with None for uncorrected amps.
918 overscans = []
919
920 for i, amp in enumerate(detector):
921 ampName = amp.getName()
922
923 ampConfig = detectorConfig.getOverscanAmpConfig(amp)
924
925 if mode == "SERIAL" and not ampConfig.doSerialOverscan:
926 self.log.debug(
927 "ISR_OSCAN: Amplifier %s/%s configured to skip serial overscan.",
928 detector.getName(),
929 ampName,
930 )
931 results = None
932 elif mode == "PARALLEL" and not ampConfig.doParallelOverscan:
933 self.log.debug(
934 "ISR_OSCAN: Amplifier %s configured to skip parallel overscan.",
935 detector.getName(),
936 ampName,
937 )
938 results = None
939 elif badAmpDict[ampName] or not ccdExposure.getBBox().contains(amp.getBBox()):
940 results = None
941 else:
942 # This check is to confirm that we are not trying to run
943 # overscan on an already trimmed image.
944 if isTrimmedExposure(ccdExposure):
945 self.log.warning(
946 "ISR_OSCAN: No overscan region for amp %s. Not performing overscan correction.",
947 ampName,
948 )
949 results = None
950 else:
951 if mode == "SERIAL":
952 # We need to set up the subtask here with a custom
953 # configuration.
954 serialOverscan = SerialOverscanCorrectionTask(config=ampConfig.serialOverscanConfig)
955 results = serialOverscan.run(ccdExposure, amp)
956 else:
957 config = ampConfig.parallelOverscanConfig
958 parallelOverscan = ParallelOverscanCorrectionTask(
959 config=config,
960 )
961
962 metadata = ccdExposure.metadata
963
964 # We need to know the saturation level that was used
965 # for the parallel overscan masking. If it isn't set
966 # then the configured parallelOverscanSaturationLevel
967 # will be used instead (assuming
968 # doParallelOverscanSaturation is True). Note that
969 # this will have the correct units (adu or electron)
970 # depending on whether the gain has been applied.
971 if self.config.doSaturation:
972 saturationLevel = metadata[f"LSST ISR SATURATION LEVEL {amp.getName()}"]
973 saturationLevel *= config.parallelOverscanSaturationLevelAdjustmentFactor
974 else:
975 saturationLevel = config.parallelOverscanSaturationLevel
976 if ccdExposure.metadata["LSST ISR UNITS"] == "electron":
977 # Need to convert to electron from adu.
978 saturationLevel *= metadata[f"LSST ISR GAIN {amp.getName()}"]
979
980 self.log.debug(
981 "Using saturation level of %.2f for parallel overscan amp %s",
982 saturationLevel,
983 amp.getName(),
984 )
985
986 parallelOverscan.maskParallelOverscanAmp(
987 ccdExposure,
988 amp,
989 saturationLevel=saturationLevel,
990 )
991
992 results = parallelOverscan.run(ccdExposure, amp)
993
994 metadata = ccdExposure.metadata
995 keyBase = "LSST ISR OVERSCAN"
996
997 # The overscan is always in adu for the serial mode,
998 # but, it may be electron in the parallel mode if
999 # doApplyGains==True. If doApplyGains==True, then the
1000 # gains are applied to the untrimmed image, so the
1001 # overscan statistics units here will always match the
1002 # units of the image at this point.
1003 metadata[f"{keyBase} {mode} UNITS"] = ccdExposure.metadata["LSST ISR UNITS"]
1004 metadata[f"{keyBase} {mode} MEAN {ampName}"] = results.overscanMean
1005 metadata[f"{keyBase} {mode} MEDIAN {ampName}"] = results.overscanMedian
1006 metadata[f"{keyBase} {mode} STDEV {ampName}"] = results.overscanSigma
1007
1008 metadata[f"{keyBase} RESIDUAL {mode} MEAN {ampName}"] = results.residualMean
1009 metadata[f"{keyBase} RESIDUAL {mode} MEDIAN {ampName}"] = results.residualMedian
1010 metadata[f"{keyBase} RESIDUAL {mode} STDEV {ampName}"] = results.residualSigma
1011
1012 overscans.append(results)
1013
1014 # Question: should this be finer grained?
1015 ccdExposure.metadata.set("OVERSCAN", "Overscan corrected")
1016
1017 return overscans
1018
1019 def correctGains(self, exposure, ptc, gains):
1020 # TODO DM 36639
1021 gains = []
1022 readNoise = []
1023
1024 return gains, readNoise
1025
1026 def maskNegativeVariance(self, exposure):
1027 """Identify and mask pixels with negative variance values.
1028
1029 Parameters
1030 ----------
1031 exposure : `lsst.afw.image.Exposure`
1032 Exposure to process.
1033
1034 See Also
1035 --------
1036 lsst.ip.isr.isrFunctions.updateVariance
1037 """
1038 maskPlane = exposure.getMask().getPlaneBitMask(self.config.negativeVarianceMaskName)
1039 bad = numpy.where(exposure.getVariance().getArray() <= 0.0)
1040 exposure.mask.array[bad] |= maskPlane
1041
1042 def addVariancePlane(self, exposure, detector):
1043 """Add the variance plane to the image.
1044
1045 The gain and read noise per amp must have been set in the
1046 exposure metadata as ``LSST ISR GAIN ampName`` and
1047 ``LSST ISR READNOISE ampName`` with the units of the image.
1048 Unit conversions for the variance plane will be done as
1049 necessary based on the exposure units.
1050
1051 The units of the variance plane will always be of the same
1052 type as the units of the input image itself
1053 (``LSST ISR UNITS``^2).
1054
1055 Parameters
1056 ----------
1057 exposure : `lsst.afw.image.Exposure`
1058 The exposure to add the variance plane.
1059 detector : `lsst.afw.cameraGeom.Detector`
1060 Detector with geometry info.
1061 """
1062 # NOTE: this will fail if the exposure is not trimmed.
1063 if not isTrimmedExposure(exposure):
1064 raise RuntimeError("Exposure must be trimmed to add variance plane.")
1065
1066 isElectrons = (exposure.metadata["LSST ISR UNITS"] == "electron")
1067
1068 for amp in detector:
1069 if exposure.getBBox().contains(amp.getBBox()):
1070 self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1071 ampExposure = exposure.Factory(exposure, amp.getBBox())
1072
1073 # The effective gain is 1.0 if we are in electron units.
1074 # The metadata read noise is in the same units as the image.
1075 gain = exposure.metadata[f"LSST ISR GAIN {amp.getName()}"] if not isElectrons else 1.0
1076 readNoise = exposure.metadata[f"LSST ISR READNOISE {amp.getName()}"]
1077
1078 isrFunctions.updateVariance(
1079 maskedImage=ampExposure.maskedImage,
1080 gain=gain,
1081 readNoise=readNoise,
1082 replace=False,
1083 )
1084
1085 if self.config.qa is not None and self.config.qa.saveStats is True:
1086 qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1087 afwMath.MEDIAN | afwMath.STDEVCLIP)
1088 self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1089 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1090 qaStats.getValue(afwMath.STDEVCLIP))
1091
1092 if self.config.maskNegativeVariance:
1093 self.maskNegativeVariance(exposure)
1094
1095 def maskDefects(self, exposure, defectBaseList):
1096 """Mask defects using mask plane "BAD", in place.
1097
1098 Parameters
1099 ----------
1100 exposure : `lsst.afw.image.Exposure`
1101 Exposure to process.
1102
1103 defectBaseList : defect-type
1104 List of defects to mask. Can be of type `lsst.ip.isr.Defects`
1105 or `list` of `lsst.afw.image.DefectBase`.
1106 """
1107 maskedImage = exposure.getMaskedImage()
1108 if not isinstance(defectBaseList, Defects):
1109 # Promotes DefectBase to Defect
1110 defectList = Defects(defectBaseList)
1111 else:
1112 defectList = defectBaseList
1113 defectList.maskPixels(maskedImage, maskName="BAD")
1114
1115 def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
1116 """Mask edge pixels with applicable mask plane.
1117
1118 Parameters
1119 ----------
1120 exposure : `lsst.afw.image.Exposure`
1121 Exposure to process.
1122 numEdgePixels : `int`, optional
1123 Number of edge pixels to mask.
1124 maskPlane : `str`, optional
1125 Mask plane name to use.
1126 level : `str`, optional
1127 Level at which to mask edges.
1128 """
1129 maskedImage = exposure.getMaskedImage()
1130 maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
1131
1132 if numEdgePixels > 0:
1133 if level == 'DETECTOR':
1134 boxes = [maskedImage.getBBox()]
1135 elif level == 'AMP':
1136 boxes = [amp.getBBox() for amp in exposure.getDetector()]
1137
1138 for box in boxes:
1139 # This makes a bbox numEdgeSuspect pixels smaller than the
1140 # image on each side
1141 subImage = maskedImage[box]
1142 box.grow(-numEdgePixels)
1143 # Mask pixels outside box
1144 SourceDetectionTask.setEdgeBits(
1145 subImage,
1146 box,
1147 maskBitMask)
1148
1149 def maskNan(self, exposure):
1150 """Mask NaNs using mask plane "UNMASKEDNAN", in place.
1151
1152 Parameters
1153 ----------
1154 exposure : `lsst.afw.image.Exposure`
1155 Exposure to process.
1156
1157 Notes
1158 -----
1159 We mask over all non-finite values (NaN, inf), including those
1160 that are masked with other bits (because those may or may not be
1161 interpolated over later, and we want to remove all NaN/infs).
1162 Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
1163 preserve the historical name.
1164 """
1165 maskedImage = exposure.getMaskedImage()
1166
1167 # Find and mask NaNs
1168 maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
1169 numNans = maskNans(maskedImage, maskVal)
1170 self.metadata["NUMNANS"] = numNans
1171 if numNans > 0:
1172 self.log.warning("There were %d unmasked NaNs.", numNans)
1173
1174 def setBadRegions(self, exposure):
1175 """Set bad regions from large contiguous regions.
1176
1177 Parameters
1178 ----------
1179 exposure : `lsst.afw.Exposure`
1180 Exposure to set bad regions.
1181
1182 Notes
1183 -----
1184 Reset and interpolate bad pixels.
1185
1186 Large contiguous bad regions (which should have the BAD mask
1187 bit set) should have their values set to the image median.
1188 This group should include defects and bad amplifiers. As the
1189 area covered by these defects are large, there's little
1190 reason to expect that interpolation would provide a more
1191 useful value.
1192
1193 Smaller defects can be safely interpolated after the larger
1194 regions have had their pixel values reset. This ensures
1195 that the remaining defects adjacent to bad amplifiers (as an
1196 example) do not attempt to interpolate extreme values.
1197 """
1198 badPixelCount, badPixelValue = isrFunctions.setBadRegions(exposure)
1199 if badPixelCount > 0:
1200 self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1201
1202 @contextmanager
1203 def flatContext(self, exp, flat, dark=None):
1204 """Context manager that applies and removes flats and darks,
1205 if the task is configured to apply them.
1206
1207 Parameters
1208 ----------
1209 exp : `lsst.afw.image.Exposure`
1210 Exposure to process.
1211 flat : `lsst.afw.image.Exposure`
1212 Flat exposure the same size as ``exp``.
1213 dark : `lsst.afw.image.Exposure`, optional
1214 Dark exposure the same size as ``exp``.
1215
1216 Yields
1217 ------
1218 exp : `lsst.afw.image.Exposure`
1219 The flat and dark corrected exposure.
1220 """
1221 if self.config.doDark and dark is not None:
1222 self.darkCorrection(exp, dark)
1223 if self.config.doFlat and flat is not None:
1224 self.flatCorrection(exp, flat)
1225 try:
1226 yield exp
1227 finally:
1228 if self.config.doFlat and flat is not None:
1229 self.flatCorrection(exp, flat, invert=True)
1230 if self.config.doDark and dark is not None:
1231 self.darkCorrection(exp, dark, invert=True)
1232
1233 def getBrighterFatterKernel(self, detector, bfKernel):
1234 detName = detector.getName()
1235
1236 # This is expected to be a dictionary of amp-wise gains.
1237 bfGains = bfKernel.gain
1238 if bfKernel.level == 'DETECTOR':
1239 if detName in bfKernel.detKernels:
1240 bfKernelOut = bfKernel.detKernels[detName]
1241 return bfKernelOut, bfGains
1242 else:
1243 raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1244 elif bfKernel.level == 'AMP':
1245 self.log.info("Making DETECTOR level kernel from AMP based brighter "
1246 "fatter kernels.")
1247 bfKernel.makeDetectorKernelFromAmpwiseKernels(detName)
1248 bfKernelOut = bfKernel.detKernels[detName]
1249 return bfKernelOut, bfGains
1250
1251 def applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, brighterFatterApplyGain,
1252 bfGains):
1253 """Apply a brighter fatter correction to the image using the
1254 method defined in Coulton et al. 2019.
1255
1256 Note that this correction requires that the image is in units
1257 electrons.
1258
1259 Parameters
1260 ----------
1261 ccdExposure : `lsst.afw.image.Exposure`
1262 Exposure to process.
1263 flat : `lsst.afw.image.Exposure`
1264 Flat exposure the same size as ``exp``.
1265 dark : `lsst.afw.image.Exposure`, optional
1266 Dark exposure the same size as ``exp``.
1267 bfKernel : `lsst.ip.isr.BrighterFatterKernel`
1268 The brighter-fatter kernel.
1269 brighterFatterApplyGain : `bool`
1270 Apply the gain to convert the image to electrons?
1271 bfGains : `dict`
1272 The gains to use if brighterFatterApplyGain = True.
1273
1274 Yields
1275 ------
1276 exp : `lsst.afw.image.Exposure`
1277 The flat and dark corrected exposure.
1278 """
1279 interpExp = ccdExposure.clone()
1280
1281 # We need to interpolate before we do B-F. Note that
1282 # brighterFatterFwhmForInterpolation is currently unused.
1283 isrFunctions.interpolateFromMask(
1284 maskedImage=interpExp.getMaskedImage(),
1285 fwhm=self.config.brighterFatterFwhmForInterpolation,
1286 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1287 maskNameList=list(self.config.brighterFatterMaskListToInterpolate),
1288 useLegacyInterp=self.config.useLegacyInterp,
1289 )
1290 bfExp = interpExp.clone()
1291 bfResults = isrFunctions.brighterFatterCorrection(bfExp, bfKernel,
1292 self.config.brighterFatterMaxIter,
1293 self.config.brighterFatterThreshold,
1294 brighterFatterApplyGain,
1295 bfGains)
1296 bfCorrIters = bfResults[1]
1297 if bfCorrIters == self.config.brighterFatterMaxIter:
1298 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1299 bfResults[0])
1300 else:
1301 self.log.info("Finished brighter-fatter correction in %d iterations.",
1302 bfResults[1])
1303
1304 image = ccdExposure.getMaskedImage().getImage()
1305 bfCorr = bfExp.getMaskedImage().getImage()
1306 bfCorr -= interpExp.getMaskedImage().getImage()
1307 image += bfCorr
1308
1309 # Applying the brighter-fatter correction applies a
1310 # convolution to the science image. At the edges this
1311 # convolution may not have sufficient valid pixels to
1312 # produce a valid correction. Mark pixels within the size
1313 # of the brighter-fatter kernel as EDGE to warn of this
1314 # fact.
1315 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1316 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1317 maskPlane="EDGE")
1318
1319 if self.config.brighterFatterMaskGrowSize > 0:
1320 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1321 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1322 isrFunctions.growMasks(ccdExposure.getMask(),
1323 radius=self.config.brighterFatterMaskGrowSize,
1324 maskNameList=maskPlane,
1325 maskValue=maskPlane)
1326
1327 return ccdExposure, bfCorrIters
1328
1329 def darkCorrection(self, exposure, darkExposure, invert=False):
1330 """Apply dark correction in place.
1331
1332 Parameters
1333 ----------
1334 exposure : `lsst.afw.image.Exposure`
1335 Exposure to process.
1336 darkExposure : `lsst.afw.image.Exposure`
1337 Dark exposure of the same size as ``exposure``.
1338 invert : `Bool`, optional
1339 If True, re-add the dark to an already corrected image.
1340
1341 Raises
1342 ------
1343 RuntimeError
1344 Raised if either ``exposure`` or ``darkExposure`` do not
1345 have their dark time defined.
1346
1347 See Also
1348 --------
1349 lsst.ip.isr.isrFunctions.darkCorrection
1350 """
1351 expScale = exposure.visitInfo.darkTime
1352 if math.isnan(expScale):
1353 raise RuntimeError("Exposure darktime is NAN.")
1354 if darkExposure.visitInfo is not None \
1355 and not math.isnan(darkExposure.visitInfo.darkTime):
1356 darkScale = darkExposure.visitInfo.darkTime
1357 else:
1358 # DM-17444: darkExposure.visitInfo is None
1359 # so darkTime does not exist.
1360 self.log.warning("darkExposure.visitInfo does not exist. Using darkScale = 1.0.")
1361 darkScale = 1.0
1362
1363 isrFunctions.darkCorrection(
1364 maskedImage=exposure.maskedImage,
1365 darkMaskedImage=darkExposure.maskedImage,
1366 expScale=expScale,
1367 darkScale=darkScale,
1368 invert=invert,
1369 )
1370
1371 @staticmethod
1373 """Extract common calibration metadata values that will be written to
1374 output header.
1375
1376 Parameters
1377 ----------
1378 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1379 Calibration to pull date information from.
1380
1381 Returns
1382 -------
1383 dateString : `str`
1384 Calibration creation date string to add to header.
1385 """
1386 if hasattr(calib, "getMetadata"):
1387 if 'CALIB_CREATION_DATE' in calib.metadata:
1388 return " ".join((calib.metadata.get("CALIB_CREATION_DATE", "Unknown"),
1389 calib.metadata.get("CALIB_CREATION_TIME", "Unknown")))
1390 else:
1391 return " ".join((calib.metadata.get("CALIB_CREATE_DATE", "Unknown"),
1392 calib.metadata.get("CALIB_CREATE_TIME", "Unknown")))
1393 else:
1394 return "Unknown Unknown"
1395
1396 def compareUnits(self, calibMetadata, calibName):
1397 """Compare units from calibration to ISR units.
1398
1399 This compares calibration units (adu or electron) to whether
1400 doApplyGain is set.
1401
1402 Parameters
1403 ----------
1404 calibMetadata : `lsst.daf.base.PropertyList`
1405 Calibration metadata from header.
1406 calibName : `str`
1407 Calibration name for log message.
1408 """
1409 calibUnits = calibMetadata.get("LSST ISR UNITS", "adu")
1410 isrUnits = "electron" if self.config.doApplyGains else "adu"
1411 if calibUnits != isrUnits:
1412 if self.config.doRaiseOnCalibMismatch:
1413 raise RuntimeError(
1414 "Unit mismatch: isr has %s units but %s has %s units",
1415 isrUnits,
1416 calibName,
1417 calibUnits,
1418 )
1419 else:
1420 self.log.warning(
1421 "Unit mismatch: isr has %s units but %s has %s units",
1422 isrUnits,
1423 calibName,
1424 calibUnits,
1425 )
1426
1427 def convertIntToFloat(self, exposure):
1428 """Convert exposure image from uint16 to float.
1429
1430 If the exposure does not need to be converted, the input is
1431 immediately returned. For exposures that are converted to use
1432 floating point pixels, the variance is set to unity and the
1433 mask to zero.
1434
1435 Parameters
1436 ----------
1437 exposure : `lsst.afw.image.Exposure`
1438 The raw exposure to be converted.
1439
1440 Returns
1441 -------
1442 newexposure : `lsst.afw.image.Exposure`
1443 The input ``exposure``, converted to floating point pixels.
1444
1445 Raises
1446 ------
1447 RuntimeError
1448 Raised if the exposure type cannot be converted to float.
1449
1450 """
1451 if isinstance(exposure, afwImage.ExposureF):
1452 # Nothing to be done
1453 self.log.debug("Exposure already of type float.")
1454 return exposure
1455 if not hasattr(exposure, "convertF"):
1456 raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
1457
1458 newexposure = exposure.convertF()
1459 newexposure.variance[:] = 1
1460 newexposure.mask[:] = 0x0
1461
1462 return newexposure
1463
1464 def ditherCounts(self, exposure, detectorConfig, fallbackSeed=12345):
1465 """Dither the counts in the exposure.
1466
1467 Parameters
1468 ----------
1469 exposure : `lsst.afw.image.Exposure`
1470 The raw exposure to be dithered.
1471 detectorConfig : `lsst.ip.isr.OverscanDetectorConfig`
1472 Configuration for overscan/etc for this detector.
1473 fallbackSeed : `int`, optional
1474 Random seed to fall back to if exposure.getInfo().getId() is
1475 not set.
1476 """
1477 if detectorConfig.integerDitherMode == "NONE":
1478 # Nothing to do here.
1479 return
1480
1481 # This ID is a unique combination of {exposure, detector} for a raw
1482 # image as we have here. We additionally need to take the lower
1483 # 32 bits to be used as a random seed.
1484 seed = exposure.info.id & 0xFFFFFFFF
1485 if seed is None:
1486 seed = fallbackSeed
1487 self.log.warning("No exposure ID found; using fallback random seed.")
1488
1489 self.log.info("Seeding dithering random number generator with %d.", seed)
1490 rng = numpy.random.RandomState(seed=seed)
1491
1492 if detectorConfig.integerDitherMode == "POSITIVE":
1493 low = 0.0
1494 high = 1.0
1495 elif detectorConfig.integerDitherMode == "NEGATIVE":
1496 low = -1.0
1497 high = 0.0
1498 elif detectorConfig.integerDitherMode == "SYMMETRIC":
1499 low = -0.5
1500 high = 0.5
1501 else:
1502 raise RuntimeError("Invalid config")
1503
1504 exposure.image.array[:, :] += rng.uniform(low=low, high=high, size=exposure.image.array.shape)
1505
1506 @deprecated(
1507 reason=(
1508 "makeBinnedImages is no longer used. "
1509 "Please subtask lsst.ip.isr.BinExposureTask instead."
1510 ),
1511 version="v28", category=FutureWarning
1512 )
1513 def makeBinnedImages(self, exposure):
1514 """Make visualizeVisit style binned exposures.
1515
1516 Parameters
1517 ----------
1518 exposure : `lsst.afw.image.Exposure`
1519 Exposure to bin.
1520
1521 Returns
1522 -------
1523 bin1 : `lsst.afw.image.Exposure`
1524 Binned exposure using binFactor1.
1525 bin2 : `lsst.afw.image.Exposure`
1526 Binned exposure using binFactor2.
1527 """
1528 mi = exposure.getMaskedImage()
1529
1530 bin1 = afwMath.binImage(mi, self.config.binFactor1)
1531 bin2 = afwMath.binImage(mi, self.config.binFactor2)
1532
1533 bin1 = afwImage.makeExposure(bin1)
1534 bin2 = afwImage.makeExposure(bin2)
1535
1536 bin1.setInfo(exposure.getInfo())
1537 bin2.setInfo(exposure.getInfo())
1538
1539 return bin1, bin2
1540
1541 def run(self, ccdExposure, *, dnlLUT=None, bias=None, deferredChargeCalib=None, linearizer=None,
1542 ptc=None, crosstalk=None, defects=None, bfKernel=None, bfGains=None, dark=None,
1543 flat=None, camera=None, **kwargs
1544 ):
1545
1546 detector = ccdExposure.getDetector()
1547
1548 overscanDetectorConfig = self.config.overscanCamera.getOverscanDetectorConfig(detector)
1549
1550 if self.config.doBootstrap and ptc is not None:
1551 self.log.warning("Task configured with doBootstrap=True. Ignoring provided PTC.")
1552 ptc = None
1553
1554 if not self.config.doBootstrap:
1555 if ptc is None:
1556 raise RuntimeError("A PTC must be supplied if config.doBootstrap is False.")
1557
1558 # Validation step: check inputs match exposure configuration.
1559 exposureMetadata = ccdExposure.metadata
1560 doRaise = self.config.doRaiseOnCalibMismatch
1561 keywords = self.config.cameraKeywordsToCompare
1562 if not self.config.doBootstrap:
1563 compareCameraKeywords(doRaise, keywords, exposureMetadata, ptc, "PTC", log=self.log)
1564 else:
1565 if self.config.doCorrectGains:
1566 raise RuntimeError("doCorrectGains is True but no ptc provided.")
1567 if self.config.doDiffNonLinearCorrection:
1568 if dnlLUT is None:
1569 raise RuntimeError("doDiffNonLinearCorrection is True but no dnlLUT provided.")
1570 compareCameraKeywords(doRaise, keywords, exposureMetadata, dnlLUT, "dnlLUT", log=self.log)
1571 if self.config.doLinearize:
1572 if linearizer is None:
1573 raise RuntimeError("doLinearize is True but no linearizer provided.")
1574 compareCameraKeywords(doRaise, keywords, exposureMetadata, linearizer, "linearizer", log=self.log)
1575 if self.config.doBias:
1576 if bias is None:
1577 raise RuntimeError("doBias is True but no bias provided.")
1578 compareCameraKeywords(doRaise, keywords, exposureMetadata, bias, "bias", log=self.log)
1579 self.compareUnits(bias.metadata, "bias")
1580 if self.config.doCrosstalk:
1581 if crosstalk is None:
1582 raise RuntimeError("doCrosstalk is True but no crosstalk provided.")
1583 compareCameraKeywords(doRaise, keywords, exposureMetadata, crosstalk, "crosstalk", log=self.log)
1584 if self.config.doDeferredCharge:
1585 if deferredChargeCalib is None:
1586 raise RuntimeError("doDeferredCharge is True but no deferredChargeCalib provided.")
1588 doRaise,
1589 keywords,
1590 exposureMetadata,
1591 deferredChargeCalib,
1592 "CTI",
1593 log=self.log,
1594 )
1595 if self.config.doDefect:
1596 if defects is None:
1597 raise RuntimeError("doDefect is True but no defects provided.")
1598 compareCameraKeywords(doRaise, keywords, exposureMetadata, defects, "defects", log=self.log)
1599 if self.config.doDark:
1600 if dark is None:
1601 raise RuntimeError("doDark is True but no dark frame provided.")
1602 compareCameraKeywords(doRaise, keywords, exposureMetadata, dark, "dark", log=self.log)
1603 self.compareUnits(bias.metadata, "dark")
1604 if self.config.doBrighterFatter:
1605 if bfKernel is None:
1606 raise RuntimeError("doBrighterFatter is True not no bfKernel provided.")
1607 compareCameraKeywords(doRaise, keywords, exposureMetadata, bfKernel, "bf", log=self.log)
1608 if self.config.doFlat:
1609 if flat is None:
1610 raise RuntimeError("doFlat is True but no flat provided.")
1611 compareCameraKeywords(doRaise, keywords, exposureMetadata, flat, "flat", log=self.log)
1612
1613 if self.config.doSaturation:
1614 if self.config.defaultSaturationSource in ["PTCTURNOFF",]:
1615 if ptc is None:
1616 raise RuntimeError(
1617 "doSaturation is True and defaultSaturationSource is "
1618 f"{self.config.defaultSaturationSource}, but no ptc provided."
1619 )
1620 if self.config.doSuspect:
1621 if self.config.defaultSuspectSource in ["PTCTURNOFF",]:
1622 if ptc is None:
1623 raise RuntimeError(
1624 "doSuspect is True and defaultSuspectSource is "
1625 f"{self.config.defaultSuspectSource}, but no ptc provided."
1626 )
1627
1628 # FIXME: Make sure that if linearity is done then it is matched
1629 # with the right PTC.
1630
1631 # We keep track of units: start in adu.
1632 exposureMetadata["LSST ISR UNITS"] = "adu"
1633 exposureMetadata["LSST ISR CROSSTALK APPLIED"] = False
1634 exposureMetadata["LSST ISR LINEARIZER APPLIED"] = False
1635 exposureMetadata["LSST ISR CTI APPLIED"] = False
1636 exposureMetadata["LSST ISR BIAS APPLIED"] = False
1637 exposureMetadata["LSST ISR DARK APPLIED"] = False
1638 exposureMetadata["LSST ISR BF APPLIED"] = False
1639 exposureMetadata["LSST ISR FLAT APPLIED"] = False
1640 exposureMetadata["LSST ISR DEFECTS APPLIED"] = False
1641
1642 if self.config.doBootstrap:
1643 self.log.info("Configured using doBootstrap=True; using gain of 1.0 (adu units)")
1644 ptc = PhotonTransferCurveDataset([amp.getName() for amp in detector], "NOMINAL_PTC", 1)
1645 for amp in detector:
1646 ptc.gain[amp.getName()] = 1.0
1647 ptc.noise[amp.getName()] = 0.0
1648
1649 exposureMetadata["LSST ISR BOOTSTRAP"] = self.config.doBootstrap
1650
1651 # Set which gains to use
1652 gains = ptc.gain
1653
1654 # And check if we have configured gains to override. This is
1655 # also a warning, since it should not be typical usage.
1656 for amp in detector:
1657 if not math.isnan(gain := overscanDetectorConfig.getOverscanAmpConfig(amp).gain):
1658 gains[amp.getName()] = gain
1659 self.log.warning(
1660 "Overriding gain for amp %s with configured value of %.3f.",
1661 amp.getName(),
1662 gain,
1663 )
1664
1665 # First we convert the exposure to floating point values
1666 # (if necessary).
1667 self.log.debug("Converting exposure to floating point values.")
1668 ccdExposure = self.convertIntToFloat(ccdExposure)
1669
1670 # Then we mark which amplifiers are completely bad from defects.
1671 badAmpDict = self.maskFullAmplifiers(ccdExposure, detector, defects, gains=gains)
1672
1673 # Now we go through ISR steps.
1674
1675 # Differential non-linearity correction.
1676 # Units: adu
1677 if self.config.doDiffNonLinearCorrection:
1678 self.diffNonLinearCorrection(ccdExposure, dnlLUT)
1679
1680 # Dither the integer counts.
1681 # Input units: integerized adu
1682 # Output units: floating-point adu
1683 self.ditherCounts(ccdExposure, overscanDetectorConfig)
1684
1685 # Serial overscan correction.
1686 # Input units: adu
1687 # Output units: adu
1688 if overscanDetectorConfig.doAnySerialOverscan:
1689 serialOverscans = self.overscanCorrection(
1690 "SERIAL",
1691 overscanDetectorConfig,
1692 detector,
1693 badAmpDict,
1694 ccdExposure,
1695 )
1696
1697 if self.config.doBootstrap:
1698 # Get the empirical read noise
1699 for amp, serialOverscan in zip(detector, serialOverscans):
1700 if serialOverscan is None:
1701 ptc.noise[amp.getName()] = 0.0
1702 else:
1703 # All PhotonTransferCurveDataset objects should contain
1704 # noise attributes in units of electrons. The read
1705 # noise measured from overscans is always in adu, so we
1706 # scale it by the gain.
1707 # Note that in bootstrap mode, these gains will always
1708 # be 1.0, but we put this conversion here for clarity.
1709 ptc.noise[amp.getName()] = serialOverscan.residualSigma * gains[amp.getName()]
1710 else:
1711 serialOverscans = [None]*len(detector)
1712
1713 # After serial overscan correction, we can mask SATURATED and
1714 # SUSPECT pixels. This updates badAmpDict if any amplifier
1715 # is fully saturated after serial overscan correction.
1716
1717 # The saturation is currently assumed to be recorded in
1718 # overscan-corrected adu.
1719 badAmpDict = self.maskSaturatedPixels(
1720 badAmpDict,
1721 ccdExposure,
1722 detector,
1723 overscanDetectorConfig,
1724 ptc=ptc,
1725 )
1726
1727 if self.config.doCorrectGains:
1728 # TODO: DM-36639
1729 # This requires the PTC (tbd) with the temperature dependence.
1730 self.log.info("Apply temperature dependence to the gains.")
1731 gains, readNoise = self.correctGains(ccdExposure, ptc, gains)
1732
1733 # Do gain normalization.
1734 # Input units: adu
1735 # Output units: electron
1736 if self.config.doApplyGains:
1737 self.log.info("Using gain values to convert from adu to electron units.")
1738 isrFunctions.applyGains(ccdExposure, normalizeGains=False, ptcGains=gains, isTrimmed=False)
1739 # The units are now electron.
1740 exposureMetadata["LSST ISR UNITS"] = "electron"
1741
1742 # Update the saturation units in the metadata if there.
1743 # These will always have the same units as the image.
1744 for amp in detector:
1745 ampName = amp.getName()
1746 if (key := f"LSST ISR SATURATION LEVEL {ampName}") in exposureMetadata:
1747 exposureMetadata[key] *= gains[ampName]
1748 if (key := f"LSST ISR SUSPECT LEVEL {ampName}") in exposureMetadata:
1749 exposureMetadata[key] *= gains[ampName]
1750
1751 # Record gain and read noise in header.
1752 metadata = ccdExposure.metadata
1753 metadata["LSST ISR READNOISE UNITS"] = "electron"
1754 for amp in detector:
1755 # This includes any gain correction (if applied).
1756 metadata[f"LSST ISR GAIN {amp.getName()}"] = gains[amp.getName()]
1757
1758 # At this stage, the read noise is always in electrons.
1759 noise = ptc.noise[amp.getName()]
1760 metadata[f"LSST ISR READNOISE {amp.getName()}"] = noise
1761
1762 # Do crosstalk correction in the full region.
1763 # Output units: electron (adu if doBootstrap=True)
1764 if self.config.doCrosstalk:
1765 self.log.info("Applying crosstalk corrections to full amplifier region.")
1766 if self.config.doBootstrap and numpy.any(crosstalk.fitGains != 0):
1767 crosstalkGains = None
1768 else:
1769 crosstalkGains = gains
1770 self.crosstalk.run(
1771 ccdExposure,
1772 crosstalk=crosstalk,
1773 gains=crosstalkGains,
1774 fullAmplifier=True,
1775 badAmpDict=badAmpDict,
1776 ignoreVariance=True,
1777 )
1778 ccdExposure.metadata["LSST ISR CROSSTALK APPLIED"] = True
1779
1780 # Parallel overscan correction.
1781 # Output units: electron (adu if doBootstrap=True)
1782 parallelOverscans = None
1783 if overscanDetectorConfig.doAnyParallelOverscan:
1784 # At the moment we do not use the return values from this task.
1785 parallelOverscans = self.overscanCorrection(
1786 "PARALLEL",
1787 overscanDetectorConfig,
1788 detector,
1789 badAmpDict,
1790 ccdExposure,
1791 )
1792
1793 # Linearity correction
1794 # Output units: electron (adu if doBootstrap=True)
1795 if self.config.doLinearize:
1796 self.log.info("Applying linearizer.")
1797 # The linearizer is in units of adu.
1798 # If our units are electron, then pass in the gains
1799 # for conversion.
1800 if exposureMetadata["LSST ISR UNITS"] == "electron":
1801 linearityGains = gains
1802 else:
1803 linearityGains = None
1804 linearizer.applyLinearity(
1805 image=ccdExposure.image,
1806 detector=detector,
1807 log=self.log,
1808 gains=linearityGains,
1809 )
1810 ccdExposure.metadata["LSST ISR LINEARIZER APPLIED"] = True
1811
1812 # Serial CTI (deferred charge) correction
1813 # This will be performed in electron units
1814 # Output units: same as input units
1815 if self.config.doDeferredCharge:
1816 self.deferredChargeCorrection.run(
1817 ccdExposure,
1818 deferredChargeCalib,
1819 gains=gains,
1820 )
1821 ccdExposure.metadata["LSST ISR CTI APPLIED"] = True
1822
1823 # Save the untrimmed version for later statistics,
1824 # which still contains the overscan information
1825 untrimmedCcdExposure = None
1826 if self.config.isrStats.doCtiStatistics:
1827 untrimmedCcdExposure = ccdExposure.clone()
1828
1829 # Assemble/trim
1830 # Output units: electron (adu if doBootstrap=True)
1831 if self.config.doAssembleCcd:
1832 self.log.info("Assembling CCD from amplifiers.")
1833 ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
1834
1835 if self.config.expectWcs and not ccdExposure.getWcs():
1836 self.log.warning("No WCS found in input exposure.")
1837
1838 # ITL Dip Masking
1839 for maskPlane in self.config.itlDipMaskPlanes:
1840 if maskPlane not in ccdExposure.mask.getMaskPlaneDict():
1841 self.log.info("Adding %s mask plane to image.", maskPlane)
1842 ccdExposure.mask.addMaskPlane(maskPlane)
1843
1844 if self.config.doITLDipMask:
1845 isrFunctions.maskITLDip(
1846 exposure=ccdExposure,
1847 detectorConfig=overscanDetectorConfig,
1848 log=self.log,
1849 maskPlaneNames=self.config.itlDipMaskPlanes,
1850 )
1851
1852 if (self.config.doITLSatSagMask or self.config.doITLEdgeBleedMask) \
1853 and detector.getPhysicalType() == 'ITL':
1854 self.maskITLSatEdgesAndColumns(exposure=ccdExposure,
1855 badAmpDict=badAmpDict)
1856
1857 # Bias subtraction
1858 # Output units: electron (adu if doBootstrap=True)
1859 if self.config.doBias:
1860 self.log.info("Applying bias correction.")
1861 # Bias frame and ISR unit consistency is checked at the top of
1862 # the run method.
1863 isrFunctions.biasCorrection(ccdExposure.maskedImage, bias.maskedImage)
1864 ccdExposure.metadata["LSST ISR BIAS APPLIED"] = True
1865
1866 # Dark subtraction
1867 # Output units: electron (adu if doBootstrap=True)
1868 if self.config.doDark:
1869 self.log.info("Applying dark subtraction.")
1870 # Dark frame and ISR unit consistency is checked at the top of
1871 # the run method.
1872 self.darkCorrection(ccdExposure, dark)
1873 ccdExposure.metadata["LSST ISR DARK APPLIED"] = True
1874
1875 # Defect masking
1876 # Masking block (defects, NAN pixels and trails).
1877 # Saturated and suspect pixels have already been masked.
1878 # Output units: electron (adu if doBootstrap=True)
1879 if self.config.doDefect:
1880 self.log.info("Applying defect masking.")
1881 self.maskDefects(ccdExposure, defects)
1882 ccdExposure.metadata["LSST ISR DEFECTS APPLIED"] = True
1883
1884 self.log.info("Adding UNMASKEDNAN mask plane to image.")
1885 ccdExposure.mask.addMaskPlane("UNMASKEDNAN")
1886 if self.config.doNanMasking:
1887 self.log.info("Masking non-finite (NAN, inf) value pixels.")
1888 self.maskNan(ccdExposure)
1889
1890 if self.config.doWidenSaturationTrails:
1891 self.log.info("Widening saturation trails.")
1892 isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
1893
1894 # Brighter-Fatter
1895 # Output units: electron (adu if doBootstrap=True)
1896 if self.config.doBrighterFatter:
1897 self.log.info("Applying brighter-fatter correction.")
1898
1899 bfKernelOut, bfGains = self.getBrighterFatterKernel(detector, bfKernel)
1900
1901 # Needs to be done in electrons; applyBrighterFatterCorrection
1902 # will convert the image if necessary.
1903 if exposureMetadata["LSST ISR UNITS"] == "electron":
1904 brighterFatterApplyGain = False
1905 else:
1906 brighterFatterApplyGain = True
1907
1908 if brighterFatterApplyGain and (ptc is not None) and (bfGains != gains):
1909 # The supplied ptc should be the same as the ptc used to
1910 # generate the bfKernel, in which case they will have the
1911 # same stored amp-keyed dictionary of gains. If not, there
1912 # is a mismatch in the calibrations being used. This should
1913 # not be always be a fatal error, but ideally, everything
1914 # should to be consistent.
1915 self.log.warning("Need to apply gain for brighter-fatter, but the stored"
1916 "gains in the kernel are not the same as the gains stored"
1917 "in the PTC. Using the kernel gains.")
1918
1919 ccdExposure, bfCorrIters = self.applyBrighterFatterCorrection(
1920 ccdExposure,
1921 flat,
1922 dark,
1923 bfKernelOut,
1924 brighterFatterApplyGain,
1925 bfGains,
1926 )
1927
1928 ccdExposure.metadata["LSST ISR BF APPLIED"] = True
1929 metadata["LSST ISR BF ITERS"] = bfCorrIters
1930
1931 # Variance plane creation
1932 # Output units: electron (adu if doBootstrap=True)
1933 if self.config.doVariance:
1934 self.addVariancePlane(ccdExposure, detector)
1935
1936 # Flat-fielding
1937 # This may move elsewhere, but this is the most convenient
1938 # location for simple flat-fielding for attractive backgrounds.
1939 # Output units: electron (adu if doBootstrap=True)
1940 if self.config.doFlat:
1941 self.log.info("Applying flat correction.")
1942 isrFunctions.flatCorrection(
1943 maskedImage=ccdExposure.maskedImage,
1944 flatMaskedImage=flat.maskedImage,
1945 scalingType=self.config.flatScalingType,
1946 userScale=self.config.flatUserScale,
1947 )
1948 ccdExposure.metadata["LSST ISR FLAT APPLIED"] = True
1949 # TODO: DM-49159
1950 # Add metadata re: type of flat.
1951
1952 # Pixel values for masked regions are set here
1953 preInterpExp = None
1954 if self.config.doSaveInterpPixels:
1955 preInterpExp = ccdExposure.clone()
1956
1957 if self.config.doSetBadRegions:
1958 self.log.info('Setting values in large contiguous bad regions.')
1959 self.setBadRegions(ccdExposure)
1960
1961 if self.config.doInterpolate:
1962 self.log.info("Interpolating masked pixels.")
1963 isrFunctions.interpolateFromMask(
1964 maskedImage=ccdExposure.getMaskedImage(),
1965 fwhm=self.config.brighterFatterFwhmForInterpolation,
1966 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1967 maskNameList=list(self.config.maskListToInterpolate),
1968 useLegacyInterp=self.config.useLegacyInterp,
1969 )
1970
1971 # Calculate amp offset corrections within the CCD.
1972 if self.config.doAmpOffset:
1973 if self.config.ampOffset.doApplyAmpOffset:
1974 self.log.info("Measuring and applying amp offset corrections.")
1975 else:
1976 self.log.info("Measuring amp offset corrections only, without applying them.")
1977 self.ampOffset.run(ccdExposure)
1978
1979 # Calculate standard image quality statistics
1980 if self.config.doStandardStatistics:
1981 metadata = ccdExposure.metadata
1982 for amp in detector:
1983 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1984 ampName = amp.getName()
1985 metadata[f"LSST ISR MASK SAT {ampName}"] = isrFunctions.countMaskedPixels(
1986 ampExposure.getMaskedImage(),
1987 [self.config.saturatedMaskName]
1988 )
1989 metadata[f"LSST ISR MASK BAD {ampName}"] = isrFunctions.countMaskedPixels(
1990 ampExposure.getMaskedImage(),
1991 ["BAD"]
1992 )
1993 metadata[f"LSST ISR MASK SUSPECT {ampName}"] = isrFunctions.countMaskedPixels(
1994 ampExposure.getMaskedImage(),
1995 ["SUSPECT"],
1996 )
1997 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
1998 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP)
1999
2000 metadata[f"LSST ISR FINAL MEAN {ampName}"] = qaStats.getValue(afwMath.MEAN)
2001 metadata[f"LSST ISR FINAL MEDIAN {ampName}"] = qaStats.getValue(afwMath.MEDIAN)
2002 metadata[f"LSST ISR FINAL STDEV {ampName}"] = qaStats.getValue(afwMath.STDEVCLIP)
2003
2004 k1 = f"LSST ISR FINAL MEDIAN {ampName}"
2005 k2 = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
2006 if overscanDetectorConfig.doAnySerialOverscan and k1 in metadata and k2 in metadata:
2007 metadata[f"LSST ISR LEVEL {ampName}"] = metadata[k1] - metadata[k2]
2008 else:
2009 metadata[f"LSST ISR LEVEL {ampName}"] = numpy.nan
2010
2011 # calculate additional statistics.
2012 outputStatistics = None
2013 if self.config.doCalculateStatistics:
2014 outputStatistics = self.isrStats.run(ccdExposure,
2015 untrimmedInputExposure=untrimmedCcdExposure,
2016 serialOverscanResults=serialOverscans,
2017 parallelOverscanResults=parallelOverscans,
2018 bias=bias, dark=dark, flat=flat,
2019 ptc=ptc, defects=defects).results
2020
2021 # do image binning.
2022 outputBin1Exposure = None
2023 outputBin2Exposure = None
2024 if self.config.doBinnedExposures:
2025 self.log.info("Creating binned exposures.")
2026 outputBin1Exposure = self.binning.run(
2027 ccdExposure,
2028 binFactor=self.config.binFactor1,
2029 ).binnedExposure
2030 outputBin2Exposure = self.binning.run(
2031 ccdExposure,
2032 binFactor=self.config.binFactor2,
2033 ).binnedExposure
2034
2035 return pipeBase.Struct(
2036 exposure=ccdExposure,
2037
2038 outputBin1Exposure=outputBin1Exposure,
2039 outputBin2Exposure=outputBin2Exposure,
2040
2041 preInterpExposure=preInterpExp,
2042 outputExposure=ccdExposure,
2043 outputStatistics=outputStatistics,
2044 )
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:82
runQuantum(self, butlerQC, inputRefs, outputRefs)
ditherCounts(self, exposure, detectorConfig, fallbackSeed=12345)
flatContext(self, exp, flat, dark=None)
correctGains(self, exposure, ptc, gains)
run(self, ccdExposure, *, dnlLUT=None, bias=None, deferredChargeCalib=None, linearizer=None, ptc=None, crosstalk=None, defects=None, bfKernel=None, bfGains=None, dark=None, flat=None, camera=None, **kwargs)
getBrighterFatterKernel(self, detector, bfKernel)
darkCorrection(self, exposure, darkExposure, invert=False)
compareUnits(self, calibMetadata, calibName)
overscanCorrection(self, mode, detectorConfig, detector, badAmpDict, ccdExposure)
applyBrighterFatterCorrection(self, ccdExposure, flat, dark, bfKernel, brighterFatterApplyGain, bfGains)
addVariancePlane(self, exposure, detector)
diffNonLinearCorrection(self, ccdExposure, dnlLUT, **kwargs)
maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR')
maskDefects(self, exposure, defectBaseList)
maskITLSatEdgesAndColumns(self, exposure, badAmpDict)
maskSaturatedPixels(self, badAmpDict, ccdExposure, detector, detectorConfig, ptc=None)
maskFullAmplifiers(self, ccdExposure, detector, defects, gains=None)
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition Exposure.h:484
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361
std::shared_ptr< ImageT > binImage(ImageT const &inImage, int const binX, int const binY, lsst::afw::math::Property const flags=lsst::afw::math::MEAN)
Definition binImage.cc:44
compareCameraKeywords(doRaiseOnCalibMismatch, cameraKeywordsToCompare, exposureMetadata, calib, calibName, log=None)
size_t maskNans(afw::image::MaskedImage< PixelT > const &mi, afw::image::MaskPixel maskVal, afw::image::MaskPixel allow=0)
Mask NANs in an image.
Definition Isr.cc:35