LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
isr.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import math
23 
24 import numpy
25 
26 import lsst.afw.geom as afwGeom
27 import lsst.afw.image as afwImage
28 import lsst.afw.detection as afwDetection
29 import lsst.afw.math as afwMath
30 import lsst.meas.algorithms as measAlg
31 import lsst.pex.exceptions as pexExcept
32 
33 def createPsf(fwhm):
34  """Make a double Gaussian PSF
35 
36  @param[in] fwhm FWHM of double Gaussian smoothing kernel
37  @return measAlg.DoubleGaussianPsf
38  """
39  ksize = 4*int(fwhm) + 1
40  return measAlg.DoubleGaussianPsf(ksize, ksize, fwhm/(2*math.sqrt(2*math.log(2))))
41 
42 def calcEffectiveGain(maskedImage):
43  """Calculate effective gain
44 
45  @param[in] maskedImage afw.image.MaskedImage to process
46  @return (median gain, mean gain) in e-/ADU
47  """
48  im = afwImage.ImageF(maskedImage.getImage(), True)
49  var = maskedImage.getVariance()
50  im /= var
51  medgain = afwMath.makeStatistics(im, afwMath.MEDIAN).getValue()
52  meangain = afwMath.makeStatistics(im, afwMath.MEANCLIP).getValue()
53  return medgain, meangain
54 
55 def transposeMaskedImage(maskedImage):
56  """Make a transposed copy of a masked image
57 
58  @param[in] maskedImage afw.image.MaskedImage to process
59  @return transposed masked image
60  """
61  transposed = maskedImage.Factory(afwGeom.Extent2I(maskedImage.getHeight(), maskedImage.getWidth()))
62  transposed.getImage().getArray()[:] = maskedImage.getImage().getArray().T
63  transposed.getMask().getArray()[:] = maskedImage.getMask().getArray().T
64  transposed.getVariance().getArray()[:] = maskedImage.getVariance().getArray().T
65  return transposed
66 
67 def interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=None):
68  """Interpolate over defects specified in a defect list
69 
70  @param[in,out] maskedImage masked image to process
71  @param[in] defectList defect list
72  @param[in] fwhm FWHM of double Gaussian smoothing kernel
73  @param[in] fallbackValue fallback value if an interpolated value cannot be determined;
74  if None then use clipped mean image value
75  """
76  psf = createPsf(fwhm)
77  if fallbackValue is None:
78  fallbackValue = afwMath.makeStatistics(maskedImage.getImage(), afwMath.MEANCLIP).getValue()
79  if 'INTRP' not in maskedImage.getMask().getMaskPlaneDict().keys():
80  maskedImage.getMask.addMaskPlane('INTRP')
81  measAlg.interpolateOverDefects(maskedImage, psf, defectList, fallbackValue, True)
82 
83 def defectListFromFootprintList(fpList, growFootprints=1):
84  """Compute a defect list from a footprint list, optionally growing the footprints
85 
86  @param[in] fpList footprint list
87  @param[in] growFootprints amount by which to grow footprints of detected regions
88  @return meas.algorithms.DefectListT
89  """
90  defectList = measAlg.DefectListT()
91  for fp in fpList:
92  if growFootprints > 0:
93  # if "True", growing requires a convolution
94  # if "False", its faster
95  fpGrow = afwDetection.growFootprint(fp, growFootprints, False)
96  else:
97  fpGrow = fp
98  for bbox in afwDetection.footprintToBBoxList(fpGrow):
99  defect = measAlg.Defect(bbox)
100  defectList.push_back(defect)
101  return defectList
102 
103 def transposeDefectList(defectList):
104  """Make a transposed copy of a defect list
105 
106  @param[in] defectList defect list
107  @return meas.algorithms.DefectListT with transposed defects
108  """
109  retDefectList = measAlg.DefectListT()
110  for defect in defectList:
111  bbox = defect.getBBox()
112  nbbox = afwGeom.Box2I(afwGeom.Point2I(bbox.getMinY(), bbox.getMinX()),
113  afwGeom.Extent2I(bbox.getDimensions()[1], bbox.getDimensions()[0]))
114  retDefectList.push_back(measAlg.Defect(nbbox))
115  return retDefectList
116 
117 def maskPixelsFromDefectList(maskedImage, defectList, maskName='BAD'):
118  """Set mask plane based on a defect list
119 
120  @param[in,out] maskedImage afw.image.MaskedImage to process; mask plane is updated
121  @param[in] defectList meas.algorithms.DefectListT
122  @param[in] maskName mask plane name
123  """
124  # mask bad pixels
125  mask = maskedImage.getMask()
126  bitmask = mask.getPlaneBitMask(maskName)
127  for defect in defectList:
128  bbox = defect.getBBox()
130 
131 def getDefectListFromMask(maskedImage, maskName, growFootprints=1):
132  """Compute a defect list from a specified mask plane
133 
134  @param[in] maskedImage masked image to process
135  @param[in] maskName mask plane name, or list of names
136  @param[in] growFootprints amount by which to grow footprints of detected regions
137  @return meas.algrithms.DefectListT of regions in mask
138  """
139  mask = maskedImage.getMask()
140  workmask = afwImage.MaskU(mask, True)
141  workmask &= mask.getPlaneBitMask(maskName)
142  thresh = afwDetection.Threshold(0.5)
143  maskimg = afwImage.ImageU(workmask.getBBox())
144  maskimg <<= workmask
145  ds = afwDetection.FootprintSet(maskimg, thresh)
146  fpList = ds.getFootprints()
147  return defectListFromFootprintList(fpList, growFootprints)
148 
149 def makeThresholdMask(maskedImage, threshold, growFootprints=1, maskName = 'SAT'):
150  """Mask pixels based on threshold detection
151 
152  @param[in,out] maskedImage afw.image.MaskedImage to process; the mask is altered
153  @param[in] threshold detection threshold
154  @param[in] growFootprints amount by which to grow footprints of detected regions
155  @param[in] maskName mask plane name
156  @return meas.algorihtms.DefectListT of regions set in the mask.
157  """
158  # find saturated regions
159  thresh = afwDetection.Threshold(threshold)
160  fs = afwDetection.FootprintSet(maskedImage, thresh)
161 
162  if growFootprints > 0:
163  fs = afwDetection.FootprintSet(fs, growFootprints)
164 
165  fpList = fs.getFootprints()
166  # set mask
167  mask = maskedImage.getMask()
168  bitmask = mask.getPlaneBitMask(maskName)
169  afwDetection.setMaskFromFootprintList(mask, fpList, bitmask)
170 
171  return defectListFromFootprintList(fpList, growFootprints=0)
172 
173 def interpolateFromMask(maskedImage, fwhm, growFootprints=1, maskName='SAT', fallbackValue=None):
174  """Interpolate over defects identified by a particular mask plane
175 
176  @param[in,out] maskedImage afw.image.MaskedImage to process
177  @param[in] fwhm FWHM of double Gaussian smoothing kernel
178  @param[in] growFootprints amount by which to grow footprints of detected regions
179  @param[in] maskName mask plane name
180  @param[in] fallbackValue value of last resort for interpolation
181  """
182  defectList = getDefectListFromMask(maskedImage, maskName, growFootprints)
183  interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue)
184 
185 def saturationCorrection(maskedImage, saturation, fwhm, growFootprints=1, interpolate=True, maskName='SAT',
186  fallbackValue=None):
187  """Mark saturated pixels and optionally interpolate over them
188 
189  @param[in,out] maskedImage afw.image.MaskedImage to process
190  @param[in] saturation saturation level (used as a detection threshold)
191  @param[in] fwhm FWHM of double Gaussian smoothing kernel
192  @param[in] growFootprints amount by which to grow footprints of detected regions
193  @param[in] interpolate interpolate over saturated pixels?
194  @param[in] maskName mask plane name
195  @param[in] fallbackValue value of last resort for interpolation
196  """
197  defectList = makeThresholdMask(
198  maskedImage = maskedImage,
199  threshold = saturation,
200  growFootprints = growFootprints,
201  maskName = maskName,
202  )
203  if interpolate:
204  interpolateDefectList(maskedImage, defectList, fwhm, fallbackValue=fallbackValue)
205 
206 def biasCorrection(maskedImage, biasMaskedImage):
207  """Apply bias correction in place
208 
209  @param[in,out] maskedImage masked image to correct
210  @param[in] biasMaskedImage bias, as a masked image
211  """
212  if maskedImage.getBBox(afwImage.LOCAL) != biasMaskedImage.getBBox(afwImage.LOCAL):
213  raise RuntimeError("maskedImage bbox %s != biasMaskedImage bbox %s" % \
214  (maskedImage.getBBox(afwImage.LOCAL), biasMaskedImage.getBBox(afwImage.LOCAL)))
215  maskedImage -= biasMaskedImage
216 
217 def darkCorrection(maskedImage, darkMaskedImage, expScale, darkScale):
218  """Apply dark correction in place
219 
220  maskedImage -= dark * expScaling / darkScaling
221 
222  @param[in,out] maskedImage afw.image.MaskedImage to correct
223  @param[in] darkMaskedImage dark afw.image.MaskedImage
224  @param[in] expScale exposure scale
225  @param[in] darkScale dark scale
226  """
227  if maskedImage.getBBox(afwImage.LOCAL) != darkMaskedImage.getBBox(afwImage.LOCAL):
228  raise RuntimeError("maskedImage bbox %s != darkMaskedImage bbox %s" % \
229  (maskedImage.getBBox(afwImage.LOCAL), darkMaskedImage.getBBox(afwImage.LOCAL)))
230 
231  scale = expScale / darkScale
232  maskedImage.scaledMinus(scale, darkMaskedImage)
233 
234 def updateVariance(maskedImage, gain, readNoise):
235  """Set the variance plane based on the image plane
236 
237  @param[in,out] maskedImage afw.image.MaskedImage; image plane is read and variance plane is written
238  @param[in] gain amplifier gain (e-/ADU)
239  @param[in] readNoise amplifier read noise (ADU/pixel)
240  """
241  var = maskedImage.getVariance()
242  var <<= maskedImage.getImage()
243  var /= gain
244  var += readNoise**2
245 
246 def flatCorrection(maskedImage, flatMaskedImage, scalingType, userScale=1.0):
247  """Apply flat correction in place
248 
249  @param[in,out] maskedImage afw.image.MaskedImage to correct
250  @param[in] flatMaskedImage flat field afw.image.MaskedImage
251  @param[in] scalingType how to compute flat scale; one of 'MEAN', 'MEDIAN' or 'USER'
252  @param[in] userScale scale to use if scalingType is 'USER', else ignored
253  """
254  if maskedImage.getBBox(afwImage.LOCAL) != flatMaskedImage.getBBox(afwImage.LOCAL):
255  raise RuntimeError("maskedImage bbox %s != flatMaskedImage bbox %s" % \
256  (maskedImage.getBBox(afwImage.LOCAL), flatMaskedImage.getBBox(afwImage.LOCAL)))
257 
258  # Figure out scale from the data
259  # Ideally the flats are normalized by the calibration product pipelin, but this allows some flexibility
260  # in the case that the flat is created by some other mechanism.
261  if scalingType == 'MEAN':
262  flatScale = afwMath.makeStatistics(flatMaskedImage.getImage(), afwMath.MEAN).getValue(afwMath.MEAN)
263  elif scalingType == 'MEDIAN':
264  flatScale = afwMath.makeStatistics(flatMaskedImage.getImage(), afwMath.MEDIAN).getValue(afwMath.MEDIAN)
265  elif scalingType == 'USER':
266  flatScale = userScale
267  else:
268  raise pexExcept.Exception, '%s : %s not implemented' % ("flatCorrection", scalingType)
269 
270  maskedImage.scaledDivides(1.0/flatScale, flatMaskedImage)
271 
272 def illuminationCorrection(maskedImage, illumMaskedImage, illumScale):
273  """Apply illumination correction in place
274 
275  @param[in,out] maskedImage afw.image.MaskedImage to correct
276  @param[in] illumMaskedImage illumination correction masked image
277  @param[in] illumScale scale value for illumination correction
278  """
279  if maskedImage.getBBox(afwImage.LOCAL) != illumMaskedImage.getBBox(afwImage.LOCAL):
280  raise RuntimeError("maskedImage bbox %s != illumMaskedImage bbox %s" % \
281  (maskedImage.getBBox(afwImage.LOCAL), illumMaskedImage.getBBox(afwImage.LOCAL)))
282 
283  maskedImage.scaledDivides(1./illumScale, illumMaskedImage)
284 
285 def overscanCorrection(ampMaskedImage, overscanImage, fitType='MEDIAN', order=1, collapseRej=3.0,
286  statControl=None):
287  """Apply overscan correction in place
288 
289  @param[in,out] ampMaskedImage masked image to correct
290  @param[in] overscanImage overscan data as an afw.image.Image or afw.image.MaskedImage.
291  If a masked image is passed in the mask plane will be used
292  to constrain the fit of the bias level.
293  @param[in] fitType type of fit for overscan correction; one of:
294  - 'MEAN'
295  - 'MEDIAN'
296  - 'POLY' (ordinary polynomial)
297  - 'CHEB' (Chebyshev polynomial)
298  - 'LEG' (Legendre polynomial)
299  - 'NATURAL_SPLINE', 'CUBIC_SPLINE', 'AKIMA_SPLINE' (splines)
300  @param[in] order polynomial order or spline knots (ignored unless fitType
301  indicates a polynomial or spline)
302  @param[in] collapseRej Rejection threshold (sigma) for collapsing dimension of overscan
303  @param[in] statControl Statistics control object
304  """
305  ampImage = ampMaskedImage.getImage()
306  if statControl is None:
307  statControl = afwMath.StatisticsControl()
308  if fitType == 'MEAN':
309  offImage = afwMath.makeStatistics(overscanImage, afwMath.MEAN, statControl).getValue(afwMath.MEAN)
310  elif fitType == 'MEDIAN':
311  offImage = afwMath.makeStatistics(overscanImage, afwMath.MEDIAN, statControl).getValue(afwMath.MEDIAN)
312  elif fitType in ('POLY', 'CHEB', 'LEG', 'NATURAL_SPLINE', 'CUBIC_SPLINE', 'AKIMA_SPLINE'):
313  if hasattr(overscanImage, "getImage"):
314  biasArray = overscanImage.getImage().getArray()
315  biasArray = numpy.ma.masked_where(overscanImage.getMask().getArray() & statControl.getAndMask(),
316  biasArray)
317  else:
318  biasArray = overscanImage.getArray()
319  # Fit along the long axis, so collapse along each short row and fit the resulting array
320  shortInd = numpy.argmin(biasArray.shape)
321  if shortInd == 0:
322  # Convert to some 'standard' representation to make things easier
323  biasArray = numpy.transpose(biasArray)
324 
325  # Do a single round of clipping to weed out CR hits and signal leaking into the overscan
326  percentiles = numpy.percentile(biasArray, [25.0, 50.0, 75.0], axis=1)
327  medianBiasArr = percentiles[1]
328  stdevBiasArr = 0.74*(percentiles[2] - percentiles[0]) # robust stdev
329  diff = numpy.abs(biasArray - medianBiasArr[:,numpy.newaxis])
330  biasMaskedArr = numpy.ma.masked_where(diff > collapseRej*stdevBiasArr[:,numpy.newaxis], biasArray)
331  collapsed = numpy.mean(biasMaskedArr, axis=1)
332  del biasArray, percentiles, stdevBiasArr, diff, biasMaskedArr
333 
334  if shortInd == 0:
335  collapsed = numpy.transpose(collapsed)
336 
337  num = len(collapsed)
338  indices = 2.0*numpy.arange(num)/float(num) - 1.0
339 
340  if fitType in ('POLY', 'CHEB', 'LEG'):
341  # A numpy polynomial
342  poly = numpy.polynomial
343  fitter, evaler = {"POLY": (poly.polynomial.polyfit, poly.polynomial.polyval),
344  "CHEB": (poly.chebyshev.chebfit, poly.chebyshev.chebval),
345  "LEG": (poly.legendre.legfit, poly.legendre.legval),
346  }[fitType]
347 
348  coeffs = fitter(indices, collapsed, order)
349  fitBiasArr = evaler(indices, coeffs)
350  elif 'SPLINE' in fitType:
351  # An afw interpolation
352  numBins = order
353  #
354  # numpy.histogram needs a real array for the mask, but numpy.ma "optimises" the case
355  # no-values-are-masked by replacing the mask array by a scalar, numpy.ma.nomask
356  #
357  # Issue DM-415
358  #
359  collapsedMask = collapsed.mask
360  try:
361  if collapsedMask == numpy.ma.nomask:
362  collapsedMask = numpy.array(len(collapsed)*[numpy.ma.nomask])
363  except ValueError: # If collapsedMask is an array the test fails [needs .all()]
364  pass
365 
366  numPerBin, binEdges = numpy.histogram(indices, bins=numBins,
367  weights=1-collapsedMask.astype(int))
368  # Binning is just a histogram, with weights equal to the values.
369  # Use a similar trick to get the bin centers (this deals with different numbers per bin).
370  values = numpy.histogram(indices, bins=numBins, weights=collapsed)[0]/numPerBin
371  binCenters = numpy.histogram(indices, bins=numBins, weights=indices)[0]/numPerBin
372  interp = afwMath.makeInterpolate(binCenters.astype(float)[numPerBin > 0],
373  values.astype(float)[numPerBin > 0],
375  fitBiasArr = numpy.array([interp.interpolate(i) for i in indices])
376 
377  import lsstDebug
378  if lsstDebug.Info(__name__).display:
379  import matplotlib.pyplot as plot
380  figure = plot.figure(1)
381  figure.clear()
382  axes = figure.add_axes((0.1, 0.1, 0.8, 0.8))
383  axes.plot(indices, collapsed, 'k+')
384  axes.plot(indices, fitBiasArr, 'r-')
385  figure.show()
386  prompt = "Press Enter or c to continue [chp]... "
387  while True:
388  ans = raw_input(prompt).lower()
389  if ans in ("", "c",):
390  break
391  if ans in ("p",):
392  import pdb; pdb.set_trace()
393  elif ans in ("h", ):
394  print "h[elp] c[ontinue] p[db]"
395  figure.close()
396 
397  offImage = ampImage.Factory(ampImage.getDimensions())
398  offArray = offImage.getArray()
399  if shortInd == 1:
400  offArray[:,:] = fitBiasArr[:,numpy.newaxis]
401  else:
402  offArray[:,:] = fitBiasArr[numpy.newaxis,:]
403 
404  # We don't trust any extrapolation: mask those pixels as SUSPECT
405  # This will occur when the top and or bottom edges of the overscan
406  # contain saturated values. The values will be extrapolated from
407  # the surrounding pixels, but we cannot entirely trust the value of
408  # the extrapolation, and will mark the image mask plane to flag the
409  # image as such.
410  mask = ampMaskedImage.getMask()
411  maskArray = mask.getArray() if shortInd == 1 else mask.getArray().transpose()
412  suspect = mask.getPlaneBitMask("SUSPECT")
413  try:
414  if collapsed.mask == numpy.ma.nomask:
415  # There is no mask, so the whole array is fine
416  pass
417  except ValueError: # If collapsed.mask is an array the test fails [needs .all()]
418  for low in xrange(num):
419  if not collapsed.mask[low]:
420  break
421  if low > 0:
422  maskArray[:low,:] |= suspect
423  for high in xrange(1, num):
424  if not collapsed.mask[-high]:
425  break
426  if high > 1:
427  maskArray[-high:,:] |= suspect
428 
429  else:
430  raise pexExcept.Exception, '%s : %s an invalid overscan type' % \
431  ("overscanCorrection", fitType)
432  ampImage -= offImage
433 
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
Definition: Interpolate.cc:262
def updateVariance
Definition: isr.py:234
boost::shared_ptr< Footprint > growFootprint(Footprint const &foot, int nGrow, bool left, bool right, bool up, bool down)
Grow a Footprint in at least one of the cardinal directions, returning a new Footprint.
def flatCorrection
Definition: isr.py:246
boost::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style=Interpolate::AKIMA_SPLINE)
Definition: Interpolate.cc:353
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:44
def transposeDefectList
Definition: isr.py:103
An integer coordinate rectangle.
Definition: Box.h:53
def maskPixelsFromDefectList
Definition: isr.py:117
def saturationCorrection
Definition: isr.py:186
void interpolateOverDefects(MaskedImageT &image, lsst::afw::detection::Psf const &psf, std::vector< Defect::Ptr > &badList, double fallbackValue=0.0, bool useFallbackValueAtEdge=false)
Process a set of known bad pixels in an image.
Definition: Interp.cc:2058
Represent a Psf as a circularly symmetrical double Gaussian.
def interpolateDefectList
Definition: isr.py:67
Pass parameters to a Statistics objectA class to pass parameters which control how the stats are calc...
Definition: Statistics.h:92
def illuminationCorrection
Definition: isr.py:272
def calcEffectiveGain
Definition: isr.py:42
def createPsf
Definition: isr.py:33
def makeThresholdMask
Definition: isr.py:149
A set of pixels in an Image.
Definition: Footprint.h:62
def getDefectListFromMask
Definition: isr.py:131
def transposeMaskedImage
Definition: isr.py:55
def overscanCorrection
Definition: isr.py:286
std::vector< lsst::afw::geom::Box2I > footprintToBBoxList(Footprint const &foot)
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
MaskT setMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels that are in the Footprint.
def interpolateFromMask
Definition: isr.py:173
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
MaskT setMaskFromFootprintList(lsst::afw::image::Mask< MaskT > *mask, boost::shared_ptr< std::vector< boost::shared_ptr< Footprint >> const > const &footprints, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels which are in the set of Footprints.
def darkCorrection
Definition: isr.py:217
Encapsulate information about a bad portion of a detector.
Definition: Interp.h:70
def biasCorrection
Definition: isr.py:206
def defectListFromFootprintList
Definition: isr.py:83