LSST Applications 26.0.0,g0265f82a02+6660c170cc,g07994bdeae+30b05a742e,g0a0026dc87+17526d298f,g0a60f58ba1+17526d298f,g0e4bf8285c+96dd2c2ea9,g0ecae5effc+c266a536c8,g1e7d6db67d+6f7cb1f4bb,g26482f50c6+6346c0633c,g2bbee38e9b+6660c170cc,g2cc88a2952+0a4e78cd49,g3273194fdb+f6908454ef,g337abbeb29+6660c170cc,g337c41fc51+9a8f8f0815,g37c6e7c3d5+7bbafe9d37,g44018dc512+6660c170cc,g4a941329ef+4f7594a38e,g4c90b7bd52+5145c320d2,g58be5f913a+bea990ba40,g635b316a6c+8d6b3a3e56,g67924a670a+bfead8c487,g6ae5381d9b+81bc2a20b4,g93c4d6e787+26b17396bd,g98cecbdb62+ed2cb6d659,g98ffbb4407+81bc2a20b4,g9ddcbc5298+7f7571301f,ga1e77700b3+99e9273977,gae46bcf261+6660c170cc,gb2715bf1a1+17526d298f,gc86a011abf+17526d298f,gcf0d15dbbd+96dd2c2ea9,gdaeeff99f8+0d8dbea60f,gdb4ec4c597+6660c170cc,ge23793e450+96dd2c2ea9,gf041782ebf+171108ac67
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
25
26from scipy.signal.windows import hamming, hann, gaussian
27
28import lsst.afw.math as afwMath
29import lsst.afw.image as afwImage
30import lsst.pipe.base as pipeBase
31import lsst.pex.config as pexConfig
32
33from lsst.afw.cameraGeom import ReadoutCorner
34
35
36class IsrStatisticsTaskConfig(pexConfig.Config):
37 """Image statistics options.
38 """
39 doCtiStatistics = pexConfig.Field(
40 dtype=bool,
41 doc="Measure CTI statistics from image and overscans?",
42 default=False,
43 )
44
45 doBandingStatistics = pexConfig.Field(
46 dtype=bool,
47 doc="Measure image banding metric?",
48 default=False,
49 )
50 bandingKernelSize = pexConfig.Field(
51 dtype=int,
52 doc="Width of box for boxcar smoothing for banding metric.",
53 default=3,
54 check=lambda x: x == 0 or x % 2 != 0,
55 )
56 bandingFractionLow = pexConfig.Field(
57 dtype=float,
58 doc="Fraction of values to exclude from low samples.",
59 default=0.1,
60 check=lambda x: x >= 0.0 and x <= 1.0
61 )
62 bandingFractionHigh = pexConfig.Field(
63 dtype=float,
64 doc="Fraction of values to exclude from high samples.",
65 default=0.9,
66 check=lambda x: x >= 0.0 and x <= 1.0,
67 )
68 bandingUseHalfDetector = pexConfig.Field(
69 dtype=float,
70 doc="Use only the first half set of amplifiers.",
71 default=True,
72 )
73
74 doProjectionStatistics = pexConfig.Field(
75 dtype=bool,
76 doc="Measure projection metric?",
77 default=False,
78 )
79 projectionKernelSize = pexConfig.Field(
80 dtype=int,
81 doc="Width of box for boxcar smoothing of projections.",
82 default=0,
83 check=lambda x: x == 0 or x % 2 != 0,
84 )
85 doProjectionFft = pexConfig.Field(
86 dtype=bool,
87 doc="Generate FFTs from the image projections?",
88 default=False,
89 )
90 projectionFftWindow = pexConfig.ChoiceField(
91 dtype=str,
92 doc="Type of windowing to use prior to calculating FFT.",
93 default="HAMMING",
94 allowed={
95 "HAMMING": "Hamming window.",
96 "HANN": "Hann window.",
97 "GAUSSIAN": "Gaussian window.",
98 "NONE": "No window."
99 }
100 )
101
102 stat = pexConfig.Field(
103 dtype=str,
104 default='MEANCLIP',
105 doc="Statistic name to use to measure regions.",
106 )
107 nSigmaClip = pexConfig.Field(
108 dtype=float,
109 default=3.0,
110 doc="Clipping threshold for background",
111 )
112 nIter = pexConfig.Field(
113 dtype=int,
114 default=3,
115 doc="Clipping iterations for background",
116 )
117 badMask = pexConfig.ListField(
118 dtype=str,
119 default=["BAD", "INTRP", "SAT"],
120 doc="Mask planes to ignore when identifying source pixels."
121 )
122
123
124class IsrStatisticsTask(pipeBase.Task):
125 """Task to measure arbitrary statistics on ISR processed exposures.
126
127 The goal is to wrap a number of optional measurements that are
128 useful for calibration production and detector stability.
129 """
130 ConfigClass = IsrStatisticsTaskConfig
131 _DefaultName = "isrStatistics"
132
133 def __init__(self, statControl=None, **kwargs):
134 super().__init__(**kwargs)
135 self.statControl = afwMath.StatisticsControl(self.config.nSigmaClip, self.config.nIter,
136 afwImage.Mask.getPlaneBitMask(self.config.badMask))
138
139 def run(self, inputExp, ptc=None, overscanResults=None, **kwargs):
140 """Task to run arbitrary statistics.
141
142 The statistics should be measured by individual methods, and
143 add to the dictionary in the return struct.
144
145 Parameters
146 ----------
147 inputExp : `lsst.afw.image.Exposure`
148 The exposure to measure.
149 ptc : `lsst.ip.isr.PtcDataset`, optional
150 A PTC object containing gains to use.
151 overscanResults : `list` [`lsst.pipe.base.Struct`], optional
152 List of overscan results. Expected fields are:
153
154 ``imageFit``
155 Value or fit subtracted from the amplifier image data
156 (scalar or `lsst.afw.image.Image`).
157 ``overscanFit``
158 Value or fit subtracted from the overscan image data
159 (scalar or `lsst.afw.image.Image`).
160 ``overscanImage``
161 Image of the overscan region with the overscan
162 correction applied (`lsst.afw.image.Image`). This
163 quantity is used to estimate the amplifier read noise
164 empirically.
165
166 Returns
167 -------
168 resultStruct : `lsst.pipe.base.Struct`
169 Contains the measured statistics as a dict stored in a
170 field named ``results``.
171
172 Raises
173 ------
174 RuntimeError
175 Raised if the amplifier gains could not be found.
176 """
177 # Find gains.
178 detector = inputExp.getDetector()
179 if ptc is not None:
180 gains = ptc.gain
181 elif detector is not None:
182 gains = {amp.getName(): amp.getGain() for amp in detector.getAmplifiers()}
183 else:
184 raise RuntimeError("No source of gains provided.")
185
186 ctiResults = None
187 if self.config.doCtiStatistics:
188 ctiResults = self.measureCti(inputExp, overscanResults, gains)
189
190 bandingResults = None
191 if self.config.doBandingStatistics:
192 bandingResults = self.measureBanding(inputExp, overscanResults)
193
194 projectionResults = None
195 if self.config.doProjectionStatistics:
196 projectionResults = self.measureProjectionStatistics(inputExp, overscanResults)
197
198 return pipeBase.Struct(
199 results={'CTI': ctiResults,
200 'BANDING': bandingResults,
201 'PROJECTION': projectionResults,
202 },
203 )
204
205 def measureCti(self, inputExp, overscans, gains):
206 """Task to measure CTI statistics.
207
208 Parameters
209 ----------
210 inputExp : `lsst.afw.image.Exposure`
211 Exposure to measure.
212 overscans : `list` [`lsst.pipe.base.Struct`]
213 List of overscan results. Expected fields are:
214
215 ``imageFit``
216 Value or fit subtracted from the amplifier image data
217 (scalar or `lsst.afw.image.Image`).
218 ``overscanFit``
219 Value or fit subtracted from the overscan image data
220 (scalar or `lsst.afw.image.Image`).
221 ``overscanImage``
222 Image of the overscan region with the overscan
223 correction applied (`lsst.afw.image.Image`). This
224 quantity is used to estimate the amplifier read noise
225 empirically.
226 gains : `dict` [`str` `float`]
227 Dictionary of per-amplifier gains, indexed by amplifier name.
228
229 Returns
230 -------
231 outputStats : `dict` [`str`, [`dict` [`str`,`float]]
232 Dictionary of measurements, keyed by amplifier name and
233 statistics segment.
234 """
235 outputStats = {}
236
237 detector = inputExp.getDetector()
238 image = inputExp.image
239
240 # Ensure we have the same number of overscans as amplifiers.
241 assert len(overscans) == len(detector.getAmplifiers())
242
243 for ampIter, amp in enumerate(detector.getAmplifiers()):
244 ampStats = {}
245 gain = gains[amp.getName()]
246 readoutCorner = amp.getReadoutCorner()
247 # Full data region.
248 dataRegion = image[amp.getBBox()]
249 ampStats['IMAGE_MEAN'] = afwMath.makeStatistics(dataRegion, self.statType,
250 self.statControl).getValue()
251
252 # First and last image columns.
253 pixelA = afwMath.makeStatistics(dataRegion.array[:, 0],
254 self.statType,
255 self.statControl).getValue()
256 pixelZ = afwMath.makeStatistics(dataRegion.array[:, -1],
257 self.statType,
258 self.statControl).getValue()
259
260 # We want these relative to the readout corner. If that's
261 # on the right side, we need to swap them.
262 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
263 ampStats['FIRST_MEAN'] = pixelZ
264 ampStats['LAST_MEAN'] = pixelA
265 else:
266 ampStats['FIRST_MEAN'] = pixelA
267 ampStats['LAST_MEAN'] = pixelZ
268
269 # Measure the columns of the overscan.
270 if overscans[ampIter] is None:
271 # The amplifier is likely entirely bad, and needs to
272 # be skipped.
273 self.log.warn("No overscan information available for ISR statistics for amp %s.",
274 amp.getName())
275 nCols = amp.getSerialOverscanBBox().getWidth()
276 ampStats['OVERSCAN_COLUMNS'] = np.full((nCols, ), np.nan)
277 ampStats['OVERSCAN_VALUES'] = np.full((nCols, ), np.nan)
278 else:
279 overscanImage = overscans[ampIter].overscanImage
280 columns = []
281 values = []
282 for column in range(0, overscanImage.getWidth()):
283 osMean = afwMath.makeStatistics(overscanImage.image.array[:, column],
284 self.statType, self.statControl).getValue()
285 columns.append(column)
286 values.append(gain * osMean)
287
288 # We want these relative to the readout corner. If that's
289 # on the right side, we need to swap them.
290 if readoutCorner in (ReadoutCorner.LR, ReadoutCorner.UR):
291 ampStats['OVERSCAN_COLUMNS'] = list(reversed(columns))
292 ampStats['OVERSCAN_VALUES'] = list(reversed(values))
293 else:
294 ampStats['OVERSCAN_COLUMNS'] = columns
295 ampStats['OVERSCAN_VALUES'] = values
296
297 outputStats[amp.getName()] = ampStats
298
299 return outputStats
300
301 @staticmethod
302 def makeKernel(kernelSize):
303 """Make a boxcar smoothing kernel.
304
305 Parameters
306 ----------
307 kernelSize : `int`
308 Size of the kernel in pixels.
309
310 Returns
311 -------
312 kernel : `np.array`
313 Kernel for boxcar smoothing.
314 """
315 if kernelSize > 0:
316 kernel = np.full(kernelSize, 1.0 / kernelSize)
317 else:
318 kernel = np.array([1.0])
319 return kernel
320
321 def measureBanding(self, inputExp, overscans):
322 """Task to measure banding statistics.
323
324 Parameters
325 ----------
326 inputExp : `lsst.afw.image.Exposure`
327 Exposure to measure.
328 overscans : `list` [`lsst.pipe.base.Struct`]
329 List of overscan results. Expected fields are:
330
331 ``imageFit``
332 Value or fit subtracted from the amplifier image data
333 (scalar or `lsst.afw.image.Image`).
334 ``overscanFit``
335 Value or fit subtracted from the overscan image data
336 (scalar or `lsst.afw.image.Image`).
337 ``overscanImage``
338 Image of the overscan region with the overscan
339 correction applied (`lsst.afw.image.Image`). This
340 quantity is used to estimate the amplifier read noise
341 empirically.
342
343 Returns
344 -------
345 outputStats : `dict` [`str`, [`dict` [`str`,`float]]
346 Dictionary of measurements, keyed by amplifier name and
347 statistics segment.
348 """
349 outputStats = {}
350
351 detector = inputExp.getDetector()
352 kernel = self.makeKernel(self.config.bandingKernelSize)
353
354 outputStats['AMP_BANDING'] = []
355 for amp, overscanData in zip(detector.getAmplifiers(), overscans):
356 overscanFit = np.array(overscanData.overscanFit)
357 overscanArray = overscanData.overscanImage.image.array
358 rawOverscan = np.mean(overscanArray + overscanFit, axis=1)
359
360 smoothedOverscan = np.convolve(rawOverscan, kernel, mode='valid')
361
362 low, high = np.quantile(smoothedOverscan, [self.config.bandingFractionLow,
363 self.config.bandingFractionHigh])
364 outputStats['AMP_BANDING'].append(float(high - low))
365
366 if self.config.bandingUseHalfDetector:
367 fullLength = len(outputStats['AMP_BANDING'])
368 outputStats['DET_BANDING'] = float(np.nanmedian(outputStats['AMP_BANDING'][0:fullLength//2]))
369 else:
370 outputStats['DET_BANDING'] = float(np.nanmedian(outputStats['AMP_BANDING']))
371
372 return outputStats
373
374 def measureProjectionStatistics(self, inputExp, overscans):
375 """Task to measure metrics from image slicing.
376
377 Parameters
378 ----------
379 inputExp : `lsst.afw.image.Exposure`
380 Exposure to measure.
381 overscans : `list` [`lsst.pipe.base.Struct`]
382 List of overscan results. Expected fields are:
383
384 ``imageFit``
385 Value or fit subtracted from the amplifier image data
386 (scalar or `lsst.afw.image.Image`).
387 ``overscanFit``
388 Value or fit subtracted from the overscan image data
389 (scalar or `lsst.afw.image.Image`).
390 ``overscanImage``
391 Image of the overscan region with the overscan
392 correction applied (`lsst.afw.image.Image`). This
393 quantity is used to estimate the amplifier read noise
394 empirically.
395
396 Returns
397 -------
398 outputStats : `dict` [`str`, [`dict` [`str`,`float]]
399 Dictionary of measurements, keyed by amplifier name and
400 statistics segment.
401 """
402 outputStats = {}
403
404 detector = inputExp.getDetector()
405 kernel = self.makeKernel(self.config.projectionKernelSize)
406
407 outputStats['AMP_VPROJECTION'] = {}
408 outputStats['AMP_HPROJECTION'] = {}
409 convolveMode = 'valid'
410 if self.config.doProjectionFft:
411 outputStats['AMP_VFFT_REAL'] = {}
412 outputStats['AMP_VFFT_IMAG'] = {}
413 outputStats['AMP_HFFT_REAL'] = {}
414 outputStats['AMP_HFFT_IMAG'] = {}
415 convolveMode = 'same'
416
417 for amp in detector.getAmplifiers():
418 ampArray = inputExp.image[amp.getBBox()].array
419
420 horizontalProjection = np.mean(ampArray, axis=0)
421 verticalProjection = np.mean(ampArray, axis=1)
422
423 horizontalProjection = np.convolve(horizontalProjection, kernel, mode=convolveMode)
424 verticalProjection = np.convolve(verticalProjection, kernel, mode=convolveMode)
425
426 outputStats['AMP_HPROJECTION'][amp.getName()] = horizontalProjection.tolist()
427 outputStats['AMP_VPROJECTION'][amp.getName()] = verticalProjection.tolist()
428
429 if self.config.doProjectionFft:
430 horizontalWindow = np.ones_like(horizontalProjection)
431 verticalWindow = np.ones_like(verticalProjection)
432 if self.config.projectionFftWindow == "NONE":
433 pass
434 elif self.config.projectionFftWindow == "HAMMING":
435 horizontalWindow = hamming(len(horizontalProjection))
436 verticalWindow = hamming(len(verticalProjection))
437 elif self.config.projectionFftWindow == "HANN":
438 horizontalWindow = hann(len(horizontalProjection))
439 verticalWindow = hann(len(verticalProjection))
440 elif self.config.projectionFftWindow == "GAUSSIAN":
441 horizontalWindow = gaussian(len(horizontalProjection))
442 verticalWindow = gaussian(len(verticalProjection))
443 else:
444 raise RuntimeError(f"Invalid window function: {self.config.projectionFftWindow}")
445
446 horizontalFFT = np.fft.rfft(np.multiply(horizontalProjection, horizontalWindow))
447 verticalFFT = np.fft.rfft(np.multiply(verticalProjection, verticalWindow))
448 outputStats['AMP_HFFT_REAL'][amp.getName()] = np.real(horizontalFFT).tolist()
449 outputStats['AMP_HFFT_IMAG'][amp.getName()] = np.imag(horizontalFFT).tolist()
450 outputStats['AMP_VFFT_REAL'][amp.getName()] = np.real(verticalFFT).tolist()
451 outputStats['AMP_VFFT_IMAG'][amp.getName()] = np.imag(verticalFFT).tolist()
452
453 return outputStats
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition Exposure.h:72
A class to represent a 2-dimensional array of pixels.
Definition Image.h:51
Pass parameters to a Statistics object.
Definition Statistics.h:83
measureCti(self, inputExp, overscans, gains)
measureProjectionStatistics(self, inputExp, overscans)
__init__(self, statControl=None, **kwargs)
daf::base::PropertyList * list
Definition fits.cc:928
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)