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