22__all__ = [
"IsrStatisticsTaskConfig",
"IsrStatisticsTask"]
26from scipy.signal.windows
import hamming, hann, gaussian
37 """Image statistics options.
39 doCtiStatistics = pexConfig.Field(
41 doc="Measure CTI statistics from image and overscans?",
45 doBandingStatistics = pexConfig.Field(
47 doc=
"Measure image banding metric?",
50 bandingKernelSize = pexConfig.Field(
52 doc=
"Width of box for boxcar smoothing for banding metric.",
54 check=
lambda x: x == 0
or x % 2 != 0,
56 bandingFractionLow = pexConfig.Field(
58 doc=
"Fraction of values to exclude from low samples.",
60 check=
lambda x: x >= 0.0
and x <= 1.0
62 bandingFractionHigh = pexConfig.Field(
64 doc=
"Fraction of values to exclude from high samples.",
66 check=
lambda x: x >= 0.0
and x <= 1.0,
68 bandingUseHalfDetector = pexConfig.Field(
70 doc=
"Use only the first half set of amplifiers.",
74 doProjectionStatistics = pexConfig.Field(
76 doc=
"Measure projection metric?",
79 projectionKernelSize = pexConfig.Field(
81 doc=
"Width of box for boxcar smoothing of projections.",
83 check=
lambda x: x == 0
or x % 2 != 0,
85 doProjectionFft = pexConfig.Field(
87 doc=
"Generate FFTs from the image projections?",
90 projectionFftWindow = pexConfig.ChoiceField(
92 doc=
"Type of windowing to use prior to calculating FFT.",
95 "HAMMING":
"Hamming window.",
96 "HANN":
"Hann window.",
97 "GAUSSIAN":
"Gaussian window.",
102 stat = pexConfig.Field(
105 doc=
"Statistic name to use to measure regions.",
107 nSigmaClip = pexConfig.Field(
110 doc=
"Clipping threshold for background",
112 nIter = pexConfig.Field(
115 doc=
"Clipping iterations for background",
117 badMask = pexConfig.ListField(
119 default=[
"BAD",
"INTRP",
"SAT"],
120 doc=
"Mask planes to ignore when identifying source pixels."
125 """Task to measure arbitrary statistics on ISR processed exposures.
127 The goal is to wrap a number of optional measurements that are
128 useful
for calibration production
and detector stability.
130 ConfigClass = IsrStatisticsTaskConfig
131 _DefaultName = "isrStatistics"
136 afwImage.Mask.getPlaneBitMask(self.config.badMask))
139 def run(self, inputExp, ptc=None, overscanResults=None, **kwargs):
140 """Task to run arbitrary statistics.
142 The statistics should be measured by individual methods, and
143 add to the dictionary
in the
return struct.
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:
155 Value
or fit subtracted
from the amplifier image data
158 Value
or fit subtracted
from the overscan image data
161 Image of the overscan region
with the overscan
163 quantity
is used to estimate the amplifier read noise
168 resultStruct : `lsst.pipe.base.Struct`
169 Contains the measured statistics
as a dict stored
in a
170 field named ``results``.
175 Raised
if the amplifier gains could
not be found.
178 detector = inputExp.getDetector()
181 elif detector
is not None:
182 gains = {amp.getName(): amp.getGain()
for amp
in detector.getAmplifiers()}
184 raise RuntimeError(
"No source of gains provided.")
187 if self.config.doCtiStatistics:
188 ctiResults = self.
measureCti(inputExp, overscanResults, gains)
190 bandingResults =
None
191 if self.config.doBandingStatistics:
194 projectionResults =
None
195 if self.config.doProjectionStatistics:
198 return pipeBase.Struct(
199 results={
'CTI': ctiResults,
200 'BANDING': bandingResults,
201 'PROJECTION': projectionResults,
206 """Task to measure CTI statistics.
212 overscans : `list` [`lsst.pipe.base.Struct`]
213 List of overscan results. Expected fields are:
216 Value or fit subtracted
from the amplifier image data
219 Value
or fit subtracted
from the overscan image data
222 Image of the overscan region
with the overscan
224 quantity
is used to estimate the amplifier read noise
226 gains : `dict` [`str` `float`]
227 Dictionary of per-amplifier gains, indexed by amplifier name.
231 outputStats : `dict` [`str`, [`dict` [`str`,`float]]
232 Dictionary of measurements, keyed by amplifier name
and
237 detector = inputExp.getDetector()
238 image = inputExp.image
241 assert len(overscans) == len(detector.getAmplifiers())
243 for ampIter, amp
in enumerate(detector.getAmplifiers()):
245 gain = gains[amp.getName()]
246 readoutCorner = amp.getReadoutCorner()
248 dataRegion = image[amp.getBBox()]
262 if readoutCorner
in (ReadoutCorner.LR, ReadoutCorner.UR):
263 ampStats[
'FIRST_MEAN'] = pixelZ
264 ampStats[
'LAST_MEAN'] = pixelA
266 ampStats[
'FIRST_MEAN'] = pixelA
267 ampStats[
'LAST_MEAN'] = pixelZ
270 if overscans[ampIter]
is None:
273 self.log.warn(
"No overscan information available for ISR statistics for amp %s.",
275 nCols = amp.getSerialOverscanBBox().getWidth()
276 ampStats[
'OVERSCAN_COLUMNS'] = np.full((nCols, ), np.nan)
277 ampStats[
'OVERSCAN_VALUES'] = np.full((nCols, ), np.nan)
279 overscanImage = overscans[ampIter].overscanImage
282 for column
in range(0, overscanImage.getWidth()):
285 columns.append(column)
286 values.append(gain * osMean)
290 if readoutCorner
in (ReadoutCorner.LR, ReadoutCorner.UR):
291 ampStats[
'OVERSCAN_COLUMNS'] =
list(reversed(columns))
292 ampStats[
'OVERSCAN_VALUES'] =
list(reversed(values))
294 ampStats[
'OVERSCAN_COLUMNS'] = columns
295 ampStats[
'OVERSCAN_VALUES'] = values
297 outputStats[amp.getName()] = ampStats
303 """Make a boxcar smoothing kernel.
308 Size of the kernel in pixels.
313 Kernel
for boxcar smoothing.
316 kernel = np.full(kernelSize, 1.0 / kernelSize)
318 kernel = np.array([1.0])
322 """Task to measure banding statistics.
328 overscans : `list` [`lsst.pipe.base.Struct`]
329 List of overscan results. Expected fields are:
332 Value or fit subtracted
from the amplifier image data
335 Value
or fit subtracted
from the overscan image data
338 Image of the overscan region
with the overscan
340 quantity
is used to estimate the amplifier read noise
345 outputStats : `dict` [`str`, [`dict` [`str`,`float]]
346 Dictionary of measurements, keyed by amplifier name
and
351 detector = inputExp.getDetector()
352 kernel = self.makeKernel(self.config.bandingKernelSize)
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)
360 smoothedOverscan = np.convolve(rawOverscan, kernel, mode=
'valid')
362 low, high = np.quantile(smoothedOverscan, [self.config.bandingFractionLow,
363 self.config.bandingFractionHigh])
364 outputStats[
'AMP_BANDING'].append(float(high - low))
366 if self.config.bandingUseHalfDetector:
367 fullLength = len(outputStats[
'AMP_BANDING'])
368 outputStats[
'DET_BANDING'] = float(np.nanmedian(outputStats[
'AMP_BANDING'][0:fullLength//2]))
370 outputStats[
'DET_BANDING'] = float(np.nanmedian(outputStats[
'AMP_BANDING']))
375 """Task to measure metrics from image slicing.
381 overscans : `list` [`lsst.pipe.base.Struct`]
382 List of overscan results. Expected fields are:
385 Value or fit subtracted
from the amplifier image data
388 Value
or fit subtracted
from the overscan image data
391 Image of the overscan region
with the overscan
393 quantity
is used to estimate the amplifier read noise
398 outputStats : `dict` [`str`, [`dict` [`str`,`float]]
399 Dictionary of measurements, keyed by amplifier name
and
404 detector = inputExp.getDetector()
405 kernel = self.makeKernel(self.config.projectionKernelSize)
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'
417 for amp
in detector.getAmplifiers():
418 ampArray = inputExp.image[amp.getBBox()].array
420 horizontalProjection = np.mean(ampArray, axis=0)
421 verticalProjection = np.mean(ampArray, axis=1)
423 horizontalProjection = np.convolve(horizontalProjection, kernel, mode=convolveMode)
424 verticalProjection = np.convolve(verticalProjection, kernel, mode=convolveMode)
426 outputStats[
'AMP_HPROJECTION'][amp.getName()] = horizontalProjection.tolist()
427 outputStats[
'AMP_VPROJECTION'][amp.getName()] = verticalProjection.tolist()
429 if self.config.doProjectionFft:
430 horizontalWindow = np.ones_like(horizontalProjection)
431 verticalWindow = np.ones_like(verticalProjection)
432 if self.config.projectionFftWindow ==
"NONE":
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))
444 raise RuntimeError(f
"Invalid window function: {self.config.projectionFftWindow}")
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()
A class to contain the data, WCS, and other information needed to describe an image of the sky.
A class to represent a 2-dimensional array of pixels.
Pass parameters to a Statistics object.
measureBanding(self, inputExp, overscans)
measureCti(self, inputExp, overscans, gains)
measureProjectionStatistics(self, inputExp, overscans)
__init__(self, statControl=None, **kwargs)
daf::base::PropertyList * list
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)
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)