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