LSST Applications g034a557a3c+dd8dd8f11d,g0afe43252f+b86e4b8053,g11f7dcd041+017865fdd3,g1cd03abf6b+8446defddb,g1ce3e0751c+f991eae79d,g28da252d5a+ca8a1a9fb3,g2bbee38e9b+b6588ad223,g2bc492864f+b6588ad223,g2cdde0e794+8523d0dbb4,g347aa1857d+b6588ad223,g35bb328faa+b86e4b8053,g3a166c0a6a+b6588ad223,g461a3dce89+b86e4b8053,g52b1c1532d+b86e4b8053,g7f3b0d46df+ad13c1b82d,g80478fca09+f29c5d6c70,g858d7b2824+293f439f82,g8cd86fa7b1+af721d2595,g965a9036f2+293f439f82,g979bb04a14+51ed57f74c,g9ddcbc5298+f24b38b85a,gae0086650b+b86e4b8053,gbb886bcc26+b97e247655,gc28159a63d+b6588ad223,gc30aee3386+a2f0f6cab9,gcaf7e4fdec+293f439f82,gcd45df26be+293f439f82,gcdd4ae20e8+70b5def7e6,gce08ada175+da9c58a417,gcf0d15dbbd+70b5def7e6,gdaeeff99f8+006e14e809,gdbce86181e+6a170ce272,ge3d4d395c2+224150c836,ge5f7162a3a+bb2241c923,ge6cb8fbbf7+d119aed356,ge79ae78c31+b6588ad223,gf048a9a2f4+40ffced2b8,gf0baf85859+b4cca3d10f,w.2024.30
LSST Data Management Base Package
Loading...
Searching...
No Matches
isrTask.py
Go to the documentation of this file.
1# This file is part of ip_isr.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22__all__ = ["IsrTask", "IsrTaskConfig"]
23
24import math
25import numpy
26import numbers
27
28import lsst.geom
29import lsst.afw.image as afwImage
30import lsst.afw.math as afwMath
31import lsst.pex.config as pexConfig
32import lsst.pipe.base as pipeBase
33import lsst.pipe.base.connectionTypes as cT
34
35from contextlib import contextmanager
36from lsstDebug import getDebugFrame
37
38from lsst.afw.cameraGeom import NullLinearityType
39from lsst.afw.display import getDisplay
40from lsst.meas.algorithms.detection import SourceDetectionTask
41from lsst.utils.timer import timeMethod
42
43from . import isrFunctions
44from . import isrQa
45from . import linearize
46from .defects import Defects
47
48from .assembleCcdTask import AssembleCcdTask
49from .crosstalk import CrosstalkTask, CrosstalkCalib
50from .fringe import FringeTask
51from .isr import maskNans
52from .masking import MaskingTask
53from .overscan import OverscanCorrectionTask
54from .straylight import StrayLightTask
55from .vignette import VignetteTask
56from .ampOffset import AmpOffsetTask
57from .deferredCharge import DeferredChargeTask
58from .isrStatistics import IsrStatisticsTask
59from .ptcDataset import PhotonTransferCurveDataset
60
61
62def crosstalkSourceLookup(datasetType, registry, quantumDataId, collections):
63 """Lookup function to identify crosstalkSource entries.
64
65 This should return an empty list under most circumstances. Only
66 when inter-chip crosstalk has been identified should this be
67 populated.
68
69 Parameters
70 ----------
71 datasetType : `str`
72 Dataset to lookup.
73 registry : `lsst.daf.butler.Registry`
74 Butler registry to query.
75 quantumDataId : `lsst.daf.butler.DataCoordinate`
76 Expanded data id to transform to identify crosstalkSources. The
77 ``detector`` entry will be stripped.
78 collections : `lsst.daf.butler.CollectionSearch`
79 Collections to search through.
80
81 Returns
82 -------
83 results : `list` [`lsst.daf.butler.DatasetRef`]
84 List of datasets that match the query that will be used as
85 crosstalkSources.
86 """
87 newDataId = quantumDataId.subset(registry.dimensions.conform(["instrument", "exposure"]))
88 results = set(registry.queryDatasets(datasetType, collections=collections, dataId=newDataId,
89 findFirst=True))
90 # In some contexts, calling `.expanded()` to expand all data IDs in the
91 # query results can be a lot faster because it vectorizes lookups. But in
92 # this case, expandDataId shouldn't need to hit the database at all in the
93 # steady state, because only the detector record is unknown and those are
94 # cached in the registry.
95 records = {k: newDataId.records[k] for k in newDataId.dimensions.elements}
96 return [ref.expanded(registry.expandDataId(ref.dataId, records=records)) for ref in results]
97
98
99class IsrTaskConnections(pipeBase.PipelineTaskConnections,
100 dimensions={"instrument", "exposure", "detector"},
101 defaultTemplates={}):
102 ccdExposure = cT.Input(
103 name="raw",
104 doc="Input exposure to process.",
105 storageClass="Exposure",
106 dimensions=["instrument", "exposure", "detector"],
107 )
108 camera = cT.PrerequisiteInput(
109 name="camera",
110 storageClass="Camera",
111 doc="Input camera to construct complete exposures.",
112 dimensions=["instrument"],
113 isCalibration=True,
114 )
115
116 crosstalk = cT.PrerequisiteInput(
117 name="crosstalk",
118 doc="Input crosstalk object",
119 storageClass="CrosstalkCalib",
120 dimensions=["instrument", "detector"],
121 isCalibration=True,
122 minimum=0, # can fall back to cameraGeom
123 )
124 crosstalkSources = cT.PrerequisiteInput(
125 name="isrOverscanCorrected",
126 doc="Overscan corrected input images.",
127 storageClass="Exposure",
128 dimensions=["instrument", "exposure", "detector"],
129 deferLoad=True,
130 multiple=True,
131 lookupFunction=crosstalkSourceLookup,
132 minimum=0, # not needed for all instruments, no config to control this
133 )
134 bias = cT.PrerequisiteInput(
135 name="bias",
136 doc="Input bias calibration.",
137 storageClass="ExposureF",
138 dimensions=["instrument", "detector"],
139 isCalibration=True,
140 )
141 dark = cT.PrerequisiteInput(
142 name='dark',
143 doc="Input dark calibration.",
144 storageClass="ExposureF",
145 dimensions=["instrument", "detector"],
146 isCalibration=True,
147 )
148 flat = cT.PrerequisiteInput(
149 name="flat",
150 doc="Input flat calibration.",
151 storageClass="ExposureF",
152 dimensions=["instrument", "physical_filter", "detector"],
153 isCalibration=True,
154 )
155 ptc = cT.PrerequisiteInput(
156 name="ptc",
157 doc="Input Photon Transfer Curve dataset",
158 storageClass="PhotonTransferCurveDataset",
159 dimensions=["instrument", "detector"],
160 isCalibration=True,
161 )
162 fringes = cT.PrerequisiteInput(
163 name="fringe",
164 doc="Input fringe calibration.",
165 storageClass="ExposureF",
166 dimensions=["instrument", "physical_filter", "detector"],
167 isCalibration=True,
168 minimum=0, # only needed for some bands, even when enabled
169 )
170 strayLightData = cT.PrerequisiteInput(
171 name='yBackground',
172 doc="Input stray light calibration.",
173 storageClass="StrayLightData",
174 dimensions=["instrument", "physical_filter", "detector"],
175 deferLoad=True,
176 isCalibration=True,
177 minimum=0, # only needed for some bands, even when enabled
178 )
179 bfKernel = cT.PrerequisiteInput(
180 name='bfKernel',
181 doc="Input brighter-fatter kernel.",
182 storageClass="NumpyArray",
183 dimensions=["instrument"],
184 isCalibration=True,
185 minimum=0, # can use either bfKernel or newBFKernel
186 )
187 newBFKernel = cT.PrerequisiteInput(
188 name='brighterFatterKernel',
189 doc="Newer complete kernel + gain solutions.",
190 storageClass="BrighterFatterKernel",
191 dimensions=["instrument", "detector"],
192 isCalibration=True,
193 minimum=0, # can use either bfKernel or newBFKernel
194 )
195 defects = cT.PrerequisiteInput(
196 name='defects',
197 doc="Input defect tables.",
198 storageClass="Defects",
199 dimensions=["instrument", "detector"],
200 isCalibration=True,
201 )
202 linearizer = cT.PrerequisiteInput(
203 name='linearizer',
204 storageClass="Linearizer",
205 doc="Linearity correction calibration.",
206 dimensions=["instrument", "detector"],
207 isCalibration=True,
208 minimum=0, # can fall back to cameraGeom
209 )
210 opticsTransmission = cT.PrerequisiteInput(
211 name="transmission_optics",
212 storageClass="TransmissionCurve",
213 doc="Transmission curve due to the optics.",
214 dimensions=["instrument"],
215 isCalibration=True,
216 )
217 filterTransmission = cT.PrerequisiteInput(
218 name="transmission_filter",
219 storageClass="TransmissionCurve",
220 doc="Transmission curve due to the filter.",
221 dimensions=["instrument", "physical_filter"],
222 isCalibration=True,
223 )
224 sensorTransmission = cT.PrerequisiteInput(
225 name="transmission_sensor",
226 storageClass="TransmissionCurve",
227 doc="Transmission curve due to the sensor.",
228 dimensions=["instrument", "detector"],
229 isCalibration=True,
230 )
231 atmosphereTransmission = cT.PrerequisiteInput(
232 name="transmission_atmosphere",
233 storageClass="TransmissionCurve",
234 doc="Transmission curve due to the atmosphere.",
235 dimensions=["instrument"],
236 isCalibration=True,
237 )
238 illumMaskedImage = cT.PrerequisiteInput(
239 name="illum",
240 doc="Input illumination correction.",
241 storageClass="MaskedImageF",
242 dimensions=["instrument", "physical_filter", "detector"],
243 isCalibration=True,
244 )
245 deferredChargeCalib = cT.PrerequisiteInput(
246 name="cpCtiCalib",
247 doc="Deferred charge/CTI correction dataset.",
248 storageClass="IsrCalib",
249 dimensions=["instrument", "detector"],
250 isCalibration=True,
251 )
252
253 outputExposure = cT.Output(
254 name='postISRCCD',
255 doc="Output ISR processed exposure.",
256 storageClass="Exposure",
257 dimensions=["instrument", "exposure", "detector"],
258 )
259 preInterpExposure = cT.Output(
260 name='preInterpISRCCD',
261 doc="Output ISR processed exposure, with pixels left uninterpolated.",
262 storageClass="ExposureF",
263 dimensions=["instrument", "exposure", "detector"],
264 )
265 outputBin1Exposure = cT.Output(
266 name="postIsrBin1",
267 doc="First binned image.",
268 storageClass="ExposureF",
269 dimensions=["instrument", "exposure", "detector"],
270 )
271 outputBin2Exposure = cT.Output(
272 name="postIsrBin2",
273 doc="Second binned image.",
274 storageClass="ExposureF",
275 dimensions=["instrument", "exposure", "detector"],
276 )
277
278 outputOssThumbnail = cT.Output(
279 name="OssThumb",
280 doc="Output Overscan-subtracted thumbnail image.",
281 storageClass="Thumbnail",
282 dimensions=["instrument", "exposure", "detector"],
283 )
284 outputFlattenedThumbnail = cT.Output(
285 name="FlattenedThumb",
286 doc="Output flat-corrected thumbnail image.",
287 storageClass="Thumbnail",
288 dimensions=["instrument", "exposure", "detector"],
289 )
290 outputStatistics = cT.Output(
291 name="isrStatistics",
292 doc="Output of additional statistics table.",
293 storageClass="StructuredDataDict",
294 dimensions=["instrument", "exposure", "detector"],
295 )
296
297 def __init__(self, *, config=None):
298 super().__init__(config=config)
299
300 if config.doBias is not True:
301 self.prerequisiteInputs.remove("bias")
302 if config.doLinearize is not True:
303 self.prerequisiteInputs.remove("linearizer")
304 if config.doCrosstalk is not True:
305 self.prerequisiteInputs.remove("crosstalkSources")
306 self.prerequisiteInputs.remove("crosstalk")
307 if config.doBrighterFatter is not True:
308 self.prerequisiteInputs.remove("bfKernel")
309 self.prerequisiteInputs.remove("newBFKernel")
310 if config.doDefect is not True:
311 self.prerequisiteInputs.remove("defects")
312 if config.doDark is not True:
313 self.prerequisiteInputs.remove("dark")
314 if config.doFlat is not True:
315 self.prerequisiteInputs.remove("flat")
316 if config.doFringe is not True:
317 self.prerequisiteInputs.remove("fringes")
318 if config.doStrayLight is not True:
319 self.prerequisiteInputs.remove("strayLightData")
320 if config.usePtcGains is not True and config.usePtcReadNoise is not True:
321 self.prerequisiteInputs.remove("ptc")
322 if config.doAttachTransmissionCurve is not True:
323 self.prerequisiteInputs.remove("opticsTransmission")
324 self.prerequisiteInputs.remove("filterTransmission")
325 self.prerequisiteInputs.remove("sensorTransmission")
326 self.prerequisiteInputs.remove("atmosphereTransmission")
327 else:
328 if config.doUseOpticsTransmission is not True:
329 self.prerequisiteInputs.remove("opticsTransmission")
330 if config.doUseFilterTransmission is not True:
331 self.prerequisiteInputs.remove("filterTransmission")
332 if config.doUseSensorTransmission is not True:
333 self.prerequisiteInputs.remove("sensorTransmission")
334 if config.doUseAtmosphereTransmission is not True:
335 self.prerequisiteInputs.remove("atmosphereTransmission")
336 if config.doIlluminationCorrection is not True:
337 self.prerequisiteInputs.remove("illumMaskedImage")
338 if config.doDeferredCharge is not True:
339 self.prerequisiteInputs.remove("deferredChargeCalib")
340
341 if config.doWrite is not True:
342 self.outputs.remove("outputExposure")
343 self.outputs.remove("preInterpExposure")
344 self.outputs.remove("outputFlattenedThumbnail")
345 self.outputs.remove("outputOssThumbnail")
346 self.outputs.remove("outputStatistics")
347 self.outputs.remove("outputBin1Exposure")
348 self.outputs.remove("outputBin2Exposure")
349 else:
350 if config.doBinnedExposures is not True:
351 self.outputs.remove("outputBin1Exposure")
352 self.outputs.remove("outputBin2Exposure")
353 if config.doSaveInterpPixels is not True:
354 self.outputs.remove("preInterpExposure")
355 if config.qa.doThumbnailOss is not True:
356 self.outputs.remove("outputOssThumbnail")
357 if config.qa.doThumbnailFlattened is not True:
358 self.outputs.remove("outputFlattenedThumbnail")
359 if config.doCalculateStatistics is not True:
360 self.outputs.remove("outputStatistics")
361
362
363class IsrTaskConfig(pipeBase.PipelineTaskConfig,
364 pipelineConnections=IsrTaskConnections):
365 """Configuration parameters for IsrTask.
366
367 Items are grouped in the order in which they are executed by the task.
368 """
369 datasetType = pexConfig.Field(
370 dtype=str,
371 doc="Dataset type for input data; users will typically leave this alone, "
372 "but camera-specific ISR tasks will override it",
373 default="raw",
374 )
375
376 fallbackFilterName = pexConfig.Field(
377 dtype=str,
378 doc="Fallback default filter name for calibrations.",
379 optional=True
380 )
381 useFallbackDate = pexConfig.Field(
382 dtype=bool,
383 doc="Pass observation date when using fallback filter.",
384 default=False,
385 )
386 expectWcs = pexConfig.Field(
387 dtype=bool,
388 default=True,
389 doc="Expect input science images to have a WCS (set False for e.g. spectrographs)."
390 )
391 fwhm = pexConfig.Field(
392 dtype=float,
393 doc="FWHM of PSF in arcseconds (currently unused).",
394 default=1.0,
395 )
396 qa = pexConfig.ConfigField(
397 dtype=isrQa.IsrQaConfig,
398 doc="QA related configuration options.",
399 )
400 doHeaderProvenance = pexConfig.Field(
401 dtype=bool,
402 default=True,
403 doc="Write calibration identifiers into output exposure header?",
404 )
405
406 # Calib checking configuration:
407 doRaiseOnCalibMismatch = pexConfig.Field(
408 dtype=bool,
409 default=False,
410 doc="Should IsrTask halt if exposure and calibration header values do not match?",
411 )
412 cameraKeywordsToCompare = pexConfig.ListField(
413 dtype=str,
414 doc="List of header keywords to compare between exposure and calibrations.",
415 default=[],
416 )
417
418 # Image conversion configuration
419 doConvertIntToFloat = pexConfig.Field(
420 dtype=bool,
421 doc="Convert integer raw images to floating point values?",
422 default=True,
423 )
424
425 # Saturated pixel handling.
426 doSaturation = pexConfig.Field(
427 dtype=bool,
428 doc="Mask saturated pixels? NB: this is totally independent of the"
429 " interpolation option - this is ONLY setting the bits in the mask."
430 " To have them interpolated make sure doSaturationInterpolation=True",
431 default=True,
432 )
433 saturatedMaskName = pexConfig.Field(
434 dtype=str,
435 doc="Name of mask plane to use in saturation detection and interpolation",
436 default="SAT",
437 )
438 saturation = pexConfig.Field(
439 dtype=float,
440 doc="The saturation level to use if no Detector is present in the Exposure (ignored if NaN)",
441 default=float("NaN"),
442 )
443 growSaturationFootprintSize = pexConfig.Field(
444 dtype=int,
445 doc="Number of pixels by which to grow the saturation footprints",
446 default=1,
447 )
448
449 # Suspect pixel handling.
450 doSuspect = pexConfig.Field(
451 dtype=bool,
452 doc="Mask suspect pixels?",
453 default=False,
454 )
455 suspectMaskName = pexConfig.Field(
456 dtype=str,
457 doc="Name of mask plane to use for suspect pixels",
458 default="SUSPECT",
459 )
460 numEdgeSuspect = pexConfig.Field(
461 dtype=int,
462 doc="Number of edge pixels to be flagged as untrustworthy.",
463 default=0,
464 )
465 edgeMaskLevel = pexConfig.ChoiceField(
466 dtype=str,
467 doc="Mask edge pixels in which coordinate frame: DETECTOR or AMP?",
468 default="DETECTOR",
469 allowed={
470 'DETECTOR': 'Mask only the edges of the full detector.',
471 'AMP': 'Mask edges of each amplifier.',
472 },
473 )
474
475 # Initial masking options.
476 doSetBadRegions = pexConfig.Field(
477 dtype=bool,
478 doc="Should we set the level of all BAD patches of the chip to the chip's average value?",
479 default=True,
480 )
481 badStatistic = pexConfig.ChoiceField(
482 dtype=str,
483 doc="How to estimate the average value for BAD regions.",
484 default='MEANCLIP',
485 allowed={
486 "MEANCLIP": "Correct using the (clipped) mean of good data",
487 "MEDIAN": "Correct using the median of the good data",
488 },
489 )
490
491 # Overscan subtraction configuration.
492 doOverscan = pexConfig.Field(
493 dtype=bool,
494 doc="Do overscan subtraction?",
495 default=True,
496 )
497 overscan = pexConfig.ConfigurableField(
498 target=OverscanCorrectionTask,
499 doc="Overscan subtraction task for image segments.",
500 )
501
502 # Amplifier to CCD assembly configuration
503 doAssembleCcd = pexConfig.Field(
504 dtype=bool,
505 default=True,
506 doc="Assemble amp-level exposures into a ccd-level exposure?"
507 )
508 assembleCcd = pexConfig.ConfigurableField(
509 target=AssembleCcdTask,
510 doc="CCD assembly task",
511 )
512
513 # General calibration configuration.
514 doAssembleIsrExposures = pexConfig.Field(
515 dtype=bool,
516 default=False,
517 doc="Assemble amp-level calibration exposures into ccd-level exposure?"
518 )
519 doTrimToMatchCalib = pexConfig.Field(
520 dtype=bool,
521 default=False,
522 doc="Trim raw data to match calibration bounding boxes?"
523 )
524
525 # Bias subtraction.
526 doBias = pexConfig.Field(
527 dtype=bool,
528 doc="Apply bias frame correction?",
529 default=True,
530 )
531 biasDataProductName = pexConfig.Field(
532 dtype=str,
533 doc="Name of the bias data product",
534 default="bias",
535 )
536 doBiasBeforeOverscan = pexConfig.Field(
537 dtype=bool,
538 doc="Reverse order of overscan and bias correction.",
539 default=False
540 )
541
542 # Deferred charge correction.
543 doDeferredCharge = pexConfig.Field(
544 dtype=bool,
545 doc="Apply deferred charge correction?",
546 default=False,
547 )
548 deferredChargeCorrection = pexConfig.ConfigurableField(
549 target=DeferredChargeTask,
550 doc="Deferred charge correction task.",
551 )
552
553 # Variance construction
554 doVariance = pexConfig.Field(
555 dtype=bool,
556 doc="Calculate variance?",
557 default=True
558 )
559 gain = pexConfig.Field(
560 dtype=float,
561 doc="The gain to use if no Detector is present in the Exposure (ignored if NaN)",
562 default=float("NaN"),
563 )
564 readNoise = pexConfig.Field(
565 dtype=float,
566 doc="The read noise to use if no Detector is present in the Exposure",
567 default=0.0,
568 )
569 doEmpiricalReadNoise = pexConfig.Field(
570 dtype=bool,
571 default=False,
572 doc="Calculate empirical read noise instead of value from AmpInfo data?"
573 )
574 usePtcReadNoise = pexConfig.Field(
575 dtype=bool,
576 default=False,
577 doc="Use readnoise values from the Photon Transfer Curve?"
578 )
579 maskNegativeVariance = pexConfig.Field(
580 dtype=bool,
581 default=True,
582 doc="Mask pixels that claim a negative variance? This likely indicates a failure "
583 "in the measurement of the overscan at an edge due to the data falling off faster "
584 "than the overscan model can account for it."
585 )
586 negativeVarianceMaskName = pexConfig.Field(
587 dtype=str,
588 default="BAD",
589 doc="Mask plane to use to mark pixels with negative variance, if `maskNegativeVariance` is True.",
590 )
591 # Linearization.
592 doLinearize = pexConfig.Field(
593 dtype=bool,
594 doc="Correct for nonlinearity of the detector's response?",
595 default=True,
596 )
597
598 # Crosstalk.
599 doCrosstalk = pexConfig.Field(
600 dtype=bool,
601 doc="Apply intra-CCD crosstalk correction?",
602 default=False,
603 )
604 doCrosstalkBeforeAssemble = pexConfig.Field(
605 dtype=bool,
606 doc="Apply crosstalk correction before CCD assembly, and before trimming?",
607 default=False,
608 )
609 crosstalk = pexConfig.ConfigurableField(
610 target=CrosstalkTask,
611 doc="Intra-CCD crosstalk correction",
612 )
613
614 # Masking options.
615 doDefect = pexConfig.Field(
616 dtype=bool,
617 doc="Apply correction for CCD defects, e.g. hot pixels?",
618 default=True,
619 )
620 doNanMasking = pexConfig.Field(
621 dtype=bool,
622 doc="Mask non-finite (NAN, inf) pixels?",
623 default=True,
624 )
625 doWidenSaturationTrails = pexConfig.Field(
626 dtype=bool,
627 doc="Widen bleed trails based on their width?",
628 default=True
629 )
630
631 # Brighter-Fatter correction.
632 doBrighterFatter = pexConfig.Field(
633 dtype=bool,
634 default=False,
635 doc="Apply the brighter-fatter correction?"
636 )
637 doFluxConservingBrighterFatterCorrection = pexConfig.Field(
638 dtype=bool,
639 default=False,
640 doc="Apply the flux-conserving BFE correction by Miller et al.?"
641 )
642 brighterFatterLevel = pexConfig.ChoiceField(
643 dtype=str,
644 default="DETECTOR",
645 doc="The level at which to correct for brighter-fatter.",
646 allowed={
647 "AMP": "Every amplifier treated separately.",
648 "DETECTOR": "One kernel per detector",
649 }
650 )
651 brighterFatterMaxIter = pexConfig.Field(
652 dtype=int,
653 default=10,
654 doc="Maximum number of iterations for the brighter-fatter correction"
655 )
656 brighterFatterThreshold = pexConfig.Field(
657 dtype=float,
658 default=1000,
659 doc="Threshold used to stop iterating the brighter-fatter correction. It is the "
660 "absolute value of the difference between the current corrected image and the one "
661 "from the previous iteration summed over all the pixels."
662 )
663 brighterFatterApplyGain = pexConfig.Field(
664 dtype=bool,
665 default=True,
666 doc="Should the gain be applied when applying the brighter-fatter correction?"
667 )
668 brighterFatterMaskListToInterpolate = pexConfig.ListField(
669 dtype=str,
670 doc="List of mask planes that should be interpolated over when applying the brighter-fatter "
671 "correction.",
672 default=["SAT", "BAD", "NO_DATA", "UNMASKEDNAN"],
673 )
674 brighterFatterMaskGrowSize = pexConfig.Field(
675 dtype=int,
676 default=0,
677 doc="Number of pixels to grow the masks listed in config.brighterFatterMaskListToInterpolate "
678 "when brighter-fatter correction is applied."
679 )
680
681 # Dark subtraction.
682 doDark = pexConfig.Field(
683 dtype=bool,
684 doc="Apply dark frame correction?",
685 default=True,
686 )
687 darkDataProductName = pexConfig.Field(
688 dtype=str,
689 doc="Name of the dark data product",
690 default="dark",
691 )
692
693 # Camera-specific stray light removal.
694 doStrayLight = pexConfig.Field(
695 dtype=bool,
696 doc="Subtract stray light in the y-band (due to encoder LEDs)?",
697 default=False,
698 )
699 strayLight = pexConfig.ConfigurableField(
700 target=StrayLightTask,
701 doc="y-band stray light correction"
702 )
703
704 # Flat correction.
705 doFlat = pexConfig.Field(
706 dtype=bool,
707 doc="Apply flat field correction?",
708 default=True,
709 )
710 flatDataProductName = pexConfig.Field(
711 dtype=str,
712 doc="Name of the flat data product",
713 default="flat",
714 )
715 flatScalingType = pexConfig.ChoiceField(
716 dtype=str,
717 doc="The method for scaling the flat on the fly.",
718 default='USER',
719 allowed={
720 "USER": "Scale by flatUserScale",
721 "MEAN": "Scale by the inverse of the mean",
722 "MEDIAN": "Scale by the inverse of the median",
723 },
724 )
725 flatUserScale = pexConfig.Field(
726 dtype=float,
727 doc="If flatScalingType is 'USER' then scale flat by this amount; ignored otherwise",
728 default=1.0,
729 )
730 doTweakFlat = pexConfig.Field(
731 dtype=bool,
732 doc="Tweak flats to match observed amplifier ratios?",
733 default=False
734 )
735
736 # Amplifier normalization based on gains instead of using flats
737 # configuration.
738 doApplyGains = pexConfig.Field(
739 dtype=bool,
740 doc="Correct the amplifiers for their gains instead of applying flat correction",
741 default=False,
742 )
743 usePtcGains = pexConfig.Field(
744 dtype=bool,
745 doc="Use the gain values from the input Photon Transfer Curve?",
746 default=False,
747 )
748 normalizeGains = pexConfig.Field(
749 dtype=bool,
750 doc="Normalize all the amplifiers in each CCD to have the same median value.",
751 default=False,
752 )
753
754 # Fringe correction.
755 doFringe = pexConfig.Field(
756 dtype=bool,
757 doc="Apply fringe correction?",
758 default=True,
759 )
760 fringe = pexConfig.ConfigurableField(
761 target=FringeTask,
762 doc="Fringe subtraction task",
763 )
764 fringeAfterFlat = pexConfig.Field(
765 dtype=bool,
766 doc="Do fringe subtraction after flat-fielding?",
767 default=True,
768 )
769
770 # Amp offset correction.
771 doAmpOffset = pexConfig.Field(
772 doc="Calculate and apply amp offset corrections?",
773 dtype=bool,
774 default=False,
775 )
776 ampOffset = pexConfig.ConfigurableField(
777 doc="Amp offset correction task.",
778 target=AmpOffsetTask,
779 )
780
781 # Initial CCD-level background statistics options.
782 doMeasureBackground = pexConfig.Field(
783 dtype=bool,
784 doc="Measure the background level on the reduced image?",
785 default=False,
786 )
787
788 # Camera-specific masking configuration.
789 doCameraSpecificMasking = pexConfig.Field(
790 dtype=bool,
791 doc="Mask camera-specific bad regions?",
792 default=False,
793 )
794 masking = pexConfig.ConfigurableField(
795 target=MaskingTask,
796 doc="Masking task."
797 )
798
799 # Interpolation options.
800 doInterpolate = pexConfig.Field(
801 dtype=bool,
802 doc="Interpolate masked pixels?",
803 default=True,
804 )
805 doSaturationInterpolation = pexConfig.Field(
806 dtype=bool,
807 doc="Perform interpolation over pixels masked as saturated?"
808 " NB: This is independent of doSaturation; if that is False this plane"
809 " will likely be blank, resulting in a no-op here.",
810 default=True,
811 )
812 doNanInterpolation = pexConfig.Field(
813 dtype=bool,
814 doc="Perform interpolation over pixels masked as NaN?"
815 " NB: This is independent of doNanMasking; if that is False this plane"
816 " will likely be blank, resulting in a no-op here.",
817 default=True,
818 )
819 doNanInterpAfterFlat = pexConfig.Field(
820 dtype=bool,
821 doc=("If True, ensure we interpolate NaNs after flat-fielding, even if we "
822 "also have to interpolate them before flat-fielding."),
823 default=False,
824 )
825 maskListToInterpolate = pexConfig.ListField(
826 dtype=str,
827 doc="List of mask planes that should be interpolated.",
828 default=['SAT', 'BAD'],
829 )
830 doSaveInterpPixels = pexConfig.Field(
831 dtype=bool,
832 doc="Save a copy of the pre-interpolated pixel values?",
833 default=False,
834 )
835
836 # Default photometric calibration options.
837 fluxMag0T1 = pexConfig.DictField(
838 keytype=str,
839 itemtype=float,
840 doc="The approximate flux of a zero-magnitude object in a one-second exposure, per filter.",
841 default=dict((f, pow(10.0, 0.4*m)) for f, m in (("Unknown", 28.0),
842 ))
843 )
844 defaultFluxMag0T1 = pexConfig.Field(
845 dtype=float,
846 doc="Default value for fluxMag0T1 (for an unrecognized filter).",
847 default=pow(10.0, 0.4*28.0)
848 )
849
850 # Vignette correction configuration.
851 doVignette = pexConfig.Field(
852 dtype=bool,
853 doc=("Compute and attach the validPolygon defining the unvignetted region to the exposure "
854 "according to vignetting parameters?"),
855 default=False,
856 )
857 doMaskVignettePolygon = pexConfig.Field(
858 dtype=bool,
859 doc=("Add a mask bit for pixels within the vignetted region. Ignored if doVignette "
860 "is False"),
861 default=True,
862 )
863 vignetteValue = pexConfig.Field(
864 dtype=float,
865 doc="Value to replace image array pixels with in the vignetted region? Ignored if None.",
866 optional=True,
867 default=None,
868 )
869 vignette = pexConfig.ConfigurableField(
870 target=VignetteTask,
871 doc="Vignetting task.",
872 )
873
874 # Transmission curve configuration.
875 doAttachTransmissionCurve = pexConfig.Field(
876 dtype=bool,
877 default=False,
878 doc="Construct and attach a wavelength-dependent throughput curve for this CCD image?"
879 )
880 doUseOpticsTransmission = pexConfig.Field(
881 dtype=bool,
882 default=True,
883 doc="Load and use transmission_optics (if doAttachTransmissionCurve is True)?"
884 )
885 doUseFilterTransmission = pexConfig.Field(
886 dtype=bool,
887 default=True,
888 doc="Load and use transmission_filter (if doAttachTransmissionCurve is True)?"
889 )
890 doUseSensorTransmission = pexConfig.Field(
891 dtype=bool,
892 default=True,
893 doc="Load and use transmission_sensor (if doAttachTransmissionCurve is True)?"
894 )
895 doUseAtmosphereTransmission = pexConfig.Field(
896 dtype=bool,
897 default=True,
898 doc="Load and use transmission_atmosphere (if doAttachTransmissionCurve is True)?"
899 )
900
901 # Illumination correction.
902 doIlluminationCorrection = pexConfig.Field(
903 dtype=bool,
904 default=False,
905 doc="Perform illumination correction?"
906 )
907 illuminationCorrectionDataProductName = pexConfig.Field(
908 dtype=str,
909 doc="Name of the illumination correction data product.",
910 default="illumcor",
911 )
912 illumScale = pexConfig.Field(
913 dtype=float,
914 doc="Scale factor for the illumination correction.",
915 default=1.0,
916 )
917 illumFilters = pexConfig.ListField(
918 dtype=str,
919 default=[],
920 doc="Only perform illumination correction for these filters."
921 )
922
923 # Calculate image quality statistics?
924 doStandardStatistics = pexConfig.Field(
925 dtype=bool,
926 doc="Should standard image quality statistics be calculated?",
927 default=True,
928 )
929 # Calculate additional statistics?
930 doCalculateStatistics = pexConfig.Field(
931 dtype=bool,
932 doc="Should additional ISR statistics be calculated?",
933 default=False,
934 )
935 isrStats = pexConfig.ConfigurableField(
936 target=IsrStatisticsTask,
937 doc="Task to calculate additional statistics.",
938 )
939
940 # Make binned images?
941 doBinnedExposures = pexConfig.Field(
942 dtype=bool,
943 doc="Should binned exposures be calculated?",
944 default=False,
945 )
946 binFactor1 = pexConfig.Field(
947 dtype=int,
948 doc="Binning factor for first binned exposure. This is intended for a finely binned output.",
949 default=8,
950 check=lambda x: x > 1,
951 )
952 binFactor2 = pexConfig.Field(
953 dtype=int,
954 doc="Binning factor for second binned exposure. This is intended for a coarsely binned output.",
955 default=64,
956 check=lambda x: x > 1,
957 )
958
959 # Write the outputs to disk. If ISR is run as a subtask, this may not
960 # be needed.
961 doWrite = pexConfig.Field(
962 dtype=bool,
963 doc="Persist postISRCCD?",
964 default=True,
965 )
966
967 def validate(self):
968 super().validate()
969 if self.doFlat and self.doApplyGains:
970 raise ValueError("You may not specify both doFlat and doApplyGains")
972 raise ValueError("You may not specify both doBiasBeforeOverscan and doTrimToMatchCalib")
977 if self.doNanInterpolation and "UNMASKEDNAN" not in self.maskListToInterpolate:
978 self.maskListToInterpolate.append("UNMASKEDNAN")
979 if self.doCalculateStatistics and self.isrStats.doCtiStatistics:
980 if self.doApplyGains != self.isrStats.doApplyGainsForCtiStatistics:
981 raise ValueError("doApplyGains must match isrStats.applyGainForCtiStatistics.")
982
983
984class IsrTask(pipeBase.PipelineTask):
985 """Apply common instrument signature correction algorithms to a raw frame.
986
987 The process for correcting imaging data is very similar from
988 camera to camera. This task provides a vanilla implementation of
989 doing these corrections, including the ability to turn certain
990 corrections off if they are not needed. The inputs to the primary
991 method, `run()`, are a raw exposure to be corrected and the
992 calibration data products. The raw input is a single chip sized
993 mosaic of all amps including overscans and other non-science
994 pixels.
995
996 The __init__ method sets up the subtasks for ISR processing, using
997 the defaults from `lsst.ip.isr`.
998
999 Parameters
1000 ----------
1001 args : `list`
1002 Positional arguments passed to the Task constructor.
1003 None used at this time.
1004 kwargs : `dict`, optional
1005 Keyword arguments passed on to the Task constructor.
1006 None used at this time.
1007 """
1008 ConfigClass = IsrTaskConfig
1009 _DefaultName = "isr"
1010
1011 def __init__(self, **kwargs):
1012 super().__init__(**kwargs)
1013 self.makeSubtask("assembleCcd")
1014 self.makeSubtask("crosstalk")
1015 self.makeSubtask("strayLight")
1016 self.makeSubtask("fringe")
1017 self.makeSubtask("masking")
1018 self.makeSubtask("overscan")
1019 self.makeSubtask("vignette")
1020 self.makeSubtask("ampOffset")
1021 self.makeSubtask("deferredChargeCorrection")
1022 self.makeSubtask("isrStats")
1023
1024 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1025 inputs = butlerQC.get(inputRefs)
1026
1027 try:
1028 inputs['detectorNum'] = inputRefs.ccdExposure.dataId['detector']
1029 except Exception as e:
1030 raise ValueError("Failure to find valid detectorNum value for Dataset %s: %s." %
1031 (inputRefs, e))
1032
1033 detector = inputs['ccdExposure'].getDetector()
1034
1035 # This is use for header provenance.
1036 additionalInputDates = {}
1037
1038 if self.config.doCrosstalk is True:
1039 # Crosstalk sources need to be defined by the pipeline
1040 # yaml if they exist.
1041 if 'crosstalk' in inputs and inputs['crosstalk'] is not None:
1042 if not isinstance(inputs['crosstalk'], CrosstalkCalib):
1043 inputs['crosstalk'] = CrosstalkCalib.fromTable(inputs['crosstalk'])
1044 else:
1045 coeffVector = (self.config.crosstalk.crosstalkValues
1046 if self.config.crosstalk.useConfigCoefficients else None)
1047 crosstalkCalib = CrosstalkCalib().fromDetector(detector, coeffVector=coeffVector)
1048 inputs['crosstalk'] = crosstalkCalib
1049 if inputs['crosstalk'].interChip and len(inputs['crosstalk'].interChip) > 0:
1050 if 'crosstalkSources' not in inputs:
1051 self.log.warning("No crosstalkSources found for chip with interChip terms!")
1052
1053 if self.doLinearize(detector) is True:
1054 if 'linearizer' in inputs:
1055 if isinstance(inputs['linearizer'], dict):
1056 linearizer = linearize.Linearizer(detector=detector, log=self.log)
1057 linearizer.fromYaml(inputs['linearizer'])
1058 self.log.warning("Dictionary linearizers will be deprecated in DM-28741.")
1059 elif isinstance(inputs['linearizer'], numpy.ndarray):
1060 linearizer = linearize.Linearizer(table=inputs.get('linearizer', None),
1061 detector=detector,
1062 log=self.log)
1063 self.log.warning("Bare lookup table linearizers will be deprecated in DM-28741.")
1064 else:
1065 linearizer = inputs['linearizer']
1066 self.log.info("Loading linearizer from the Butler.")
1067 linearizer.log = self.log
1068 inputs['linearizer'] = linearizer
1069 else:
1070 inputs['linearizer'] = linearize.Linearizer(detector=detector, log=self.log)
1071 self.log.info("Constructing linearizer from cameraGeom information.")
1072
1073 if self.config.doDefect is True:
1074 if "defects" in inputs and inputs['defects'] is not None:
1075 # defects is loaded as a BaseCatalog with columns
1076 # x0, y0, width, height. Masking expects a list of defects
1077 # defined by their bounding box
1078 if not isinstance(inputs["defects"], Defects):
1079 inputs["defects"] = Defects.fromTable(inputs["defects"])
1080
1081 # Load the correct style of brighter-fatter kernel, and repack
1082 # the information as a numpy array.
1083 brighterFatterSource = None
1084 if self.config.doBrighterFatter:
1085 brighterFatterKernel = inputs.pop('newBFKernel', None)
1086 if brighterFatterKernel is None:
1087 # This type of kernel must be in (y, x) index
1088 # ordering, as it used directly as the .array
1089 # component of the afwImage kernel.
1090 brighterFatterKernel = inputs.get('bfKernel', None)
1091 brighterFatterSource = 'bfKernel'
1092 additionalInputDates[brighterFatterSource] = self.extractCalibDate(brighterFatterKernel)
1093
1094 if brighterFatterKernel is None:
1095 # This was requested by the config, but none were found.
1096 raise RuntimeError("No brighter-fatter kernel was supplied.")
1097 elif not isinstance(brighterFatterKernel, numpy.ndarray):
1098 # This is a ISR calib kernel. These kernels are
1099 # generated in (x, y) index ordering, and need to be
1100 # transposed to be used directly as the .array
1101 # component of the afwImage kernel. This is done
1102 # explicitly below when setting the ``bfKernel``
1103 # input.
1104 brighterFatterSource = 'newBFKernel'
1105 additionalInputDates[brighterFatterSource] = self.extractCalibDate(brighterFatterKernel)
1106
1107 detName = detector.getName()
1108 level = brighterFatterKernel.level
1109
1110 # This is expected to be a dictionary of amp-wise gains.
1111 inputs['bfGains'] = brighterFatterKernel.gain
1112 if self.config.brighterFatterLevel == 'DETECTOR':
1113 kernel = None
1114 if level == 'DETECTOR':
1115 if detName in brighterFatterKernel.detKernels:
1116 kernel = brighterFatterKernel.detKernels[detName]
1117 else:
1118 raise RuntimeError("Failed to extract kernel from new-style BF kernel.")
1119 elif level == 'AMP':
1120 self.log.warning("Making DETECTOR level kernel from AMP based brighter "
1121 "fatter kernels.")
1122 brighterFatterKernel.makeDetectorKernelFromAmpwiseKernels(detName)
1123 kernel = brighterFatterKernel.detKernels[detName]
1124 if kernel is None:
1125 raise RuntimeError("Could not identify brighter-fatter kernel!")
1126 # Do the one single transpose here so the kernel
1127 # can be directly loaded into the afwImage .array
1128 # component.
1129 inputs['bfKernel'] = numpy.transpose(kernel)
1130 elif self.config.brighterFatterLevel == 'AMP':
1131 raise NotImplementedError("Per-amplifier brighter-fatter correction not implemented")
1132
1133 if self.config.doFringe is True and self.fringe.checkFilter(inputs['ccdExposure']):
1134 expId = inputs['ccdExposure'].info.id
1135 inputs['fringes'] = self.fringe.loadFringes(inputs['fringes'],
1136 expId=expId,
1137 assembler=self.assembleCcd
1138 if self.config.doAssembleIsrExposures else None)
1139 else:
1140 inputs['fringes'] = pipeBase.Struct(fringes=None)
1141
1142 if self.config.doStrayLight is True and self.strayLight.checkFilter(inputs['ccdExposure']):
1143 if 'strayLightData' not in inputs:
1144 inputs['strayLightData'] = None
1145
1146 if self.config.doHeaderProvenance:
1147 # Add calibration provenanace info to header.
1148 exposureMetadata = inputs['ccdExposure'].getMetadata()
1149
1150 # These inputs change name during this step. These should
1151 # have matching entries in the additionalInputDates dict.
1152 additionalInputs = []
1153 if self.config.doBrighterFatter:
1154 additionalInputs.append(brighterFatterSource)
1155
1156 for inputName in sorted(list(inputs.keys()) + additionalInputs):
1157 reference = getattr(inputRefs, inputName, None)
1158 if reference is not None and hasattr(reference, "run"):
1159 runKey = f"LSST CALIB RUN {inputName.upper()}"
1160 runValue = reference.run
1161 idKey = f"LSST CALIB UUID {inputName.upper()}"
1162 idValue = str(reference.id)
1163 dateKey = f"LSST CALIB DATE {inputName.upper()}"
1164
1165 if inputName in additionalInputDates:
1166 dateValue = additionalInputDates[inputName]
1167 else:
1168 dateValue = self.extractCalibDate(inputs[inputName])
1169
1170 exposureMetadata[runKey] = runValue
1171 exposureMetadata[idKey] = idValue
1172 exposureMetadata[dateKey] = dateValue
1173
1174 outputs = self.run(**inputs)
1175 butlerQC.put(outputs, outputRefs)
1176
1177 @timeMethod
1178 def run(self, ccdExposure, *, camera=None, bias=None, linearizer=None,
1179 crosstalk=None, crosstalkSources=None,
1180 dark=None, flat=None, ptc=None, bfKernel=None, bfGains=None, defects=None,
1181 fringes=pipeBase.Struct(fringes=None), opticsTransmission=None, filterTransmission=None,
1182 sensorTransmission=None, atmosphereTransmission=None,
1183 detectorNum=None, strayLightData=None, illumMaskedImage=None,
1184 deferredChargeCalib=None,
1185 ):
1186 """Perform instrument signature removal on an exposure.
1187
1188 Steps included in the ISR processing, in order performed, are:
1189
1190 - saturation and suspect pixel masking
1191 - overscan subtraction
1192 - CCD assembly of individual amplifiers
1193 - bias subtraction
1194 - variance image construction
1195 - linearization of non-linear response
1196 - crosstalk masking
1197 - brighter-fatter correction
1198 - dark subtraction
1199 - fringe correction
1200 - stray light subtraction
1201 - flat correction
1202 - masking of known defects and camera specific features
1203 - vignette calculation
1204 - appending transmission curve and distortion model
1205
1206 Parameters
1207 ----------
1208 ccdExposure : `lsst.afw.image.Exposure`
1209 The raw exposure that is to be run through ISR. The
1210 exposure is modified by this method.
1211 camera : `lsst.afw.cameraGeom.Camera`, optional
1212 The camera geometry for this exposure. Required if
1213 one or more of ``ccdExposure``, ``bias``, ``dark``, or
1214 ``flat`` does not have an associated detector.
1215 bias : `lsst.afw.image.Exposure`, optional
1216 Bias calibration frame.
1217 linearizer : `lsst.ip.isr.linearize.LinearizeBase`, optional
1218 Functor for linearization.
1219 crosstalk : `lsst.ip.isr.crosstalk.CrosstalkCalib`, optional
1220 Calibration for crosstalk.
1221 crosstalkSources : `list`, optional
1222 List of possible crosstalk sources.
1223 dark : `lsst.afw.image.Exposure`, optional
1224 Dark calibration frame.
1225 flat : `lsst.afw.image.Exposure`, optional
1226 Flat calibration frame.
1227 ptc : `lsst.ip.isr.PhotonTransferCurveDataset`, optional
1228 Photon transfer curve dataset, with, e.g., gains
1229 and read noise.
1230 bfKernel : `numpy.ndarray`, optional
1231 Brighter-fatter kernel.
1232 bfGains : `dict` of `float`, optional
1233 Gains used to override the detector's nominal gains for the
1234 brighter-fatter correction. A dict keyed by amplifier name for
1235 the detector in question.
1236 defects : `lsst.ip.isr.Defects`, optional
1237 List of defects.
1238 fringes : `lsst.pipe.base.Struct`, optional
1239 Struct containing the fringe correction data, with
1240 elements:
1241
1242 ``fringes``
1243 fringe calibration frame (`lsst.afw.image.Exposure`)
1244 ``seed``
1245 random seed derived from the ``ccdExposureId`` for random
1246 number generator (`numpy.uint32`)
1247 opticsTransmission: `lsst.afw.image.TransmissionCurve`, optional
1248 A ``TransmissionCurve`` that represents the throughput of the,
1249 optics, to be evaluated in focal-plane coordinates.
1250 filterTransmission : `lsst.afw.image.TransmissionCurve`
1251 A ``TransmissionCurve`` that represents the throughput of the
1252 filter itself, to be evaluated in focal-plane coordinates.
1253 sensorTransmission : `lsst.afw.image.TransmissionCurve`
1254 A ``TransmissionCurve`` that represents the throughput of the
1255 sensor itself, to be evaluated in post-assembly trimmed detector
1256 coordinates.
1257 atmosphereTransmission : `lsst.afw.image.TransmissionCurve`
1258 A ``TransmissionCurve`` that represents the throughput of the
1259 atmosphere, assumed to be spatially constant.
1260 detectorNum : `int`, optional
1261 The integer number for the detector to process.
1262 strayLightData : `object`, optional
1263 Opaque object containing calibration information for stray-light
1264 correction. If `None`, no correction will be performed.
1265 illumMaskedImage : `lsst.afw.image.MaskedImage`, optional
1266 Illumination correction image.
1267
1268 Returns
1269 -------
1270 result : `lsst.pipe.base.Struct`
1271 Result struct with component:
1272
1273 ``exposure``
1274 The fully ISR corrected exposure.
1275 (`lsst.afw.image.Exposure`)
1276 ``outputExposure``
1277 An alias for ``exposure``. (`lsst.afw.image.Exposure`)
1278 ``ossThumb``
1279 Thumbnail image of the exposure after overscan subtraction.
1280 (`numpy.ndarray`)
1281 ``flattenedThumb``
1282 Thumbnail image of the exposure after flat-field correction.
1283 (`numpy.ndarray`)
1284 ``outputStatistics``
1285 Values of the additional statistics calculated.
1286
1287 Raises
1288 ------
1289 RuntimeError
1290 Raised if a configuration option is set to `True`, but the
1291 required calibration data has not been specified.
1292
1293 Notes
1294 -----
1295 The current processed exposure can be viewed by setting the
1296 appropriate `lsstDebug` entries in the ``debug.display``
1297 dictionary. The names of these entries correspond to some of
1298 the `IsrTaskConfig` Boolean options, with the value denoting the
1299 frame to use. The exposure is shown inside the matching
1300 option check and after the processing of that step has
1301 finished. The steps with debug points are:
1302
1303 * doAssembleCcd
1304 * doBias
1305 * doCrosstalk
1306 * doBrighterFatter
1307 * doDark
1308 * doFringe
1309 * doStrayLight
1310 * doFlat
1311
1312 In addition, setting the ``postISRCCD`` entry displays the
1313 exposure after all ISR processing has finished.
1314 """
1315
1316 ccdExposure = self.ensureExposure(ccdExposure, camera, detectorNum)
1317 bias = self.ensureExposure(bias, camera, detectorNum)
1318 dark = self.ensureExposure(dark, camera, detectorNum)
1319 flat = self.ensureExposure(flat, camera, detectorNum)
1320
1321 ccd = ccdExposure.getDetector()
1322 filterLabel = ccdExposure.getFilter()
1323 physicalFilter = isrFunctions.getPhysicalFilter(filterLabel, self.log)
1324
1325 if not ccd:
1326 assert not self.config.doAssembleCcd, "You need a Detector to run assembleCcd."
1327 ccd = [FakeAmp(ccdExposure, self.config)]
1328
1329 # Validate Input
1330 if self.config.doBias and bias is None:
1331 raise RuntimeError("Must supply a bias exposure if config.doBias=True.")
1332 if self.doLinearize(ccd) and linearizer is None:
1333 raise RuntimeError("Must supply a linearizer if config.doLinearize=True for this detector.")
1334 if self.config.doBrighterFatter and bfKernel is None:
1335 raise RuntimeError("Must supply a kernel if config.doBrighterFatter=True.")
1336 if self.config.doDark and dark is None:
1337 raise RuntimeError("Must supply a dark exposure if config.doDark=True.")
1338 if self.config.doFlat and flat is None:
1339 raise RuntimeError("Must supply a flat exposure if config.doFlat=True.")
1340 if self.config.doDefect and defects is None:
1341 raise RuntimeError("Must supply defects if config.doDefect=True.")
1342 if (self.config.doFringe and physicalFilter in self.fringe.config.filters
1343 and fringes.fringes is None):
1344 # The `fringes` object needs to be a pipeBase.Struct, as
1345 # we use it as a `dict` for the parameters of
1346 # `FringeTask.run()`. The `fringes.fringes` `list` may
1347 # not be `None` if `doFringe=True`. Otherwise, raise.
1348 raise RuntimeError("Must supply fringe exposure as a pipeBase.Struct.")
1349 if (self.config.doIlluminationCorrection and physicalFilter in self.config.illumFilters
1350 and illumMaskedImage is None):
1351 raise RuntimeError("Must supply an illumcor if config.doIlluminationCorrection=True.")
1352 if (self.config.doDeferredCharge and deferredChargeCalib is None):
1353 raise RuntimeError("Must supply a deferred charge calibration if config.doDeferredCharge=True.")
1354 if (self.config.usePtcGains and ptc is None):
1355 raise RuntimeError("No ptcDataset provided to use PTC gains.")
1356 if (self.config.usePtcReadNoise and ptc is None):
1357 raise RuntimeError("No ptcDataset provided to use PTC read noise.")
1358
1359 # Validate that the inputs match the exposure configuration.
1360 exposureMetadata = ccdExposure.getMetadata()
1361 if self.config.doBias:
1362 self.compareCameraKeywords(exposureMetadata, bias, "bias")
1363 if self.config.doBrighterFatter:
1364 self.compareCameraKeywords(exposureMetadata, bfKernel, "brighter-fatter")
1365 if self.config.doCrosstalk:
1366 self.compareCameraKeywords(exposureMetadata, crosstalk, "crosstalk")
1367 if self.config.doDark:
1368 self.compareCameraKeywords(exposureMetadata, dark, "dark")
1369 if self.config.doDefect:
1370 self.compareCameraKeywords(exposureMetadata, defects, "defects")
1371 if self.config.doDeferredCharge:
1372 self.compareCameraKeywords(exposureMetadata, deferredChargeCalib, "CTI")
1373 if self.config.doFlat:
1374 self.compareCameraKeywords(exposureMetadata, flat, "flat")
1375 if (self.config.doFringe and physicalFilter in self.fringe.config.filters):
1376 self.compareCameraKeywords(exposureMetadata, fringes.fringes, "fringe")
1377 if (self.config.doIlluminationCorrection and physicalFilter in self.config.illumFilters):
1378 self.compareCameraKeywords(exposureMetadata, illumMaskedImage, "illumination")
1379 if self.doLinearize(ccd):
1380 self.compareCameraKeywords(exposureMetadata, linearizer, "linearizer")
1381 if self.config.usePtcGains or self.config.usePtcReadNoise:
1382 self.compareCameraKeywords(exposureMetadata, ptc, "PTC")
1383 if self.config.doStrayLight:
1384 self.compareCameraKeywords(exposureMetadata, strayLightData, "straylight")
1385
1386 # Start in ADU. Update units to electrons when gain is applied:
1387 # updateVariance, applyGains
1388 # Check if needed during/after BFE correction, CTI correction.
1389 exposureMetadata["LSST ISR UNITS"] = "ADU"
1390
1391 # Begin ISR processing.
1392 if self.config.doConvertIntToFloat:
1393 self.log.info("Converting exposure to floating point values.")
1394 ccdExposure = self.convertIntToFloat(ccdExposure)
1395
1396 if self.config.doBias and self.config.doBiasBeforeOverscan:
1397 self.log.info("Applying bias correction.")
1398 isrFunctions.biasCorrection(ccdExposure.getMaskedImage(), bias.getMaskedImage(),
1399 trimToFit=self.config.doTrimToMatchCalib)
1400 self.debugView(ccdExposure, "doBias")
1401
1402 # Amplifier level processing.
1403 overscans = []
1404
1405 if self.config.doOverscan and self.config.overscan.doParallelOverscan:
1406 # This will attempt to mask bleed pixels across all amplifiers.
1407 self.overscan.maskParallelOverscan(ccdExposure, ccd)
1408
1409 for amp in ccd:
1410 # if ccdExposure is one amp,
1411 # check for coverage to prevent performing ops multiple times
1412 if ccdExposure.getBBox().contains(amp.getBBox()):
1413 # Check for fully masked bad amplifiers,
1414 # and generate masks for SUSPECT and SATURATED values.
1415 badAmp = self.maskAmplifier(ccdExposure, amp, defects)
1416
1417 if self.config.doOverscan and not badAmp:
1418 # Overscan correction on amp-by-amp basis.
1419 overscanResults = self.overscanCorrection(ccdExposure, amp)
1420 self.log.debug("Corrected overscan for amplifier %s.", amp.getName())
1421 if overscanResults is not None and \
1422 self.config.qa is not None and self.config.qa.saveStats is True:
1423 if isinstance(overscanResults.overscanMean, float):
1424 # Only serial overscan was run
1425 mean = overscanResults.overscanMean
1426 sigma = overscanResults.overscanSigma
1427 residMean = overscanResults.residualMean
1428 residSigma = overscanResults.residualSigma
1429 else:
1430 # Both serial and parallel overscan were
1431 # run. Only report serial here.
1432 mean = overscanResults.overscanMean[0]
1433 sigma = overscanResults.overscanSigma[0]
1434 residMean = overscanResults.residualMean[0]
1435 residSigma = overscanResults.residualSigma[0]
1436
1437 self.metadata[f"FIT MEDIAN {amp.getName()}"] = mean
1438 self.metadata[f"FIT STDEV {amp.getName()}"] = sigma
1439 self.log.debug(" Overscan stats for amplifer %s: %f +/- %f",
1440 amp.getName(), mean, sigma)
1441
1442 self.metadata[f"RESIDUAL MEDIAN {amp.getName()}"] = residMean
1443 self.metadata[f"RESIDUAL STDEV {amp.getName()}"] = residSigma
1444 self.log.debug(" Overscan stats for amplifer %s after correction: %f +/- %f",
1445 amp.getName(), residMean, residSigma)
1446
1447 ccdExposure.getMetadata().set('OVERSCAN', "Overscan corrected")
1448 else:
1449 if badAmp:
1450 self.log.warning("Amplifier %s is bad.", amp.getName())
1451 overscanResults = None
1452
1453 overscans.append(overscanResults if overscanResults is not None else None)
1454 else:
1455 self.log.info("Skipped OSCAN for %s.", amp.getName())
1456
1457 # Define an effective PTC that will contain the gain and readout
1458 # noise to be used throughout the ISR task.
1459 ptc = self.defineEffectivePtc(ptc, ccd, bfGains, overscans, exposureMetadata)
1460
1461 if self.config.doDeferredCharge:
1462 self.log.info("Applying deferred charge/CTI correction.")
1463 self.deferredChargeCorrection.run(ccdExposure, deferredChargeCalib)
1464 self.debugView(ccdExposure, "doDeferredCharge")
1465
1466 if self.config.doCrosstalk and self.config.doCrosstalkBeforeAssemble:
1467 self.log.info("Applying crosstalk correction.")
1468 self.crosstalk.run(ccdExposure, crosstalk=crosstalk,
1469 crosstalkSources=crosstalkSources, camera=camera)
1470 self.debugView(ccdExposure, "doCrosstalk")
1471
1472 if self.config.doAssembleCcd:
1473 self.log.info("Assembling CCD from amplifiers.")
1474 ccdExposure = self.assembleCcd.assembleCcd(ccdExposure)
1475
1476 if self.config.expectWcs and not ccdExposure.getWcs():
1477 self.log.warning("No WCS found in input exposure.")
1478 self.debugView(ccdExposure, "doAssembleCcd")
1479
1480 ossThumb = None
1481 if self.config.qa.doThumbnailOss:
1482 ossThumb = isrQa.makeThumbnail(ccdExposure, isrQaConfig=self.config.qa)
1483
1484 if self.config.doBias and not self.config.doBiasBeforeOverscan:
1485 self.log.info("Applying bias correction.")
1486 isrFunctions.biasCorrection(ccdExposure.getMaskedImage(), bias.getMaskedImage(),
1487 trimToFit=self.config.doTrimToMatchCalib)
1488 self.debugView(ccdExposure, "doBias")
1489
1490 if self.config.doVariance:
1491 for amp in ccd:
1492 if ccdExposure.getBBox().contains(amp.getBBox()):
1493 self.log.debug("Constructing variance map for amplifer %s.", amp.getName())
1494 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1495 self.updateVariance(ampExposure, amp, ptc)
1496
1497 if self.config.qa is not None and self.config.qa.saveStats is True:
1498 qaStats = afwMath.makeStatistics(ampExposure.getVariance(),
1499 afwMath.MEDIAN | afwMath.STDEVCLIP)
1500 self.metadata[f"ISR VARIANCE {amp.getName()} MEDIAN"] = \
1501 qaStats.getValue(afwMath.MEDIAN)
1502 self.metadata[f"ISR VARIANCE {amp.getName()} STDEV"] = \
1503 qaStats.getValue(afwMath.STDEVCLIP)
1504 self.log.debug(" Variance stats for amplifer %s: %f +/- %f.",
1505 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1506 qaStats.getValue(afwMath.STDEVCLIP))
1507 if self.config.maskNegativeVariance:
1508 self.maskNegativeVariance(ccdExposure)
1509
1510 if self.doLinearize(ccd):
1511 self.log.info("Applying linearizer.")
1512 linearizer.applyLinearity(image=ccdExposure.getMaskedImage().getImage(),
1513 detector=ccd, log=self.log)
1514
1515 if self.config.doCrosstalk and not self.config.doCrosstalkBeforeAssemble:
1516 self.log.info("Applying crosstalk correction.")
1517 self.crosstalk.run(ccdExposure, crosstalk=crosstalk,
1518 crosstalkSources=crosstalkSources, isTrimmed=True)
1519 self.debugView(ccdExposure, "doCrosstalk")
1520
1521 # Masking block. Optionally mask known defects, NAN/inf pixels,
1522 # widen trails, and do anything else the camera needs. Saturated and
1523 # suspect pixels have already been masked.
1524 if self.config.doDefect:
1525 self.log.info("Masking defects.")
1526 self.maskDefect(ccdExposure, defects)
1527
1528 if self.config.numEdgeSuspect > 0:
1529 self.log.info("Masking edges as SUSPECT.")
1530 self.maskEdges(ccdExposure, numEdgePixels=self.config.numEdgeSuspect,
1531 maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
1532
1533 if self.config.doNanMasking:
1534 self.log.info("Masking non-finite (NAN, inf) value pixels.")
1535 self.maskNan(ccdExposure)
1536
1537 if self.config.doWidenSaturationTrails:
1538 self.log.info("Widening saturation trails.")
1539 isrFunctions.widenSaturationTrails(ccdExposure.getMaskedImage().getMask())
1540
1541 if self.config.doCameraSpecificMasking:
1542 self.log.info("Masking regions for camera specific reasons.")
1543 self.masking.run(ccdExposure)
1544
1545 if self.config.doBrighterFatter:
1546 # We need to apply flats and darks before we can interpolate, and
1547 # we need to interpolate before we do B-F, but we do B-F without
1548 # the flats and darks applied so we can work in units of electrons
1549 # or holes. This context manager applies and then removes the darks
1550 # and flats.
1551 #
1552 # We also do not want to interpolate values here, so operate on
1553 # temporary images so we can apply only the BF-correction and roll
1554 # back the interpolation.
1555 interpExp = ccdExposure.clone()
1556 with self.flatContext(interpExp, flat, dark):
1557 isrFunctions.interpolateFromMask(
1558 maskedImage=interpExp.getMaskedImage(),
1559 fwhm=self.config.fwhm,
1560 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1561 maskNameList=list(self.config.brighterFatterMaskListToInterpolate)
1562 )
1563 bfExp = interpExp.clone()
1564
1565 self.log.info("Applying brighter-fatter correction using kernel type %s / gains %s.",
1566 type(bfKernel), type(bfGains))
1567 if self.config.doFluxConservingBrighterFatterCorrection:
1568 bfResults = isrFunctions.fluxConservingBrighterFatterCorrection(
1569 bfExp,
1570 bfKernel,
1571 self.config.brighterFatterMaxIter,
1572 self.config.brighterFatterThreshold,
1573 self.config.brighterFatterApplyGain,
1574 bfGains
1575 )
1576 else:
1577 bfResults = isrFunctions.brighterFatterCorrection(
1578 bfExp,
1579 bfKernel,
1580 self.config.brighterFatterMaxIter,
1581 self.config.brighterFatterThreshold,
1582 self.config.brighterFatterApplyGain,
1583 bfGains
1584 )
1585 if bfResults[1] == self.config.brighterFatterMaxIter - 1:
1586 self.log.warning("Brighter-fatter correction did not converge, final difference %f.",
1587 bfResults[0])
1588 else:
1589 self.log.info("Finished brighter-fatter correction in %d iterations.",
1590 bfResults[1])
1591 image = ccdExposure.getMaskedImage().getImage()
1592 bfCorr = bfExp.getMaskedImage().getImage()
1593 bfCorr -= interpExp.getMaskedImage().getImage()
1594 image += bfCorr
1595
1596 # Applying the brighter-fatter correction applies a
1597 # convolution to the science image. At the edges this
1598 # convolution may not have sufficient valid pixels to
1599 # produce a valid correction. Mark pixels within the size
1600 # of the brighter-fatter kernel as EDGE to warn of this
1601 # fact.
1602 self.log.info("Ensuring image edges are masked as EDGE to the brighter-fatter kernel size.")
1603 self.maskEdges(ccdExposure, numEdgePixels=numpy.max(bfKernel.shape) // 2,
1604 maskPlane="EDGE")
1605
1606 if self.config.brighterFatterMaskGrowSize > 0:
1607 self.log.info("Growing masks to account for brighter-fatter kernel convolution.")
1608 for maskPlane in self.config.brighterFatterMaskListToInterpolate:
1609 isrFunctions.growMasks(ccdExposure.getMask(),
1610 radius=self.config.brighterFatterMaskGrowSize,
1611 maskNameList=maskPlane,
1612 maskValue=maskPlane)
1613
1614 self.debugView(ccdExposure, "doBrighterFatter")
1615
1616 if self.config.doDark:
1617 self.log.info("Applying dark correction.")
1618 self.darkCorrection(ccdExposure, dark)
1619 self.debugView(ccdExposure, "doDark")
1620
1621 if self.config.doFringe and not self.config.fringeAfterFlat:
1622 self.log.info("Applying fringe correction before flat.")
1623 self.fringe.run(ccdExposure, **fringes.getDict())
1624 self.debugView(ccdExposure, "doFringe")
1625
1626 if self.config.doStrayLight and self.strayLight.check(ccdExposure):
1627 self.log.info("Checking strayLight correction.")
1628 self.strayLight.run(ccdExposure, strayLightData)
1629 self.debugView(ccdExposure, "doStrayLight")
1630
1631 if self.config.doFlat:
1632 self.log.info("Applying flat correction.")
1633 self.flatCorrection(ccdExposure, flat)
1634 self.debugView(ccdExposure, "doFlat")
1635
1636 if self.config.doApplyGains:
1637 self.log.info("Applying gain correction instead of flat.")
1638 isrFunctions.applyGains(ccdExposure, self.config.normalizeGains,
1639 ptcGains=ptc.gain)
1640 exposureMetadata["LSST ISR UNITS"] = "electrons"
1641
1642 if self.config.doFringe and self.config.fringeAfterFlat:
1643 self.log.info("Applying fringe correction after flat.")
1644 self.fringe.run(ccdExposure, **fringes.getDict())
1645
1646 if self.config.doVignette:
1647 if self.config.doMaskVignettePolygon:
1648 self.log.info("Constructing, attaching, and masking vignette polygon.")
1649 else:
1650 self.log.info("Constructing and attaching vignette polygon.")
1651 self.vignettePolygon = self.vignette.run(
1652 exposure=ccdExposure, doUpdateMask=self.config.doMaskVignettePolygon,
1653 vignetteValue=self.config.vignetteValue, log=self.log)
1654
1655 if self.config.doAttachTransmissionCurve:
1656 self.log.info("Adding transmission curves.")
1657 isrFunctions.attachTransmissionCurve(ccdExposure, opticsTransmission=opticsTransmission,
1658 filterTransmission=filterTransmission,
1659 sensorTransmission=sensorTransmission,
1660 atmosphereTransmission=atmosphereTransmission)
1661
1662 flattenedThumb = None
1663 if self.config.qa.doThumbnailFlattened:
1664 flattenedThumb = isrQa.makeThumbnail(ccdExposure, isrQaConfig=self.config.qa)
1665
1666 if self.config.doIlluminationCorrection and physicalFilter in self.config.illumFilters:
1667 self.log.info("Performing illumination correction.")
1668 isrFunctions.illuminationCorrection(ccdExposure.getMaskedImage(),
1669 illumMaskedImage, illumScale=self.config.illumScale,
1670 trimToFit=self.config.doTrimToMatchCalib)
1671
1672 preInterpExp = None
1673 if self.config.doSaveInterpPixels:
1674 preInterpExp = ccdExposure.clone()
1675
1676 # Reset and interpolate bad pixels.
1677 #
1678 # Large contiguous bad regions (which should have the BAD mask
1679 # bit set) should have their values set to the image median.
1680 # This group should include defects and bad amplifiers. As the
1681 # area covered by these defects are large, there's little
1682 # reason to expect that interpolation would provide a more
1683 # useful value.
1684 #
1685 # Smaller defects can be safely interpolated after the larger
1686 # regions have had their pixel values reset. This ensures
1687 # that the remaining defects adjacent to bad amplifiers (as an
1688 # example) do not attempt to interpolate extreme values.
1689 if self.config.doSetBadRegions:
1690 badPixelCount, badPixelValue = isrFunctions.setBadRegions(ccdExposure)
1691 if badPixelCount > 0:
1692 self.log.info("Set %d BAD pixels to %f.", badPixelCount, badPixelValue)
1693
1694 if self.config.doInterpolate:
1695 self.log.info("Interpolating masked pixels.")
1696 isrFunctions.interpolateFromMask(
1697 maskedImage=ccdExposure.getMaskedImage(),
1698 fwhm=self.config.fwhm,
1699 growSaturatedFootprints=self.config.growSaturationFootprintSize,
1700 maskNameList=list(self.config.maskListToInterpolate)
1701 )
1702
1703 self.roughZeroPoint(ccdExposure)
1704
1705 # correct for amp offsets within the CCD
1706 if self.config.doAmpOffset:
1707 self.log.info("Correcting amp offsets.")
1708 self.ampOffset.run(ccdExposure)
1709
1710 if self.config.doMeasureBackground:
1711 self.log.info("Measuring background level.")
1712 self.measureBackground(ccdExposure, self.config.qa)
1713
1714 if self.config.qa is not None and self.config.qa.saveStats is True:
1715 for amp in ccd:
1716 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1717 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
1718 afwMath.MEDIAN | afwMath.STDEVCLIP)
1719 self.metadata[f"ISR BACKGROUND {amp.getName()} MEDIAN"] = qaStats.getValue(afwMath.MEDIAN)
1720 self.metadata[f"ISR BACKGROUND {amp.getName()} STDEV"] = \
1721 qaStats.getValue(afwMath.STDEVCLIP)
1722 self.log.debug(" Background stats for amplifer %s: %f +/- %f",
1723 amp.getName(), qaStats.getValue(afwMath.MEDIAN),
1724 qaStats.getValue(afwMath.STDEVCLIP))
1725
1726 # Calculate standard image quality statistics
1727 if self.config.doStandardStatistics:
1728 metadata = ccdExposure.getMetadata()
1729 for amp in ccd:
1730 ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
1731 ampName = amp.getName()
1732 metadata[f"LSST ISR MASK SAT {ampName}"] = isrFunctions.countMaskedPixels(
1733 ampExposure.getMaskedImage(),
1734 [self.config.saturatedMaskName]
1735 )
1736 metadata[f"LSST ISR MASK BAD {ampName}"] = isrFunctions.countMaskedPixels(
1737 ampExposure.getMaskedImage(),
1738 ["BAD"]
1739 )
1740 qaStats = afwMath.makeStatistics(ampExposure.getImage(),
1741 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP)
1742
1743 metadata[f"LSST ISR FINAL MEAN {ampName}"] = qaStats.getValue(afwMath.MEAN)
1744 metadata[f"LSST ISR FINAL MEDIAN {ampName}"] = qaStats.getValue(afwMath.MEDIAN)
1745 metadata[f"LSST ISR FINAL STDEV {ampName}"] = qaStats.getValue(afwMath.STDEVCLIP)
1746
1747 k1 = f"LSST ISR FINAL MEDIAN {ampName}"
1748 k2 = f"LSST ISR OVERSCAN SERIAL MEDIAN {ampName}"
1749 if self.config.doOverscan and k1 in metadata and k2 in metadata:
1750 metadata[f"LSST ISR LEVEL {ampName}"] = metadata[k1] - metadata[k2]
1751 else:
1752 metadata[f"LSST ISR LEVEL {ampName}"] = numpy.nan
1753
1754 # calculate additional statistics.
1755 outputStatistics = None
1756 if self.config.doCalculateStatistics:
1757 outputStatistics = self.isrStats.run(ccdExposure, overscanResults=overscans,
1758 bias=bias, dark=dark, flat=flat, ptc=ptc,
1759 defects=defects).results
1760
1761 # do any binning.
1762 outputBin1Exposure = None
1763 outputBin2Exposure = None
1764 if self.config.doBinnedExposures:
1765 outputBin1Exposure, outputBin2Exposure = self.makeBinnedImages(ccdExposure)
1766
1767 self.debugView(ccdExposure, "postISRCCD")
1768
1769 return pipeBase.Struct(
1770 exposure=ccdExposure,
1771 ossThumb=ossThumb,
1772 flattenedThumb=flattenedThumb,
1773
1774 outputBin1Exposure=outputBin1Exposure,
1775 outputBin2Exposure=outputBin2Exposure,
1776
1777 preInterpExposure=preInterpExp,
1778 outputExposure=ccdExposure,
1779 outputOssThumbnail=ossThumb,
1780 outputFlattenedThumbnail=flattenedThumb,
1781 outputStatistics=outputStatistics,
1782 )
1783
1784 def defineEffectivePtc(self, ptcDataset, detector, bfGains, overScans, metadata):
1785 """Define an effective Photon Transfer Curve dataset
1786 with nominal gains and noise.
1787
1788 Parameters
1789 ----------
1790 ptcDataset : `lsst.ip.isr.PhotonTransferCurveDataset`
1791 Input Photon Transfer Curve dataset.
1792 detector : `lsst.afw.cameraGeom.Detector`
1793 Detector object.
1794 bfGains : `dict`
1795 Gains from running the brighter-fatter code.
1796 A dict keyed by amplifier name for the detector
1797 in question.
1798 ovserScans : `list` [`lsst.pipe.base.Struct`]
1799 List of overscanResults structures
1800 metadata : `lsst.daf.base.PropertyList`
1801 Exposure metadata to update gain and noise provenance.
1802
1803 Returns
1804 -------
1805 effectivePtc : `lsst.ip.isr.PhotonTransferCurveDataset`
1806 PTC dataset containing gains and readout noise
1807 values to be used throughout
1808 Instrument Signature Removal.
1809 """
1810 amps = detector.getAmplifiers()
1811 ampNames = [amp.getName() for amp in amps]
1812 detName = detector.getName()
1813 effectivePtc = PhotonTransferCurveDataset(ampNames, 'EFFECTIVE_PTC', 1)
1814 boolGainMismatch = False
1815 doWarningPtcValidation = True
1816
1817 for amp, overscanResults in zip(amps, overScans):
1818 ampName = amp.getName()
1819 # Gain:
1820 # Try first with the PTC gains.
1821 gainProvenanceString = "amp"
1822 if self.config.usePtcGains:
1823 gain = ptcDataset.gain[ampName]
1824 gainProvenanceString = "ptc"
1825 self.log.debug("Using gain from Photon Transfer Curve.")
1826 else:
1827 # Try then with the amplifier gain.
1828 # We already have a detector at this point. If there was no
1829 # detector to begin with, one would have been created with
1830 # self.config.gain and self.config.noise. Same comment
1831 # applies for the noise block below.
1832 gain = amp.getGain()
1833
1834 # Check if the gain up to this point differs from the
1835 # gain in bfGains. If so, raise or warn, accordingly.
1836 if not boolGainMismatch and bfGains is not None and ampName in bfGains:
1837 bfGain = bfGains[ampName]
1838 if not math.isclose(gain, bfGain, rel_tol=1e-4):
1839 if self.config.doRaiseOnCalibMismatch:
1840 raise RuntimeError("Gain mismatch for det %s amp %s: "
1841 "(gain (%s): %s, bfGain: %s)",
1842 detName, ampName, gainProvenanceString,
1843 gain, bfGain)
1844 else:
1845 self.log.warning("Gain mismatch for det %s amp %s: "
1846 "(gain (%s): %s, bfGain: %s)",
1847 detName, ampName, gainProvenanceString,
1848 gain, bfGain)
1849 boolGainMismatch = True
1850
1851 # Noise:
1852 # Try first with the empirical noise from the overscan.
1853 noiseProvenanceString = "amp"
1854 if self.config.doEmpiricalReadNoise and overscanResults is not None:
1855 noiseProvenanceString = "serial overscan"
1856 if isinstance(overscanResults.residualSigma, float):
1857 # Only serial overscan was run
1858 noise = overscanResults.residualSigma
1859 else:
1860 # Both serial and parallel overscan were
1861 # run. Only report noise from serial here.
1862 noise = overscanResults.residualSigma[0]
1863 elif self.config.usePtcReadNoise:
1864 # Try then with the PTC noise.
1865 noise = ptcDataset.noise[ampName]
1866 noiseProvenanceString = "ptc"
1867 self.log.debug("Using noise from Photon Transfer Curve.")
1868 else:
1869 # Finally, try with the amplifier noise.
1870 # We already have a detector at this point. If there
1871 # was no detector to begin with, one would have
1872 # been created with self.config.gain and
1873 # self.config.noise.
1874 noise = amp.getReadNoise()
1875
1876 if math.isnan(gain):
1877 gain = 1.0
1878 self.log.warning("Gain for amp %s set to NAN! Updating to"
1879 " 1.0 to generate Poisson variance.", ampName)
1880 elif gain <= 0:
1881 patchedGain = 1.0
1882 self.log.warning("Gain for amp %s == %g <= 0; setting to %f.",
1883 ampName, gain, patchedGain)
1884 gain = patchedGain
1885
1886 # PTC Turnoff:
1887 # Copy it over from the input PTC if it's positive. If it's a nan
1888 # set it to a high value.
1889 if ptcDataset is not None:
1890 ptcTurnoff = ptcDataset.ptcTurnoff[ampName]
1891 else:
1892 ptcTurnoff = 2e19
1893
1894 if (isinstance(ptcTurnoff, numbers.Real) and ptcTurnoff > 0):
1895 effectivePtc.ptcTurnoff[ampName] = ptcTurnoff
1896 elif math.isnan(ptcTurnoff):
1897 effectivePtc.ptcTurnoff[ampName] = 2e19
1898
1899 effectivePtc.gain[ampName] = gain
1900 effectivePtc.noise[ampName] = noise
1901 # Make sure noise,turnoff, and gain make sense
1902 effectivePtc.validateGainNoiseTurnoffValues(ampName, doWarn=doWarningPtcValidation)
1903 doWarningPtcValidation = False
1904
1905 metadata[f"LSST GAIN {amp.getName()}"] = effectivePtc.gain[ampName]
1906 metadata[f"LSST READNOISE {amp.getName()}"] = effectivePtc.noise[ampName]
1907
1908 self.log.info("Det: %s - Noise provenance: %s, Gain provenance: %s",
1909 detName,
1910 noiseProvenanceString,
1911 gainProvenanceString)
1912 metadata["LSST ISR GAIN SOURCE"] = gainProvenanceString
1913 metadata["LSST ISR NOISE SOURCE"] = noiseProvenanceString
1914
1915 return effectivePtc
1916
1917 def ensureExposure(self, inputExp, camera=None, detectorNum=None):
1918 """Ensure that the data returned by Butler is a fully constructed exp.
1919
1920 ISR requires exposure-level image data for historical reasons, so if we
1921 did not recieve that from Butler, construct it from what we have,
1922 modifying the input in place.
1923
1924 Parameters
1925 ----------
1926 inputExp : `lsst.afw.image` image-type.
1927 The input data structure obtained from Butler.
1928 Can be `lsst.afw.image.Exposure`,
1929 `lsst.afw.image.DecoratedImageU`,
1930 or `lsst.afw.image.ImageF`
1931 camera : `lsst.afw.cameraGeom.camera`, optional
1932 The camera associated with the image. Used to find the appropriate
1933 detector if detector is not already set.
1934 detectorNum : `int`, optional
1935 The detector in the camera to attach, if the detector is not
1936 already set.
1937
1938 Returns
1939 -------
1940 inputExp : `lsst.afw.image.Exposure`
1941 The re-constructed exposure, with appropriate detector parameters.
1942
1943 Raises
1944 ------
1945 TypeError
1946 Raised if the input data cannot be used to construct an exposure.
1947 """
1948 if isinstance(inputExp, afwImage.DecoratedImageU):
1950 elif isinstance(inputExp, afwImage.ImageF):
1952 elif isinstance(inputExp, afwImage.MaskedImageF):
1953 inputExp = afwImage.makeExposure(inputExp)
1954 elif isinstance(inputExp, afwImage.Exposure):
1955 pass
1956 elif inputExp is None:
1957 # Assume this will be caught by the setup if it is a problem.
1958 return inputExp
1959 else:
1960 raise TypeError("Input Exposure is not known type in isrTask.ensureExposure: %s." %
1961 (type(inputExp), ))
1962
1963 if inputExp.getDetector() is None:
1964 if camera is None or detectorNum is None:
1965 raise RuntimeError('Must supply both a camera and detector number when using exposures '
1966 'without a detector set.')
1967 inputExp.setDetector(camera[detectorNum])
1968
1969 return inputExp
1970
1971 @staticmethod
1973 """Extract common calibration metadata values that will be written to
1974 output header.
1975
1976 Parameters
1977 ----------
1978 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
1979 Calibration to pull date information from.
1980
1981 Returns
1982 -------
1983 dateString : `str`
1984 Calibration creation date string to add to header.
1985 """
1986 if hasattr(calib, "getMetadata"):
1987 if 'CALIB_CREATION_DATE' in calib.getMetadata():
1988 return " ".join((calib.getMetadata().get("CALIB_CREATION_DATE", "Unknown"),
1989 calib.getMetadata().get("CALIB_CREATION_TIME", "Unknown")))
1990 else:
1991 return " ".join((calib.getMetadata().get("CALIB_CREATE_DATE", "Unknown"),
1992 calib.getMetadata().get("CALIB_CREATE_TIME", "Unknown")))
1993 else:
1994 return "Unknown Unknown"
1995
1996 def compareCameraKeywords(self, exposureMetadata, calib, calibName):
1997 """Compare header keywords to confirm camera states match.
1998
1999 Parameters
2000 ----------
2001 exposureMetadata : `lsst.daf.base.PropertySet`
2002 Header for the exposure being processed.
2003 calib : `lsst.afw.image.Exposure` or `lsst.ip.isr.IsrCalib`
2004 Calibration to be applied.
2005 calibName : `str`
2006 Calib type for log message.
2007 """
2008 try:
2009 calibMetadata = calib.getMetadata()
2010 except AttributeError:
2011 return
2012 for keyword in self.config.cameraKeywordsToCompare:
2013 if keyword in exposureMetadata and keyword in calibMetadata:
2014 if exposureMetadata[keyword] != calibMetadata[keyword]:
2015 if self.config.doRaiseOnCalibMismatch:
2016 raise RuntimeError("Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
2017 calibName, keyword,
2018 exposureMetadata[keyword], calibMetadata[keyword])
2019 else:
2020 self.log.warning("Sequencer mismatch for %s [%s]: exposure: %s calib: %s",
2021 calibName, keyword,
2022 exposureMetadata[keyword], calibMetadata[keyword])
2023 else:
2024 self.log.debug("Sequencer keyword %s not found.", keyword)
2025
2026 def convertIntToFloat(self, exposure):
2027 """Convert exposure image from uint16 to float.
2028
2029 If the exposure does not need to be converted, the input is
2030 immediately returned. For exposures that are converted to use
2031 floating point pixels, the variance is set to unity and the
2032 mask to zero.
2033
2034 Parameters
2035 ----------
2036 exposure : `lsst.afw.image.Exposure`
2037 The raw exposure to be converted.
2038
2039 Returns
2040 -------
2041 newexposure : `lsst.afw.image.Exposure`
2042 The input ``exposure``, converted to floating point pixels.
2043
2044 Raises
2045 ------
2046 RuntimeError
2047 Raised if the exposure type cannot be converted to float.
2048
2049 """
2050 if isinstance(exposure, afwImage.ExposureF):
2051 # Nothing to be done
2052 self.log.debug("Exposure already of type float.")
2053 return exposure
2054 if not hasattr(exposure, "convertF"):
2055 raise RuntimeError("Unable to convert exposure (%s) to float." % type(exposure))
2056
2057 newexposure = exposure.convertF()
2058 newexposure.variance[:] = 1
2059 newexposure.mask[:] = 0x0
2060
2061 return newexposure
2062
2063 def maskAmplifier(self, ccdExposure, amp, defects):
2064 """Identify bad amplifiers, saturated and suspect pixels.
2065
2066 Parameters
2067 ----------
2068 ccdExposure : `lsst.afw.image.Exposure`
2069 Input exposure to be masked.
2070 amp : `lsst.afw.cameraGeom.Amplifier`
2071 Catalog of parameters defining the amplifier on this
2072 exposure to mask.
2073 defects : `lsst.ip.isr.Defects`
2074 List of defects. Used to determine if the entire
2075 amplifier is bad.
2076
2077 Returns
2078 -------
2079 badAmp : `Bool`
2080 If this is true, the entire amplifier area is covered by
2081 defects and unusable.
2082
2083 """
2084 maskedImage = ccdExposure.getMaskedImage()
2085
2086 badAmp = False
2087
2088 # Check if entire amp region is defined as a defect
2089 # NB: need to use amp.getBBox() for correct comparison with current
2090 # defects definition.
2091 if defects is not None:
2092 badAmp = bool(sum([v.getBBox().contains(amp.getBBox()) for v in defects]))
2093
2094 # In the case of a bad amp, we will set mask to "BAD"
2095 # (here use amp.getRawBBox() for correct association with pixels in
2096 # current ccdExposure).
2097 if badAmp:
2098 dataView = afwImage.MaskedImageF(maskedImage, amp.getRawBBox(),
2099 afwImage.PARENT)
2100 maskView = dataView.getMask()
2101 maskView |= maskView.getPlaneBitMask("BAD")
2102 del maskView
2103 return badAmp
2104
2105 # Mask remaining defects after assembleCcd() to allow for defects that
2106 # cross amplifier boundaries. Saturation and suspect pixels can be
2107 # masked now, though.
2108 limits = dict()
2109 if self.config.doSaturation and not badAmp:
2110 # Set to the default from the camera model.
2111 limits.update({self.config.saturatedMaskName: amp.getSaturation()})
2112 # And update if it is set in the config.
2113 if math.isfinite(self.config.saturation):
2114 limits.update({self.config.saturatedMaskName: self.config.saturation})
2115 if self.config.doSuspect and not badAmp:
2116 limits.update({self.config.suspectMaskName: amp.getSuspectLevel()})
2117
2118 for maskName, maskThreshold in limits.items():
2119 if not math.isnan(maskThreshold):
2120 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2121 isrFunctions.makeThresholdMask(
2122 maskedImage=dataView,
2123 threshold=maskThreshold,
2124 growFootprints=0,
2125 maskName=maskName
2126 )
2127
2128 # Determine if we've fully masked this amplifier with SUSPECT and
2129 # SAT pixels.
2130 maskView = afwImage.Mask(maskedImage.getMask(), amp.getRawDataBBox(),
2131 afwImage.PARENT)
2132 maskVal = maskView.getPlaneBitMask([self.config.saturatedMaskName,
2133 self.config.suspectMaskName])
2134 if numpy.all(maskView.getArray() & maskVal > 0):
2135 badAmp = True
2136 maskView |= maskView.getPlaneBitMask("BAD")
2137
2138 return badAmp
2139
2140 def overscanCorrection(self, ccdExposure, amp):
2141 """Apply overscan correction in place.
2142
2143 This method does initial pixel rejection of the overscan
2144 region. The overscan can also be optionally segmented to
2145 allow for discontinuous overscan responses to be fit
2146 separately. The actual overscan subtraction is performed by
2147 the `lsst.ip.isr.overscan.OverscanTask`, which is called here
2148 after the amplifier is preprocessed.
2149
2150 Parameters
2151 ----------
2152 ccdExposure : `lsst.afw.image.Exposure`
2153 Exposure to have overscan correction performed.
2154 amp : `lsst.afw.cameraGeom.Amplifer`
2155 The amplifier to consider while correcting the overscan.
2156
2157 Returns
2158 -------
2159 overscanResults : `lsst.pipe.base.Struct`
2160 Result struct with components:
2161
2162 ``imageFit``
2163 Value or fit subtracted from the amplifier image data.
2164 (scalar or `lsst.afw.image.Image`)
2165 ``overscanFit``
2166 Value or fit subtracted from the overscan image data.
2167 (scalar or `lsst.afw.image.Image`)
2168 ``overscanImage``
2169 Image of the overscan region with the overscan
2170 correction applied. This quantity is used to estimate
2171 the amplifier read noise empirically.
2172 (`lsst.afw.image.Image`)
2173 ``edgeMask``
2174 Mask of the suspect pixels. (`lsst.afw.image.Mask`)
2175 ``overscanMean``
2176 Median overscan fit value. (`float`)
2177 ``overscanSigma``
2178 Clipped standard deviation of the overscan after
2179 correction. (`float`)
2180
2181 Raises
2182 ------
2183 RuntimeError
2184 Raised if the ``amp`` does not contain raw pixel information.
2185
2186 See Also
2187 --------
2188 lsst.ip.isr.overscan.OverscanTask
2189 """
2190 if amp.getRawHorizontalOverscanBBox().isEmpty():
2191 self.log.info("ISR_OSCAN: No overscan region. Not performing overscan correction.")
2192 return None
2193
2194 # Perform overscan correction on subregions.
2195 overscanResults = self.overscan.run(ccdExposure, amp)
2196
2197 metadata = ccdExposure.getMetadata()
2198 ampName = amp.getName()
2199
2200 keyBase = "LSST ISR OVERSCAN"
2201 # Updated quantities
2202 if isinstance(overscanResults.overscanMean, float):
2203 # Serial overscan correction only:
2204 metadata[f"{keyBase} SERIAL MEAN {ampName}"] = overscanResults.overscanMean
2205 metadata[f"{keyBase} SERIAL MEDIAN {ampName}"] = overscanResults.overscanMedian
2206 metadata[f"{keyBase} SERIAL STDEV {ampName}"] = overscanResults.overscanSigma
2207
2208 metadata[f"{keyBase} RESIDUAL SERIAL MEAN {ampName}"] = overscanResults.residualMean
2209 metadata[f"{keyBase} RESIDUAL SERIAL MEDIAN {ampName}"] = overscanResults.residualMedian
2210 metadata[f"{keyBase} RESIDUAL SERIAL STDEV {ampName}"] = overscanResults.residualSigma
2211 elif isinstance(overscanResults.overscanMean, tuple):
2212 # Both serial and parallel overscan have run:
2213 metadata[f"{keyBase} SERIAL MEAN {ampName}"] = overscanResults.overscanMean[0]
2214 metadata[f"{keyBase} SERIAL MEDIAN {ampName}"] = overscanResults.overscanMedian[0]
2215 metadata[f"{keyBase} SERIAL STDEV {ampName}"] = overscanResults.overscanSigma[0]
2216
2217 metadata[f"{keyBase} PARALLEL MEAN {ampName}"] = overscanResults.overscanMean[1]
2218 metadata[f"{keyBase} PARALLEL MEDIAN {ampName}"] = overscanResults.overscanMedian[1]
2219 metadata[f"{keyBase} PARALLEL STDEV {ampName}"] = overscanResults.overscanSigma[1]
2220
2221 metadata[f"{keyBase} RESIDUAL SERIAL MEAN {ampName}"] = overscanResults.residualMean[0]
2222 metadata[f"{keyBase} RESIDUAL SERIAL MEDIAN {ampName}"] = overscanResults.residualMedian[0]
2223 metadata[f"{keyBase} RESIDUAL SERIAL STDEV {ampName}"] = overscanResults.residualSigma[0]
2224
2225 metadata[f"{keyBase} RESIDUAL PARALLEL MEAN {ampName}"] = overscanResults.residualMean[1]
2226 metadata[f"{keyBase} RESIDUAL PARALLEL MEDIAN {ampName}"] = overscanResults.residualMedian[1]
2227 metadata[f"{keyBase} RESIDUAL PARALLEL STDEV {ampName}"] = overscanResults.residualSigma[1]
2228 else:
2229 self.log.warning("Unexpected type for overscan values; none added to header.")
2230
2231 return overscanResults
2232
2233 def updateVariance(self, ampExposure, amp, ptcDataset):
2234 """Set the variance plane using the gain and read noise
2235
2236 The read noise is calculated from the ``overscanImage`` if the
2237 ``doEmpiricalReadNoise`` option is set in the configuration; otherwise
2238 the value from the amplifier data is used.
2239
2240 Parameters
2241 ----------
2242 ampExposure : `lsst.afw.image.Exposure`
2243 Exposure to process.
2244 amp : `lsst.afw.cameraGeom.Amplifier` or `FakeAmp`
2245 Amplifier detector data.
2246 ptcDataset : `lsst.ip.isr.PhotonTransferCurveDataset`
2247 Effective PTC dataset containing the gains and read noise.
2248
2249 See also
2250 --------
2251 lsst.ip.isr.isrFunctions.updateVariance
2252 """
2253 ampName = amp.getName()
2254 # At this point, the effective PTC should have
2255 # gain and noise values.
2256 gain = ptcDataset.gain[ampName]
2257 readNoise = ptcDataset.noise[ampName]
2258
2259 isrFunctions.updateVariance(
2260 maskedImage=ampExposure.getMaskedImage(),
2261 gain=gain,
2262 readNoise=readNoise,
2263 )
2264
2265 def maskNegativeVariance(self, exposure):
2266 """Identify and mask pixels with negative variance values.
2267
2268 Parameters
2269 ----------
2270 exposure : `lsst.afw.image.Exposure`
2271 Exposure to process.
2272
2273 See Also
2274 --------
2275 lsst.ip.isr.isrFunctions.updateVariance
2276 """
2277 maskPlane = exposure.getMask().getPlaneBitMask(self.config.negativeVarianceMaskName)
2278 bad = numpy.where(exposure.getVariance().getArray() <= 0.0)
2279 exposure.mask.array[bad] |= maskPlane
2280
2281 def darkCorrection(self, exposure, darkExposure, invert=False):
2282 """Apply dark correction in place.
2283
2284 Parameters
2285 ----------
2286 exposure : `lsst.afw.image.Exposure`
2287 Exposure to process.
2288 darkExposure : `lsst.afw.image.Exposure`
2289 Dark exposure of the same size as ``exposure``.
2290 invert : `Bool`, optional
2291 If True, re-add the dark to an already corrected image.
2292
2293 Raises
2294 ------
2295 RuntimeError
2296 Raised if either ``exposure`` or ``darkExposure`` do not
2297 have their dark time defined.
2298
2299 See Also
2300 --------
2301 lsst.ip.isr.isrFunctions.darkCorrection
2302 """
2303 expScale = exposure.getInfo().getVisitInfo().getDarkTime()
2304 if math.isnan(expScale):
2305 raise RuntimeError("Exposure darktime is NAN.")
2306 if darkExposure.getInfo().getVisitInfo() is not None \
2307 and not math.isnan(darkExposure.getInfo().getVisitInfo().getDarkTime()):
2308 darkScale = darkExposure.getInfo().getVisitInfo().getDarkTime()
2309 else:
2310 # DM-17444: darkExposure.getInfo.getVisitInfo() is None
2311 # so getDarkTime() does not exist.
2312 self.log.warning("darkExposure.getInfo().getVisitInfo() does not exist. Using darkScale = 1.0.")
2313 darkScale = 1.0
2314
2315 isrFunctions.darkCorrection(
2316 maskedImage=exposure.getMaskedImage(),
2317 darkMaskedImage=darkExposure.getMaskedImage(),
2318 expScale=expScale,
2319 darkScale=darkScale,
2320 invert=invert,
2321 trimToFit=self.config.doTrimToMatchCalib
2322 )
2323
2324 def doLinearize(self, detector):
2325 """Check if linearization is needed for the detector cameraGeom.
2326
2327 Checks config.doLinearize and the linearity type of the first
2328 amplifier.
2329
2330 Parameters
2331 ----------
2332 detector : `lsst.afw.cameraGeom.Detector`
2333 Detector to get linearity type from.
2334
2335 Returns
2336 -------
2337 doLinearize : `Bool`
2338 If True, linearization should be performed.
2339 """
2340 return self.config.doLinearize and \
2341 detector.getAmplifiers()[0].getLinearityType() != NullLinearityType
2342
2343 def flatCorrection(self, exposure, flatExposure, invert=False):
2344 """Apply flat correction in place.
2345
2346 Parameters
2347 ----------
2348 exposure : `lsst.afw.image.Exposure`
2349 Exposure to process.
2350 flatExposure : `lsst.afw.image.Exposure`
2351 Flat exposure of the same size as ``exposure``.
2352 invert : `Bool`, optional
2353 If True, unflatten an already flattened image.
2354
2355 See Also
2356 --------
2357 lsst.ip.isr.isrFunctions.flatCorrection
2358 """
2359 isrFunctions.flatCorrection(
2360 maskedImage=exposure.getMaskedImage(),
2361 flatMaskedImage=flatExposure.getMaskedImage(),
2362 scalingType=self.config.flatScalingType,
2363 userScale=self.config.flatUserScale,
2364 invert=invert,
2365 trimToFit=self.config.doTrimToMatchCalib
2366 )
2367
2368 def saturationDetection(self, exposure, amp):
2369 """Detect and mask saturated pixels in config.saturatedMaskName.
2370
2371 Parameters
2372 ----------
2373 exposure : `lsst.afw.image.Exposure`
2374 Exposure to process. Only the amplifier DataSec is processed.
2375 amp : `lsst.afw.cameraGeom.Amplifier`
2376 Amplifier detector data.
2377
2378 See Also
2379 --------
2380 lsst.ip.isr.isrFunctions.makeThresholdMask
2381 """
2382 if not math.isnan(amp.getSaturation()):
2383 maskedImage = exposure.getMaskedImage()
2384 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2385 isrFunctions.makeThresholdMask(
2386 maskedImage=dataView,
2387 threshold=amp.getSaturation(),
2388 growFootprints=0,
2389 maskName=self.config.saturatedMaskName,
2390 )
2391
2392 def saturationInterpolation(self, exposure):
2393 """Interpolate over saturated pixels, in place.
2394
2395 This method should be called after `saturationDetection`, to
2396 ensure that the saturated pixels have been identified in the
2397 SAT mask. It should also be called after `assembleCcd`, since
2398 saturated regions may cross amplifier boundaries.
2399
2400 Parameters
2401 ----------
2402 exposure : `lsst.afw.image.Exposure`
2403 Exposure to process.
2404
2405 See Also
2406 --------
2407 lsst.ip.isr.isrTask.saturationDetection
2408 lsst.ip.isr.isrFunctions.interpolateFromMask
2409 """
2410 isrFunctions.interpolateFromMask(
2411 maskedImage=exposure.getMaskedImage(),
2412 fwhm=self.config.fwhm,
2413 growSaturatedFootprints=self.config.growSaturationFootprintSize,
2414 maskNameList=list(self.config.saturatedMaskName),
2415 )
2416
2417 def suspectDetection(self, exposure, amp):
2418 """Detect and mask suspect pixels in config.suspectMaskName.
2419
2420 Parameters
2421 ----------
2422 exposure : `lsst.afw.image.Exposure`
2423 Exposure to process. Only the amplifier DataSec is processed.
2424 amp : `lsst.afw.cameraGeom.Amplifier`
2425 Amplifier detector data.
2426
2427 See Also
2428 --------
2429 lsst.ip.isr.isrFunctions.makeThresholdMask
2430
2431 Notes
2432 -----
2433 Suspect pixels are pixels whose value is greater than
2434 amp.getSuspectLevel(). This is intended to indicate pixels that may be
2435 affected by unknown systematics; for example if non-linearity
2436 corrections above a certain level are unstable then that would be a
2437 useful value for suspectLevel. A value of `nan` indicates that no such
2438 level exists and no pixels are to be masked as suspicious.
2439 """
2440 suspectLevel = amp.getSuspectLevel()
2441 if math.isnan(suspectLevel):
2442 return
2443
2444 maskedImage = exposure.getMaskedImage()
2445 dataView = maskedImage.Factory(maskedImage, amp.getRawBBox())
2446 isrFunctions.makeThresholdMask(
2447 maskedImage=dataView,
2448 threshold=suspectLevel,
2449 growFootprints=0,
2450 maskName=self.config.suspectMaskName,
2451 )
2452
2453 def maskDefect(self, exposure, defectBaseList):
2454 """Mask defects using mask plane "BAD", in place.
2455
2456 Parameters
2457 ----------
2458 exposure : `lsst.afw.image.Exposure`
2459 Exposure to process.
2460 defectBaseList : defect-type
2461 List of defects to mask. Can be of type `lsst.ip.isr.Defects`
2462 or `list` of `lsst.afw.image.DefectBase`.
2463
2464 Notes
2465 -----
2466 Call this after CCD assembly, since defects may cross amplifier
2467 boundaries.
2468 """
2469 maskedImage = exposure.getMaskedImage()
2470 if not isinstance(defectBaseList, Defects):
2471 # Promotes DefectBase to Defect
2472 defectList = Defects(defectBaseList)
2473 else:
2474 defectList = defectBaseList
2475 defectList.maskPixels(maskedImage, maskName="BAD")
2476
2477 def maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR'):
2478 """Mask edge pixels with applicable mask plane.
2479
2480 Parameters
2481 ----------
2482 exposure : `lsst.afw.image.Exposure`
2483 Exposure to process.
2484 numEdgePixels : `int`, optional
2485 Number of edge pixels to mask.
2486 maskPlane : `str`, optional
2487 Mask plane name to use.
2488 level : `str`, optional
2489 Level at which to mask edges.
2490 """
2491 maskedImage = exposure.getMaskedImage()
2492 maskBitMask = maskedImage.getMask().getPlaneBitMask(maskPlane)
2493
2494 if numEdgePixels > 0:
2495 if level == 'DETECTOR':
2496 boxes = [maskedImage.getBBox()]
2497 elif level == 'AMP':
2498 boxes = [amp.getBBox() for amp in exposure.getDetector()]
2499
2500 for box in boxes:
2501 # This makes a bbox numEdgeSuspect pixels smaller than the
2502 # image on each side
2503 subImage = maskedImage[box]
2504 box.grow(-numEdgePixels)
2505 # Mask pixels outside box
2506 SourceDetectionTask.setEdgeBits(
2507 subImage,
2508 box,
2509 maskBitMask)
2510
2511 def maskAndInterpolateDefects(self, exposure, defectBaseList):
2512 """Mask and interpolate defects using mask plane "BAD", in place.
2513
2514 Parameters
2515 ----------
2516 exposure : `lsst.afw.image.Exposure`
2517 Exposure to process.
2518 defectBaseList : defects-like
2519 List of defects to mask and interpolate. Can be
2520 `lsst.ip.isr.Defects` or `list` of `lsst.afw.image.DefectBase`.
2521
2522 See Also
2523 --------
2524 lsst.ip.isr.isrTask.maskDefect
2525 """
2526 self.maskDefect(exposure, defectBaseList)
2527 self.maskEdges(exposure, numEdgePixels=self.config.numEdgeSuspect,
2528 maskPlane="SUSPECT", level=self.config.edgeMaskLevel)
2529 isrFunctions.interpolateFromMask(
2530 maskedImage=exposure.getMaskedImage(),
2531 fwhm=self.config.fwhm,
2532 growSaturatedFootprints=0,
2533 maskNameList=["BAD"],
2534 )
2535
2536 def maskNan(self, exposure):
2537 """Mask NaNs using mask plane "UNMASKEDNAN", in place.
2538
2539 Parameters
2540 ----------
2541 exposure : `lsst.afw.image.Exposure`
2542 Exposure to process.
2543
2544 Notes
2545 -----
2546 We mask over all non-finite values (NaN, inf), including those
2547 that are masked with other bits (because those may or may not be
2548 interpolated over later, and we want to remove all NaN/infs).
2549 Despite this behaviour, the "UNMASKEDNAN" mask plane is used to
2550 preserve the historical name.
2551 """
2552 maskedImage = exposure.getMaskedImage()
2553
2554 # Find and mask NaNs
2555 maskedImage.getMask().addMaskPlane("UNMASKEDNAN")
2556 maskVal = maskedImage.getMask().getPlaneBitMask("UNMASKEDNAN")
2557 numNans = maskNans(maskedImage, maskVal)
2558 self.metadata["NUMNANS"] = numNans
2559 if numNans > 0:
2560 self.log.warning("There were %d unmasked NaNs.", numNans)
2561
2562 def maskAndInterpolateNan(self, exposure):
2563 """"Mask and interpolate NaN/infs using mask plane "UNMASKEDNAN",
2564 in place.
2565
2566 Parameters
2567 ----------
2568 exposure : `lsst.afw.image.Exposure`
2569 Exposure to process.
2570
2571 See Also
2572 --------
2573 lsst.ip.isr.isrTask.maskNan
2574 """
2575 self.maskNan(exposure)
2576 isrFunctions.interpolateFromMask(
2577 maskedImage=exposure.getMaskedImage(),
2578 fwhm=self.config.fwhm,
2579 growSaturatedFootprints=0,
2580 maskNameList=["UNMASKEDNAN"],
2581 )
2582
2583 def measureBackground(self, exposure, IsrQaConfig=None):
2584 """Measure the image background in subgrids, for quality control.
2585
2586 Parameters
2587 ----------
2588 exposure : `lsst.afw.image.Exposure`
2589 Exposure to process.
2590 IsrQaConfig : `lsst.ip.isr.isrQa.IsrQaConfig`
2591 Configuration object containing parameters on which background
2592 statistics and subgrids to use.
2593 """
2594 if IsrQaConfig is not None:
2595 statsControl = afwMath.StatisticsControl(IsrQaConfig.flatness.clipSigma,
2596 IsrQaConfig.flatness.nIter)
2597 maskVal = exposure.getMaskedImage().getMask().getPlaneBitMask(["BAD", "SAT", "DETECTED"])
2598 statsControl.setAndMask(maskVal)
2599 maskedImage = exposure.getMaskedImage()
2600 stats = afwMath.makeStatistics(maskedImage, afwMath.MEDIAN | afwMath.STDEVCLIP, statsControl)
2601 skyLevel = stats.getValue(afwMath.MEDIAN)
2602 skySigma = stats.getValue(afwMath.STDEVCLIP)
2603 self.log.info("Flattened sky level: %f +/- %f.", skyLevel, skySigma)
2604 metadata = exposure.getMetadata()
2605 metadata["SKYLEVEL"] = skyLevel
2606 metadata["SKYSIGMA"] = skySigma
2607
2608 # calcluating flatlevel over the subgrids
2609 stat = afwMath.MEANCLIP if IsrQaConfig.flatness.doClip else afwMath.MEAN
2610 meshXHalf = int(IsrQaConfig.flatness.meshX/2.)
2611 meshYHalf = int(IsrQaConfig.flatness.meshY/2.)
2612 nX = int((exposure.getWidth() + meshXHalf) / IsrQaConfig.flatness.meshX)
2613 nY = int((exposure.getHeight() + meshYHalf) / IsrQaConfig.flatness.meshY)
2614 skyLevels = numpy.zeros((nX, nY))
2615
2616 for j in range(nY):
2617 yc = meshYHalf + j * IsrQaConfig.flatness.meshY
2618 for i in range(nX):
2619 xc = meshXHalf + i * IsrQaConfig.flatness.meshX
2620
2621 xLLC = xc - meshXHalf
2622 yLLC = yc - meshYHalf
2623 xURC = xc + meshXHalf - 1
2624 yURC = yc + meshYHalf - 1
2625
2626 bbox = lsst.geom.Box2I(lsst.geom.Point2I(xLLC, yLLC), lsst.geom.Point2I(xURC, yURC))
2627 miMesh = maskedImage.Factory(exposure.getMaskedImage(), bbox, afwImage.LOCAL)
2628
2629 skyLevels[i, j] = afwMath.makeStatistics(miMesh, stat, statsControl).getValue()
2630
2631 good = numpy.where(numpy.isfinite(skyLevels))
2632 skyMedian = numpy.median(skyLevels[good])
2633 flatness = (skyLevels[good] - skyMedian) / skyMedian
2634 flatness_rms = numpy.std(flatness)
2635 flatness_pp = flatness.max() - flatness.min() if len(flatness) > 0 else numpy.nan
2636
2637 self.log.info("Measuring sky levels in %dx%d grids: %f.", nX, nY, skyMedian)
2638 self.log.info("Sky flatness in %dx%d grids - pp: %f rms: %f.",
2639 nX, nY, flatness_pp, flatness_rms)
2640
2641 metadata["FLATNESS_PP"] = float(flatness_pp)
2642 metadata["FLATNESS_RMS"] = float(flatness_rms)
2643 metadata["FLATNESS_NGRIDS"] = '%dx%d' % (nX, nY)
2644 metadata["FLATNESS_MESHX"] = IsrQaConfig.flatness.meshX
2645 metadata["FLATNESS_MESHY"] = IsrQaConfig.flatness.meshY
2646
2647 def roughZeroPoint(self, exposure):
2648 """Set an approximate magnitude zero point for the exposure.
2649
2650 Parameters
2651 ----------
2652 exposure : `lsst.afw.image.Exposure`
2653 Exposure to process.
2654 """
2655 filterLabel = exposure.getFilter()
2656 physicalFilter = isrFunctions.getPhysicalFilter(filterLabel, self.log)
2657
2658 if physicalFilter in self.config.fluxMag0T1:
2659 fluxMag0 = self.config.fluxMag0T1[physicalFilter]
2660 else:
2661 self.log.warning("No rough magnitude zero point defined for filter %s.", physicalFilter)
2662 fluxMag0 = self.config.defaultFluxMag0T1
2663
2664 expTime = exposure.getInfo().getVisitInfo().getExposureTime()
2665 if not expTime > 0: # handle NaN as well as <= 0
2666 self.log.warning("Non-positive exposure time; skipping rough zero point.")
2667 return
2668
2669 self.log.info("Setting rough magnitude zero point for filter %s: %f",
2670 physicalFilter, 2.5*math.log10(fluxMag0*expTime))
2671 exposure.setPhotoCalib(afwImage.makePhotoCalibFromCalibZeroPoint(fluxMag0*expTime, 0.0))
2672
2673 @contextmanager
2674 def flatContext(self, exp, flat, dark=None):
2675 """Context manager that applies and removes flats and darks,
2676 if the task is configured to apply them.
2677
2678 Parameters
2679 ----------
2680 exp : `lsst.afw.image.Exposure`
2681 Exposure to process.
2682 flat : `lsst.afw.image.Exposure`
2683 Flat exposure the same size as ``exp``.
2684 dark : `lsst.afw.image.Exposure`, optional
2685 Dark exposure the same size as ``exp``.
2686
2687 Yields
2688 ------
2689 exp : `lsst.afw.image.Exposure`
2690 The flat and dark corrected exposure.
2691 """
2692 if self.config.doDark and dark is not None:
2693 self.darkCorrection(exp, dark)
2694 if self.config.doFlat:
2695 self.flatCorrection(exp, flat)
2696 try:
2697 yield exp
2698 finally:
2699 if self.config.doFlat:
2700 self.flatCorrection(exp, flat, invert=True)
2701 if self.config.doDark and dark is not None:
2702 self.darkCorrection(exp, dark, invert=True)
2703
2704 def makeBinnedImages(self, exposure):
2705 """Make visualizeVisit style binned exposures.
2706
2707 Parameters
2708 ----------
2709 exposure : `lsst.afw.image.Exposure`
2710 Exposure to bin.
2711
2712 Returns
2713 -------
2714 bin1 : `lsst.afw.image.Exposure`
2715 Binned exposure using binFactor1.
2716 bin2 : `lsst.afw.image.Exposure`
2717 Binned exposure using binFactor2.
2718 """
2719 mi = exposure.getMaskedImage()
2720
2721 bin1 = afwMath.binImage(mi, self.config.binFactor1)
2722 bin2 = afwMath.binImage(mi, self.config.binFactor2)
2723
2724 bin1 = afwImage.makeExposure(bin1)
2725 bin2 = afwImage.makeExposure(bin2)
2726
2727 bin1.setInfo(exposure.getInfo())
2728 bin2.setInfo(exposure.getInfo())
2729
2730 return bin1, bin2
2731
2732 def debugView(self, exposure, stepname):
2733 """Utility function to examine ISR exposure at different stages.
2734
2735 Parameters
2736 ----------
2737 exposure : `lsst.afw.image.Exposure`
2738 Exposure to view.
2739 stepname : `str`
2740 State of processing to view.
2741 """
2742 frame = getDebugFrame(self._display, stepname)
2743 if frame:
2744 display = getDisplay(frame)
2745 display.scale('asinh', 'zscale')
2746 display.mtv(exposure)
2747 prompt = "Press Enter to continue [c]... "
2748 while True:
2749 ans = input(prompt).lower()
2750 if ans in ("", "c",):
2751 break
2752
2753
2754class FakeAmp(object):
2755 """A Detector-like object that supports returning gain and saturation level
2756
2757 This is used when the input exposure does not have a detector.
2758
2759 Parameters
2760 ----------
2761 exposure : `lsst.afw.image.Exposure`
2762 Exposure to generate a fake amplifier for.
2763 config : `lsst.ip.isr.isrTaskConfig`
2764 Configuration to apply to the fake amplifier.
2765 """
2766
2767 def __init__(self, exposure, config):
2768 self._bbox = exposure.getBBox(afwImage.LOCAL)
2770 self._gain = config.gain
2771 self._readNoise = config.readNoise
2772 self._saturation = config.saturation
2773
2774 def getBBox(self):
2775 return self._bbox
2776
2777 def getRawBBox(self):
2778 return self._bbox
2779
2783 def getGain(self):
2784 return self._gain
2785
2786 def getReadNoise(self):
2787 return self._readNoise
2788
2789 def getSaturation(self):
2790 return self._saturation
2791
2793 return float("NaN")
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:77
Pass parameters to a Statistics object.
Definition Statistics.h:83
An integer coordinate rectangle.
Definition Box.h:55
__init__(self, exposure, config)
Definition isrTask.py:2767
run(self, ccdExposure, *camera=None, bias=None, linearizer=None, crosstalk=None, crosstalkSources=None, dark=None, flat=None, ptc=None, bfKernel=None, bfGains=None, defects=None, fringes=pipeBase.Struct(fringes=None), opticsTransmission=None, filterTransmission=None, sensorTransmission=None, atmosphereTransmission=None, detectorNum=None, strayLightData=None, illumMaskedImage=None, deferredChargeCalib=None)
Definition isrTask.py:1185
maskEdges(self, exposure, numEdgePixels=0, maskPlane="SUSPECT", level='DETECTOR')
Definition isrTask.py:2477
updateVariance(self, ampExposure, amp, ptcDataset)
Definition isrTask.py:2233
maskAndInterpolateDefects(self, exposure, defectBaseList)
Definition isrTask.py:2511
roughZeroPoint(self, exposure)
Definition isrTask.py:2647
defineEffectivePtc(self, ptcDataset, detector, bfGains, overScans, metadata)
Definition isrTask.py:1784
runQuantum(self, butlerQC, inputRefs, outputRefs)
Definition isrTask.py:1024
ensureExposure(self, inputExp, camera=None, detectorNum=None)
Definition isrTask.py:1917
maskDefect(self, exposure, defectBaseList)
Definition isrTask.py:2453
convertIntToFloat(self, exposure)
Definition isrTask.py:2026
maskAndInterpolateNan(self, exposure)
Definition isrTask.py:2562
flatCorrection(self, exposure, flatExposure, invert=False)
Definition isrTask.py:2343
debugView(self, exposure, stepname)
Definition isrTask.py:2732
measureBackground(self, exposure, IsrQaConfig=None)
Definition isrTask.py:2583
compareCameraKeywords(self, exposureMetadata, calib, calibName)
Definition isrTask.py:1996
maskAmplifier(self, ccdExposure, amp, defects)
Definition isrTask.py:2063
saturationInterpolation(self, exposure)
Definition isrTask.py:2392
saturationDetection(self, exposure, amp)
Definition isrTask.py:2368
doLinearize(self, detector)
Definition isrTask.py:2324
darkCorrection(self, exposure, darkExposure, invert=False)
Definition isrTask.py:2281
suspectDetection(self, exposure, amp)
Definition isrTask.py:2417
overscanCorrection(self, ccdExposure, amp)
Definition isrTask.py:2140
makeBinnedImages(self, exposure)
Definition isrTask.py:2704
flatContext(self, exp, flat, dark=None)
Definition isrTask.py:2674
maskNegativeVariance(self, exposure)
Definition isrTask.py:2265
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT > > image, typename std::shared_ptr< Mask< MaskPixelT > > mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT > > variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
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
crosstalkSourceLookup(datasetType, registry, quantumDataId, collections)
Definition isrTask.py:62