LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
background.py
Go to the documentation of this file.
1 import numpy
2 import itertools
3 from scipy.ndimage import gaussian_filter
4 
5 import lsst.afw.math as afwMath
6 import lsst.afw.image as afwImage
7 import lsst.afw.geom as afwGeom
8 import lsst.afw.cameraGeom as afwCameraGeom
9 import lsst.geom as geom
10 import lsst.meas.algorithms as measAlg
11 import lsst.afw.table as afwTable
12 
13 from lsst.pex.config import Config, Field, ListField, ChoiceField, ConfigField, RangeField, ConfigurableField
14 from lsst.pipe.base import Task
15 
16 
17 def robustMean(array, rej=3.0):
18  """Measure a robust mean of an array
19 
20  Parameters
21  ----------
22  array : `numpy.ndarray`
23  Array for which to measure the mean.
24  rej : `float`
25  k-sigma rejection threshold.
26 
27  Returns
28  -------
29  mean : `array.dtype`
30  Robust mean of `array`.
31  """
32  q1, median, q3 = numpy.percentile(array, [25.0, 50.0, 100.0])
33  good = numpy.abs(array - median) < rej*0.74*(q3 - q1)
34  return array[good].mean()
35 
36 
38  """Configuration for background measurement"""
39  statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points",
40  allowed={"MEANCLIP": "clipped mean",
41  "MEAN": "unclipped mean",
42  "MEDIAN": "median"})
43  xBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in x")
44  yBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in y")
45  algorithm = ChoiceField(dtype=str, default="NATURAL_SPLINE", optional=True,
46  doc="How to interpolate the background values. "
47  "This maps to an enum; see afw::math::Background",
48  allowed={
49  "CONSTANT": "Use a single constant value",
50  "LINEAR": "Use linear interpolation",
51  "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
52  "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust"
53  " to outliers",
54  "NONE": "No background estimation is to be attempted",
55  })
56  mask = ListField(dtype=str, default=["SAT", "BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA"],
57  doc="Names of mask planes to ignore while estimating the background")
58 
59 
61  """Parameters controlling the measurement of sky statistics"""
62  statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points",
63  allowed={"MEANCLIP": "clipped mean",
64  "MEAN": "unclipped mean",
65  "MEDIAN": "median"})
66  clip = Field(doc="Clipping threshold for background", dtype=float, default=3.0)
67  nIter = Field(doc="Clipping iterations for background", dtype=int, default=3)
68  mask = ListField(doc="Mask planes to reject", dtype=str,
69  default=["SAT", "DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA"])
70 
71 
73  """Configuration for SkyMeasurementTask"""
74  skyIter = Field(dtype=int, default=3, doc="k-sigma rejection iterations for sky scale")
75  skyRej = Field(dtype=float, default=3.0, doc="k-sigma rejection threshold for sky scale")
76  background = ConfigField(dtype=BackgroundConfig, doc="Background measurement")
77  xNumSamples = Field(dtype=int, default=4, doc="Number of samples in x for scaling sky frame")
78  yNumSamples = Field(dtype=int, default=4, doc="Number of samples in y for scaling sky frame")
79  stats = ConfigField(dtype=SkyStatsConfig, doc="Measurement of sky statistics in the samples")
80 
81 
82 class SkyMeasurementTask(Task):
83  """Task for creating, persisting and using sky frames
84 
85  A sky frame is like a fringe frame (the sum of many exposures of the night sky,
86  combined with rejection to remove astrophysical objects) except the structure
87  is on larger scales, and hence we bin the images and represent them as a
88  background model (a `lsst.afw.math.BackgroundMI`). The sky frame represents
89  the dominant response of the camera to the sky background.
90  """
91  ConfigClass = SkyMeasurementConfig
92 
93  def getSkyData(self, butler, calibId):
94  """Retrieve sky frame from the butler
95 
96  Parameters
97  ----------
98  butler : `lsst.daf.persistence.Butler`
99  Data butler
100  calibId : `dict`
101  Data identifier for calib
102 
103  Returns
104  -------
105  sky : `lsst.afw.math.BackgroundList`
106  Sky frame
107  """
108  exp = butler.get("sky", calibId)
109  return self.exposureToBackgroundexposureToBackground(exp)
110 
111  @staticmethod
113  """Convert an exposure to background model
114 
115  Calibs need to be persisted as an Exposure, so we need to convert
116  the persisted Exposure to a background model.
117 
118  Parameters
119  ----------
120  bgExp : `lsst.afw.image.Exposure`
121  Background model in Exposure format.
122 
123  Returns
124  -------
125  bg : `lsst.afw.math.BackgroundList`
126  Background model
127  """
128  header = bgExp.getMetadata()
129  xMin = header.getScalar("BOX.MINX")
130  yMin = header.getScalar("BOX.MINY")
131  xMax = header.getScalar("BOX.MAXX")
132  yMax = header.getScalar("BOX.MAXY")
133  algorithm = header.getScalar("ALGORITHM")
134  bbox = geom.Box2I(geom.Point2I(xMin, yMin), geom.Point2I(xMax, yMax))
135  return afwMath.BackgroundList(
136  (afwMath.BackgroundMI(bbox, bgExp.getMaskedImage()),
137  afwMath.stringToInterpStyle(algorithm),
138  afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
139  afwMath.ApproximateControl.UNKNOWN,
140  0, 0, False))
141 
142  def backgroundToExposure(self, statsImage, bbox):
143  """Convert a background model to an exposure
144 
145  Calibs need to be persisted as an Exposure, so we need to convert
146  the background model to an Exposure.
147 
148  Parameters
149  ----------
150  statsImage : `lsst.afw.image.MaskedImageF`
151  Background model's statistics image.
152  bbox : `lsst.geom.Box2I`
153  Bounding box for image.
154 
155  Returns
156  -------
157  exp : `lsst.afw.image.Exposure`
158  Background model in Exposure format.
159  """
160  exp = afwImage.makeExposure(statsImage)
161  header = exp.getMetadata()
162  header.set("BOX.MINX", bbox.getMinX())
163  header.set("BOX.MINY", bbox.getMinY())
164  header.set("BOX.MAXX", bbox.getMaxX())
165  header.set("BOX.MAXY", bbox.getMaxY())
166  header.set("ALGORITHM", self.config.background.algorithm)
167  return exp
168 
169  def measureBackground(self, image):
170  """Measure a background model for image
171 
172  This doesn't use a full-featured background model (e.g., no Chebyshev
173  approximation) because we just want the binning behaviour. This will
174  allow us to average the bins later (`averageBackgrounds`).
175 
176  The `BackgroundMI` is wrapped in a `BackgroundList` so it can be
177  pickled and persisted.
178 
179  Parameters
180  ----------
181  image : `lsst.afw.image.MaskedImage`
182  Image for which to measure background.
183 
184  Returns
185  -------
186  bgModel : `lsst.afw.math.BackgroundList`
187  Background model.
188  """
189  stats = afwMath.StatisticsControl()
190  stats.setAndMask(image.getMask().getPlaneBitMask(self.config.background.mask))
191  stats.setNanSafe(True)
193  self.config.background.algorithm,
194  max(int(image.getWidth()/self.config.background.xBinSize + 0.5), 1),
195  max(int(image.getHeight()/self.config.background.yBinSize + 0.5), 1),
196  "REDUCE_INTERP_ORDER",
197  stats,
198  self.config.background.statistic
199  )
200 
201  bg = afwMath.makeBackground(image, ctrl)
202 
203  return afwMath.BackgroundList((
204  bg,
205  afwMath.stringToInterpStyle(self.config.background.algorithm),
206  afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
207  afwMath.ApproximateControl.UNKNOWN,
208  0, 0, False
209  ))
210 
211  def averageBackgrounds(self, bgList):
212  """Average multiple background models
213 
214  The input background models should be a `BackgroundList` consisting
215  of a single `BackgroundMI`.
216 
217  Parameters
218  ----------
219  bgList : `list` of `lsst.afw.math.BackgroundList`
220  Background models to average.
221 
222  Returns
223  -------
224  bgExp : `lsst.afw.image.Exposure`
225  Background model in Exposure format.
226  """
227  assert all(len(bg) == 1 for bg in bgList), "Mixed bgList: %s" % ([len(bg) for bg in bgList],)
228  images = [bg[0][0].getStatsImage() for bg in bgList]
229  boxes = [bg[0][0].getImageBBox() for bg in bgList]
230  assert len(set((box.getMinX(), box.getMinY(), box.getMaxX(), box.getMaxY()) for box in boxes)) == 1, \
231  "Bounding boxes not all equal"
232  bbox = boxes.pop(0)
233 
234  # Ensure bad pixels are masked
235  maskVal = afwImage.Mask.getPlaneBitMask("BAD")
236  for img in images:
237  bad = numpy.isnan(img.getImage().getArray())
238  img.getMask().getArray()[bad] = maskVal
239 
240  stats = afwMath.StatisticsControl()
241  stats.setAndMask(maskVal)
242  stats.setNanSafe(True)
243  combined = afwMath.statisticsStack(images, afwMath.MEANCLIP, stats)
244 
245  # Set bad pixels to the median
246  # Specifically NOT going to attempt to interpolate the bad values because we're only working on a
247  # single CCD here and can't use the entire field-of-view to do the interpolation (which is what we
248  # would need to avoid introducing problems at the edge of CCDs).
249  array = combined.getImage().getArray()
250  bad = numpy.isnan(array)
251  median = numpy.median(array[~bad])
252  array[bad] = median
253 
254  # Put it into an exposure, which is required for calibs
255  return self.backgroundToExposurebackgroundToExposure(combined, bbox)
256 
257  def measureScale(self, image, skyBackground):
258  """Measure scale of background model in image
259 
260  We treat the sky frame much as we would a fringe frame
261  (except the length scale of the variations is different):
262  we measure samples on the input image and the sky frame,
263  which we will use to determine the scaling factor in the
264  'solveScales` method.
265 
266  Parameters
267  ----------
268  image : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage`
269  Science image for which to measure scale.
270  skyBackground : `lsst.afw.math.BackgroundList`
271  Sky background model.
272 
273  Returns
274  -------
275  imageSamples : `numpy.ndarray`
276  Sample measurements on image.
277  skySamples : `numpy.ndarray`
278  Sample measurements on sky frame.
279  """
280  if isinstance(image, afwImage.Exposure):
281  image = image.getMaskedImage()
282  # Ensure more samples than pixels
283  xNumSamples = min(self.config.xNumSamples, image.getWidth())
284  yNumSamples = min(self.config.yNumSamples, image.getHeight())
285  xLimits = numpy.linspace(0, image.getWidth(), xNumSamples + 1, dtype=int)
286  yLimits = numpy.linspace(0, image.getHeight(), yNumSamples + 1, dtype=int)
287  sky = skyBackground.getImage()
288  maskVal = image.getMask().getPlaneBitMask(self.config.stats.mask)
289  ctrl = afwMath.StatisticsControl(self.config.stats.clip, self.config.stats.nIter, maskVal)
290  statistic = afwMath.stringToStatisticsProperty(self.config.stats.statistic)
291  imageSamples = []
292  skySamples = []
293  for xIndex, yIndex in itertools.product(range(xNumSamples), range(yNumSamples)):
294  # -1 on the stop because Box2I is inclusive of the end point and we don't want to overlap boxes
295  xStart, xStop = xLimits[xIndex], xLimits[xIndex + 1] - 1
296  yStart, yStop = yLimits[yIndex], yLimits[yIndex + 1] - 1
297  box = geom.Box2I(geom.Point2I(xStart, yStart), geom.Point2I(xStop, yStop))
298  subImage = image.Factory(image, box)
299  subSky = sky.Factory(sky, box)
300  imageSamples.append(afwMath.makeStatistics(subImage, statistic, ctrl).getValue())
301  skySamples.append(afwMath.makeStatistics(subSky, statistic, ctrl).getValue())
302  return imageSamples, skySamples
303 
304  def solveScales(self, scales):
305  """Solve multiple scales for a single scale factor
306 
307  Having measured samples from the image and sky frame, we
308  fit for the scaling factor.
309 
310  Parameters
311  ----------
312  scales : `list` of a `tuple` of two `numpy.ndarray` arrays
313  A `list` of the results from `measureScale` method.
314 
315  Returns
316  -------
317  scale : `float`
318  Scale factor.
319  """
320  imageSamples = []
321  skySamples = []
322  for ii, ss in scales:
323  imageSamples.extend(ii)
324  skySamples.extend(ss)
325  assert len(imageSamples) == len(skySamples)
326  imageSamples = numpy.array(imageSamples)
327  skySamples = numpy.array(skySamples)
328 
329  def solve(mask):
330  return afwMath.LeastSquares.fromDesignMatrix(skySamples[mask].reshape(mask.sum(), 1),
331  imageSamples[mask],
332  afwMath.LeastSquares.DIRECT_SVD).getSolution()
333 
334  mask = numpy.isfinite(imageSamples) & numpy.isfinite(skySamples)
335  for ii in range(self.config.skyIter):
336  solution = solve(mask)
337  residuals = imageSamples - solution*skySamples
338  lq, uq = numpy.percentile(residuals[mask], [25, 75])
339  stdev = 0.741*(uq - lq) # Robust stdev from IQR
340  with numpy.errstate(invalid="ignore"): # suppress NAN warnings
341  bad = numpy.abs(residuals) > self.config.skyRej*stdev
342  mask[bad] = False
343 
344  return solve(mask)
345 
346  def subtractSkyFrame(self, image, skyBackground, scale, bgList=None):
347  """Subtract sky frame from science image
348 
349  Parameters
350  ----------
351  image : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage`
352  Science image.
353  skyBackground : `lsst.afw.math.BackgroundList`
354  Sky background model.
355  scale : `float`
356  Scale to apply to background model.
357  bgList : `lsst.afw.math.BackgroundList`
358  List of backgrounds applied to image
359  """
360  if isinstance(image, afwImage.Exposure):
361  image = image.getMaskedImage()
362  if isinstance(image, afwImage.MaskedImage):
363  image = image.getImage()
364  image.scaledMinus(scale, skyBackground.getImage())
365  if bgList is not None:
366  # Append the sky frame to the list of applied background models
367  bgData = list(skyBackground[0])
368  bg = bgData[0]
369  statsImage = bg.getStatsImage().clone()
370  statsImage *= scale
371  newBg = afwMath.BackgroundMI(bg.getImageBBox(), statsImage)
372  newBgData = [newBg] + bgData[1:]
373  bgList.append(newBgData)
374 
375 
376 def interpolate1D(method, xSample, ySample, xInterp):
377  """Interpolate in one dimension
378 
379  Interpolates the curve provided by `xSample` and `ySample` at
380  the positions of `xInterp`. Automatically backs off the
381  interpolation method to achieve successful interpolation.
382 
383  Parameters
384  ----------
385  method : `lsst.afw.math.Interpolate.Style`
386  Interpolation method to use.
387  xSample : `numpy.ndarray`
388  Vector of ordinates.
389  ySample : `numpy.ndarray`
390  Vector of coordinates.
391  xInterp : `numpy.ndarray`
392  Vector of ordinates to which to interpolate.
393 
394  Returns
395  -------
396  yInterp : `numpy.ndarray`
397  Vector of interpolated coordinates.
398 
399  """
400  if len(xSample) == 0:
401  return numpy.ones_like(xInterp)*numpy.nan
402  try:
403  return afwMath.makeInterpolate(xSample.astype(float), ySample.astype(float),
404  method).interpolate(xInterp.astype(float))
405  except Exception:
406  if method == afwMath.Interpolate.CONSTANT:
407  # We've already tried the most basic interpolation and it failed
408  return numpy.ones_like(xInterp)*numpy.nan
409  newMethod = afwMath.lookupMaxInterpStyle(len(xSample))
410  if newMethod == method:
411  newMethod = afwMath.Interpolate.CONSTANT
412  return interpolate1D(newMethod, xSample, ySample, xInterp)
413 
414 
415 def interpolateBadPixels(array, isBad, interpolationStyle):
416  """Interpolate bad pixels in an image array
417 
418  The bad pixels are modified in the array.
419 
420  Parameters
421  ----------
422  array : `numpy.ndarray`
423  Image array with bad pixels.
424  isBad : `numpy.ndarray` of type `bool`
425  Boolean array indicating which pixels are bad.
426  interpolationStyle : `str`
427  Style for interpolation (see `lsst.afw.math.Background`);
428  supported values are CONSTANT, LINEAR, NATURAL_SPLINE,
429  AKIMA_SPLINE.
430  """
431  if numpy.all(isBad):
432  raise RuntimeError("No good pixels in image array")
433  height, width = array.shape
434  xIndices = numpy.arange(width, dtype=float)
435  yIndices = numpy.arange(height, dtype=float)
436  method = afwMath.stringToInterpStyle(interpolationStyle)
437  isGood = ~isBad
438  for y in range(height):
439  if numpy.any(isBad[y, :]) and numpy.any(isGood[y, :]):
440  array[y][isBad[y]] = interpolate1D(method, xIndices[isGood[y]], array[y][isGood[y]],
441  xIndices[isBad[y]])
442 
443  isBad = numpy.isnan(array)
444  isGood = ~isBad
445  for x in range(width):
446  if numpy.any(isBad[:, x]) and numpy.any(isGood[:, x]):
447  array[:, x][isBad[:, x]] = interpolate1D(method, yIndices[isGood[:, x]],
448  array[:, x][isGood[:, x]], yIndices[isBad[:, x]])
449 
450 
452  """Configuration for FocalPlaneBackground
453 
454  Note that `xSize` and `ySize` are floating-point values, as
455  the focal plane frame is usually defined in units of microns
456  or millimetres rather than pixels. As such, their values will
457  need to be revised according to each particular camera. For
458  this reason, no defaults are set for those.
459  """
460  xSize = Field(dtype=float, doc="Bin size in x")
461  ySize = Field(dtype=float, doc="Bin size in y")
462  pixelSize = Field(dtype=float, default=1.0, doc="Pixel size in same units as xSize/ySize")
463  minFrac = Field(dtype=float, default=0.1, doc="Minimum fraction of bin size for good measurement")
464  mask = ListField(dtype=str, doc="Mask planes to treat as bad",
465  default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"])
466  interpolation = ChoiceField(
467  doc="how to interpolate the background values. This maps to an enum; see afw::math::Background",
468  dtype=str, default="AKIMA_SPLINE", optional=True,
469  allowed={
470  "CONSTANT": "Use a single constant value",
471  "LINEAR": "Use linear interpolation",
472  "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
473  "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
474  "NONE": "No background estimation is to be attempted",
475  },
476  )
477  doSmooth = Field(dtype=bool, default=False, doc="Do smoothing?")
478  smoothScale = Field(dtype=float, default=2.0, doc="Smoothing scale, as a multiple of the bin size")
479  binning = Field(dtype=int, default=64, doc="Binning to use for CCD background model (pixels)")
480 
481 
483  """Background model for a focal plane camera
484 
485  We model the background empirically with the "superpixel" method: we
486  measure the background in each superpixel and interpolate between
487  superpixels to yield the model.
488 
489  The principal difference between this and `lsst.afw.math.BackgroundMI`
490  is that here the superpixels are defined in the frame of the focal
491  plane of the camera which removes discontinuities across detectors.
492 
493  The constructor you probably want to use is the `fromCamera` classmethod.
494 
495  There are two use patterns for building a background model:
496 
497  * Serial: create a `FocalPlaneBackground`, then `addCcd` for each of the
498  CCDs in an exposure.
499 
500  * Parallel: create a `FocalPlaneBackground`, then `clone` it for each
501  of the CCDs in an exposure and use those to `addCcd` their respective
502  CCD image. Finally, `merge` all the clones into the original.
503 
504  Once you've built the background model, you can apply it to individual
505  CCDs with the `toCcdBackground` method.
506  """
507  @classmethod
508  def fromCamera(cls, config, camera):
509  """Construct from a camera object
510 
511  Parameters
512  ----------
513  config : `FocalPlaneBackgroundConfig`
514  Configuration for measuring backgrounds.
515  camera : `lsst.afw.cameraGeom.Camera`
516  Camera for which to measure backgrounds.
517  """
518  cameraBox = geom.Box2D()
519  for ccd in camera:
520  for point in ccd.getCorners(afwCameraGeom.FOCAL_PLANE):
521  cameraBox.include(point)
522 
523  width, height = cameraBox.getDimensions()
524  # Offset so that we run from zero
525  offset = geom.Extent2D(cameraBox.getMin())*-1
526  # Add an extra pixel buffer on either side
527  dims = geom.Extent2I(int(numpy.ceil(width/config.xSize)) + 2,
528  int(numpy.ceil(height/config.ySize)) + 2)
529  # Transform takes us from focal plane coordinates --> sample coordinates
530  transform = (geom.AffineTransform.makeTranslation(geom.Extent2D(1, 1))*
531  geom.AffineTransform.makeScaling(1.0/config.xSize, 1.0/config.ySize)*
532  geom.AffineTransform.makeTranslation(offset))
533 
534  return cls(config, dims, afwGeom.makeTransform(transform))
535 
536  def __init__(self, config, dims, transform, values=None, numbers=None):
537  """Constructor
538 
539  Developers should note that changes to the signature of this method
540  require coordinated changes to the `__reduce__` and `clone` methods.
541 
542  Parameters
543  ----------
544  config : `FocalPlaneBackgroundConfig`
545  Configuration for measuring backgrounds.
546  dims : `lsst.geom.Extent2I`
547  Dimensions for background samples.
548  transform : `lsst.afw.geom.TransformPoint2ToPoint2`
549  Transformation from focal plane coordinates to sample coordinates.
550  values : `lsst.afw.image.ImageF`
551  Measured background values.
552  numbers : `lsst.afw.image.ImageF`
553  Number of pixels in each background measurement.
554  """
555  self.configconfig = config
556  self.dimsdims = dims
557  self.transformtransform = transform
558 
559  if values is None:
560  values = afwImage.ImageF(self.dimsdims)
561  values.set(0.0)
562  else:
563  values = values.clone()
564  assert(values.getDimensions() == self.dimsdims)
565  self._values_values = values
566  if numbers is None:
567  numbers = afwImage.ImageF(self.dimsdims) # float for dynamic range and convenience
568  numbers.set(0.0)
569  else:
570  numbers = numbers.clone()
571  assert(numbers.getDimensions() == self.dimsdims)
572  self._numbers_numbers = numbers
573 
574  def __reduce__(self):
575  return self.__class__, (self.configconfig, self.dimsdims, self.transformtransform, self._values_values, self._numbers_numbers)
576 
577  def clone(self):
578  return self.__class__(self.configconfig, self.dimsdims, self.transformtransform, self._values_values, self._numbers_numbers)
579 
580  def addCcd(self, exposure):
581  """Add CCD to model
582 
583  We measure the background on the CCD (clipped mean), and record
584  the results in the model. For simplicity, measurements are made
585  in a box on the CCD corresponding to the warped coordinates of the
586  superpixel rather than accounting for little rotations, etc.
587  We also record the number of pixels used in the measurement so we
588  can have a measure of confidence in each bin's value.
589 
590  Parameters
591  ----------
592  exposure : `lsst.afw.image.Exposure`
593  CCD exposure to measure
594  """
595  detector = exposure.getDetector()
596  transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS),
597  detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE))
598  image = exposure.getMaskedImage()
599  maskVal = image.getMask().getPlaneBitMask(self.configconfig.mask)
600 
601  # Warp the binned image to the focal plane
602  toSample = transform.then(self.transformtransform) # CCD pixels --> focal plane --> sample
603 
604  warped = afwImage.ImageF(self._values_values.getBBox())
605  warpedCounts = afwImage.ImageF(self._numbers_numbers.getBBox())
606  width, height = warped.getDimensions()
607 
608  stats = afwMath.StatisticsControl()
609  stats.setAndMask(maskVal)
610  stats.setNanSafe(True)
611  # Iterating over individual pixels in python is usually bad because it's slow, but there aren't many.
612  pixels = itertools.product(range(width), range(height))
613  for xx, yy in pixels:
614  llc = toSample.applyInverse(geom.Point2D(xx - 0.5, yy - 0.5))
615  urc = toSample.applyInverse(geom.Point2D(xx + 0.5, yy + 0.5))
616  bbox = geom.Box2I(geom.Point2I(llc), geom.Point2I(urc))
617  bbox.clip(image.getBBox())
618  if bbox.isEmpty():
619  continue
620  subImage = image.Factory(image, bbox)
621  result = afwMath.makeStatistics(subImage, afwMath.MEANCLIP | afwMath.NPOINT, stats)
622  mean = result.getValue(afwMath.MEANCLIP)
623  num = result.getValue(afwMath.NPOINT)
624  if not numpy.isfinite(mean) or not numpy.isfinite(num):
625  continue
626  warped[xx, yy, afwImage.LOCAL] = mean*num
627  warpedCounts[xx, yy, afwImage.LOCAL] = num
628 
629  self._values_values += warped
630  self._numbers_numbers += warpedCounts
631 
632  def toCcdBackground(self, detector, bbox):
633  """Produce a background model for a CCD
634 
635  The superpixel background model is warped back to the
636  CCD frame, for application to the individual CCD.
637 
638  Parameters
639  ----------
640  detector : `lsst.afw.cameraGeom.Detector`
641  CCD for which to produce background model.
642  bbox : `lsst.geom.Box2I`
643  Bounding box of CCD exposure.
644 
645  Returns
646  -------
647  bg : `lsst.afw.math.BackgroundList`
648  Background model for CCD.
649  """
650  transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS),
651  detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE))
652  binTransform = (geom.AffineTransform.makeScaling(self.configconfig.binning)*
653  geom.AffineTransform.makeTranslation(geom.Extent2D(0.5, 0.5)))
654 
655  # Binned image on CCD --> unbinned image on CCD --> focal plane --> binned focal plane
656  toSample = afwGeom.makeTransform(binTransform).then(transform).then(self.transformtransform)
657 
658  focalPlane = self.getStatsImagegetStatsImage()
659  fpNorm = afwImage.ImageF(focalPlane.getBBox())
660  fpNorm.set(1.0)
661 
662  image = afwImage.ImageF(bbox.getDimensions()//self.configconfig.binning)
663  norm = afwImage.ImageF(image.getBBox())
664  ctrl = afwMath.WarpingControl("bilinear")
665  afwMath.warpImage(image, focalPlane, toSample.inverted(), ctrl)
666  afwMath.warpImage(norm, fpNorm, toSample.inverted(), ctrl)
667  image /= norm
668 
669  mask = afwImage.Mask(image.getBBox())
670  isBad = numpy.isnan(image.getArray())
671  mask.getArray()[isBad] = mask.getPlaneBitMask("BAD")
672  image.getArray()[isBad] = image.getArray()[~isBad].mean()
673 
674  return afwMath.BackgroundList(
675  (afwMath.BackgroundMI(bbox, afwImage.makeMaskedImage(image, mask)),
676  afwMath.stringToInterpStyle(self.configconfig.interpolation),
677  afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
678  afwMath.ApproximateControl.UNKNOWN,
679  0, 0, False)
680  )
681 
682  def merge(self, other):
683  """Merge with another FocalPlaneBackground
684 
685  This allows multiple background models to be constructed from
686  different CCDs, and then merged to form a single consistent
687  background model for the entire focal plane.
688 
689  Parameters
690  ----------
691  other : `FocalPlaneBackground`
692  Another background model to merge.
693 
694  Returns
695  -------
696  self : `FocalPlaneBackground`
697  The merged background model.
698  """
699  if (self.configconfig.xSize, self.configconfig.ySize) != (other.config.xSize, other.config.ySize):
700  raise RuntimeError("Size mismatch: %s vs %s" % ((self.configconfig.xSize, self.configconfig.ySize),
701  (other.config.xSize, other.config.ySize)))
702  if self.dimsdims != other.dims:
703  raise RuntimeError("Dimensions mismatch: %s vs %s" % (self.dimsdims, other.dims))
704  self._values_values += other._values
705  self._numbers_numbers += other._numbers
706  return self
707 
708  def __iadd__(self, other):
709  """Merge with another FocalPlaneBackground
710 
711  Parameters
712  ----------
713  other : `FocalPlaneBackground`
714  Another background model to merge.
715 
716  Returns
717  -------
718  self : `FocalPlaneBackground`
719  The merged background model.
720  """
721  return self.mergemerge(other)
722 
723  def getStatsImage(self):
724  """Return the background model data
725 
726  This is the measurement of the background for each of the superpixels.
727  """
728  values = self._values_values.clone()
729  values /= self._numbers_numbers
730  thresh = (self.configconfig.minFrac*
731  (self.configconfig.xSize/self.configconfig.pixelSize)*(self.configconfig.ySize/self.configconfig.pixelSize))
732  isBad = self._numbers_numbers.getArray() < thresh
733  if self.configconfig.doSmooth:
734  array = values.getArray()
735  array[:] = smoothArray(array, isBad, self.configconfig.smoothScale)
736  isBad = numpy.isnan(values.array)
737  if numpy.any(isBad):
738  interpolateBadPixels(values.getArray(), isBad, self.configconfig.interpolation)
739  return values
740 
741 
743  """Configuration for MaskObjectsTask"""
744  nIter = Field(dtype=int, default=3, doc="Number of iterations")
745  subtractBackground = ConfigurableField(target=measAlg.SubtractBackgroundTask,
746  doc="Background subtraction")
747  detection = ConfigurableField(target=measAlg.SourceDetectionTask, doc="Source detection")
748  detectSigma = Field(dtype=float, default=5.0, doc="Detection threshold (standard deviations)")
749  doInterpolate = Field(dtype=bool, default=True, doc="Interpolate when removing objects?")
750  interpolate = ConfigurableField(target=measAlg.SubtractBackgroundTask, doc="Interpolation")
751 
752  def setDefaults(self):
753  self.detectiondetection.reEstimateBackground = False
754  self.detectiondetection.doTempLocalBackground = False
755  self.detectiondetection.doTempWideBackground = False
756  self.detectiondetection.thresholdValue = 2.5
757  self.subtractBackgroundsubtractBackground.binSize = 1024
758  self.subtractBackgroundsubtractBackground.useApprox = False
759  self.interpolateinterpolate.binSize = 256
760  self.interpolateinterpolate.useApprox = False
761 
762  def validate(self):
763  if (self.detectiondetection.reEstimateBackground or
764  self.detectiondetection.doTempLocalBackground or
765  self.detectiondetection.doTempWideBackground):
766  raise RuntimeError("Incorrect settings for object masking: reEstimateBackground, "
767  "doTempLocalBackground and doTempWideBackground must be False")
768 
769 
770 class MaskObjectsTask(Task):
771  """Iterative masking of objects on an Exposure
772 
773  This task makes more exhaustive object mask by iteratively doing detection
774  and background-subtraction. The purpose of this task is to get true
775  background removing faint tails of large objects. This is useful to get a
776  clean sky estimate from relatively small number of visits.
777 
778  We deliberately use the specified ``detectSigma`` instead of the PSF,
779  in order to better pick up the faint wings of objects.
780  """
781  ConfigClass = MaskObjectsConfig
782 
783  def __init__(self, *args, **kwargs):
784  super().__init__(*args, **kwargs)
785  # Disposable schema suppresses warning from SourceDetectionTask.__init__
786  self.makeSubtask("detection", schema=afwTable.Schema())
787  self.makeSubtask("interpolate")
788  self.makeSubtask("subtractBackground")
789 
790  def run(self, exposure, maskPlanes=None):
791  """Mask objects on Exposure
792 
793  Objects are found and removed.
794 
795  Parameters
796  ----------
797  exposure : `lsst.afw.image.Exposure`
798  Exposure on which to mask objects.
799  maskPlanes : iterable of `str`, optional
800  List of mask planes to remove.
801  """
802  self.findObjectsfindObjects(exposure)
803  self.removeObjectsremoveObjects(exposure, maskPlanes)
804 
805  def findObjects(self, exposure):
806  """Iteratively find objects on an exposure
807 
808  Objects are masked with the ``DETECTED`` mask plane.
809 
810  Parameters
811  ----------
812  exposure : `lsst.afw.image.Exposure`
813  Exposure on which to mask objects.
814  """
815  for _ in range(self.config.nIter):
816  bg = self.subtractBackground.run(exposure).background
817  self.detection.detectFootprints(exposure, sigma=self.config.detectSigma, clearMask=True)
818  exposure.maskedImage += bg.getImage()
819 
820  def removeObjects(self, exposure, maskPlanes=None):
821  """Remove objects from exposure
822 
823  We interpolate over using a background model if ``doInterpolate`` is
824  set; otherwise we simply replace everything with the median.
825 
826  Parameters
827  ----------
828  exposure : `lsst.afw.image.Exposure`
829  Exposure on which to mask objects.
830  maskPlanes : iterable of `str`, optional
831  List of mask planes to remove. ``DETECTED`` will be added as well.
832  """
833  image = exposure.image
834  mask = exposure.mask
835  maskVal = mask.getPlaneBitMask("DETECTED")
836  if maskPlanes is not None:
837  maskVal |= mask.getPlaneBitMask(maskPlanes)
838  isBad = mask.array & maskVal > 0
839 
840  if self.config.doInterpolate:
841  smooth = self.interpolate.fitBackground(exposure.maskedImage)
842  replace = smooth.getImageF().array[isBad]
843  mask.array &= ~mask.getPlaneBitMask(["DETECTED"])
844  else:
845  replace = numpy.median(image.array[~isBad])
846  image.array[isBad] = replace
847 
848 
849 def smoothArray(array, bad, sigma):
850  """Gaussian-smooth an array while ignoring bad pixels
851 
852  It's not sufficient to set the bad pixels to zero, as then they're treated
853  as if they are zero, rather than being ignored altogether. We need to apply
854  a correction to that image that removes the effect of the bad pixels.
855 
856  Parameters
857  ----------
858  array : `numpy.ndarray` of floating-point
859  Array to smooth.
860  bad : `numpy.ndarray` of `bool`
861  Flag array indicating bad pixels.
862  sigma : `float`
863  Gaussian sigma.
864 
865  Returns
866  -------
867  convolved : `numpy.ndarray`
868  Smoothed image.
869  """
870  convolved = gaussian_filter(numpy.where(bad, 0.0, array), sigma, mode="constant", cval=0.0)
871  denominator = gaussian_filter(numpy.where(bad, 0.0, 1.0), sigma, mode="constant", cval=0.0)
872  return convolved/denominator
int min
int max
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition: Exposure.h:72
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
Pass parameters to a Background object.
Definition: Background.h:56
A class to evaluate image background levels.
Definition: Background.h:434
Pass parameters to a Statistics object.
Definition: Statistics.h:92
Parameters to control convolution.
Definition: warpExposure.h:276
Defines the fields and offsets for a table.
Definition: Schema.h:51
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An integer coordinate rectangle.
Definition: Box.h:55
def __init__(self, config, dims, transform, values=None, numbers=None)
Definition: background.py:536
def __init__(self, *args, **kwargs)
Definition: background.py:783
def removeObjects(self, exposure, maskPlanes=None)
Definition: background.py:820
def run(self, exposure, maskPlanes=None)
Definition: background.py:790
def backgroundToExposure(self, statsImage, bbox)
Definition: background.py:142
def subtractSkyFrame(self, image, skyBackground, scale, bgList=None)
Definition: background.py:346
def measureScale(self, image, skyBackground)
Definition: background.py:257
daf::base::PropertyList * list
Definition: fits.cc:913
daf::base::PropertySet * set
Definition: fits.cc:912
std::shared_ptr< TransformPoint2ToPoint2 > makeTransform(lsst::geom::AffineTransform const &affine)
Wrap an lsst::geom::AffineTransform as a Transform.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition: Exposure.h:462
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT >> image, typename std::shared_ptr< Mask< MaskPixelT >> mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT >> variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
Definition: MaskedImage.h:1240
Interpolate::Style lookupMaxInterpStyle(int const n)
Get the highest order Interpolation::Style available for 'n' points.
Definition: Interpolate.cc:274
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:359
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
Definition: Statistics.cc:738
std::shared_ptr< lsst::afw::image::Image< PixelT > > statisticsStack(std::vector< std::shared_ptr< lsst::afw::image::Image< PixelT >>> &images, Property flags, StatisticsControl const &sctrl=StatisticsControl(), std::vector< lsst::afw::image::VariancePixel > const &wvector=std::vector< lsst::afw::image::VariancePixel >(0))
A function to compute some statistics of a stack of Images.
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
Definition: Interpolate.cc:256
int warpImage(DestImageT &destImage, geom::SkyWcs const &destWcs, SrcImageT const &srcImage, geom::SkyWcs const &srcWcs, WarpingControl const &control, typename DestImageT::SinglePixel padValue=lsst::afw::math::edgePixel< DestImageT >(typename lsst::afw::image::detail::image_traits< DestImageT >::image_category()))
Warp an Image or MaskedImage to a new Wcs.
UndersampleStyle stringToUndersampleStyle(std::string const &style)
Conversion function to switch a string to an UndersampleStyle.
Definition: Background.cc:117
std::shared_ptr< Background > makeBackground(ImageT const &img, BackgroundControl const &bgCtrl)
A convenience function that uses function overloading to make the correct type of Background.
Definition: Background.h:526
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:342
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def interpolate1D(method, xSample, ySample, xInterp)
Definition: background.py:376
def interpolateBadPixels(array, isBad, interpolationStyle)
Definition: background.py:415
def smoothArray(array, bad, sigma)
Definition: background.py:849
def robustMean(array, rej=3.0)
Definition: background.py:17