From 4260e5d2c5b4e9bbcdfc31f53bb20926b57be238 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Fri, 28 Jun 2024 05:49:42 -0700 Subject: [PATCH 01/13] Update MatchBackgroundsTask to Gen3, first pass Initial version of this task is written using old architecture, and so needs updating. In MatchBackgroundsTask, and its method selectRefExposure, required parameters were equally outdated: DataId, DatasetType, and ImageScaler. All of these now seem consolidated under lsst.afw.image.Exposure, so separate calls to DataId and DatasetType are now single calls for Exposure objects. ImageScaler calls were replaced in-line with Exposure.getPhotoCalib() calls, to scale all image flux to the same zeropoint (nJy). Also, we want to process visit-level images using this, so a MatchBackgroundsConnections class was created, MatchBackgroundsConfig was updated to inherit from PipelineTaskConfig (and those connections), and a rudimentary runQuantum method was added to MatchBackgroundsTask. --- python/lsst/pipe/tasks/matchBackgrounds.py | 196 ++++++++++++--------- 1 file changed, 108 insertions(+), 88 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 210ecc209..405339777 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -30,8 +30,42 @@ import lsstDebug from lsst.utils.timer import timeMethod +from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections +import lsst.pipe.base.connectionTypes as cT -class MatchBackgroundsConfig(pexConfig.Config): +import pdb + + +class MatchBackgroundsConnections(PipelineTaskConnections, + defaultTemplates={}, + dimensions={"instrument",}): + + calExpList = cT.Input( + doc="Input exposures to process", + dimensions=("visit", "detector",), + storageClass="ExposureF", + name="calexp", + multiple=True, # Probably this will be needed for the List format + deferLoad=True, # All image calls currently use .get(), so... + ) + # refCalExp = cT.Input( + # doc="Reference exposure", + # dimensions=("visit", "detector", "band"), + # storageClass="ExposureF", + # name="calexp", + # ) + backgroundInfoList = cT.Output( + doc="List of differential backgrounds, w/goodness of fit params", + dimensions=("visit", "detector",), + storageClass="Background", + name="calexpBackground_diff", + multiple=True, + ) + + +# class MatchBackgroundsConfig(pexConfig.Config): +class MatchBackgroundsConfig(PipelineTaskConfig, + pipelineConnections=MatchBackgroundsConnections): usePolynomial = pexConfig.Field( dtype=bool, @@ -139,40 +173,40 @@ class MatchBackgroundsConfig(pexConfig.Config): ) -class MatchBackgroundsTask(pipeBase.Task): +class MatchBackgroundsTask(pipeBase.PipelineTask): ConfigClass = MatchBackgroundsConfig _DefaultName = "matchBackgrounds" - def __init__(self, *args, **kwargs): - pipeBase.Task.__init__(self, *args, **kwargs) + # def __init__(self, *args, **kwargs): + # pipeBase.Task.__init__(self, *args, **kwargs) - self.sctrl = afwMath.StatisticsControl() - self.sctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)) - self.sctrl.setNanSafe(True) + # self.sctrl = afwMath.StatisticsControl() + # self.sctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)) + # self.sctrl.setNanSafe(True) + + def runQuantum(self, butlerQC, inputRefs, outputRefs): + ''' + Boilerplate for now, until bug-testing commences + ''' + inputs = butlerQC.get(inputRefs) + outputs = self.run(**inputs) + butlerQC.put(outputs, outputRefs) @timeMethod - def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=None, refImageScaler=None): - """Match the backgrounds of a list of coadd temp exposures to a reference coadd temp exposure. + def run(self, calExpList): + """Match the backgrounds of a list of science exposures to a reference science exposure. - Choose a refExpDataRef automatically if none supplied. + Choose a refCalExp automatically if none supplied. Parameters ---------- - expRefList : `list` - List of data references to science exposures to be background-matched; + calExpList : `list` of `lsst.afw.image.Exposure` + List of science exposures to be background-matched; all exposures must exist. - expDatasetType : `str` - Dataset type of exposures, e.g. 'goodSeeingCoadd_tempExp'. - imageScalerList : `list`, optional - List of image scalers (coaddUtils.ImageScaler); - if None then the images are not scaled. - refExpDataRef : `Unknown`, optional - Data reference for the reference exposure. - If None, then this task selects the best exposures from expRefList. - If not None then must be one of the exposures in expRefList. - refImageScaler : `Unknown`, optional - Image scaler for reference image; - ignored if refExpDataRef is None, else scaling is not performed if None. + refCalExp : `lsst.afw.image.Exposure`, optional + Reference science exposure. + If None, then this task selects the best exposures from calExpList. + If not None then must be one of the exposures in calExpList. Returns ------- @@ -200,53 +234,46 @@ def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=No RuntimeError Raised if an exposure does not exist on disk. """ - numExp = len(expRefList) + numExp = len(calExpList) if numExp < 1: raise pipeBase.TaskError("No exposures to match") - if expDatasetType is None: - raise pipeBase.TaskError("Must specify expDatasetType") - - if imageScalerList is None: - self.log.info("imageScalerList is None; no scaling will be performed") - imageScalerList = [None] * numExp - - if len(expRefList) != len(imageScalerList): - raise RuntimeError("len(expRefList) = %s != %s = len(imageScalerList)" % - (len(expRefList), len(imageScalerList))) + # Removed several raise statements that no longer apply + breakpoint() + # The selectRefExposure() method now only requires one input + # refInd is the index in calExpList corresponding to refCalExp refInd = None - if refExpDataRef is None: - # select the best reference exposure from expRefList - refInd = self.selectRefExposure( - expRefList=expRefList, - imageScalerList=imageScalerList, - expDatasetType=expDatasetType, - ) - refExpDataRef = expRefList[refInd] - refImageScaler = imageScalerList[refInd] + #if refCalExp is None: + # select the best reference exposure from calExpList + refInd = self.selectRefExposure( + calExpList=calExpList, + ) + refCalExp = calExpList[refInd] # refIndSet is the index of all exposures in expDataList that match the reference. # It is used to avoid background-matching an exposure to itself. It is a list # because it is possible (though unlikely) that expDataList will contain duplicates. - expKeyList = refExpDataRef.butlerSubset.butler.getKeys(expDatasetType) - refMatcher = DataRefMatcher(refExpDataRef.butlerSubset.butler, expDatasetType) - refIndSet = set(refMatcher.matchList(ref0=refExpDataRef, refList=expRefList)) + refId = refCalExp.getInfo().id # Must be a better ID to use + expIds = numpy.array([exp.getInfo().id for exp in calExpList]) + refIndSet = list(numpy.where(expIds == refId)[0]) if refInd is not None and refInd not in refIndSet: - raise RuntimeError("Internal error: selected reference %s not found in expRefList") + raise RuntimeError("Internal error: selected reference %s not found in calExpList") - refExposure = refExpDataRef.get() - if refImageScaler is not None: - refMI = refExposure.getMaskedImage() - refImageScaler.scaleMaskedImage(refMI) + # Images must be scaled to a common ZP + # Converting everything to nJy to accomplish this + refZp = refCalExp.getPhotoCalib().instFluxToNanojansky(1) + refCalExp.image.array *= refZp - debugIdKeyList = tuple(set(expKeyList) - set(['tract', 'patch'])) + # Will determine later what this is supposed to be + # All unknown quantities being logged currently commented out + # debugIdKeyList = tuple(set(expKeyList) - set(['tract', 'patch'])) self.log.info("Matching %d Exposures", numExp) backgroundInfoList = [] - for ind, (toMatchRef, imageScaler) in enumerate(zip(expRefList, imageScalerList)): + for ind, exp in enumerate(calExpList): if ind in refIndSet: backgroundInfoStruct = pipeBase.Struct( isReference=True, @@ -256,21 +283,23 @@ def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=No diffImVar=None, ) else: - self.log.info("Matching background of %s to %s", toMatchRef.dataId, refExpDataRef.dataId) + # For L245: find the Gen3 equivalent to these IDs? + # self.log.info("Matching background of %s to %s", toMatchRef.dataId, refExpDataRef.dataId) try: - toMatchExposure = toMatchRef.get() - if imageScaler is not None: - toMatchMI = toMatchExposure.getMaskedImage() - imageScaler.scaleMaskedImage(toMatchMI) + toMatchExposure = exp.get() + toMatchMI = toMatchExposure.getMaskedImage() + fluxZp = toMatchExposure.getPhotoCalib().instFluxToNanojansky(1) + toMatchMI.image.array *= fluxZp # Left off here # store a string specifying the visit to label debug plot - self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList]) + # self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList]) + backgroundInfoStruct = self.matchBackgrounds( - refExposure=refExposure, + refExposure=refCalExp, sciExposure=toMatchExposure, ) backgroundInfoStruct.isReference = False except Exception as e: - self.log.warning("Failed to fit background %s: %s", toMatchRef.dataId, e) + # self.log.warning("Failed to fit background %s: %s", toMatchRef.dataId, e) backgroundInfoStruct = pipeBase.Struct( isReference=False, backgroundModel=None, @@ -285,7 +314,7 @@ def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=No backgroundInfoList=backgroundInfoList) @timeMethod - def selectRefExposure(self, expRefList, imageScalerList, expDatasetType): + def selectRefExposure(self, calExpList): """Find best exposure to use as the reference exposure. Calculate an appropriate reference exposure by minimizing a cost function that penalizes @@ -296,15 +325,9 @@ def selectRefExposure(self, expRefList, imageScalerList, expDatasetType): Parameters ---------- - expRefList : `list` - List of data references to exposures. - Retrieves dataset type specified by expDatasetType. + calExpList : `list` of `lsst.afw.image.Exposure` + List of science exposures. If an exposure is not found, it is skipped with a warning. - imageScalerList : `list` - List of image scalers (coaddUtils.ImageScaler); - must be the same length as expRefList. - expDatasetType : `str` - Dataset type of exposure: e.g. 'goodSeeingCoadd_tempExp'. Returns ------- @@ -314,29 +337,26 @@ def selectRefExposure(self, expRefList, imageScalerList, expDatasetType): Raises ------ RuntimeError - Raised if none of the exposures in expRefList are found. + Raised if none of the exposures in calExpList are found. """ self.log.info("Calculating best reference visit") varList = [] meanBkgdLevelList = [] coverageList = [] - if len(expRefList) != len(imageScalerList): - raise RuntimeError("len(expRefList) = %s != %s = len(imageScalerList)" % - (len(expRefList), len(imageScalerList))) - - for expRef, imageScaler in zip(expRefList, imageScalerList): - exposure = expRef.get() + for exp in calExpList: + exposure = exp.get() maskedImage = exposure.getMaskedImage() - if imageScaler is not None: - try: - imageScaler.scaleMaskedImage(maskedImage) - except Exception: - # need to put a place holder in Arr - varList.append(numpy.nan) - meanBkgdLevelList.append(numpy.nan) - coverageList.append(numpy.nan) - continue + # Convert images to nJy before doing statistics + try: + fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) + exposure.image.array *= fluxZp + except Exception: + # need to put a place holder in Arr + varList.append(numpy.nan) + meanBkgdLevelList.append(numpy.nan) + coverageList.append(numpy.nan) + continue statObjIm = afwMath.makeStatistics(maskedImage.getImage(), maskedImage.getMask(), afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, self.sctrl) meanVar, meanVarErr = statObjIm.getResult(afwMath.VARIANCE) @@ -347,7 +367,7 @@ def selectRefExposure(self, expRefList, imageScalerList, expDatasetType): coverageList.append(npoints) if not coverageList: raise pipeBase.TaskError( - "None of the candidate %s exist; cannot select best reference exposure" % (expDatasetType,)) + "None of the candidates exist; cannot select best reference exposure") # Normalize metrics to range from 0 to 1 varArr = numpy.array(varList)/numpy.nanmax(varList) From 5ed37bf11dc0c4b4b6da23b295c2d66845805445 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Mon, 8 Jul 2024 06:49:48 -0700 Subject: [PATCH 02/13] Change input data type to deepCoadd_psfMatchedWarp Code now runs without complaint through self.matchBackgrounds. Also added a self._fluxScale method to replace repeat code blocks. Will decide later if scaling to nJy is the best way to do this. --- python/lsst/pipe/tasks/matchBackgrounds.py | 192 ++++++++++++--------- 1 file changed, 108 insertions(+), 84 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 405339777..2a78aa49e 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -21,7 +21,7 @@ __all__ = ["MatchBackgroundsConfig", "MatchBackgroundsTask"] -import numpy +import numpy as np import lsst.afw.image as afwImage import lsst.afw.math as afwMath import lsst.geom as geom @@ -37,28 +37,34 @@ class MatchBackgroundsConnections(PipelineTaskConnections, - defaultTemplates={}, - dimensions={"instrument",}): - - calExpList = cT.Input( - doc="Input exposures to process", - dimensions=("visit", "detector",), + dimensions=("tract", "patch", "band", "skymap"), + defaultTemplates={"inputCoaddName": "deep", + "outputCoaddName": "deep", + "warpType": "direct", + "warpTypeSuffix": ""}): + + psfMatchedWarps = pipeBase.connectionTypes.Input( + doc=("PSF-matched warps to be subtracted from the reference warp"), + name="{inputCoaddName}Coadd_psfMatchedWarp", storageClass="ExposureF", - name="calexp", - multiple=True, # Probably this will be needed for the List format - deferLoad=True, # All image calls currently use .get(), so... + dimensions=("tract", "patch", "skymap", "visit"), + deferLoad=True, + multiple=True, ) - # refCalExp = cT.Input( + # The reference exposure needs to be another psfMatchedWarp + # Let the software choose it automatically for now + # refMatchedWarp = cT.Input( # doc="Reference exposure", # dimensions=("visit", "detector", "band"), # storageClass="ExposureF", # name="calexp", # ) + # This needs to be the models of each differential BG in warped coords backgroundInfoList = cT.Output( doc="List of differential backgrounds, w/goodness of fit params", + name="calexpBackground_diff", # This needs to change dimensions=("visit", "detector",), storageClass="Background", - name="calexpBackground_diff", multiple=True, ) @@ -177,12 +183,12 @@ class MatchBackgroundsTask(pipeBase.PipelineTask): ConfigClass = MatchBackgroundsConfig _DefaultName = "matchBackgrounds" - # def __init__(self, *args, **kwargs): - # pipeBase.Task.__init__(self, *args, **kwargs) + def __init__(self, *args, **kwargs): + super().__init__(**kwargs) - # self.sctrl = afwMath.StatisticsControl() - # self.sctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)) - # self.sctrl.setNanSafe(True) + self.sctrl = afwMath.StatisticsControl() + self.sctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)) + self.sctrl.setNanSafe(True) def runQuantum(self, butlerQC, inputRefs, outputRefs): ''' @@ -193,20 +199,20 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): butlerQC.put(outputs, outputRefs) @timeMethod - def run(self, calExpList): + def run(self, psfMatchedWarps): """Match the backgrounds of a list of science exposures to a reference science exposure. - Choose a refCalExp automatically if none supplied. + Choose a refMatchedWarp automatically if none supplied. Parameters ---------- - calExpList : `list` of `lsst.afw.image.Exposure` - List of science exposures to be background-matched; + psfMatchedWarps : `list` of `lsst.afw.image.Exposure` + List of warped science exposures to be background-matched; all exposures must exist. - refCalExp : `lsst.afw.image.Exposure`, optional + refMatchedWarp : `lsst.afw.image.Exposure`, optional Reference science exposure. - If None, then this task selects the best exposures from calExpList. - If not None then must be one of the exposures in calExpList. + If None, then this task selects the best exposures from psfMatchedWarps. + If not None then must be one of the exposures in psfMatchedWarps. Returns ------- @@ -214,7 +220,7 @@ def run(self, calExpList): Results as a struct with attributes: ``backgroundInfoList`` - A `list` of `pipeBase.Struct`, one per exposure in expRefList, + A `list` of `pipeBase.Struct`, one per exposure in psfMatchedWarps\s, each of which contains these fields: - ``isReference``: This is the reference exposure (only one returned Struct will contain True for this @@ -234,37 +240,39 @@ def run(self, calExpList): RuntimeError Raised if an exposure does not exist on disk. """ - numExp = len(calExpList) + numExp = len(psfMatchedWarps) if numExp < 1: raise pipeBase.TaskError("No exposures to match") - # Removed several raise statements that no longer apply - breakpoint() - # The selectRefExposure() method now only requires one input - # refInd is the index in calExpList corresponding to refCalExp - refInd = None - #if refCalExp is None: - # select the best reference exposure from calExpList + # refInd is the index in psfMatchedWarps corresponding to refMatchedWarp + # refInd = None + # if refMatchedWarp is None: + # select the best reference exposure from psfMatchedWarps + # note: just selecting from the set for testing purposes refInd = self.selectRefExposure( - calExpList=calExpList, + psfMatchedWarps=psfMatchedWarps, ) - refCalExp = calExpList[refInd] + refMatchedWarp = psfMatchedWarps[refInd] - # refIndSet is the index of all exposures in expDataList that match the reference. + # refIndSet is the index of all exposures in psfMatchedWarps that match the reference. # It is used to avoid background-matching an exposure to itself. It is a list - # because it is possible (though unlikely) that expDataList will contain duplicates. - refId = refCalExp.getInfo().id # Must be a better ID to use - expIds = numpy.array([exp.getInfo().id for exp in calExpList]) - refIndSet = list(numpy.where(expIds == refId)[0]) + # because it is possible (though unlikely) that psfMatchedWarps will contain duplicates. + refId = refMatchedWarp.dataId + expIds = [exp.dataId for exp in psfMatchedWarps] + try: + refIndSet = [expIds.index(i) for i in [refId]] + except ValueError: + raise RuntimeError("Internal error: selected reference %s not found in psfMatchedWarps") - if refInd is not None and refInd not in refIndSet: - raise RuntimeError("Internal error: selected reference %s not found in calExpList") + # This is the original check, probably no longer needed. + # if refInd is not None and refInd not in refIndSet: + # raise RuntimeError("Internal error: selected reference %s not found in psfMatchedWarps") # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this - refZp = refCalExp.getPhotoCalib().instFluxToNanojansky(1) - refCalExp.image.array *= refZp + refExposure = refMatchedWarp.get() + refMI = self._fluxScale(refExposure) # Also modifies refExposure # Will determine later what this is supposed to be # All unknown quantities being logged currently commented out @@ -273,7 +281,7 @@ def run(self, calExpList): self.log.info("Matching %d Exposures", numExp) backgroundInfoList = [] - for ind, exp in enumerate(calExpList): + for ind, exp in enumerate(psfMatchedWarps): if ind in refIndSet: backgroundInfoStruct = pipeBase.Struct( isReference=True, @@ -283,18 +291,17 @@ def run(self, calExpList): diffImVar=None, ) else: - # For L245: find the Gen3 equivalent to these IDs? - # self.log.info("Matching background of %s to %s", toMatchRef.dataId, refExpDataRef.dataId) + self.log.info("Matching background of %s to %s", exp.dataId, refMatchedWarp.dataId) + toMatchExposure = exp.get() try: - toMatchExposure = exp.get() - toMatchMI = toMatchExposure.getMaskedImage() - fluxZp = toMatchExposure.getPhotoCalib().instFluxToNanojansky(1) - toMatchMI.image.array *= fluxZp # Left off here + # Why use exposure and not maskedImage? + toMatchMI = self._fluxScale(toMatchExposure) + # store a string specifying the visit to label debug plot # self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList]) backgroundInfoStruct = self.matchBackgrounds( - refExposure=refCalExp, + refExposure=refExposure, sciExposure=toMatchExposure, ) backgroundInfoStruct.isReference = False @@ -314,7 +321,7 @@ def run(self, calExpList): backgroundInfoList=backgroundInfoList) @timeMethod - def selectRefExposure(self, calExpList): + def selectRefExposure(self, psfMatchedWarps): """Find best exposure to use as the reference exposure. Calculate an appropriate reference exposure by minimizing a cost function that penalizes @@ -325,7 +332,7 @@ def selectRefExposure(self, calExpList): Parameters ---------- - calExpList : `list` of `lsst.afw.image.Exposure` + psfMatchedWarps : `list` of `lsst.afw.image.Exposure` List of science exposures. If an exposure is not found, it is skipped with a warning. @@ -337,25 +344,23 @@ def selectRefExposure(self, calExpList): Raises ------ RuntimeError - Raised if none of the exposures in calExpList are found. + Raised if none of the exposures in psfMatchedWarps are found. """ self.log.info("Calculating best reference visit") varList = [] meanBkgdLevelList = [] coverageList = [] - for exp in calExpList: + for exp in psfMatchedWarps: exposure = exp.get() - maskedImage = exposure.getMaskedImage() # Convert images to nJy before doing statistics try: - fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) - exposure.image.array *= fluxZp + maskedImage = self._fluxScale(exposure) except Exception: # need to put a place holder in Arr - varList.append(numpy.nan) - meanBkgdLevelList.append(numpy.nan) - coverageList.append(numpy.nan) + varList.append(np.nan) + meanBkgdLevelList.append(np.nan) + coverageList.append(np.nan) continue statObjIm = afwMath.makeStatistics(maskedImage.getImage(), maskedImage.getMask(), afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, self.sctrl) @@ -370,14 +375,14 @@ def selectRefExposure(self, calExpList): "None of the candidates exist; cannot select best reference exposure") # Normalize metrics to range from 0 to 1 - varArr = numpy.array(varList)/numpy.nanmax(varList) - meanBkgdLevelArr = numpy.array(meanBkgdLevelList)/numpy.nanmax(meanBkgdLevelList) - coverageArr = numpy.nanmin(coverageList)/numpy.array(coverageList) + varArr = np.array(varList)/np.nanmax(varList) + meanBkgdLevelArr = np.array(meanBkgdLevelList)/np.nanmax(meanBkgdLevelList) + coverageArr = np.nanmin(coverageList)/np.array(coverageList) costFunctionArr = self.config.bestRefWeightVariance * varArr costFunctionArr += self.config.bestRefWeightLevel * meanBkgdLevelArr costFunctionArr += self.config.bestRefWeightCoverage * coverageArr - return numpy.nanargmin(costFunctionArr) + return np.nanargmin(costFunctionArr) @timeMethod def matchBackgrounds(self, refExposure, sciExposure): @@ -489,7 +494,7 @@ def matchBackgrounds(self, refExposure, sciExposure): self.log.warning("Decreasing binsize to %d", newBinSize) # If there is no variance in any image pixels, do not weight bins by inverse variance - isUniformImageDiff = not numpy.any(bgdZ > self.config.gridStdevEpsilon) + isUniformImageDiff = not np.any(bgdZ > self.config.gridStdevEpsilon) weightByInverseVariance = False if isUniformImageDiff else self.config.approxWeighting # Add offset to sciExposure @@ -514,11 +519,11 @@ def matchBackgrounds(self, refExposure, sciExposure): rms = 0.0 X, Y, Z, dZ = self._gridImage(diffMI, self.config.binSize, statsFlag) x0, y0 = diffMI.getXY0() - modelValueArr = numpy.empty(len(Z)) + modelValueArr = np.empty(len(Z)) for i in range(len(X)): modelValueArr[i] = bkgdImage[int(X[i]-x0), int(Y[i]-y0), afwImage.LOCAL] resids = Z - modelValueArr - rms = numpy.sqrt(numpy.mean(resids[~numpy.isnan(resids)]**2)) + rms = np.sqrt(np.mean(resids[~np.isnan(resids)]**2)) if lsstDebug.Info(__name__).savefits: sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciMatchedExposure.fits') @@ -546,6 +551,25 @@ def matchBackgrounds(self, refExposure, sciExposure): matchedMSE=mse, diffImVar=meanVar) + def _fluxScale(self, exposure): + """Scales image to nJy flux using photometric calibration. + + Parameters + ---------- + exposure: `` + Exposure to scale. + + Returns + ------- + maskedImage: `` + Flux-scaled masked exposure. + """ + maskedImage = exposure.getMaskedImage() + fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) + exposure.image.array *= fluxZp + + return maskedImage + def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): """Generate a plot showing the background fit and residuals. @@ -555,17 +579,17 @@ def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): Parameters ---------- - X : `numpy.ndarray`, (N,) + X : `np.ndarray`, (N,) Array of x positions. - Y : `numpy.ndarray`, (N,) + Y : `np.ndarray`, (N,) Array of y positions. - Z : `numpy.ndarray` + Z : `np.ndarray` Array of the grid values that were interpolated. - dZ : `numpy.ndarray`, (len(Z),) + dZ : `np.ndarray`, (len(Z),) Array of the error on the grid values. modelImage : `Unknown` Image of the model of the fit. - model : `numpy.ndarray`, (len(Z),) + model : `np.ndarray`, (len(Z),) Array of len(Z) containing the grid values predicted by the model. resids : `Unknown` Z - model. @@ -580,13 +604,13 @@ def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): if len(Z) == 0: self.log.warning("No grid. Skipping plot generation.") else: - max, min = numpy.max(Z), numpy.min(Z) + max, min = np.max(Z), np.min(Z) norm = matplotlib.colors.normalize(vmax=max, vmin=min) - maxdiff = numpy.max(numpy.abs(resids)) + maxdiff = np.max(np.abs(resids)) diffnorm = matplotlib.colors.normalize(vmax=maxdiff, vmin=-maxdiff) - rms = numpy.sqrt(numpy.mean(resids**2)) + rms = np.sqrt(np.mean(resids**2)) fig = plt.figure(1, (8, 6)) - meanDz = numpy.mean(dZ) + meanDz = np.mean(dZ) grid = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.1, share_all=True, label_mode="L", cbar_mode="each", cbar_size="7%", cbar_pad="2%", cbar_location="top") @@ -613,10 +637,10 @@ def _gridImage(self, maskedImage, binsize, statsFlag): """Private method to grid an image for debugging.""" width, height = maskedImage.getDimensions() x0, y0 = maskedImage.getXY0() - xedges = numpy.arange(0, width, binsize) - yedges = numpy.arange(0, height, binsize) - xedges = numpy.hstack((xedges, width)) # add final edge - yedges = numpy.hstack((yedges, height)) # add final edge + xedges = np.arange(0, width, binsize) + yedges = np.arange(0, height, binsize) + xedges = np.hstack((xedges, width)) # add final edge + yedges = np.hstack((yedges, height)) # add final edge # Use lists/append to protect against the case where # a bin has no valid pixels and should not be included in the fit @@ -641,11 +665,11 @@ def _gridImage(self, maskedImage, binsize, statsFlag): stdev = self.config.gridStdevEpsilon bgX.append(0.5 * (x0 + xmin + x0 + xmax)) bgY.append(0.5 * (y0 + ymin + y0 + ymax)) - bgdZ.append(stdev/numpy.sqrt(npoints)) + bgdZ.append(stdev/np.sqrt(npoints)) est, _ = stats.getResult(statsFlag) bgZ.append(est) - return numpy.array(bgX), numpy.array(bgY), numpy.array(bgZ), numpy.array(bgdZ) + return np.array(bgX), np.array(bgY), np.array(bgZ), np.array(bgdZ) class DataRefMatcher: From 73f93725fc83177f59836d5b1f0dbd51373b96d3 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Wed, 10 Jul 2024 01:10:25 -0700 Subject: [PATCH 03/13] Change output from Struct to BackgroundList Code is now functional, in that it accepts images and returns difference image background models as "psfMatchedWarpBackground_diff" (name likely to be altered later). Uses a fit to a blank image for that corresponding to the reference image. --- python/lsst/pipe/tasks/matchBackgrounds.py | 58 ++++++++++++++-------- 1 file changed, 37 insertions(+), 21 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 2a78aa49e..a9bfc84a2 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -62,8 +62,8 @@ class MatchBackgroundsConnections(PipelineTaskConnections, # This needs to be the models of each differential BG in warped coords backgroundInfoList = cT.Output( doc="List of differential backgrounds, w/goodness of fit params", - name="calexpBackground_diff", # This needs to change - dimensions=("visit", "detector",), + name="psfMatchedWarpBackground_diff", # This needs to change + dimensions=("tract", "patch", "skymap", "visit"), storageClass="Background", multiple=True, ) @@ -280,16 +280,35 @@ def run(self, psfMatchedWarps): self.log.info("Matching %d Exposures", numExp) + # Creating a null BackgroundList object by fitting a blank image + statsFlag = getattr(afwMath, self.config.gridStatistic) + self.sctrl.setNumSigmaClip(self.config.numSigmaClip) + self.sctrl.setNumIter(self.config.numIter) + + # TODO: refactor below to construct blank bg model + im = refExposure.getMaskedImage() + blankIm = im.Factory(im, True) # Don't do this + blankIm.image.array *= 0 + + width = blankIm.getWidth() + height = blankIm.getHeight() + nx = width // self.config.binSize + if width % self.config.binSize != 0: + nx += 1 + ny = height // self.config.binSize + if height % self.config.binSize != 0: + ny += 1 + + bctrl = afwMath.BackgroundControl(nx, ny, self.sctrl, statsFlag) + bctrl.setUndersampleStyle(self.config.undersampleStyle) + + bkgd = afwMath.makeBackground(blankIm, bctrl) + + backgroundInfoList = [] for ind, exp in enumerate(psfMatchedWarps): if ind in refIndSet: - backgroundInfoStruct = pipeBase.Struct( - isReference=True, - backgroundModel=None, - fitRMS=0.0, - matchedMSE=None, - diffImVar=None, - ) + backgroundInfoStruct = afwMath.BackgroundList(bkgd,) else: self.log.info("Matching background of %s to %s", exp.dataId, refMatchedWarp.dataId) toMatchExposure = exp.get() @@ -307,13 +326,7 @@ def run(self, psfMatchedWarps): backgroundInfoStruct.isReference = False except Exception as e: # self.log.warning("Failed to fit background %s: %s", toMatchRef.dataId, e) - backgroundInfoStruct = pipeBase.Struct( - isReference=False, - backgroundModel=None, - fitRMS=None, - matchedMSE=None, - diffImVar=None, - ) + backgroundInfoStruct = afwMath.BackgroundList(bkgd,) backgroundInfoList.append(backgroundInfoStruct) @@ -545,11 +558,14 @@ def matchBackgrounds(self, refExposure, sciExposure): outBkgd = approx if self.config.usePolynomial else bkgd - return pipeBase.Struct( - backgroundModel=outBkgd, - fitRMS=rms, - matchedMSE=mse, - diffImVar=meanVar) + # Type `Background` can't use a struct. Should fitRMS &c. be added to + # a log instead of output? + # return pipeBase.Struct( + # backgroundModel=afwMath.BackgroundList(outBkgd), + # fitRMS=rms, + # matchedMSE=mse, + # diffImVar=meanVar) + return afwMath.BackgroundList(outBkgd,) def _fluxScale(self, exposure): """Scales image to nJy flux using photometric calibration. From c59c4a90af89ee5ca1727b239d8b7551623aa049 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Tue, 16 Jul 2024 02:58:46 -0700 Subject: [PATCH 04/13] Apply spline parameters to BackgroundList Difference background models are now formatted properly, to allow for image creation from the spline parameters. Also did some adjustments to documentation for Flake8 formatting. --- python/lsst/pipe/tasks/matchBackgrounds.py | 332 ++++++++++++--------- 1 file changed, 184 insertions(+), 148 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index a9bfc84a2..91db387af 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -21,27 +21,30 @@ __all__ = ["MatchBackgroundsConfig", "MatchBackgroundsTask"] -import numpy as np +import pdb + import lsst.afw.image as afwImage import lsst.afw.math as afwMath import lsst.geom as geom import lsst.pex.config as pexConfig import lsst.pipe.base as pipeBase +import lsst.pipe.base.connectionTypes as cT import lsstDebug +import numpy as np +from lsst.pipe.base import PipelineTaskConfig, PipelineTaskConnections from lsst.utils.timer import timeMethod -from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections -import lsst.pipe.base.connectionTypes as cT - -import pdb - -class MatchBackgroundsConnections(PipelineTaskConnections, - dimensions=("tract", "patch", "band", "skymap"), - defaultTemplates={"inputCoaddName": "deep", - "outputCoaddName": "deep", - "warpType": "direct", - "warpTypeSuffix": ""}): +class MatchBackgroundsConnections( + PipelineTaskConnections, + dimensions=("tract", "patch", "band", "skymap"), + defaultTemplates={ + "inputCoaddName": "deep", + "outputCoaddName": "deep", + "warpType": "direct", + "warpTypeSuffix": "", + }, +): psfMatchedWarps = pipeBase.connectionTypes.Input( doc=("PSF-matched warps to be subtracted from the reference warp"), @@ -70,34 +73,30 @@ class MatchBackgroundsConnections(PipelineTaskConnections, # class MatchBackgroundsConfig(pexConfig.Config): -class MatchBackgroundsConfig(PipelineTaskConfig, - pipelineConnections=MatchBackgroundsConnections): +class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgroundsConnections): usePolynomial = pexConfig.Field( dtype=bool, doc="Fit background difference with Chebychev polynomial interpolation " "(using afw.math.Approximate)? If False, fit with spline interpolation using afw.math.Background", - default=False + default=False, ) order = pexConfig.Field( dtype=int, doc="Order of Chebyshev polynomial background model. Ignored if usePolynomial False", - default=8 + default=8, ) badMaskPlanes = pexConfig.ListField( doc="Names of mask planes to ignore while estimating the background", - dtype=str, default=["NO_DATA", "DETECTED", "DETECTED_NEGATIVE", "SAT", "BAD", "INTRP", "CR"], + dtype=str, + default=["NO_DATA", "DETECTED", "DETECTED_NEGATIVE", "SAT", "BAD", "INTRP", "CR"], itemCheck=lambda x: x in afwImage.Mask().getMaskPlaneDict(), ) gridStatistic = pexConfig.ChoiceField( dtype=str, doc="Type of statistic to estimate pixel value for the grid points", default="MEAN", - allowed={ - "MEAN": "mean", - "MEDIAN": "median", - "MEANCLIP": "clipped mean" - } + allowed={"MEAN": "mean", "MEDIAN": "median", "MEANCLIP": "clipped mean"}, ) undersampleStyle = pexConfig.ChoiceField( doc="Behaviour if there are too few points in grid for requested interpolation style. " @@ -108,12 +107,10 @@ class MatchBackgroundsConfig(PipelineTaskConfig, "THROW_EXCEPTION": "throw an exception if there are too few points", "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.", "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.", - } + }, ) binSize = pexConfig.Field( - doc="Bin size for gridding the difference image and fitting a spatial model", - dtype=int, - default=256 + doc="Bin size for gridding the difference image and fitting a spatial model", dtype=int, default=256 ) interpStyle = pexConfig.ChoiceField( dtype=str, @@ -126,17 +123,15 @@ class MatchBackgroundsConfig(PipelineTaskConfig, "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints", "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers", "NONE": "No background estimation is to be attempted", - } + }, ) numSigmaClip = pexConfig.Field( - dtype=int, - doc="Sigma for outlier rejection; ignored if gridStatistic != 'MEANCLIP'.", - default=3 + dtype=int, doc="Sigma for outlier rejection; ignored if gridStatistic != 'MEANCLIP'.", default=3 ) numIter = pexConfig.Field( dtype=int, doc="Number of iterations of outlier rejection; ignored if gridStatistic != 'MEANCLIP'.", - default=2 + default=2, ) bestRefWeightCoverage = pexConfig.RangeField( dtype=float, @@ -144,14 +139,16 @@ class MatchBackgroundsConfig(PipelineTaskConfig, "when calculating best reference exposure. Higher weight prefers exposures with high coverage." "Ignored when reference visit is supplied", default=0.4, - min=0., max=1. + min=0.0, + max=1.0, ) bestRefWeightVariance = pexConfig.RangeField( dtype=float, doc="Weight given to image variance when calculating best reference exposure. " "Higher weight prefers exposures with low image variance. Ignored when reference visit is supplied", default=0.4, - min=0., max=1. + min=0.0, + max=1.0, ) bestRefWeightLevel = pexConfig.RangeField( dtype=float, @@ -159,14 +156,17 @@ class MatchBackgroundsConfig(PipelineTaskConfig, "Higher weight prefers exposures with low mean background level. " "Ignored when reference visit is supplied.", default=0.2, - min=0., max=1. + min=0.0, + max=1.0, ) approxWeighting = pexConfig.Field( dtype=bool, - doc=("Use inverse-variance weighting when approximating background offset model? " - "This will fail when the background offset is constant " - "(this is usually only the case in testing with artificial images)." - "(usePolynomial=True)"), + doc=( + "Use inverse-variance weighting when approximating background offset model? " + "This will fail when the background offset is constant " + "(this is usually only the case in testing with artificial images)." + "(usePolynomial=True)" + ), default=True, ) gridStdevEpsilon = pexConfig.RangeField( @@ -175,7 +175,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, "If all bins have a standard deviation below this value, the background offset model " "is approximated without inverse-variance weighting. (usePolynomial=True)", default=1e-8, - min=0. + min=0.0, ) @@ -191,18 +191,17 @@ def __init__(self, *args, **kwargs): self.sctrl.setNanSafe(True) def runQuantum(self, butlerQC, inputRefs, outputRefs): - ''' - Boilerplate for now, until bug-testing commences - ''' inputs = butlerQC.get(inputRefs) outputs = self.run(**inputs) butlerQC.put(outputs, outputRefs) @timeMethod def run(self, psfMatchedWarps): - """Match the backgrounds of a list of science exposures to a reference science exposure. + """Match the backgrounds of a list of science exposures to a reference + science exposure. Choose a refMatchedWarp automatically if none supplied. + (above is legacy: right now this only does it automatically) Parameters ---------- @@ -211,29 +210,15 @@ def run(self, psfMatchedWarps): all exposures must exist. refMatchedWarp : `lsst.afw.image.Exposure`, optional Reference science exposure. - If None, then this task selects the best exposures from psfMatchedWarps. + If None, then this task selects the best exposures from + psfMatchedWarps. If not None then must be one of the exposures in psfMatchedWarps. Returns ------- - result : `lsst.pipe.base.Struct` - Results as a struct with attributes: - - ``backgroundInfoList`` - A `list` of `pipeBase.Struct`, one per exposure in psfMatchedWarps\s, - each of which contains these fields: - - ``isReference``: This is the reference exposure (only one - returned Struct will contain True for this - value, unless the ref exposure is listed multiple times). - - ``backgroundModel``: Differential background model - (afw.Math.Background or afw.Math.Approximate). - Add this to the science exposure to match the reference exposure. - - ``fitRMS``: The RMS of the fit. This is the sqrt(mean(residuals**2)). - - ``matchedMSE``: The MSE of the reference and matched images: - mean((refImage - matchedSciImage)**2); - should be comparable to difference image's mean variance. - - ``diffImVar``: The mean variance of the difference image. - All fields except isReference will be None if isReference True or the fit failed. + result : `lsst.afw.math._backgroundList.BackgroundList` + Differential background model + Add this to the science exposure to match the reference exposure. Raises ------ @@ -244,8 +229,7 @@ def run(self, psfMatchedWarps): if numExp < 1: raise pipeBase.TaskError("No exposures to match") - # The selectRefExposure() method now only requires one input - # refInd is the index in psfMatchedWarps corresponding to refMatchedWarp + # refInd: index in psfMatchedWarps corresponding to refMatchedWarp # refInd = None # if refMatchedWarp is None: # select the best reference exposure from psfMatchedWarps @@ -255,9 +239,10 @@ def run(self, psfMatchedWarps): ) refMatchedWarp = psfMatchedWarps[refInd] - # refIndSet is the index of all exposures in psfMatchedWarps that match the reference. - # It is used to avoid background-matching an exposure to itself. It is a list - # because it is possible (though unlikely) that psfMatchedWarps will contain duplicates. + # refIndSet: indices of exposures in psfMatchedWarps matching ref. + # It is used to avoid background-matching an exposure to itself. + # It is a list because it is possible (though unlikely) that + # psfMatchedWarps will contain duplicates. refId = refMatchedWarp.dataId expIds = [exp.dataId for exp in psfMatchedWarps] try: @@ -265,17 +250,12 @@ def run(self, psfMatchedWarps): except ValueError: raise RuntimeError("Internal error: selected reference %s not found in psfMatchedWarps") - # This is the original check, probably no longer needed. - # if refInd is not None and refInd not in refIndSet: - # raise RuntimeError("Internal error: selected reference %s not found in psfMatchedWarps") - # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this refExposure = refMatchedWarp.get() refMI = self._fluxScale(refExposure) # Also modifies refExposure - # Will determine later what this is supposed to be - # All unknown quantities being logged currently commented out + # TODO: figure out what this was creating and why # debugIdKeyList = tuple(set(expKeyList) - set(['tract', 'patch'])) self.log.info("Matching %d Exposures", numExp) @@ -304,16 +284,26 @@ def run(self, psfMatchedWarps): bkgd = afwMath.makeBackground(blankIm, bctrl) - backgroundInfoList = [] for ind, exp in enumerate(psfMatchedWarps): if ind in refIndSet: - backgroundInfoStruct = afwMath.BackgroundList(bkgd,) + backgroundInfoStruct = afwMath.BackgroundList( + ( + bkgd, + afwMath.stringToInterpStyle(self.config.interpStyle), + afwMath.stringToUndersampleStyle(self.config.undersampleStyle), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) else: + # TODO: simplify this maybe, using only visit IDs? self.log.info("Matching background of %s to %s", exp.dataId, refMatchedWarp.dataId) toMatchExposure = exp.get() try: - # Why use exposure and not maskedImage? + # Seems to require exposure, not masked exposure. toMatchMI = self._fluxScale(toMatchExposure) # store a string specifying the visit to label debug plot @@ -325,20 +315,30 @@ def run(self, psfMatchedWarps): ) backgroundInfoStruct.isReference = False except Exception as e: - # self.log.warning("Failed to fit background %s: %s", toMatchRef.dataId, e) - backgroundInfoStruct = afwMath.BackgroundList(bkgd,) + self.log.warning("Failed to fit background %s: %s", exp.dataId, e) + backgroundInfoStruct = afwMath.BackgroundList( + ( + bkgd, + afwMath.stringToInterpStyle(self.config.interpStyle), + afwMath.stringToUndersampleStyle(self.config.undersampleStyle), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) backgroundInfoList.append(backgroundInfoStruct) - return pipeBase.Struct( - backgroundInfoList=backgroundInfoList) + return pipeBase.Struct(backgroundInfoList=backgroundInfoList) @timeMethod def selectRefExposure(self, psfMatchedWarps): """Find best exposure to use as the reference exposure. - Calculate an appropriate reference exposure by minimizing a cost function that penalizes - high variance, high background level, and low coverage. Use the following config parameters: + Calculate an appropriate reference exposure by minimizing a cost + function that penalizes high variance, high background level, and low + coverage. Use the following config parameters: - bestRefWeightCoverage - bestRefWeightVariance - bestRefWeightLevel @@ -375,8 +375,12 @@ def selectRefExposure(self, psfMatchedWarps): meanBkgdLevelList.append(np.nan) coverageList.append(np.nan) continue - statObjIm = afwMath.makeStatistics(maskedImage.getImage(), maskedImage.getMask(), - afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, self.sctrl) + statObjIm = afwMath.makeStatistics( + maskedImage.getImage(), + maskedImage.getMask(), + afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, + self.sctrl, + ) meanVar, meanVarErr = statObjIm.getResult(afwMath.VARIANCE) meanBkgdLevel, meanBkgdLevelErr = statObjIm.getResult(afwMath.MEAN) npoints, npointsErr = statObjIm.getResult(afwMath.NPOINT) @@ -384,13 +388,12 @@ def selectRefExposure(self, psfMatchedWarps): meanBkgdLevelList.append(meanBkgdLevel) coverageList.append(npoints) if not coverageList: - raise pipeBase.TaskError( - "None of the candidates exist; cannot select best reference exposure") + raise pipeBase.TaskError("None of the candidates exist; cannot select best reference exposure") # Normalize metrics to range from 0 to 1 - varArr = np.array(varList)/np.nanmax(varList) - meanBkgdLevelArr = np.array(meanBkgdLevelList)/np.nanmax(meanBkgdLevelList) - coverageArr = np.nanmin(coverageList)/np.array(coverageList) + varArr = np.array(varList) / np.nanmax(varList) + meanBkgdLevelArr = np.array(meanBkgdLevelList) / np.nanmax(meanBkgdLevelList) + coverageArr = np.nanmin(coverageList) / np.array(coverageList) costFunctionArr = self.config.bestRefWeightVariance * varArr costFunctionArr += self.config.bestRefWeightLevel * meanBkgdLevelArr @@ -399,16 +402,19 @@ def selectRefExposure(self, psfMatchedWarps): @timeMethod def matchBackgrounds(self, refExposure, sciExposure): - """Match science exposure's background level to that of reference exposure. + """Match science exposure's background level to that of reference + exposure. - Process creates a difference image of the reference exposure minus the science exposure, and then - generates an afw.math.Background object. It assumes (but does not require/check) that the mask plane - already has detections set. If detections have not been set/masked, sources will bias the - background estimation. + Process creates a difference image of the reference exposure minus the + science exposure, and then generates an afw.math.Background object. It + assumes (but does not require/check) that the mask plane already has + detections set. If detections have not been set/masked, sources will + bias the background estimation. - The 'background' of the difference image is smoothed by spline interpolation (by the Background class) - or by polynomial interpolation by the Approximate class. This model of difference image - is added to the science exposure in memory. + The 'background' of the difference image is smoothed by spline + interpolation (by the Background class) or by polynomial interpolation + by the Approximate class. This model of difference image is added to + the science exposure in memory. Fit diagnostics are also calculated and returned. @@ -428,24 +434,27 @@ def matchBackgrounds(self, refExposure, sciExposure): ``backgroundModel`` An afw.math.Approximate or an afw.math.Background. ``fitRMS`` - RMS of the fit. This is the sqrt(mean(residuals**2)), (`float`). + RMS of the fit = sqrt(mean(residuals**2)), (`float`). ``matchedMSE`` - The MSE of the reference and matched images: mean((refImage - matchedSciImage)**2); - should be comparable to difference image's mean variance (`float`). + The MSE of the reference and matched images: + mean((refImage - matchedSciImage)**2); + Should be comparable to difference image's mean variance + (`float`). ``diffImVar`` The mean variance of the difference image (`float`). """ if lsstDebug.Info(__name__).savefits: - refExposure.writeFits(lsstDebug.Info(__name__).figpath + 'refExposure.fits') - sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciExposure.fits') + refExposure.writeFits(lsstDebug.Info(__name__).figpath + "refExposure.fits") + sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciExposure.fits") # Check Configs for polynomials: if self.config.usePolynomial: x, y = sciExposure.getDimensions() shortSideLength = min(x, y) if shortSideLength < self.config.binSize: - raise ValueError("%d = config.binSize > shorter dimension = %d" % (self.config.binSize, - shortSideLength)) + raise ValueError( + "%d = config.binSize > shorter dimension = %d" % (self.config.binSize, shortSideLength) + ) npoints = shortSideLength // self.config.binSize if shortSideLength % self.config.binSize != 0: npoints += 1 @@ -454,12 +463,12 @@ def matchBackgrounds(self, refExposure, sciExposure): raise ValueError("%d = config.order > npoints - 1 = %d" % (self.config.order, npoints - 1)) # Check that exposures are same shape - if (sciExposure.getDimensions() != refExposure.getDimensions()): + if sciExposure.getDimensions() != refExposure.getDimensions(): wSci, hSci = sciExposure.getDimensions() wRef, hRef = refExposure.getDimensions() raise RuntimeError( - "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" % - (wSci, hSci, wRef, hRef)) + "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" % (wSci, hSci, wRef, hRef) + ) statsFlag = getattr(afwMath, self.config.gridStatistic) self.sctrl.setNumSigmaClip(self.config.numSigmaClip) @@ -500,7 +509,7 @@ def matchBackgrounds(self, refExposure, sciExposure): self.log.warning("Reducing order to %d", (minNumberGridPoints - 1)) order = minNumberGridPoints - 1 elif self.config.undersampleStyle == "INCREASE_NXNYSAMPLE": - newBinSize = (minNumberGridPoints*self.config.binSize) // (self.config.order + 1) + newBinSize = (minNumberGridPoints * self.config.binSize) // (self.config.order + 1) bctrl.setNxSample(newBinSize) bctrl.setNySample(newBinSize) bkgd = afwMath.makeBackground(diffMI, bctrl) # do over @@ -513,16 +522,18 @@ def matchBackgrounds(self, refExposure, sciExposure): # Add offset to sciExposure try: if self.config.usePolynomial: - actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV, - order, order, weightByInverseVariance) + actrl = afwMath.ApproximateControl( + afwMath.ApproximateControl.CHEBYSHEV, order, order, weightByInverseVariance + ) undersampleStyle = getattr(afwMath, self.config.undersampleStyle) approx = bkgd.getApproximate(actrl, undersampleStyle) bkgdImage = approx.getImage() else: bkgdImage = bkgd.getImageF(self.config.interpStyle, self.config.undersampleStyle) except Exception as e: - raise RuntimeError("Background/Approximation failed to interp image %s: %s" % ( - self.debugDataIdString, e)) + raise RuntimeError( + "Background/Approximation failed to interp image %s: %s" % (self.debugDataIdString, e) + ) sciMI = sciExposure.getMaskedImage() sciMI += bkgdImage @@ -534,22 +545,23 @@ def matchBackgrounds(self, refExposure, sciExposure): x0, y0 = diffMI.getXY0() modelValueArr = np.empty(len(Z)) for i in range(len(X)): - modelValueArr[i] = bkgdImage[int(X[i]-x0), int(Y[i]-y0), afwImage.LOCAL] + modelValueArr[i] = bkgdImage[int(X[i] - x0), int(Y[i] - y0), afwImage.LOCAL] resids = Z - modelValueArr - rms = np.sqrt(np.mean(resids[~np.isnan(resids)]**2)) + rms = np.sqrt(np.mean(resids[~np.isnan(resids)] ** 2)) if lsstDebug.Info(__name__).savefits: - sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciMatchedExposure.fits') + sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciMatchedExposure.fits") if lsstDebug.Info(__name__).savefig: bbox = geom.Box2D(refExposure.getMaskedImage().getBBox()) try: self._debugPlot(X, Y, Z, dZ, bkgdImage, bbox, modelValueArr, resids) except Exception as e: - self.log.warning('Debug plot not generated: %s', e) + self.log.warning("Debug plot not generated: %s", e) - meanVar = afwMath.makeStatistics(diffMI.getVariance(), diffMI.getMask(), - afwMath.MEANCLIP, self.sctrl).getValue() + meanVar = afwMath.makeStatistics( + diffMI.getVariance(), diffMI.getMask(), afwMath.MEANCLIP, self.sctrl + ).getValue() diffIm = diffMI.getImage() diffIm -= bkgdImage # diffMI should now have a mean ~ 0 @@ -558,14 +570,25 @@ def matchBackgrounds(self, refExposure, sciExposure): outBkgd = approx if self.config.usePolynomial else bkgd - # Type `Background` can't use a struct. Should fitRMS &c. be added to - # a log instead of output? - # return pipeBase.Struct( - # backgroundModel=afwMath.BackgroundList(outBkgd), - # fitRMS=rms, - # matchedMSE=mse, - # diffImVar=meanVar) - return afwMath.BackgroundList(outBkgd,) + self.log.info( + "Visit %d; difference BG fit RMS=%.1f cts, matched MSE=%.1f cts, mean variance=%.1f cts", + sciExposure.getInfo().getVisitInfo().id, + rms, + mse, + meanVar, + ) + # TODO: verify this is correct (borrowed from background.py) + return afwMath.BackgroundList( + ( + outBkgd, + afwMath.stringToInterpStyle(self.config.interpStyle), + afwMath.stringToUndersampleStyle(self.config.undersampleStyle), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) def _fluxScale(self, exposure): """Scales image to nJy flux using photometric calibration. @@ -610,9 +633,10 @@ def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): resids : `Unknown` Z - model. """ - import matplotlib.pyplot as plt import matplotlib.colors + import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import ImageGrid + zeroIm = afwImage.MaskedImageF(geom.Box2I(bbox)) zeroIm += modelImage x0, y0 = zeroIm.getXY0() @@ -627,24 +651,33 @@ def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): rms = np.sqrt(np.mean(resids**2)) fig = plt.figure(1, (8, 6)) meanDz = np.mean(dZ) - grid = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.1, - share_all=True, label_mode="L", cbar_mode="each", - cbar_size="7%", cbar_pad="2%", cbar_location="top") - im = grid[0].imshow(zeroIm.getImage().getArray(), - extent=(x0, x0+dx, y0+dy, y0), norm=norm, - cmap='Spectral') - im = grid[0].scatter(X, Y, c=Z, s=15.*meanDz/dZ, edgecolor='none', norm=norm, - marker='o', cmap='Spectral') - im2 = grid[1].scatter(X, Y, c=resids, edgecolor='none', norm=diffnorm, - marker='s', cmap='seismic') + grid = ImageGrid( + fig, + 111, + nrows_ncols=(1, 2), + axes_pad=0.1, + share_all=True, + label_mode="L", + cbar_mode="each", + cbar_size="7%", + cbar_pad="2%", + cbar_location="top", + ) + im = grid[0].imshow( + zeroIm.getImage().getArray(), extent=(x0, x0 + dx, y0 + dy, y0), norm=norm, cmap="Spectral" + ) + im = grid[0].scatter( + X, Y, c=Z, s=15.0 * meanDz / dZ, edgecolor="none", norm=norm, marker="o", cmap="Spectral" + ) + im2 = grid[1].scatter(X, Y, c=resids, edgecolor="none", norm=diffnorm, marker="s", cmap="seismic") grid.cbar_axes[0].colorbar(im) grid.cbar_axes[1].colorbar(im2) - grid[0].axis([x0, x0+dx, y0+dy, y0]) - grid[1].axis([x0, x0+dx, y0+dy, y0]) + grid[0].axis([x0, x0 + dx, y0 + dy, y0]) + grid[1].axis([x0, x0 + dx, y0 + dy, y0]) grid[0].set_xlabel("model and grid") - grid[1].set_xlabel("residuals. rms = %0.3f"%(rms)) + grid[1].set_xlabel("residuals. rms = %0.3f" % (rms)) if lsstDebug.Info(__name__).savefig: - fig.savefig(lsstDebug.Info(__name__).figpath + self.debugDataIdString + '.png') + fig.savefig(lsstDebug.Info(__name__).figpath + self.debugDataIdString + ".png") if lsstDebug.Info(__name__).display: plt.show() plt.clf() @@ -667,13 +700,16 @@ def _gridImage(self, maskedImage, binsize, statsFlag): for ymin, ymax in zip(yedges[0:-1], yedges[1:]): for xmin, xmax in zip(xedges[0:-1], xedges[1:]): - subBBox = geom.Box2I(geom.PointI(int(x0 + xmin), int(y0 + ymin)), - geom.PointI(int(x0 + xmax-1), int(y0 + ymax-1))) + subBBox = geom.Box2I( + geom.PointI(int(x0 + xmin), int(y0 + ymin)), + geom.PointI(int(x0 + xmax - 1), int(y0 + ymax - 1)), + ) subIm = afwImage.MaskedImageF(maskedImage, subBBox, afwImage.PARENT, False) - stats = afwMath.makeStatistics(subIm, - afwMath.MEAN | afwMath.MEANCLIP | afwMath.MEDIAN - | afwMath.NPOINT | afwMath.STDEV, - self.sctrl) + stats = afwMath.makeStatistics( + subIm, + afwMath.MEAN | afwMath.MEANCLIP | afwMath.MEDIAN | afwMath.NPOINT | afwMath.STDEV, + self.sctrl, + ) npoints, _ = stats.getResult(afwMath.NPOINT) if npoints >= 2: stdev, _ = stats.getResult(afwMath.STDEV) @@ -681,7 +717,7 @@ def _gridImage(self, maskedImage, binsize, statsFlag): stdev = self.config.gridStdevEpsilon bgX.append(0.5 * (x0 + xmin + x0 + xmax)) bgY.append(0.5 * (y0 + ymin + y0 + ymax)) - bgdZ.append(stdev/np.sqrt(npoints)) + bgdZ.append(stdev / np.sqrt(npoints)) est, _ = stats.getResult(statsFlag) bgZ.append(est) From 835ca9eb4c8337c0b152e686b720ba199ef25dce Mon Sep 17 00:00:00 2001 From: Lee Kelvin Date: Fri, 19 Jul 2024 06:58:20 -0700 Subject: [PATCH 05/13] Refactor by LSK --- python/lsst/pipe/tasks/matchBackgrounds.py | 532 ++++++++++----------- 1 file changed, 265 insertions(+), 267 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 91db387af..82ee3e6d9 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -19,176 +19,205 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -__all__ = ["MatchBackgroundsConfig", "MatchBackgroundsTask"] +__all__ = ["MatchBackgroundsConnections", "MatchBackgroundsConfig", "MatchBackgroundsTask"] -import pdb - -import lsst.afw.image as afwImage -import lsst.afw.math as afwMath -import lsst.geom as geom -import lsst.pex.config as pexConfig -import lsst.pipe.base as pipeBase -import lsst.pipe.base.connectionTypes as cT import lsstDebug import numpy as np -from lsst.pipe.base import PipelineTaskConfig, PipelineTaskConnections +from lsst.afw.image import LOCAL, PARENT, Mask, MaskedImageF +from lsst.afw.math import ( + MEAN, + MEANCLIP, + MEANSQUARE, + MEDIAN, + NPOINT, + STDEV, + VARIANCE, + ApproximateControl, + BackgroundControl, + BackgroundList, + StatisticsControl, + makeBackground, + makeStatistics, + stringToInterpStyle, + stringToStatisticsProperty, + stringToUndersampleStyle, +) +from lsst.geom import Box2D, Box2I, PointI +from lsst.pex.config import ChoiceField, Field, ListField, RangeField +from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct, TaskError +from lsst.pipe.base.connectionTypes import Input, Output from lsst.utils.timer import timeMethod class MatchBackgroundsConnections( PipelineTaskConnections, - dimensions=("tract", "patch", "band", "skymap"), + dimensions=("skymap", "tract", "patch", "band"), defaultTemplates={ "inputCoaddName": "deep", "outputCoaddName": "deep", - "warpType": "direct", + "warpType": "psfMatched", "warpTypeSuffix": "", }, ): - psfMatchedWarps = pipeBase.connectionTypes.Input( - doc=("PSF-matched warps to be subtracted from the reference warp"), - name="{inputCoaddName}Coadd_psfMatchedWarp", + warps = Input( + doc=("Warps used to construct a list of matched backgrounds."), + name="{inputCoaddName}Coadd_{warpType}Warp", storageClass="ExposureF", - dimensions=("tract", "patch", "skymap", "visit"), + dimensions=("skymap", "tract", "patch", "visit"), deferLoad=True, multiple=True, ) - # The reference exposure needs to be another psfMatchedWarp - # Let the software choose it automatically for now - # refMatchedWarp = cT.Input( - # doc="Reference exposure", - # dimensions=("visit", "detector", "band"), - # storageClass="ExposureF", - # name="calexp", - # ) - # This needs to be the models of each differential BG in warped coords - backgroundInfoList = cT.Output( - doc="List of differential backgrounds, w/goodness of fit params", + backgroundInfoList = Output( + doc="List of differential backgrounds, with goodness of fit params.", name="psfMatchedWarpBackground_diff", # This needs to change - dimensions=("tract", "patch", "skymap", "visit"), + dimensions=("skymap", "tract", "patch", "visit"), storageClass="Background", multiple=True, ) -# class MatchBackgroundsConfig(pexConfig.Config): class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgroundsConnections): - usePolynomial = pexConfig.Field( - dtype=bool, - doc="Fit background difference with Chebychev polynomial interpolation " - "(using afw.math.Approximate)? If False, fit with spline interpolation using afw.math.Background", + # Reference warp selection + refWarpVisit = Field[int]( + doc="Visit ID of the reference warp. If None, the best warp is chosen from the list of warps.", + optional=True, + ) + bestRefWeightCoverage = RangeField( + dtype=float, + doc="Coverage weight (Number of pixels overlapping the patch) when calculating the best reference " + "exposure. Higher weights prefer exposures with high coverage. Ignored when a ref visit supplied.", + default=0.4, + min=0.0, + max=1.0, + ) + bestRefWeightVariance = RangeField( + dtype=float, + doc="Image variance weight when calculating the best reference exposure. " + "Higher weights prefers exposures with low image variances. Ignored when ref visit supplied", + default=0.4, + min=0.0, + max=1.0, + ) + bestRefWeightLevel = RangeField( + dtype=float, + doc="Mean background level weight when calculating the best reference exposure. " + "Higher weights prefer exposures with low mean background levels. Ignored when ref visit supplied.", + default=0.2, + min=0.0, + max=1.0, + ) + + # Background matching + usePolynomial = Field[bool]( + doc="Fit background difference with a Chebychev polynomial interpolation? " + "If False, fit with spline interpolation instead.", default=False, ) - order = pexConfig.Field( - dtype=int, - doc="Order of Chebyshev polynomial background model. Ignored if usePolynomial False", + order = Field[int]( + doc="Order of Chebyshev polynomial background model. Ignored if ``usePolynomial=False``.", default=8, ) - badMaskPlanes = pexConfig.ListField( - doc="Names of mask planes to ignore while estimating the background", - dtype=str, + badMaskPlanes = ListField[str]( + doc="Names of mask planes to ignore while estimating the background.", default=["NO_DATA", "DETECTED", "DETECTED_NEGATIVE", "SAT", "BAD", "INTRP", "CR"], - itemCheck=lambda x: x in afwImage.Mask().getMaskPlaneDict(), + itemCheck=lambda x: x in Mask().getMaskPlaneDict(), ) - gridStatistic = pexConfig.ChoiceField( + gridStatistic = ChoiceField( dtype=str, - doc="Type of statistic to estimate pixel value for the grid points", + doc="Type of statistic to estimate pixel value for the grid points.", default="MEAN", allowed={"MEAN": "mean", "MEDIAN": "median", "MEANCLIP": "clipped mean"}, ) - undersampleStyle = pexConfig.ChoiceField( - doc="Behaviour if there are too few points in grid for requested interpolation style. " - "Note: INCREASE_NXNYSAMPLE only allowed for usePolynomial=True.", + undersampleStyle = ChoiceField( dtype=str, + doc="Behaviour if there are too few points in the grid for requested interpolation style. " + "Note: choice ``INCREASE_NXNYSAMPLE`` only allowed for ``usePolynomial=True``.", default="REDUCE_INTERP_ORDER", allowed={ - "THROW_EXCEPTION": "throw an exception if there are too few points", + "THROW_EXCEPTION": "throw an exception if there are too few points.", "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.", "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.", }, ) - binSize = pexConfig.Field( - doc="Bin size for gridding the difference image and fitting a spatial model", dtype=int, default=256 + binSize = Field[int]( + doc="Bin size for gridding the difference image and fitting a spatial model.", + default=256, ) - interpStyle = pexConfig.ChoiceField( + interpStyle = ChoiceField( dtype=str, - doc="Algorithm to interpolate the background values; ignored if usePolynomial is True" - "Maps to an enum; see afw.math.Background", + doc="Algorithm to interpolate the background values; ignored if ``usePolynomial=True``." + "Maps to an enum; see afw.math.Background for more information.", default="AKIMA_SPLINE", allowed={ - "CONSTANT": "Use a single constant value", - "LINEAR": "Use linear interpolation", - "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints", - "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers", - "NONE": "No background estimation is to be attempted", + "CONSTANT": "Use a single constant value.", + "LINEAR": "Use linear interpolation.", + "NATURAL_SPLINE": "A cubic spline with zero second derivative at endpoints.", + "AKIMA_SPLINE": "A higher-level non-linear spline that is more robust to outliers.", + "NONE": "No background estimation is to be attempted.", }, ) - numSigmaClip = pexConfig.Field( - dtype=int, doc="Sigma for outlier rejection; ignored if gridStatistic != 'MEANCLIP'.", default=3 + numSigmaClip = Field[int]( + doc="Sigma for outlier rejection. Ignored if ``gridStatistic != 'MEANCLIP'``.", + default=3, ) - numIter = pexConfig.Field( - dtype=int, - doc="Number of iterations of outlier rejection; ignored if gridStatistic != 'MEANCLIP'.", + numIter = Field[int]( + doc="Number of iterations of outlier rejection. Ignored if ``gridStatistic != 'MEANCLIP'``.", default=2, ) - bestRefWeightCoverage = pexConfig.RangeField( - dtype=float, - doc="Weight given to coverage (number of pixels that overlap with patch), " - "when calculating best reference exposure. Higher weight prefers exposures with high coverage." - "Ignored when reference visit is supplied", - default=0.4, - min=0.0, - max=1.0, - ) - bestRefWeightVariance = pexConfig.RangeField( - dtype=float, - doc="Weight given to image variance when calculating best reference exposure. " - "Higher weight prefers exposures with low image variance. Ignored when reference visit is supplied", - default=0.4, - min=0.0, - max=1.0, - ) - bestRefWeightLevel = pexConfig.RangeField( - dtype=float, - doc="Weight given to mean background level when calculating best reference exposure. " - "Higher weight prefers exposures with low mean background level. " - "Ignored when reference visit is supplied.", - default=0.2, - min=0.0, - max=1.0, - ) - approxWeighting = pexConfig.Field( - dtype=bool, - doc=( - "Use inverse-variance weighting when approximating background offset model? " - "This will fail when the background offset is constant " - "(this is usually only the case in testing with artificial images)." - "(usePolynomial=True)" - ), + + approxWeighting = Field[bool]( + doc="Use inverse-variance weighting when approximating the background offset model? This will fail " + "when the background offset is constant (usually only the case in testing with artificial images)." + "Only applied if ``usePolynomial=True``.", default=True, ) - gridStdevEpsilon = pexConfig.RangeField( + gridStdevEpsilon = RangeField( dtype=float, - doc="Tolerance on almost zero standard deviation in a background-offset grid bin. " - "If all bins have a standard deviation below this value, the background offset model " - "is approximated without inverse-variance weighting. (usePolynomial=True)", + doc="Tolerance on almost zero standard deviation in a background-offset grid bin. If all bins have a " + "standard deviation below this value, the background offset model is approximated without " + "inverse-variance weighting. Only applied if ``usePolynomial=True``.", default=1e-8, min=0.0, ) -class MatchBackgroundsTask(pipeBase.PipelineTask): +class MatchBackgroundsTask(PipelineTask): + """Match the backgrounds of a list of warped exposures to a reference. + + This task is a part of the background subtraction pipeline. + It matches the backgrounds of a list of science exposures to a reference + science exposure. + The reference exposure is chosen from the list of science exposures by + minimizing a cost function that penalizes high variance, high background + level, and low coverage. + The cost function is a weighted sum of these three metrics. + The weights are set by the config parameters: + - ``bestRefWeightCoverage`` + - ``bestRefWeightVariance`` + - ``bestRefWeightLevel`` + + + Attributes + ---------- + config : `MatchBackgroundsConfig` + Configuration for this task. + statsCtrl : `~lsst.afw.math.StatisticsControl` + Statistics control object. + """ + ConfigClass = MatchBackgroundsConfig + config: MatchBackgroundsConfig _DefaultName = "matchBackgrounds" def __init__(self, *args, **kwargs): super().__init__(**kwargs) - - self.sctrl = afwMath.StatisticsControl() - self.sctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes)) - self.sctrl.setNanSafe(True) + self.statsCtrl = StatisticsControl() + # TODO: Check that setting the mask planes here work - these planes + # can vary from exposure to exposure, I think? + self.statsCtrl.setAndMask(Mask.getPlaneBitMask(self.config.badMaskPlanes)) + self.statsCtrl.setNanSafe(True) def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) @@ -196,27 +225,19 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): butlerQC.put(outputs, outputRefs) @timeMethod - def run(self, psfMatchedWarps): - """Match the backgrounds of a list of science exposures to a reference - science exposure. + def run(self, warps): + """Match the backgrounds of a list of warped exposures to a reference. - Choose a refMatchedWarp automatically if none supplied. - (above is legacy: right now this only does it automatically) + A reference warp will be chosen automatically if none is supplied. Parameters ---------- - psfMatchedWarps : `list` of `lsst.afw.image.Exposure` - List of warped science exposures to be background-matched; - all exposures must exist. - refMatchedWarp : `lsst.afw.image.Exposure`, optional - Reference science exposure. - If None, then this task selects the best exposures from - psfMatchedWarps. - If not None then must be one of the exposures in psfMatchedWarps. + warps : `list`[`~lsst.afw.image.Exposure`] + List of warped science exposures to be background-matched. Returns ------- - result : `lsst.afw.math._backgroundList.BackgroundList` + result : `~lsst.afw.math.BackgroundList` Differential background model Add this to the science exposure to match the reference exposure. @@ -225,30 +246,12 @@ def run(self, psfMatchedWarps): RuntimeError Raised if an exposure does not exist on disk. """ - numExp = len(psfMatchedWarps) - if numExp < 1: - raise pipeBase.TaskError("No exposures to match") - - # refInd: index in psfMatchedWarps corresponding to refMatchedWarp - # refInd = None - # if refMatchedWarp is None: - # select the best reference exposure from psfMatchedWarps - # note: just selecting from the set for testing purposes - refInd = self.selectRefExposure( - psfMatchedWarps=psfMatchedWarps, - ) - refMatchedWarp = psfMatchedWarps[refInd] - - # refIndSet: indices of exposures in psfMatchedWarps matching ref. - # It is used to avoid background-matching an exposure to itself. - # It is a list because it is possible (though unlikely) that - # psfMatchedWarps will contain duplicates. - refId = refMatchedWarp.dataId - expIds = [exp.dataId for exp in psfMatchedWarps] - try: - refIndSet = [expIds.index(i) for i in [refId]] - except ValueError: - raise RuntimeError("Internal error: selected reference %s not found in psfMatchedWarps") + if (numExp := len(warps)) < 1: + raise TaskError("No exposures to match") + + # Define a reference warp; 'warps' is modified in-place to exclude it + refWarp = self._defineWarps(warps=warps, refWarpVisit=self.config.refWarpVisit) + breakpoint() # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this @@ -261,9 +264,9 @@ def run(self, psfMatchedWarps): self.log.info("Matching %d Exposures", numExp) # Creating a null BackgroundList object by fitting a blank image - statsFlag = getattr(afwMath, self.config.gridStatistic) - self.sctrl.setNumSigmaClip(self.config.numSigmaClip) - self.sctrl.setNumIter(self.config.numIter) + statsFlag = stringToStatisticsProperty(self.config.gridStatistic) + self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) + self.statsCtrl.setNumIter(self.config.numIter) # TODO: refactor below to construct blank bg model im = refExposure.getMaskedImage() @@ -279,20 +282,20 @@ def run(self, psfMatchedWarps): if height % self.config.binSize != 0: ny += 1 - bctrl = afwMath.BackgroundControl(nx, ny, self.sctrl, statsFlag) + bctrl = BackgroundControl(nx, ny, self.statsCtrl, statsFlag) bctrl.setUndersampleStyle(self.config.undersampleStyle) - bkgd = afwMath.makeBackground(blankIm, bctrl) + bkgd = makeBackground(blankIm, bctrl) backgroundInfoList = [] - for ind, exp in enumerate(psfMatchedWarps): + for ind, exp in enumerate(warps): if ind in refIndSet: - backgroundInfoStruct = afwMath.BackgroundList( + backgroundInfoStruct = BackgroundList( ( bkgd, - afwMath.stringToInterpStyle(self.config.interpStyle), - afwMath.stringToUndersampleStyle(self.config.undersampleStyle), - afwMath.ApproximateControl.UNKNOWN, + stringToInterpStyle(self.config.interpStyle), + stringToUndersampleStyle(self.config.undersampleStyle), + ApproximateControl.UNKNOWN, 0, 0, False, @@ -316,12 +319,12 @@ def run(self, psfMatchedWarps): backgroundInfoStruct.isReference = False except Exception as e: self.log.warning("Failed to fit background %s: %s", exp.dataId, e) - backgroundInfoStruct = afwMath.BackgroundList( + backgroundInfoStruct = BackgroundList( ( bkgd, - afwMath.stringToInterpStyle(self.config.interpStyle), - afwMath.stringToUndersampleStyle(self.config.undersampleStyle), - afwMath.ApproximateControl.UNKNOWN, + stringToInterpStyle(self.config.interpStyle), + stringToUndersampleStyle(self.config.undersampleStyle), + ApproximateControl.UNKNOWN, 0, 0, False, @@ -330,75 +333,91 @@ def run(self, psfMatchedWarps): backgroundInfoList.append(backgroundInfoStruct) - return pipeBase.Struct(backgroundInfoList=backgroundInfoList) + return Struct(backgroundInfoList=backgroundInfoList) @timeMethod - def selectRefExposure(self, psfMatchedWarps): - """Find best exposure to use as the reference exposure. + def _defineWarps(self, warps, refWarpVisit=None): + """Define the reference warp and list of comparison warps. - Calculate an appropriate reference exposure by minimizing a cost - function that penalizes high variance, high background level, and low - coverage. Use the following config parameters: - - bestRefWeightCoverage - - bestRefWeightVariance - - bestRefWeightLevel + If no reference warp ID is supplied, this method calculates an + appropriate reference exposure from the supplied list of warps by + minimizing a cost function that penalizes high variance, high + background level, and low coverage. + + To find a reference warp, the following config parameters are used: + - ``bestRefWeightCoverage`` + - ``bestRefWeightVariance`` + - ``bestRefWeightLevel`` Parameters ---------- - psfMatchedWarps : `list` of `lsst.afw.image.Exposure` - List of science exposures. - If an exposure is not found, it is skipped with a warning. + warps : `list`[`~lsst.afw.image.Exposure`] + List of warped science exposures. Returns ------- - bestIdx : `int` - Index of best exposure. - - Raises - ------ - RuntimeError - Raised if none of the exposures in psfMatchedWarps are found. + refWarp : `~lsst.afw.image.Exposure` + Reference warped exposure. + compWarps : `list`[`~lsst.afw.image.Exposure`] + List of warped science exposures to compare to the reference. """ - self.log.info("Calculating best reference visit") - varList = [] - meanBkgdLevelList = [] - coverageList = [] - - for exp in psfMatchedWarps: - exposure = exp.get() - # Convert images to nJy before doing statistics + # User a reference visit, if one has been supplied + if refWarpVisit: + warpVisits = [warp.dataId["visit"] for warp in warps] try: - maskedImage = self._fluxScale(exposure) - except Exception: - # need to put a place holder in Arr - varList.append(np.nan) - meanBkgdLevelList.append(np.nan) - coverageList.append(np.nan) - continue - statObjIm = afwMath.makeStatistics( - maskedImage.getImage(), - maskedImage.getMask(), - afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, - self.sctrl, - ) - meanVar, meanVarErr = statObjIm.getResult(afwMath.VARIANCE) - meanBkgdLevel, meanBkgdLevelErr = statObjIm.getResult(afwMath.MEAN) - npoints, npointsErr = statObjIm.getResult(afwMath.NPOINT) - varList.append(meanVar) - meanBkgdLevelList.append(meanBkgdLevel) - coverageList.append(npoints) - if not coverageList: - raise pipeBase.TaskError("None of the candidates exist; cannot select best reference exposure") - - # Normalize metrics to range from 0 to 1 - varArr = np.array(varList) / np.nanmax(varList) - meanBkgdLevelArr = np.array(meanBkgdLevelList) / np.nanmax(meanBkgdLevelList) - coverageArr = np.nanmin(coverageList) / np.array(coverageList) - - costFunctionArr = self.config.bestRefWeightVariance * varArr - costFunctionArr += self.config.bestRefWeightLevel * meanBkgdLevelArr - costFunctionArr += self.config.bestRefWeightCoverage * coverageArr - return np.nanargmin(costFunctionArr) + refWarp = warps.pop(warpVisits.index(refWarpVisit)) + self.log.info("Using user-supplied reference visit %d", refWarpVisit) + return refWarp + except ValueError: + raise TaskError(f"Reference visit {refWarpVisit} is not found in the list of warps.") + + # Extract mean/var/npoints for each warp + warpMeans = [] + warpVars = [] + warpNPoints = [] + for warpDDH in warps: + warp = warpDDH.get() + warp.image.array *= warp.getPhotoCalib().instFluxToNanojansky(1) # Convert image plane to nJy + warpStats = makeStatistics(warp.image, warp.mask, MEAN | VARIANCE | NPOINT, self.statsCtrl) + warpMean, _ = warpStats.getResult(MEAN) + warpVar, _ = warpStats.getResult(VARIANCE) + warpNPoint, _ = warpStats.getResult(NPOINT) + warpMeans.append(warpMean) + warpVars.append(warpVar) + warpNPoints.append(int(warpNPoint)) + + # Normalize mean/var/npoints to range from 0 to 1 + warpMeansFrac = np.array(warpMeans) / np.nanmax(warpMeans) + warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) + warpNPointsFrac = np.nanmin(warpNPoints) / np.array(warpNPoints) + + # Calculate cost function values + costFunctionVals = self.config.bestRefWeightLevel * warpMeansFrac + costFunctionVals += self.config.bestRefWeightVariance * warpVarsFrac + costFunctionVals += self.config.bestRefWeightCoverage * warpNPointsFrac + + refWarp = warps.pop(np.nanargmin(costFunctionVals)) + self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) + return refWarp + + def _fluxScale(self, exposure): + """Scales image to nJy flux using photometric calibration. + + Parameters + ---------- + exposure: `` + Exposure to scale. + + Returns + ------- + maskedImage: `` + Flux-scaled masked exposure. + """ + maskedImage = exposure.getMaskedImage() + fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) + exposure.image.array *= fluxZp + + return maskedImage @timeMethod def matchBackgrounds(self, refExposure, sciExposure): @@ -470,9 +489,9 @@ def matchBackgrounds(self, refExposure, sciExposure): "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" % (wSci, hSci, wRef, hRef) ) - statsFlag = getattr(afwMath, self.config.gridStatistic) - self.sctrl.setNumSigmaClip(self.config.numSigmaClip) - self.sctrl.setNumIter(self.config.numIter) + statsFlag = stringToStatisticsProperty(self.config.gridStatistic) + self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) + self.statsCtrl.setNumIter(self.config.numIter) im = refExposure.getMaskedImage() diffMI = im.Factory(im, True) @@ -487,10 +506,10 @@ def matchBackgrounds(self, refExposure, sciExposure): if height % self.config.binSize != 0: ny += 1 - bctrl = afwMath.BackgroundControl(nx, ny, self.sctrl, statsFlag) + bctrl = BackgroundControl(nx, ny, self.statsCtrl, statsFlag) bctrl.setUndersampleStyle(self.config.undersampleStyle) - bkgd = afwMath.makeBackground(diffMI, bctrl) + bkgd = makeBackground(diffMI, bctrl) # Some config and input checks if config.usePolynomial: # 1) Check that order/bin size make sense: @@ -512,7 +531,7 @@ def matchBackgrounds(self, refExposure, sciExposure): newBinSize = (minNumberGridPoints * self.config.binSize) // (self.config.order + 1) bctrl.setNxSample(newBinSize) bctrl.setNySample(newBinSize) - bkgd = afwMath.makeBackground(diffMI, bctrl) # do over + bkgd = makeBackground(diffMI, bctrl) # do over self.log.warning("Decreasing binsize to %d", newBinSize) # If there is no variance in any image pixels, do not weight bins by inverse variance @@ -522,10 +541,10 @@ def matchBackgrounds(self, refExposure, sciExposure): # Add offset to sciExposure try: if self.config.usePolynomial: - actrl = afwMath.ApproximateControl( - afwMath.ApproximateControl.CHEBYSHEV, order, order, weightByInverseVariance + actrl = ApproximateControl( + ApproximateControl.CHEBYSHEV, order, order, weightByInverseVariance ) - undersampleStyle = getattr(afwMath, self.config.undersampleStyle) + undersampleStyle = stringToUndersampleStyle(self.config.undersampleStyle) approx = bkgd.getApproximate(actrl, undersampleStyle) bkgdImage = approx.getImage() else: @@ -541,32 +560,30 @@ def matchBackgrounds(self, refExposure, sciExposure): # Need RMS from fit: 2895 will replace this: rms = 0.0 - X, Y, Z, dZ = self._gridImage(diffMI, self.config.binSize, statsFlag) + bgX, bgY, bgZ, bgdZ = self._gridImage(diffMI, self.config.binSize, statsFlag) x0, y0 = diffMI.getXY0() - modelValueArr = np.empty(len(Z)) - for i in range(len(X)): - modelValueArr[i] = bkgdImage[int(X[i] - x0), int(Y[i] - y0), afwImage.LOCAL] - resids = Z - modelValueArr + modelValueArr = np.empty(len(bgZ)) + for i in range(len(bgX)): + modelValueArr[i] = bkgdImage[int(bgX[i] - x0), int(bgY[i] - y0), LOCAL] + resids = bgZ - modelValueArr rms = np.sqrt(np.mean(resids[~np.isnan(resids)] ** 2)) if lsstDebug.Info(__name__).savefits: sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciMatchedExposure.fits") if lsstDebug.Info(__name__).savefig: - bbox = geom.Box2D(refExposure.getMaskedImage().getBBox()) + bbox = Box2D(refExposure.getMaskedImage().getBBox()) try: - self._debugPlot(X, Y, Z, dZ, bkgdImage, bbox, modelValueArr, resids) + self._debugPlot(bgX, bgY, bgZ, bgdZ, bkgdImage, bbox, modelValueArr, resids) except Exception as e: self.log.warning("Debug plot not generated: %s", e) - meanVar = afwMath.makeStatistics( - diffMI.getVariance(), diffMI.getMask(), afwMath.MEANCLIP, self.sctrl - ).getValue() + meanVar = makeStatistics(diffMI.getVariance(), diffMI.getMask(), MEANCLIP, self.statsCtrl).getValue() diffIm = diffMI.getImage() diffIm -= bkgdImage # diffMI should now have a mean ~ 0 del diffIm - mse = afwMath.makeStatistics(diffMI, afwMath.MEANSQUARE, self.sctrl).getValue() + mse = makeStatistics(diffMI, MEANSQUARE, self.statsCtrl).getValue() outBkgd = approx if self.config.usePolynomial else bkgd @@ -578,37 +595,18 @@ def matchBackgrounds(self, refExposure, sciExposure): meanVar, ) # TODO: verify this is correct (borrowed from background.py) - return afwMath.BackgroundList( + return BackgroundList( ( outBkgd, - afwMath.stringToInterpStyle(self.config.interpStyle), - afwMath.stringToUndersampleStyle(self.config.undersampleStyle), - afwMath.ApproximateControl.UNKNOWN, + stringToInterpStyle(self.config.interpStyle), + stringToUndersampleStyle(self.config.undersampleStyle), + ApproximateControl.UNKNOWN, 0, 0, False, ) ) - def _fluxScale(self, exposure): - """Scales image to nJy flux using photometric calibration. - - Parameters - ---------- - exposure: `` - Exposure to scale. - - Returns - ------- - maskedImage: `` - Flux-scaled masked exposure. - """ - maskedImage = exposure.getMaskedImage() - fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) - exposure.image.array *= fluxZp - - return maskedImage - def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): """Generate a plot showing the background fit and residuals. @@ -637,7 +635,7 @@ def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import ImageGrid - zeroIm = afwImage.MaskedImageF(geom.Box2I(bbox)) + zeroIm = MaskedImageF(Box2I(bbox)) zeroIm += modelImage x0, y0 = zeroIm.getXY0() dx, dy = zeroIm.getDimensions() @@ -700,19 +698,19 @@ def _gridImage(self, maskedImage, binsize, statsFlag): for ymin, ymax in zip(yedges[0:-1], yedges[1:]): for xmin, xmax in zip(xedges[0:-1], xedges[1:]): - subBBox = geom.Box2I( - geom.PointI(int(x0 + xmin), int(y0 + ymin)), - geom.PointI(int(x0 + xmax - 1), int(y0 + ymax - 1)), + subBBox = Box2I( + PointI(int(x0 + xmin), int(y0 + ymin)), + PointI(int(x0 + xmax - 1), int(y0 + ymax - 1)), ) - subIm = afwImage.MaskedImageF(maskedImage, subBBox, afwImage.PARENT, False) - stats = afwMath.makeStatistics( + subIm = MaskedImageF(maskedImage, subBBox, PARENT, False) + stats = makeStatistics( subIm, - afwMath.MEAN | afwMath.MEANCLIP | afwMath.MEDIAN | afwMath.NPOINT | afwMath.STDEV, - self.sctrl, + MEAN | MEANCLIP | MEDIAN | NPOINT | STDEV, + self.statsCtrl, ) - npoints, _ = stats.getResult(afwMath.NPOINT) + npoints, _ = stats.getResult(NPOINT) if npoints >= 2: - stdev, _ = stats.getResult(afwMath.STDEV) + stdev, _ = stats.getResult(STDEV) if stdev < self.config.gridStdevEpsilon: stdev = self.config.gridStdevEpsilon bgX.append(0.5 * (x0 + xmin + x0 + xmax)) From feaff2c81d9f4387475142d037b0a2dd0b4977e9 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Thu, 25 Jul 2024 01:40:14 -0700 Subject: [PATCH 06/13] Use sky-subtracted warps for ref-image selection _defineWarps() now rejects any image with all NaNs along any image edge, and creates the cost function using a sky-subtracted image. This sky-subtraction fits a 1st order Chebyshev polynomial to the masked image background. Also fixed a bug from LSK refactor by inserting a blank sky model into the background model list at the chosen reference image index. --- python/lsst/pipe/tasks/matchBackgrounds.py | 185 +++++++++++++-------- 1 file changed, 114 insertions(+), 71 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 82ee3e6d9..6a6578adc 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -23,7 +23,7 @@ import lsstDebug import numpy as np -from lsst.afw.image import LOCAL, PARENT, Mask, MaskedImageF +from lsst.afw.image import LOCAL, PARENT, ImageF, Mask, MaskedImageF from lsst.afw.math import ( MEAN, MEANCLIP, @@ -250,12 +250,11 @@ def run(self, warps): raise TaskError("No exposures to match") # Define a reference warp; 'warps' is modified in-place to exclude it - refWarp = self._defineWarps(warps=warps, refWarpVisit=self.config.refWarpVisit) - breakpoint() + refWarp, refInd = self._defineWarps(warps=warps, refWarpVisit=self.config.refWarpVisit) # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this - refExposure = refMatchedWarp.get() + refExposure = refWarp.get() refMI = self._fluxScale(refExposure) # Also modifies refExposure # TODO: figure out what this was creating and why @@ -270,7 +269,7 @@ def run(self, warps): # TODO: refactor below to construct blank bg model im = refExposure.getMaskedImage() - blankIm = im.Factory(im, True) # Don't do this + blankIm = im.Factory(im, True) blankIm.image.array *= 0 width = blankIm.getWidth() @@ -286,53 +285,43 @@ def run(self, warps): bctrl.setUndersampleStyle(self.config.undersampleStyle) bkgd = makeBackground(blankIm, bctrl) + blank = BackgroundList( + ( + bkgd, + stringToInterpStyle(self.config.interpStyle), + stringToUndersampleStyle(self.config.undersampleStyle), + ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) backgroundInfoList = [] for ind, exp in enumerate(warps): - if ind in refIndSet: - backgroundInfoStruct = BackgroundList( - ( - bkgd, - stringToInterpStyle(self.config.interpStyle), - stringToUndersampleStyle(self.config.undersampleStyle), - ApproximateControl.UNKNOWN, - 0, - 0, - False, - ) + # TODO: simplify this maybe, using only visit IDs? + self.log.info("Matching background of %s to %s", exp.dataId, refWarp.dataId) + toMatchExposure = exp.get() + try: + # TODO: adjust logic to avoid creating spurious MIs like this + toMatchMI = self._fluxScale(toMatchExposure) + + # store a string specifying the visit to label debug plot + # self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList]) + + backgroundInfoStruct = self.matchBackgrounds( + refExposure=refExposure, + sciExposure=toMatchExposure, ) - else: - # TODO: simplify this maybe, using only visit IDs? - self.log.info("Matching background of %s to %s", exp.dataId, refMatchedWarp.dataId) - toMatchExposure = exp.get() - try: - # Seems to require exposure, not masked exposure. - toMatchMI = self._fluxScale(toMatchExposure) - - # store a string specifying the visit to label debug plot - # self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList]) - - backgroundInfoStruct = self.matchBackgrounds( - refExposure=refExposure, - sciExposure=toMatchExposure, - ) - backgroundInfoStruct.isReference = False - except Exception as e: - self.log.warning("Failed to fit background %s: %s", exp.dataId, e) - backgroundInfoStruct = BackgroundList( - ( - bkgd, - stringToInterpStyle(self.config.interpStyle), - stringToUndersampleStyle(self.config.undersampleStyle), - ApproximateControl.UNKNOWN, - 0, - 0, - False, - ) - ) + backgroundInfoStruct.isReference = False + except Exception as e: + self.log.warning("Failed to fit background %s: %s", exp.dataId, e) + backgroundInfoStruct = blank backgroundInfoList.append(backgroundInfoStruct) + # TODO: more elegant solution than inserting blank model at ref ind? + backgroundInfoList.insert(refInd, blank) return Struct(backgroundInfoList=backgroundInfoList) @timeMethod @@ -377,8 +366,28 @@ def _defineWarps(self, warps, refWarpVisit=None): warpNPoints = [] for warpDDH in warps: warp = warpDDH.get() + + # First check if any image edge is all NaN + # If so, don't use + leftBool = np.all(np.isnan(warp.image.array[:, 0])) + rightBool = np.all(np.isnan(warp.image.array[:, warp.image.getHeight() - 1])) + bottomBool = np.all(np.isnan(warp.image.array[0, :])) + topBool = np.all(np.isnan(warp.image.array[warp.image.getWidth() - 1, :])) + if np.any([leftBool, rightBool, bottomBool, topBool]): + continue + warp.image.array *= warp.getPhotoCalib().instFluxToNanojansky(1) # Convert image plane to nJy - warpStats = makeStatistics(warp.image, warp.mask, MEAN | VARIANCE | NPOINT, self.statsCtrl) + + # Select reference based on BG of plane sky-subtracted images + bkgd, __, __, __ = self._setupBackground(warp) + + weightByInverseVariance = self.config.approxWeighting + actrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, weightByInverseVariance) + undersampleStyle = stringToUndersampleStyle(self.config.undersampleStyle) + approx = bkgd.getApproximate(actrl, undersampleStyle) + bgSubImage = ImageF(warp.image.array - approx.getImage().array) + + warpStats = makeStatistics(bgSubImage, warp.mask, MEAN | VARIANCE | NPOINT, self.statsCtrl) warpMean, _ = warpStats.getResult(MEAN) warpVar, _ = warpStats.getResult(VARIANCE) warpNPoint, _ = warpStats.getResult(NPOINT) @@ -386,6 +395,9 @@ def _defineWarps(self, warps, refWarpVisit=None): warpVars.append(warpVar) warpNPoints.append(int(warpNPoint)) + if len(warpNPoints) == 0: + raise TaskError("No suitable reference visit found in list of warps.") + # Normalize mean/var/npoints to range from 0 to 1 warpMeansFrac = np.array(warpMeans) / np.nanmax(warpMeans) warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) @@ -396,21 +408,22 @@ def _defineWarps(self, warps, refWarpVisit=None): costFunctionVals += self.config.bestRefWeightVariance * warpVarsFrac costFunctionVals += self.config.bestRefWeightCoverage * warpNPointsFrac - refWarp = warps.pop(np.nanargmin(costFunctionVals)) + ind = np.nanargmin(costFunctionVals) + refWarp = warps.pop(ind) self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) - return refWarp + return refWarp, ind def _fluxScale(self, exposure): """Scales image to nJy flux using photometric calibration. Parameters ---------- - exposure: `` + exposure: `lsst.afw.image._exposure.ExposureF` Exposure to scale. Returns ------- - maskedImage: `` + maskedImage: `lsst.afw.image._maskedImage.MaskedImageF` Flux-scaled masked exposure. """ maskedImage = exposure.getMaskedImage() @@ -419,6 +432,56 @@ def _fluxScale(self, exposure): return maskedImage + def _setupBackground(self, warp): + """Set up and return a background model container and associated images + and controllers. + + Uses the following config parameters: + - ``gridStatistic`` + - ``numSigmaClip`` + - ``numIter`` + - ``binSize`` + - ``undersampleStyle`` + + Parameters + ---------- + warp: `lsst.afw.image._exposure.ExposureF` + Warped exposure or difference image for which to estimate + background. + + Returns + ------- + bkgd: `lsst.afw.math.BackgroundMI` + Background model container. + bctrl: `lsst.afw.math.BackgroundControl` + Background model control + warpMI: `lsst.afw.image._maskedImage.MaskedImageF` + Masked image associated with warp + statsFlag: `lsst.afw.math.Property` + Flag for grid statistic property (self.config.gridStatistic) + """ + statsFlag = stringToStatisticsProperty(self.config.gridStatistic) + self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) + self.statsCtrl.setNumIter(self.config.numIter) + + warpMI = warp.getMaskedImage() + + width = warpMI.getWidth() + height = warpMI.getHeight() + nx = width // self.config.binSize + if width % self.config.binSize != 0: + nx += 1 + ny = height // self.config.binSize + if height % self.config.binSize != 0: + ny += 1 + + bctrl = BackgroundControl(nx, ny, self.statsCtrl, statsFlag) + bctrl.setUndersampleStyle(self.config.undersampleStyle) + + bkgd = makeBackground(warpMI, bctrl) + + return bkgd, bctrl, warpMI, statsFlag + @timeMethod def matchBackgrounds(self, refExposure, sciExposure): """Match science exposure's background level to that of reference @@ -489,27 +552,7 @@ def matchBackgrounds(self, refExposure, sciExposure): "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" % (wSci, hSci, wRef, hRef) ) - statsFlag = stringToStatisticsProperty(self.config.gridStatistic) - self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) - self.statsCtrl.setNumIter(self.config.numIter) - - im = refExposure.getMaskedImage() - diffMI = im.Factory(im, True) - diffMI -= sciExposure.getMaskedImage() - - width = diffMI.getWidth() - height = diffMI.getHeight() - nx = width // self.config.binSize - if width % self.config.binSize != 0: - nx += 1 - ny = height // self.config.binSize - if height % self.config.binSize != 0: - ny += 1 - - bctrl = BackgroundControl(nx, ny, self.statsCtrl, statsFlag) - bctrl.setUndersampleStyle(self.config.undersampleStyle) - - bkgd = makeBackground(diffMI, bctrl) + bkgd, bctrl, diffMI, statsFlag = self._setupBackground(refExposure) # Some config and input checks if config.usePolynomial: # 1) Check that order/bin size make sense: From 76b73760c8132b3e07efd3fd6bcafaed2621f3aa Mon Sep 17 00:00:00 2001 From: Lee Kelvin Date: Fri, 26 Jul 2024 06:01:08 -0700 Subject: [PATCH 07/13] Refactor by LSK, round 2 --- python/lsst/pipe/tasks/matchBackgrounds.py | 207 ++++++++++----------- 1 file changed, 98 insertions(+), 109 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 6a6578adc..117c8467a 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -23,7 +23,7 @@ import lsstDebug import numpy as np -from lsst.afw.image import LOCAL, PARENT, ImageF, Mask, MaskedImageF +from lsst.afw.image import LOCAL, PARENT, ExposureF, ImageF, Mask, MaskedImageF from lsst.afw.math import ( MEAN, MEANCLIP, @@ -35,6 +35,7 @@ ApproximateControl, BackgroundControl, BackgroundList, + BackgroundMI, StatisticsControl, makeBackground, makeStatistics, @@ -84,11 +85,11 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr doc="Visit ID of the reference warp. If None, the best warp is chosen from the list of warps.", optional=True, ) - bestRefWeightCoverage = RangeField( + bestRefWeightChi2 = RangeField( dtype=float, - doc="Coverage weight (Number of pixels overlapping the patch) when calculating the best reference " - "exposure. Higher weights prefer exposures with high coverage. Ignored when a ref visit supplied.", - default=0.4, + doc="Mean background goodness of fit statistic weight when calculating the best reference exposure. " + "Higher weights prefer exposures with flatter backgrounds. Ignored when ref visit supplied.", + default=0.2, min=0.0, max=1.0, ) @@ -100,10 +101,18 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr min=0.0, max=1.0, ) - bestRefWeightLevel = RangeField( + bestRefWeightGlobalCoverage = RangeField( dtype=float, - doc="Mean background level weight when calculating the best reference exposure. " - "Higher weights prefer exposures with low mean background levels. Ignored when ref visit supplied.", + doc="Global coverage weight (total number of valid pixels) when calculating the best reference " + "exposure. Higher weights prefer exposures with high coverage. Ignored when a ref visit supplied.", + default=0.2, + min=0.0, + max=1.0, + ) + bestRefWeightEdgeCoverage = RangeField( + dtype=float, + doc="Edge coverage weight (number of valid edge pixels) when calculating the best reference " + "exposure. Higher weights prefer exposures with high coverage. Ignored when a ref visit supplied.", default=0.2, min=0.0, max=1.0, @@ -143,7 +152,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr ) binSize = Field[int]( doc="Bin size for gridding the difference image and fitting a spatial model.", - default=256, + default=1024, ) interpStyle = ChoiceField( dtype=str, @@ -213,11 +222,16 @@ class MatchBackgroundsTask(PipelineTask): def __init__(self, *args, **kwargs): super().__init__(**kwargs) + self.statsFlag = stringToStatisticsProperty(self.config.gridStatistic) self.statsCtrl = StatisticsControl() # TODO: Check that setting the mask planes here work - these planes # can vary from exposure to exposure, I think? self.statsCtrl.setAndMask(Mask.getPlaneBitMask(self.config.badMaskPlanes)) self.statsCtrl.setNanSafe(True) + self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) + self.statsCtrl.setNumIter(self.config.numIter) + self.stringToInterpStyle = stringToInterpStyle(self.config.interpStyle) + self.undersampleStyle = stringToUndersampleStyle(self.config.undersampleStyle) def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) @@ -333,154 +347,129 @@ def _defineWarps(self, warps, refWarpVisit=None): minimizing a cost function that penalizes high variance, high background level, and low coverage. - To find a reference warp, the following config parameters are used: - - ``bestRefWeightCoverage`` - - ``bestRefWeightVariance`` - - ``bestRefWeightLevel`` - Parameters ---------- - warps : `list`[`~lsst.afw.image.Exposure`] - List of warped science exposures. + warps : `list`[`~lsst.daf.butler.DeferredDatasetHandle`] + List of warped exposures (of type `~lsst.afw.image.ExposureF`). + refWarpVisit : `int`, optional + Visit ID of the reference warp. + If None, the best warp is chosen from the list of warps. Returns ------- - refWarp : `~lsst.afw.image.Exposure` + refWarp : `~lsst.afw.image.ExposureF` Reference warped exposure. - compWarps : `list`[`~lsst.afw.image.Exposure`] - List of warped science exposures to compare to the reference. + refWarpIndex : `int` + Index of the reference removed from the list of warps. + + Notes + ----- + This method modifies the input list of warps in place by removing the + reference warp from it. """ - # User a reference visit, if one has been supplied + # User-defined reference visit, if one has been supplied if refWarpVisit: - warpVisits = [warp.dataId["visit"] for warp in warps] + warpVisits = [warpDDH.dataId["visit"] for warpDDH in warps] try: - refWarp = warps.pop(warpVisits.index(refWarpVisit)) + refWarpIndex = warpVisits.index(refWarpVisit) + refWarpDDH = warps.pop(refWarpIndex) self.log.info("Using user-supplied reference visit %d", refWarpVisit) - return refWarp + return refWarpDDH.get(), refWarpIndex except ValueError: raise TaskError(f"Reference visit {refWarpVisit} is not found in the list of warps.") # Extract mean/var/npoints for each warp - warpMeans = [] - warpVars = [] - warpNPoints = [] + warpChi2s = [] # Background goodness of fit + warpVars = [] # Variance + warpNPointsGlobal = [] # Global coverage + warpNPointsEdge = [] # Edge coverage for warpDDH in warps: warp = warpDDH.get() - - # First check if any image edge is all NaN - # If so, don't use - leftBool = np.all(np.isnan(warp.image.array[:, 0])) - rightBool = np.all(np.isnan(warp.image.array[:, warp.image.getHeight() - 1])) - bottomBool = np.all(np.isnan(warp.image.array[0, :])) - topBool = np.all(np.isnan(warp.image.array[warp.image.getWidth() - 1, :])) - if np.any([leftBool, rightBool, bottomBool, topBool]): - continue - - warp.image.array *= warp.getPhotoCalib().instFluxToNanojansky(1) # Convert image plane to nJy - - # Select reference based on BG of plane sky-subtracted images - bkgd, __, __, __ = self._setupBackground(warp) - - weightByInverseVariance = self.config.approxWeighting - actrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, weightByInverseVariance) - undersampleStyle = stringToUndersampleStyle(self.config.undersampleStyle) - approx = bkgd.getApproximate(actrl, undersampleStyle) - bgSubImage = ImageF(warp.image.array - approx.getImage().array) - - warpStats = makeStatistics(bgSubImage, warp.mask, MEAN | VARIANCE | NPOINT, self.statsCtrl) - warpMean, _ = warpStats.getResult(MEAN) + instFluxToNanojansky = warp.getPhotoCalib().instFluxToNanojansky(1) + warp.image *= instFluxToNanojansky # Images in nJy to facilitate difference imaging + warp.variance *= instFluxToNanojansky + warpBg, _ = self._makeBackground(warp) + + # Return an approximation to the background + approxCtrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, self.config.approxWeighting) + warpApprox = warpBg.getApproximate(approxCtrl, self.undersampleStyle) + warpBgSub = ImageF(warp.image.array - warpApprox.getImage().array) + + warpStats = makeStatistics(warpBgSub, warp.mask, VARIANCE | NPOINT, self.statsCtrl) + # TODO: need to modify this to account for the background mask + warpChi2 = np.nansum(warpBgSub.array**2 / warp.variance.array) warpVar, _ = warpStats.getResult(VARIANCE) - warpNPoint, _ = warpStats.getResult(NPOINT) - warpMeans.append(warpMean) + warpNPointGlobal, _ = warpStats.getResult(NPOINT) + warpNPointEdge = ( + np.sum(~np.isnan(warp.image.array[:, 0])) # Left edge + + np.sum(~np.isnan(warp.image.array[:, -1])) # Right edge + + np.sum(~np.isnan(warp.image.array[0, :])) # Bottom edge + + np.sum(~np.isnan(warp.image.array[-1, :])) # Top edge + ) + warpChi2s.append(warpChi2) warpVars.append(warpVar) - warpNPoints.append(int(warpNPoint)) - - if len(warpNPoints) == 0: - raise TaskError("No suitable reference visit found in list of warps.") + warpNPointsGlobal.append(int(warpNPointGlobal)) + warpNPointsEdge.append(warpNPointEdge) # Normalize mean/var/npoints to range from 0 to 1 - warpMeansFrac = np.array(warpMeans) / np.nanmax(warpMeans) + warpChi2sFrac = np.array(warpChi2s) / np.nanmax(warpChi2s) warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) - warpNPointsFrac = np.nanmin(warpNPoints) / np.array(warpNPoints) + warpNPointsGlobalFrac = np.nanmin(warpNPointsGlobal) / np.array(warpNPointsGlobal) + warpNPointsEdgeFrac = np.nanmin(warpNPointsEdge) / np.array(warpNPointsEdge) # Calculate cost function values - costFunctionVals = self.config.bestRefWeightLevel * warpMeansFrac + costFunctionVals = self.config.bestRefWeightChi2 * warpChi2sFrac costFunctionVals += self.config.bestRefWeightVariance * warpVarsFrac - costFunctionVals += self.config.bestRefWeightCoverage * warpNPointsFrac + costFunctionVals += self.config.bestRefWeightGlobalCoverage * warpNPointsGlobalFrac + costFunctionVals += self.config.bestRefWeightEdgeCoverage * warpNPointsEdgeFrac ind = np.nanargmin(costFunctionVals) refWarp = warps.pop(ind) self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) return refWarp, ind - def _fluxScale(self, exposure): - """Scales image to nJy flux using photometric calibration. + def _makeBackground(self, warp: ExposureF) -> tuple[BackgroundMI, BackgroundControl]: + """Generate a simple binned background masked image for warped data. Parameters ---------- - exposure: `lsst.afw.image._exposure.ExposureF` - Exposure to scale. + warp: `~lsst.afw.image.ExposureF` + Warped exposure for which to estimate background. Returns ------- - maskedImage: `lsst.afw.image._maskedImage.MaskedImageF` - Flux-scaled masked exposure. + warpBgMI: `~lsst.afw.math.BackgroundMI` + Background-subtracted masked image. + bgCtrl: `~lsst.afw.math.BackgroundControl` + Background control object. """ - maskedImage = exposure.getMaskedImage() - fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) - exposure.image.array *= fluxZp + nx = warp.getWidth() // self.config.binSize + ny = warp.getHeight() // self.config.binSize - return maskedImage + bgCtrl = BackgroundControl(nx, ny, self.statsCtrl, self.statsFlag) + bgCtrl.setUndersampleStyle(self.config.undersampleStyle) + warpBgMI = makeBackground(warp.getMaskedImage(), bgCtrl) - def _setupBackground(self, warp): - """Set up and return a background model container and associated images - and controllers. + return warpBgMI, bgCtrl - Uses the following config parameters: - - ``gridStatistic`` - - ``numSigmaClip`` - - ``numIter`` - - ``binSize`` - - ``undersampleStyle`` + def _fluxScale(self, exposure): + """Scales image to nJy flux using photometric calibration. Parameters ---------- - warp: `lsst.afw.image._exposure.ExposureF` - Warped exposure or difference image for which to estimate - background. + exposure: `lsst.afw.image._exposure.ExposureF` + Exposure to scale. Returns ------- - bkgd: `lsst.afw.math.BackgroundMI` - Background model container. - bctrl: `lsst.afw.math.BackgroundControl` - Background model control - warpMI: `lsst.afw.image._maskedImage.MaskedImageF` - Masked image associated with warp - statsFlag: `lsst.afw.math.Property` - Flag for grid statistic property (self.config.gridStatistic) + maskedImage: `lsst.afw.image._maskedImage.MaskedImageF` + Flux-scaled masked exposure. """ - statsFlag = stringToStatisticsProperty(self.config.gridStatistic) - self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) - self.statsCtrl.setNumIter(self.config.numIter) - - warpMI = warp.getMaskedImage() - - width = warpMI.getWidth() - height = warpMI.getHeight() - nx = width // self.config.binSize - if width % self.config.binSize != 0: - nx += 1 - ny = height // self.config.binSize - if height % self.config.binSize != 0: - ny += 1 - - bctrl = BackgroundControl(nx, ny, self.statsCtrl, statsFlag) - bctrl.setUndersampleStyle(self.config.undersampleStyle) - - bkgd = makeBackground(warpMI, bctrl) + maskedImage = exposure.getMaskedImage() + fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) + exposure.image.array *= fluxZp - return bkgd, bctrl, warpMI, statsFlag + return maskedImage @timeMethod def matchBackgrounds(self, refExposure, sciExposure): From 836ac7977d8989d869303a01b4b59ce68ef4e837 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Mon, 30 Sep 2024 02:09:48 -0700 Subject: [PATCH 08/13] Add background-matched image data type Otherwise, changes are clean-up from previous refactoring to restore functionality, plus a bug fix. Bug fix was the restoration of two lines of code in MatchBackgroundsTask.matchBackgrounds() which produced a difference image to work from. --- python/lsst/pipe/tasks/matchBackgrounds.py | 33 +++++++++++++++++----- 1 file changed, 26 insertions(+), 7 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 117c8467a..dfc75224b 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -76,6 +76,13 @@ class MatchBackgroundsConnections( storageClass="Background", multiple=True, ) + matchedImageList = Output( + doc="List of background-matched warps.", + name="{inputCoaddName}Coadd_{warpType}Warp_bgMatched", + storageClass="ExposureF", + dimensions=("skymap", "tract", "patch", "visit"), + multiple=True, + ) class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgroundsConnections): @@ -312,6 +319,7 @@ def run(self, warps): ) backgroundInfoList = [] + matchedImageList = [] for ind, exp in enumerate(warps): # TODO: simplify this maybe, using only visit IDs? self.log.info("Matching background of %s to %s", exp.dataId, refWarp.dataId) @@ -333,10 +341,13 @@ def run(self, warps): backgroundInfoStruct = blank backgroundInfoList.append(backgroundInfoStruct) + matchedImageList.append(toMatchExposure) # TODO: more elegant solution than inserting blank model at ref ind? backgroundInfoList.insert(refInd, blank) - return Struct(backgroundInfoList=backgroundInfoList) + matchedImageList.insert(refInd, refWarp.get()) + return Struct(backgroundInfoList=backgroundInfoList, + matchedImageList=matchedImageList) @timeMethod def _defineWarps(self, warps, refWarpVisit=None): @@ -439,7 +450,7 @@ def _makeBackground(self, warp: ExposureF) -> tuple[BackgroundMI, BackgroundCont Returns ------- warpBgMI: `~lsst.afw.math.BackgroundMI` - Background-subtracted masked image. + Background model of masked warp. bgCtrl: `~lsst.afw.math.BackgroundControl` Background control object. """ @@ -448,7 +459,11 @@ def _makeBackground(self, warp: ExposureF) -> tuple[BackgroundMI, BackgroundCont bgCtrl = BackgroundControl(nx, ny, self.statsCtrl, self.statsFlag) bgCtrl.setUndersampleStyle(self.config.undersampleStyle) - warpBgMI = makeBackground(warp.getMaskedImage(), bgCtrl) + # Difference image not in ExposureF format! And no reason it should be. + try: + warpBgMI = makeBackground(warp.getMaskedImage(), bgCtrl) + except AttributeError: + warpBgMI = makeBackground(warp, bgCtrl) return warpBgMI, bgCtrl @@ -541,14 +556,18 @@ def matchBackgrounds(self, refExposure, sciExposure): "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" % (wSci, hSci, wRef, hRef) ) - bkgd, bctrl, diffMI, statsFlag = self._setupBackground(refExposure) + im = refExposure.getMaskedImage() + diffMI = im.clone() + diffMI -= sciExposure.getMaskedImage() + + bkgd, bctrl = self._makeBackground(diffMI) # Some config and input checks if config.usePolynomial: # 1) Check that order/bin size make sense: # 2) Change binsize or order if underconstrained. if self.config.usePolynomial: order = self.config.order - bgX, bgY, bgZ, bgdZ = self._gridImage(diffMI, self.config.binSize, statsFlag) + bgX, bgY, bgZ, bgdZ = self._gridImage(diffMI, self.config.binSize, self.statsFlag) minNumberGridPoints = min(len(set(bgX)), len(set(bgY))) if len(bgZ) == 0: raise ValueError("No overlap with reference. Nothing to match") @@ -588,11 +607,11 @@ def matchBackgrounds(self, refExposure, sciExposure): sciMI = sciExposure.getMaskedImage() sciMI += bkgdImage - del sciMI + del sciMI # sciExposure is now a BG-matched image # Need RMS from fit: 2895 will replace this: rms = 0.0 - bgX, bgY, bgZ, bgdZ = self._gridImage(diffMI, self.config.binSize, statsFlag) + bgX, bgY, bgZ, bgdZ = self._gridImage(diffMI, self.config.binSize, self.statsFlag) x0, y0 = diffMI.getXY0() modelValueArr = np.empty(len(bgZ)) for i in range(len(bgX)): From 4add9cc9e68974c7c8b930ef7f91cd67c8d78906 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Wed, 2 Oct 2024 04:28:02 -0700 Subject: [PATCH 09/13] Refactor and bug-fix, AEW All images and background models now returned in counts, not nJy. --- python/lsst/pipe/tasks/matchBackgrounds.py | 209 +++++++-------------- 1 file changed, 72 insertions(+), 137 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index dfc75224b..a7f352959 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -143,7 +143,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr gridStatistic = ChoiceField( dtype=str, doc="Type of statistic to estimate pixel value for the grid points.", - default="MEAN", + default="MEANCLIP", allowed={"MEAN": "mean", "MEDIAN": "median", "MEANCLIP": "clipped mean"}, ) undersampleStyle = ChoiceField( @@ -180,7 +180,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr ) numIter = Field[int]( doc="Number of iterations of outlier rejection. Ignored if ``gridStatistic != 'MEANCLIP'``.", - default=2, + default=3, ) approxWeighting = Field[bool]( @@ -206,13 +206,15 @@ class MatchBackgroundsTask(PipelineTask): It matches the backgrounds of a list of science exposures to a reference science exposure. The reference exposure is chosen from the list of science exposures by - minimizing a cost function that penalizes high variance, high background - level, and low coverage. - The cost function is a weighted sum of these three metrics. + minimizing a cost function that penalizes high background complexity + (divergence from a plane), high variance, low global coverage, and low + coverage along image edges. + The cost function is a weighted sum of these four metrics. The weights are set by the config parameters: - - ``bestRefWeightCoverage`` + - ``bestRefWeightChi2`` - ``bestRefWeightVariance`` - - ``bestRefWeightLevel`` + - ``bestRefWeightGlobalCoverage`` + - ``bestRefWeightEdgeCoverage`` Attributes @@ -233,6 +235,7 @@ def __init__(self, *args, **kwargs): self.statsCtrl = StatisticsControl() # TODO: Check that setting the mask planes here work - these planes # can vary from exposure to exposure, I think? + # Aaron: I think only the bit values vary, not the names, which this is referencing. self.statsCtrl.setAndMask(Mask.getPlaneBitMask(self.config.badMaskPlanes)) self.statsCtrl.setNanSafe(True) self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) @@ -276,10 +279,7 @@ def run(self, warps): # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this refExposure = refWarp.get() - refMI = self._fluxScale(refExposure) # Also modifies refExposure - - # TODO: figure out what this was creating and why - # debugIdKeyList = tuple(set(expKeyList) - set(['tract', 'patch'])) + instFluxToNanojanskyRef = self._fluxScale(refExposure) self.log.info("Matching %d Exposures", numExp) @@ -290,7 +290,7 @@ def run(self, warps): # TODO: refactor below to construct blank bg model im = refExposure.getMaskedImage() - blankIm = im.Factory(im, True) + blankIm = im.clone() blankIm.image.array *= 0 width = blankIm.getWidth() @@ -320,17 +320,16 @@ def run(self, warps): backgroundInfoList = [] matchedImageList = [] - for ind, exp in enumerate(warps): - # TODO: simplify this maybe, using only visit IDs? - self.log.info("Matching background of %s to %s", exp.dataId, refWarp.dataId) + for exp in warps: + # TODO: simplify what this prints? + self.log.info( + "Matching background of %s to %s", + exp.dataId, + refWarp.dataId, + ) toMatchExposure = exp.get() + instFluxToNanojansky = self._fluxScale(toMatchExposure) try: - # TODO: adjust logic to avoid creating spurious MIs like this - toMatchMI = self._fluxScale(toMatchExposure) - - # store a string specifying the visit to label debug plot - # self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList]) - backgroundInfoStruct = self.matchBackgrounds( refExposure=refExposure, sciExposure=toMatchExposure, @@ -341,13 +340,14 @@ def run(self, warps): backgroundInfoStruct = blank backgroundInfoList.append(backgroundInfoStruct) + toMatchExposure.image /= instFluxToNanojansky # Back to cts matchedImageList.append(toMatchExposure) # TODO: more elegant solution than inserting blank model at ref ind? backgroundInfoList.insert(refInd, blank) - matchedImageList.insert(refInd, refWarp.get()) - return Struct(backgroundInfoList=backgroundInfoList, - matchedImageList=matchedImageList) + refExposure.image /= instFluxToNanojanskyRef # Back to cts + matchedImageList.insert(refInd, refExposure) + return Struct(backgroundInfoList=backgroundInfoList, matchedImageList=matchedImageList) @timeMethod def _defineWarps(self, warps, refWarpVisit=None): @@ -355,8 +355,9 @@ def _defineWarps(self, warps, refWarpVisit=None): If no reference warp ID is supplied, this method calculates an appropriate reference exposure from the supplied list of warps by - minimizing a cost function that penalizes high variance, high - background level, and low coverage. + minimizing a cost function that penalizes high background complexity + (divergence from a plane), high variance, low global coverage, and low + edge coverage. Parameters ---------- @@ -399,7 +400,7 @@ def _defineWarps(self, warps, refWarpVisit=None): instFluxToNanojansky = warp.getPhotoCalib().instFluxToNanojansky(1) warp.image *= instFluxToNanojansky # Images in nJy to facilitate difference imaging warp.variance *= instFluxToNanojansky - warpBg, _ = self._makeBackground(warp) + warpBg, _ = self._makeBackground(warp.getMaskedImage()) # Return an approximation to the background approxCtrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, self.config.approxWeighting) @@ -407,8 +408,11 @@ def _defineWarps(self, warps, refWarpVisit=None): warpBgSub = ImageF(warp.image.array - warpApprox.getImage().array) warpStats = makeStatistics(warpBgSub, warp.mask, VARIANCE | NPOINT, self.statsCtrl) - # TODO: need to modify this to account for the background mask - warpChi2 = np.nansum(warpBgSub.array**2 / warp.variance.array) + + bad_mask_bit_mask = warp.mask.getPlaneBitMask(self.config.badMaskPlanes) + good = (warp.mask.array.astype(int) & bad_mask_bit_mask) == 0 + dof = len(good[good]) - 6 # Assuming eq. of plane + warpChi2 = np.nansum(warpBgSub.array[good] ** 2 / warp.variance.array[good]) / dof warpVar, _ = warpStats.getResult(VARIANCE) warpNPointGlobal, _ = warpStats.getResult(NPOINT) warpNPointEdge = ( @@ -422,6 +426,15 @@ def _defineWarps(self, warps, refWarpVisit=None): warpNPointsGlobal.append(int(warpNPointGlobal)) warpNPointsEdge.append(warpNPointEdge) + self.log.info( + "Sci exp. visit %d; BG fit Chi^2=%.1f, var=%.1f nJy, nPoints global=%d, nPoints edge=%d", + warp.getInfo().getVisitInfo().id, + warpChi2, + warpVar, + warpNPointGlobal, + warpNPointEdge, + ) + # Normalize mean/var/npoints to range from 0 to 1 warpChi2sFrac = np.array(warpChi2s) / np.nanmax(warpChi2s) warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) @@ -439,17 +452,17 @@ def _defineWarps(self, warps, refWarpVisit=None): self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) return refWarp, ind - def _makeBackground(self, warp: ExposureF) -> tuple[BackgroundMI, BackgroundControl]: + def _makeBackground(self, warp: MaskedImageF) -> tuple[BackgroundMI, BackgroundControl]: """Generate a simple binned background masked image for warped data. Parameters ---------- - warp: `~lsst.afw.image.ExposureF` + warp: `~lsst.afw.image.MaskedImageF` Warped exposure for which to estimate background. Returns ------- - warpBgMI: `~lsst.afw.math.BackgroundMI` + bkgd: `~lsst.afw.math.BackgroundMI` Background model of masked warp. bgCtrl: `~lsst.afw.math.BackgroundControl` Background control object. @@ -459,13 +472,9 @@ def _makeBackground(self, warp: ExposureF) -> tuple[BackgroundMI, BackgroundCont bgCtrl = BackgroundControl(nx, ny, self.statsCtrl, self.statsFlag) bgCtrl.setUndersampleStyle(self.config.undersampleStyle) - # Difference image not in ExposureF format! And no reason it should be. - try: - warpBgMI = makeBackground(warp.getMaskedImage(), bgCtrl) - except AttributeError: - warpBgMI = makeBackground(warp, bgCtrl) + bkgd = makeBackground(warp, bgCtrl) - return warpBgMI, bgCtrl + return bkgd, bgCtrl def _fluxScale(self, exposure): """Scales image to nJy flux using photometric calibration. @@ -477,14 +486,12 @@ def _fluxScale(self, exposure): Returns ------- - maskedImage: `lsst.afw.image._maskedImage.MaskedImageF` - Flux-scaled masked exposure. + fluxZp: `float` + Counts to nanojanskies conversion factor """ - maskedImage = exposure.getMaskedImage() fluxZp = exposure.getPhotoCalib().instFluxToNanojansky(1) - exposure.image.array *= fluxZp - - return maskedImage + exposure.image *= fluxZp + return fluxZp @timeMethod def matchBackgrounds(self, refExposure, sciExposure): @@ -509,26 +516,15 @@ def matchBackgrounds(self, refExposure, sciExposure): refExposure : `lsst.afw.image.Exposure` Reference exposure. sciExposure : `lsst.afw.image.Exposure` - Science exposure; modified by changing the background level - to match that of the reference exposure. + Science exposure; ultimately modified by changing the background + level to match that of the reference exposure. Returns ------- - model : `lsst.pipe.base.Struct` - Background model as a struct with attributes: - - ``backgroundModel`` - An afw.math.Approximate or an afw.math.Background. - ``fitRMS`` - RMS of the fit = sqrt(mean(residuals**2)), (`float`). - ``matchedMSE`` - The MSE of the reference and matched images: - mean((refImage - matchedSciImage)**2); - Should be comparable to difference image's mean variance - (`float`). - ``diffImVar`` - The mean variance of the difference image (`float`). + model : `~lsst.afw.math.BackgroundMI` + Background model of difference image, reference - science """ + # TODO: this is deprecated if lsstDebug.Info(__name__).savefits: refExposure.writeFits(lsstDebug.Info(__name__).figpath + "refExposure.fits") sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciExposure.fits") @@ -545,6 +541,7 @@ def matchBackgrounds(self, refExposure, sciExposure): if shortSideLength % self.config.binSize != 0: npoints += 1 + # If order of polynomial to be fit > number of bins to fit, error if self.config.order > npoints - 1: raise ValueError("%d = config.order > npoints - 1 = %d" % (self.config.order, npoints - 1)) @@ -572,7 +569,7 @@ def matchBackgrounds(self, refExposure, sciExposure): if len(bgZ) == 0: raise ValueError("No overlap with reference. Nothing to match") elif minNumberGridPoints <= self.config.order: - # must either lower order or raise number of bins or throw exception + # must lower order or raise number of bins, or throw exception if self.config.undersampleStyle == "THROW_EXCEPTION": raise ValueError("Image does not cover enough of ref image for order and binsize") elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER": @@ -585,7 +582,8 @@ def matchBackgrounds(self, refExposure, sciExposure): bkgd = makeBackground(diffMI, bctrl) # do over self.log.warning("Decreasing binsize to %d", newBinSize) - # If there is no variance in any image pixels, do not weight bins by inverse variance + # If there is no variance in any image pixels, + # do not weight bins by inverse variance isUniformImageDiff = not np.any(bgdZ > self.config.gridStdevEpsilon) weightByInverseVariance = False if isUniformImageDiff else self.config.approxWeighting @@ -602,9 +600,10 @@ def matchBackgrounds(self, refExposure, sciExposure): bkgdImage = bkgd.getImageF(self.config.interpStyle, self.config.undersampleStyle) except Exception as e: raise RuntimeError( - "Background/Approximation failed to interp image %s: %s" % (self.debugDataIdString, e) + "Background/Approximation failed to interp image %s: %s" % (sciExposure.dataId, e) ) + instFluxToNanojansky = sciExposure.getPhotoCalib().instFluxToNanojansky(1) sciMI = sciExposure.getMaskedImage() sciMI += bkgdImage del sciMI # sciExposure is now a BG-matched image @@ -619,6 +618,7 @@ def matchBackgrounds(self, refExposure, sciExposure): resids = bgZ - modelValueArr rms = np.sqrt(np.mean(resids[~np.isnan(resids)] ** 2)) + # TODO: also deprecated; _gridImage() maybe can go? if lsstDebug.Info(__name__).savefits: sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciMatchedExposure.fits") @@ -637,6 +637,12 @@ def matchBackgrounds(self, refExposure, sciExposure): mse = makeStatistics(diffMI, MEANSQUARE, self.statsCtrl).getValue() outBkgd = approx if self.config.usePolynomial else bkgd + # Convert this back into counts + # TODO: is there a one-line way to do this? + statsIm = outBkgd.getStatsImage() + statsIm /= instFluxToNanojansky + bkgdIm = outBkgd.getImageF() + bkgdIm /= instFluxToNanojansky self.log.info( "Visit %d; difference BG fit RMS=%.1f cts, matched MSE=%.1f cts, mean variance=%.1f cts", @@ -659,7 +665,9 @@ def matchBackgrounds(self, refExposure, sciExposure): ) def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): - """Generate a plot showing the background fit and residuals. + """ + Consider deleting this entirely + Generate a plot showing the background fit and residuals. It is called when lsstDebug.Info(__name__).savefig = True. Saves the fig to lsstDebug.Info(__name__).figpath. @@ -771,76 +779,3 @@ def _gridImage(self, maskedImage, binsize, statsFlag): bgZ.append(est) return np.array(bgX), np.array(bgY), np.array(bgZ), np.array(bgdZ) - - -class DataRefMatcher: - """Match data references for a specified dataset type. - - Note that this is not exact, but should suffice for this task - until there is better support for this kind of thing in the butler. - - Parameters - ---------- - butler : `lsst.daf.butler.Butler` - Butler to search for maches in. - datasetType : `str` - Dataset type to match. - """ - - def __init__(self, butler, datasetType): - self._datasetType = datasetType # for diagnostics - self._keyNames = butler.getKeys(datasetType) - - def _makeKey(self, ref): - """Return a tuple of values for the specified keyNames. - - Parameters - ---------- - ref : `Unknown` - Data reference. - - Raises - ------ - KeyError - Raised if ref.dataId is missing a key in keyNames. - """ - return tuple(ref.dataId[key] for key in self._keyNames) - - def isMatch(self, ref0, ref1): - """Return True if ref0 == ref1. - - Parameters - ---------- - ref0 : `Unknown` - Data for ref 0. - ref1 : `Unknown` - Data for ref 1. - - Raises - ------ - KeyError - Raised if either ID is missing a key in keyNames. - """ - return self._makeKey(ref0) == self._makeKey(ref1) - - def matchList(self, ref0, refList): - """Return a list of indices of matches. - - Parameters - ---------- - ref0 : `Unknown` - Data for ref 0. - `refList` : `list` - - Returns - ------- - matches : `tuple` - Tuple of indices of matches. - - Raises - ------ - KeyError - Raised if any ID is missing a key in keyNames. - """ - key0 = self._makeKey(ref0) - return tuple(ind for ind, ref in enumerate(refList) if self._makeKey(ref) == key0) From f7405ade38f69ca8c8eb776c04f8daf55d869dbb Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Mon, 7 Oct 2024 07:18:17 -0700 Subject: [PATCH 10/13] Add config for reference image selection bin size --- python/lsst/pipe/tasks/matchBackgrounds.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index a7f352959..d2e0cf2be 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -104,7 +104,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr dtype=float, doc="Image variance weight when calculating the best reference exposure. " "Higher weights prefers exposures with low image variances. Ignored when ref visit supplied", - default=0.4, + default=0.2, min=0.0, max=1.0, ) @@ -120,7 +120,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr dtype=float, doc="Edge coverage weight (number of valid edge pixels) when calculating the best reference " "exposure. Higher weights prefer exposures with high coverage. Ignored when a ref visit supplied.", - default=0.2, + default=0.4, min=0.0, max=1.0, ) @@ -159,6 +159,10 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr ) binSize = Field[int]( doc="Bin size for gridding the difference image and fitting a spatial model.", + default=256, + ) + chi2BinSize = Field[int]( + doc="Bin size for gridding images when choosing best reference exposure.", default=1024, ) interpStyle = ChoiceField( @@ -400,7 +404,7 @@ def _defineWarps(self, warps, refWarpVisit=None): instFluxToNanojansky = warp.getPhotoCalib().instFluxToNanojansky(1) warp.image *= instFluxToNanojansky # Images in nJy to facilitate difference imaging warp.variance *= instFluxToNanojansky - warpBg, _ = self._makeBackground(warp.getMaskedImage()) + warpBg, _ = self._makeBackground(warp.getMaskedImage(), binSize=self.config.chi2BinSize) # Return an approximation to the background approxCtrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, self.config.approxWeighting) @@ -452,7 +456,7 @@ def _defineWarps(self, warps, refWarpVisit=None): self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) return refWarp, ind - def _makeBackground(self, warp: MaskedImageF) -> tuple[BackgroundMI, BackgroundControl]: + def _makeBackground(self, warp: MaskedImageF, binSize) -> tuple[BackgroundMI, BackgroundControl]: """Generate a simple binned background masked image for warped data. Parameters @@ -467,8 +471,8 @@ def _makeBackground(self, warp: MaskedImageF) -> tuple[BackgroundMI, BackgroundC bgCtrl: `~lsst.afw.math.BackgroundControl` Background control object. """ - nx = warp.getWidth() // self.config.binSize - ny = warp.getHeight() // self.config.binSize + nx = warp.getWidth() // binSize + ny = warp.getHeight() // binSize bgCtrl = BackgroundControl(nx, ny, self.statsCtrl, self.statsFlag) bgCtrl.setUndersampleStyle(self.config.undersampleStyle) @@ -557,7 +561,7 @@ def matchBackgrounds(self, refExposure, sciExposure): diffMI = im.clone() diffMI -= sciExposure.getMaskedImage() - bkgd, bctrl = self._makeBackground(diffMI) + bkgd, bctrl = self._makeBackground(diffMI, binSize=self.config.binSize) # Some config and input checks if config.usePolynomial: # 1) Check that order/bin size make sense: From 6050e87a0af9418518c43be821e7f626024f44bd Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Wed, 20 Nov 2024 01:33:50 -0800 Subject: [PATCH 11/13] Refactor AEW --- python/lsst/pipe/tasks/matchBackgrounds.py | 190 +++------------------ 1 file changed, 20 insertions(+), 170 deletions(-) diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index d2e0cf2be..8906dba5c 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -21,16 +21,12 @@ __all__ = ["MatchBackgroundsConnections", "MatchBackgroundsConfig", "MatchBackgroundsTask"] -import lsstDebug import numpy as np -from lsst.afw.image import LOCAL, PARENT, ExposureF, ImageF, Mask, MaskedImageF +from lsst.afw.image import LOCAL, ImageF, Mask, MaskedImageF from lsst.afw.math import ( - MEAN, MEANCLIP, MEANSQUARE, - MEDIAN, NPOINT, - STDEV, VARIANCE, ApproximateControl, BackgroundControl, @@ -43,7 +39,6 @@ stringToStatisticsProperty, stringToUndersampleStyle, ) -from lsst.geom import Box2D, Box2I, PointI from lsst.pex.config import ChoiceField, Field, ListField, RangeField from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct, TaskError from lsst.pipe.base.connectionTypes import Input, Output @@ -137,7 +132,16 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr ) badMaskPlanes = ListField[str]( doc="Names of mask planes to ignore while estimating the background.", - default=["NO_DATA", "DETECTED", "DETECTED_NEGATIVE", "SAT", "BAD", "INTRP", "CR"], + default=[ + "NO_DATA", + "DETECTED", + "DETECTED_NEGATIVE", + "SAT", + "BAD", + "INTRP", + "CR", + "NOT_DEBLENDED", + ], itemCheck=lambda x: x in Mask().getMaskPlaneDict(), ) gridStatistic = ChoiceField( @@ -237,9 +241,6 @@ def __init__(self, *args, **kwargs): super().__init__(**kwargs) self.statsFlag = stringToStatisticsProperty(self.config.gridStatistic) self.statsCtrl = StatisticsControl() - # TODO: Check that setting the mask planes here work - these planes - # can vary from exposure to exposure, I think? - # Aaron: I think only the bit values vary, not the names, which this is referencing. self.statsCtrl.setAndMask(Mask.getPlaneBitMask(self.config.badMaskPlanes)) self.statsCtrl.setNanSafe(True) self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) @@ -278,7 +279,7 @@ def run(self, warps): raise TaskError("No exposures to match") # Define a reference warp; 'warps' is modified in-place to exclude it - refWarp, refInd = self._defineWarps(warps=warps, refWarpVisit=self.config.refWarpVisit) + refWarp, refInd, bkgd = self._defineWarps(warps=warps, refWarpVisit=self.config.refWarpVisit) # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this @@ -287,29 +288,11 @@ def run(self, warps): self.log.info("Matching %d Exposures", numExp) - # Creating a null BackgroundList object by fitting a blank image - statsFlag = stringToStatisticsProperty(self.config.gridStatistic) - self.statsCtrl.setNumSigmaClip(self.config.numSigmaClip) - self.statsCtrl.setNumIter(self.config.numIter) - - # TODO: refactor below to construct blank bg model - im = refExposure.getMaskedImage() - blankIm = im.clone() - blankIm.image.array *= 0 - - width = blankIm.getWidth() - height = blankIm.getHeight() - nx = width // self.config.binSize - if width % self.config.binSize != 0: - nx += 1 - ny = height // self.config.binSize - if height % self.config.binSize != 0: - ny += 1 - - bctrl = BackgroundControl(nx, ny, self.statsCtrl, statsFlag) - bctrl.setUndersampleStyle(self.config.undersampleStyle) - - bkgd = makeBackground(blankIm, bctrl) + # Blank ref warp background as reference background + bkgdIm = bkgd.getImageF() + bkgdStatsIm = bkgd.getStatsImage() + bkgdIm *= 0 + bkgdStatsIm *= 0 blank = BackgroundList( ( bkgd, @@ -325,7 +308,6 @@ def run(self, warps): backgroundInfoList = [] matchedImageList = [] for exp in warps: - # TODO: simplify what this prints? self.log.info( "Matching background of %s to %s", exp.dataId, @@ -347,7 +329,6 @@ def run(self, warps): toMatchExposure.image /= instFluxToNanojansky # Back to cts matchedImageList.append(toMatchExposure) - # TODO: more elegant solution than inserting blank model at ref ind? backgroundInfoList.insert(refInd, blank) refExposure.image /= instFluxToNanojanskyRef # Back to cts matchedImageList.insert(refInd, refExposure) @@ -377,6 +358,8 @@ def _defineWarps(self, warps, refWarpVisit=None): Reference warped exposure. refWarpIndex : `int` Index of the reference removed from the list of warps. + warpBg : `~lsst.afw.math.BackgroundMI` + Temporary background model, used to make a blank BG for the ref Notes ----- @@ -454,7 +437,7 @@ def _defineWarps(self, warps, refWarpVisit=None): ind = np.nanargmin(costFunctionVals) refWarp = warps.pop(ind) self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) - return refWarp, ind + return refWarp, ind, warpBg def _makeBackground(self, warp: MaskedImageF, binSize) -> tuple[BackgroundMI, BackgroundControl]: """Generate a simple binned background masked image for warped data. @@ -528,11 +511,6 @@ def matchBackgrounds(self, refExposure, sciExposure): model : `~lsst.afw.math.BackgroundMI` Background model of difference image, reference - science """ - # TODO: this is deprecated - if lsstDebug.Info(__name__).savefits: - refExposure.writeFits(lsstDebug.Info(__name__).figpath + "refExposure.fits") - sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciExposure.fits") - # Check Configs for polynomials: if self.config.usePolynomial: x, y = sciExposure.getDimensions() @@ -622,17 +600,6 @@ def matchBackgrounds(self, refExposure, sciExposure): resids = bgZ - modelValueArr rms = np.sqrt(np.mean(resids[~np.isnan(resids)] ** 2)) - # TODO: also deprecated; _gridImage() maybe can go? - if lsstDebug.Info(__name__).savefits: - sciExposure.writeFits(lsstDebug.Info(__name__).figpath + "sciMatchedExposure.fits") - - if lsstDebug.Info(__name__).savefig: - bbox = Box2D(refExposure.getMaskedImage().getBBox()) - try: - self._debugPlot(bgX, bgY, bgZ, bgdZ, bkgdImage, bbox, modelValueArr, resids) - except Exception as e: - self.log.warning("Debug plot not generated: %s", e) - meanVar = makeStatistics(diffMI.getVariance(), diffMI.getMask(), MEANCLIP, self.statsCtrl).getValue() diffIm = diffMI.getImage() @@ -642,7 +609,6 @@ def matchBackgrounds(self, refExposure, sciExposure): outBkgd = approx if self.config.usePolynomial else bkgd # Convert this back into counts - # TODO: is there a one-line way to do this? statsIm = outBkgd.getStatsImage() statsIm /= instFluxToNanojansky bkgdIm = outBkgd.getImageF() @@ -667,119 +633,3 @@ def matchBackgrounds(self, refExposure, sciExposure): False, ) ) - - def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids): - """ - Consider deleting this entirely - Generate a plot showing the background fit and residuals. - - It is called when lsstDebug.Info(__name__).savefig = True. - Saves the fig to lsstDebug.Info(__name__).figpath. - Displays on screen if lsstDebug.Info(__name__).display = True. - - Parameters - ---------- - X : `np.ndarray`, (N,) - Array of x positions. - Y : `np.ndarray`, (N,) - Array of y positions. - Z : `np.ndarray` - Array of the grid values that were interpolated. - dZ : `np.ndarray`, (len(Z),) - Array of the error on the grid values. - modelImage : `Unknown` - Image of the model of the fit. - model : `np.ndarray`, (len(Z),) - Array of len(Z) containing the grid values predicted by the model. - resids : `Unknown` - Z - model. - """ - import matplotlib.colors - import matplotlib.pyplot as plt - from mpl_toolkits.axes_grid1 import ImageGrid - - zeroIm = MaskedImageF(Box2I(bbox)) - zeroIm += modelImage - x0, y0 = zeroIm.getXY0() - dx, dy = zeroIm.getDimensions() - if len(Z) == 0: - self.log.warning("No grid. Skipping plot generation.") - else: - max, min = np.max(Z), np.min(Z) - norm = matplotlib.colors.normalize(vmax=max, vmin=min) - maxdiff = np.max(np.abs(resids)) - diffnorm = matplotlib.colors.normalize(vmax=maxdiff, vmin=-maxdiff) - rms = np.sqrt(np.mean(resids**2)) - fig = plt.figure(1, (8, 6)) - meanDz = np.mean(dZ) - grid = ImageGrid( - fig, - 111, - nrows_ncols=(1, 2), - axes_pad=0.1, - share_all=True, - label_mode="L", - cbar_mode="each", - cbar_size="7%", - cbar_pad="2%", - cbar_location="top", - ) - im = grid[0].imshow( - zeroIm.getImage().getArray(), extent=(x0, x0 + dx, y0 + dy, y0), norm=norm, cmap="Spectral" - ) - im = grid[0].scatter( - X, Y, c=Z, s=15.0 * meanDz / dZ, edgecolor="none", norm=norm, marker="o", cmap="Spectral" - ) - im2 = grid[1].scatter(X, Y, c=resids, edgecolor="none", norm=diffnorm, marker="s", cmap="seismic") - grid.cbar_axes[0].colorbar(im) - grid.cbar_axes[1].colorbar(im2) - grid[0].axis([x0, x0 + dx, y0 + dy, y0]) - grid[1].axis([x0, x0 + dx, y0 + dy, y0]) - grid[0].set_xlabel("model and grid") - grid[1].set_xlabel("residuals. rms = %0.3f" % (rms)) - if lsstDebug.Info(__name__).savefig: - fig.savefig(lsstDebug.Info(__name__).figpath + self.debugDataIdString + ".png") - if lsstDebug.Info(__name__).display: - plt.show() - plt.clf() - - def _gridImage(self, maskedImage, binsize, statsFlag): - """Private method to grid an image for debugging.""" - width, height = maskedImage.getDimensions() - x0, y0 = maskedImage.getXY0() - xedges = np.arange(0, width, binsize) - yedges = np.arange(0, height, binsize) - xedges = np.hstack((xedges, width)) # add final edge - yedges = np.hstack((yedges, height)) # add final edge - - # Use lists/append to protect against the case where - # a bin has no valid pixels and should not be included in the fit - bgX = [] - bgY = [] - bgZ = [] - bgdZ = [] - - for ymin, ymax in zip(yedges[0:-1], yedges[1:]): - for xmin, xmax in zip(xedges[0:-1], xedges[1:]): - subBBox = Box2I( - PointI(int(x0 + xmin), int(y0 + ymin)), - PointI(int(x0 + xmax - 1), int(y0 + ymax - 1)), - ) - subIm = MaskedImageF(maskedImage, subBBox, PARENT, False) - stats = makeStatistics( - subIm, - MEAN | MEANCLIP | MEDIAN | NPOINT | STDEV, - self.statsCtrl, - ) - npoints, _ = stats.getResult(NPOINT) - if npoints >= 2: - stdev, _ = stats.getResult(STDEV) - if stdev < self.config.gridStdevEpsilon: - stdev = self.config.gridStdevEpsilon - bgX.append(0.5 * (x0 + xmin + x0 + xmax)) - bgY.append(0.5 * (y0 + ymin + y0 + ymax)) - bgdZ.append(stdev / np.sqrt(npoints)) - est, _ = stats.getResult(statsFlag) - bgZ.append(est) - - return np.array(bgX), np.array(bgY), np.array(bgZ), np.array(bgdZ) From 27b7f908f975c35cd3d67ffb135561671021f4e3 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Mon, 27 Jan 2025 05:15:18 -0800 Subject: [PATCH 12/13] Add full tract background functionality. `matchBackgrounds` in its original form matches by warps, i.e. single detectors. A more sensible thing to do is to match backgrounds across the whole focal plane. This functionality needed to be added to do this using warped exposures, in tract coordinates, so this is now added to `backgrounds.py`. This commit also includes a partial revision to `matchBackgrounds.py` using this new functionality, choosing a reference visit ID based on these background models instead of on individual warps. Full functionality has yet to be restored in light of these changes. --- python/lsst/pipe/tasks/background.py | 533 +++++++++++++++++---- python/lsst/pipe/tasks/matchBackgrounds.py | 240 +++++++--- 2 files changed, 604 insertions(+), 169 deletions(-) diff --git a/python/lsst/pipe/tasks/background.py b/python/lsst/pipe/tasks/background.py index 532a2fba3..ae90d4167 100644 --- a/python/lsst/pipe/tasks/background.py +++ b/python/lsst/pipe/tasks/background.py @@ -30,22 +30,21 @@ "SkyStatsConfig", ] -import sys -import numpy import importlib import itertools -from scipy.ndimage import gaussian_filter +import sys -import lsst.afw.math as afwMath -import lsst.afw.image as afwImage -import lsst.afw.geom as afwGeom import lsst.afw.cameraGeom as afwCameraGeom +import lsst.afw.geom as afwGeom +import lsst.afw.image as afwImage +import lsst.afw.math as afwMath +import lsst.afw.table as afwTable import lsst.geom as geom import lsst.meas.algorithms as measAlg -import lsst.afw.table as afwTable - -from lsst.pex.config import Config, Field, ListField, ChoiceField, ConfigField, RangeField, ConfigurableField +import numpy +from lsst.pex.config import ChoiceField, Config, ConfigField, ConfigurableField, Field, ListField, RangeField from lsst.pipe.base import Task +from scipy.ndimage import gaussian_filter def robustMean(array, rej=3.0): @@ -64,47 +63,62 @@ def robustMean(array, rej=3.0): Robust mean of `array`. """ q1, median, q3 = numpy.percentile(array, [25.0, 50.0, 100.0]) - good = numpy.abs(array - median) < rej*0.74*(q3 - q1) + good = numpy.abs(array - median) < rej * 0.74 * (q3 - q1) return array[good].mean() class BackgroundConfig(Config): """Configuration for background measurement""" - statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points", - allowed={"MEANCLIP": "clipped mean", - "MEAN": "unclipped mean", - "MEDIAN": "median"}) + + statistic = ChoiceField( + dtype=str, + default="MEANCLIP", + doc="type of statistic to use for grid points", + allowed={"MEANCLIP": "clipped mean", "MEAN": "unclipped mean", "MEDIAN": "median"}, + ) xBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in x") yBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in y") - algorithm = ChoiceField(dtype=str, default="NATURAL_SPLINE", optional=True, - doc="How to interpolate the background values. " - "This maps to an enum; see afw::math::Background", - allowed={ - "CONSTANT": "Use a single constant value", - "LINEAR": "Use linear interpolation", - "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints", - "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust" - " to outliers", - "NONE": "No background estimation is to be attempted", - }) - mask = ListField(dtype=str, default=["SAT", "BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA"], - doc="Names of mask planes to ignore while estimating the background") + algorithm = ChoiceField( + dtype=str, + default="NATURAL_SPLINE", + optional=True, + doc="How to interpolate the background values. " "This maps to an enum; see afw::math::Background", + allowed={ + "CONSTANT": "Use a single constant value", + "LINEAR": "Use linear interpolation", + "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints", + "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust" " to outliers", + "NONE": "No background estimation is to be attempted", + }, + ) + mask = ListField( + dtype=str, + default=["SAT", "BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA"], + doc="Names of mask planes to ignore while estimating the background", + ) class SkyStatsConfig(Config): """Parameters controlling the measurement of sky statistics""" - statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points", - allowed={"MEANCLIP": "clipped mean", - "MEAN": "unclipped mean", - "MEDIAN": "median"}) + + statistic = ChoiceField( + dtype=str, + default="MEANCLIP", + doc="type of statistic to use for grid points", + allowed={"MEANCLIP": "clipped mean", "MEAN": "unclipped mean", "MEDIAN": "median"}, + ) clip = Field(doc="Clipping threshold for background", dtype=float, default=3.0) nIter = Field(doc="Clipping iterations for background", dtype=int, default=3) - mask = ListField(doc="Mask planes to reject", dtype=str, - default=["SAT", "DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA"]) + mask = ListField( + doc="Mask planes to reject", + dtype=str, + default=["SAT", "DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA"], + ) class SkyMeasurementConfig(Config): """Configuration for SkyMeasurementTask""" + skyIter = Field(dtype=int, default=3, doc="k-sigma rejection iterations for sky scale") skyRej = Field(dtype=float, default=3.0, doc="k-sigma rejection threshold for sky scale") background = ConfigField(dtype=BackgroundConfig, doc="Background measurement") @@ -122,6 +136,7 @@ class SkyMeasurementTask(Task): background model (a `lsst.afw.math.BackgroundMI`). The sky frame represents the dominant response of the camera to the sky background. """ + ConfigClass = SkyMeasurementConfig @staticmethod @@ -149,11 +164,16 @@ def exposureToBackground(bgExp): algorithm = header.getScalar("ALGORITHM") bbox = geom.Box2I(geom.Point2I(xMin, yMin), geom.Point2I(xMax, yMax)) return afwMath.BackgroundList( - (afwMath.BackgroundMI(bbox, bgExp.getMaskedImage()), - afwMath.stringToInterpStyle(algorithm), - afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), - afwMath.ApproximateControl.UNKNOWN, - 0, 0, False)) + ( + afwMath.BackgroundMI(bbox, bgExp.getMaskedImage()), + afwMath.stringToInterpStyle(algorithm), + afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) def backgroundToExposure(self, statsImage, bbox): """Convert a background model to an exposure @@ -207,22 +227,26 @@ def measureBackground(self, image): stats.setNanSafe(True) ctrl = afwMath.BackgroundControl( self.config.background.algorithm, - max(int(image.getWidth()/self.config.background.xBinSize + 0.5), 1), - max(int(image.getHeight()/self.config.background.yBinSize + 0.5), 1), + max(int(image.getWidth() / self.config.background.xBinSize + 0.5), 1), + max(int(image.getHeight() / self.config.background.yBinSize + 0.5), 1), "REDUCE_INTERP_ORDER", stats, - self.config.background.statistic + self.config.background.statistic, ) bg = afwMath.makeBackground(image, ctrl) - return afwMath.BackgroundList(( - bg, - afwMath.stringToInterpStyle(self.config.background.algorithm), - afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), - afwMath.ApproximateControl.UNKNOWN, - 0, 0, False - )) + return afwMath.BackgroundList( + ( + bg, + afwMath.stringToInterpStyle(self.config.background.algorithm), + afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) def averageBackgrounds(self, bgList): """Average multiple background models @@ -243,8 +267,9 @@ def averageBackgrounds(self, bgList): assert all(len(bg) == 1 for bg in bgList), "Mixed bgList: %s" % ([len(bg) for bg in bgList],) images = [bg[0][0].getStatsImage() for bg in bgList] boxes = [bg[0][0].getImageBBox() for bg in bgList] - assert len(set((box.getMinX(), box.getMinY(), box.getMaxX(), box.getMaxY()) for box in boxes)) == 1, \ - "Bounding boxes not all equal" + assert ( + len(set((box.getMinX(), box.getMinY(), box.getMaxX(), box.getMaxY()) for box in boxes)) == 1 + ), "Bounding boxes not all equal" bbox = boxes.pop(0) # Ensure bad pixels are masked @@ -343,19 +368,25 @@ def solveScales(self, scales): skySamples = numpy.array(skySamples) def solve(mask): +<<<<<<< HEAD # Make sure we return a float, not an array. return afwMath.LeastSquares.fromDesignMatrix(skySamples[mask].reshape(mask.sum(), 1), imageSamples[mask], afwMath.LeastSquares.DIRECT_SVD).getSolution()[0] +======= + return afwMath.LeastSquares.fromDesignMatrix( + skySamples[mask].reshape(mask.sum(), 1), imageSamples[mask], afwMath.LeastSquares.DIRECT_SVD + ).getSolution() +>>>>>>> 8cf1e691 (Add full tract background functionality.) mask = numpy.isfinite(imageSamples) & numpy.isfinite(skySamples) for ii in range(self.config.skyIter): solution = solve(mask) - residuals = imageSamples - solution*skySamples + residuals = imageSamples - solution * skySamples lq, uq = numpy.percentile(residuals[mask], [25, 75]) - stdev = 0.741*(uq - lq) # Robust stdev from IQR + stdev = 0.741 * (uq - lq) # Robust stdev from IQR with numpy.errstate(invalid="ignore"): # suppress NAN warnings - bad = numpy.abs(residuals) > self.config.skyRej*stdev + bad = numpy.abs(residuals) > self.config.skyRej * stdev mask[bad] = False return solve(mask) @@ -415,14 +446,15 @@ def interpolate1D(method, xSample, ySample, xInterp): """ if len(xSample) == 0: - return numpy.ones_like(xInterp)*numpy.nan + return numpy.ones_like(xInterp) * numpy.nan try: - return afwMath.makeInterpolate(xSample.astype(float), ySample.astype(float), - method).interpolate(xInterp.astype(float)) + return afwMath.makeInterpolate(xSample.astype(float), ySample.astype(float), method).interpolate( + xInterp.astype(float) + ) except Exception: if method == afwMath.Interpolate.CONSTANT: # We've already tried the most basic interpolation and it failed - return numpy.ones_like(xInterp)*numpy.nan + return numpy.ones_like(xInterp) * numpy.nan newMethod = afwMath.lookupMaxInterpStyle(len(xSample)) if newMethod == method: newMethod = afwMath.Interpolate.CONSTANT @@ -454,15 +486,17 @@ def interpolateBadPixels(array, isBad, interpolationStyle): isGood = ~isBad for y in range(height): if numpy.any(isBad[y, :]) and numpy.any(isGood[y, :]): - array[y][isBad[y]] = interpolate1D(method, xIndices[isGood[y]], array[y][isGood[y]], - xIndices[isBad[y]]) + array[y][isBad[y]] = interpolate1D( + method, xIndices[isGood[y]], array[y][isGood[y]], xIndices[isBad[y]] + ) isBad = numpy.isnan(array) isGood = ~isBad for x in range(width): if numpy.any(isBad[:, x]) and numpy.any(isGood[:, x]): - array[:, x][isBad[:, x]] = interpolate1D(method, yIndices[isGood[:, x]], - array[:, x][isGood[:, x]], yIndices[isBad[:, x]]) + array[:, x][isBad[:, x]] = interpolate1D( + method, yIndices[isGood[:, x]], array[:, x][isGood[:, x]], yIndices[isBad[:, x]] + ) class FocalPlaneBackgroundConfig(Config): @@ -474,15 +508,21 @@ class FocalPlaneBackgroundConfig(Config): need to be revised according to each particular camera. For this reason, no defaults are set for those. """ + xSize = Field(dtype=float, doc="Bin size in x") ySize = Field(dtype=float, doc="Bin size in y") pixelSize = Field(dtype=float, default=1.0, doc="Pixel size in same units as xSize/ySize") minFrac = Field(dtype=float, default=0.1, doc="Minimum fraction of bin size for good measurement") - mask = ListField(dtype=str, doc="Mask planes to treat as bad", - default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"]) + mask = ListField( + dtype=str, + doc="Mask planes to treat as bad", + default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"], + ) interpolation = ChoiceField( doc="how to interpolate the background values. This maps to an enum; see afw::math::Background", - dtype=str, default="AKIMA_SPLINE", optional=True, + dtype=str, + default="AKIMA_SPLINE", + optional=True, allowed={ "CONSTANT": "Use a single constant value", "LINEAR": "Use linear interpolation", @@ -521,6 +561,7 @@ class FocalPlaneBackground: Once you've built the background model, you can apply it to individual CCDs with the `toCcdBackground` method. """ + @classmethod def fromCamera(cls, config, camera): """Construct from a camera object @@ -539,14 +580,17 @@ def fromCamera(cls, config, camera): width, height = cameraBox.getDimensions() # Offset so that we run from zero - offset = geom.Extent2D(cameraBox.getMin())*-1 + offset = geom.Extent2D(cameraBox.getMin()) * -1 # Add an extra pixel buffer on either side - dims = geom.Extent2I(int(numpy.ceil(width/config.xSize)) + 2, - int(numpy.ceil(height/config.ySize)) + 2) + dims = geom.Extent2I( + int(numpy.ceil(width / config.xSize)) + 2, int(numpy.ceil(height / config.ySize)) + 2 + ) # Transform takes us from focal plane coordinates --> sample coordinates - transform = (geom.AffineTransform.makeTranslation(geom.Extent2D(1, 1)) - * geom.AffineTransform.makeScaling(1.0/config.xSize, 1.0/config.ySize) - * geom.AffineTransform.makeTranslation(offset)) + transform = ( + geom.AffineTransform.makeTranslation(geom.Extent2D(1, 1)) + * geom.AffineTransform.makeScaling(1.0 / config.xSize, 1.0 / config.ySize) + * geom.AffineTransform.makeTranslation(offset) + ) return cls(config, dims, afwGeom.makeTransform(transform)) @@ -627,8 +671,9 @@ def addCcd(self, exposure): CCD exposure to measure """ detector = exposure.getDetector() - transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS), - detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE)) + transform = detector.getTransformMap().getTransform( + detector.makeCameraSys(afwCameraGeom.PIXELS), detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE) + ) image = exposure.getMaskedImage() maskVal = image.getMask().getPlaneBitMask(self.config.mask) @@ -657,7 +702,7 @@ def addCcd(self, exposure): num = result.getValue(afwMath.NPOINT) if not numpy.isfinite(mean) or not numpy.isfinite(num): continue - warped[xx, yy, afwImage.LOCAL] = mean*num + warped[xx, yy, afwImage.LOCAL] = mean * num warpedCounts[xx, yy, afwImage.LOCAL] = num self._values += warped @@ -681,10 +726,12 @@ def toCcdBackground(self, detector, bbox): bg : `lsst.afw.math.BackgroundList` Background model for CCD. """ - transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS), - detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE)) - binTransform = (geom.AffineTransform.makeScaling(self.config.binning) - * geom.AffineTransform.makeTranslation(geom.Extent2D(0.5, 0.5))) + transform = detector.getTransformMap().getTransform( + detector.makeCameraSys(afwCameraGeom.PIXELS), detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE) + ) + binTransform = geom.AffineTransform.makeScaling( + self.config.binning + ) * geom.AffineTransform.makeTranslation(geom.Extent2D(0.5, 0.5)) # Binned image on CCD --> unbinned image on CCD --> focal plane --> binned focal plane toSample = afwGeom.makeTransform(binTransform).then(transform).then(self.transform) @@ -693,7 +740,7 @@ def toCcdBackground(self, detector, bbox): fpNorm = afwImage.ImageF(focalPlane.getBBox()) fpNorm.set(1.0) - image = afwImage.ImageF(bbox.getDimensions()//self.config.binning) + image = afwImage.ImageF(bbox.getDimensions() // self.config.binning) norm = afwImage.ImageF(image.getBBox()) ctrl = afwMath.WarpingControl("bilinear") afwMath.warpImage(image, focalPlane, toSample.inverted(), ctrl) @@ -706,11 +753,15 @@ def toCcdBackground(self, detector, bbox): image.getArray()[isBad] = image.getArray()[~isBad].mean() return afwMath.BackgroundList( - (afwMath.BackgroundMI(bbox, afwImage.makeMaskedImage(image, mask)), - afwMath.stringToInterpStyle(self.config.interpolation), - afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), - afwMath.ApproximateControl.UNKNOWN, - 0, 0, False) + ( + afwMath.BackgroundMI(bbox, afwImage.makeMaskedImage(image, mask)), + afwMath.stringToInterpStyle(self.config.interpolation), + afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) ) def merge(self, other): @@ -731,8 +782,10 @@ def merge(self, other): The merged background model. """ if (self.config.xSize, self.config.ySize) != (other.config.xSize, other.config.ySize): - raise RuntimeError("Size mismatch: %s vs %s" % ((self.config.xSize, self.config.ySize), - (other.config.xSize, other.config.ySize))) + raise RuntimeError( + "Size mismatch: %s vs %s" + % ((self.config.xSize, self.config.ySize), (other.config.xSize, other.config.ySize)) + ) if self.dims != other.dims: raise RuntimeError("Dimensions mismatch: %s vs %s" % (self.dims, other.dims)) self._values += other._values @@ -761,8 +814,11 @@ def getStatsImage(self): """ values = self._values.clone() values /= self._numbers - thresh = (self.config.minFrac - * (self.config.xSize/self.config.pixelSize)*(self.config.ySize/self.config.pixelSize)) + thresh = ( + self.config.minFrac + * (self.config.xSize / self.config.pixelSize) + * (self.config.ySize / self.config.pixelSize) + ) isBad = self._numbers.getArray() < thresh if self.config.doSmooth: array = values.getArray() @@ -773,11 +829,293 @@ def getStatsImage(self): return values +class TractBackgroundConfig(Config): + """Configuration for TractBackground + + Note that `xBin` and `yBin` are in pixels, as unlike FocalPlaneBackground, + translation from warps to tract and back only requires geometric + transformations in the warped pixel plane. + """ + + xBin = Field(dtype=float, default=500, doc="Bin size in x") + yBin = Field(dtype=float, default=500, doc="Bin size in y") + minFrac = Field(dtype=float, default=0.1, doc="Minimum fraction of bin size for good measurement") + mask = ListField( + dtype=str, + doc="Mask planes to treat as bad", + default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"], + ) + interpolation = ChoiceField( + doc="how to interpolate the background values. This maps to an enum; see afw::math::Background", + dtype=str, + default="AKIMA_SPLINE", + optional=True, + allowed={ + "CONSTANT": "Use a single constant value", + "LINEAR": "Use linear interpolation", + "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints", + "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers", + "NONE": "No background estimation is to be attempted", + }, + ) + doSmooth = Field(dtype=bool, default=False, doc="Do smoothing?") + smoothScale = Field(dtype=float, default=2.0, doc="Smoothing scale, as a multiple of the bin size") + binning = Field(dtype=int, default=64, doc="Binning to use for warp background model (pixels)") + + +class TractBackground: + """ + As FocalPlaneBackground, but works in warped tract coordinates + """ + + @classmethod + def fromSimilar(cls, other): + """Construct from an object that has the same interface. + + Parameters + ---------- + other : `TractBackground`-like + An object that matches the interface of `TractBackground` + but which may be different. + + Returns + ------- + background : `TractBackground` + Something guaranteed to be a `TractBackground`. + """ + return cls(other.config, other.dims, other.transform, other._values, other._numbers) + + def __init__(self, config, values=None, numbers=None): + """Constructor + + Developers should note that changes to the signature of this method + require coordinated changes to the `__reduce__` and `clone` methods. + + Parameters + ---------- + config : `TractBackgroundConfig` + Configuration for measuring tract backgrounds. + transform : `lsst.afw.geom.TransformPoint2ToPoint2` + Transformation from tract coordinates to warp coordinates. + values : `lsst.afw.image.ImageF` + Measured background values. + numbers : `lsst.afw.image.ImageF` + Number of pixels in each background measurement. + """ + self.config = config + # TODO: dynamic tract dimensions? + self.dims = geom.Extent2I(36000 / self.config.xBin, 36000 / self.config.yBin) + + if values is None: + values = afwImage.ImageF(self.dims) + values.set(0.0) + else: + values = values.clone() + assert values.getDimensions() == self.dims + self._values = values + if numbers is None: + numbers = afwImage.ImageF(self.dims) # float for dynamic range and convenience + numbers.set(0.0) + else: + numbers = numbers.clone() + assert numbers.getDimensions() == self.dims + self._numbers = numbers + + def __reduce__(self): + return self.__class__, (self.config, self._values, self._numbers) + + def clone(self): + return self.__class__(self.config, self._values, self._numbers) + + def addWarp(self, warp): + """ + Equivalent to FocalPlaneBackground.addCcd(), but on warps instead. + Bins masked images of warps and adds these values into a blank image + with the binned tract dimensions at the location of the warp in the + tract. + + Parameters + ---------- + warp : `lsst.afw.image.ExposureF` + Warped image corresponding to a single patch in a single visit + """ + image = warp.getMaskedImage() + maskVal = image.getMask().getPlaneBitMask(self.config.mask) + # Photometric scaling necessary for contiguous background across tract + image.image.array *= warp.getPhotoCalib().instFluxToNanojansky(1) + + warped = afwImage.ImageF(self._values.getBBox()) + warpedCounts = afwImage.ImageF(self._numbers.getBBox()) + + stats = afwMath.StatisticsControl() + stats.setAndMask(maskVal) + stats.setNanSafe(True) + + # Pixel locations in binned tract-scale image + pixels = itertools.product( + numpy.arange(warped.getBBox().getMinX(), warped.getBBox().getMaxX() + 1), + numpy.arange(warped.getBBox().getMinY(), warped.getBBox().getMaxY() + 1), + ) + for xx, yy in pixels: + llc = geom.Point2D((xx - 0.5) * self.config.xBin, (yy - 0.5) * self.config.yBin) + urc = geom.Point2D( + (xx + 0.5) * self.config.xBin + self.config.xBin - 1, + (yy + 0.5) * self.config.yBin + self.config.yBin - 1, + ) + bbox = geom.Box2I(geom.Point2I(llc), geom.Point2I(urc)) + bbox.clip(image.getBBox()) # Works in tract coordinates + if bbox.isEmpty(): + continue + subImage = image.Factory(image, bbox) + result = afwMath.makeStatistics(subImage, afwMath.MEANCLIP | afwMath.NPOINT, stats) + mean = result.getValue(afwMath.MEANCLIP) + num = result.getValue(afwMath.NPOINT) + if not numpy.isfinite(mean) or not numpy.isfinite(num): + continue + warped[xx, yy, afwImage.LOCAL] = mean * num + warpedCounts[xx, yy, afwImage.LOCAL] = num + + self._values += warped + self._numbers += warpedCounts + + def merge(self, other): + """Merge with another TractBackground + + This allows multiple background models to be constructed from + different warps, and then merged to form a single consistent + background model for the entire tract. + + Parameters + ---------- + other : `TractBackground` + Another background model to merge. + + Returns + ------- + self : `TractBackground` + The merged background model. + """ + if (self.config.xBin, self.config.yBin) != (other.config.xBin, other.config.yBin): + raise RuntimeError( + "Size mismatch: %s vs %s" + % ((self.config.xBin, self.config.yBin), (other.config.xBin, other.config.yBin)) + ) + if self.dims != other.dims: + raise RuntimeError("Dimensions mismatch: %s vs %s" % (self.dims, other.dims)) + self._values += other._values + self._numbers += other._numbers + return self + + def __iadd__(self, other): + """Merge with another TractBackground + + Parameters + ---------- + other : `TractBackground` + Another background model to merge. + + Returns + ------- + self : `TractBackground` + The merged background model. + """ + return self.merge(other) + + def toWarpBackground(self, warp): + """ + Equivalent of FocalPlaneBackground.toCcdBackground(), but creates a + background model for a warp using a full tract model. + + Parameters + ---------- + warp : `lsst.afw.image.ExposureF` + Warped image corresponding to a single patch in a single visit + + Returns + ------- + bg : `lsst.afw.math.BackgroundList` + Background model for warp + """ + # Transform to binned warp plane + binTransform = geom.AffineTransform.makeScaling(self.config.binning) + + # Transform from binned tract plane to tract plane + # Start at the patch corner, not the warp corner overlap region + corner = warp.getBBox().getMin() + if corner[0] % 4000 != 0: # TODO: hard-coded patch dimensions are bad + corner[0] += 100 + corner[1] += 100 + offset = geom.Extent2D(corner[0], corner[1]) + tractTransform = ( + geom.AffineTransform.makeTranslation(geom.Extent2D(-0.5, -0.5)) + * geom.AffineTransform.makeScaling(1.0 / self.config.xBin, 1.0 / self.config.yBin) + * geom.AffineTransform.makeTranslation(offset) + ) + transform = afwGeom.makeTransform(tractTransform) + + # Full transform + toSample = afwGeom.makeTransform(binTransform).then(transform) + + # Full tract sky model and normalization array + tractPlane = self.getStatsImage() + tpNorm = afwImage.ImageF(tractPlane.getBBox()) + tpNorm.set(1.0) + + # Binned warp image and normalization array + bbox = warp.getBBox() + image = afwImage.ImageF(bbox.getDimensions() // self.config.binning) + norm = afwImage.ImageF(image.getBBox()) + + ctrl = afwMath.WarpingControl("bilinear") + afwMath.warpImage(image, tractPlane, toSample.inverted(), ctrl) + afwMath.warpImage(norm, tpNorm, toSample.inverted(), ctrl) + image /= norm + # Convert back to counts so the model can be subtracted w/o conversion + image /= warp.getPhotoCalib().instFluxToNanojansky(1) + + mask = afwImage.Mask(image.getBBox()) + isBad = numpy.isnan(image.getArray()) + mask.getArray()[isBad] = mask.getPlaneBitMask("BAD") + image.getArray()[isBad] = image.getArray()[~isBad].mean() + + return afwMath.BackgroundList( + ( + afwMath.BackgroundMI(warp.getBBox(), afwImage.makeMaskedImage(image, mask)), + afwMath.stringToInterpStyle(self.config.interpolation), + afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), + afwMath.ApproximateControl.UNKNOWN, + 0, + 0, + False, + ) + ) + + def getStatsImage(self): + """Return the background model data + + This is the measurement of the background for each of the superpixels. + """ + values = self._values.clone() + values /= self._numbers + # This old logic doesn't work because everything outside the FP is bad + # thresh = self.config.minFrac * (self.config.xBin) * (self.config.yBin) + # isBad = self._numbers.getArray() < thresh + # if self.config.doSmooth: + # array = values.getArray() + # array[:] = smoothArray(array, isBad, self.config.smoothScale) + # isBad = numpy.isnan(values.array) + # if numpy.any(isBad): + # interpolateBadPixels(values.getArray(), isBad, self.config.interpolation) + return values + + class MaskObjectsConfig(Config): """Configuration for MaskObjectsTask""" + nIter = Field(dtype=int, default=3, doc="Number of iterations") - subtractBackground = ConfigurableField(target=measAlg.SubtractBackgroundTask, - doc="Background subtraction") + subtractBackground = ConfigurableField( + target=measAlg.SubtractBackgroundTask, doc="Background subtraction" + ) detection = ConfigurableField(target=measAlg.SourceDetectionTask, doc="Source detection") detectSigma = Field(dtype=float, default=5.0, doc="Detection threshold (standard deviations)") doInterpolate = Field(dtype=bool, default=True, doc="Interpolate when removing objects?") @@ -794,11 +1132,15 @@ def setDefaults(self): self.interpolate.useApprox = False def validate(self): - if (self.detection.reEstimateBackground - or self.detection.doTempLocalBackground - or self.detection.doTempWideBackground): - raise RuntimeError("Incorrect settings for object masking: reEstimateBackground, " - "doTempLocalBackground and doTempWideBackground must be False") + if ( + self.detection.reEstimateBackground + or self.detection.doTempLocalBackground + or self.detection.doTempWideBackground + ): + raise RuntimeError( + "Incorrect settings for object masking: reEstimateBackground, " + "doTempLocalBackground and doTempWideBackground must be False" + ) class MaskObjectsTask(Task): @@ -812,6 +1154,7 @@ class MaskObjectsTask(Task): We deliberately use the specified ``detectSigma`` instead of the PSF, in order to better pick up the faint wings of objects. """ + ConfigClass = MaskObjectsConfig def __init__(self, *args, **kwargs): @@ -903,7 +1246,7 @@ def smoothArray(array, bad, sigma): """ convolved = gaussian_filter(numpy.where(bad, 0.0, array), sigma, mode="constant", cval=0.0) denominator = gaussian_filter(numpy.where(bad, 0.0, 1.0), sigma, mode="constant", cval=0.0) - return convolved/denominator + return convolved / denominator def _create_module_child(name): diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 8906dba5c..4f45b67eb 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -39,15 +39,20 @@ stringToStatisticsProperty, stringToUndersampleStyle, ) -from lsst.pex.config import ChoiceField, Field, ListField, RangeField +from lsst.pex.config import ChoiceField, ConfigField, Field, ListField, RangeField from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct, TaskError from lsst.pipe.base.connectionTypes import Input, Output +from lsst.pipe.tasks.background import TractBackground, TractBackgroundConfig from lsst.utils.timer import timeMethod class MatchBackgroundsConnections( PipelineTaskConnections, - dimensions=("skymap", "tract", "patch", "band"), + # How to get it to do visit, warp? + # Need to kill all collections w/new dataset types to try other combos here... + # https://pipelines.lsst.io/v/weekly/modules/lsst.pipe.base/creating-a-pipelinetask.html#pipelinetask-processing-multiple-datasets + # dimensions=("skymap", "tract", "patch", "band"), + dimensions=("skymap", "tract", "band"), # Don't want to mix bands... defaultTemplates={ "inputCoaddName": "deep", "outputCoaddName": "deep", @@ -67,7 +72,8 @@ class MatchBackgroundsConnections( backgroundInfoList = Output( doc="List of differential backgrounds, with goodness of fit params.", name="psfMatchedWarpBackground_diff", # This needs to change - dimensions=("skymap", "tract", "patch", "visit"), + # dimensions=("skymap", "tract", "patch", "visit"), + dimensions=("skymap", "tract", "visit", "patch"), storageClass="Background", multiple=True, ) @@ -75,7 +81,8 @@ class MatchBackgroundsConnections( doc="List of background-matched warps.", name="{inputCoaddName}Coadd_{warpType}Warp_bgMatched", storageClass="ExposureF", - dimensions=("skymap", "tract", "patch", "visit"), + # dimensions=("skymap", "tract", "patch", "visit"), + dimensions=("skymap", "tract", "visit", "patch"), multiple=True, ) @@ -121,6 +128,10 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr ) # Background matching + tractBgModel = ConfigField( + dtype=TractBackgroundConfig, + doc="Background model for the entire tract", + ) usePolynomial = Field[bool]( doc="Fit background difference with a Chebychev polynomial interpolation? " "If False, fit with spline interpolation instead.", @@ -140,7 +151,6 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr "BAD", "INTRP", "CR", - "NOT_DEBLENDED", ], itemCheck=lambda x: x in Mask().getMaskPlaneDict(), ) @@ -277,9 +287,13 @@ def run(self, warps): """ if (numExp := len(warps)) < 1: raise TaskError("No exposures to match") + # TODO: store ref visit ID between runs, to skip selection process? + + # First, build FFP BG models of each visit + visitTractBgs = self._makeTractBackgrounds(warps) # Define a reference warp; 'warps' is modified in-place to exclude it - refWarp, refInd, bkgd = self._defineWarps(warps=warps, refWarpVisit=self.config.refWarpVisit) + refBg, refVis, bkgd = self._defineWarps(visitTractBgs, refVisitId=self.config.refWarpVisit) # Images must be scaled to a common ZP # Converting everything to nJy to accomplish this @@ -335,109 +349,187 @@ def run(self, warps): return Struct(backgroundInfoList=backgroundInfoList, matchedImageList=matchedImageList) @timeMethod - def _defineWarps(self, warps, refWarpVisit=None): - """Define the reference warp and list of comparison warps. - - If no reference warp ID is supplied, this method calculates an - appropriate reference exposure from the supplied list of warps by - minimizing a cost function that penalizes high background complexity - (divergence from a plane), high variance, low global coverage, and low - edge coverage. + def _makeTractBackgrounds(self, warps, refWarpDDF=None): + """Create full tract model of the backgrounds of all visits. Parameters ---------- warps : `list`[`~lsst.daf.butler.DeferredDatasetHandle`] List of warped exposures (of type `~lsst.afw.image.ExposureF`). - refWarpVisit : `int`, optional - Visit ID of the reference warp. - If None, the best warp is chosen from the list of warps. + This is ordered by patch ID, then by visit ID + + refWarpDDF : ``~lsst.daf.butler.DeferredDatasetHandle` optional + Chosen reference warp to match to, if supplied + + Returns + ------- + visitTractBackrounds : `dict`{`TractBackground`} + Models of full tract backgrounds for all visits, accessed by visit + IDs + """ + # First, separate warps by visit + visits = np.unique([i.dataId["visit"] for i in warps]) + + # Then build background models for each visit, and store + visitTractBackgrounds = {} + for i in range(len(visits)): + # Find all the warps for the visit + visitWarps = [j for j in warps if j.dataId["visit"] == visits[i]] + # Set up empty full tract background model object + bgModelBase = TractBackground(self.config.tractBgModel) + + bgModels = [] + for warp in visitWarps: + msg = "Constructing FFP background model for visit %d using %d patches" + self.log.debug( + msg, + visits[i], + len(warps), + ) + if refWarpDDF is not None: + msg = "Doing difference imaging: reference warp visit ID: %d" + self.log.debug(msg, refWarpDDF.dataId["visit"]) + visitWarp = warp.get() + # If a reference image is supplied, make a BG of the difference + if refWarpDDF is not None: + refWarp = refWarpDDF.get() + visitWarp = refWarp - visitWarp + bgModel = bgModelBase.clone() + bgModel.addWarp(visitWarp) + bgModels.append(bgModel) + + # Merge warp models to make a single full tract background model + for bgModel, warp in zip(bgModels, visitWarps): + msg = ( + "Patch %d: Merging %d unmasked pixels (%.1f%s of detector area) into tract plane BG " + "model" + ) + self.log.debug( + msg, + warp.dataId["patch"], + bgModel._numbers.getArray().sum(), + 100 * bgModel._numbers.getArray().sum() / visitWarp.getBBox().getArea(), + "%", + ) + bgModelBase.merge(bgModel) + visitTractBackgrounds[visits[i]] = bgModelBase + return visitTractBackgrounds + + @timeMethod + def _defineWarps(self, visitTractBackgrounds, refVisitId=None): + """Define the reference visit and list of comparison visits. + + If no reference visit ID is supplied, this method calculates an + appropriate reference exposure from the supplied list of visit + backgrounds by minimizing a cost function that penalizes high + background complexity (divergence from a plane), high variance, and low + global coverage. + + Parameters + ---------- + visitTractBackgrounds : `dict`{`TractBackground`} + Models of full tract backgrounds for all visits, accessed by visit + IDs + refVisitId : `int`, optional + ID of the reference visit. + If None, the best visit is chosen using the dictionary of existing + backgrounds. Returns ------- - refWarp : `~lsst.afw.image.ExposureF` - Reference warped exposure. - refWarpIndex : `int` - Index of the reference removed from the list of warps. - warpBg : `~lsst.afw.math.BackgroundMI` + refBg : `~lsst.afw.math.BackgroundMI` + Reference background to match to. + refVis : `int` + Index of the reference visit removed from the dictionary. + fitBg : `~lsst.afw.math.BackgroundMI` Temporary background model, used to make a blank BG for the ref Notes ----- This method modifies the input list of warps in place by removing the reference warp from it. + """ # User-defined reference visit, if one has been supplied - if refWarpVisit: - warpVisits = [warpDDH.dataId["visit"] for warpDDH in warps] + if refVisitId: + visits = [visId for visId in visitTractBackgrounds.keys()] try: - refWarpIndex = warpVisits.index(refWarpVisit) - refWarpDDH = warps.pop(refWarpIndex) - self.log.info("Using user-supplied reference visit %d", refWarpVisit) - return refWarpDDH.get(), refWarpIndex + refTractBackground = visitTractBackgrounds.pop(refVisitId) + self.log.info("Using user-supplied reference visit %d", refVisitId) + # TODO: need to return a background object here! + return refTractBackground, refVisitId except ValueError: - raise TaskError(f"Reference visit {refWarpVisit} is not found in the list of warps.") + raise TaskError(f"Reference visit {refVisitId} is not found in the list of warps.") # Extract mean/var/npoints for each warp - warpChi2s = [] # Background goodness of fit - warpVars = [] # Variance - warpNPointsGlobal = [] # Global coverage - warpNPointsEdge = [] # Edge coverage - for warpDDH in warps: - warp = warpDDH.get() - instFluxToNanojansky = warp.getPhotoCalib().instFluxToNanojansky(1) - warp.image *= instFluxToNanojansky # Images in nJy to facilitate difference imaging - warp.variance *= instFluxToNanojansky - warpBg, _ = self._makeBackground(warp.getMaskedImage(), binSize=self.config.chi2BinSize) + fitChi2s = [] # Background goodness of fit + # warpVars = [] # Variance + fitNPointsGlobal = [] # Global coverage + # warpNPointsEdge = [] # Edge coverage + visits = [] # To ensure dictionary key order is correct + for vis in visitTractBackgrounds: + visits.append(vis) + # Fit a model to the FFP + # TODO: need a variance plane in the tractBg as well + tractBg = visitTractBackgrounds[vis].getStatsImage() + fitBg, _ = self._makeBackground(tractBg, binSize=1) + fitBg.getStatsImage().variance = ImageF(np.sqrt(fitBg.getStatsImage().image.array)) # Return an approximation to the background approxCtrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, self.config.approxWeighting) - warpApprox = warpBg.getApproximate(approxCtrl, self.undersampleStyle) - warpBgSub = ImageF(warp.image.array - warpApprox.getImage().array) + fitApprox = fitBg.getApproximate(approxCtrl, self.undersampleStyle) - warpStats = makeStatistics(warpBgSub, warp.mask, VARIANCE | NPOINT, self.statsCtrl) + fitBgSub = MaskedImageF(ImageF(tractBg.array - fitApprox.getImage().array)) + bad_mask_bit_mask = fitBgSub.mask.getPlaneBitMask("BAD") + fitBgSub.mask.array[np.isnan(fitBgSub.image.array)] = bad_mask_bit_mask - bad_mask_bit_mask = warp.mask.getPlaneBitMask(self.config.badMaskPlanes) - good = (warp.mask.array.astype(int) & bad_mask_bit_mask) == 0 + fitStats = makeStatistics(fitBgSub.image, fitBgSub.mask, VARIANCE | NPOINT, self.statsCtrl) + + good = (fitBgSub.mask.array.astype(int) & bad_mask_bit_mask) == 0 dof = len(good[good]) - 6 # Assuming eq. of plane - warpChi2 = np.nansum(warpBgSub.array[good] ** 2 / warp.variance.array[good]) / dof - warpVar, _ = warpStats.getResult(VARIANCE) - warpNPointGlobal, _ = warpStats.getResult(NPOINT) - warpNPointEdge = ( - np.sum(~np.isnan(warp.image.array[:, 0])) # Left edge - + np.sum(~np.isnan(warp.image.array[:, -1])) # Right edge - + np.sum(~np.isnan(warp.image.array[0, :])) # Bottom edge - + np.sum(~np.isnan(warp.image.array[-1, :])) # Top edge + fitChi2 = ( + np.nansum(fitBgSub.image.array[good] ** 2 / fitBg.getStatsImage().variance.array[good]) / dof ) - warpChi2s.append(warpChi2) - warpVars.append(warpVar) - warpNPointsGlobal.append(int(warpNPointGlobal)) - warpNPointsEdge.append(warpNPointEdge) + # fitVar, _ = fitStats.getResult(VARIANCE) # Will add this back later + fitNPointGlobal, _ = fitStats.getResult(NPOINT) + # This becomes pointless: it's a circle within a square. Coverage is what's needed for that. + # warpNPointEdge = ( + # np.sum(~np.isnan(warp.image.array[:, 0])) # Left edge + # + np.sum(~np.isnan(warp.image.array[:, -1])) # Right edge + # + np.sum(~np.isnan(warp.image.array[0, :])) # Bottom edge + # + np.sum(~np.isnan(warp.image.array[-1, :])) # Top edge + # ) + fitChi2s.append(fitChi2) + # warpVars.append(warpVar) + fitNPointsGlobal.append(int(fitNPointGlobal)) + # warpNPointsEdge.append(warpNPointEdge) self.log.info( - "Sci exp. visit %d; BG fit Chi^2=%.1f, var=%.1f nJy, nPoints global=%d, nPoints edge=%d", - warp.getInfo().getVisitInfo().id, - warpChi2, - warpVar, - warpNPointGlobal, - warpNPointEdge, + # "Sci exp. visit %d; BG fit Chi^2=%.1f, var=%.1f nJy, nPoints global=%d, nPoints edge=%d", + "Sci exp. visit %d; BG fit Chi^2=%.1f, nPoints global=%d", + vis, + fitChi2, + # warpVar, + fitNPointGlobal, + # warpNPointEdge, ) - # Normalize mean/var/npoints to range from 0 to 1 - warpChi2sFrac = np.array(warpChi2s) / np.nanmax(warpChi2s) - warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) - warpNPointsGlobalFrac = np.nanmin(warpNPointsGlobal) / np.array(warpNPointsGlobal) - warpNPointsEdgeFrac = np.nanmin(warpNPointsEdge) / np.array(warpNPointsEdge) + fitChi2sFrac = np.array(fitChi2s) / np.nanmax(fitChi2s) + # warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) + fitNPointsGlobalFrac = np.nanmin(fitNPointsGlobal) / np.array(fitNPointsGlobal) + # warpNPointsEdgeFrac = np.nanmin(warpNPointsEdge) / np.array(warpNPointsEdge) # Calculate cost function values - costFunctionVals = self.config.bestRefWeightChi2 * warpChi2sFrac - costFunctionVals += self.config.bestRefWeightVariance * warpVarsFrac - costFunctionVals += self.config.bestRefWeightGlobalCoverage * warpNPointsGlobalFrac - costFunctionVals += self.config.bestRefWeightEdgeCoverage * warpNPointsEdgeFrac + costFunctionVals = self.config.bestRefWeightChi2 * fitChi2sFrac + # costFunctionVals += self.config.bestRefWeightVariance * warpVarsFrac + costFunctionVals += self.config.bestRefWeightGlobalCoverage * fitNPointsGlobalFrac + # costFunctionVals += self.config.bestRefWeightEdgeCoverage * warpNPointsEdgeFrac ind = np.nanargmin(costFunctionVals) - refWarp = warps.pop(ind) - self.log.info("Using best reference visit %d", refWarp.dataId["visit"]) - return refWarp, ind, warpBg + refVis = visits[ind] + refBg = visitTractBackgrounds.pop(refVis) + self.log.info("Using best reference visit %d", refVis) + return refBg, refVis, fitBg def _makeBackground(self, warp: MaskedImageF, binSize) -> tuple[BackgroundMI, BackgroundControl]: """Generate a simple binned background masked image for warped data. From c652ce09a7eb9df170df192d0b7c2c5612fc63b8 Mon Sep 17 00:00:00 2001 From: Aaron Watkins Date: Thu, 17 Apr 2025 07:36:33 -0700 Subject: [PATCH 13/13] Full tract background refactor AEW All methods, including run, updated to function properly with tract backgrounds rather than warps. Task completes without error when run on three full visits, and output appears roughly correct. --- python/lsst/pipe/tasks/background.py | 9 +- python/lsst/pipe/tasks/matchBackgrounds.py | 424 ++++++++++++--------- 2 files changed, 253 insertions(+), 180 deletions(-) diff --git a/python/lsst/pipe/tasks/background.py b/python/lsst/pipe/tasks/background.py index ae90d4167..643a61f68 100644 --- a/python/lsst/pipe/tasks/background.py +++ b/python/lsst/pipe/tasks/background.py @@ -860,7 +860,7 @@ class TractBackgroundConfig(Config): ) doSmooth = Field(dtype=bool, default=False, doc="Do smoothing?") smoothScale = Field(dtype=float, default=2.0, doc="Smoothing scale, as a multiple of the bin size") - binning = Field(dtype=int, default=64, doc="Binning to use for warp background model (pixels)") + binning = Field(dtype=int, default=200, doc="Binning to use for warp background model (pixels)") class TractBackground: @@ -895,8 +895,6 @@ def __init__(self, config, values=None, numbers=None): ---------- config : `TractBackgroundConfig` Configuration for measuring tract backgrounds. - transform : `lsst.afw.geom.TransformPoint2ToPoint2` - Transformation from tract coordinates to warp coordinates. values : `lsst.afw.image.ImageF` Measured background values. numbers : `lsst.afw.image.ImageF` @@ -1097,7 +1095,10 @@ def getStatsImage(self): """ values = self._values.clone() values /= self._numbers - # This old logic doesn't work because everything outside the FP is bad + # TODO: this logic smoothes over bad parts of the image, including NaN + # values. But it doesn't work here because NaN pixels are always found + # at the image edges. Could ignore it, or devise an alternative. + # tract have no overlap with the visit? # thresh = self.config.minFrac * (self.config.xBin) * (self.config.yBin) # isBad = self._numbers.getArray() < thresh # if self.config.doSmooth: diff --git a/python/lsst/pipe/tasks/matchBackgrounds.py b/python/lsst/pipe/tasks/matchBackgrounds.py index 4f45b67eb..17cf53c65 100644 --- a/python/lsst/pipe/tasks/matchBackgrounds.py +++ b/python/lsst/pipe/tasks/matchBackgrounds.py @@ -23,22 +23,10 @@ import numpy as np from lsst.afw.image import LOCAL, ImageF, Mask, MaskedImageF -from lsst.afw.math import ( - MEANCLIP, - MEANSQUARE, - NPOINT, - VARIANCE, - ApproximateControl, - BackgroundControl, - BackgroundList, - BackgroundMI, - StatisticsControl, - makeBackground, - makeStatistics, - stringToInterpStyle, - stringToStatisticsProperty, - stringToUndersampleStyle, -) +from lsst.afw.math import (MEANCLIP, MEANSQUARE, NPOINT, STDEV, VARIANCE, ApproximateControl, + BackgroundControl, BackgroundList, BackgroundMI, StatisticsControl, makeBackground, + makeStatistics, stringToInterpStyle, stringToStatisticsProperty, + stringToUndersampleStyle) from lsst.pex.config import ChoiceField, ConfigField, Field, ListField, RangeField from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct, TaskError from lsst.pipe.base.connectionTypes import Input, Output @@ -48,11 +36,7 @@ class MatchBackgroundsConnections( PipelineTaskConnections, - # How to get it to do visit, warp? - # Need to kill all collections w/new dataset types to try other combos here... - # https://pipelines.lsst.io/v/weekly/modules/lsst.pipe.base/creating-a-pipelinetask.html#pipelinetask-processing-multiple-datasets - # dimensions=("skymap", "tract", "patch", "band"), - dimensions=("skymap", "tract", "band"), # Don't want to mix bands... + dimensions=("skymap", "tract", "band"), defaultTemplates={ "inputCoaddName": "deep", "outputCoaddName": "deep", @@ -71,8 +55,7 @@ class MatchBackgroundsConnections( ) backgroundInfoList = Output( doc="List of differential backgrounds, with goodness of fit params.", - name="psfMatchedWarpBackground_diff", # This needs to change - # dimensions=("skymap", "tract", "patch", "visit"), + name="psfMatchedWarpBackground_diff", # TODO: come up with better name dimensions=("skymap", "tract", "visit", "patch"), storageClass="Background", multiple=True, @@ -85,6 +68,7 @@ class MatchBackgroundsConnections( dimensions=("skymap", "tract", "visit", "patch"), multiple=True, ) + # TODO: should we also output the full tract background images? class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgroundsConnections): @@ -110,6 +94,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr min=0.0, max=1.0, ) + # TODO: removed edge coverage, so rename this? bestRefWeightGlobalCoverage = RangeField( dtype=float, doc="Global coverage weight (total number of valid pixels) when calculating the best reference " @@ -118,14 +103,6 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr min=0.0, max=1.0, ) - bestRefWeightEdgeCoverage = RangeField( - dtype=float, - doc="Edge coverage weight (number of valid edge pixels) when calculating the best reference " - "exposure. Higher weights prefer exposures with high coverage. Ignored when a ref visit supplied.", - default=0.4, - min=0.0, - max=1.0, - ) # Background matching tractBgModel = ConfigField( @@ -139,8 +116,10 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr ) order = Field[int]( doc="Order of Chebyshev polynomial background model. Ignored if ``usePolynomial=False``.", - default=8, + default=1, ) + # TODO: because the binning has moved to TractBackground, the mask planes + # are also now defined there. I think only NaNSafe is necessary now. badMaskPlanes = ListField[str]( doc="Names of mask planes to ignore while estimating the background.", default=[ @@ -171,6 +150,8 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.", }, ) + # TODO: bin sizes have been moved to background.py TractBackground configs + # so both of these have no purpose anymore. binSize = Field[int]( doc="Bin size for gridding the difference image and fitting a spatial model.", default=256, @@ -201,6 +182,7 @@ class MatchBackgroundsConfig(PipelineTaskConfig, pipelineConnections=MatchBackgr default=3, ) + # TODO: this might have been made redundant due to how binning proceeds now approxWeighting = Field[bool]( doc="Use inverse-variance weighting when approximating the background offset model? This will fail " "when the background offset is constant (usually only the case in testing with artificial images)." @@ -225,15 +207,12 @@ class MatchBackgroundsTask(PipelineTask): science exposure. The reference exposure is chosen from the list of science exposures by minimizing a cost function that penalizes high background complexity - (divergence from a plane), high variance, low global coverage, and low - coverage along image edges. - The cost function is a weighted sum of these four metrics. + (divergence from a plane), high variance, and low global coverage. + The cost function is a weighted sum of these three metrics. The weights are set by the config parameters: - ``bestRefWeightChi2`` - ``bestRefWeightVariance`` - ``bestRefWeightGlobalCoverage`` - - ``bestRefWeightEdgeCoverage`` - Attributes ---------- @@ -265,9 +244,10 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): @timeMethod def run(self, warps): - """Match the backgrounds of a list of warped exposures to a reference. + """Match the backgrounds of a list of warped exposures to the same + patches in a reference visit. - A reference warp will be chosen automatically if none is supplied. + A reference visit ID will be chosen automatically if none is supplied. Parameters ---------- @@ -276,9 +256,9 @@ def run(self, warps): Returns ------- - result : `~lsst.afw.math.BackgroundList` - Differential background model - Add this to the science exposure to match the reference exposure. + result : `~lsst.afw.math.BackgroundList`, `~lsst.afw.image.Exposure` + Differential background model and associated background-matched + image. Raises ------ @@ -287,70 +267,32 @@ def run(self, warps): """ if (numExp := len(warps)) < 1: raise TaskError("No exposures to match") - # TODO: store ref visit ID between runs, to skip selection process? - - # First, build FFP BG models of each visit - visitTractBgs = self._makeTractBackgrounds(warps) - # Define a reference warp; 'warps' is modified in-place to exclude it - refBg, refVis, bkgd = self._defineWarps(visitTractBgs, refVisitId=self.config.refWarpVisit) + if self.config.refWarpVisit is None: + # First, build FFP BG models of each visit + # Dictionary of background models (in nJy) accessed by visit ID + visitTractBgs = self._makeTractBackgrounds(warps) - # Images must be scaled to a common ZP - # Converting everything to nJy to accomplish this - refExposure = refWarp.get() - instFluxToNanojanskyRef = self._fluxScale(refExposure) + refVisId = self._defineWarps(visitTractBgs) + else: + self.log.info("Using user-supplied reference visit %d", + self.config.refWarpVisit) + refVisId = self.config.refWarpVisit self.log.info("Matching %d Exposures", numExp) - # Blank ref warp background as reference background - bkgdIm = bkgd.getImageF() - bkgdStatsIm = bkgd.getStatsImage() - bkgdIm *= 0 - bkgdStatsIm *= 0 - blank = BackgroundList( - ( - bkgd, - stringToInterpStyle(self.config.interpStyle), - stringToUndersampleStyle(self.config.undersampleStyle), - ApproximateControl.UNKNOWN, - 0, - 0, - False, - ) - ) + backgroundInfoList, matchedImageList = self.matchBackgrounds(warps, refVisId) + + # TODO: should we include an optional check on the differential BG + # models reconstructed at the warp resolution from the full tract BG? - backgroundInfoList = [] - matchedImageList = [] - for exp in warps: - self.log.info( - "Matching background of %s to %s", - exp.dataId, - refWarp.dataId, - ) - toMatchExposure = exp.get() - instFluxToNanojansky = self._fluxScale(toMatchExposure) - try: - backgroundInfoStruct = self.matchBackgrounds( - refExposure=refExposure, - sciExposure=toMatchExposure, - ) - backgroundInfoStruct.isReference = False - except Exception as e: - self.log.warning("Failed to fit background %s: %s", exp.dataId, e) - backgroundInfoStruct = blank - - backgroundInfoList.append(backgroundInfoStruct) - toMatchExposure.image /= instFluxToNanojansky # Back to cts - matchedImageList.append(toMatchExposure) - - backgroundInfoList.insert(refInd, blank) - refExposure.image /= instFluxToNanojanskyRef # Back to cts - matchedImageList.insert(refInd, refExposure) return Struct(backgroundInfoList=backgroundInfoList, matchedImageList=matchedImageList) @timeMethod - def _makeTractBackgrounds(self, warps, refWarpDDF=None): - """Create full tract model of the backgrounds of all visits. + def _makeTractBackgrounds(self, warps, refVisitId=None): + """Create full tract models of the backgrounds of all visits. If a + reference visit ID is supplied, create full tract models of difference + image backgrounds of all visits. Parameters ---------- @@ -358,48 +300,64 @@ def _makeTractBackgrounds(self, warps, refWarpDDF=None): List of warped exposures (of type `~lsst.afw.image.ExposureF`). This is ordered by patch ID, then by visit ID - refWarpDDF : ``~lsst.daf.butler.DeferredDatasetHandle` optional - Chosen reference warp to match to, if supplied + refVisitId : `int` optional + Chosen reference visit ID to match to, if supplied Returns ------- visitTractBackrounds : `dict`{`TractBackground`} - Models of full tract backgrounds for all visits, accessed by visit - IDs + Models of full tract backgrounds for all visits, in nanojanskies. + Accessed by visit ID. + + Notes + ----- + Input warps, including reference if ID is supplied, are converted in- + place to nanojanskies. """ # First, separate warps by visit visits = np.unique([i.dataId["visit"] for i in warps]) - # Then build background models for each visit, and store + # Then build background models for each visit and store visitTractBackgrounds = {} for i in range(len(visits)): - # Find all the warps for the visit - visitWarps = [j for j in warps if j.dataId["visit"] == visits[i]] + visitWarpDDFs = [j for j in warps if j.dataId["visit"] == visits[i]] # Set up empty full tract background model object bgModelBase = TractBackground(self.config.tractBgModel) bgModels = [] - for warp in visitWarps: + for warp in visitWarpDDFs: msg = "Constructing FFP background model for visit %d using %d patches" self.log.debug( msg, visits[i], - len(warps), + len(visitWarpDDFs), ) - if refWarpDDF is not None: + if refVisitId is not None: msg = "Doing difference imaging: reference warp visit ID: %d" - self.log.debug(msg, refWarpDDF.dataId["visit"]) + self.log.debug(msg, refVisitId) visitWarp = warp.get() - # If a reference image is supplied, make a BG of the difference - if refWarpDDF is not None: - refWarp = refWarpDDF.get() - visitWarp = refWarp - visitWarp + self._fluxScale(visitWarp) + + # If a reference visit is supplied, makes model of difference + # image backgrounds + if refVisitId is not None: + patchId = warp.dataId["patch"] + refWarpDDFs = [j for j in warps if j.dataId["visit"] == refVisitId] + refPatches = [j.dataId["patch"] for j in refWarpDDFs] + try: + idx = refPatches.index(patchId) + refWarp = refWarpDDFs[idx].get() + except ValueError: + refWarp = visitWarp.clone() + refWarp.image += np.nan + self._fluxScale(refWarp) + visitWarp.image.array = refWarp.image.array - visitWarp.image.array bgModel = bgModelBase.clone() bgModel.addWarp(visitWarp) bgModels.append(bgModel) # Merge warp models to make a single full tract background model - for bgModel, warp in zip(bgModels, visitWarps): + for bgModel, warp in zip(bgModels, visitWarpDDFs): msg = ( "Patch %d: Merging %d unmasked pixels (%.1f%s of detector area) into tract plane BG " "model" @@ -412,68 +370,114 @@ def _makeTractBackgrounds(self, warps, refWarpDDF=None): "%", ) bgModelBase.merge(bgModel) + + if refVisitId is not None and visits[i] != refVisitId: + # If using a diff image, fit that to extrapolate into regions + # of no overlap between visit and reference visit + + # Some config and input checks if config.usePolynomial: + # 1) Check that order/bin size make sense + # 2) Change order if underconstrained + bgModelImage = bgModelBase.getStatsImage() + if self.config.usePolynomial: + order = self.config.order + dimX, dimY = bgModelImage.array.shape + stats = makeStatistics(bgModelImage, NPOINT | STDEV, self.statsCtrl) + npoints, _ = stats.getResult(NPOINT) + stdev, _ = stats.getResult(STDEV) + if stdev < self.config.gridStdevEpsilon: + stdev = self.config.gridStdevEpsilon + minNumberGridPoints = min(dimX, dimY) + if npoints == 0: + raise ValueError("No overlap with reference. Nothing to match") + elif minNumberGridPoints <= order: + # Must lower order or throw exception + if self.config.undersampleStyle == "THROW_EXCEPTION": + raise ValueError("Image does not cover enough of ref image for order and binsize") + elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER": + self.log.warning("Reducing order to %d", (minNumberGridPoints - 1)) + order = minNumberGridPoints - 1 + # TODO: INCREASE_NXNYSAMPLE, i.e. change bin size, is + # harder to pull off now as it involves remaking the + # full-tract background models from scratch. Remove + # that config option? + + # TODO: see below, but using binSize=1 results in all 0s in + # the output variance image, meaning no inverse variance + # weighting is possible. Remove related config? + weightByInverseVariance = False + + bkgd, _ = self._makeBackground(bgModelImage, binSize=1) # Already binned + try: + if self.config.usePolynomial: + actrl = ApproximateControl( + ApproximateControl.CHEBYSHEV, order, order, weightByInverseVariance + ) + undersampleStyle = stringToUndersampleStyle(self.config.undersampleStyle) + approx = bkgd.getApproximate(actrl, undersampleStyle) + bkgdImage = approx.getImage() + else: + bkgdImage = bkgd.getImageF(self.config.interpStyle, self.config.undersampleStyle) + except Exception as e: + raise RuntimeError( + "Background/Approximation failed to interp image %s: %s" % (warp.dataId, e) + ) + # Calculate RMS of fit and log here, too + resids = ImageF(bgModelImage.array - bkgdImage.array) + rms = np.sqrt(np.nanmean(resids.array ** 2)) + mse = makeStatistics(resids, MEANSQUARE, self.statsCtrl).getValue() + + # TODO: because there's no variance image, there's no mean + # variance parameter to measure and report here. + self.log.info( + "Visit %d; difference BG fit RMS=%.2f nJy, matched MSE=%.2f nJy", + visits[i], + rms, + mse, + ) + # Replace statsImage w/model fit + bgModelBase._numbers /= bgModelBase._numbers + bgModelBase._numbers.array[np.isnan(bgModelBase._numbers.array)] = 1. + bgModelBase._values = bkgdImage + visitTractBackgrounds[visits[i]] = bgModelBase return visitTractBackgrounds @timeMethod - def _defineWarps(self, visitTractBackgrounds, refVisitId=None): + def _defineWarps(self, visitTractBackgrounds): """Define the reference visit and list of comparison visits. - If no reference visit ID is supplied, this method calculates an - appropriate reference exposure from the supplied list of visit - backgrounds by minimizing a cost function that penalizes high - background complexity (divergence from a plane), high variance, and low - global coverage. + This method calculates an appropriate reference exposure from the + supplied list of visit backgrounds by minimizing a cost function that + penalizes high background complexity (divergence from a plane), high + variance, and low global coverage. Parameters ---------- visitTractBackgrounds : `dict`{`TractBackground`} Models of full tract backgrounds for all visits, accessed by visit IDs - refVisitId : `int`, optional - ID of the reference visit. - If None, the best visit is chosen using the dictionary of existing - backgrounds. Returns ------- - refBg : `~lsst.afw.math.BackgroundMI` - Reference background to match to. - refVis : `int` + refVisId : `int` Index of the reference visit removed from the dictionary. - fitBg : `~lsst.afw.math.BackgroundMI` - Temporary background model, used to make a blank BG for the ref - - Notes - ----- - This method modifies the input list of warps in place by removing the - reference warp from it. - """ - # User-defined reference visit, if one has been supplied - if refVisitId: - visits = [visId for visId in visitTractBackgrounds.keys()] - try: - refTractBackground = visitTractBackgrounds.pop(refVisitId) - self.log.info("Using user-supplied reference visit %d", refVisitId) - # TODO: need to return a background object here! - return refTractBackground, refVisitId - except ValueError: - raise TaskError(f"Reference visit {refVisitId} is not found in the list of warps.") + # TODO: note that supplying a refID here is pointless, because the only + # thing visitTractBackgrounds is used for is to choose one. # Extract mean/var/npoints for each warp fitChi2s = [] # Background goodness of fit - # warpVars = [] # Variance + fitVars = [] # Variance fitNPointsGlobal = [] # Global coverage - # warpNPointsEdge = [] # Edge coverage visits = [] # To ensure dictionary key order is correct for vis in visitTractBackgrounds: visits.append(vis) # Fit a model to the FFP - # TODO: need a variance plane in the tractBg as well tractBg = visitTractBackgrounds[vis].getStatsImage() fitBg, _ = self._makeBackground(tractBg, binSize=1) - fitBg.getStatsImage().variance = ImageF(np.sqrt(fitBg.getStatsImage().image.array)) + # Use the input image as the variance image + fitBg.getStatsImage().variance = ImageF(tractBg.array) # Return an approximation to the background approxCtrl = ApproximateControl(ApproximateControl.CHEBYSHEV, 1, 1, self.config.approxWeighting) @@ -490,49 +494,37 @@ def _defineWarps(self, visitTractBackgrounds, refVisitId=None): fitChi2 = ( np.nansum(fitBgSub.image.array[good] ** 2 / fitBg.getStatsImage().variance.array[good]) / dof ) - # fitVar, _ = fitStats.getResult(VARIANCE) # Will add this back later + fitVar, _ = fitStats.getResult(VARIANCE) fitNPointGlobal, _ = fitStats.getResult(NPOINT) - # This becomes pointless: it's a circle within a square. Coverage is what's needed for that. - # warpNPointEdge = ( - # np.sum(~np.isnan(warp.image.array[:, 0])) # Left edge - # + np.sum(~np.isnan(warp.image.array[:, -1])) # Right edge - # + np.sum(~np.isnan(warp.image.array[0, :])) # Bottom edge - # + np.sum(~np.isnan(warp.image.array[-1, :])) # Top edge - # ) fitChi2s.append(fitChi2) - # warpVars.append(warpVar) + fitVars.append(fitVar) fitNPointsGlobal.append(int(fitNPointGlobal)) - # warpNPointsEdge.append(warpNPointEdge) self.log.info( - # "Sci exp. visit %d; BG fit Chi^2=%.1f, var=%.1f nJy, nPoints global=%d, nPoints edge=%d", - "Sci exp. visit %d; BG fit Chi^2=%.1f, nPoints global=%d", + "Sci exp. visit %d; BG fit Chi^2=%.2f, var=%.2f nJy, nPoints global=%d", vis, fitChi2, - # warpVar, + fitVar, fitNPointGlobal, - # warpNPointEdge, ) # Normalize mean/var/npoints to range from 0 to 1 fitChi2sFrac = np.array(fitChi2s) / np.nanmax(fitChi2s) - # warpVarsFrac = np.array(warpVars) / np.nanmax(warpVars) + fitVarsFrac = np.array(fitVars) / np.nanmax(fitVars) fitNPointsGlobalFrac = np.nanmin(fitNPointsGlobal) / np.array(fitNPointsGlobal) - # warpNPointsEdgeFrac = np.nanmin(warpNPointsEdge) / np.array(warpNPointsEdge) # Calculate cost function values costFunctionVals = self.config.bestRefWeightChi2 * fitChi2sFrac - # costFunctionVals += self.config.bestRefWeightVariance * warpVarsFrac + costFunctionVals += self.config.bestRefWeightVariance * fitVarsFrac costFunctionVals += self.config.bestRefWeightGlobalCoverage * fitNPointsGlobalFrac - # costFunctionVals += self.config.bestRefWeightEdgeCoverage * warpNPointsEdgeFrac ind = np.nanargmin(costFunctionVals) - refVis = visits[ind] - refBg = visitTractBackgrounds.pop(refVis) - self.log.info("Using best reference visit %d", refVis) - return refBg, refVis, fitBg + refVisitId = visits[ind] + self.log.info("Using best reference visit %d", refVisitId) + return refVisitId def _makeBackground(self, warp: MaskedImageF, binSize) -> tuple[BackgroundMI, BackgroundControl]: - """Generate a simple binned background masked image for warped data. + """Generate a simple binned background masked image for warped or other + data. Parameters ---------- @@ -546,6 +538,9 @@ def _makeBackground(self, warp: MaskedImageF, binSize) -> tuple[BackgroundMI, Ba bgCtrl: `~lsst.afw.math.BackgroundControl` Background control object. """ + # TODO: the only thing this is used for now is to model the + # TractBackground statsImage(), which has no mask and is already + # binned. Simplify this? nx = warp.getWidth() // binSize ny = warp.getHeight() // binSize @@ -573,8 +568,87 @@ def _fluxScale(self, exposure): return fluxZp @timeMethod - def matchBackgrounds(self, refExposure, sciExposure): - """Match science exposure's background level to that of reference + def matchBackgrounds(self, warps, refVisitId): + """Match science exposures' background level to that of reference + exposure. + + Process creates binned images of the full focal plane (in tract + coordinates) for all visit IDs, subtracts each from a similarly + binned FFP reference image, then generates TractBackground + objects. It assumes (but does not require/check) that the mask planes + already have detections set. If detections have not been set/masked, + sources will bias the difference image background estimation. + + The TractBackground objects representing the difference image + backgrounds are then used to generate 'background' models for each warp + comprising the full science exposure visit, which are then added to + each warp. + + Fit diagnostics are also calculated and returned. + + Parameters + ---------- + warps : `list`[`~lsst.daf.butler.DeferredDatasetHandle`] + List of warped exposures (of type `~lsst.afw.image.ExposureF`). + This is ordered by patch ID, then by visit ID + refVisitId : `int` + Chosen reference visit ID to match to + + Returns + ------- + backgroundInfoList : `list`[`TractBackground`] + List of all difference image backgrounds used to match to reference + visit warps, in counts + matchedImageList : `list`[`~lsst.afw.image.ExposureF`] + List of all background-matched warps, in counts + """ + visits = np.unique([i.dataId["visit"] for i in warps]) + self.log.info("Processing %d visits", len(visits)) + + backgroundInfoList = [] + matchedImageList = [] + diffTractBackgrounds = self._makeTractBackgrounds(warps, refVisitId) + + # Blank ref warp background as reference backgrounds + im = warps[0].get() # Use arbitrary image as base + bkgd = diffTractBackgrounds[refVisitId].toWarpBackground(im) + blank = bkgd.getImage() + blank *= 0 + + for warp in warps: + visId = warp.dataId["visit"] + if visId == refVisitId: + backgroundInfoList.append(bkgd) + matchedImageList.append(warp.get()) + continue + self.log.info( + "Matching background of %s to same patch in visit %s", + warp.dataId, + refVisitId, + ) + im = warp.get() + maskIm = im.getMaskedImage() + # Matching must be done at common zeropoint + instFluxToNanojansky = self._fluxScale(im) + tractBg = diffTractBackgrounds[visId] + diffModel = tractBg.toWarpBackground(im) + bkgdIm = diffModel.getImage() + maskIm.image += bkgdIm + # Then convert everything back to counts + maskIm.image /= instFluxToNanojansky + bkgdIm /= instFluxToNanojansky + + backgroundInfoList.append(diffModel) + matchedImageList.append(im) + + return backgroundInfoList, matchedImageList + + @timeMethod + def matchBackgroundsWarps(self, refExposure, sciExposure): + """ + DEPRECATED, CURRENTLY ONLY HERE AS A QUICK REFERENCE + + Match science exposure's background level to that of reference exposure. Process creates a difference image of the reference exposure minus the @@ -603,7 +677,6 @@ def matchBackgrounds(self, refExposure, sciExposure): model : `~lsst.afw.math.BackgroundMI` Background model of difference image, reference - science """ - # Check Configs for polynomials: if self.config.usePolynomial: x, y = sciExposure.getDimensions() shortSideLength = min(x, y) @@ -619,7 +692,6 @@ def matchBackgrounds(self, refExposure, sciExposure): if self.config.order > npoints - 1: raise ValueError("%d = config.order > npoints - 1 = %d" % (self.config.order, npoints - 1)) - # Check that exposures are same shape if sciExposure.getDimensions() != refExposure.getDimensions(): wSci, hSci = sciExposure.getDimensions() wRef, hRef = refExposure.getDimensions() @@ -713,7 +785,7 @@ def matchBackgrounds(self, refExposure, sciExposure): mse, meanVar, ) - # TODO: verify this is correct (borrowed from background.py) + return BackgroundList( ( outBkgd,