DM-52108: Add simplified image differencing task

python/lsst/ip/diffim/makeKernel.py

@@ -80,7 +80,8 @@ def setDefaults(self):
         # Minimal set of measurments for star selection
         self.selectMeasurement.algorithms.names.clear()
         self.selectMeasurement.algorithms.names = ('base_SdssCentroid', 'base_PsfFlux', 'base_PixelFlags',
-                                                   'base_SdssShape', 'base_GaussianFlux', 'base_SkyCoord')
+                                                   'base_SdssShape', 'base_GaussianFlux', 'base_SkyCoord',
+                                                   'base_ClassificationSizeExtendedness')
         self.selectMeasurement.slots.modelFlux = None
         self.selectMeasurement.slots.apFlux = None
         self.selectMeasurement.slots.calibFlux = None

python/lsst/ip/diffim/subtractImages.py

@@ -37,6 +37,7 @@
 
 __all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask",
            "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask",
+           "SimplifiedSubtractConfig", "SimplifiedSubtractTask",
            "InsufficientKernelSourcesError"]
 
 _dimensions = ("instrument", "visit", "detector")
@@ -154,6 +155,24 @@ class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutp
     pass
 
 
+class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections):
+    inputPsfMatchingKernel = connectionTypes.Input(
+        doc="Kernel used to PSF match the science and template images.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="MatchingKernel",
+        name="{fakesType}{coaddName}Diff_psfMatchKernel",
+    )
+
+    def __init__(self, *, config=None):
+        super().__init__(config=config)
+        del self.sources
+        if config.useExistingKernel:
+            del self.psfMatchingKernel
+            del self.kernelSources
+        else:
+            del self.inputPsfMatchingKernel
+
+
 class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
     makeKernel = lsst.pex.config.ConfigurableField(
         target=MakeKernelTask,
@@ -209,6 +228,14 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
         target=ScienceSourceSelectorTask,
         doc="Task to select sources to be used for PSF matching.",
     )
+    fallbackSourceSelector = lsst.pex.config.ConfigurableField(
+        target=ScienceSourceSelectorTask,
+        doc="Task to select sources to be used for PSF matching."
+        "Used only if the kernel calculation fails and"
+        "`allowKernelSourceDetection` is set. The fallback source detection"
+        " will not include all of the same plugins as the original source "
+        " detection, so not all of the same flags can be used.",
+    )
     detectionThreshold = lsst.pex.config.Field(
         dtype=float,
         default=10,
@@ -282,6 +309,14 @@ def setDefaults(self):
         self.sourceSelector.doSignalToNoise = True  # apply signal to noise filter
         self.sourceSelector.signalToNoise.minimum = 10
         self.sourceSelector.signalToNoise.maximum = 500
+        self.fallbackSourceSelector.doSkySources = False  # Do not include sky sources
+        self.fallbackSourceSelector.doSignalToNoise = True  # apply signal to noise filter
+        self.fallbackSourceSelector.signalToNoise.minimum = 10
+        # The following two configs should not be necessary to be turned on for
+        # PSF-matching, and the fallback kernel source selection will fail if
+        # they are set since it does not run deblending.
+        self.fallbackSourceSelector.doIsolated = False  # Do not apply isolated star selection
+        self.fallbackSourceSelector.doRequirePrimary = False  # Do not apply primary flag selection
 
 
 class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
@@ -308,6 +343,7 @@ def __init__(self, **kwargs):
         self.makeSubtask("decorrelate")
         self.makeSubtask("makeKernel")
         self.makeSubtask("sourceSelector")
+        self.makeSubtask("fallbackSourceSelector")
         if self.config.doScaleVariance:
             self.makeSubtask("scaleVariance")
 
@@ -541,17 +577,7 @@ def runMakeKernel(self, template, science, sources, convolveTemplate=True):
                 if self.config.allowKernelSourceDetection and convolveTemplate:
                     self.log.warning("Error encountered trying to construct the matching kernel"
                                      f" Running source detection and retrying. {e}")
-                    kernelSize = self.makeKernel.makeKernelBasisList(
-                        referenceFwhmPix, targetFwhmPix)[0].getWidth()
-                    sigmaToFwhm = 2*np.log(2*np.sqrt(2))
-                    candidateList = self.makeKernel.makeCandidateList(reference, target, kernelSize,
-                                                                      candidateList=None,
-                                                                      sigma=targetFwhmPix/sigmaToFwhm)
-                    kernelSources = self.makeKernel.selectKernelSources(reference, target,
-                                                                        candidateList=candidateList,
-                                                                        preconvolved=False,
-                                                                        templateFwhmPix=referenceFwhmPix,
-                                                                        scienceFwhmPix=targetFwhmPix)
+                    kernelSources = self.runKernelSourceDetection(template, science)
                     kernelResult = self.makeKernel.run(reference, target, kernelSources,
                                                        preconvolved=False,
                                                        templateFwhmPix=referenceFwhmPix,
@@ -571,6 +597,38 @@ def runMakeKernel(self, template, science, sources, convolveTemplate=True):
                                      psfMatchingKernel=kernelResult.psfMatchingKernel,
                                      kernelSources=kernelSources)
 
+    def runKernelSourceDetection(self, template, science):
+        """Run detection on the science image and use the template mask plane
+        to reject candidate sources.
+
+        Parameters
+        ----------
+        template : `lsst.afw.image.ExposureF`
+            Template exposure, warped to match the science exposure.
+        science : `lsst.afw.image.ExposureF`
+            Science exposure to subtract from the template.
+
+        Returns
+        -------
+        kernelSources : `lsst.afw.table.SourceCatalog`
+                Sources from the input catalog to use to construct the
+                PSF-matching kernel.
+        """
+        kernelSize = self.makeKernel.makeKernelBasisList(
+            self.templatePsfSize, self.sciencePsfSize)[0].getWidth()
+        sigmaToFwhm = 2*np.log(2*np.sqrt(2))
+        candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize,
+                                                          candidateList=None,
+                                                          sigma=self.sciencePsfSize/sigmaToFwhm)
+        sources = self.makeKernel.selectKernelSources(template, science,
+                                                      candidateList=candidateList,
+                                                      preconvolved=False,
+                                                      templateFwhmPix=self.templatePsfSize,
+                                                      scienceFwhmPix=self.sciencePsfSize)
+
+        # return sources
+        return self._sourceSelector(sources, science.getBBox(), fallback=True)
+
     def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None):
         """Convolve the template image with a PSF-matching kernel and subtract
         from the science image.
@@ -850,7 +908,7 @@ def _convolveExposure(self, exposure, kernel, convolutionControl,
         else:
             return convolvedExposure[bbox]
 
-    def _sourceSelector(self, sources, bbox):
+    def _sourceSelector(self, sources, bbox, fallback=False):
         """Select sources from a catalog that meet the selection criteria.
         The selection criteria include any configured parameters of the
         `sourceSelector` subtask, as well as distance from the edge if
