Coverage for python / lsst / ip / isr / overscan.py: 14%
426 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-23 08:35 +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__ = ["OverscanCorrectionTaskConfig", "OverscanCorrectionTask",
23 "SerialOverscanCorrectionTaskConfig", "SerialOverscanCorrectionTask",
24 "ParallelOverscanCorrectionTaskConfig", "ParallelOverscanCorrectionTask",
25 ]
27import numpy as np
28from scipy.signal import medfilt
30import lsst.afw.math as afwMath
31import lsst.afw.image as afwImage
32import lsst.geom as geom
33import lsst.pipe.base as pipeBase
34import lsst.pex.config as pexConfig
36from .isr import fitOverscanImage, fitOverscanImageMean
37from .isrFunctions import makeThresholdMask
40class OverscanCorrectionTaskConfigBase(pexConfig.Config):
41 """Overscan correction options.
42 """
43 fitType = pexConfig.ChoiceField(
44 dtype=str,
45 doc="The method for fitting the overscan bias level.",
46 default='MEDIAN',
47 allowed={
48 "POLY": "Fit ordinary polynomial to the longest axis of the overscan region",
49 "CHEB": "Fit Chebyshev polynomial to the longest axis of the overscan region",
50 "LEG": "Fit Legendre polynomial to the longest axis of the overscan region",
51 "NATURAL_SPLINE": "Fit natural spline to the longest axis of the overscan region",
52 "CUBIC_SPLINE": "Fit cubic spline to the longest axis of the overscan region",
53 "AKIMA_SPLINE": "Fit Akima spline to the longest axis of the overscan region",
54 "MEAN": "Correct using the mean of the overscan region",
55 "MEANCLIP": "Correct using a clipped mean of the overscan region",
56 "MEDIAN": "Correct using the median of the overscan region",
57 "MEDIAN_PER_ROW": "Correct using the median per row of the overscan region",
58 "MEAN_PER_ROW": "Correct using the mean per row of the overscan region",
59 },
60 )
61 order = pexConfig.Field(
62 dtype=int,
63 doc=("Order of polynomial to fit if overscan fit type is a polynomial, "
64 "or number of spline knots if overscan fit type is a spline."),
65 default=1,
66 )
67 numSigmaClip = pexConfig.Field(
68 dtype=float,
69 doc="Rejection threshold (sigma) for collapsing overscan before fit",
70 default=3.0,
71 )
72 maskPlanes = pexConfig.ListField(
73 dtype=str,
74 doc="Mask planes to reject when measuring overscan",
75 default=['BAD', 'SAT'],
76 )
77 overscanIsInt = pexConfig.Field(
78 dtype=bool,
79 doc="Treat overscan as an integer image for purposes of fitType=MEDIAN"
80 " and fitType=MEDIAN_PER_ROW.",
81 default=True,
82 )
83 maxDeviation = pexConfig.Field(
84 dtype=float,
85 doc="Maximum deviation from median (in ADU) to mask in overscan correction; "
86 "Will be applied to the absolute deviation if doAbsoluteMaxDeviation=True.",
87 default=1000.0, check=lambda x: x > 0,
88 )
89 doAbsoluteMaxDeviation = pexConfig.Field(
90 dtype=bool,
91 doc="Apply the maxDeviation to the absolute value of the deviation? If "
92 "False, this will be a one-sided cut for positive-only deviations "
93 "(typically for parallel overscan subtraction.",
94 default=True,
95 )
98class OverscanCorrectionTaskConfig(OverscanCorrectionTaskConfigBase):
99 doParallelOverscan = pexConfig.Field(
100 dtype=bool,
101 doc="Correct using parallel overscan after serial overscan correction?",
102 default=False,
103 )
104 parallelOverscanMaskThreshold = pexConfig.Field(
105 dtype=int,
106 doc="Threshold above which pixels in the parallel overscan are masked as bleeds.",
107 default=100000,
108 )
109 parallelOverscanMaskGrowSize = pexConfig.Field(
110 dtype=int,
111 doc="Grow the SAT mask in the parallel overscan region by this many pixels. "
112 "This value was determined from the ITL chip in the LATISS camera.",
113 default=7,
114 )
115 leadingColumnsToSkip = pexConfig.Field(
116 dtype=int,
117 doc="Number of leading columns to skip in serial overscan correction.",
118 default=0,
119 )
120 trailingColumnsToSkip = pexConfig.Field(
121 dtype=int,
122 doc="Number of trailing columns to skip in serial overscan correction.",
123 default=0,
124 )
125 leadingRowsToSkip = pexConfig.Field(
126 dtype=int,
127 doc="Number of leading rows to skip in parallel overscan correction.",
128 default=0,
129 )
130 trailingRowsToSkip = pexConfig.Field(
131 dtype=int,
132 doc="Number of trailing rows to skip in parallel overscan correction.",
133 default=0,
134 )
137class OverscanCorrectionTaskBase(pipeBase.Task):
138 """Base Correction task for overscan.
140 This class contains a number of utilities that are easier to
141 understand and use when they are not embedded in nested if/else
142 loops.
144 Parameters
145 ----------
146 statControl : `lsst.afw.math.StatisticsControl`, optional
147 Statistics control object.
148 """
149 ConfigClass = OverscanCorrectionTaskConfigBase
150 _DefaultName = "overscanBase"
152 def __init__(self, statControl=None, **kwargs):
153 super().__init__(**kwargs)
154 self.allowDebug = True
156 if statControl:
157 self.statControl = statControl
158 else:
159 self.statControl = afwMath.StatisticsControl()
160 self.statControl.setNumSigmaClip(self.config.numSigmaClip)
161 self.statControl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.maskPlanes))
163 def run(self, exposure, amp, isTransposed=False):
164 raise NotImplementedError("run method is not defined for OverscanCorrectionTaskBase")
166 @staticmethod
167 def _maskRowsOrColumns(
168 exposure,
169 overscanBBox,
170 overscanMaskedImage,
171 overscanMask,
172 maxDeviation,
173 maskedRowColumnGrowSize,
174 medianSmoothingKernel,
175 medianSmoothingOutlierThreshold,
176 doAbsoluteMaxDeviation,
177 isTransposed,
178 ):
179 """Mask overscan rows (~serial) or columns (~parallel).
181 Parameters
182 ----------
183 exposure : `lsst.afw.image.Exposure`
184 Exposure containing the data.
185 overscanBBox: `lsst.geom.Box2I`
186 Bounding box for the overscan data.
187 overscanMaskedImage : `lsst.afw.image.MaskedImage`
188 Masked image containing the overscan data.
189 overscanMask : `np.ndarray`
190 Numpy array of the mask bits, anded with appropriate
191 mask planes.
192 maxDeviation : `float`
193 Maximum deviation from median (overscan units) to mask in overscan
194 correction. For parallel overscan this is a one-sided (positive
195 only) cut.
196 maskedRowColumnGrowSize : `int`
197 If a column (parallel overscan) or row (serial overscan) is
198 completely masked, then grow the mask by this radius. If the
199 value is <=0 then this will not be checked.
200 medianSmoothingKernel : `int`
201 Kernel (pixels) to smooth rows/columns. If <=0, median smoothing
202 is skipped. Otherwise must be odd.
203 medianSmoothingOutlierThreshold : `float`
204 Outlier threshold after median smoothing (overscan units). This
205 is applied only to positive outliers.
206 doAbsoluteMaxDeviation : `bool`
207 If true, deviation comparisons will use the absolute value;
208 otherwise it will cut positive outliers only.
209 isTransposed : `bool`
210 If true, then the data will be transposed before fitting
211 the overscan.
213 Returns
214 -------
215 badRowsOrColumns : `np.ndarray`
216 Array of bad rows (serial) or columns (parallel) that were
217 found, prior to dilation by maskedRowColumnGrowSize.
218 """
219 overscanArray = overscanMaskedImage.image.array
221 badRowsOrColumns = np.zeros(0, dtype=np.int64)
223 median = np.ma.median(np.ma.masked_where(overscanMask, overscanArray))
224 if doAbsoluteMaxDeviation:
225 delta = np.abs(overscanArray - median)
226 else:
227 delta = overscanArray - median
229 bad = np.where((delta > maxDeviation) & (~overscanMask))
230 # Mark the bad pixels as BAD
231 overscanMaskedImage.mask.array[bad] |= overscanMaskedImage.mask.getPlaneBitMask("BAD")
233 if isTransposed:
234 axis = 0
235 nComp = overscanArray.shape[0]
236 else:
237 axis = 1
238 nComp = overscanArray.shape[1]
240 # Check for completely masked row/column (from maxDeviation or
241 # previously applied SAT flag.)
242 if len(bad) > 0:
243 # We only need to look at the bad pixels set here for this
244 # mask growth.
245 overscanMaskTemp = np.zeros_like(overscanMask)
246 overscanMaskTemp[bad] = True
248 nMaskedArray = np.sum(overscanMaskTemp, axis=axis, dtype=np.int32)
249 badRowsOrColumns, = np.where(nMaskedArray == nComp)
251 # Perform median-smoothing outlier rejection if desired.
252 if medianSmoothingKernel > 0:
253 # We do a straight numpy median ignoring the mask.
254 # This will be fine because it avoids missing values,
255 # and very large deviations have already been flagged by
256 # maxDeviation or SAT.
257 rowsCols = np.median(overscanArray, axis=axis)
258 filtered = medfilt(rowsCols, kernel_size=medianSmoothingKernel)
259 delta = rowsCols - filtered
261 # We cannot reliably look for outliers within a kernel length
262 # of the edges.
263 high, = np.where(delta[medianSmoothingKernel: -medianSmoothingKernel]
264 >= medianSmoothingOutlierThreshold)
265 high += medianSmoothingKernel
267 if len(high) > 0:
268 badRowsOrColumns = np.unique(np.append(badRowsOrColumns, high))
270 # If we have any bad rows/columns, we need to dilate them
271 # and apply the mask to the parent overscan image.
272 if len(badRowsOrColumns) > 0:
273 dataView = afwImage.MaskedImageF(exposure.maskedImage,
274 overscanBBox,
275 afwImage.PARENT)
276 if isTransposed:
277 pixelsCopy = dataView.image.array[:, badRowsOrColumns].copy()
278 dataView.image.array[:, badRowsOrColumns] = 1e30
279 else:
280 pixelsCopy = dataView.image.array[badRowsOrColumns, :].copy()
281 dataView.image.array[badRowsOrColumns, :] = 1e30
283 makeThresholdMask(
284 maskedImage=dataView,
285 threshold=1e30,
286 growFootprints=maskedRowColumnGrowSize,
287 maskName="BAD",
288 )
290 if isTransposed:
291 dataView.image.array[:, badRowsOrColumns] = pixelsCopy
292 else:
293 dataView.image.array[badRowsOrColumns, :] = pixelsCopy
295 return badRowsOrColumns
297 def correctOverscan(
298 self,
299 exposure,
300 amp,
301 imageBBox,
302 overscanBBox,
303 isTransposed=True,
304 leadingToSkip=0,
305 trailingToSkip=0,
306 overscanFraction=1.0,
307 imageThreshold=np.inf,
308 maskedRowColumnGrowSize=0,
309 medianSmoothingKernel=0,
310 medianSmoothingOutlierThreshold=np.inf,
311 ):
312 """Trim the exposure, fit the overscan, subtract the fit, and
313 calculate statistics.
315 Parameters
316 ----------
317 exposure : `lsst.afw.image.Exposure`
318 Exposure containing the data.
319 amp : `lsst.afw.cameraGeom.Amplifier`
320 The amplifier that is to be corrected.
321 imageBBox: `lsst.geom.Box2I`
322 Bounding box of the image data that will have the overscan
323 subtracted. If parallel overscan will be performed, that
324 area is added to the image bounding box during serial
325 overscan correction.
326 overscanBBox: `lsst.geom.Box2I`
327 Bounding box for the overscan data.
328 isTransposed: `bool`
329 If true, then the data will be transposed before fitting
330 the overscan.
331 leadingToSkip : `int`, optional
332 Leading rows/columns to skip.
333 trailingToSkip : `int`, optional
334 Leading rows/columns to skip.
335 overscanFraction : `float`, optional
336 If the overscan region median is greater than overscanFraction
337 and the imaging region median is greater than imageThreshold
338 then overscan correction will be skipped.
339 maxLevel : `float`, optional
340 If the overscan region median is greater than overscanFraction
341 and the imaging region median is greater than imageThreshold
342 then overscan correction will be skipped.
343 maskedRowColumnGrowSize : `int`, optional
344 If a column (parallel overscan) or row (serial overscan) is
345 completely masked, then grow the mask by this radius. If the
346 value is <=0 then this will not be checked.
347 medianSmoothingKernel : `int`, optional
348 Kernel (pixels) to use to smooth rows/columns for row/column
349 outlier rejection. Must be odd if positive; if <=0 median
350 smoothing will not be used to find outliers.
351 medianSmoothingOutlierThreshold : `float`, optional
352 Threshold to look for outliers after median smoothing (adu).
354 Returns
355 -------
356 results : `lsst.pipe.base.Struct`
357 ``ampOverscanModel``
358 Overscan model broadcast to the full image size.
359 (`lsst.afw.image.Exposure`)
360 ``overscanOverscanModel``
361 Overscan model broadcast to the full overscan image
362 size. (`lsst.afw.image.Exposure`)
363 ``overscanImage``
364 Overscan image with the overscan fit subtracted.
365 (`lsst.afw.image.Exposure`)
366 ``overscanValue``
367 Overscan model. (`float` or `np.array`)
368 ``overscanMean``
369 Mean value of the overscan fit. (`float`)
370 ``overscanMedian``
371 Median value of the overscan fit. (`float`)
372 ``overscanSigma``
373 Standard deviation of the overscan fit. (`float`)
374 ``overscanMeanResidual``
375 Mean value of the overscan region after overscan
376 subtraction. (`float`)
377 ``overscanMedianResidual``
378 Median value of the overscan region after overscan
379 subtraction. (`float`)
380 ``overscanSigmaResidual``
381 Standard deviation of the overscan region after
382 overscan subtraction. (`float`)
383 """
384 overscanBox = self.trimOverscan(exposure, amp, overscanBBox,
385 leadingToSkip,
386 trailingToSkip,
387 transpose=isTransposed)
388 overscanImage = exposure[overscanBox].getMaskedImage()
390 # Record the gain value if necessary to convert configs from
391 # electron to adu.
392 if exposure.metadata.get("LSST ISR UNITS", "adu") == "electron":
393 gain = exposure.metadata[f"LSST ISR GAIN {amp.getName()}"]
394 else:
395 gain = 1.0
397 # Mask pixels.
398 maskVal = overscanImage.mask.getPlaneBitMask(self.config.maskPlanes)
399 overscanMask = ~((overscanImage.mask.array & maskVal) == 0)
401 badResults = False
402 overscanMedian = np.nanmedian(overscanImage.image.array)
403 imageMedian = np.nanmedian(exposure[imageBBox].image.array)
405 if np.all(overscanMask):
406 self.log.warning(
407 "All overscan pixels masked when attempting overscan correction for %s",
408 amp.getName(),
409 )
410 badResults = True
411 elif overscanMedian/imageMedian > overscanFraction and imageMedian > imageThreshold:
412 self.log.warning(
413 "The level in the overscan region (%.2f) compared to the image region (%.2f) is "
414 "greater than the maximum fraction (%.2f) for %s",
415 overscanMedian,
416 imageMedian,
417 overscanFraction,
418 amp.getName(),
419 )
420 badResults = True
422 if badResults:
423 # Do not do overscan subtraction at all.
424 badRowsOrColumns = np.zeros(0, dtype=np.int64)
425 overscanResults = pipeBase.Struct(
426 overscanValue=0.0,
427 overscanMean=0.0,
428 overscanMedian=0.0,
429 overscanSigma=0.0,
430 )
431 else:
432 badRowsOrColumns = self._maskRowsOrColumns(
433 exposure,
434 overscanBBox,
435 overscanImage,
436 overscanMask,
437 gain * self.config.maxDeviation,
438 maskedRowColumnGrowSize,
439 medianSmoothingKernel,
440 gain * medianSmoothingOutlierThreshold,
441 self.config.doAbsoluteMaxDeviation,
442 isTransposed,
443 )
444 # Do overscan fit.
445 # CZW: Handle transposed correctly.
446 overscanResults = self.fitOverscan(overscanImage, isTransposed=isTransposed)
448 # Correct image region (and possibly parallel-overscan region).
449 ampImage = exposure[imageBBox]
450 ampOverscanModel = self.broadcastFitToImage(overscanResults.overscanValue,
451 ampImage.image.array,
452 transpose=isTransposed)
453 ampImage.image.array -= ampOverscanModel
455 # Correct overscan region (and possibly doubly-overscaned
456 # region).
457 overscanImage = exposure[overscanBBox]
458 # CZW: Transposed?
459 overscanOverscanModel = self.broadcastFitToImage(overscanResults.overscanValue,
460 overscanImage.image.array)
461 self.debugView(overscanImage, overscanResults.overscanValue, amp, isTransposed=isTransposed)
462 overscanImage.image.array -= overscanOverscanModel
464 # Find residual fit statistics.
465 stats = afwMath.makeStatistics(overscanImage.getMaskedImage(),
466 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP, self.statControl)
467 residualMean = stats.getValue(afwMath.MEAN)
468 residualMedian = stats.getValue(afwMath.MEDIAN)
469 residualSigma = stats.getValue(afwMath.STDEVCLIP)
471 return pipeBase.Struct(
472 ampOverscanModel=ampOverscanModel,
473 overscanOverscanModel=overscanOverscanModel,
474 overscanImage=overscanImage,
475 overscanValue=overscanResults.overscanValue,
476 overscanMean=overscanResults.overscanMean,
477 overscanMedian=overscanResults.overscanMedian,
478 overscanSigma=overscanResults.overscanSigma,
479 overscanMeanResidual=residualMean,
480 overscanMedianResidual=residualMedian,
481 overscanSigmaResidual=residualSigma,
482 overscanBadRowsOrColumns=badRowsOrColumns,
483 )
485 def broadcastFitToImage(self, overscanValue, imageArray, transpose=False):
486 """Broadcast 0 or 1 dimension fit to appropriate shape.
488 Parameters
489 ----------
490 overscanValue : `numpy.ndarray`, (Nrows, ) or scalar
491 Overscan fit to broadcast.
492 imageArray : `numpy.ndarray`, (Nrows, Ncols)
493 Image array that we want to match.
494 transpose : `bool`, optional
495 Switch order to broadcast along the other axis.
497 Returns
498 -------
499 overscanModel : `numpy.ndarray`, (Nrows, Ncols) or scalar
500 Expanded overscan fit.
502 Raises
503 ------
504 RuntimeError
505 Raised if no axis has the appropriate dimension.
506 """
507 if isinstance(overscanValue, np.ndarray):
508 overscanModel = np.zeros_like(imageArray)
510 if transpose is False:
511 if imageArray.shape[0] == overscanValue.shape[0]:
512 overscanModel[:, :] = overscanValue[:, np.newaxis]
513 elif imageArray.shape[1] == overscanValue.shape[0]:
514 overscanModel[:, :] = overscanValue[np.newaxis, :]
515 elif imageArray.shape[0] == overscanValue.shape[1]:
516 overscanModel[:, :] = overscanValue[np.newaxis, :]
517 else:
518 raise RuntimeError(f"Could not broadcast {overscanValue.shape} to "
519 f"match {imageArray.shape}")
520 else:
521 if imageArray.shape[1] == overscanValue.shape[0]:
522 overscanModel[:, :] = overscanValue[np.newaxis, :]
523 elif imageArray.shape[0] == overscanValue.shape[0]:
524 overscanModel[:, :] = overscanValue[:, np.newaxis]
525 elif imageArray.shape[1] == overscanValue.shape[1]:
526 overscanModel[:, :] = overscanValue[:, np.newaxis]
527 else:
528 raise RuntimeError(f"Could not broadcast {overscanValue.shape} to "
529 f"match {imageArray.shape}")
530 else:
531 overscanModel = overscanValue
533 return overscanModel
535 def trimOverscan(self, exposure, amp, bbox, skipLeading, skipTrailing, transpose=False):
536 """Trim overscan region to remove edges.
538 Parameters
539 ----------
540 exposure : `lsst.afw.image.Exposure`
541 Exposure containing data.
542 amp : `lsst.afw.cameraGeom.Amplifier`
543 Amplifier containing geometry information.
544 bbox : `lsst.geom.Box2I`
545 Bounding box of the overscan region.
546 skipLeading : `int`
547 Number of leading (towards data region) rows/columns to skip.
548 skipTrailing : `int`
549 Number of trailing (away from data region) rows/columns to skip.
550 transpose : `bool`, optional
551 Operate on the transposed array.
553 Returns
554 -------
555 overscanArray : `numpy.array`, (N, M)
556 Data array to fit.
557 overscanMask : `numpy.array`, (N, M)
558 Data mask.
559 """
560 dx0, dy0, dx1, dy1 = (0, 0, 0, 0)
561 dataBBox = amp.getRawDataBBox()
562 if transpose:
563 if dataBBox.getBeginY() < bbox.getBeginY():
564 dy0 += skipLeading
565 dy1 -= skipTrailing
566 else:
567 dy0 += skipTrailing
568 dy1 -= skipLeading
569 else:
570 if dataBBox.getBeginX() < bbox.getBeginX():
571 dx0 += skipLeading
572 dx1 -= skipTrailing
573 else:
574 dx0 += skipTrailing
575 dx1 -= skipLeading
577 overscanBBox = geom.Box2I(bbox.getBegin() + geom.Extent2I(dx0, dy0),
578 geom.Extent2I(bbox.getWidth() - dx0 + dx1,
579 bbox.getHeight() - dy0 + dy1))
580 return overscanBBox
582 def fitOverscan(self, overscanImage, isTransposed=False):
583 if self.config.fitType in ('MEAN', 'MEANCLIP', 'MEDIAN'):
584 # Transposition has no effect here.
585 overscanResult = self.measureConstantOverscan(overscanImage)
586 overscanValue = overscanResult.overscanValue
587 overscanMean = overscanValue
588 overscanMedian = overscanValue
589 overscanSigma = 0.0
590 elif self.config.fitType in ('MEDIAN_PER_ROW', 'MEAN_PER_ROW', 'POLY', 'CHEB', 'LEG',
591 'NATURAL_SPLINE', 'CUBIC_SPLINE', 'AKIMA_SPLINE'):
592 # Force transposes as needed
593 overscanResult = self.measureVectorOverscan(overscanImage, isTransposed)
594 overscanValue = overscanResult.overscanValue
596 stats = afwMath.makeStatistics(overscanResult.overscanValue,
597 afwMath.MEAN | afwMath.MEDIAN | afwMath.STDEVCLIP,
598 self.statControl)
599 overscanMean = stats.getValue(afwMath.MEAN)
600 overscanMedian = stats.getValue(afwMath.MEDIAN)
601 overscanSigma = stats.getValue(afwMath.STDEVCLIP)
602 else:
603 raise ValueError('%s : %s an invalid overscan type' %
604 ("overscanCorrection", self.config.fitType))
606 return pipeBase.Struct(overscanValue=overscanValue,
607 overscanMean=overscanMean,
608 overscanMedian=overscanMedian,
609 overscanSigma=overscanSigma,
610 )
612 def maskParallelOverscan(self, exposure, detector):
613 """Mask the union of high values on all amplifiers in the parallel
614 overscan.
616 This operates on the image in-place.
618 Parameters
619 ----------
620 exposure : `lsst.afw.image.Exposure`
621 An untrimmed raw exposure.
622 detector : `lsst.afw.cameraGeom.Detector`
623 The detetor to use for amplifier geometry.
624 """
625 parallelMask = None
627 for amp in detector:
628 dataView = afwImage.MaskedImageF(exposure.getMaskedImage(),
629 amp.getRawParallelOverscanBBox(),
630 afwImage.PARENT)
631 # This should mark all the saturated pixels as SAT.
632 makeThresholdMask(
633 maskedImage=dataView,
634 threshold=self.config.parallelOverscanMaskThreshold,
635 growFootprints=self.config.parallelOverscanMaskGrowSize,
636 maskName="SAT"
637 )
638 if parallelMask is None:
639 parallelMask = dataView.mask.array
640 else:
641 parallelMask |= dataView.mask.array
642 for amp in detector:
643 dataView = afwImage.MaskedImageF(exposure.getMaskedImage(),
644 amp.getRawParallelOverscanBBox(),
645 afwImage.PARENT)
646 dataView.mask.array |= parallelMask
648 # Constant methods
649 def measureConstantOverscan(self, image):
650 """Measure a constant overscan value.
652 Parameters
653 ----------
654 image : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage`
655 Image data to measure the overscan from.
657 Returns
658 -------
659 results : `lsst.pipe.base.Struct`
660 Overscan result with entries:
661 - ``overscanValue``: Overscan value to subtract (`float`)
662 - ``isTransposed``: Orientation of the overscan (`bool`)
663 """
664 fitType = afwMath.stringToStatisticsProperty(self.config.fitType)
665 overscanValue = afwMath.makeStatistics(image, fitType, self.statControl).getValue()
667 return pipeBase.Struct(overscanValue=overscanValue,
668 isTransposed=False)
670 # Vector correction utilities
671 def getImageArray(self, image):
672 """Extract the numpy array from the input image.
674 Parameters
675 ----------
676 image : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage`
677 Image data to pull array from.
679 calcImage : `numpy.ndarray`
680 Image data array for numpy operating.
681 """
682 if hasattr(image, "getImage"):
683 calcImage = image.getImage().getArray()
684 calcImage = np.ma.masked_where(image.getMask().getArray() & self.statControl.getAndMask(),
685 calcImage)
686 else:
687 calcImage = image.getArray()
688 return calcImage
690 def maskOutliers(self, imageArray):
691 """Mask outliers in a row of overscan data from a robust sigma
692 clipping procedure.
694 Parameters
695 ----------
696 imageArray : `numpy.ndarray`
697 Image to filter along numpy axis=1.
699 Returns
700 -------
701 maskedArray : `numpy.ma.masked_array`
702 Masked image marking outliers.
703 """
704 lq, median, uq = np.percentile(np.ma.getdata(imageArray),
705 [25.0, 50.0, 75.0], axis=1)
706 axisMedians = median
707 axisStdev = 0.74*(uq - lq) # robust stdev
709 # Replace pixels that have excessively large stdev values
710 # with the median of stdev values. A large stdev likely
711 # indicates a bleed is spilling into the overscan.
712 axisStdev = np.where(axisStdev > 2.0 * np.median(axisStdev),
713 np.median(axisStdev), axisStdev)
715 # Mask pixels that are N-sigma away from their array medians.
716 diff = np.abs(imageArray - axisMedians[:, np.newaxis])
717 masked = np.ma.masked_where(diff > self.statControl.getNumSigmaClip()
718 * axisStdev[:, np.newaxis], imageArray)
720 return masked
722 def fillMaskedPixels(self, overscanVector):
723 """Fill masked/NaN pixels in the overscan.
725 Parameters
726 ----------
727 overscanVector : `np.array` or `np.ma.masked_array`
728 Overscan vector to fill.
730 Returns
731 -------
732 overscanVector : `np.ma.masked_array`
733 Filled vector.
735 Notes
736 -----
737 Each maskSlice is a section of overscan with contiguous masks.
738 Ideally this adds 5 pixels from the left and right of that
739 mask slice, and takes the median of those values to fill the
740 slice. If this isn't possible, the median of all non-masked
741 values is used. The mask is removed for the pixels filled.
742 """
743 workingCopy = overscanVector
744 if not isinstance(overscanVector, np.ma.MaskedArray):
745 workingCopy = np.ma.masked_array(overscanVector,
746 mask=~np.isfinite(overscanVector))
748 defaultValue = np.median(workingCopy.data[~workingCopy.mask])
749 for maskSlice in np.ma.clump_masked(workingCopy):
750 neighborhood = []
751 if maskSlice.start > 5:
752 neighborhood.extend(workingCopy[maskSlice.start - 5:maskSlice.start].data)
753 if maskSlice.stop < workingCopy.size - 5:
754 neighborhood.extend(workingCopy[maskSlice.stop:maskSlice.stop+5].data)
755 if len(neighborhood) > 0:
756 workingCopy.data[maskSlice] = np.nanmedian(neighborhood)
757 workingCopy.mask[maskSlice] = False
758 else:
759 workingCopy.data[maskSlice] = defaultValue
760 workingCopy.mask[maskSlice] = False
761 return workingCopy
763 def collapseArray(self, maskedArray, fillMasked=True):
764 """Collapse overscan array (and mask) to a 1-D vector of values.
766 Parameters
767 ----------
768 maskedArray : `numpy.ma.masked_array`
769 Masked array of input overscan data.
770 fillMasked : `bool`, optional
771 If true, fill any pixels that are masked with a median of
772 neighbors.
774 Returns
775 -------
776 collapsed : `numpy.ma.masked_array`
777 Single dimensional overscan data, combined with the mean.
779 """
780 collapsed = np.mean(maskedArray, axis=1)
781 if collapsed.mask.sum() > 0 and fillMasked:
782 collapsed = self.fillMaskedPixels(collapsed)
784 return collapsed
786 def splineFit(self, indices, collapsed, numBins):
787 """Wrapper function to match spline fit API to polynomial fit API.
789 Parameters
790 ----------
791 indices : `numpy.ndarray`
792 Locations to evaluate the spline.
793 collapsed : `numpy.ndarray`
794 Collapsed overscan values corresponding to the spline
795 evaluation points.
796 numBins : `int`
797 Number of bins to use in constructing the spline.
799 Returns
800 -------
801 interp : `lsst.afw.math.Interpolate`
802 Interpolation object for later evaluation.
803 """
804 if not np.ma.is_masked(collapsed):
805 collapsed.mask = np.array(len(collapsed)*[np.ma.nomask])
807 numPerBin, binEdges = np.histogram(indices, bins=numBins,
808 weights=1 - collapsed.mask.astype(int))
809 with np.errstate(invalid="ignore"):
810 values = np.histogram(indices, bins=numBins,
811 weights=collapsed.data*~collapsed.mask)[0]/numPerBin
812 binCenters = np.histogram(indices, bins=numBins,
813 weights=indices*~collapsed.mask)[0]/numPerBin
815 if len(binCenters[numPerBin > 0]) < 5:
816 self.log.warning("Cannot do spline fitting for overscan: %s valid points.",
817 len(binCenters[numPerBin > 0]))
818 # Return a scalar value if we have one, otherwise
819 # return zero. This amplifier is hopefully already
820 # masked.
821 if len(values[numPerBin > 0]) != 0:
822 return float(values[numPerBin > 0][0])
823 else:
824 return 0.0
826 interp = afwMath.makeInterpolate(binCenters.astype(float)[numPerBin > 0],
827 values.astype(float)[numPerBin > 0],
828 afwMath.stringToInterpStyle(self.config.fitType))
829 return interp
831 @staticmethod
832 def splineEval(indices, interp):
833 """Wrapper function to match spline evaluation API to polynomial fit
834 API.
836 Parameters
837 ----------
838 indices : `numpy.ndarray`
839 Locations to evaluate the spline.
840 interp : `lsst.afw.math.interpolate`
841 Interpolation object to use.
843 Returns
844 -------
845 values : `numpy.ndarray`
846 Evaluated spline values at each index.
847 """
849 return interp.interpolate(indices.astype(float))
851 @staticmethod
852 def maskExtrapolated(collapsed):
853 """Create mask if edges are extrapolated.
855 Parameters
856 ----------
857 collapsed : `numpy.ma.masked_array`
858 Masked array to check the edges of.
860 Returns
861 -------
862 maskArray : `numpy.ndarray`
863 Boolean numpy array of pixels to mask.
864 """
865 maskArray = np.full_like(collapsed, False, dtype=bool)
866 if np.ma.is_masked(collapsed):
867 num = len(collapsed)
868 for low in range(num):
869 if not collapsed.mask[low]:
870 break
871 if low > 0:
872 maskArray[:low] = True
873 for high in range(1, num):
874 if not collapsed.mask[-high]:
875 break
876 if high > 1:
877 maskArray[-high:] = True
878 return maskArray
880 def measureVectorOverscan(self, image, isTransposed=False):
881 """Calculate the 1-d vector overscan from the input overscan image.
883 Parameters
884 ----------
885 image : `lsst.afw.image.MaskedImage`
886 Image containing the overscan data.
887 isTransposed : `bool`
888 If true, the image has been transposed.
890 Returns
891 -------
892 results : `lsst.pipe.base.Struct`
893 Overscan result with entries:
895 ``overscanValue``
896 Overscan value to subtract (`float`)
897 ``maskArray``
898 List of rows that should be masked as ``SUSPECT`` when the
899 overscan solution is applied. (`list` [ `bool` ])
900 ``isTransposed``
901 Indicates if the overscan data was transposed during
902 calcuation, noting along which axis the overscan should be
903 subtracted. (`bool`)
904 """
905 calcImage = self.getImageArray(image)
907 # operate on numpy-arrays from here
908 if isTransposed:
909 calcImage = np.transpose(calcImage)
910 masked = self.maskOutliers(calcImage)
912 if self.config.fitType in ('MEDIAN_PER_ROW', "MEAN_PER_ROW"):
913 if self.config.overscanIsInt:
914 mi = afwImage.MaskedImageI(image.getBBox())
915 masked = masked.astype(int)
916 else:
917 mi = image.clone()
919 if isTransposed:
920 masked = masked.transpose()
922 mi.image.array[:, :] = masked.data[:, :]
923 if bool(masked.mask.shape):
924 mi.mask.array[:, :] = masked.mask[:, :]
926 if self.config.fitType == "MEDIAN_PER_ROW":
927 overscanVector = fitOverscanImage(mi, self.config.maskPlanes, isTransposed)
928 else:
929 overscanVector = fitOverscanImageMean(mi, self.config.maskPlanes, isTransposed)
931 overscanVector = self.fillMaskedPixels(overscanVector)
932 maskArray = self.maskExtrapolated(overscanVector)
933 else:
934 collapsed = self.collapseArray(masked)
936 num = len(collapsed)
937 indices = 2.0*np.arange(num)/float(num) - 1.0
939 poly = np.polynomial
940 fitter, evaler = {
941 'POLY': (poly.polynomial.polyfit, poly.polynomial.polyval),
942 'CHEB': (poly.chebyshev.chebfit, poly.chebyshev.chebval),
943 'LEG': (poly.legendre.legfit, poly.legendre.legval),
944 'NATURAL_SPLINE': (self.splineFit, self.splineEval),
945 'CUBIC_SPLINE': (self.splineFit, self.splineEval),
946 'AKIMA_SPLINE': (self.splineFit, self.splineEval)
947 }[self.config.fitType]
949 # These are the polynomial coefficients, or an
950 # interpolation object.
951 coeffs = fitter(indices, collapsed, self.config.order)
953 if isinstance(coeffs, float):
954 self.log.warning("Using fallback value %f due to fitter failure. Amplifier will be masked.",
955 coeffs)
956 overscanVector = np.full_like(indices, coeffs)
957 maskArray = np.full_like(collapsed, True, dtype=bool)
958 else:
959 # Otherwise we can just use things as normal.
960 overscanVector = evaler(indices, coeffs)
961 maskArray = self.maskExtrapolated(collapsed)
963 return pipeBase.Struct(overscanValue=np.array(overscanVector),
964 maskArray=maskArray,
965 isTransposed=isTransposed)
967 def debugView(self, image, model, amp=None, isTransposed=True):
968 """Debug display for the final overscan solution.
970 Parameters
971 ----------
972 image : `lsst.afw.image.Image`
973 Input image the overscan solution was determined from.
974 model : `numpy.ndarray` or `float`
975 Overscan model determined for the image.
976 amp : `lsst.afw.cameraGeom.Amplifier`, optional
977 Amplifier to extract diagnostic information.
978 isTransposed : `bool`, optional
979 Does the data need to be transposed before display?
980 """
981 import lsstDebug
982 if not lsstDebug.Info(__name__).display:
983 return
984 if not self.allowDebug:
985 return
987 calcImage = self.getImageArray(image)
988 # CZW: Check that this is ok
989 if isTransposed:
990 calcImage = np.transpose(calcImage)
991 masked = self.maskOutliers(calcImage)
992 collapsed = self.collapseArray(masked, fillMasked=False)
994 num = len(collapsed)
995 indices = 2.0 * np.arange(num)/float(num) - 1.0
996 indices = np.arange(num)
998 if np.ma.is_masked(collapsed):
999 collapsedMask = collapsed.mask
1000 else:
1001 collapsedMask = np.array(num*[np.ma.nomask])
1003 import matplotlib.pyplot as plot
1004 figure = plot.figure(1)
1005 figure.clear()
1006 axes = figure.add_axes((0.1, 0.1, 0.8, 0.8))
1007 axes.plot(indices[~collapsedMask], collapsed[~collapsedMask], 'k+')
1008 if collapsedMask.sum() > 0:
1009 axes.plot(indices[collapsedMask], collapsed.data[collapsedMask], 'b+')
1010 if isinstance(model, np.ndarray):
1011 plotModel = model
1012 else:
1013 plotModel = np.zeros_like(indices)
1014 plotModel += model
1016 axes.plot(indices, plotModel, 'r-')
1017 plot.xlabel("position along overscan region")
1018 plot.ylabel("pixel value/fit value")
1019 if amp:
1020 plot.title(f"{amp.getName()} DataX: "
1021 f"[{amp.getRawDataBBox().getBeginX()}:{amp.getRawBBox().getEndX()}]"
1022 f"OscanX: [{amp.getRawHorizontalOverscanBBox().getBeginX()}:"
1023 f"{amp.getRawHorizontalOverscanBBox().getEndX()}] {self.config.fitType}")
1024 else:
1025 plot.title("No amp supplied.")
1026 figure.show()
1027 prompt = "Press Enter or c to continue [chp]..."
1028 while True:
1029 ans = input(prompt).lower()
1030 if ans in ("", " ", "c",):
1031 break
1032 elif ans in ("p", ):
1033 import pdb
1034 pdb.set_trace()
1035 elif ans in ('x', ):
1036 self.allowDebug = False
1037 break
1038 elif ans in ("h", ):
1039 print("[h]elp [c]ontinue [p]db e[x]itDebug")
1040 plot.close()
1043class OverscanCorrectionTask(OverscanCorrectionTaskBase):
1044 """Correction task for serial/parallel overscan.
1046 (Will be deprecated)
1048 This class contains a number of utilities that are easier to
1049 understand and use when they are not embedded in nested if/else
1050 loops.
1052 Parameters
1053 ----------
1054 statControl : `lsst.afw.math.StatisticsControl`, optional
1055 Statistics control object.
1056 """
1057 ConfigClass = OverscanCorrectionTaskConfig
1058 _DefaultName = "overscan"
1060 def run(self, exposure, amp, isTransposed=False):
1061 """Measure and remove serial/parallel overscan from an amplifier image.
1063 This will be deprecated.
1065 Parameters
1066 ----------
1067 exposure : `lsst.afw.image.Exposure`
1068 Image data that will have the overscan corrections applied.
1069 amp : `lsst.afw.cameraGeom.Amplifier`
1070 Amplifier to use for debugging purposes.
1071 isTransposed : `bool`, optional
1072 Is the image transposed, such that serial and parallel
1073 overscan regions are reversed? Default is False.
1075 Returns
1076 -------
1077 overscanResults : `lsst.pipe.base.Struct`
1078 Result struct with components:
1080 ``imageFit``
1081 Value or fit subtracted from the amplifier image data
1082 (scalar or `lsst.afw.image.Image`).
1083 ``overscanFit``
1084 Value or fit subtracted from the serial overscan image
1085 data (scalar or `lsst.afw.image.Image`).
1086 ``overscanImage``
1087 Image of the serial overscan region with the serial
1088 overscan correction applied
1089 (`lsst.afw.image.Image`). This quantity is used to
1090 estimate the amplifier read noise empirically.
1091 ``parallelOverscanFit``
1092 Value or fit subtracted from the parallel overscan
1093 image data (scalar, `lsst.afw.image.Image`, or None).
1094 ``parallelOverscanImage``
1095 Image of the parallel overscan region with the
1096 parallel overscan correction applied
1097 (`lsst.afw.image.Image` or None).
1098 ``overscanMean``
1099 Mean of the fit serial overscan region.
1100 This and the following values will be tuples of
1101 (serial, parallel) if doParallelOverscan=True.
1102 ``overscanMedian``
1103 Median of the fit serial overscan region.
1104 ``overscanSigma``
1105 Sigma of the fit serial overscan region.
1106 ``residualMean``
1107 Mean of the residual of the serial overscan region after
1108 correction.
1109 ``residualMedian``
1110 Median of the residual of the serial overscan region after
1111 correction.
1112 ``residualSigma``
1113 Mean of the residual of the serial overscan region after
1114 correction.
1117 Raises
1118 ------
1119 RuntimeError
1120 Raised if an invalid overscan type is set.
1121 """
1122 # Do Serial overscan first.
1123 serialOverscanBBox = amp.getRawSerialOverscanBBox()
1124 imageBBox = amp.getRawDataBBox()
1126 if self.config.doParallelOverscan:
1127 # We need to extend the serial overscan BBox to the full
1128 # size of the detector.
1129 parallelOverscanBBox = amp.getRawParallelOverscanBBox()
1130 imageBBox = imageBBox.expandedTo(parallelOverscanBBox)
1132 if isTransposed:
1133 serialOverscanBBox = geom.Box2I(
1134 geom.Point2I(serialOverscanBBox.getMinX(), imageBBox.getEndY()),
1135 geom.Extent2I(imageBBox.getWidth(), serialOverscanBBox.getHeight()),
1136 )
1137 else:
1138 serialOverscanBBox = geom.Box2I(
1139 geom.Point2I(serialOverscanBBox.getMinX(),
1140 imageBBox.getMinY()),
1141 geom.Extent2I(serialOverscanBBox.getWidth(),
1142 imageBBox.getHeight()),
1143 )
1145 serialResults = self.correctOverscan(
1146 exposure,
1147 amp,
1148 imageBBox,
1149 serialOverscanBBox,
1150 isTransposed=isTransposed,
1151 leadingToSkip=self.config.leadingColumnsToSkip,
1152 trailingToSkip=self.config.trailingColumnsToSkip,
1153 )
1154 overscanMean = serialResults.overscanMean
1155 overscanMedian = serialResults.overscanMedian
1156 overscanSigma = serialResults.overscanSigma
1157 residualMean = serialResults.overscanMeanResidual
1158 residualMedian = serialResults.overscanMedianResidual
1159 residualSigma = serialResults.overscanSigmaResidual
1161 # Do Parallel Overscan
1162 parallelResults = None
1163 if self.config.doParallelOverscan:
1164 # This does not need any extensions, as we'll only
1165 # subtract it from the data region.
1166 parallelOverscanBBox = amp.getRawParallelOverscanBBox()
1167 imageBBox = amp.getRawDataBBox()
1169 # The serial overscan correction has removed some signal
1170 # from the parallel overscan region, but that is largely a
1171 # constant offset. The collapseArray method now attempts
1172 # to fill fully masked columns with the median of
1173 # neighboring values, with a fallback to the median of the
1174 # correction in all other columns. Filling with neighbor
1175 # values ensures that large variations in the parallel
1176 # overscan do not create new outlier points. The
1177 # MEDIAN_PER_ROW method does this filling as a separate
1178 # operation, using the same method.
1179 parallelResults = self.correctOverscan(
1180 exposure,
1181 amp,
1182 imageBBox,
1183 parallelOverscanBBox,
1184 isTransposed=not isTransposed,
1185 leadingToSkip=self.config.leadingRowsToSkip,
1186 trailingToSkip=self.config.trailingRowsToSkip,
1187 )
1188 overscanMean = (overscanMean, parallelResults.overscanMean)
1189 overscanMedian = (overscanMedian, parallelResults.overscanMedian)
1190 overscanSigma = (overscanSigma, parallelResults.overscanSigma)
1191 residualMean = (residualMean, parallelResults.overscanMeanResidual)
1192 residualMedian = (residualMedian, parallelResults.overscanMedianResidual)
1193 residualSigma = (residualSigma, parallelResults.overscanSigmaResidual)
1195 parallelOverscanFit = parallelResults.overscanOverscanModel if parallelResults else None
1196 parallelOverscanImage = parallelResults.overscanImage if parallelResults else None
1198 return pipeBase.Struct(imageFit=serialResults.ampOverscanModel,
1199 overscanFit=serialResults.overscanOverscanModel,
1200 overscanImage=serialResults.overscanImage,
1202 parallelOverscanFit=parallelOverscanFit,
1203 parallelOverscanImage=parallelOverscanImage,
1204 overscanMean=overscanMean,
1205 overscanMedian=overscanMedian,
1206 overscanSigma=overscanSigma,
1207 residualMean=residualMean,
1208 residualMedian=residualMedian,
1209 residualSigma=residualSigma)
1212class SerialOverscanCorrectionTaskConfig(OverscanCorrectionTaskConfigBase):
1213 leadingToSkip = pexConfig.Field(
1214 dtype=int,
1215 doc="Number of leading values to skip in serial overscan correction.",
1216 default=0,
1217 )
1218 trailingToSkip = pexConfig.Field(
1219 dtype=int,
1220 doc="Number of trailing values to skip in serial overscan correction.",
1221 default=0,
1222 )
1225class SerialOverscanCorrectionTask(OverscanCorrectionTaskBase):
1226 """Correction task for serial overscan.
1228 Parameters
1229 ----------
1230 statControl : `lsst.afw.math.StatisticsControl`, optional
1231 Statistics control object.
1232 """
1233 ConfigClass = SerialOverscanCorrectionTaskConfig
1234 _DefaultName = "serialOverscan"
1236 def run(self, exposure, amp, isTransposed=False):
1237 """Measure and remove serial overscan from an amplifier image.
1239 Parameters
1240 ----------
1241 exposure : `lsst.afw.image.Exposure`
1242 Image data that will have the overscan corrections applied.
1243 amp : `lsst.afw.cameraGeom.Amplifier`
1244 Amplifier to use for debugging purposes.
1245 isTransposed : `bool`, optional
1246 Is the image transposed, such that serial and parallel
1247 overscan regions are reversed? Default is False.
1249 Returns
1250 -------
1251 overscanResults : `lsst.pipe.base.Struct`
1252 Result struct with components:
1254 ``imageFit``
1255 Value or fit subtracted from the amplifier image data
1256 (scalar or `lsst.afw.image.Image`).
1257 ``overscanFit``
1258 Value or fit subtracted from the serial overscan image
1259 data (scalar or `lsst.afw.image.Image`).
1260 ``overscanImage``
1261 Image of the serial overscan region with the serial
1262 overscan correction applied
1263 (`lsst.afw.image.Image`). This quantity is used to
1264 estimate the amplifier read noise empirically.
1265 ``overscanMean``
1266 Mean of the fit serial overscan region.
1267 ``overscanMedian``
1268 Median of the fit serial overscan region.
1269 ``overscanSigma``
1270 Sigma of the fit serial overscan region.
1271 ``residualMean``
1272 Mean of the residual of the serial overscan region after
1273 correction.
1274 ``residualMedian``
1275 Median of the residual of the serial overscan region after
1276 correction.
1277 ``residualSigma``
1278 Mean of the residual of the serial overscan region after
1279 correction.
1281 Raises
1282 ------
1283 RuntimeError
1284 Raised if an invalid overscan type is set.
1285 """
1286 serialOverscanBBox = amp.getRawSerialOverscanBBox()
1287 imageBBox = amp.getRawDataBBox()
1289 # We always want to extend the serial overscan bounding box to
1290 # the full size of the detector.
1291 parallelOverscanBBox = amp.getRawParallelOverscanBBox()
1292 imageBBox = imageBBox.expandedTo(parallelOverscanBBox)
1294 if isTransposed:
1295 serialOverscanBBox = geom.Box2I(
1296 geom.Point2I(serialOverscanBBox.getMinX(), imageBBox.getEndY()),
1297 geom.Extent2I(imageBBox.getWidth(), serialOverscanBBox.getHeight()),
1298 )
1299 else:
1300 serialOverscanBBox = geom.Box2I(
1301 geom.Point2I(serialOverscanBBox.getMinX(),
1302 imageBBox.getMinY()),
1303 geom.Extent2I(serialOverscanBBox.getWidth(),
1304 imageBBox.getHeight()),
1305 )
1307 results = self.correctOverscan(
1308 exposure,
1309 amp,
1310 imageBBox,
1311 serialOverscanBBox,
1312 isTransposed=isTransposed,
1313 leadingToSkip=self.config.leadingToSkip,
1314 trailingToSkip=self.config.trailingToSkip,
1315 )
1316 overscanMean = results.overscanMean
1317 overscanMedian = results.overscanMedian
1318 overscanSigma = results.overscanSigma
1319 residualMean = results.overscanMeanResidual
1320 residualMedian = results.overscanMedianResidual
1321 residualSigma = results.overscanSigmaResidual
1322 badRowsOrColumns = results.overscanBadRowsOrColumns
1324 return pipeBase.Struct(
1325 imageFit=results.ampOverscanModel,
1326 overscanFit=results.overscanOverscanModel,
1327 overscanImage=results.overscanImage,
1328 overscanMean=overscanMean,
1329 overscanMedian=overscanMedian,
1330 overscanSigma=overscanSigma,
1331 residualMean=residualMean,
1332 residualMedian=residualMedian,
1333 residualSigma=residualSigma,
1334 overscanBadRowsOrColumns=badRowsOrColumns,
1335 )
1338class ParallelOverscanCorrectionTaskConfig(OverscanCorrectionTaskConfigBase):
1339 doParallelOverscanSaturation = pexConfig.Field(
1340 dtype=bool,
1341 doc="Mask saturated pixels in parallel overscan region?",
1342 default=True,
1343 )
1344 parallelOverscanSaturationLevel = pexConfig.Field(
1345 dtype=float,
1346 doc="The saturation level (adu) to use if not specified in call to "
1347 "maskParallelOverscanAmp. This should be low enough to capture "
1348 "all possible amplifiers for defect detection.",
1349 default=20000.,
1350 )
1351 parallelOverscanSaturationLevelAdjustmentFactor = pexConfig.Field(
1352 dtype=float,
1353 doc="The parallel overscan saturation level may be below that of "
1354 "the data region. This factor is applied to the amplifier "
1355 "saturation value when evaluating saturation in the parallel "
1356 "overscan region.",
1357 default=0.75,
1358 )
1359 parallelOverscanMaskGrowSize = pexConfig.Field(
1360 dtype=int,
1361 doc="Grow the SAT mask in the parallel overscan region by this many pixels. "
1362 "This value was determined from the ITL chip in the LATISS camera.",
1363 default=7,
1364 )
1365 parallelOverscanMaskedColumnGrowSize = pexConfig.Field(
1366 dtype=int,
1367 doc="When a full column is masked in the parallel overscan (at less "
1368 "than saturation) the mask should be grown by this many pixels. "
1369 "This value is determined from ITL chips in LATISS and LSSTCam.",
1370 default=2,
1371 )
1372 leadingToSkip = pexConfig.Field(
1373 dtype=int,
1374 doc="Number of leading values to skip in parallel overscan correction.",
1375 default=0,
1376 )
1377 trailingToSkip = pexConfig.Field(
1378 dtype=int,
1379 doc="Number of trailing values to skip in parallel overscan correction.",
1380 default=0,
1381 )
1382 parallelOverscanFraction = pexConfig.Field(
1383 dtype=float,
1384 doc="When the parallel overscan region median is greater than parallelOverscanFraction "
1385 "and the imaging region median is greater than parallelOverscanImageThreshold "
1386 "then parallel overscan subtraction will be turned off, as this is usually "
1387 "due to the region being flooded with spillover from a super-saturated flat.",
1388 default=0.5,
1389 )
1390 parallelOverscanImageThreshold = pexConfig.Field(
1391 dtype=float,
1392 doc="When the parallel overscan region median is greater than parallelOverscanFraction "
1393 "and the imaging region median is greater than parallelOverscanImageThreshold "
1394 "then parallel overscan subtraction will be turned off, as this is usually "
1395 "due to the region being flooded with spillover from a super-saturated flat.",
1396 default=10000.0,
1397 )
1398 doMedianSmoothingOutlierRejection = pexConfig.Field(
1399 dtype=bool,
1400 doc="Do column-by-column median smoothing outlier rejection? Columns that are rejected "
1401 "in this way will be grown by parallelOverscanMaskedColumnGrowSize.",
1402 default=True,
1403 )
1404 medianSmoothingKernel = pexConfig.Field(
1405 dtype=int,
1406 doc="Kernel (pixels) to use to smooth the parallel overscan columns. Must be odd.",
1407 default=5,
1408 check=lambda x: x // 2 != x / 2,
1409 )
1410 medianSmoothingOutlierThreshold = pexConfig.Field(
1411 dtype=float,
1412 doc="Outlier threshold after parallel median smoothing (adu). This is applied only "
1413 "to positive outliers.",
1414 default=5.0,
1415 check=lambda x: x > 0.0,
1416 )
1418 def setDefaults(self):
1419 super().setDefaults()
1421 # For parallel overscan, this should be a one-sided cut.
1422 self.doAbsoluteMaxDeviation = False
1425class ParallelOverscanCorrectionTask(OverscanCorrectionTaskBase):
1426 """Correction task for parallel overscan.
1428 Parameters
1429 ----------
1430 statControl : `lsst.afw.math.StatisticsControl`, optional
1431 Statistics control object.
1432 """
1433 ConfigClass = ParallelOverscanCorrectionTaskConfig
1434 _DefaultName = "parallelOverscan"
1436 def run(self, exposure, amp, isTransposed=False):
1437 """Measure and remove parallel overscan from an amplifier image.
1439 This method assumes that serial overscan has already been
1440 removed from the amplifier.
1442 Parameters
1443 ----------
1444 exposure : `lsst.afw.image.Exposure`
1445 Image data that will have the overscan corrections applied.
1446 amp : `lsst.afw.cameraGeom.Amplifier`
1447 Amplifier to use for debugging purposes.
1448 isTransposed : `bool`, optional
1449 Is the image transposed, such that serial and parallel
1450 overscan regions are reversed? Default is False.
1452 Returns
1453 -------
1454 overscanResults : `lsst.pipe.base.Struct`
1455 Result struct with components:
1457 ``imageFit``
1458 Value or fit subtracted from the amplifier image data
1459 (scalar or `lsst.afw.image.Image`).
1460 ``overscanFit``
1461 Value or fit subtracted from the parallel overscan image
1462 data (scalar or `lsst.afw.image.Image`).
1463 ``overscanImage``
1464 Image of the parallel overscan region with the parallel
1465 overscan correction applied
1466 (`lsst.afw.image.Image`). This quantity is used to
1467 estimate the amplifier read noise empirically.
1468 ``overscanMean``
1469 Mean of the fit parallel overscan region.
1470 ``overscanMedian``
1471 Median of the fit parallel overscan region.
1472 ``overscanSigma``
1473 Sigma of the fit parallel overscan region.
1474 ``residualMean``
1475 Mean of the residual of the parallel overscan region after
1476 correction.
1477 ``residualMedian``
1478 Median of the residual of the parallel overscan region after
1479 correction.
1480 ``residualSigma``
1481 Mean of the residual of the parallel overscan region after
1482 correction.
1484 Raises
1485 ------
1486 RuntimeError
1487 Raised if an invalid overscan type is set.
1488 """
1489 # This does not need any extending, as we only subtract
1490 # from the data region.
1491 parallelOverscanBBox = amp.getRawParallelOverscanBBox()
1492 imageBBox = amp.getRawDataBBox()
1494 medianSmoothingKernel = self.config.medianSmoothingKernel if \
1495 self.config.doMedianSmoothingOutlierRejection else 0
1497 # The serial overscan correction has removed some signal
1498 # from the parallel overscan region, but that is largely a
1499 # constant offset. The collapseArray method now attempts
1500 # to fill fully masked columns with the median of
1501 # neighboring values, with a fallback to the median of the
1502 # correction in all other columns. Filling with neighbor
1503 # values ensures that large variations in the parallel
1504 # overscan do not create new outlier points. The
1505 # MEDIAN_PER_ROW method does this filling as a separate
1506 # operation, using the same method.
1507 results = self.correctOverscan(
1508 exposure,
1509 amp,
1510 imageBBox,
1511 parallelOverscanBBox,
1512 isTransposed=not isTransposed,
1513 leadingToSkip=self.config.leadingToSkip,
1514 trailingToSkip=self.config.trailingToSkip,
1515 overscanFraction=self.config.parallelOverscanFraction,
1516 imageThreshold=self.config.parallelOverscanImageThreshold,
1517 maskedRowColumnGrowSize=self.config.parallelOverscanMaskedColumnGrowSize,
1518 medianSmoothingKernel=medianSmoothingKernel,
1519 medianSmoothingOutlierThreshold=self.config.medianSmoothingOutlierThreshold,
1520 )
1521 overscanMean = results.overscanMean
1522 overscanMedian = results.overscanMedian
1523 overscanSigma = results.overscanSigma
1524 residualMean = results.overscanMeanResidual
1525 residualMedian = results.overscanMedianResidual
1526 residualSigma = results.overscanSigmaResidual
1527 badRowsOrColumns = results.overscanBadRowsOrColumns
1529 return pipeBase.Struct(
1530 imageFit=results.ampOverscanModel,
1531 overscanFit=results.overscanOverscanModel,
1532 overscanImage=results.overscanImage,
1533 overscanMean=overscanMean,
1534 overscanMedian=overscanMedian,
1535 overscanSigma=overscanSigma,
1536 residualMean=residualMean,
1537 residualMedian=residualMedian,
1538 residualSigma=residualSigma,
1539 badRowsOrColumns=badRowsOrColumns,
1540 )
1542 def maskParallelOverscanAmp(self, exposure, amp, saturationLevel=None):
1543 """Mask parallel overscan, growing saturated pixels.
1545 This operates on the image in-place.
1547 Parameters
1548 ----------
1549 exposure : `lsst.afw.image.Exposure`
1550 An untrimmed raw exposure.
1551 amp : `lsst.afw.cameraGeom.Amplifier`
1552 The amplifier to use for masking.
1553 saturationLevel : `float`, optional
1554 Saturation level to use for masking.
1555 """
1556 if not self.config.doParallelOverscanSaturation:
1557 # This is a no-op.
1558 return
1560 if saturationLevel is None:
1561 saturationLevel = self.config.parallelOverscanSaturationLevel
1563 dataView = afwImage.MaskedImageF(exposure.getMaskedImage(),
1564 amp.getRawParallelOverscanBBox(),
1565 afwImage.PARENT)
1566 # This should mark all of these saturated pixels as SAT.
1567 makeThresholdMask(
1568 maskedImage=dataView,
1569 threshold=saturationLevel,
1570 growFootprints=self.config.parallelOverscanMaskGrowSize,
1571 maskName="SAT",
1572 )