332 def matchBackgrounds(self, refExposure, sciExposure):
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.sctrl.setNumSigmaClip(self.config.numSigmaClip)
384 self.sctrl.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(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" % (
443 self.debugDataIdString, e))
445 sciMI = sciExposure.getMaskedImage()
451 X, Y, Z, dZ = self._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(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.sctrl).getValue()
472 diffIm = diffMI.getImage()
477 outBkgd = approx
if self.config.usePolynomial
else bkgd
479 return pipeBase.Struct(
480 backgroundModel=outBkgd,
Control how to make an approximation.
Pass parameters to a Background object.
A floating-point coordinate rectangle geometry.
daf::base::PropertySet * set
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.