LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
isrStatistics.py
Go to the documentation of this file.
1# This file is part of ip_isr.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22__all__ = ["IsrStatisticsTaskConfig", "IsrStatisticsTask"]
23
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
29
30import lsst.afw.math as afwMath
31import lsst.afw.image as afwImage
32import lsst.pipe.base as pipeBase
33import lsst.pex.config as pexConfig
34
35from lsst.afw.cameraGeom import ReadoutCorner
36
37
38class IsrStatisticsTaskConfig(pexConfig.Config):
39 """Image statistics options.
40 """
41 doCtiStatistics = pexConfig.Field(
42 dtype=bool,
43 doc="Measure CTI statistics from image and overscans?",
44 default=False,
45 )
46 doApplyGainsForCtiStatistics = pexConfig.Field(
47 dtype=bool,
48 doc="Apply gain to the overscan region when measuring CTI statistics?",
49 default=True,
50 )
51
52 doBandingStatistics = pexConfig.Field(
53 dtype=bool,
54 doc="Measure image banding metric?",
55 default=False,
56 )
57 bandingKernelSize = pexConfig.Field(
58 dtype=int,
59 doc="Width of box for boxcar smoothing for banding metric.",
60 default=3,
61 check=lambda x: x == 0 or x % 2 != 0,
62 )
63 bandingFractionLow = pexConfig.Field(
64 dtype=float,
65 doc="Fraction of values to exclude from low samples.",
66 default=0.1,
67 check=lambda x: x >= 0.0 and x <= 1.0
68 )
69 bandingFractionHigh = pexConfig.Field(
70 dtype=float,
71 doc="Fraction of values to exclude from high samples.",
72 default=0.9,
73 check=lambda x: x >= 0.0 and x <= 1.0,
74 )
75 bandingUseHalfDetector = pexConfig.Field(
76 dtype=float,
77 doc="Use only the first half set of amplifiers.",
78 default=True,
79 )
80
81 doProjectionStatistics = pexConfig.Field(
82 dtype=bool,
83 doc="Measure projection metric?",
84 default=False,
85 )
86 projectionKernelSize = pexConfig.Field(
87 dtype=int,
88 doc="Width of box for boxcar smoothing of projections.",
89 default=0,
90 check=lambda x: x == 0 or x % 2 != 0,
91 )
92 doProjectionFft = pexConfig.Field(
93 dtype=bool,
94 doc="Generate FFTs from the image projections?",
95 default=False,
96 )
97 projectionFftWindow = pexConfig.ChoiceField(
98 dtype=str,
99 doc="Type of windowing to use prior to calculating FFT.",
100 default="HAMMING",
101 allowed={
102 "HAMMING": "Hamming window.",
103 "HANN": "Hann window.",
104 "GAUSSIAN": "Gaussian window.",
105 "NONE": "No window."
106 }
107 )
108 doDivisaderoStatistics = pexConfig.Field(
109 dtype=bool,
110 doc="Measure divisadero tearing statistics?",
111 default=False,
112 )
113 divisaderoEdgePixels = pexConfig.Field(
114 dtype=int,
115 doc="Number of edge pixels excluded from divisadero linear fit.",
116 default=25,
117 )
118 divisaderoNumImpactPixels = pexConfig.Field(
119 dtype=int,
120 doc="Number of edge pixels to examine for divisadero tearing.",
121 default=2,
122 )
123 divisaderoProjectionMinimum = pexConfig.Field(
124 dtype=int,
125 doc="Minimum row to consider when taking robust mean of columns.",
126 default=10,
127 )
128 divisaderoProjectionMaximum = pexConfig.Field(
129 dtype=int,
130 doc="Maximum row to consider when taking robust mean of columns",
131 default=210,
132 )
133 doCopyCalibDistributionStatistics = pexConfig.Field(
134 dtype=bool,
135 doc="Copy calibration distribution statistics to output?",
136 default=False,
137 )
138 expectedDistributionLevels = pexConfig.ListField(
139 dtype=float,
140 doc="Percentile levels expected in the calibration header.",
141 default=[0, 5, 16, 50, 84, 95, 100],
142 )
143
144 doBiasShiftStatistics = pexConfig.Field(
145 dtype=bool,
146 doc="Measure number of image shifts in overscan?",
147 default=False,
148 )
149 biasShiftFilterOrder = pexConfig.Field(
150 dtype=int,
151 doc="Filter order for Butterworth highpass filter.",
152 default=5,
153 )
154 biasShiftCutoff = pexConfig.Field(
155 dtype=float,
156 doc="Cutoff frequency for highpass filter.",
157 default=1.0/15.0,
158 )
159 biasShiftWindow = pexConfig.Field(
160 dtype=int,
161 doc="Filter window size in pixels for highpass filter.",
162 default=30,
163 )
164 biasShiftThreshold = pexConfig.Field(
165 dtype=float,
166 doc="S/N threshold for bias shift detection.",
167 default=3.0,
168 )
169 biasShiftRowSkip = pexConfig.Field(
170 dtype=int,
171 doc="Number of rows to skip for the bias shift detection.",
172 default=30,
173 )
174 biasShiftColumnSkip = pexConfig.Field(
175 dtype=int,
176 doc="Number of columns to skip when averaging the overscan region.",
177 default=3,
178 )
179
180 doAmplifierCorrelationStatistics = pexConfig.Field(
181 dtype=bool,
182 doc="Measure amplifier correlations?",
183 default=False,
184 )
185
186 stat = pexConfig.Field(
187 dtype=str,
188 default="MEANCLIP",
189 doc="Statistic name to use to measure regions.",
190 )
191 nSigmaClip = pexConfig.Field(
192 dtype=float,
193 default=3.0,
194 doc="Clipping threshold for background",
195 )
196 nIter = pexConfig.Field(
197 dtype=int,
198 default=3,
199 doc="Clipping iterations for background",
200 )
201 badMask = pexConfig.ListField(
202 dtype=str,
203 default=["BAD", "INTRP", "SAT"],
204 doc="Mask planes to ignore when identifying source pixels."
205 )
206
207
208class IsrStatisticsTask(pipeBase.Task):
209 """Task to measure arbitrary statistics on ISR processed exposures.
210
211 The goal is to wrap a number of optional measurements that are
212 useful for calibration production and detector stability.
213 """
214 ConfigClass = IsrStatisticsTaskConfig
215 _DefaultName = "isrStatistics"
216
217 def __init__(self, statControl=None, **kwargs):
218 super().__init__(**kwargs)
219 self.statControl = afwMath.StatisticsControl(self.config.nSigmaClip, self.config.nIter,
220 afwImage.Mask.getPlaneBitMask(self.config.badMask))
222
223 def run(self, inputExp, ptc=None, overscanResults=None, **kwargs):
224 """Task to run arbitrary statistics.
225
226 The statistics should be measured by individual methods, and
227 add to the dictionary in the return struct.
228
229 Parameters
230 ----------
231 inputExp : `lsst.afw.image.Exposure`
232 The exposure to measure.
233 ptc : `lsst.ip.isr.PtcDataset`, optional
234 A PTC object containing gains to use.
235 overscanResults : `list` [`lsst.pipe.base.Struct`], optional
236 List of overscan results. Expected fields are:
237
238 ``imageFit``
239 Value or fit subtracted from the amplifier image data
240 (scalar or `lsst.afw.image.Image`).
241 ``overscanFit``
242 Value or fit subtracted from the overscan image data
243 (scalar or `lsst.afw.image.Image`).
244 ``overscanImage``
245 Image of the overscan region with the overscan
246 correction applied (`lsst.afw.image.Image`). This
247 quantity is used to estimate the amplifier read noise
248 empirically.
249 **kwargs :
250 Keyword arguments. Calibrations being passed in should
251 have an entry here.
252
253 Returns
254 -------
255 resultStruct : `lsst.pipe.base.Struct`
256 Contains the measured statistics as a dict stored in a
257 field named ``results``.
258
259 Raises
260 ------
261 RuntimeError
262 Raised if the amplifier gains could not be found.
263 """
264 # Find gains.
265 detector = inputExp.getDetector()
266 if ptc is not None:
267 gains = ptc.gain
268 elif detector is not None:
269 gains = {amp.getName(): amp.getGain() for amp in detector.getAmplifiers()}
270 else:
271 raise RuntimeError("No source of gains provided.")
272
273 ctiResults = None
274 if self.config.doCtiStatistics:
275 ctiResults = self.measureCti(inputExp, overscanResults, gains)
276
277 bandingResults = None
278 if self.config.doBandingStatistics:
279 bandingResults = self.measureBanding(inputExp, overscanResults)
280
281 projectionResults = None
282 if self.config.doProjectionStatistics:
283 projectionResults = self.measureProjectionStatistics(inputExp, overscanResults)
284
285 divisaderoResults = None
286 if self.config.doDivisaderoStatistics:
287 divisaderoResults = self.measureDivisaderoStatistics(inputExp, **kwargs)
288
289 calibDistributionResults = None
290 if self.config.doCopyCalibDistributionStatistics:
291 calibDistributionResults = self.copyCalibDistributionStatistics(inputExp, **kwargs)
292
293 biasShiftResults = None
294 if self.config.doBiasShiftStatistics:
295 biasShiftResults = self.measureBiasShifts(inputExp, overscanResults)
296
297 ampCorrelationResults = None
298 if self.config.doAmplifierCorrelationStatistics:
299 ampCorrelationResults = self.measureAmpCorrelations(inputExp, overscanResults)
300
301 mjd = inputExp.getMetadata().get("MJD", None)
302
303 return pipeBase.Struct(
304 results={"CTI": ctiResults,
305 "BANDING": bandingResults,
306 "PROJECTION": projectionResults,
307 "CALIBDIST": calibDistributionResults,
308 "BIASSHIFT": biasShiftResults,
309 "AMPCORR": ampCorrelationResults,
310 "MJD": mjd,
311 'DIVISADERO': divisaderoResults,
312 },
313 )
314
315 def measureCti(self, inputExp, overscans, gains):
316 """Task to measure CTI statistics.
317
318 Parameters
319 ----------
320 inputExp : `lsst.afw.image.Exposure`
321 Exposure to measure.
322 overscans : `list` [`lsst.pipe.base.Struct`]
323 List of overscan results. Expected fields are:
324
325 ``imageFit``
326 Value or fit subtracted from the amplifier image data
327 (scalar or `lsst.afw.image.Image`).
328 ``overscanFit``
329 Value or fit subtracted from the overscan image data
330 (scalar or `lsst.afw.image.Image`).
331 ``overscanImage``
332 Image of the overscan region with the overscan
333 correction applied (`lsst.afw.image.Image`). This
334 quantity is used to estimate the amplifier read noise
335 empirically.
336 gains : `dict` [`str` `float`]
337 Dictionary of per-amplifier gains, indexed by amplifier name.
338
339 Returns
340 -------
341 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
342 Dictionary of measurements, keyed by amplifier name and
343 statistics segment.
344 """
345 outputStats = {}
346
347 detector = inputExp.getDetector()
348 image = inputExp.image
349
350 # Ensure we have the same number of overscans as amplifiers.
351 assert len(overscans) == len(detector.getAmplifiers())
352
353 for ampIter, amp in enumerate(detector.getAmplifiers()):
354 ampStats = {}
355 gain = gains[amp.getName()]
356 readoutCorner = amp.getReadoutCorner()
357 # Full data region.
358 dataRegion = image[amp.getBBox()]
359 ampStats["IMAGE_MEAN"] = afwMath.makeStatistics(dataRegion, self.statType,
360 self.statControl).getValue()
361
362 # First and last image columns.
363 pixelA = afwMath.makeStatistics(dataRegion.array[:, 0],
364 self.statType,
365 self.statControl).getValue()
366 pixelZ = afwMath.makeStatistics(dataRegion.array[:, -1],
367 self.statType,
368 self.statControl).getValue()
369
370 # We want these relative to the readout corner. If that's
371 # on the right side, we need to swap them.
372 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
373 ampStats["FIRST_MEAN"] = pixelZ
374 ampStats["LAST_MEAN"] = pixelA
375 else:
376 ampStats["FIRST_MEAN"] = pixelA
377 ampStats["LAST_MEAN"] = pixelZ
378
379 # Measure the columns of the overscan.
380 if overscans[ampIter] is None:
381 # The amplifier is likely entirely bad, and needs to
382 # be skipped.
383 self.log.warning("No overscan information available for ISR statistics for amp %s.",
384 amp.getName())
385 nCols = amp.getRawSerialOverscanBBox().getWidth()
386 ampStats["OVERSCAN_COLUMNS"] = np.full((nCols, ), np.nan)
387 ampStats["OVERSCAN_VALUES"] = np.full((nCols, ), np.nan)
388 else:
389 overscanImage = overscans[ampIter].overscanImage
390 columns = []
391 values = []
392 for column in range(0, overscanImage.getWidth()):
393 # If overscan.doParallelOverscan=True, the overscanImage
394 # will contain both the serial and parallel overscan
395 # regions.
396 # Only the serial overscan correction is implemented,
397 # so we must select only the serial overscan rows
398 # for a given column.
399 nRows = amp.getRawSerialOverscanBBox().getHeight()
400 osMean = afwMath.makeStatistics(overscanImage.image.array[:nRows, column],
401 self.statType, self.statControl).getValue()
402 columns.append(column)
403 if self.config.doApplyGainsForCtiStatistics:
404 values.append(gain * osMean)
405 else:
406 values.append(osMean)
407
408 # We want these relative to the readout corner. If that's
409 # on the right side, we need to swap them.
410 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
411 ampStats["OVERSCAN_COLUMNS"] = list(reversed(columns))
412 ampStats["OVERSCAN_VALUES"] = list(reversed(values))
413 else:
414 ampStats["OVERSCAN_COLUMNS"] = columns
415 ampStats["OVERSCAN_VALUES"] = values
416
417 outputStats[amp.getName()] = ampStats
418
419 return outputStats
420
421 @staticmethod
422 def makeKernel(kernelSize):
423 """Make a boxcar smoothing kernel.
424
425 Parameters
426 ----------
427 kernelSize : `int`
428 Size of the kernel in pixels.
429
430 Returns
431 -------
432 kernel : `np.array`
433 Kernel for boxcar smoothing.
434 """
435 if kernelSize > 0:
436 kernel = np.full(kernelSize, 1.0 / kernelSize)
437 else:
438 kernel = np.array([1.0])
439 return kernel
440
441 def measureBanding(self, inputExp, overscans):
442 """Task to measure banding statistics.
443
444 Parameters
445 ----------
446 inputExp : `lsst.afw.image.Exposure`
447 Exposure to measure.
448 overscans : `list` [`lsst.pipe.base.Struct`]
449 List of overscan results. Expected fields are:
450
451 ``imageFit``
452 Value or fit subtracted from the amplifier image data
453 (scalar or `lsst.afw.image.Image`).
454 ``overscanFit``
455 Value or fit subtracted from the overscan image data
456 (scalar or `lsst.afw.image.Image`).
457 ``overscanImage``
458 Image of the overscan region with the overscan
459 correction applied (`lsst.afw.image.Image`). This
460 quantity is used to estimate the amplifier read noise
461 empirically.
462
463 Returns
464 -------
465 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
466 Dictionary of measurements, keyed by amplifier name and
467 statistics segment.
468 """
469 outputStats = {}
470
471 detector = inputExp.getDetector()
472 kernel = self.makeKernel(self.config.bandingKernelSize)
473
474 outputStats["AMP_BANDING"] = []
475 for amp, overscanData in zip(detector.getAmplifiers(), overscans):
476 overscanFit = np.array(overscanData.overscanFit)
477 overscanArray = overscanData.overscanImage.image.array
478 rawOverscan = np.mean(overscanArray + overscanFit, axis=1)
479
480 smoothedOverscan = np.convolve(rawOverscan, kernel, mode="valid")
481
482 low, high = np.quantile(smoothedOverscan, [self.config.bandingFractionLow,
483 self.config.bandingFractionHigh])
484 outputStats["AMP_BANDING"].append(float(high - low))
485
486 if self.config.bandingUseHalfDetector:
487 fullLength = len(outputStats["AMP_BANDING"])
488 outputStats["DET_BANDING"] = float(np.nanmedian(outputStats["AMP_BANDING"][0:fullLength//2]))
489 else:
490 outputStats["DET_BANDING"] = float(np.nanmedian(outputStats["AMP_BANDING"]))
491
492 return outputStats
493
494 def measureProjectionStatistics(self, inputExp, overscans):
495 """Task to measure metrics from image slicing.
496
497 Parameters
498 ----------
499 inputExp : `lsst.afw.image.Exposure`
500 Exposure to measure.
501 overscans : `list` [`lsst.pipe.base.Struct`]
502 List of overscan results. Expected fields are:
503
504 ``imageFit``
505 Value or fit subtracted from the amplifier image data
506 (scalar or `lsst.afw.image.Image`).
507 ``overscanFit``
508 Value or fit subtracted from the overscan image data
509 (scalar or `lsst.afw.image.Image`).
510 ``overscanImage``
511 Image of the overscan region with the overscan
512 correction applied (`lsst.afw.image.Image`). This
513 quantity is used to estimate the amplifier read noise
514 empirically.
515
516 Returns
517 -------
518 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
519 Dictionary of measurements, keyed by amplifier name and
520 statistics segment.
521 """
522 outputStats = {}
523
524 detector = inputExp.getDetector()
525 kernel = self.makeKernel(self.config.projectionKernelSize)
526
527 outputStats["AMP_VPROJECTION"] = {}
528 outputStats["AMP_HPROJECTION"] = {}
529 convolveMode = "valid"
530 if self.config.doProjectionFft:
531 outputStats["AMP_VFFT_REAL"] = {}
532 outputStats["AMP_VFFT_IMAG"] = {}
533 outputStats["AMP_HFFT_REAL"] = {}
534 outputStats["AMP_HFFT_IMAG"] = {}
535 convolveMode = "same"
536
537 for amp in detector.getAmplifiers():
538 ampArray = inputExp.image[amp.getBBox()].array
539
540 horizontalProjection = np.mean(ampArray, axis=0)
541 verticalProjection = np.mean(ampArray, axis=1)
542
543 horizontalProjection = np.convolve(horizontalProjection, kernel, mode=convolveMode)
544 verticalProjection = np.convolve(verticalProjection, kernel, mode=convolveMode)
545
546 outputStats["AMP_HPROJECTION"][amp.getName()] = horizontalProjection.tolist()
547 outputStats["AMP_VPROJECTION"][amp.getName()] = verticalProjection.tolist()
548
549 if self.config.doProjectionFft:
550 horizontalWindow = np.ones_like(horizontalProjection)
551 verticalWindow = np.ones_like(verticalProjection)
552 if self.config.projectionFftWindow == "NONE":
553 pass
554 elif self.config.projectionFftWindow == "HAMMING":
555 horizontalWindow = hamming(len(horizontalProjection))
556 verticalWindow = hamming(len(verticalProjection))
557 elif self.config.projectionFftWindow == "HANN":
558 horizontalWindow = hann(len(horizontalProjection))
559 verticalWindow = hann(len(verticalProjection))
560 elif self.config.projectionFftWindow == "GAUSSIAN":
561 horizontalWindow = gaussian(len(horizontalProjection))
562 verticalWindow = gaussian(len(verticalProjection))
563 else:
564 raise RuntimeError(f"Invalid window function: {self.config.projectionFftWindow}")
565
566 horizontalFFT = np.fft.rfft(np.multiply(horizontalProjection, horizontalWindow))
567 verticalFFT = np.fft.rfft(np.multiply(verticalProjection, verticalWindow))
568
569 outputStats["AMP_HFFT_REAL"][amp.getName()] = np.real(horizontalFFT).tolist()
570 outputStats["AMP_HFFT_IMAG"][amp.getName()] = np.imag(horizontalFFT).tolist()
571 outputStats["AMP_VFFT_REAL"][amp.getName()] = np.real(verticalFFT).tolist()
572 outputStats["AMP_VFFT_IMAG"][amp.getName()] = np.imag(verticalFFT).tolist()
573
574 return outputStats
575
576 def copyCalibDistributionStatistics(self, inputExp, **kwargs):
577 """Copy calibration statistics for this exposure.
578
579 Parameters
580 ----------
581 inputExp : `lsst.afw.image.Exposure`
582 The exposure being processed.
583 **kwargs :
584 Keyword arguments with calibrations.
585
586 Returns
587 -------
588 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
589 Dictionary of measurements, keyed by amplifier name and
590 statistics segment.
591 """
592 outputStats = {}
593
594 # Amp level elements
595 for amp in inputExp.getDetector():
596 ampStats = {}
597
598 for calibType in ("bias", "dark", "flat"):
599 if kwargs.get(calibType, None) is not None:
600 metadata = kwargs[calibType].getMetadata()
601 for pct in self.config.expectedDistributionLevels:
602 key = f"LSST CALIB {calibType.upper()} {amp.getName()} DISTRIBUTION {pct}-PCT"
603 ampStats[key] = metadata.get(key, np.nan)
604
605 for calibType in ("defects"):
606 if kwargs.get(calibType, None) is not None:
607 metadata = kwargs[calibType].getMetadata()
608 for key in (f"LSST CALIB {calibType.upper()} {amp.getName()} N_HOT",
609 f"LSST CALIB {calibType.upper()} {amp.getName()} N_COLD"):
610 ampStats[key] = metadata.get(key, np.nan)
611 outputStats[amp.getName()] = ampStats
612
613 # Detector level elements
614 for calibType in ("defects"):
615 if kwargs.get(calibType, None) is not None:
616 metadata = kwargs[calibType].getMetadata()
617 for key in (f"LSST CALIB {calibType.upper()} N_BAD_COLUMNS"):
618 outputStats["detector"][key] = metadata.get(key, np.nan)
619
620 return outputStats
621
622 def measureBiasShifts(self, inputExp, overscanResults):
623 """Measure number of bias shifts from overscan data.
624
625 Parameters
626 ----------
627 inputExp : `lsst.afw.image.Exposure`
628 Exposure to measure.
629 overscans : `list` [`lsst.pipe.base.Struct`]
630 List of overscan results. Expected fields are:
631
632 ``imageFit``
633 Value or fit subtracted from the amplifier image data
634 (scalar or `lsst.afw.image.Image`).
635 ``overscanFit``
636 Value or fit subtracted from the overscan image data
637 (scalar or `lsst.afw.image.Image`).
638 ``overscanImage``
639 Image of the overscan region with the overscan
640 correction applied (`lsst.afw.image.Image`). This
641 quantity is used to estimate the amplifier read noise
642 empirically.
643
644 Returns
645 -------
646 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
647 Dictionary of measurements, keyed by amplifier name and
648 statistics segment.
649
650 Notes
651 -----
652 Based on eop_pipe implementation:
653 https://github.com/lsst-camera-dh/eo_pipe/blob/main/python/lsst/eo/pipe/biasShiftsTask.py # noqa: E501 W505
654 """
655 outputStats = {}
656
657 detector = inputExp.getDetector()
658 for amp, overscans in zip(detector, overscanResults):
659 ampStats = {}
660 # Add fit back to data
661 rawOverscan = overscans.overscanImage.image.array + overscans.overscanFit
662
663 # Collapse array, skipping first three columns
664 rawOverscan = np.mean(rawOverscan[:, self.config.biasShiftColumnSkip:], axis=1)
665
666 # Scan for shifts
667 noise, shift_peaks = self._scan_for_shifts(rawOverscan)
668 ampStats["LOCAL_NOISE"] = float(noise)
669 ampStats["BIAS_SHIFTS"] = shift_peaks
670
671 outputStats[amp.getName()] = ampStats
672 return outputStats
673
674 def _scan_for_shifts(self, overscanData):
675 """Scan overscan data for shifts.
676
677 Parameters
678 ----------
679 overscanData : `list` [`float`]
680 Overscan data to search for shifts.
681
682 Returns
683 -------
684 noise : `float`
685 Noise estimated from Butterworth filtered overscan data.
686 peaks : `list` [`float`, `float`, `int`, `int`]
687 Shift peak information, containing the convolved peak
688 value, the raw peak value, and the lower and upper bounds
689 of the region checked.
690 """
691 numerator, denominator = butter(self.config.biasShiftFilterOrder,
692 self.config.biasShiftCutoff,
693 btype="high", analog=False)
694 noise = np.std(filtfilt(numerator, denominator, overscanData))
695 kernel = np.concatenate([np.arange(self.config.biasShiftWindow),
696 np.arange(-self.config.biasShiftWindow + 1, 0)])
697 kernel = kernel/np.sum(kernel[:self.config.biasShiftWindow])
698
699 convolved = np.convolve(overscanData, kernel, mode="valid")
700 convolved = np.pad(convolved, (self.config.biasShiftWindow - 1, self.config.biasShiftWindow))
701
702 shift_check = np.abs(convolved)/noise
703 shift_mask = shift_check > self.config.biasShiftThreshold
704 shift_mask[:self.config.biasShiftRowSkip] = False
705
706 shift_regions = np.flatnonzero(np.diff(np.r_[np.int8(0),
707 shift_mask.view(np.int8),
708 np.int8(0)])).reshape(-1, 2)
709 shift_peaks = []
710 for region in shift_regions:
711 region_peak = np.argmax(shift_check[region[0]:region[1]]) + region[0]
712 if self._satisfies_flatness(region_peak, convolved[region_peak], overscanData):
713 shift_peaks.append(
714 [float(convolved[region_peak]), float(region_peak),
715 int(region[0]), int(region[1])])
716 return noise, shift_peaks
717
718 def _satisfies_flatness(self, shiftRow, shiftPeak, overscanData):
719 """Determine if a region is flat.
720
721 Parameters
722 ----------
723 shiftRow : `int`
724 Row with possible peak.
725 shiftPeak : `float`
726 Value at the possible peak.
727 overscanData : `list` [`float`]
728 Overscan data used to fit around the possible peak.
729
730 Returns
731 -------
732 isFlat : `bool`
733 Indicates if the region is flat, and so the peak is valid.
734 """
735 prerange = np.arange(shiftRow - self.config.biasShiftWindow, shiftRow)
736 postrange = np.arange(shiftRow, shiftRow + self.config.biasShiftWindow)
737
738 preFit = linregress(prerange, overscanData[prerange])
739 postFit = linregress(postrange, overscanData[postrange])
740
741 if shiftPeak > 0:
742 preTrend = (2*preFit[0]*len(prerange) < shiftPeak)
743 postTrend = (2*postFit[0]*len(postrange) < shiftPeak)
744 else:
745 preTrend = (2*preFit[0]*len(prerange) > shiftPeak)
746 postTrend = (2*postFit[0]*len(postrange) > shiftPeak)
747
748 return (preTrend and postTrend)
749
750 def measureAmpCorrelations(self, inputExp, overscanResults):
751 """Measure correlations between amplifier segments.
752
753 Parameters
754 ----------
755 inputExp : `lsst.afw.image.Exposure`
756 Exposure to measure.
757 overscans : `list` [`lsst.pipe.base.Struct`]
758 List of overscan results. Expected fields are:
759
760 ``imageFit``
761 Value or fit subtracted from the amplifier image data
762 (scalar or `lsst.afw.image.Image`).
763 ``overscanFit``
764 Value or fit subtracted from the overscan image data
765 (scalar or `lsst.afw.image.Image`).
766 ``overscanImage``
767 Image of the overscan region with the overscan
768 correction applied (`lsst.afw.image.Image`). This
769 quantity is used to estimate the amplifier read noise
770 empirically.
771
772 Returns
773 -------
774 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
775 Dictionary of measurements, keyed by amplifier name and
776 statistics segment.
777
778 Notes
779 -----
780 Based on eo_pipe implementation:
781 https://github.com/lsst-camera-dh/eo_pipe/blob/main/python/lsst/eo/pipe/raft_level_correlations.py # noqa: E501 W505
782 """
783 detector = inputExp.getDetector()
784 rows = []
785
786 for ampId, overscan in enumerate(overscanResults):
787 rawOverscan = overscan.overscanImage.image.array + overscan.overscanFit
788 rawOverscan = rawOverscan.ravel()
789
790 ampImage = inputExp[detector[ampId].getBBox()]
791 ampImage = ampImage.image.array.ravel()
792
793 for ampId2, overscan2 in enumerate(overscanResults):
794 osC = 1.0
795 imC = 1.0
796 if ampId2 != ampId:
797 rawOverscan2 = overscan2.overscanImage.image.array + overscan2.overscanFit
798 rawOverscan2 = rawOverscan2.ravel()
799
800 osC = np.corrcoef(rawOverscan, rawOverscan2)[0, 1]
801
802 ampImage2 = inputExp[detector[ampId2].getBBox()]
803 ampImage2 = ampImage2.image.array.ravel()
804
805 imC = np.corrcoef(ampImage, ampImage2)[0, 1]
806 rows.append(
807 {'detector': detector.getId(),
808 'detectorComp': detector.getId(),
809 'ampName': detector[ampId].getName(),
810 'ampComp': detector[ampId2].getName(),
811 'imageCorr': float(imC),
812 'overscanCorr': float(osC),
813 }
814 )
815 return rows
816
817 def measureDivisaderoStatistics(self, inputExp, **kwargs):
818 """Measure Max Divisadero Tearing effect per amp.
819
820 Parameters
821 ----------
822 inputExp : `lsst.afw.image.Exposure`
823 Exposure to measure. Usually a flat.
824 **kwargs :
825 The flat will be selected from here.
826
827 Returns
828 -------
829 outputStats : `dict` [`str`, [`dict` [`str`, `float`]]]
830 Dictionary of measurements, keyed by amplifier name and
831 statistics segment.
832 Measurements include
833
834 - DIVISADERO_PROFILE: Robust mean of rows between
835 divisaderoProjection<Maximum|Minumum> on readout edge of ccd
836 normalized by a linear fit to the same rows.
837 - DIVISADERO_MAX_PAIR: Tuple of maximum of the absolute values of
838 the DIVISADERO_PROFILE, for number of pixels (specified by
839 divisaderoNumImpactPixels on left and right side of amp.
840 - DIVISADERO_MAX: Maximum of the absolute values of the
841 the DIVISADERO_PROFILE, for the divisaderoNumImpactPixels on
842 boundaries of neighboring amps (including the pixels in those
843 neighborboring amps).
844 """
845 outputStats = {}
846
847 for amp in inputExp.getDetector():
848 # Copy unneeded if we do not ever modify the array by flipping
849 ampArray = inputExp.image[amp.getBBox()].array
850 # slice the outer top or bottom of the amp: the readout side
851 if amp.getReadoutCorner().name in ('UL', 'UR'):
852 minRow = amp.getBBox().getHeight() - self.config.divisaderoProjectionMaximum
853 maxRow = amp.getBBox().getHeight() - self.config.divisaderoProjectionMinimum
854 else:
855 minRow = self.config.divisaderoProjectionMinimum
856 maxRow = self.config.divisaderoProjectionMaximum
857
858 segment = slice(minRow, maxRow)
859 projection, _, _ = astropy.stats.sigma_clipped_stats(ampArray[segment, :], axis=0)
860
861 ampStats = {}
862 projection /= np.median(projection)
863 columns = np.arange(len(projection))
864
865 segment = slice(self.config.divisaderoEdgePixels, -self.config.divisaderoEdgePixels)
866 model = np.polyfit(columns[segment], projection[segment], 1)
867 modelProjection = model[0] * columns + model[1]
868 divisaderoProfile = projection / modelProjection
869
870 # look for max at the edges:
871 leftMax = np.nanmax(np.abs(divisaderoProfile[0:self.config.divisaderoNumImpactPixels] - 1.0))
872 rightMax = np.nanmax(np.abs(divisaderoProfile[-self.config.divisaderoNumImpactPixels:] - 1.0))
873
874 ampStats['DIVISADERO_PROFILE'] = np.array(divisaderoProfile).tolist()
875 ampStats['DIVISADERO_MAX_PAIR'] = [leftMax, rightMax]
876 outputStats[amp.getName()] = ampStats
877
878 detector = inputExp.getDetector()
879 xCenters = [amp.getBBox().getCenterX() for amp in detector]
880 yCenters = [amp.getBBox().getCenterY() for amp in detector]
881 xIndices = np.ceil(xCenters / np.min(xCenters) / 2).astype(int) - 1
882 yIndices = np.ceil(yCenters / np.min(yCenters) / 2).astype(int) - 1
883 ampIds = np.zeros((len(set(yIndices)), len(set(xIndices))), dtype=int)
884 for ampId, xIndex, yIndex in zip(np.arange(len(detector)), xIndices, yIndices):
885 ampIds[yIndex, xIndex] = ampId
886
887 # Loop over amps again because the DIVISIDERO_MAX will be the max
888 # of the profile on its boundary with its neighboring amps
889 for i, amp in enumerate(detector):
890 y, x = np.where(ampIds == i)
891 end = ampIds.shape[1] - 1
892 xInd = x[0]
893 yInd = y[0]
894 thisAmpsPair = outputStats[amp.getName()]['DIVISADERO_MAX_PAIR']
895
896 if x == 0:
897 # leftmost amp: take the max of your right side and
898 myMax = thisAmpsPair[1]
899 # your neighbor's left side
900 neighborMax = outputStats[detector[ampIds[yInd, 1]].getName()]['DIVISADERO_MAX_PAIR'][0]
901 elif x == end:
902 # rightmost amp: take the max of your left side and
903 myMax = thisAmpsPair[0]
904 # your neighbor's right side
905 neighborMax = outputStats[detector[ampIds[yInd, end - 1]].getName()]['DIVISADERO_MAX_PAIR'][1]
906 else:
907 # Middle amp: take the max of both your own sides and the
908 myMax = max(thisAmpsPair)
909 leftName = detector[ampIds[yInd, max(xInd - 1, 0)]].getName()
910 rightName = detector[ampIds[yInd, min(xInd + 1, ampIds.shape[1] - 1)]].getName()
911 # right side of the neighbor to your left
912 # and left side of your neighbor to your right
913 neighborMax = max(outputStats[leftName]['DIVISADERO_MAX_PAIR'][1],
914 outputStats[rightName]['DIVISADERO_MAX_PAIR'][0])
915
916 divisaderoMax = max([myMax, neighborMax])
917 outputStats[amp.getName()]['DIVISADERO_MAX'] = divisaderoMax
918
919 return outputStats
int min
int max
Pass parameters to a Statistics object.
Definition Statistics.h:83
measureDivisaderoStatistics(self, inputExp, **kwargs)
measureCti(self, inputExp, overscans, gains)
copyCalibDistributionStatistics(self, inputExp, **kwargs)
measureAmpCorrelations(self, inputExp, overscanResults)
measureProjectionStatistics(self, inputExp, overscans)
run(self, inputExp, ptc=None, overscanResults=None, **kwargs)
_satisfies_flatness(self, shiftRow, shiftPeak, overscanData)
measureBiasShifts(self, inputExp, overscanResults)
__init__(self, statControl=None, **kwargs)
daf::base::PropertySet * set
Definition fits.cc:931
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)