28 from lsst.utils.timer
import timeMethod
33 usePolynomial = pexConfig.Field(
35 doc=
"Fit background difference with Chebychev polynomial interpolation "
36 "(using afw.math.Approximate)? If False, fit with spline interpolation using afw.math.Background",
39 order = pexConfig.Field(
41 doc=
"Order of Chebyshev polynomial background model. Ignored if usePolynomial False",
44 badMaskPlanes = pexConfig.ListField(
45 doc=
"Names of mask planes to ignore while estimating the background",
46 dtype=str, default=[
"NO_DATA",
"DETECTED",
"DETECTED_NEGATIVE",
"SAT",
"BAD",
"INTRP",
"CR"],
49 gridStatistic = pexConfig.ChoiceField(
51 doc=
"Type of statistic to estimate pixel value for the grid points",
56 "MEANCLIP":
"clipped mean"
59 undersampleStyle = pexConfig.ChoiceField(
60 doc=
"Behaviour if there are too few points in grid for requested interpolation style. "
61 "Note: INCREASE_NXNYSAMPLE only allowed for usePolynomial=True.",
63 default=
"REDUCE_INTERP_ORDER",
65 "THROW_EXCEPTION":
"throw an exception if there are too few points",
66 "REDUCE_INTERP_ORDER":
"use an interpolation style with a lower order.",
67 "INCREASE_NXNYSAMPLE":
"Increase the number of samples used to make the interpolation grid.",
70 binSize = pexConfig.Field(
71 doc=
"Bin size for gridding the difference image and fitting a spatial model",
75 interpStyle = pexConfig.ChoiceField(
77 doc=
"Algorithm to interpolate the background values; ignored if usePolynomial is True"
78 "Maps to an enum; see afw.math.Background",
79 default=
"AKIMA_SPLINE",
81 "CONSTANT":
"Use a single constant value",
82 "LINEAR":
"Use linear interpolation",
83 "NATURAL_SPLINE":
"cubic spline with zero second derivative at endpoints",
84 "AKIMA_SPLINE":
"higher-level nonlinear spline that is more robust to outliers",
85 "NONE":
"No background estimation is to be attempted",
88 numSigmaClip = pexConfig.Field(
90 doc=
"Sigma for outlier rejection; ignored if gridStatistic != 'MEANCLIP'.",
93 numIter = pexConfig.Field(
95 doc=
"Number of iterations of outlier rejection; ignored if gridStatistic != 'MEANCLIP'.",
98 bestRefWeightCoverage = pexConfig.RangeField(
100 doc=
"Weight given to coverage (number of pixels that overlap with patch), "
101 "when calculating best reference exposure. Higher weight prefers exposures with high coverage."
102 "Ignored when reference visit is supplied",
106 bestRefWeightVariance = pexConfig.RangeField(
108 doc=
"Weight given to image variance when calculating best reference exposure. "
109 "Higher weight prefers exposures with low image variance. Ignored when reference visit is supplied",
113 bestRefWeightLevel = pexConfig.RangeField(
115 doc=
"Weight given to mean background level when calculating best reference exposure. "
116 "Higher weight prefers exposures with low mean background level. "
117 "Ignored when reference visit is supplied.",
121 approxWeighting = pexConfig.Field(
123 doc=(
"Use inverse-variance weighting when approximating background offset model? "
124 "This will fail when the background offset is constant "
125 "(this is usually only the case in testing with artificial images)."
126 "(usePolynomial=True)"),
129 gridStdevEpsilon = pexConfig.RangeField(
131 doc=
"Tolerance on almost zero standard deviation in a background-offset grid bin. "
132 "If all bins have a standard deviation below this value, the background offset model "
133 "is approximated without inverse-variance weighting. (usePolynomial=True)",
140 ConfigClass = MatchBackgroundsConfig
141 _DefaultName =
"matchBackgrounds"
144 pipeBase.Task.__init__(self, *args, **kwargs)
147 self.
sctrlsctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
148 self.
sctrlsctrl.setNanSafe(
True)
151 def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=None, refImageScaler=None):
152 """Match the backgrounds of a list of coadd temp exposures to a reference coadd temp exposure.
154 Choose a refExpDataRef automatically if none supplied.
156 @param[in] expRefList: list of data references to science exposures to be background-matched;
157 all exposures must exist.
158 @param[in] expDatasetType: dataset type of exposures, e.g. 'goodSeeingCoadd_tempExp'
159 @param[in] imageScalerList: list of image scalers (coaddUtils.ImageScaler);
160 if None then the images are not scaled
161 @param[in] refExpDataRef: data reference for the reference exposure.
162 If None, then this task selects the best exposures from expRefList.
163 if not None then must be one of the exposures in expRefList.
164 @param[in] refImageScaler: image scaler for reference image;
165 ignored if refExpDataRef is None, else scaling is not performed if None
167 @return: a pipBase.Struct containing these fields:
168 - backgroundInfoList: a list of pipeBase.Struct, one per exposure in expRefList,
169 each of which contains these fields:
170 - isReference: this is the reference exposure (only one returned Struct will
171 contain True for this value, unless the ref exposure is listed multiple times)
172 - backgroundModel: differential background model (afw.Math.Background or afw.Math.Approximate).
173 Add this to the science exposure to match the reference exposure.
174 - fitRMS: rms of the fit. This is the sqrt(mean(residuals**2)).
175 - matchedMSE: the MSE of the reference and matched images: mean((refImage - matchedSciImage)**2);
176 should be comparable to difference image's mean variance.
177 - diffImVar: the mean variance of the difference image.
178 All fields except isReference will be None if isReference True or the fit failed.
180 @warning: all exposures must exist on disk
183 numExp = len(expRefList)
185 raise pipeBase.TaskError(
"No exposures to match")
187 if expDatasetType
is None:
188 raise pipeBase.TaskError(
"Must specify expDatasetType")
190 if imageScalerList
is None:
191 self.log.
info(
"imageScalerList is None; no scaling will be performed")
192 imageScalerList = [
None] * numExp
194 if len(expRefList) != len(imageScalerList):
195 raise RuntimeError(
"len(expRefList) = %s != %s = len(imageScalerList)" %
196 (len(expRefList), len(imageScalerList)))
199 if refExpDataRef
is None:
202 expRefList=expRefList,
203 imageScalerList=imageScalerList,
204 expDatasetType=expDatasetType,
206 refExpDataRef = expRefList[refInd]
207 refImageScaler = imageScalerList[refInd]
212 expKeyList = refExpDataRef.butlerSubset.butler.getKeys(expDatasetType)
213 refMatcher =
DataRefMatcher(refExpDataRef.butlerSubset.butler, expDatasetType)
214 refIndSet =
set(refMatcher.matchList(ref0=refExpDataRef, refList=expRefList))
216 if refInd
is not None and refInd
not in refIndSet:
217 raise RuntimeError(
"Internal error: selected reference %s not found in expRefList")
219 refExposure = refExpDataRef.get(expDatasetType, immediate=
True)
220 if refImageScaler
is not None:
221 refMI = refExposure.getMaskedImage()
222 refImageScaler.scaleMaskedImage(refMI)
224 debugIdKeyList = tuple(
set(expKeyList) -
set([
'tract',
'patch']))
226 self.log.
info(
"Matching %d Exposures", numExp)
228 backgroundInfoList = []
229 for ind, (toMatchRef, imageScaler)
in enumerate(zip(expRefList, imageScalerList)):
231 backgroundInfoStruct = pipeBase.Struct(
233 backgroundModel=
None,
239 self.log.
info(
"Matching background of %s to %s", toMatchRef.dataId, refExpDataRef.dataId)
241 toMatchExposure = toMatchRef.get(expDatasetType, immediate=
True)
242 if imageScaler
is not None:
243 toMatchMI = toMatchExposure.getMaskedImage()
244 imageScaler.scaleMaskedImage(toMatchMI)
246 self.
debugDataIdStringdebugDataIdString =
''.join([str(toMatchRef.dataId[vk])
for vk
in debugIdKeyList])
248 refExposure=refExposure,
249 sciExposure=toMatchExposure,
251 backgroundInfoStruct.isReference =
False
252 except Exception
as e:
253 self.log.
warning(
"Failed to fit background %s: %s", toMatchRef.dataId, e)
254 backgroundInfoStruct = pipeBase.Struct(
256 backgroundModel=
None,
262 backgroundInfoList.append(backgroundInfoStruct)
264 return pipeBase.Struct(
265 backgroundInfoList=backgroundInfoList)
269 """Find best exposure to use as the reference exposure
271 Calculate an appropriate reference exposure by minimizing a cost function that penalizes
272 high variance, high background level, and low coverage. Use the following config parameters:
273 - bestRefWeightCoverage
274 - bestRefWeightVariance
277 @param[in] expRefList: list of data references to exposures.
278 Retrieves dataset type specified by expDatasetType.
279 If an exposure is not found, it is skipped with a warning.
280 @param[in] imageScalerList: list of image scalers (coaddUtils.ImageScaler);
281 must be the same length as expRefList
282 @param[in] expDatasetType: dataset type of exposure: e.g. 'goodSeeingCoadd_tempExp'
284 @return: index of best exposure
286 @raise pipeBase.TaskError if none of the exposures in expRefList are found.
288 self.log.
info(
"Calculating best reference visit")
290 meanBkgdLevelList = []
293 if len(expRefList) != len(imageScalerList):
294 raise RuntimeError(
"len(expRefList) = %s != %s = len(imageScalerList)" %
295 (len(expRefList), len(imageScalerList)))
297 for expRef, imageScaler
in zip(expRefList, imageScalerList):
298 exposure = expRef.get(expDatasetType, immediate=
True)
299 maskedImage = exposure.getMaskedImage()
300 if imageScaler
is not None:
302 imageScaler.scaleMaskedImage(maskedImage)
305 varList.append(numpy.nan)
306 meanBkgdLevelList.append(numpy.nan)
307 coverageList.append(numpy.nan)
310 afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, self.
sctrlsctrl)
311 meanVar, meanVarErr = statObjIm.getResult(afwMath.VARIANCE)
312 meanBkgdLevel, meanBkgdLevelErr = statObjIm.getResult(afwMath.MEAN)
313 npoints, npointsErr = statObjIm.getResult(afwMath.NPOINT)
314 varList.append(meanVar)
315 meanBkgdLevelList.append(meanBkgdLevel)
316 coverageList.append(npoints)
318 raise pipeBase.TaskError(
319 "None of the candidate %s exist; cannot select best reference exposure" % (expDatasetType,))
322 varArr = numpy.array(varList)/numpy.nanmax(varList)
323 meanBkgdLevelArr = numpy.array(meanBkgdLevelList)/numpy.nanmax(meanBkgdLevelList)
324 coverageArr = numpy.nanmin(coverageList)/numpy.array(coverageList)
326 costFunctionArr = self.config.bestRefWeightVariance * varArr
327 costFunctionArr += self.config.bestRefWeightLevel * meanBkgdLevelArr
328 costFunctionArr += self.config.bestRefWeightCoverage * coverageArr
329 return numpy.nanargmin(costFunctionArr)
334 Match science exposure's background level to that of reference exposure.
336 Process creates a difference image of the reference exposure minus the science exposure, and then
337 generates an afw.math.Background object. It assumes (but does not require/check) that the mask plane
338 already has detections set. If detections have not been set/masked, sources will bias the
339 background estimation.
340 The 'background' of the difference image is smoothed by spline interpolation (by the Background class)
341 or by polynomial interpolation by the Approximate class. This model of difference image
342 is added to the science exposure in memory.
343 Fit diagnostics are also calculated and returned.
345 @param[in] refExposure: reference exposure
346 @param[in,out] sciExposure: science exposure; modified by changing the background level
347 to match that of the reference exposure
348 @returns a pipBase.Struct with fields:
349 - backgroundModel: an afw.math.Approximate or an afw.math.Background.
350 - fitRMS: rms of the fit. This is the sqrt(mean(residuals**2)).
351 - matchedMSE: the MSE of the reference and matched images: mean((refImage - matchedSciImage)**2);
352 should be comparable to difference image's mean variance.
353 - diffImVar: the mean variance of the difference image.
357 refExposure.writeFits(
lsstDebug.Info(__name__).figpath +
'refExposure.fits')
358 sciExposure.writeFits(
lsstDebug.Info(__name__).figpath +
'sciExposure.fits')
361 if self.config.usePolynomial:
362 x, y = sciExposure.getDimensions()
363 shortSideLength =
min(x, y)
364 if shortSideLength < self.config.binSize:
365 raise ValueError(
"%d = config.binSize > shorter dimension = %d" % (self.config.binSize,
367 npoints = shortSideLength // self.config.binSize
368 if shortSideLength % self.config.binSize != 0:
371 if self.config.order > npoints - 1:
372 raise ValueError(
"%d = config.order > npoints - 1 = %d" % (self.config.order, npoints - 1))
375 if (sciExposure.getDimensions() != refExposure.getDimensions()):
376 wSci, hSci = sciExposure.getDimensions()
377 wRef, hRef = refExposure.getDimensions()
379 "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" %
380 (wSci, hSci, wRef, hRef))
382 statsFlag = getattr(afwMath, self.config.gridStatistic)
383 self.
sctrlsctrl.setNumSigmaClip(self.config.numSigmaClip)
384 self.
sctrlsctrl.setNumIter(self.config.numIter)
386 im = refExposure.getMaskedImage()
387 diffMI = im.Factory(im,
True)
388 diffMI -= sciExposure.getMaskedImage()
390 width = diffMI.getWidth()
391 height = diffMI.getHeight()
392 nx = width // self.config.binSize
393 if width % self.config.binSize != 0:
395 ny = height // self.config.binSize
396 if height % self.config.binSize != 0:
400 bctrl.setUndersampleStyle(self.config.undersampleStyle)
407 if self.config.usePolynomial:
408 order = self.config.order
409 bgX, bgY, bgZ, bgdZ = self.
_gridImage_gridImage(diffMI, self.config.binSize, statsFlag)
410 minNumberGridPoints =
min(len(
set(bgX)), len(
set(bgY)))
412 raise ValueError(
"No overlap with reference. Nothing to match")
413 elif minNumberGridPoints <= self.config.order:
415 if self.config.undersampleStyle ==
"THROW_EXCEPTION":
416 raise ValueError(
"Image does not cover enough of ref image for order and binsize")
417 elif self.config.undersampleStyle ==
"REDUCE_INTERP_ORDER":
418 self.log.
warning(
"Reducing order to %d", (minNumberGridPoints - 1))
419 order = minNumberGridPoints - 1
420 elif self.config.undersampleStyle ==
"INCREASE_NXNYSAMPLE":
421 newBinSize = (minNumberGridPoints*self.config.binSize) // (self.config.order + 1)
422 bctrl.setNxSample(newBinSize)
423 bctrl.setNySample(newBinSize)
425 self.log.
warning(
"Decreasing binsize to %d", newBinSize)
428 isUniformImageDiff =
not numpy.any(bgdZ > self.config.gridStdevEpsilon)
429 weightByInverseVariance =
False if isUniformImageDiff
else self.config.approxWeighting
433 if self.config.usePolynomial:
435 order, order, weightByInverseVariance)
436 undersampleStyle = getattr(afwMath, self.config.undersampleStyle)
437 approx = bkgd.getApproximate(actrl, undersampleStyle)
438 bkgdImage = approx.getImage()
440 bkgdImage = bkgd.getImageF(self.config.interpStyle, self.config.undersampleStyle)
441 except Exception
as e:
442 raise RuntimeError(
"Background/Approximation failed to interp image %s: %s" % (
445 sciMI = sciExposure.getMaskedImage()
451 X, Y, Z, dZ = self.
_gridImage_gridImage(diffMI, self.config.binSize, statsFlag)
452 x0, y0 = diffMI.getXY0()
453 modelValueArr = numpy.empty(len(Z))
454 for i
in range(len(X)):
455 modelValueArr[i] = bkgdImage[int(X[i]-x0), int(Y[i]-y0), afwImage.LOCAL]
456 resids = Z - modelValueArr
457 rms = numpy.sqrt(numpy.mean(resids[~numpy.isnan(resids)]**2))
460 sciExposure.writeFits(
lsstDebug.Info(__name__).figpath +
'sciMatchedExposure.fits')
463 bbox =
geom.Box2D(refExposure.getMaskedImage().getBBox())
465 self.
_debugPlot_debugPlot(X, Y, Z, dZ, bkgdImage, bbox, modelValueArr, resids)
466 except Exception
as e:
467 self.log.
warning(
'Debug plot not generated: %s', e)
470 afwMath.MEANCLIP, self.
sctrlsctrl).getValue()
472 diffIm = diffMI.getImage()
477 outBkgd = approx
if self.config.usePolynomial
else bkgd
479 return pipeBase.Struct(
480 backgroundModel=outBkgd,
485 def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids):
486 """Generate a plot showing the background fit and residuals.
488 It is called when lsstDebug.Info(__name__).savefig = True
489 Saves the fig to lsstDebug.Info(__name__).figpath
490 Displays on screen if lsstDebug.Info(__name__).display = True
492 @param X: array of x positions
493 @param Y: array of y positions
494 @param Z: array of the grid values that were interpolated
495 @param dZ: array of the error on the grid values
496 @param modelImage: image ofthe model of the fit
497 @param model: array of len(Z) containing the grid values predicted by the model
498 @param resids: Z - model
500 import matplotlib.pyplot
as plt
501 import matplotlib.colors
502 from mpl_toolkits.axes_grid1
import ImageGrid
503 zeroIm = afwImage.MaskedImageF(
geom.Box2I(bbox))
505 x0, y0 = zeroIm.getXY0()
506 dx, dy = zeroIm.getDimensions()
508 self.log.
warning(
"No grid. Skipping plot generation.")
510 max, min = numpy.max(Z), numpy.min(Z)
511 norm = matplotlib.colors.normalize(vmax=max, vmin=min)
512 maxdiff = numpy.max(numpy.abs(resids))
513 diffnorm = matplotlib.colors.normalize(vmax=maxdiff, vmin=-maxdiff)
514 rms = numpy.sqrt(numpy.mean(resids**2))
515 fig = plt.figure(1, (8, 6))
516 meanDz = numpy.mean(dZ)
517 grid = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.1,
518 share_all=
True, label_mode=
"L", cbar_mode=
"each",
519 cbar_size=
"7%", cbar_pad=
"2%", cbar_location=
"top")
520 im = grid[0].imshow(zeroIm.getImage().getArray(),
521 extent=(x0, x0+dx, y0+dy, y0), norm=norm,
523 im = grid[0].scatter(X, Y, c=Z, s=15.*meanDz/dZ, edgecolor=
'none', norm=norm,
524 marker=
'o', cmap=
'Spectral')
525 im2 = grid[1].scatter(X, Y, c=resids, edgecolor=
'none', norm=diffnorm,
526 marker=
's', cmap=
'seismic')
527 grid.cbar_axes[0].colorbar(im)
528 grid.cbar_axes[1].colorbar(im2)
529 grid[0].axis([x0, x0+dx, y0+dy, y0])
530 grid[1].axis([x0, x0+dx, y0+dy, y0])
531 grid[0].set_xlabel(
"model and grid")
532 grid[1].set_xlabel(
"residuals. rms = %0.3f"%(rms))
539 def _gridImage(self, maskedImage, binsize, statsFlag):
540 """Private method to grid an image for debugging"""
541 width, height = maskedImage.getDimensions()
542 x0, y0 = maskedImage.getXY0()
543 xedges = numpy.arange(0, width, binsize)
544 yedges = numpy.arange(0, height, binsize)
545 xedges = numpy.hstack((xedges, width))
546 yedges = numpy.hstack((yedges, height))
555 for ymin, ymax
in zip(yedges[0:-1], yedges[1:]):
556 for xmin, xmax
in zip(xedges[0:-1], xedges[1:]):
559 subIm = afwImage.MaskedImageF(maskedImage, subBBox, afwImage.PARENT,
False)
561 afwMath.MEAN | afwMath.MEANCLIP | afwMath.MEDIAN
562 | afwMath.NPOINT | afwMath.STDEV,
564 npoints, _ = stats.getResult(afwMath.NPOINT)
566 stdev, _ = stats.getResult(afwMath.STDEV)
567 if stdev < self.config.gridStdevEpsilon:
568 stdev = self.config.gridStdevEpsilon
569 bgX.append(0.5 * (x0 + xmin + x0 + xmax))
570 bgY.append(0.5 * (y0 + ymin + y0 + ymax))
571 bgdZ.append(stdev/numpy.sqrt(npoints))
572 est, _ = stats.getResult(statsFlag)
575 return numpy.array(bgX), numpy.array(bgY), numpy.array(bgZ), numpy.array(bgdZ)
579 """Match data references for a specified dataset type
581 Note that this is not exact, but should suffice for this task
582 until there is better support for this kind of thing in the butler.
586 """Construct a DataRefMatcher
589 @param[in] datasetType: dataset type to match
592 self.
_keyNames_keyNames = butler.getKeys(datasetType)
594 def _makeKey(self, ref):
595 """Return a tuple of values for the specified keyNames
597 @param[in] ref: data reference
599 @raise KeyError if ref.dataId is missing a key in keyNames
601 return tuple(ref.dataId[key]
for key
in self.
_keyNames_keyNames)
604 """Return True if ref0 == ref1
606 @param[in] ref0: data ref 0
607 @param[in] ref1: data ref 1
609 @raise KeyError if either ID is missing a key in keyNames
614 """Return a list of indices of matches
616 @return tuple of indices of matches
618 @raise KeyError if any ID is missing a key in keyNames
621 return tuple(ind
for ind, ref
in enumerate(refList)
if self.
_makeKey_makeKey(ref) == key0)
Represent a 2-dimensional array of bitmask pixels.
Control how to make an approximation.
Pass parameters to a Background object.
Pass parameters to a Statistics object.
A floating-point coordinate rectangle geometry.
An integer coordinate rectangle.
def __init__(self, butler, datasetType)
def matchList(self, ref0, refList)
def isMatch(self, ref0, ref1)
def selectRefExposure(self, expRefList, imageScalerList, expDatasetType)
def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=None, refImageScaler=None)
def _gridImage(self, maskedImage, binsize, statsFlag)
def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids)
def matchBackgrounds(self, refExposure, sciExposure)
def __init__(self, *args, **kwargs)
daf::base::PropertySet * set
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
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.