@@ -875,8 +933,10 @@ def _sourceSelector(self, sources, bbox):
             If there are too few sources to compute the PSF matching kernel
             remaining after source selection.
         """
-
-        selected = self.sourceSelector.selectSources(sources).selected
+        if fallback:
+            selected = self.fallbackSourceSelector.selectSources(sources).selected
+        else:
+            selected = self.sourceSelector.selectSources(sources).selected
         if self.config.restrictKernelEdgeSources:
             rejectRadius = 2*self.config.makeKernel.kernel.active.kernelSize
             bbox.grow(-rejectRadius)
@@ -1356,6 +1416,159 @@ def _shapeTest(exp1, exp2, fwhmExposureBuffer, fwhmExposureGrid):
     return xTest | yTest
 
 
+class SimplifiedSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig,
+                               pipelineConnections=SimplifiedSubtractConnections):
+    mode = lsst.pex.config.ChoiceField(
+        dtype=str,
+        default="convolveTemplate",
+        allowed={"auto": "Choose which image to convolve at runtime.",
+                 "convolveScience": "Only convolve the science image.",
+                 "convolveTemplate": "Only convolve the template image."},
+        doc="Choose which image to convolve at runtime, or require that a specific image is convolved."
+    )
+    useExistingKernel = lsst.pex.config.Field(
+        dtype=bool,
+        default=True,
+        doc="Use a pre-existing PSF matching kernel?"
+        "If False, source detection and measurement will be run."
+    )
+
+
+class SimplifiedSubtractTask(AlardLuptonSubtractTask):
+    """Compute the image difference of a science and template image using
+    the Alard & Lupton (1998) algorithm.
+    """
+    ConfigClass = SimplifiedSubtractConfig
+    _DefaultName = "simplifiedSubtract"
+
+    @timeMethod
+    def run(self, template, science, visitSummary=None, inputPsfMatchingKernel=None):
+        """PSF match, subtract, and decorrelate two images.
+
+        Parameters
+        ----------
+        template : `lsst.afw.image.ExposureF`
+            Template exposure, warped to match the science exposure.
+        science : `lsst.afw.image.ExposureF`
+            Science exposure to subtract from the template.
+        visitSummary : `lsst.afw.table.ExposureCatalog`, optional
+            Exposure catalog with external calibrations to be applied. Catalog
+            uses the detector id for the catalog id, sorted on id for fast
+            lookup.
+
+        Returns
+        -------
+        results : `lsst.pipe.base.Struct`
+            ``difference`` : `lsst.afw.image.ExposureF`
+                Result of subtracting template and science.
+            ``matchedTemplate`` : `lsst.afw.image.ExposureF`
+                Warped and PSF-matched template exposure.
+            ``backgroundModel`` : `lsst.afw.math.Function2D`
+                Background model that was fit while solving for the
+                PSF-matching kernel
+            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
+                Kernel used to PSF-match the convolved image.
+            ``kernelSources` : `lsst.afw.table.SourceCatalog`
+                Sources detected on the science image that were used to
+                construct the PSF-matching kernel.
+
+        Raises
+        ------
+        RuntimeError
+            If an unsupported convolution mode is supplied.
+        RuntimeError
+            If there are too few sources to calculate the PSF matching kernel.
+        lsst.pipe.base.NoWorkFound
+            Raised if fraction of good pixels, defined as not having NO_DATA
+            set, is less then the configured requiredTemplateFraction
+        """
+        self._prepareInputs(template, science, visitSummary=visitSummary)
+
+        convolveTemplate = self.chooseConvolutionMethod(template, science)
+
+        if self.config.useExistingKernel:
+            psfMatchingKernel = inputPsfMatchingKernel
+            backgroundModel = None
+            kernelSources = None
+        else:
+            kernelResult = self.runMakeKernel(template, science, convolveTemplate=convolveTemplate)
+            psfMatchingKernel = kernelResult.psfMatchingKernel
+            kernelSources = kernelResult.kernelSources
+            if self.config.doSubtractBackground:
+                backgroundModel = kernelResult.backgroundModel
+            else:
+                backgroundModel = None
+        if convolveTemplate:
+            subtractResults = self.runConvolveTemplate(template, science, psfMatchingKernel,
+                                                       backgroundModel=backgroundModel)
+        else:
+            subtractResults = self.runConvolveScience(template, science, psfMatchingKernel,
+                                                      backgroundModel=backgroundModel)
+        if kernelSources is not None:
+            subtractResults.kernelSources = kernelSources
+        return subtractResults
+
+    def runMakeKernel(self, template, science, convolveTemplate=True):
+        """Construct the PSF-matching kernel.
+
+        Parameters
+        ----------
+        template : `lsst.afw.image.ExposureF`
+            Template exposure, warped to match the science exposure.
+        science : `lsst.afw.image.ExposureF`
+            Science exposure to subtract from the template.
+        sources : `lsst.afw.table.SourceCatalog`
+            Identified sources on the science exposure. This catalog is used to
+            select sources in order to perform the AL PSF matching on stamp
+            images around them.
+        convolveTemplate : `bool`, optional
+            Construct the matching kernel to convolve the template?
+
+        Returns
+        -------
+        results : `lsst.pipe.base.Struct`
+            ``backgroundModel`` : `lsst.afw.math.Function2D`
+                Background model that was fit while solving for the
+                PSF-matching kernel
+            ``psfMatchingKernel`` : `lsst.afw.math.Kernel`
+                Kernel used to PSF-match the convolved image.
+            ``kernelSources` : `lsst.afw.table.SourceCatalog`
+                Sources from the input catalog that were used to construct the
+                PSF-matching kernel.
+        """
+        if convolveTemplate:
+            reference = template
+            target = science
+            referenceFwhmPix = self.templatePsfSize
+            targetFwhmPix = self.sciencePsfSize
+        else:
+            reference = science
+            target = template
+            referenceFwhmPix = self.sciencePsfSize
+            targetFwhmPix = self.templatePsfSize
+        try:
+            # The try..except block catches any error, and raises
+            # NoWorkFound if the template coverage is insufficient. Otherwise,
+            # the original error is raised.
+            kernelSources = self.runKernelSourceDetection(template, science)
+            kernelResult = self.makeKernel.run(reference, target, kernelSources,
+                                               preconvolved=False,
+                                               templateFwhmPix=referenceFwhmPix,
+                                               scienceFwhmPix=targetFwhmPix)
+        except (RuntimeError, lsst.pex.exceptions.Exception) as e:
+            self.log.warning("Failed to match template. Checking coverage")
+            #  Raise NoWorkFound if template fraction is insufficient
+            checkTemplateIsSufficient(template[science.getBBox()], science, self.log,
+                                      self.config.minTemplateFractionForExpectedSuccess,
+                                      exceptionMessage="Template coverage lower than expected to succeed."
+                                      f" Failure is tolerable: {e}")
+            #  checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception
+            raise e
+        return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel,
+                                     psfMatchingKernel=kernelResult.psfMatchingKernel,
+                                     kernelSources=kernelSources)
+
+
 def _interpolateImage(maskedImage, badMaskPlanes, fallbackValue=None):
     """Replace masked image pixels with interpolated values.