Coverage for python/lsst/ip/isr/isrStatistics.py: 11%
402 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-10 01:49 -0700
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 # Full data region.
401 dataRegion = image[amp.getBBox()]
402 serialOverscanImage = untrimmedImage[amp.getRawSerialOverscanBBox()]
403 parallelOverscanImage = untrimmedImage[amp.getRawSerialOverscanBBox()]
405 # Get the mean of the image
406 ampStats["IMAGE_MEAN"] = afwMath.makeStatistics(dataRegion, self.statType,
407 self.statControl).getValue()
409 # First and last image columns.
410 colA = afwMath.makeStatistics(
411 dataRegion.array[:, 0],
412 self.statType,
413 self.statControl,
414 ).getValue()
415 colZ = afwMath.makeStatistics(
416 dataRegion.array[:, -1],
417 self.statType,
418 self.statControl,
419 ).getValue()
421 # First and last image rows.
422 rowA = afwMath.makeStatistics(
423 dataRegion.array[0, :],
424 self.statType,
425 self.statControl,
426 ).getValue()
427 rowZ = afwMath.makeStatistics(
428 dataRegion.array[-1, :],
429 self.statType,
430 self.statControl,
431 ).getValue()
433 # We want these relative to the readout corner. If that's
434 # on the right side, we need to swap them.
435 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
436 ampStats["FIRST_COLUMN_MEAN"] = colZ
437 ampStats["LAST_COLUMN_MEAN"] = colA
438 else:
439 ampStats["FIRST_COLUMN_MEAN"] = colA
440 ampStats["LAST_COLUMN_MEAN"] = colZ
442 # We want these relative to the readout corner. If that's
443 # on the top, we need to swap them.
444 if readoutCorner in (ReadoutCorner.UR, ReadoutCorner.UL):
445 ampStats["FIRST_ROW_MEAN"] = rowZ
446 ampStats["LAST_ROW_MEAN"] = rowA
447 else:
448 ampStats["FIRST_ROW_MEAN"] = rowA
449 ampStats["LAST_ROW_MEAN"] = rowZ
451 # Measure the columns of the serial overscan and the rows
452 # of the parallel overscan.
453 nSerialOverscanCols = amp.getRawSerialOverscanBBox().getWidth()
454 nParallelOverscanRows = amp.getRawParallelOverscanBBox().getHeight()
456 # Calculate serial overscan statistics
457 columns = []
458 columnValues = []
459 for idx in range(0, nSerialOverscanCols):
460 # If overscan.doParallelOverscan=True, the
461 # overscanImage will contain both the serial
462 # and parallel overscan regions.
463 serialOverscanColMean = afwMath.makeStatistics(
464 serialOverscanImage.array[:, idx],
465 self.statType,
466 self.statControl,
467 ).getValue()
468 columns.append(idx)
469 # The overscan input is always in adu, but it only
470 # makes sense to measure CTI in electron units.
471 columnValues.append(serialOverscanColMean)
473 # We want these relative to the readout corner. If that's
474 # on the right side, we need to swap them.
475 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
476 ampStats["SERIAL_OVERSCAN_COLUMNS"] = list(reversed(columns))
477 ampStats["SERIAL_OVERSCAN_VALUES"] = list(reversed(columnValues))
478 else:
479 ampStats["SERIAL_OVERSCAN_COLUMNS"] = columns
480 ampStats["SERIAL_OVERSCAN_VALUES"] = columnValues
482 # Calculate parallel overscan statistics
483 rows = []
484 rowValues = []
485 for idx in range(0, nParallelOverscanRows):
486 parallelOverscanRowMean = afwMath.makeStatistics(
487 parallelOverscanImage.array[idx, :],
488 self.statType,
489 self.statControl,
490 ).getValue()
491 rows.append(idx)
492 # The overscan input is always in adu, but it only
493 # makes sense to measure CTI in electron units.
494 rowValues.append(parallelOverscanRowMean)
496 # We want these relative to the readout corner. If that's
497 # on the top, we need to swap them.
498 if readoutCorner in (ReadoutCorner.UR, ReadoutCorner.UL):
499 ampStats["PARALLEL_OVERSCAN_ROWS"] = list(reversed(rows))
500 ampStats["PARALLEL_OVERSCAN_VALUES"] = list(reversed(rowValues))
501 else:
502 ampStats["PARALLEL_OVERSCAN_ROWS"] = rows
503 ampStats["PARALLEL_OVERSCAN_VALUES"] = rowValues
505 outputStats[amp.getName()] = ampStats
507 return outputStats
509 def measureCtiLegacy(self, inputExp, serialOverscans, gains):
510 """Task to measure CTI statistics.
512 Parameters
513 ----------
514 inputExp : `lsst.afw.image.Exposure`
515 Exposure to measure.
516 serialOverscans : `list` [`lsst.pipe.base.Struct`]
517 List of serial overscan results (expects base units of adu).
518 Expected fields are:
520 ``imageFit``
521 Value or fit subtracted from the amplifier image data
522 (scalar or `lsst.afw.image.Image`).
523 ``overscanFit``
524 Value or fit subtracted from the overscan image data
525 (scalar or `lsst.afw.image.Image`).
526 ``overscanImage``
527 Image of the overscan region with the overscan
528 correction applied (`lsst.afw.image.Image`). This
529 quantity is used to estimate the amplifier read noise
530 empirically.
531 gains : `dict` [`str` `float`]
532 Dictionary of per-amplifier gains, indexed by amplifier name.
534 Returns
535 -------
536 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
537 Dictionary of measurements, keyed by amplifier name and
538 statistics segment. Everything in units based on electron.
539 """
540 self.log.warning("Using the legacy version of CTI statistics may give wrong CTI")
542 outputStats = {}
544 detector = inputExp.getDetector()
545 image = inputExp.image
547 # It only makes sense to measure CTI in electron units.
548 # Make it so.
549 imageUnits = inputExp.getMetadata().get("LSST ISR UNITS")
550 applyGain = False
551 if imageUnits == "adu":
552 applyGain = True
554 # Check if the image is trimmed.
555 isTrimmed = isTrimmedExposure(inputExp)
557 # Ensure we have the same number of overscans as amplifiers.
558 assert len(serialOverscans) == len(detector.getAmplifiers())
560 with gainContext(inputExp, image, applyGain, gains, isTrimmed=isTrimmed):
561 for ampIter, amp in enumerate(detector.getAmplifiers()):
562 ampStats = {}
563 readoutCorner = amp.getReadoutCorner()
565 # Full data region.
566 dataRegion = image[amp.getBBox() if isTrimmed else amp.getRawDataBBox()]
568 # Get the mean of the image
569 ampStats["IMAGE_MEAN"] = afwMath.makeStatistics(dataRegion, self.statType,
570 self.statControl).getValue()
572 # First and last image columns.
573 pixelA = afwMath.makeStatistics(dataRegion.array[:, 0],
574 self.statType,
575 self.statControl).getValue()
576 pixelZ = afwMath.makeStatistics(dataRegion.array[:, -1],
577 self.statType,
578 self.statControl).getValue()
580 # We want these relative to the readout corner. If that's
581 # on the right side, we need to swap them.
582 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
583 ampStats["FIRST_MEAN"] = pixelZ
584 ampStats["LAST_MEAN"] = pixelA
585 else:
586 ampStats["FIRST_MEAN"] = pixelA
587 ampStats["LAST_MEAN"] = pixelZ
589 # Measure the columns of the overscan.
590 if serialOverscans[ampIter] is None:
591 # The amplifier is likely entirely bad, and needs to
592 # be skipped.
593 self.log.warning("No overscan information available for ISR statistics for amp %s.",
594 amp.getName())
595 nCols = amp.getRawSerialOverscanBBox().getWidth()
596 ampStats["OVERSCAN_COLUMNS"] = np.full((nCols, ), np.nan)
597 ampStats["OVERSCAN_VALUES"] = np.full((nCols, ), np.nan)
598 else:
599 overscanImage = serialOverscans[ampIter].overscanImage
601 columns = []
602 values = []
603 for column in range(0, overscanImage.getWidth()):
604 # If overscan.doParallelOverscan=True, the
605 # overscanImage will contain both the serial
606 # and parallel overscan regions.
607 # Only the serial CTI correction has been
608 # implemented, so we must select only the
609 # serial overscan rows for a given column.
610 nRows = amp.getRawSerialOverscanBBox().getHeight()
611 osMean = afwMath.makeStatistics(overscanImage.image.array[:nRows, column],
612 self.statType, self.statControl).getValue()
613 columns.append(column)
614 # The overscan input is always in adu, but it only
615 # makes sense to measure CTI in electron units.
616 values.append(osMean * gains[amp.getName()])
618 # We want these relative to the readout corner. If that's
619 # on the right side, we need to swap them.
620 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
621 ampStats["OVERSCAN_COLUMNS"] = list(reversed(columns))
622 ampStats["OVERSCAN_VALUES"] = list(reversed(values))
623 else:
624 ampStats["OVERSCAN_COLUMNS"] = columns
625 ampStats["OVERSCAN_VALUES"] = values
627 outputStats[amp.getName()] = ampStats
629 return outputStats
631 @staticmethod
632 def makeKernel(kernelSize):
633 """Make a boxcar smoothing kernel.
635 Parameters
636 ----------
637 kernelSize : `int`
638 Size of the kernel in pixels.
640 Returns
641 -------
642 kernel : `np.array`
643 Kernel for boxcar smoothing.
644 """
645 if kernelSize > 0:
646 kernel = np.full(kernelSize, 1.0 / kernelSize)
647 else:
648 kernel = np.array([1.0])
649 return kernel
651 def measureBanding(self, inputExp, overscans):
652 """Task to measure banding statistics.
654 Parameters
655 ----------
656 inputExp : `lsst.afw.image.Exposure`
657 Exposure to measure.
658 overscans : `list` [`lsst.pipe.base.Struct`]
659 List of overscan results. Expected fields are:
661 ``imageFit``
662 Value or fit subtracted from the amplifier image data
663 (scalar or `lsst.afw.image.Image`).
664 ``overscanFit``
665 Value or fit subtracted from the overscan image data
666 (scalar or `lsst.afw.image.Image`).
667 ``overscanImage``
668 Image of the overscan region with the overscan
669 correction applied (`lsst.afw.image.Image`). This
670 quantity is used to estimate the amplifier read noise
671 empirically.
673 Returns
674 -------
675 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
676 Dictionary of measurements, keyed by amplifier name and
677 statistics segment.
678 """
679 outputStats = {}
681 detector = inputExp.getDetector()
682 kernel = self.makeKernel(self.config.bandingKernelSize)
684 outputStats["AMP_BANDING"] = []
685 for amp, overscanData in zip(detector.getAmplifiers(), overscans):
686 overscanFit = np.array(overscanData.overscanFit)
687 overscanArray = overscanData.overscanImage.image.array
688 rawOverscan = np.mean(overscanArray + overscanFit, axis=1)
690 smoothedOverscan = np.convolve(rawOverscan, kernel, mode="valid")
692 low, high = np.quantile(smoothedOverscan, [self.config.bandingFractionLow,
693 self.config.bandingFractionHigh])
694 outputStats["AMP_BANDING"].append(float(high - low))
696 if self.config.bandingUseHalfDetector:
697 fullLength = len(outputStats["AMP_BANDING"])
698 outputStats["DET_BANDING"] = float(np.nanmedian(outputStats["AMP_BANDING"][0:fullLength//2]))
699 else:
700 outputStats["DET_BANDING"] = float(np.nanmedian(outputStats["AMP_BANDING"]))
702 return outputStats
704 def measureProjectionStatistics(self, inputExp, overscans):
705 """Task to measure metrics from image slicing.
707 Parameters
708 ----------
709 inputExp : `lsst.afw.image.Exposure`
710 Exposure to measure.
711 overscans : `list` [`lsst.pipe.base.Struct`]
712 List of overscan results. Expected fields are:
714 ``imageFit``
715 Value or fit subtracted from the amplifier image data
716 (scalar or `lsst.afw.image.Image`).
717 ``overscanFit``
718 Value or fit subtracted from the overscan image data
719 (scalar or `lsst.afw.image.Image`).
720 ``overscanImage``
721 Image of the overscan region with the overscan
722 correction applied (`lsst.afw.image.Image`). This
723 quantity is used to estimate the amplifier read noise
724 empirically.
726 Returns
727 -------
728 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
729 Dictionary of measurements, keyed by amplifier name and
730 statistics segment.
731 """
732 outputStats = {}
734 detector = inputExp.getDetector()
735 kernel = self.makeKernel(self.config.projectionKernelSize)
737 outputStats["AMP_VPROJECTION"] = {}
738 outputStats["AMP_HPROJECTION"] = {}
739 convolveMode = "valid"
740 if self.config.doProjectionFft:
741 outputStats["AMP_VFFT_REAL"] = {}
742 outputStats["AMP_VFFT_IMAG"] = {}
743 outputStats["AMP_HFFT_REAL"] = {}
744 outputStats["AMP_HFFT_IMAG"] = {}
745 convolveMode = "same"
747 for amp in detector.getAmplifiers():
748 ampArray = inputExp.image[amp.getBBox()].array
750 horizontalProjection = np.mean(ampArray, axis=0)
751 verticalProjection = np.mean(ampArray, axis=1)
753 horizontalProjection = np.convolve(horizontalProjection, kernel, mode=convolveMode)
754 verticalProjection = np.convolve(verticalProjection, kernel, mode=convolveMode)
756 outputStats["AMP_HPROJECTION"][amp.getName()] = horizontalProjection.tolist()
757 outputStats["AMP_VPROJECTION"][amp.getName()] = verticalProjection.tolist()
759 if self.config.doProjectionFft:
760 horizontalWindow = np.ones_like(horizontalProjection)
761 verticalWindow = np.ones_like(verticalProjection)
762 if self.config.projectionFftWindow == "NONE":
763 pass
764 elif self.config.projectionFftWindow == "HAMMING":
765 horizontalWindow = hamming(len(horizontalProjection))
766 verticalWindow = hamming(len(verticalProjection))
767 elif self.config.projectionFftWindow == "HANN":
768 horizontalWindow = hann(len(horizontalProjection))
769 verticalWindow = hann(len(verticalProjection))
770 elif self.config.projectionFftWindow == "GAUSSIAN":
771 horizontalWindow = gaussian(len(horizontalProjection))
772 verticalWindow = gaussian(len(verticalProjection))
773 else:
774 raise RuntimeError(f"Invalid window function: {self.config.projectionFftWindow}")
776 horizontalFFT = np.fft.rfft(np.multiply(horizontalProjection, horizontalWindow))
777 verticalFFT = np.fft.rfft(np.multiply(verticalProjection, verticalWindow))
779 outputStats["AMP_HFFT_REAL"][amp.getName()] = np.real(horizontalFFT).tolist()
780 outputStats["AMP_HFFT_IMAG"][amp.getName()] = np.imag(horizontalFFT).tolist()
781 outputStats["AMP_VFFT_REAL"][amp.getName()] = np.real(verticalFFT).tolist()
782 outputStats["AMP_VFFT_IMAG"][amp.getName()] = np.imag(verticalFFT).tolist()
784 return outputStats
786 def copyCalibDistributionStatistics(self, inputExp, **kwargs):
787 """Copy calibration statistics for this exposure.
789 Parameters
790 ----------
791 inputExp : `lsst.afw.image.Exposure`
792 The exposure being processed.
793 **kwargs :
794 Keyword arguments with calibrations.
796 Returns
797 -------
798 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
799 Dictionary of measurements, keyed by amplifier name and
800 statistics segment.
801 """
802 outputStats = {}
804 # Amp level elements
805 for amp in inputExp.getDetector():
806 ampStats = {}
808 for calibType in ("bias", "dark", "flat"):
809 if kwargs.get(calibType, None) is not None:
810 metadata = kwargs[calibType].getMetadata()
811 for pct in self.config.expectedDistributionLevels:
812 key = f"LSST CALIB {calibType.upper()} {amp.getName()} DISTRIBUTION {pct}-PCT"
813 ampStats[key] = metadata.get(key, np.nan)
815 for calibType in ("defects"):
816 if kwargs.get(calibType, None) is not None:
817 metadata = kwargs[calibType].getMetadata()
818 for key in (f"LSST CALIB {calibType.upper()} {amp.getName()} N_HOT",
819 f"LSST CALIB {calibType.upper()} {amp.getName()} N_COLD"):
820 ampStats[key] = metadata.get(key, np.nan)
821 outputStats[amp.getName()] = ampStats
823 # Detector level elements
824 for calibType in ("defects"):
825 if kwargs.get(calibType, None) is not None:
826 metadata = kwargs[calibType].getMetadata()
827 for key in (f"LSST CALIB {calibType.upper()} N_BAD_COLUMNS"):
828 outputStats["detector"][key] = metadata.get(key, np.nan)
830 return outputStats
832 def measureBiasShifts(self, inputExp, serialOverscanResults):
833 """Measure number of bias shifts from overscan data.
835 Parameters
836 ----------
837 inputExp : `lsst.afw.image.Exposure`
838 Exposure to measure.
839 overscans : `list` [`lsst.pipe.base.Struct`]
840 List of overscan results. Expected fields are:
842 ``imageFit``
843 Value or fit subtracted from the amplifier image data
844 (scalar or `lsst.afw.image.Image`).
845 ``overscanFit``
846 Value or fit subtracted from the overscan image data
847 (scalar or `lsst.afw.image.Image`).
848 ``overscanImage``
849 Image of the overscan region with the overscan
850 correction applied (`lsst.afw.image.Image`). This
851 quantity is used to estimate the amplifier read noise
852 empirically.
854 Returns
855 -------
856 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
857 Dictionary of measurements, keyed by amplifier name and
858 statistics segment.
860 Notes
861 -----
862 Based on eo_pipe implementation:
863 https://github.com/lsst-camera-dh/eo_pipe/blob/main/python/lsst/eo/pipe/biasShiftsTask.py # noqa: E501 W505
864 """
865 outputStats = {}
867 detector = inputExp.getDetector()
868 for amp, overscans in zip(detector, serialOverscanResults):
869 ampStats = {}
871 # Check if the overscan results are None
872 # This can happen if the amp is in the
873 # input PTC's badAmp list, for instance.
874 if overscans is None:
875 ampStats["LOCAL_NOISE"] = 0.0
876 ampStats["BIAS_SHIFTS"] = []
877 outputStats[amp.getName()] = ampStats
878 continue
880 # Add fit back to data
881 rawOverscan = overscans.overscanImage.image.array + overscans.overscanFit
883 # Collapse array, skipping first three columns
884 rawOverscan = np.mean(rawOverscan[:, self.config.biasShiftColumnSkip:], axis=1)
886 # Scan for shifts
887 noise, shift_peaks = self._scan_for_shifts(rawOverscan)
888 ampStats["LOCAL_NOISE"] = float(noise)
889 ampStats["BIAS_SHIFTS"] = shift_peaks
891 outputStats[amp.getName()] = ampStats
892 return outputStats
894 def _scan_for_shifts(self, overscanData):
895 """Scan overscan data for shifts.
897 Parameters
898 ----------
899 overscanData : `list` [`float`]
900 Overscan data to search for shifts.
902 Returns
903 -------
904 noise : `float`
905 Noise estimated from Butterworth filtered overscan data.
906 peaks : `list` [`float`, `float`, `int`, `int`]
907 Shift peak information, containing the convolved peak
908 value, the raw peak value, and the lower and upper bounds
909 of the region checked.
910 """
911 numerator, denominator = butter(self.config.biasShiftFilterOrder,
912 self.config.biasShiftCutoff,
913 btype="high", analog=False)
914 noise = np.std(filtfilt(numerator, denominator, overscanData))
915 kernel = np.concatenate([np.arange(self.config.biasShiftWindow),
916 np.arange(-self.config.biasShiftWindow + 1, 0)])
917 kernel = kernel/np.sum(kernel[:self.config.biasShiftWindow])
919 convolved = np.convolve(overscanData, kernel, mode="valid")
920 convolved = np.pad(convolved, (self.config.biasShiftWindow - 1, self.config.biasShiftWindow))
922 shift_check = np.abs(convolved)/noise
923 shift_mask = shift_check > self.config.biasShiftThreshold
924 shift_mask[:self.config.biasShiftRowSkip] = False
926 shift_regions = np.flatnonzero(np.diff(np.r_[np.int8(0),
927 shift_mask.view(np.int8),
928 np.int8(0)])).reshape(-1, 2)
929 shift_peaks = []
930 for region in shift_regions:
931 region_peak = np.argmax(shift_check[region[0]:region[1]]) + region[0]
932 if self._satisfies_flatness(region_peak, convolved[region_peak], overscanData):
933 shift_peaks.append(
934 [float(convolved[region_peak]), float(region_peak),
935 int(region[0]), int(region[1])])
936 return noise, shift_peaks
938 def _satisfies_flatness(self, shiftRow, shiftPeak, overscanData):
939 """Determine if a region is flat.
941 Parameters
942 ----------
943 shiftRow : `int`
944 Row with possible peak.
945 shiftPeak : `float`
946 Value at the possible peak.
947 overscanData : `list` [`float`]
948 Overscan data used to fit around the possible peak.
950 Returns
951 -------
952 isFlat : `bool`
953 Indicates if the region is flat, and so the peak is valid.
954 """
955 prerange = np.arange(shiftRow - self.config.biasShiftWindow, shiftRow)
956 postrange = np.arange(shiftRow, shiftRow + self.config.biasShiftWindow)
958 preFit = linregress(prerange, overscanData[prerange])
959 postFit = linregress(postrange, overscanData[postrange])
961 if shiftPeak > 0:
962 preTrend = (2*preFit[0]*len(prerange) < shiftPeak)
963 postTrend = (2*postFit[0]*len(postrange) < shiftPeak)
964 else:
965 preTrend = (2*preFit[0]*len(prerange) > shiftPeak)
966 postTrend = (2*postFit[0]*len(postrange) > shiftPeak)
968 return (preTrend and postTrend)
970 def measureAmpCorrelations(self, inputExp, serialOverscanResults):
971 """Measure correlations between amplifier segments.
973 Parameters
974 ----------
975 inputExp : `lsst.afw.image.Exposure`
976 Exposure to measure.
977 overscans : `list` [`lsst.pipe.base.Struct`]
978 List of overscan results. Expected fields are:
980 ``imageFit``
981 Value or fit subtracted from the amplifier image data
982 (scalar or `lsst.afw.image.Image`).
983 ``overscanFit``
984 Value or fit subtracted from the overscan image data
985 (scalar or `lsst.afw.image.Image`).
986 ``overscanImage``
987 Image of the overscan region with the overscan
988 correction applied (`lsst.afw.image.Image`). This
989 quantity is used to estimate the amplifier read noise
990 empirically.
992 Returns
993 -------
994 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
995 Dictionary of measurements, keyed by amplifier name and
996 statistics segment.
998 Notes
999 -----
1000 Based on eo_pipe implementation:
1001 https://github.com/lsst-camera-dh/eo_pipe/blob/main/python/lsst/eo/pipe/raft_level_correlations.py # noqa: E501 W505
1002 """
1003 detector = inputExp.getDetector()
1004 rows = []
1006 for ampId, overscan in enumerate(serialOverscanResults):
1007 # Check if the overscan results are None
1008 # This can happen if the amp is in the
1009 # input PTC's badAmp list, for instance.
1010 if overscan is None:
1011 for ampId2, overscan2 in enumerate(serialOverscanResults):
1012 if ampId2 != ampId:
1013 osC = 0.0
1014 imC = 0.0
1015 row = {
1016 'detector': detector.getId(),
1017 'detectorComp': detector.getId(),
1018 'ampName': detector[ampId].getName(),
1019 'ampComp': detector[ampId2].getName(),
1020 'imageCorr': float(imC),
1021 'overscanCorr': float(osC),
1022 }
1023 rows.append(row)
1024 else:
1025 rawOverscan = overscan.overscanImage.image.array + overscan.overscanFit
1026 rawOverscan = rawOverscan.ravel()
1028 ampImage = inputExp[detector[ampId].getBBox()]
1029 ampImage = ampImage.image.array.ravel()
1031 for ampId2, overscan2 in enumerate(serialOverscanResults):
1032 osC = 1.0
1033 imC = 1.0
1034 if ampId2 != ampId:
1035 # Check if the overscan results are None
1036 # This can happen if the amp is in the
1037 # input PTC's badAmp list, for instance.
1038 if overscan2 is None:
1039 osC = 0.0
1040 imC = 0.0
1041 else:
1042 rawOverscan2 = overscan2.overscanImage.image.array + overscan2.overscanFit
1043 rawOverscan2 = rawOverscan2.ravel()
1045 osC = np.corrcoef(rawOverscan, rawOverscan2)[0, 1]
1047 ampImage2 = inputExp[detector[ampId2].getBBox()]
1048 ampImage2 = ampImage2.image.array.ravel()
1050 imC = np.corrcoef(ampImage, ampImage2)[0, 1]
1051 row = {
1052 'detector': detector.getId(),
1053 'detectorComp': detector.getId(),
1054 'ampName': detector[ampId].getName(),
1055 'ampComp': detector[ampId2].getName(),
1056 'imageCorr': float(imC),
1057 'overscanCorr': float(osC),
1058 }
1059 rows.append(row)
1060 return rows
1062 def measureDivisaderoStatistics(self, inputExp, **kwargs):
1063 """Measure Max Divisadero Tearing effect per amp.
1065 Parameters
1066 ----------
1067 inputExp : `lsst.afw.image.Exposure`
1068 Exposure to measure. Usually a flat.
1069 **kwargs :
1070 The flat will be selected from here.
1072 Returns
1073 -------
1074 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
1075 Dictionary of measurements, keyed by amplifier name and
1076 statistics segment.
1077 Measurements include
1079 - DIVISADERO_PROFILE: Robust mean of rows between
1080 divisaderoProjection<Maximum|Minumum> on readout edge of ccd
1081 normalized by a linear fit to the same rows.
1082 - DIVISADERO_MAX_PAIR: Tuple of maximum of the absolute values of
1083 the DIVISADERO_PROFILE, for number of pixels (specified by
1084 divisaderoNumImpactPixels on left and right side of amp.
1085 - DIVISADERO_MAX: Maximum of the absolute values of the
1086 the DIVISADERO_PROFILE, for the divisaderoNumImpactPixels on
1087 boundaries of neighboring amps (including the pixels in those
1088 neighborboring amps).
1089 """
1090 outputStats = {}
1092 for amp in inputExp.getDetector():
1093 # Copy unneeded if we do not ever modify the array by flipping
1094 ampArray = inputExp.image[amp.getBBox()].array
1095 # slice the outer top or bottom of the amp: the readout side
1096 if amp.getReadoutCorner().name in ('UL', 'UR'):
1097 minRow = amp.getBBox().getHeight() - self.config.divisaderoProjectionMaximum
1098 maxRow = amp.getBBox().getHeight() - self.config.divisaderoProjectionMinimum
1099 else:
1100 minRow = self.config.divisaderoProjectionMinimum
1101 maxRow = self.config.divisaderoProjectionMaximum
1103 segment = slice(minRow, maxRow)
1104 projection, _, _ = astropy.stats.sigma_clipped_stats(ampArray[segment, :], axis=0)
1106 ampStats = {}
1107 projection /= np.median(projection)
1108 columns = np.arange(len(projection))
1110 segment = slice(self.config.divisaderoEdgePixels, -self.config.divisaderoEdgePixels)
1111 model = np.polyfit(columns[segment], projection[segment], 1)
1112 modelProjection = model[0] * columns + model[1]
1113 divisaderoProfile = projection / modelProjection
1115 # look for max at the edges:
1116 leftMax = np.nanmax(np.abs(divisaderoProfile[0:self.config.divisaderoNumImpactPixels] - 1.0))
1117 rightMax = np.nanmax(np.abs(divisaderoProfile[-self.config.divisaderoNumImpactPixels:] - 1.0))
1119 ampStats['DIVISADERO_PROFILE'] = np.array(divisaderoProfile).tolist()
1120 ampStats['DIVISADERO_MAX_PAIR'] = [leftMax, rightMax]
1121 outputStats[amp.getName()] = ampStats
1123 detector = inputExp.getDetector()
1124 xCenters = [amp.getBBox().getCenterX() for amp in detector]
1125 yCenters = [amp.getBBox().getCenterY() for amp in detector]
1126 xIndices = np.ceil(xCenters / np.min(xCenters) / 2).astype(int) - 1
1127 yIndices = np.ceil(yCenters / np.min(yCenters) / 2).astype(int) - 1
1128 ampIds = np.zeros((len(set(yIndices)), len(set(xIndices))), dtype=int)
1129 for ampId, xIndex, yIndex in zip(np.arange(len(detector)), xIndices, yIndices):
1130 ampIds[yIndex, xIndex] = ampId
1132 # Loop over amps again because the DIVISIDERO_MAX will be the max
1133 # of the profile on its boundary with its neighboring amps
1134 for i, amp in enumerate(detector):
1135 y, x = np.where(ampIds == i)
1136 end = ampIds.shape[1] - 1
1137 xInd = x[0]
1138 yInd = y[0]
1139 thisAmpsPair = outputStats[amp.getName()]['DIVISADERO_MAX_PAIR']
1141 if x == 0:
1142 # leftmost amp: take the max of your right side and
1143 myMax = thisAmpsPair[1]
1144 # your neighbor's left side
1145 neighborMax = outputStats[detector[ampIds[yInd, 1]].getName()]['DIVISADERO_MAX_PAIR'][0]
1146 elif x == end:
1147 # rightmost amp: take the max of your left side and
1148 myMax = thisAmpsPair[0]
1149 # your neighbor's right side
1150 neighborMax = outputStats[detector[ampIds[yInd, end - 1]].getName()]['DIVISADERO_MAX_PAIR'][1]
1151 else:
1152 # Middle amp: take the max of both your own sides and the
1153 myMax = max(thisAmpsPair)
1154 leftName = detector[ampIds[yInd, max(xInd - 1, 0)]].getName()
1155 rightName = detector[ampIds[yInd, min(xInd + 1, ampIds.shape[1] - 1)]].getName()
1156 # right side of the neighbor to your left
1157 # and left side of your neighbor to your right
1158 neighborMax = max(outputStats[leftName]['DIVISADERO_MAX_PAIR'][1],
1159 outputStats[rightName]['DIVISADERO_MAX_PAIR'][0])
1161 divisaderoMax = max([myMax, neighborMax])
1162 outputStats[amp.getName()]['DIVISADERO_MAX'] = divisaderoMax
1164 return outputStats