Coverage for python / lsst / ip / isr / isrStatistics.py: 11%
408 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-07 08:28 +0000
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/>.
22__all__ = ["IsrStatisticsTaskConfig", "IsrStatisticsTask"]
24import numpy as np
25import astropy.stats
26from scipy.signal.windows import hamming, hann, gaussian
27from scipy.signal import butter, filtfilt
28from scipy.stats import linregress
30import lsst.afw.math as afwMath
31import lsst.afw.image as afwImage
32import lsst.pipe.base as pipeBase
33import lsst.pex.config as pexConfig
35from lsst.afw.cameraGeom import ReadoutCorner
36from .isrFunctions import gainContext, isTrimmedExposure
39class IsrStatisticsTaskConfig(pexConfig.Config):
40 """Image statistics options.
41 """
42 doCtiStatistics = pexConfig.Field(
43 dtype=bool,
44 doc="Measure CTI statistics from image and overscans?",
45 default=False,
46 )
47 doApplyGainsForCtiStatistics = pexConfig.Field(
48 dtype=bool,
49 doc="Apply gain to the overscan region when measuring CTI statistics?",
50 default=True,
51 # TODO: DM-46721
52 deprecated="This field is no longer used. Will be removed after v28.",
53 )
55 doBandingStatistics = pexConfig.Field(
56 dtype=bool,
57 doc="Measure image banding metric?",
58 default=False,
59 )
60 bandingKernelSize = pexConfig.Field(
61 dtype=int,
62 doc="Width of box for boxcar smoothing for banding metric.",
63 default=3,
64 check=lambda x: x == 0 or x % 2 != 0,
65 )
66 bandingFractionLow = pexConfig.Field(
67 dtype=float,
68 doc="Fraction of values to exclude from low samples.",
69 default=0.1,
70 check=lambda x: x >= 0.0 and x <= 1.0
71 )
72 bandingFractionHigh = pexConfig.Field(
73 dtype=float,
74 doc="Fraction of values to exclude from high samples.",
75 default=0.9,
76 check=lambda x: x >= 0.0 and x <= 1.0,
77 )
78 bandingUseHalfDetector = pexConfig.Field(
79 dtype=float,
80 doc="Use only the first half set of amplifiers.",
81 default=True,
82 )
84 doProjectionStatistics = pexConfig.Field(
85 dtype=bool,
86 doc="Measure projection metric?",
87 default=False,
88 )
89 projectionKernelSize = pexConfig.Field(
90 dtype=int,
91 doc="Width of box for boxcar smoothing of projections.",
92 default=0,
93 check=lambda x: x == 0 or x % 2 != 0,
94 )
95 doProjectionFft = pexConfig.Field(
96 dtype=bool,
97 doc="Generate FFTs from the image projections?",
98 default=False,
99 )
100 projectionFftWindow = pexConfig.ChoiceField(
101 dtype=str,
102 doc="Type of windowing to use prior to calculating FFT.",
103 default="HAMMING",
104 allowed={
105 "HAMMING": "Hamming window.",
106 "HANN": "Hann window.",
107 "GAUSSIAN": "Gaussian window.",
108 "NONE": "No window."
109 }
110 )
112 doDivisaderoStatistics = pexConfig.Field(
113 dtype=bool,
114 doc="Measure divisadero tearing statistics?",
115 default=False,
116 )
117 divisaderoEdgePixels = pexConfig.Field(
118 dtype=int,
119 doc="Number of edge pixels excluded from divisadero linear fit.",
120 default=25,
121 )
122 divisaderoNumImpactPixels = pexConfig.Field(
123 dtype=int,
124 doc="Number of edge pixels to examine for divisadero tearing.",
125 default=2,
126 )
127 divisaderoProjectionMinimum = pexConfig.Field(
128 dtype=int,
129 doc="Minimum row to consider when taking robust mean of columns.",
130 default=10,
131 )
132 divisaderoProjectionMaximum = pexConfig.Field(
133 dtype=int,
134 doc="Maximum row to consider when taking robust mean of columns",
135 default=210,
136 )
138 doCopyCalibDistributionStatistics = pexConfig.Field(
139 dtype=bool,
140 doc="Copy calibration distribution statistics to output?",
141 default=False,
142 )
143 expectedDistributionLevels = pexConfig.ListField(
144 dtype=float,
145 doc="Percentile levels expected in the calibration header.",
146 default=[0, 5, 16, 50, 84, 95, 100],
147 )
149 doBiasShiftStatistics = pexConfig.Field(
150 dtype=bool,
151 doc="Measure number of image shifts in overscan?",
152 default=False,
153 )
154 biasShiftFilterOrder = pexConfig.Field(
155 dtype=int,
156 doc="Filter order for Butterworth highpass filter.",
157 default=5,
158 )
159 biasShiftCutoff = pexConfig.Field(
160 dtype=float,
161 doc="Cutoff frequency for highpass filter.",
162 default=1.0/15.0,
163 )
164 biasShiftWindow = pexConfig.Field(
165 dtype=int,
166 doc="Filter window size in pixels for highpass filter.",
167 default=30,
168 )
169 biasShiftThreshold = pexConfig.Field(
170 dtype=float,
171 doc="S/N threshold for bias shift detection.",
172 default=3.0,
173 )
174 biasShiftRowSkip = pexConfig.Field(
175 dtype=int,
176 doc="Number of rows to skip for the bias shift detection.",
177 default=30,
178 )
179 biasShiftColumnSkip = pexConfig.Field(
180 dtype=int,
181 doc="Number of columns to skip when averaging the overscan region.",
182 default=3,
183 )
185 doAmplifierCorrelationStatistics = pexConfig.Field(
186 dtype=bool,
187 doc="Measure amplifier correlations?",
188 default=False,
189 )
191 stat = pexConfig.Field(
192 dtype=str,
193 default="MEANCLIP",
194 doc="Statistic name to use to measure regions.",
195 )
196 nSigmaClip = pexConfig.Field(
197 dtype=float,
198 default=3.0,
199 doc="Clipping threshold for background",
200 )
201 nIter = pexConfig.Field(
202 dtype=int,
203 default=3,
204 doc="Clipping iterations for background",
205 )
206 badMask = pexConfig.ListField(
207 dtype=str,
208 default=["BAD", "INTRP", "SAT"],
209 doc="Mask planes to ignore when identifying source pixels."
210 )
213class IsrStatisticsTask(pipeBase.Task):
214 """Task to measure arbitrary statistics on ISR processed exposures.
216 The goal is to wrap a number of optional measurements that are
217 useful for calibration production and detector stability.
218 """
219 ConfigClass = IsrStatisticsTaskConfig
220 _DefaultName = "isrStatistics"
222 def __init__(self, statControl=None, **kwargs):
223 super().__init__(**kwargs)
224 self.statControl = afwMath.StatisticsControl(self.config.nSigmaClip, self.config.nIter,
225 afwImage.Mask.getPlaneBitMask(self.config.badMask))
226 self.statType = afwMath.stringToStatisticsProperty(self.config.stat)
228 def run(self, inputExp, untrimmedInputExposure=None, ptc=None, serialOverscanResults=None,
229 parallelOverscanResults=None, doLegacyCtiStatistics=False, **kwargs):
230 """Task to run arbitrary statistics.
232 The statistics should be measured by individual methods, and
233 add to the dictionary in the return struct.
235 Parameters
236 ----------
237 inputExp : `lsst.afw.image.Exposure`
238 The exposure to measure.
239 untrimmedInputExp :
240 The exposure to measure overscan statistics from.
241 ptc : `lsst.ip.isr.PtcDataset`, optional
242 A PTC object containing gains to use.
243 serialOverscanResults : `list` [`lsst.pipe.base.Struct`], optional
244 List of serial overscan results. Expected fields are:
246 ``imageFit``
247 Value or fit subtracted from the amplifier image data
248 (scalar or `lsst.afw.image.Image`).
249 ``overscanFit``
250 Value or fit subtracted from the overscan image data
251 (scalar or `lsst.afw.image.Image`).
252 ``overscanImage``
253 Image of the overscan region with the overscan
254 correction applied (`lsst.afw.image.Image`). This
255 quantity is used to estimate the amplifier read noise
256 empirically.
257 parallelOverscanResults : `list` [`lsst.pipe.base.Struct`], optional
258 List of parallel overscan results. Expected fields are:
260 ``imageFit``
261 Value or fit subtracted from the amplifier image data
262 (scalar or `lsst.afw.image.Image`).
263 ``overscanFit``
264 Value or fit subtracted from the overscan image data
265 (scalar or `lsst.afw.image.Image`).
266 ``overscanImage``
267 Image of the overscan region with the overscan
268 correction applied (`lsst.afw.image.Image`). This
269 quantity is used to estimate the amplifier read noise
270 empirically.
271 doLegacyCtiStatistics : `bool`, optional
272 Use the older version of measureCti (not recommended).
273 This should be True if and only if this task is called
274 from IsrTask. TODO: Deprecate legacy CTI + CTI correction
275 from IsrTask (DM-48757).
276 **kwargs :
277 Keyword arguments. Calibrations being passed in should
278 have an entry here.
280 Returns
281 -------
282 resultStruct : `lsst.pipe.base.Struct`
283 Contains the measured statistics as a dict stored in a
284 field named ``results``.
286 Raises
287 ------
288 RuntimeError
289 Raised if the amplifier gains could not be found.
290 """
291 # Verify inputs
292 if self.config.doCtiStatistics and not doLegacyCtiStatistics:
293 if untrimmedInputExposure is None:
294 raise RuntimeError("Must pass untrimmed exposure if doCtiStatistics==True.")
296 # Find gains.
297 detector = inputExp.getDetector()
298 if ptc is not None:
299 gains = ptc.gain
300 elif detector is not None:
301 gains = {amp.getName(): amp.getGain() for amp in detector.getAmplifiers()}
302 else:
303 raise RuntimeError("No source of gains provided.")
305 ctiResults = None
306 if self.config.doCtiStatistics:
307 if doLegacyCtiStatistics:
308 self.log.warning("Calculating the legacy CTI statistics is not recommended.")
309 ctiResults = self.measureCtiLegacy(inputExp, serialOverscanResults, gains)
310 else:
311 ctiResults = self.measureCti(inputExp, untrimmedInputExposure, gains)
313 bandingResults = None
314 if self.config.doBandingStatistics:
315 bandingResults = self.measureBanding(inputExp, serialOverscanResults)
317 projectionResults = None
318 if self.config.doProjectionStatistics:
319 projectionResults = self.measureProjectionStatistics(inputExp, serialOverscanResults)
321 divisaderoResults = None
322 if self.config.doDivisaderoStatistics:
323 divisaderoResults = self.measureDivisaderoStatistics(inputExp, **kwargs)
325 calibDistributionResults = None
326 if self.config.doCopyCalibDistributionStatistics:
327 calibDistributionResults = self.copyCalibDistributionStatistics(inputExp, **kwargs)
329 biasShiftResults = None
330 if self.config.doBiasShiftStatistics:
331 biasShiftResults = self.measureBiasShifts(inputExp, serialOverscanResults)
333 ampCorrelationResults = None
334 if self.config.doAmplifierCorrelationStatistics:
335 ampCorrelationResults = self.measureAmpCorrelations(inputExp, serialOverscanResults)
337 mjd = inputExp.getMetadata().get("MJD", None)
339 return pipeBase.Struct(
340 results={"CTI": ctiResults,
341 "BANDING": bandingResults,
342 "PROJECTION": projectionResults,
343 "CALIBDIST": calibDistributionResults,
344 "BIASSHIFT": biasShiftResults,
345 "AMPCORR": ampCorrelationResults,
346 "MJD": mjd,
347 'DIVISADERO': divisaderoResults,
348 },
349 )
351 def measureCti(self, inputExp, untrimmedInputExp, gains):
352 """Task to measure CTI statistics.
354 Parameters
355 ----------
356 inputExp : `lsst.afw.image.Exposure`
357 The exposure to measure.
358 untrimmedInputExp : `lsst.afw.image.Exposure`
359 Exposure to measure overscan from.
360 gains : `dict` [`str` `float`]
361 Dictionary of per-amplifier gains, indexed by amplifier name.
363 Returns
364 -------
365 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
366 Dictionary of measurements, keyed by amplifier name and
367 statistics segment. Everything in units based on electron.
369 Notes
370 -------
371 The input exposure is needed because it contains the last imaging
372 pixel, with defects applied. And the untrimmed input exposure is
373 needed because it contains the overscan regions. It needs to be
374 this way because the defect masking code requires that the image
375 be trimmed, but we need the image with defects masked to measure
376 the CTI from the last imaging pixel.
377 """
378 outputStats = {}
380 detector = inputExp.getDetector()
381 image = inputExp.image
382 untrimmedImage = untrimmedInputExp.image
384 # It only makes sense to measure CTI in electron units.
385 # Make it so.
386 imageUnits = inputExp.getMetadata().get("LSST ISR UNITS")
387 untrimmedImageUnits = untrimmedInputExp.getMetadata().get("LSST ISR UNITS")
389 # Ensure we always have the same units.
390 applyGainToImage = True if imageUnits == "adu" else False
391 applyGainToUntrimmedImage = True if untrimmedImageUnits == "adu" else False
393 with gainContext(inputExp, image, applyGainToImage, gains, isTrimmed=True):
394 with gainContext(untrimmedInputExp, untrimmedImage, applyGainToUntrimmedImage,
395 gains, isTrimmed=False):
396 for amp in detector.getAmplifiers():
397 ampStats = {}
398 readoutCorner = amp.getReadoutCorner()
400 ampStats["INPUT_GAIN"] = float(gains[amp.getName()])
401 turnoff = inputExp.metadata[f"LSST ISR PTCTURNOFF {amp.getName()}"]
402 ampStats["INPUT_PTCTURNOFF"] = float(turnoff)
404 # Full data region.
405 dataRegion = image[amp.getBBox()]
406 serialOverscanImage = untrimmedImage[amp.getRawSerialOverscanBBox()]
407 parallelOverscanImage = untrimmedImage[amp.getRawSerialOverscanBBox()]
409 # Get the mean of the image
410 ampStats["IMAGE_MEAN"] = afwMath.makeStatistics(dataRegion, self.statType,
411 self.statControl).getValue()
413 # First and last image columns.
414 colA = afwMath.makeStatistics(
415 dataRegion.array[:, 0],
416 self.statType,
417 self.statControl,
418 ).getValue()
419 colZ = afwMath.makeStatistics(
420 dataRegion.array[:, -1],
421 self.statType,
422 self.statControl,
423 ).getValue()
425 # First and last image rows.
426 rowA = afwMath.makeStatistics(
427 dataRegion.array[0, :],
428 self.statType,
429 self.statControl,
430 ).getValue()
431 rowZ = afwMath.makeStatistics(
432 dataRegion.array[-1, :],
433 self.statType,
434 self.statControl,
435 ).getValue()
437 # We want these relative to the readout corner. If that's
438 # on the right side, we need to swap them.
439 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
440 ampStats["FIRST_COLUMN_MEAN"] = colZ
441 ampStats["LAST_COLUMN_MEAN"] = colA
442 else:
443 ampStats["FIRST_COLUMN_MEAN"] = colA
444 ampStats["LAST_COLUMN_MEAN"] = colZ
446 # We want these relative to the readout corner. If that's
447 # on the top, we need to swap them.
448 if readoutCorner in (ReadoutCorner.UR, ReadoutCorner.UL):
449 ampStats["FIRST_ROW_MEAN"] = rowZ
450 ampStats["LAST_ROW_MEAN"] = rowA
451 else:
452 ampStats["FIRST_ROW_MEAN"] = rowA
453 ampStats["LAST_ROW_MEAN"] = rowZ
455 # Measure the columns of the serial overscan and the rows
456 # of the parallel overscan.
457 nSerialOverscanCols = amp.getRawSerialOverscanBBox().getWidth()
458 nParallelOverscanRows = amp.getRawParallelOverscanBBox().getHeight()
460 # Calculate serial overscan statistics
461 columns = []
462 columnValues = []
463 for idx in range(0, nSerialOverscanCols):
464 # If overscan.doParallelOverscan=True, the
465 # overscanImage will contain both the serial
466 # and parallel overscan regions.
467 serialOverscanColMean = afwMath.makeStatistics(
468 serialOverscanImage.array[:, idx],
469 self.statType,
470 self.statControl,
471 ).getValue()
472 columns.append(idx)
473 # The overscan input is always in adu, but it only
474 # makes sense to measure CTI in electron units.
475 columnValues.append(serialOverscanColMean)
477 # We want these relative to the readout corner. If that's
478 # on the right side, we need to swap them.
479 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
480 ampStats["SERIAL_OVERSCAN_COLUMNS"] = list(reversed(columns))
481 ampStats["SERIAL_OVERSCAN_VALUES"] = list(reversed(columnValues))
482 else:
483 ampStats["SERIAL_OVERSCAN_COLUMNS"] = columns
484 ampStats["SERIAL_OVERSCAN_VALUES"] = columnValues
486 # Calculate parallel overscan statistics
487 rows = []
488 rowValues = []
489 for idx in range(0, nParallelOverscanRows):
490 parallelOverscanRowMean = afwMath.makeStatistics(
491 parallelOverscanImage.array[idx, :],
492 self.statType,
493 self.statControl,
494 ).getValue()
495 rows.append(idx)
496 # The overscan input is always in adu, but it only
497 # makes sense to measure CTI in electron units.
498 rowValues.append(parallelOverscanRowMean)
500 # We want these relative to the readout corner. If that's
501 # on the top, we need to swap them.
502 if readoutCorner in (ReadoutCorner.UR, ReadoutCorner.UL):
503 ampStats["PARALLEL_OVERSCAN_ROWS"] = list(reversed(rows))
504 ampStats["PARALLEL_OVERSCAN_VALUES"] = list(reversed(rowValues))
505 else:
506 ampStats["PARALLEL_OVERSCAN_ROWS"] = rows
507 ampStats["PARALLEL_OVERSCAN_VALUES"] = rowValues
509 outputStats[amp.getName()] = ampStats
511 return outputStats
513 def measureCtiLegacy(self, inputExp, serialOverscans, gains):
514 """Task to measure CTI statistics.
516 Parameters
517 ----------
518 inputExp : `lsst.afw.image.Exposure`
519 Exposure to measure.
520 serialOverscans : `list` [`lsst.pipe.base.Struct`]
521 List of serial overscan results (expects base units of adu).
522 Expected fields are:
524 ``imageFit``
525 Value or fit subtracted from the amplifier image data
526 (scalar or `lsst.afw.image.Image`).
527 ``overscanFit``
528 Value or fit subtracted from the overscan image data
529 (scalar or `lsst.afw.image.Image`).
530 ``overscanImage``
531 Image of the overscan region with the overscan
532 correction applied (`lsst.afw.image.Image`). This
533 quantity is used to estimate the amplifier read noise
534 empirically.
535 gains : `dict` [`str` `float`]
536 Dictionary of per-amplifier gains, indexed by amplifier name.
538 Returns
539 -------
540 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
541 Dictionary of measurements, keyed by amplifier name and
542 statistics segment. Everything in units based on electron.
543 """
544 self.log.warning("Using the legacy version of CTI statistics may give wrong CTI")
546 outputStats = {}
548 detector = inputExp.getDetector()
549 image = inputExp.image
551 # It only makes sense to measure CTI in electron units.
552 # Make it so.
553 imageUnits = inputExp.getMetadata().get("LSST ISR UNITS")
554 applyGain = False
555 if imageUnits == "adu":
556 applyGain = True
558 # Check if the image is trimmed.
559 isTrimmed = isTrimmedExposure(inputExp)
561 # Ensure we have the same number of overscans as amplifiers.
562 assert len(serialOverscans) == len(detector.getAmplifiers())
564 with gainContext(inputExp, image, applyGain, gains, isTrimmed=isTrimmed):
565 for ampIter, amp in enumerate(detector.getAmplifiers()):
566 ampStats = {}
567 readoutCorner = amp.getReadoutCorner()
569 ampStats["INPUT_GAIN"] = float(gains[amp.getName()])
570 turnoff = inputExp.metadata[f"LSST ISR PTCTURNOFF {amp.getName()}"]
571 ampStats["INPUT_PTCTURNOFF"] = float(turnoff)
573 # Full data region.
574 dataRegion = image[amp.getBBox() if isTrimmed else amp.getRawDataBBox()]
576 # Get the mean of the image
577 ampStats["IMAGE_MEAN"] = afwMath.makeStatistics(dataRegion, self.statType,
578 self.statControl).getValue()
580 # First and last image columns.
581 pixelA = afwMath.makeStatistics(dataRegion.array[:, 0],
582 self.statType,
583 self.statControl).getValue()
584 pixelZ = afwMath.makeStatistics(dataRegion.array[:, -1],
585 self.statType,
586 self.statControl).getValue()
588 # We want these relative to the readout corner. If that's
589 # on the right side, we need to swap them.
590 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
591 ampStats["FIRST_MEAN"] = pixelZ
592 ampStats["LAST_MEAN"] = pixelA
593 else:
594 ampStats["FIRST_MEAN"] = pixelA
595 ampStats["LAST_MEAN"] = pixelZ
597 # Measure the columns of the overscan.
598 if serialOverscans[ampIter] is None:
599 # The amplifier is likely entirely bad, and needs to
600 # be skipped.
601 self.log.warning("No overscan information available for ISR statistics for amp %s.",
602 amp.getName())
603 nCols = amp.getRawSerialOverscanBBox().getWidth()
604 ampStats["OVERSCAN_COLUMNS"] = np.full((nCols, ), np.nan)
605 ampStats["OVERSCAN_VALUES"] = np.full((nCols, ), np.nan)
606 else:
607 overscanImage = serialOverscans[ampIter].overscanImage
609 columns = []
610 values = []
611 for column in range(0, overscanImage.getWidth()):
612 # If overscan.doParallelOverscan=True, the
613 # overscanImage will contain both the serial
614 # and parallel overscan regions.
615 # Only the serial CTI correction has been
616 # implemented, so we must select only the
617 # serial overscan rows for a given column.
618 nRows = amp.getRawSerialOverscanBBox().getHeight()
619 osMean = afwMath.makeStatistics(overscanImage.image.array[:nRows, column],
620 self.statType, self.statControl).getValue()
621 columns.append(column)
622 # The overscan input is always in adu, but it only
623 # makes sense to measure CTI in electron units.
624 values.append(osMean * gains[amp.getName()])
626 # We want these relative to the readout corner. If that's
627 # on the right side, we need to swap them.
628 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
629 ampStats["OVERSCAN_COLUMNS"] = list(reversed(columns))
630 ampStats["OVERSCAN_VALUES"] = list(reversed(values))
631 else:
632 ampStats["OVERSCAN_COLUMNS"] = columns
633 ampStats["OVERSCAN_VALUES"] = values
635 outputStats[amp.getName()] = ampStats
637 return outputStats
639 @staticmethod
640 def makeKernel(kernelSize):
641 """Make a boxcar smoothing kernel.
643 Parameters
644 ----------
645 kernelSize : `int`
646 Size of the kernel in pixels.
648 Returns
649 -------
650 kernel : `np.array`
651 Kernel for boxcar smoothing.
652 """
653 if kernelSize > 0:
654 kernel = np.full(kernelSize, 1.0 / kernelSize)
655 else:
656 kernel = np.array([1.0])
657 return kernel
659 def measureBanding(self, inputExp, overscans):
660 """Task to measure banding statistics.
662 Parameters
663 ----------
664 inputExp : `lsst.afw.image.Exposure`
665 Exposure to measure.
666 overscans : `list` [`lsst.pipe.base.Struct`]
667 List of overscan results. Expected fields are:
669 ``imageFit``
670 Value or fit subtracted from the amplifier image data
671 (scalar or `lsst.afw.image.Image`).
672 ``overscanFit``
673 Value or fit subtracted from the overscan image data
674 (scalar or `lsst.afw.image.Image`).
675 ``overscanImage``
676 Image of the overscan region with the overscan
677 correction applied (`lsst.afw.image.Image`). This
678 quantity is used to estimate the amplifier read noise
679 empirically.
681 Returns
682 -------
683 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
684 Dictionary of measurements, keyed by amplifier name and
685 statistics segment.
686 """
687 outputStats = {}
689 detector = inputExp.getDetector()
690 kernel = self.makeKernel(self.config.bandingKernelSize)
692 outputStats["AMP_BANDING"] = []
693 for amp, overscanData in zip(detector.getAmplifiers(), overscans):
694 overscanFit = np.array(overscanData.overscanFit)
695 overscanArray = overscanData.overscanImage.image.array
696 rawOverscan = np.mean(overscanArray + overscanFit, axis=1)
698 smoothedOverscan = np.convolve(rawOverscan, kernel, mode="valid")
700 low, high = np.quantile(smoothedOverscan, [self.config.bandingFractionLow,
701 self.config.bandingFractionHigh])
702 outputStats["AMP_BANDING"].append(float(high - low))
704 if self.config.bandingUseHalfDetector:
705 fullLength = len(outputStats["AMP_BANDING"])
706 outputStats["DET_BANDING"] = float(np.nanmedian(outputStats["AMP_BANDING"][0:fullLength//2]))
707 else:
708 outputStats["DET_BANDING"] = float(np.nanmedian(outputStats["AMP_BANDING"]))
710 return outputStats
712 def measureProjectionStatistics(self, inputExp, overscans):
713 """Task to measure metrics from image slicing.
715 Parameters
716 ----------
717 inputExp : `lsst.afw.image.Exposure`
718 Exposure to measure.
719 overscans : `list` [`lsst.pipe.base.Struct`]
720 List of overscan results. Expected fields are:
722 ``imageFit``
723 Value or fit subtracted from the amplifier image data
724 (scalar or `lsst.afw.image.Image`).
725 ``overscanFit``
726 Value or fit subtracted from the overscan image data
727 (scalar or `lsst.afw.image.Image`).
728 ``overscanImage``
729 Image of the overscan region with the overscan
730 correction applied (`lsst.afw.image.Image`). This
731 quantity is used to estimate the amplifier read noise
732 empirically.
734 Returns
735 -------
736 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
737 Dictionary of measurements, keyed by amplifier name and
738 statistics segment.
739 """
740 outputStats = {}
742 detector = inputExp.getDetector()
743 kernel = self.makeKernel(self.config.projectionKernelSize)
745 outputStats["AMP_VPROJECTION"] = {}
746 outputStats["AMP_HPROJECTION"] = {}
747 convolveMode = "valid"
748 if self.config.doProjectionFft:
749 outputStats["AMP_VFFT_REAL"] = {}
750 outputStats["AMP_VFFT_IMAG"] = {}
751 outputStats["AMP_HFFT_REAL"] = {}
752 outputStats["AMP_HFFT_IMAG"] = {}
753 convolveMode = "same"
755 for amp in detector.getAmplifiers():
756 ampArray = inputExp.image[amp.getBBox()].array
758 horizontalProjection = np.mean(ampArray, axis=0)
759 verticalProjection = np.mean(ampArray, axis=1)
761 horizontalProjection = np.convolve(horizontalProjection, kernel, mode=convolveMode)
762 verticalProjection = np.convolve(verticalProjection, kernel, mode=convolveMode)
764 outputStats["AMP_HPROJECTION"][amp.getName()] = horizontalProjection.tolist()
765 outputStats["AMP_VPROJECTION"][amp.getName()] = verticalProjection.tolist()
767 if self.config.doProjectionFft:
768 horizontalWindow = np.ones_like(horizontalProjection)
769 verticalWindow = np.ones_like(verticalProjection)
770 if self.config.projectionFftWindow == "NONE":
771 pass
772 elif self.config.projectionFftWindow == "HAMMING":
773 horizontalWindow = hamming(len(horizontalProjection))
774 verticalWindow = hamming(len(verticalProjection))
775 elif self.config.projectionFftWindow == "HANN":
776 horizontalWindow = hann(len(horizontalProjection))
777 verticalWindow = hann(len(verticalProjection))
778 elif self.config.projectionFftWindow == "GAUSSIAN":
779 horizontalWindow = gaussian(len(horizontalProjection))
780 verticalWindow = gaussian(len(verticalProjection))
781 else:
782 raise RuntimeError(f"Invalid window function: {self.config.projectionFftWindow}")
784 horizontalFFT = np.fft.rfft(np.multiply(horizontalProjection, horizontalWindow))
785 verticalFFT = np.fft.rfft(np.multiply(verticalProjection, verticalWindow))
787 outputStats["AMP_HFFT_REAL"][amp.getName()] = np.real(horizontalFFT).tolist()
788 outputStats["AMP_HFFT_IMAG"][amp.getName()] = np.imag(horizontalFFT).tolist()
789 outputStats["AMP_VFFT_REAL"][amp.getName()] = np.real(verticalFFT).tolist()
790 outputStats["AMP_VFFT_IMAG"][amp.getName()] = np.imag(verticalFFT).tolist()
792 return outputStats
794 def copyCalibDistributionStatistics(self, inputExp, **kwargs):
795 """Copy calibration statistics for this exposure.
797 Parameters
798 ----------
799 inputExp : `lsst.afw.image.Exposure`
800 The exposure being processed.
801 **kwargs :
802 Keyword arguments with calibrations.
804 Returns
805 -------
806 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
807 Dictionary of measurements, keyed by amplifier name and
808 statistics segment.
809 """
810 outputStats = {}
812 # Amp level elements
813 for amp in inputExp.getDetector():
814 ampStats = {}
816 for calibType in ("bias", "dark", "flat"):
817 if kwargs.get(calibType, None) is not None:
818 metadata = kwargs[calibType].getMetadata()
819 for pct in self.config.expectedDistributionLevels:
820 key = f"LSST CALIB {calibType.upper()} {amp.getName()} DISTRIBUTION {pct}-PCT"
821 ampStats[key] = metadata.get(key, np.nan)
823 for calibType in ("defects"):
824 if kwargs.get(calibType, None) is not None:
825 metadata = kwargs[calibType].getMetadata()
826 for key in (f"LSST CALIB {calibType.upper()} {amp.getName()} N_HOT",
827 f"LSST CALIB {calibType.upper()} {amp.getName()} N_COLD"):
828 ampStats[key] = metadata.get(key, np.nan)
829 outputStats[amp.getName()] = ampStats
831 # Detector level elements
832 for calibType in ("defects"):
833 if kwargs.get(calibType, None) is not None:
834 metadata = kwargs[calibType].getMetadata()
835 for key in (f"LSST CALIB {calibType.upper()} N_BAD_COLUMNS"):
836 outputStats["detector"][key] = metadata.get(key, np.nan)
838 return outputStats
840 def measureBiasShifts(self, inputExp, serialOverscanResults):
841 """Measure number of bias shifts from overscan data.
843 Parameters
844 ----------
845 inputExp : `lsst.afw.image.Exposure`
846 Exposure to measure.
847 overscans : `list` [`lsst.pipe.base.Struct`]
848 List of overscan results. Expected fields are:
850 ``imageFit``
851 Value or fit subtracted from the amplifier image data
852 (scalar or `lsst.afw.image.Image`).
853 ``overscanFit``
854 Value or fit subtracted from the overscan image data
855 (scalar or `lsst.afw.image.Image`).
856 ``overscanImage``
857 Image of the overscan region with the overscan
858 correction applied (`lsst.afw.image.Image`). This
859 quantity is used to estimate the amplifier read noise
860 empirically.
862 Returns
863 -------
864 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
865 Dictionary of measurements, keyed by amplifier name and
866 statistics segment.
868 Notes
869 -----
870 Based on eo_pipe implementation:
871 https://github.com/lsst-camera-dh/eo_pipe/blob/main/python/lsst/eo/pipe/biasShiftsTask.py # noqa: E501 W505
872 """
873 outputStats = {}
875 detector = inputExp.getDetector()
876 for amp, overscans in zip(detector, serialOverscanResults):
877 ampStats = {}
879 # Check if the overscan results are None
880 # This can happen if the amp is in the
881 # input PTC's badAmp list, for instance.
882 if overscans is None:
883 ampStats["LOCAL_NOISE"] = 0.0
884 ampStats["BIAS_SHIFTS"] = []
885 outputStats[amp.getName()] = ampStats
886 continue
888 # Add fit back to data
889 rawOverscan = overscans.overscanImage.image.array + overscans.overscanFit
891 # Collapse array, skipping first three columns
892 rawOverscan = np.mean(rawOverscan[:, self.config.biasShiftColumnSkip:], axis=1)
894 # Scan for shifts
895 noise, shift_peaks = self._scan_for_shifts(rawOverscan)
896 ampStats["LOCAL_NOISE"] = float(noise)
897 ampStats["BIAS_SHIFTS"] = shift_peaks
899 outputStats[amp.getName()] = ampStats
900 return outputStats
902 def _scan_for_shifts(self, overscanData):
903 """Scan overscan data for shifts.
905 Parameters
906 ----------
907 overscanData : `list` [`float`]
908 Overscan data to search for shifts.
910 Returns
911 -------
912 noise : `float`
913 Noise estimated from Butterworth filtered overscan data.
914 peaks : `list` [`float`, `float`, `int`, `int`]
915 Shift peak information, containing the convolved peak
916 value, the raw peak value, and the lower and upper bounds
917 of the region checked.
918 """
919 numerator, denominator = butter(self.config.biasShiftFilterOrder,
920 self.config.biasShiftCutoff,
921 btype="high", analog=False)
922 noise = np.std(filtfilt(numerator, denominator, overscanData))
923 kernel = np.concatenate([np.arange(self.config.biasShiftWindow),
924 np.arange(-self.config.biasShiftWindow + 1, 0)])
925 kernel = kernel/np.sum(kernel[:self.config.biasShiftWindow])
927 convolved = np.convolve(overscanData, kernel, mode="valid")
928 convolved = np.pad(convolved, (self.config.biasShiftWindow - 1, self.config.biasShiftWindow))
930 shift_check = np.abs(convolved)/noise
931 shift_mask = shift_check > self.config.biasShiftThreshold
932 shift_mask[:self.config.biasShiftRowSkip] = False
934 shift_regions = np.flatnonzero(np.diff(np.r_[np.int8(0),
935 shift_mask.view(np.int8),
936 np.int8(0)])).reshape(-1, 2)
937 shift_peaks = []
938 for region in shift_regions:
939 region_peak = np.argmax(shift_check[region[0]:region[1]]) + region[0]
940 if self._satisfies_flatness(region_peak, convolved[region_peak], overscanData):
941 shift_peaks.append(
942 [float(convolved[region_peak]), float(region_peak),
943 int(region[0]), int(region[1])])
944 return noise, shift_peaks
946 def _satisfies_flatness(self, shiftRow, shiftPeak, overscanData):
947 """Determine if a region is flat.
949 Parameters
950 ----------
951 shiftRow : `int`
952 Row with possible peak.
953 shiftPeak : `float`
954 Value at the possible peak.
955 overscanData : `list` [`float`]
956 Overscan data used to fit around the possible peak.
958 Returns
959 -------
960 isFlat : `bool`
961 Indicates if the region is flat, and so the peak is valid.
962 """
963 prerange = np.arange(shiftRow - self.config.biasShiftWindow, shiftRow)
964 postrange = np.arange(shiftRow, shiftRow + self.config.biasShiftWindow)
966 preFit = linregress(prerange, overscanData[prerange])
967 postFit = linregress(postrange, overscanData[postrange])
969 if shiftPeak > 0:
970 preTrend = (2*preFit[0]*len(prerange) < shiftPeak)
971 postTrend = (2*postFit[0]*len(postrange) < shiftPeak)
972 else:
973 preTrend = (2*preFit[0]*len(prerange) > shiftPeak)
974 postTrend = (2*postFit[0]*len(postrange) > shiftPeak)
976 return (preTrend and postTrend)
978 def measureAmpCorrelations(self, inputExp, serialOverscanResults):
979 """Measure correlations between amplifier segments.
981 Parameters
982 ----------
983 inputExp : `lsst.afw.image.Exposure`
984 Exposure to measure.
985 overscans : `list` [`lsst.pipe.base.Struct`]
986 List of overscan results. Expected fields are:
988 ``imageFit``
989 Value or fit subtracted from the amplifier image data
990 (scalar or `lsst.afw.image.Image`).
991 ``overscanFit``
992 Value or fit subtracted from the overscan image data
993 (scalar or `lsst.afw.image.Image`).
994 ``overscanImage``
995 Image of the overscan region with the overscan
996 correction applied (`lsst.afw.image.Image`). This
997 quantity is used to estimate the amplifier read noise
998 empirically.
1000 Returns
1001 -------
1002 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
1003 Dictionary of measurements, keyed by amplifier name and
1004 statistics segment.
1006 Notes
1007 -----
1008 Based on eo_pipe implementation:
1009 https://github.com/lsst-camera-dh/eo_pipe/blob/main/python/lsst/eo/pipe/raft_level_correlations.py # noqa: E501 W505
1010 """
1011 detector = inputExp.getDetector()
1012 rows = []
1014 for ampId, overscan in enumerate(serialOverscanResults):
1015 # Check if the overscan results are None
1016 # This can happen if the amp is in the
1017 # input PTC's badAmp list, for instance.
1018 if overscan is None:
1019 for ampId2, overscan2 in enumerate(serialOverscanResults):
1020 if ampId2 != ampId:
1021 osC = 0.0
1022 imC = 0.0
1023 row = {
1024 'detector': detector.getId(),
1025 'detectorComp': detector.getId(),
1026 'ampName': detector[ampId].getName(),
1027 'ampComp': detector[ampId2].getName(),
1028 'imageCorr': float(imC),
1029 'overscanCorr': float(osC),
1030 }
1031 rows.append(row)
1032 else:
1033 rawOverscan = overscan.overscanImage.image.array + overscan.overscanFit
1034 rawOverscan = rawOverscan.ravel()
1036 ampImage = inputExp[detector[ampId].getBBox()]
1037 ampImage = ampImage.image.array.ravel()
1039 for ampId2, overscan2 in enumerate(serialOverscanResults):
1040 osC = 1.0
1041 imC = 1.0
1042 if ampId2 != ampId:
1043 # Check if the overscan results are None
1044 # This can happen if the amp is in the
1045 # input PTC's badAmp list, for instance.
1046 if overscan2 is None:
1047 osC = 0.0
1048 imC = 0.0
1049 else:
1050 rawOverscan2 = overscan2.overscanImage.image.array + overscan2.overscanFit
1051 rawOverscan2 = rawOverscan2.ravel()
1053 osC = np.corrcoef(rawOverscan, rawOverscan2)[0, 1]
1055 ampImage2 = inputExp[detector[ampId2].getBBox()]
1056 ampImage2 = ampImage2.image.array.ravel()
1058 imC = np.corrcoef(ampImage, ampImage2)[0, 1]
1059 row = {
1060 'detector': detector.getId(),
1061 'detectorComp': detector.getId(),
1062 'ampName': detector[ampId].getName(),
1063 'ampComp': detector[ampId2].getName(),
1064 'imageCorr': float(imC),
1065 'overscanCorr': float(osC),
1066 }
1067 rows.append(row)
1068 return rows
1070 def measureDivisaderoStatistics(self, inputExp, **kwargs):
1071 """Measure Max Divisadero Tearing effect per amp.
1073 Parameters
1074 ----------
1075 inputExp : `lsst.afw.image.Exposure`
1076 Exposure to measure. Usually a flat.
1077 **kwargs :
1078 The flat will be selected from here.
1080 Returns
1081 -------
1082 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
1083 Dictionary of measurements, keyed by amplifier name and
1084 statistics segment.
1085 Measurements include
1087 - DIVISADERO_PROFILE: Robust mean of rows between
1088 divisaderoProjection<Maximum|Minumum> on readout edge of ccd
1089 normalized by a linear fit to the same rows.
1090 - DIVISADERO_MAX_PAIR: Tuple of maximum of the absolute values of
1091 the DIVISADERO_PROFILE, for number of pixels (specified by
1092 divisaderoNumImpactPixels on left and right side of amp.
1093 - DIVISADERO_MAX: Maximum of the absolute values of the
1094 the DIVISADERO_PROFILE, for the divisaderoNumImpactPixels on
1095 boundaries of neighboring amps (including the pixels in those
1096 neighborboring amps).
1097 """
1098 outputStats = {}
1100 for amp in inputExp.getDetector():
1101 # Copy unneeded if we do not ever modify the array by flipping
1102 ampArray = inputExp.image[amp.getBBox()].array
1103 # slice the outer top or bottom of the amp: the readout side
1104 if amp.getReadoutCorner().name in ('UL', 'UR'):
1105 minRow = amp.getBBox().getHeight() - self.config.divisaderoProjectionMaximum
1106 maxRow = amp.getBBox().getHeight() - self.config.divisaderoProjectionMinimum
1107 else:
1108 minRow = self.config.divisaderoProjectionMinimum
1109 maxRow = self.config.divisaderoProjectionMaximum
1111 segment = slice(minRow, maxRow)
1112 projection, _, _ = astropy.stats.sigma_clipped_stats(ampArray[segment, :], axis=0)
1114 ampStats = {}
1115 projection /= np.median(projection)
1116 columns = np.arange(len(projection))
1118 segment = slice(self.config.divisaderoEdgePixels, -self.config.divisaderoEdgePixels)
1119 model = np.polyfit(columns[segment], projection[segment], 1)
1120 modelProjection = model[0] * columns + model[1]
1121 divisaderoProfile = projection / modelProjection
1123 # look for max at the edges:
1124 leftMax = np.nanmax(np.abs(divisaderoProfile[0:self.config.divisaderoNumImpactPixels] - 1.0))
1125 rightMax = np.nanmax(np.abs(divisaderoProfile[-self.config.divisaderoNumImpactPixels:] - 1.0))
1127 ampStats['DIVISADERO_PROFILE'] = np.array(divisaderoProfile).tolist()
1128 ampStats['DIVISADERO_MAX_PAIR'] = [leftMax, rightMax]
1129 outputStats[amp.getName()] = ampStats
1131 detector = inputExp.getDetector()
1132 xCenters = [amp.getBBox().getCenterX() for amp in detector]
1133 yCenters = [amp.getBBox().getCenterY() for amp in detector]
1134 xIndices = np.ceil(xCenters / np.min(xCenters) / 2).astype(int) - 1
1135 yIndices = np.ceil(yCenters / np.min(yCenters) / 2).astype(int) - 1
1136 ampIds = np.zeros((len(set(yIndices)), len(set(xIndices))), dtype=int)
1137 for ampId, xIndex, yIndex in zip(np.arange(len(detector)), xIndices, yIndices):
1138 ampIds[yIndex, xIndex] = ampId
1140 # Loop over amps again because the DIVISIDERO_MAX will be the max
1141 # of the profile on its boundary with its neighboring amps
1142 for i, amp in enumerate(detector):
1143 y, x = np.where(ampIds == i)
1144 end = ampIds.shape[1] - 1
1145 xInd = x[0]
1146 yInd = y[0]
1147 thisAmpsPair = outputStats[amp.getName()]['DIVISADERO_MAX_PAIR']
1149 if x == 0:
1150 # leftmost amp: take the max of your right side and
1151 myMax = thisAmpsPair[1]
1152 # your neighbor's left side
1153 neighborMax = outputStats[detector[ampIds[yInd, 1]].getName()]['DIVISADERO_MAX_PAIR'][0]
1154 elif x == end:
1155 # rightmost amp: take the max of your left side and
1156 myMax = thisAmpsPair[0]
1157 # your neighbor's right side
1158 neighborMax = outputStats[detector[ampIds[yInd, end - 1]].getName()]['DIVISADERO_MAX_PAIR'][1]
1159 else:
1160 # Middle amp: take the max of both your own sides and the
1161 myMax = max(thisAmpsPair)
1162 leftName = detector[ampIds[yInd, max(xInd - 1, 0)]].getName()
1163 rightName = detector[ampIds[yInd, min(xInd + 1, ampIds.shape[1] - 1)]].getName()
1164 # right side of the neighbor to your left
1165 # and left side of your neighbor to your right
1166 neighborMax = max(outputStats[leftName]['DIVISADERO_MAX_PAIR'][1],
1167 outputStats[rightName]['DIVISADERO_MAX_PAIR'][0])
1169 divisaderoMax = max([myMax, neighborMax])
1170 outputStats[amp.getName()]['DIVISADERO_MAX'] = divisaderoMax
1172 return outputStats