LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
overscan.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 import numpy as np
23 import lsst.afw.math as afwMath
24 import lsst.afw.image as afwImage
25 import lsst.pipe.base as pipeBase
26 import lsst.pex.config as pexConfig
27 
28 __all__ = ["OverscanCorrectionTaskConfig", "OverscanCorrectionTask"]
29 
30 
31 class OverscanCorrectionTaskConfig(pexConfig.Config):
32  """Overscan correction options.
33  """
34  fitType = pexConfig.ChoiceField(
35  dtype=str,
36  doc="The method for fitting the overscan bias level.",
37  default='MEDIAN',
38  allowed={
39  "POLY": "Fit ordinary polynomial to the longest axis of the overscan region",
40  "CHEB": "Fit Chebyshev polynomial to the longest axis of the overscan region",
41  "LEG": "Fit Legendre polynomial to the longest axis of the overscan region",
42  "NATURAL_SPLINE": "Fit natural spline to the longest axis of the overscan region",
43  "CUBIC_SPLINE": "Fit cubic spline to the longest axis of the overscan region",
44  "AKIMA_SPLINE": "Fit Akima spline to the longest axis of the overscan region",
45  "MEAN": "Correct using the mean of the overscan region",
46  "MEANCLIP": "Correct using a clipped mean of the overscan region",
47  "MEDIAN": "Correct using the median of the overscan region",
48  "MEDIAN_PER_ROW": "Correct using the median per row of the overscan region",
49  },
50  )
51  order = pexConfig.Field(
52  dtype=int,
53  doc=("Order of polynomial to fit if overscan fit type is a polynomial, "
54  "or number of spline knots if overscan fit type is a spline."),
55  default=1,
56  )
57  numSigmaClip = pexConfig.Field(
58  dtype=float,
59  doc="Rejection threshold (sigma) for collapsing overscan before fit",
60  default=3.0,
61  )
62  maskPlanes = pexConfig.ListField(
63  dtype=str,
64  doc="Mask planes to reject when measuring overscan",
65  default=['SAT'],
66  )
67  overscanIsInt = pexConfig.Field(
68  dtype=bool,
69  doc="Treat overscan as an integer image for purposes of fitType=MEDIAN"
70  " and fitType=MEDIAN_PER_ROW.",
71  default=True,
72  )
73 
74 
75 class OverscanCorrectionTask(pipeBase.Task):
76  """Correction task for overscan.
77 
78  This class contains a number of utilities that are easier to
79  understand and use when they are not embedded in nested if/else
80  loops.
81 
82  Parameters
83  ----------
84  statControl : `lsst.afw.math.StatisticsControl`, optional
85  Statistics control object.
86  """
87  ConfigClass = OverscanCorrectionTaskConfig
88  _DefaultName = "overscan"
89 
90  def __init__(self, statControl=None, **kwargs):
91  super().__init__(**kwargs)
92  self.allowDebugallowDebug = True
93 
94  if statControl:
95  self.statControlstatControl = statControl
96  else:
97  self.statControlstatControl = afwMath.StatisticsControl()
98  self.statControlstatControl.setNumSigmaClip(self.config.numSigmaClip)
99  self.statControlstatControl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.maskPlanes))
100 
101  def run(self, ampImage, overscanImage, amp=None):
102  """Measure and remove an overscan from an amplifier image.
103 
104  Parameters
105  ----------
106  ampImage : `lsst.afw.image.Image`
107  Image data that will have the overscan removed.
108  overscanImage : `lsst.afw.image.Image`
109  Overscan data that the overscan is measured from.
110  amp : `lsst.afw.cameraGeom.Amplifier`, optional
111  Amplifier to use for debugging purposes.
112 
113  Returns
114  -------
115  overscanResults : `lsst.pipe.base.Struct`
116  Result struct with components:
117 
118  ``imageFit``
119  Value or fit subtracted from the amplifier image data
120  (scalar or `lsst.afw.image.Image`).
121  ``overscanFit``
122  Value or fit subtracted from the overscan image data
123  (scalar or `lsst.afw.image.Image`).
124  ``overscanImage``
125  Image of the overscan region with the overscan
126  correction applied (`lsst.afw.image.Image`). This
127  quantity is used to estimate the amplifier read noise
128  empirically.
129 
130  Raises
131  ------
132  RuntimeError
133  Raised if an invalid overscan type is set.
134 
135  """
136  if self.config.fitType in ('MEAN', 'MEANCLIP', 'MEDIAN'):
137  overscanResult = self.measureConstantOverscanmeasureConstantOverscan(overscanImage)
138  overscanValue = overscanResult.overscanValue
139  offImage = overscanValue
140  overscanModel = overscanValue
141  maskSuspect = None
142  elif self.config.fitType in ('MEDIAN_PER_ROW', 'POLY', 'CHEB', 'LEG',
143  'NATURAL_SPLINE', 'CUBIC_SPLINE', 'AKIMA_SPLINE'):
144  overscanResult = self.measureVectorOverscanmeasureVectorOverscan(overscanImage)
145  overscanValue = overscanResult.overscanValue
146  maskArray = overscanResult.maskArray
147  isTransposed = overscanResult.isTransposed
148 
149  offImage = afwImage.ImageF(ampImage.getDimensions())
150  offArray = offImage.getArray()
151  overscanModel = afwImage.ImageF(overscanImage.getDimensions())
152  overscanArray = overscanModel.getArray()
153 
154  if hasattr(ampImage, 'getMask'):
155  maskSuspect = afwImage.Mask(ampImage.getDimensions())
156  else:
157  maskSuspect = None
158 
159  if isTransposed:
160  offArray[:, :] = overscanValue[np.newaxis, :]
161  overscanArray[:, :] = overscanValue[np.newaxis, :]
162  if maskSuspect:
163  maskSuspect.getArray()[:, maskArray] |= ampImage.getMask().getPlaneBitMask("SUSPECT")
164  else:
165  offArray[:, :] = overscanValue[:, np.newaxis]
166  overscanArray[:, :] = overscanValue[:, np.newaxis]
167  if maskSuspect:
168  maskSuspect.getArray()[maskArray, :] |= ampImage.getMask().getPlaneBitMask("SUSPECT")
169  else:
170  raise RuntimeError('%s : %s an invalid overscan type' %
171  ("overscanCorrection", self.config.fitType))
172 
173  self.debugViewdebugView(overscanImage, overscanValue, amp)
174 
175  ampImage -= offImage
176  if maskSuspect:
177  ampImage.getMask().getArray()[:, :] |= maskSuspect.getArray()[:, :]
178  overscanImage -= overscanModel
179  return pipeBase.Struct(imageFit=offImage,
180  overscanFit=overscanModel,
181  overscanImage=overscanImage,
182  edgeMask=maskSuspect)
183 
184  @staticmethod
185  def integerConvert(image):
186  """Return an integer version of the input image.
187 
188  Parameters
189  ----------
190  image : `numpy.ndarray`, `lsst.afw.image.Image` or `MaskedImage`
191  Image to convert to integers.
192 
193  Returns
194  -------
195  outI : `numpy.ndarray`, `lsst.afw.image.Image` or `MaskedImage`
196  The integer converted image.
197 
198  Raises
199  ------
200  RuntimeError
201  Raised if the input image could not be converted.
202  """
203  if hasattr(image, "image"):
204  # Is a maskedImage:
205  imageI = image.image.convertI()
206  outI = afwImage.MaskedImageI(imageI, image.mask, image.variance)
207  elif hasattr(image, "convertI"):
208  # Is an Image:
209  outI = image.convertI()
210  elif hasattr(image, "astype"):
211  # Is a numpy array:
212  outI = image.astype(int)
213  else:
214  raise RuntimeError("Could not convert this to integers: %s %s %s",
215  image, type(image), dir(image))
216  return outI
217 
218  # Constant methods
219  def measureConstantOverscan(self, image):
220  """Measure a constant overscan value.
221 
222  Parameters
223  ----------
224  image : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage`
225  Image data to measure the overscan from.
226 
227  Returns
228  -------
229  results : `lsst.pipe.base.Struct`
230  Overscan result with entries:
231  - ``overscanValue``: Overscan value to subtract (`float`)
232  - ``maskArray``: Placeholder for a mask array (`list`)
233  - ``isTransposed``: Orientation of the overscan (`bool`)
234  """
235  if self.config.fitType == 'MEDIAN':
236  calcImage = self.integerConvertintegerConvert(image)
237  else:
238  calcImage = image
239 
240  fitType = afwMath.stringToStatisticsProperty(self.config.fitType)
241  overscanValue = afwMath.makeStatistics(calcImage, fitType, self.statControlstatControl).getValue()
242 
243  return pipeBase.Struct(overscanValue=overscanValue,
244  maskArray=None,
245  isTransposed=False)
246 
247  # Vector correction utilities
248  def getImageArray(self, image):
249  """Extract the numpy array from the input image.
250 
251  Parameters
252  ----------
253  image : `lsst.afw.image.Image` or `lsst.afw.image.MaskedImage`
254  Image data to pull array from.
255 
256  calcImage : `numpy.ndarray`
257  Image data array for numpy operating.
258  """
259  if hasattr(image, "getImage"):
260  calcImage = image.getImage().getArray()
261  calcImage = np.ma.masked_where(image.getMask().getArray() & self.statControlstatControl.getAndMask(),
262  calcImage)
263  else:
264  calcImage = image.getArray()
265  return calcImage
266 
267  @staticmethod
268  def transpose(imageArray):
269  """Transpose input numpy array if necessary.
270 
271  Parameters
272  ----------
273  imageArray : `numpy.ndarray`
274  Image data to transpose.
275 
276  Returns
277  -------
278  imageArray : `numpy.ndarray`
279  Transposed image data.
280  isTransposed : `bool`
281  Indicates whether the input data was transposed.
282  """
283  if np.argmin(imageArray.shape) == 0:
284  return np.transpose(imageArray), True
285  else:
286  return imageArray, False
287 
288  def maskOutliers(self, imageArray):
289  """Mask outliers in a row of overscan data from a robust sigma
290  clipping procedure.
291 
292  Parameters
293  ----------
294  imageArray : `numpy.ndarray`
295  Image to filter along numpy axis=1.
296 
297  Returns
298  -------
299  maskedArray : `numpy.ma.masked_array`
300  Masked image marking outliers.
301  """
302  lq, median, uq = np.percentile(imageArray, [25.0, 50.0, 75.0], axis=1)
303  axisMedians = median
304  axisStdev = 0.74*(uq - lq) # robust stdev
305 
306  diff = np.abs(imageArray - axisMedians[:, np.newaxis])
307  return np.ma.masked_where(diff > self.statControlstatControl.getNumSigmaClip()
308  * axisStdev[:, np.newaxis], imageArray)
309 
310  @staticmethod
311  def collapseArray(maskedArray):
312  """Collapse overscan array (and mask) to a 1-D vector of values.
313 
314  Parameters
315  ----------
316  maskedArray : `numpy.ma.masked_array`
317  Masked array of input overscan data.
318 
319  Returns
320  -------
321  collapsed : `numpy.ma.masked_array`
322  Single dimensional overscan data, combined with the mean.
323  """
324  collapsed = np.mean(maskedArray, axis=1)
325  if collapsed.mask.sum() > 0:
326  collapsed.data[collapsed.mask] = np.mean(maskedArray.data[collapsed.mask], axis=1)
327  return collapsed
328 
329  def collapseArrayMedian(self, maskedArray):
330  """Collapse overscan array (and mask) to a 1-D vector of using the
331  correct integer median of row-values.
332 
333  Parameters
334  ----------
335  maskedArray : `numpy.ma.masked_array`
336  Masked array of input overscan data.
337 
338  Returns
339  -------
340  collapsed : `numpy.ma.masked_array`
341  Single dimensional overscan data, combined with the afwMath median.
342  """
343  integerMI = self.integerConvertintegerConvert(maskedArray)
344 
345  collapsed = []
346  fitType = afwMath.stringToStatisticsProperty('MEDIAN')
347  for row in integerMI:
348  newRow = row.compressed()
349  if len(newRow) > 0:
350  rowMedian = afwMath.makeStatistics(newRow, fitType, self.statControlstatControl).getValue()
351  else:
352  rowMedian = np.nan
353  collapsed.append(rowMedian)
354 
355  return np.array(collapsed)
356 
357  def splineFit(self, indices, collapsed, numBins):
358  """Wrapper function to match spline fit API to polynomial fit API.
359 
360  Parameters
361  ----------
362  indices : `numpy.ndarray`
363  Locations to evaluate the spline.
364  collapsed : `numpy.ndarray`
365  Collapsed overscan values corresponding to the spline
366  evaluation points.
367  numBins : `int`
368  Number of bins to use in constructing the spline.
369 
370  Returns
371  -------
372  interp : `lsst.afw.math.Interpolate`
373  Interpolation object for later evaluation.
374  """
375  if not np.ma.is_masked(collapsed):
376  collapsed.mask = np.array(len(collapsed)*[np.ma.nomask])
377 
378  numPerBin, binEdges = np.histogram(indices, bins=numBins,
379  weights=1 - collapsed.mask.astype(int))
380  with np.errstate(invalid="ignore"):
381  values = np.histogram(indices, bins=numBins,
382  weights=collapsed.data*~collapsed.mask)[0]/numPerBin
383  binCenters = np.histogram(indices, bins=numBins,
384  weights=indices*~collapsed.mask)[0]/numPerBin
385  interp = afwMath.makeInterpolate(binCenters.astype(float)[numPerBin > 0],
386  values.astype(float)[numPerBin > 0],
387  afwMath.stringToInterpStyle(self.config.fitType))
388  return interp
389 
390  @staticmethod
391  def splineEval(indices, interp):
392  """Wrapper function to match spline evaluation API to polynomial fit API.
393 
394  Parameters
395  ----------
396  indices : `numpy.ndarray`
397  Locations to evaluate the spline.
398  interp : `lsst.afw.math.interpolate`
399  Interpolation object to use.
400 
401  Returns
402  -------
403  values : `numpy.ndarray`
404  Evaluated spline values at each index.
405  """
406 
407  return interp.interpolate(indices.astype(float))
408 
409  @staticmethod
410  def maskExtrapolated(collapsed):
411  """Create mask if edges are extrapolated.
412 
413  Parameters
414  ----------
415  collapsed : `numpy.ma.masked_array`
416  Masked array to check the edges of.
417 
418  Returns
419  -------
420  maskArray : `numpy.ndarray`
421  Boolean numpy array of pixels to mask.
422  """
423  maskArray = np.full_like(collapsed, False, dtype=bool)
424  if np.ma.is_masked(collapsed):
425  num = len(collapsed)
426  for low in range(num):
427  if not collapsed.mask[low]:
428  break
429  if low > 0:
430  maskArray[:low] = True
431  for high in range(1, num):
432  if not collapsed.mask[-high]:
433  break
434  if high > 1:
435  maskArray[-high:] = True
436  return maskArray
437 
438  def measureVectorOverscan(self, image):
439  """Calculate the 1-d vector overscan from the input overscan image.
440 
441  Parameters
442  ----------
443  image : `lsst.afw.image.MaskedImage`
444  Image containing the overscan data.
445 
446  Returns
447  -------
448  results : `lsst.pipe.base.Struct`
449  Overscan result with entries:
450  - ``overscanValue``: Overscan value to subtract (`float`)
451  - ``maskArray`` : `list` [ `bool` ]
452  List of rows that should be masked as ``SUSPECT`` when the
453  overscan solution is applied.
454  - ``isTransposed`` : `bool`
455  Indicates if the overscan data was transposed during
456  calcuation, noting along which axis the overscan should be
457  subtracted.
458  """
459  calcImage = self.getImageArraygetImageArray(image)
460 
461  # operate on numpy-arrays from here
462  calcImage, isTransposed = self.transposetranspose(calcImage)
463  masked = self.maskOutliersmaskOutliers(calcImage)
464 
465  if self.config.fitType == 'MEDIAN_PER_ROW':
466  overscanVector = self.collapseArrayMediancollapseArrayMedian(masked)
467  maskArray = self.maskExtrapolatedmaskExtrapolated(overscanVector)
468  else:
469  collapsed = self.collapseArraycollapseArray(masked)
470 
471  num = len(collapsed)
472  indices = 2.0*np.arange(num)/float(num) - 1.0
473 
474  poly = np.polynomial
475  fitter, evaler = {
476  'POLY': (poly.polynomial.polyfit, poly.polynomial.polyval),
477  'CHEB': (poly.chebyshev.chebfit, poly.chebyshev.chebval),
478  'LEG': (poly.legendre.legfit, poly.legendre.legval),
479  'NATURAL_SPLINE': (self.splineFitsplineFit, self.splineEvalsplineEval),
480  'CUBIC_SPLINE': (self.splineFitsplineFit, self.splineEvalsplineEval),
481  'AKIMA_SPLINE': (self.splineFitsplineFit, self.splineEvalsplineEval)
482  }[self.config.fitType]
483 
484  coeffs = fitter(indices, collapsed, self.config.order)
485  overscanVector = evaler(indices, coeffs)
486  maskArray = self.maskExtrapolatedmaskExtrapolated(collapsed)
487  return pipeBase.Struct(overscanValue=np.array(overscanVector),
488  maskArray=maskArray,
489  isTransposed=isTransposed)
490 
491  def debugView(self, image, model, amp=None):
492  """Debug display for the final overscan solution.
493 
494  Parameters
495  ----------
496  image : `lsst.afw.image.Image`
497  Input image the overscan solution was determined from.
498  model : `numpy.ndarray` or `float`
499  Overscan model determined for the image.
500  amp : `lsst.afw.cameraGeom.Amplifier`, optional
501  Amplifier to extract diagnostic information.
502  """
503  import lsstDebug
504  if not lsstDebug.Info(__name__).display:
505  return
506  if not self.allowDebugallowDebug:
507  return
508 
509  calcImage = self.getImageArraygetImageArray(image)
510  calcImage, isTransposed = self.transposetranspose(calcImage)
511  masked = self.maskOutliersmaskOutliers(calcImage)
512  collapsed = self.collapseArraycollapseArray(masked)
513 
514  num = len(collapsed)
515  indices = 2.0 * np.arange(num)/float(num) - 1.0
516 
517  if np.ma.is_masked(collapsed):
518  collapsedMask = collapsed.mask
519  else:
520  collapsedMask = np.array(num*[np.ma.nomask])
521 
522  import matplotlib.pyplot as plot
523  figure = plot.figure(1)
524  figure.clear()
525  axes = figure.add_axes((0.1, 0.1, 0.8, 0.8))
526  axes.plot(indices[~collapsedMask], collapsed[~collapsedMask], 'k+')
527  if collapsedMask.sum() > 0:
528  axes.plot(indices[collapsedMask], collapsed.data[collapsedMask], 'b+')
529  if isinstance(model, np.ndarray):
530  plotModel = model
531  else:
532  plotModel = np.zeros_like(indices)
533  plotModel += model
534  axes.plot(indices, plotModel, 'r-')
535  plot.xlabel("centered/scaled position along overscan region")
536  plot.ylabel("pixel value/fit value")
537  if amp:
538  plot.title(f"{amp.getName()} DataX: "
539  f"[{amp.getRawDataBBox().getBeginX()}:{amp.getRawBBox().getEndX()}]"
540  f"OscanX: [{amp.getRawHorizontalOverscanBBox().getBeginX()}:"
541  f"{amp.getRawHorizontalOverscanBBox().getEndX()}] {self.config.fitType}")
542  else:
543  plot.title("No amp supplied.")
544  figure.show()
545  prompt = "Press Enter or c to continue [chp]..."
546  while True:
547  ans = input(prompt).lower()
548  if ans in ("", " ", "c",):
549  break
550  elif ans in ("p", ):
551  import pdb
552  pdb.set_trace()
553  elif ans in ('x', ):
554  self.allowDebugallowDebug = False
555  break
556  elif ans in ("h", ):
557  print("[h]elp [c]ontinue [p]db e[x]itDebug")
558  plot.close()
table::Key< int > type
Definition: Detector.cc:163
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
Pass parameters to a Statistics object.
Definition: Statistics.h:93
def collapseArrayMedian(self, maskedArray)
Definition: overscan.py:329
def __init__(self, statControl=None, **kwargs)
Definition: overscan.py:90
def debugView(self, image, model, amp=None)
Definition: overscan.py:491
def run(self, ampImage, overscanImage, amp=None)
Definition: overscan.py:101
def splineFit(self, indices, collapsed, numBins)
Definition: overscan.py:357
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
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:354
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:747
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
Definition: Interpolate.cc:257
std::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style=Interpolate::AKIMA_SPLINE)
A factory function to make Interpolate objects.
Definition: Interpolate.cc:343