From 4de28032593a79059a8bc6b456524eb6f41431d5 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Fri, 10 Apr 2026 14:16:26 -0700
Subject: [PATCH 01/24] Support partial outputs for detection on the
 preconvolved diffim

---
 python/lsst/ip/diffim/detectAndMeasure.py | 27 +++++++++++++++--------
 1 file changed, 18 insertions(+), 9 deletions(-)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index acebadd1..19bd8c7d 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -1368,7 +1368,7 @@ class DetectAndMeasureScoreTask(DetectAndMeasureTask):
 
     @timeMethod
     def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources,
-            idFactory=None):
+            idFactory=None, measurementResults=None):
         """Detect and measure sources on a score image.
 
         Parameters
@@ -1388,6 +1388,10 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
             Generator object used to assign ids to detected sources in the
             difference image. Ids from this generator are not set until after
             deblending and merging positive/negative peaks.
+        measurementResults : `lsst.pipe.base.Struct`, optional
+            Result struct that is modified to allow saving of partial outputs
+            for some failure conditions. If the task completes successfully,
+            this is also returned.
 
         Returns
         -------
@@ -1398,6 +1402,8 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
             ``diaSources``  : `lsst.afw.table.SourceCatalog`
                 The catalog of detected sources.
         """
+        if measurementResults is None:
+            measurementResults = pipeBase.Struct()
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
@@ -1420,17 +1426,20 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
                                                           results.positive,
                                                           results.negative)
 
-            return self.processResults(science, matchedTemplate, difference,
-                                       sources, idFactory, kernelSources,
-                                       positives=positives,
-                                       negatives=negatives)
+            self.processResults(science, matchedTemplate, difference,
+                                sources, idFactory, kernelSources,
+                                positives=positives,
+                                negatives=negatives,
+                                measurementResults=measurementResults)
 
         else:
             positives = afwTable.SourceCatalog(self.schema)
             results.positive.makeSources(positives)
             negatives = afwTable.SourceCatalog(self.schema)
             results.negative.makeSources(negatives)
-            return self.processResults(science, matchedTemplate, difference,
-                                       results.sources, idFactory, kernelSources,
-                                       positives=positives,
-                                       negatives=negatives)
+            self.processResults(science, matchedTemplate, difference,
+                                results.sources, idFactory, kernelSources,
+                                positives=positives,
+                                negatives=negatives,
+                                measurementResults=measurementResults)
+        return measurementResults

From acfe806ab02caace5ed366bb1053fc8d01410294 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Tue, 21 Apr 2026 14:57:07 -0700
Subject: [PATCH 02/24] Reflect the science image PSF for preconvolution.

Also make sure to use a consistent PSF size for preconvolution and decorrelation.
---
 python/lsst/ip/diffim/makeKernel.py     | 12 +++--
 python/lsst/ip/diffim/subtractImages.py | 70 ++++++++++++++++++++-----
 2 files changed, 65 insertions(+), 17 deletions(-)

diff --git a/python/lsst/ip/diffim/makeKernel.py b/python/lsst/ip/diffim/makeKernel.py
index b20cd0ef..ad0a1d2d 100644
--- a/python/lsst/ip/diffim/makeKernel.py
+++ b/python/lsst/ip/diffim/makeKernel.py
@@ -155,9 +155,11 @@ def run(self, template, science, kernelSources, preconvolved=False,
                                                      fwhmExposureBuffer=self.config.fwhmExposureBuffer,
                                                      fwhmExposureGrid=self.config.fwhmExposureGrid
                                                      )
-
-        if preconvolved:
-            scienceFwhmPix *= np.sqrt(2)
+            # Apply the Gaussian approximation FWHM boost only when
+            # the size was auto-computed. An explicit ``scienceFwhmPix`` is
+            # assumed already to reflect the preconvolved effective width.
+            if preconvolved:
+                scienceFwhmPix *= np.sqrt(2)
         basisList = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix,
                                              metadata=self.metadata)
         spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
@@ -208,8 +210,8 @@ def selectKernelSources(self, template, science, candidateList=None, preconvolve
                                                      fwhmExposureBuffer=self.config.fwhmExposureBuffer,
                                                      fwhmExposureGrid=self.config.fwhmExposureGrid
                                                      )
-        if preconvolved:
-            scienceFwhmPix *= np.sqrt(2)
+            if preconvolved:
+                scienceFwhmPix *= np.sqrt(2)
         kernelSize = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix)[0].getWidth()
         kernelSources = self.makeCandidateList(template, science, kernelSize,
                                                candidateList=candidateList,
diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 2cebb25c..3514f6fb 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -27,7 +27,8 @@
 import lsst.afw.image
 import lsst.afw.math
 import lsst.geom
-from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, computeDifferenceImageMetrics,
+from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm,
+                                  computeDifferenceImageMetrics,
                                   checkMask, setSourceFootprints)
 from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask
 import lsst.pex.config
@@ -909,7 +910,8 @@ def _convolveExposure(self, exposure, kernel, convolutionControl,
         else:
             return convolvedExposure[bbox]
 
-    def _sourceSelector(self, template, science, sources, fallback=False, preconvolved=False):
+    def _sourceSelector(self, template, science, sources, fallback=False, preconvolved=False,
+                        scienceFwhmPix=None):
         """Select sources from a catalog that meet the selection criteria.
         The selection criteria include any configured parameters of the
         `sourceSelector` subtask, as well as checking the science and template
@@ -931,6 +933,11 @@ def _sourceSelector(self, template, science, sources, fallback=False, preconvolv
         preconvolved : `bool`, optional
             If set, exclude a wider buffer around the edge of the image to
             account for an extra convolution.
+        scienceFwhmPix : `float`, optional
+            Effective FWHM of the science image in pixels to use for sizing
+            the kernel basis and source footprints. If `None`, the
+            width is derived from ``self.sciencePsfSize`` (boosted by
+            ``sqrt(2)`` in the ``preconvolved`` case).
 
         Returns
         -------
@@ -948,8 +955,9 @@ def _sourceSelector(self, template, science, sources, fallback=False, preconvolv
             selected = self.fallbackSourceSelector.selectSources(sources).selected
         else:
             selected = self.sourceSelector.selectSources(sources).selected
-        sciencePsfSize = self.sciencePsfSize*np.sqrt(2) if preconvolved else self.sciencePsfSize
-        kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, sciencePsfSize)[0].getWidth()
+        if scienceFwhmPix is None:
+            scienceFwhmPix = self.sciencePsfSize*np.sqrt(2) if preconvolved else self.sciencePsfSize
+        kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, scienceFwhmPix)[0].getWidth()
         selectSources = sources[selected].copy(deep=True)
         # Set the footprints, to be used in `makeKernel` and `checkMask`.
         kernelSources = setSourceFootprints(selectSources, kernelSize=kSize)
@@ -1229,15 +1237,24 @@ def run(self, template, science, sources, visitSummary=None):
         """
         self._prepareInputs(template, science, visitSummary=visitSummary)
 
-        # TODO: DM-37212 we need to mirror the kernel in order to get correct cross correlation
-        scienceKernel = science.psf.getKernel()
-        matchedScience = self._convolveExposure(science, scienceKernel, self.convolutionControl,
+        # Reflect the PSF kernel in both axes so the preconvolution implements
+        # a matched-filter correlation rather than a convolution, which
+        # centers the score on the true source location even for asymmetric
+        # PSFs. The same reflected kernel is reused as ``preConvKernel`` in
+        # the decorrelation step so that both stages are self-consistent.
+        convolutionKernel = self._makePreconvolutionKernel(science.psf)
+        matchedScience = self._convolveExposure(science, convolutionKernel, self.convolutionControl,
                                                 interpolateBadMaskPlanes=True)
         self.metadata["convolvedExposure"] = "Preconvolution"
+
+        self.preconvolvedSciencePsfSize = self.sciencePsfSize*np.sqrt(2)
+        self.log.info("Preconvolved science PSF FWHM: %f pixels", self.preconvolvedSciencePsfSize)
+        self.metadata["preconvolvedSciencePsfSize"] = self.preconvolvedSciencePsfSize
         try:
-            kernelSources = self._sourceSelector(template, matchedScience, sources, preconvolved=True)
+            kernelSources = self._sourceSelector(template, matchedScience, sources, preconvolved=True,
+                                                 scienceFwhmPix=self.preconvolvedSciencePsfSize)
             subtractResults = self.runPreconvolve(template, science, matchedScience,
-                                                  kernelSources, scienceKernel)
+                                                  kernelSources, convolutionKernel)
 
         except (RuntimeError, lsst.pex.exceptions.Exception) as e:
             self.log.warning("Failed to match template. Checking coverage")
@@ -1251,6 +1268,35 @@ def run(self, template, science, sources, visitSummary=None):
 
         return subtractResults
 
+    @staticmethod
+    def _makePreconvolutionKernel(psf):
+        """Build a normalized, reflected matched-filter kernel from a PSF.
+
+        Convolving an image with this kernel is equivalent to correlating
+        the image with the PSF, so peaks in the output align with the PSF's
+        centroid — even for asymmetric PSFs. The kernel is evaluated at the
+        PSF's average position and returned as a constant
+        `~lsst.afw.math.FixedKernel`.
+
+        Parameters
+        ----------
+        psf : `~lsst.afw.detection.Psf`
+            The PSF to derive the preconvolution kernel from.
+
+        Returns
+        -------
+        kernel : `~lsst.afw.math.FixedKernel`
+            The PSF reflected about both axes, normalized to sum to one.
+        """
+        avgPos = psf.getAveragePosition()
+        localKernel = psf.getLocalKernel(avgPos)
+        kimg = lsst.afw.image.ImageD(localKernel.getDimensions())
+        localKernel.computeImage(kimg, True)  # normalize to unit sum
+        # Reflect about the kernel center. PSF kernels are odd-sized,
+        # so ``[::-1, ::-1]`` places the peak at the same pixel.
+        kimg.array[...] = kimg.array[::-1, ::-1]
+        return lsst.afw.math.FixedKernel(kimg)
+
     def runPreconvolve(self, template, science, matchedScience, kernelSources, preConvKernel):
         """Convolve the science image with its own PSF, then convolve the
         template with a matching kernel and subtract to form the Score
@@ -1269,8 +1315,8 @@ def runPreconvolve(self, template, science, matchedScience, kernelSources, preCo
             select sources in order to perform the AL PSF matching on stamp
             images around them.
         preConvKernel : `lsst.afw.math.Kernel`
-            The reflection of the kernel that was used to preconvolve the
-            `science` exposure. Must be normalized to sum to 1.
+            The kernel that was used to preconvolve the ``science``
+            exposure. Must be normalized to sum to 1.
 
         Returns
         -------
@@ -1298,7 +1344,7 @@ def runPreconvolve(self, template, science, matchedScience, kernelSources, preCo
         kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
                                            preconvolved=True,
                                            templateFwhmPix=self.templatePsfSize,
-                                           scienceFwhmPix=self.sciencePsfSize)
+                                           scienceFwhmPix=self.preconvolvedSciencePsfSize)
 
         matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
                                                  self.convolutionControl,

From 62743a9faf7988b6703186aca14795e6416746c3 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Wed, 22 Apr 2026 16:16:37 -0700
Subject: [PATCH 03/24] Raise if the kernel is not odd-sized

---
 python/lsst/ip/diffim/subtractImages.py | 17 ++++++++++++++++-
 tests/test_detectAndMeasure.py          |  4 ++--
 tests/utils.py                          | 11 +++++++----
 3 files changed, 25 insertions(+), 7 deletions(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 3514f6fb..36386cb8 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -1287,10 +1287,25 @@ def _makePreconvolutionKernel(psf):
         -------
         kernel : `~lsst.afw.math.FixedKernel`
             The PSF reflected about both axes, normalized to sum to one.
+
+        Raises
+        ------
+        RuntimeError
+            Raised if the PSF kernel has an even size along either axis.
+            The reflection ``[::-1, ::-1]`` only preserves the peak pixel
+            (and therefore the kernel center) for odd-sized kernels; an
+            even-sized kernel would be silently off-center by one pixel
+            after reflection and produce a misaligned matched filter.
         """
         avgPos = psf.getAveragePosition()
         localKernel = psf.getLocalKernel(avgPos)
-        kimg = lsst.afw.image.ImageD(localKernel.getDimensions())
+        dims = localKernel.getDimensions()
+        width, height = dims.getX(), dims.getY()
+        if width % 2 == 0 or height % 2 == 0:
+            raise RuntimeError(
+                f"Preconvolution requires an odd-sized PSF kernel, got {width}x{height}. "
+            )
+        kimg = lsst.afw.image.ImageD(dims)
         localKernel.computeImage(kimg, True)  # normalize to unit sum
         # Reflect about the kernel center. PSF kernels are odd-sized,
         # so ``[::-1, ::-1]`` places the peak at the same pixel.
diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index 70c7c5ee..38b71551 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -485,7 +485,7 @@ def test_detect_dipoles(self):
         offset = 1
         xSize = 300
         ySize = 300
-        kernelSize = 32
+        kernelSize = 31
         # Avoid placing sources near the edge for this test, so that we can
         # easily check that the correct number of sources are detected.
         templateBorderSize = kernelSize//2
@@ -1050,7 +1050,7 @@ def test_detect_dipoles(self):
         offset = 1
         xSize = 300
         ySize = 300
-        kernelSize = 32
+        kernelSize = 31
         # Avoid placing sources near the edge for this test, so that we can
         # easily check that the correct number of sources are detected.
         templateBorderSize = kernelSize//2
diff --git a/tests/utils.py b/tests/utils.py
index 80a3c0b4..47f5fce8 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -910,7 +910,7 @@ def makeFakeWcs():
 
 def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
                   noiseSeed=6, fluxLevel=500., fluxRange=2.,
-                  kernelSize=32, templateBorderSize=0,
+                  kernelSize=31, templateBorderSize=0,
                   background=None,
                   xSize=256,
                   ySize=256,
@@ -945,7 +945,7 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
     fluxRange : `float`, optional
         Range in flux amplitude of the simulated sources.
     kernelSize : `int`, optional
-        Size in pixels of the kernel for simulating sources.
+        Size in pixels of the kernel for simulating sources. Must be odd.
     templateBorderSize : `int`, optional
         Size in pixels of the image border used to pad the image.
     background : `lsst.afw.math.Chebyshev1Function2D`, optional
@@ -980,11 +980,14 @@ def makeTestImage(seed=5, nSrc=20, psfSize=2., noiseLevel=5.,
     Raises
     ------
     ValueError
-        If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths.
+        If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths,
+        or if ``kernelSize`` is even.
     """
+    if kernelSize % 2 == 0:
+        raise ValueError(f"kernelSize must be odd, got {kernelSize}.")
     # Distance from the inner edge of the bounding box to avoid placing test
     # sources in the model images.
-    bufferSize = kernelSize/2 + templateBorderSize + 1
+    bufferSize = kernelSize//2 + templateBorderSize + 1
 
     bbox = geom.Box2I(geom.Point2I(x0, y0), geom.Extent2I(xSize, ySize))
     if templateBorderSize > 0:

From 488525adef8721d8567a99e92b492b3556f5de0d Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Tue, 21 Apr 2026 15:22:19 -0700
Subject: [PATCH 04/24] Ensure that the border of the Score image has `EDGE`
 mask set

---
 python/lsst/ip/diffim/subtractImages.py | 32 +++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 36386cb8..5daf2407 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -1268,6 +1268,35 @@ def run(self, template, science, sources, visitSummary=None):
 
         return subtractResults
 
+    @staticmethod
+    def _flagScoreEdge(exposure, innerBBox):
+        """Set the ``EDGE`` mask bit on pixels outside ``innerBBox``.
+
+        Parameters
+        ----------
+        exposure : `~lsst.afw.image.Mask`
+            Exposure that will be modified in place. Must have
+            an ``EDGE`` mask plane.
+        innerBBox : `~lsst.geom.Box2I`
+            The valid inner region. Pixels
+            outside this bbox will have their ``EDGE`` bit set.
+        """
+        mask = exposure.mask
+        bbox = exposure.getBBox()
+        edgeBit = mask.getPlaneBitMask("EDGE")
+        dx0 = innerBBox.getMinX() - bbox.getMinX()
+        dx1 = bbox.getMaxX() - innerBBox.getMaxX()
+        dy0 = innerBBox.getMinY() - bbox.getMinY()
+        dy1 = bbox.getMaxY() - innerBBox.getMaxY()
+        if dy0 > 0:
+            mask.array[:dy0, :] |= edgeBit
+        if dy1 > 0:
+            mask.array[-dy1:, :] |= edgeBit
+        if dx0 > 0:
+            mask.array[:, :dx0] |= edgeBit
+        if dx1 > 0:
+            mask.array[:, -dx1:] |= edgeBit
+
     @staticmethod
     def _makePreconvolutionKernel(psf):
         """Build a normalized, reflected matched-filter kernel from a PSF.
@@ -1375,6 +1404,9 @@ def runPreconvolve(self, template, science, matchedScience, kernelSources, preCo
                                        templateMatched=True, preConvMode=True,
                                        preConvKernel=preConvKernel)
 
+        # Flag the outer ``preConvKernel/2``-wide border as EDGE.
+        self._flagScoreEdge(correctedScore, innerBBox)
+
         return lsst.pipe.base.Struct(scoreExposure=correctedScore,
                                      matchedTemplate=matchedTemplate,
                                      matchedScience=matchedScience,

From a8dfd3ffa198e69f633b61fb1934795f726e327e Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Wed, 22 Apr 2026 11:09:17 -0700
Subject: [PATCH 05/24] Add background subtraction for preconvolution

---
 python/lsst/ip/diffim/detectAndMeasure.py | 36 ++++++++++++-
 tests/test_detectAndMeasure.py            | 63 +++++++++++++++++++++++
 2 files changed, 97 insertions(+), 2 deletions(-)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 19bd8c7d..1cd3d5c9 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -1407,7 +1407,20 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
-        self._prepareInputs(scoreExposure)
+        if self.config.doSubtractBackground:
+            # Background is measured on the difference image, and then
+            # subtracted from both the difference and score images.
+            # The preconvolution kernel is normalized to 1, so the
+            # same background level applies to both.
+            background = self.subtractInitialBackground.run(difference.clone()).background
+            detectionScoreExposure = scoreExposure.clone()
+            detectionScoreExposure.image -= background.getImage()
+        else:
+            detectionScoreExposure = scoreExposure
+            background = afwMath.BackgroundList()
+
+        self._prepareInputs(detectionScoreExposure)
+
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't
@@ -1415,9 +1428,28 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
         table = afwTable.SourceTable.make(self.schema)
         results = self.detection.run(
             table=table,
-            exposure=scoreExposure,
+            exposure=detectionScoreExposure,
             doSmooth=False,
+            background=background,
         )
+
+        if self.config.doSubtractBackground:
+            # Refine the background with detected peaks flagged in the mask.
+            # The refinement is fit on the difference and also applied
+            # to the score image.
+            difference.setMask(detectionScoreExposure.mask)
+            background = self.subtractFinalBackground.run(difference).background
+            scoreExposure.image -= background.getImage()
+
+            table = afwTable.SourceTable.make(self.schema)
+            results = self.detection.run(
+                table=table,
+                exposure=scoreExposure,
+                doSmooth=False,
+                background=background,
+            )
+        measurementResults.differenceBackground = background
+
         # Copy the detection mask from the Score image to the difference image
         difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())
 
diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index 38b71551..be0a2054 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -1197,6 +1197,69 @@ def _detection_wrapper(setFlags=True):
         _detection_wrapper(setFlags=False)
         _detection_wrapper(setFlags=True)
 
+    def test_detect_transients_with_background(self):
+        """Background measured on the difference image should be subtracted
+        from the score image so that transients are still detected cleanly.
+        """
+        # Set up the simulated images
+        noiseLevel = 1.
+        staticSeed = 1
+        transientSeed = 6
+        fluxLevel = 500
+        xSize = 512
+        ySize = 512
+        x0 = 123
+        y0 = 456
+        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+                  "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
+        params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
+
+        bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
+        background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
+        scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
+                                             background=background_model, **kwargs)
+        matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
+        subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
+        scienceKernel = scienceBase.psf.getKernel()
+
+        # Configure the detection Task
+        detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
+        kwargs["seed"] = transientSeed
+        kwargs["nSrc"] = 10
+        kwargs["fluxLevel"] = 1000
+
+        def _detection_wrapper(positive=True):
+            transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
+            science = scienceBase.clone()
+            if positive:
+                science.maskedImage += transients.maskedImage
+            else:
+                science.maskedImage -= transients.maskedImage
+            difference = science.clone()
+            difference.maskedImage -= matchedTemplate.maskedImage
+            score = subtractTask._convolveExposure(difference, scienceKernel,
+                                                   subtractTask.convolutionControl)
+            # Keep a copy of the score image before detection so we can verify
+            # that the background was actually subtracted from it.
+            originalScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
+            originalDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
+            output = detectionTask.run(science, matchedTemplate, difference, score, sources)
+
+            # The background subtraction should leave a non-empty BackgroundList.
+            self.assertGreater(len(output.differenceBackground), 0)
+            # The difference image should have been modified in place by the
+            # background subtraction, so its median should shift noticeably.
+            subtractedScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
+            self.assertNotAlmostEqual(originalScoreMedian, subtractedScoreMedian, places=3)
+            subtractedDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
+            self.assertNotAlmostEqual(originalDiffimMedian, subtractedDiffimMedian, places=3)
+
+            refIds = []
+            scale = 1. if positive else -1.
+            for transient in transientSources:
+                self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
+        _detection_wrapper(positive=True)
+        _detection_wrapper(positive=False)
 
 class TestNegativePeaks(lsst.utils.tests.TestCase):
     """Tests of deblending and merging negative peaks, to test fixes for the

From 92959d84ee9a7f8d1cd41094be0612a8adb0f8c9 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Wed, 22 Apr 2026 11:09:38 -0700
Subject: [PATCH 06/24] Add cosmic ray detection for preconvolution

---
 python/lsst/ip/diffim/detectAndMeasure.py |  9 +++++
 tests/test_detectAndMeasure.py            | 47 +++++++++++++++++++++++
 2 files changed, 56 insertions(+)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 1cd3d5c9..a3e9748c 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -1421,6 +1421,11 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
 
         self._prepareInputs(detectionScoreExposure)
 
+        if self.config.doFindCosmicRays and not self.config.doSubtractBackground:
+            # Detect cosmic rays on the difference image and propagate the mask
+            # to the score image.
+            self.findAndMaskCosmicRays(difference)
+            scoreExposure.mask.array |= difference.mask.array
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't
@@ -1441,6 +1446,10 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
             background = self.subtractFinalBackground.run(difference).background
             scoreExposure.image -= background.getImage()
 
+            if self.config.doFindCosmicRays:
+                self.findAndMaskCosmicRays(difference)
+                scoreExposure.mask.array |= difference.mask.array
+
             table = afwTable.SourceTable.make(self.schema)
             results = self.detection.run(
                 table=table,
diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index be0a2054..aea15872 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -1261,6 +1261,53 @@ def _detection_wrapper(positive=True):
         _detection_wrapper(positive=True)
         _detection_wrapper(positive=False)
 
+    def test_mask_cosmic_rays(self):
+        """Cosmic rays detected on the difference image should propagate
+        to the mask of the returned (measured) exposure.
+        """
+        # Set up the simulated images
+        noiseLevel = 1.
+        staticSeed = 1
+        fluxLevel = 500
+        xSize = 400
+        ySize = 400
+        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+                  "xSize": xSize, "ySize": ySize}
+        science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
+        matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
+        subtractTask = subtractImages.AlardLuptonPreconvolveSubtractTask()
+        scienceKernel = science.psf.getKernel()
+        crMask = science.mask.getPlaneBitMask("CR")
+
+        # Configure the detection Task
+        detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
+
+        # Test that no CRs are detected when none are present.
+        transients, _ = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, nSrc=10, **kwargs)
+        science.maskedImage += transients.maskedImage
+        difference = science.clone()
+        difference.maskedImage -= matchedTemplate.maskedImage
+        score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
+        output = detectionTask.run(science, matchedTemplate, difference, score, sources)
+        crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
+        self.assertTrue(np.all(crMaskSet == 0))
+
+        crX0 = round(sources.getX()[0] - science.getX0())
+        crY0 = round(sources.getY()[0] - science.getY0())
+        crX1 = round(sources.getX()[1] - science.getX0())
+        crY1 = round(sources.getY()[1] - science.getY0())
+        # Inject CR-like shapes into the difference image. CR detection runs
+        # on the difference, and the mask should propagate to the measured
+        # exposure returned by the task.
+        difference.image.array[crX0:crX0+1, crY0:crY0+5] += 1234
+        difference.image.array[crX1:crX1+5, crY1:crY1+1] += 1234
+        score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
+        output = detectionTask.run(science, matchedTemplate, difference, score, sources)
+        crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
+        self.assertTrue(np.all(crMaskSet[crX0:crX0+1, crY0:crY0+5]))
+        self.assertTrue(np.all(crMaskSet[crX1:crX1+5, crY1:crY1+1]))
+
+
 class TestNegativePeaks(lsst.utils.tests.TestCase):
     """Tests of deblending and merging negative peaks, to test fixes for the
     various problems found on DM-48596.

From f6b826afb523af74de3ca88393fe4ff586d5339f Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Fri, 24 Apr 2026 09:59:47 -0700
Subject: [PATCH 07/24] The matchedTemplate should not be scaled a second time

---
 python/lsst/ip/diffim/subtractImages.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 5daf2407..958c7e1c 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -717,7 +717,6 @@ def runConvolveScience(self, template, science, psfMatchingKernel, backgroundMod
         # Place back on native photometric scale
         matchedScience.maskedImage /= norm
         matchedTemplate = template.clone()[bbox]
-        matchedTemplate.maskedImage /= norm
         matchedTemplate.setPhotoCalib(science.photoCalib)
 
         if backgroundModel is not None:

From 2998f82b7519521e6ab1cc5449dfef0e49daa815 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Fri, 24 Apr 2026 10:00:32 -0700
Subject: [PATCH 08/24] Reduce PSF size in tests to keep PSF fully contained in
 kernel size

---
 tests/test_subtractTask.py | 36 ++++++++++++++++++------------------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py
index dbf31271..e081c092 100644
--- a/tests/test_subtractTask.py
+++ b/tests/test_subtractTask.py
@@ -119,7 +119,7 @@ def test_incomplete_template_coverage(self):
         xSize = 400
         ySize = 400
         nSources = 80
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6, nSrc=nSources,
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6, nSrc=nSources,
                                          xSize=xSize, ySize=ySize)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, nSrc=nSources,
                                     templateBorderSize=border, doApplyCalibration=True,
@@ -155,7 +155,7 @@ def _run_and_check_coverage(template_coverage,
 
     def test_clear_template_mask(self):
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
@@ -304,7 +304,7 @@ def test_auto_convolveTemplate(self):
         the template psf is the smaller.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction(mode="convolveTemplate")
@@ -320,7 +320,7 @@ def test_auto_convolveScience(self):
         """
         noiseLevel = 1.
         science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction(mode="convolveScience")
         output = task.run(template.clone(), science.clone(), sources)
@@ -338,7 +338,7 @@ def test_science_better(self):
 
         def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
             science, sources = makeTestImage(psfSize=2.0, noiseLevel=scienceNoiseLevel, noiseSeed=6)
-            template, _ = makeTestImage(psfSize=3.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
+            template, _ = makeTestImage(psfSize=2.3, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                         templateBorderSize=20, doApplyCalibration=True)
             task = self._setup_subtraction(mode="convolveScience")
             output = task.run(template, science, sources)
@@ -372,7 +372,7 @@ def test_template_better(self):
         statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
 
         def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
-            science, sources = makeTestImage(psfSize=3.0, noiseLevel=scienceNoiseLevel, noiseSeed=6)
+            science, sources = makeTestImage(psfSize=2.4, noiseLevel=scienceNoiseLevel, noiseSeed=6)
             template, _ = makeTestImage(psfSize=2.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                         templateBorderSize=20, doApplyCalibration=True)
             task = self._setup_subtraction()
@@ -410,7 +410,7 @@ def test_symmetry(self):
         #  comparable when we swap which image is treated as the "science" image.
         science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel,
                                          noiseSeed=6, templateBorderSize=0, doApplyCalibration=True)
-        template, _ = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel,
+        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel,
                                     noiseSeed=7, templateBorderSize=0, doApplyCalibration=True)
         task = self._setup_subtraction(mode='auto')
 
@@ -640,7 +640,7 @@ def _run_and_check_images(science, template, sources, statsCtrl,
             varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
             self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
 
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=scienceNoiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=scienceNoiseLevel, noiseSeed=6)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         # Verify that the variance plane of the difference image is correct
@@ -710,7 +710,7 @@ def _run_and_check_images(science, template, sources, statsCtrl,
             self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
 
         science, sources = makeTestImage(psfSize=2.0, noiseLevel=scienceNoiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=3.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=2.4, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         # Verify that the variance plane of the difference image is correct
         #  when the template and science variance planes are correct
@@ -743,7 +743,7 @@ def test_exposure_properties_convolve_template(self):
         noiseLevel = 1.
         seed = 37
         rng = np.random.RandomState(seed)
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
         psf = science.psf
         psfAvgPos = psf.getAveragePosition()
         psfSize = getPsfFwhm(science.psf)
@@ -791,7 +791,7 @@ def test_exposure_properties_convolve_science(self):
         seed = 37
         rng = np.random.RandomState(seed)
         science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         psf = template.psf
         psfAvgPos = psf.getAveragePosition()
@@ -1044,7 +1044,7 @@ def test_incomplete_template_coverage(self):
         border = 20
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=border, doApplyCalibration=True,
@@ -1077,7 +1077,7 @@ def test_clear_template_mask(self):
         noiseLevel = 1.
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True,
@@ -1125,7 +1125,7 @@ def test_agnostic_template_psf(self):
         science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel,
                                          noiseSeed=6, templateBorderSize=0,
                                          xSize=xSize, ySize=ySize)
-        template1, _ = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel,
+        template1, _ = makeTestImage(psfSize=2.8, noiseLevel=noiseLevel,
                                      noiseSeed=7, doApplyCalibration=True,
                                      xSize=xSize, ySize=ySize)
         template2, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel,
@@ -1266,7 +1266,7 @@ def _run_and_check_images(science, template, sources, statsCtrl,
                                               statsCtrl)
             self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
 
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=scienceNoiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=scienceNoiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True,
@@ -1302,7 +1302,7 @@ def test_exposure_properties(self):
         noiseLevel = 1.
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
         psf = science.psf
         psfAvgPos = psf.getAveragePosition()
@@ -1344,7 +1344,7 @@ def test_runSimplifiedTaskWithExistingKernel(self):
         `AlardLuptonSubtractTask` if it uses the AL kernel.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         alTask = AlardLuptonSubtractTest._setup_subtraction(AlardLuptonSubtractTest())
@@ -1361,7 +1361,7 @@ def test_runSimplifiedTaskWithSourceDetection(self):
         reasonable output.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=3.0, noiseLevel=noiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
         template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction(useExistingKernel=False,

From b28cfbc7def54c50202c56c0454c0fd700e94716 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Fri, 24 Apr 2026 16:27:20 -0700
Subject: [PATCH 09/24] Fix discontinuity in kernel basis

Otherwise the results are unstable when the science image has a psfSize of 2.0 in the unit tests and the template has a psfSize of 2.3 to 2.4.  This change makes the basis a smooth transition instead of a sudden jump.
---
 python/lsst/ip/diffim/makeKernelBasisList.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/python/lsst/ip/diffim/makeKernelBasisList.py b/python/lsst/ip/diffim/makeKernelBasisList.py
index 8e9efabd..a686b85c 100644
--- a/python/lsst/ip/diffim/makeKernelBasisList.py
+++ b/python/lsst/ip/diffim/makeKernelBasisList.py
@@ -291,11 +291,14 @@ def _calculateBasisSigmas(referenceSigma, targetSigma, basisMinSigma, basisGauss
     if kernelSigma < basisMinSigma:
         kernelSigma = basisMinSigma
 
-    # If more than one gaussian is requested and kernelSigma is more than a
-    #  factor of basisGaussBeta larger than the minimum sigma, use
-    #  kernelSigma/basisGaussBeta as the size of the first gaussian.
-    if ((kernelSigma/basisGaussBeta) > basisMinSigma) & (basisNGauss > 1):
-        basisSigmaGauss = [kernelSigma/basisGaussBeta, ]
+    # The smallest basis Gaussian sigma is kernelSigma/basisGaussBeta to center
+    # the geometric basis spacing around kernelSigma, but bounded below by
+    # basisMinSigma. Using max() here (rather than a kernelSigma-vs-threshold
+    # branch) keeps the basis continuous as kernelSigma crosses
+    # basisMinSigma*basisGaussBeta, avoiding a discontinuous jump in the AL
+    # basis that produced unstable kernels.
+    if basisNGauss > 1:
+        basisSigmaGauss = [max(basisMinSigma, kernelSigma/basisGaussBeta), ]
     else:
         basisSigmaGauss = [kernelSigma, ]
 

From 8119893a1d44e61215973b5eb709fbbd94351496 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sat, 25 Apr 2026 11:43:15 -0700
Subject: [PATCH 10/24] Fix broken error message

---
 python/lsst/ip/diffim/subtractImages.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 958c7e1c..cb4ffe57 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -525,7 +525,7 @@ def chooseConvolutionMethod(self, template, science):
             self.log.info("`convolveScience` is set: convolving science image.")
             convolveTemplate = False
         else:
-            raise RuntimeError("Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
+            raise RuntimeError(f"Cannot handle AlardLuptonSubtract mode: {self.config.mode}")
         return convolveTemplate
 
     def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False):

From 61cc6fd82413a89f69b68a3154d47aa64799fe17 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sat, 25 Apr 2026 11:48:01 -0700
Subject: [PATCH 11/24] Evaluate psf matching kernel norm in the center of the
 image, not default (0,0)

---
 python/lsst/ip/diffim/psfMatch.py       | 7 ++++++-
 python/lsst/ip/diffim/subtractImages.py | 3 ++-
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/python/lsst/ip/diffim/psfMatch.py b/python/lsst/ip/diffim/psfMatch.py
index 96ccb33b..7c46f9e4 100644
--- a/python/lsst/ip/diffim/psfMatch.py
+++ b/python/lsst/ip/diffim/psfMatch.py
@@ -595,7 +595,12 @@ def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
         """
         # What is the final kernel sum
         kImage = afwImage.ImageD(spatialKernel.getDimensions())
-        kSum = spatialKernel.computeImage(kImage, False)
+        # Evaluate the kernel at the cell-set center, not at the world
+        # origin (the default (0, 0) lies outside the cellset bbox for any
+        # non-origin region and produces a meaningless kernel sum).
+        regionBBox = kernelCellSet.getBBox()
+        xcen, ycen = regionBBox.getCenter()
+        kSum = spatialKernel.computeImage(kImage, False, xcen, ycen)
         self.log.info("Final spatial kernel sum %.3f", kSum)
 
         # Look at how well conditioned the matrix is
diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index cb4ffe57..6d9b7e73 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -708,7 +708,8 @@ def runConvolveScience(self, template, science, psfMatchingKernel, backgroundMod
         bbox = science.getBBox()
 
         kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions())
-        norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False)
+        xcen, ycen = bbox.getCenter()
+        norm = psfMatchingKernel.computeImage(kernelImage, False, xcen, ycen)
 
         matchedScience = self._convolveExposure(science, psfMatchingKernel,
                                                 self.convolutionControl,

From 0fc71fd06ea37d985a8ae65e53e8775aec6a920f Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sat, 25 Apr 2026 11:49:02 -0700
Subject: [PATCH 12/24] Clone the background model before inverting

---
 python/lsst/ip/diffim/subtractImages.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 6d9b7e73..081ad3e4 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -721,9 +721,10 @@ def runConvolveScience(self, template, science, psfMatchingKernel, backgroundMod
         matchedTemplate.setPhotoCalib(science.photoCalib)
 
         if backgroundModel is not None:
-            modelParams = backgroundModel.getParameters()
             # We must invert the background model if the matching kernel is solved for the science image.
-            backgroundModel.setParameters([-p for p in modelParams])
+            invertedBackground = backgroundModel.clone()
+            invertedBackground.setParameters([-p for p in backgroundModel.getParameters()])
+            backgroundModel = invertedBackground
 
         difference = _subtractImages(matchedScience, matchedTemplate, backgroundModel=backgroundModel)
 

From 2499cd07d5d85e84d68fc55de01eda53dbe8fdc1 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sat, 25 Apr 2026 11:50:13 -0700
Subject: [PATCH 13/24] Errors in calculating limiting magnitudes were being
 silently squashed

---
 python/lsst/ip/diffim/subtractImages.py | 20 +++++++-------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 081ad3e4..9aba1087 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -809,29 +809,23 @@ def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None):
         """
         if exposure.photoCalib is None:
             return np.nan
-        # Set maglim to nan upfront in case on an unexpected RuntimeError
-        maglim = np.nan
         try:
             psf = exposure.getPsf()
             psf_shape = psf.computeShape(psf.getAveragePosition())
         except (lsst.pex.exceptions.InvalidParameterError,
                 afwDetection.InvalidPsfError,
                 lsst.pex.exceptions.RangeError):
-            if fallbackPsfSize is not None:
-                self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
-                psf_area = np.pi*(fallbackPsfSize/2)**2
-                zeropoint = exposure.photoCalib.instFluxToMagnitude(1)
-                maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
-            else:
+            if fallbackPsfSize is None:
                 self.log.info("Unable to evaluate PSF, setting maglim to nan")
-                maglim = np.nan
+                return np.nan
+            self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize)
+            psf_area = np.pi*(fallbackPsfSize/2)**2
         else:
             # Get a more accurate area than `psf_shape.getArea()` via moments
             psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy())
-            zeropoint = exposure.photoCalib.instFluxToMagnitude(1)
-            maglim = zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
-        finally:
-            return maglim
+
+        zeropoint = exposure.photoCalib.instFluxToMagnitude(1)
+        return zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area))
 
     @staticmethod
     def _validateExposures(template, science):

From 2e197dfde03b4b86f2850dc31ca6a21c77d8d564 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sat, 25 Apr 2026 12:47:25 -0700
Subject: [PATCH 14/24] Fix docstring

---
 python/lsst/ip/diffim/utils.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py
index c7bc12b3..debed678 100644
--- a/python/lsst/ip/diffim/utils.py
+++ b/python/lsst/ip/diffim/utils.py
@@ -502,17 +502,23 @@ def checkMask(mask, sources, excludeMaskPlanes):
 def setSourceFootprints(sources, kernelSize):
     """Add footprints of fixed size to a source catalog
 
+    Each source's footprint is replaced by a centered box of side
+    ``2*kernelSize + 1`` pixels (always odd), preserving the brightest peak
+    of the original footprint.
+
     Parameters
     ----------
     sources : `lsst.afw.table.SourceCatalog`
-        The source catalog to add footprints to.
+        The source catalog to add footprints to. Modified in place.
     kernelSize : `int`
-        The "radius" of the footprint, i.e half the size of the bounding box.
+        Width (in pixels) of the PSF-matching kernel. The footprint side is
+        ``2*kernelSize + 1`` so that the kernel can be convolved across the
+        source without falling off the stamp.
 
     Returns
     -------
     sources : `lsst.afw.table.SourceCatalog`
-        The modified source catalog
+        The same (modified) source catalog, returned for convenience.
     """
     size = 2*kernelSize + 1
     for source in sources:

From 217a9467b4d23ecf12505997a5e24a68eec66055 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Mon, 27 Apr 2026 13:54:41 -0700
Subject: [PATCH 15/24] Persist the Score image post-detection and BG
 subtraction

---
 python/lsst/ip/diffim/detectAndMeasure.py | 7 +++++++
 python/lsst/ip/diffim/subtractImages.py   | 2 +-
 2 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index a3e9748c..18d1594a 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -1339,6 +1339,12 @@ def _runStreakMasking(self, difference):
 
 class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
     scoreExposure = pipeBase.connectionTypes.Input(
+        doc="Maximum likelihood image for detection.",
+        dimensions=("instrument", "visit", "detector"),
+        storageClass="ExposureF",
+        name="{fakesType}{coaddName}Diff_scoreTempExp",
+    )
+    scoreMeasuredExposure = pipeBase.connectionTypes.Output(
         doc="Maximum likelihood image for detection.",
         dimensions=("instrument", "visit", "detector"),
         storageClass="ExposureF",
@@ -1458,6 +1464,7 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
                 background=background,
             )
         measurementResults.differenceBackground = background
+        measurementResults.scoreMeasuredExposure = scoreExposure
 
         # Copy the detection mask from the Score image to the difference image
         difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())
diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 9aba1087..2a4d3291 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -138,7 +138,7 @@ class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections,
         doc="The maximum likelihood image, used for the detection of diaSources.",
         dimensions=("instrument", "visit", "detector"),
         storageClass="ExposureF",
-        name="{fakesType}{coaddName}Diff_scoreExp",
+        name="{fakesType}{coaddName}Diff_scoreTempExp",
     )
     psfMatchingKernel = connectionTypes.Output(
         doc="Kernel used to PSF match the science and template images.",

From cd5766b3a47614b36899e428d714a47e1a09463b Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Mon, 27 Apr 2026 16:35:34 -0700
Subject: [PATCH 16/24] Add option to reject bad peaks prior to measurement.

This does not suffer from the same bug as setting the image plane to 0 for masked pixels, since the peak is removed after detection.
---
 python/lsst/ip/diffim/detectAndMeasure.py | 67 +++++++++++++++++++++++
 tests/test_detectAndMeasure.py            | 28 +++++++---
 2 files changed, 86 insertions(+), 9 deletions(-)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 18d1594a..4515e72b 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -346,6 +346,19 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
         doc="Mask planes to clear before running detection.",
         default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"),
     )
+    doRejectBadMaskPlaneDetections = pexConfig.Field(
+        dtype=bool,
+        default=True,
+        doc="Reject any peaks detected on ``badMaskPlanes`` before measurement."
+            "These should all be rejected downstream by ``badSourceFlags``, "
+            "but filtering earlier can save time."
+    )
+    badMaskPlanes = lsst.pex.config.ListField(
+        dtype=str,
+        doc="Detections whose footprint peak lies on a pixel with any of these "
+            "mask planes set will be rejected before measurement.",
+        default=("NO_DATA", "BAD", "SAT", "EDGE"),
+    )
     raiseOnBadSubtractionRatio = pexConfig.Field(
         dtype=bool,
         default=True,
@@ -823,6 +836,11 @@ def processResults(self, science, matchedTemplate, difference, sources, idFactor
         if self.config.doSkySources:
             self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
 
+        if self.config.doRejectBadMaskPlaneDetections:
+            # Save time by rejecting peaks on bad mask planes prior to
+            # measurement.
+            initialDiaSources = self._rejectBadMaskedDetections(initialDiaSources, difference.mask)
+
         if not initialDiaSources.isContiguous():
             initialDiaSources = initialDiaSources.copy(deep=True)
 
@@ -940,6 +958,55 @@ def deblend(footprints, negative=False):
             sources[-len(negatives):]["is_negative"] = True
         return sources, positives, negatives
 
+    def _rejectBadMaskedDetections(self, initialDiaSources, mask):
+        """Remove detections whose footprint peak lies on a bad-masked pixel.
+
+        Detections are rejected if any peak of the source footprint falls on a
+        pixel with any of the ``config.badMaskPlanes`` bits set.
+
+        Parameters
+        ----------
+        initialDiaSources : `lsst.afw.table.SourceCatalog`
+            The catalog of detected sources, before measurement.
+        mask : `lsst.afw.image.Mask`
+            Mask of the image used for detection.
+
+        Returns
+        -------
+        initialDiaSources : `lsst.afw.table.SourceCatalog`
+            The catalog with bad-masked detections removed.
+        """
+        if not self.config.badMaskPlanes or len(initialDiaSources) == 0:
+            self.metadata["nRejectedBadMaskPlanes"] = 0
+            return initialDiaSources
+
+        presentPlanes = [mp for mp in self.config.badMaskPlanes
+                         if mp in mask.getMaskPlaneDict()]
+        if not presentPlanes:
+            self.metadata["nRejectedBadMaskPlanes"] = 0
+            return initialDiaSources
+
+        badBitMask = mask.getPlaneBitMask(presentPlanes)
+        x0 = mask.getX0()
+        y0 = mask.getY0()
+
+        selector = np.ones(len(initialDiaSources), dtype=bool)
+        for i, source in enumerate(initialDiaSources):
+            peaks = source.getFootprint().getPeaks()
+            for peak in peaks:
+                ix = peak.getIx() - x0
+                iy = peak.getIy() - y0
+                if mask.array[iy, ix] & badBitMask:
+                    selector[i] = False
+                    break
+        nRejected = np.count_nonzero(~selector)
+        self.metadata["nRejectedBadMaskPlanes"] = nRejected
+        if nRejected > 0:
+            self.log.info("Rejected %d detections on pixels with bad mask "
+                          "planes %s before measurement.",
+                          nRejected, list(presentPlanes))
+        return initialDiaSources[selector].copy(deep=True)
+
     def _removeBadSources(self, diaSources):
         """Remove unphysical diaSources from the catalog.
 
diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index aea15872..ab32f4da 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -903,7 +903,13 @@ def test_detection_xy0(self):
         noiseLevel = 1.
         staticSeed = 1
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
+        kernelSize = 31
+        # Buffer source positions away from the score image's EDGE strip.
+        # Preconvolution adds ~kernelSize//2 pixels of EDGE on top of the
+        # EDGE already set by detection smoothing on the science image.
+        templateBorderSize = kernelSize//2
+        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
+                  "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
         difference = science.clone()
@@ -930,12 +936,9 @@ def test_detection_xy0(self):
         self.assertEqual(len(output.diaSources), len(sources))
         # no sources should be flagged as negative
         self.assertEqual(len(~output.diaSources["is_negative"]), len(output.diaSources))
-        # TODO DM-41496: restore this block once we handle detections on edge
-        # pixels better; at least one of these sources currently has a bad
-        # centroid because most of the source is rejected as EDGE.
-        # refIds = []
-        # for source in sources:
-        #     self._check_diaSource(output.diaSources, source, refIds=refIds)
+        refIds = []
+        for source in sources:
+            self._check_diaSource(output.diaSources, source, refIds=refIds)
 
     def test_measurements_finite(self):
         """Measured fluxes and centroids should always be finite.
@@ -1210,11 +1213,18 @@ def test_detect_transients_with_background(self):
         ySize = 512
         x0 = 123
         y0 = 456
+        kernelSize = 31
+        templateBorderSize = kernelSize//2
         kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
-                  "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
+                  "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0,
+                  "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
         params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
 
-        bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
+        # The background model must cover the grown image bbox, otherwise it
+        # would be extrapolated outside its declared domain.
+        bbox2D = lsst.geom.Box2D(
+            lsst.geom.Point2D(x0 - templateBorderSize, y0 - templateBorderSize),
+            lsst.geom.Extent2D(xSize + 2*templateBorderSize, ySize + 2*templateBorderSize))
         background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
         scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
                                              background=background_model, **kwargs)

From 9ce1be20dbba30f55ab15d26200e8e58080c6712 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Tue, 28 Apr 2026 17:37:23 -0700
Subject: [PATCH 17/24] Use smaller template coverage in test that expects
 kernel failure.

The larger coverage fraction could actually sometimes succeed on x86 machines, leading to a test failure since it unexpectedly did not raise `NoWorkFound`.
This change is necessary following the tightening of the template PSF size used in the tests, which changed the kernel behavior slightly.
---
 tests/test_subtractTask.py | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py
index e081c092..728ae362 100644
--- a/tests/test_subtractTask.py
+++ b/tests/test_subtractTask.py
@@ -150,7 +150,7 @@ def _run_and_check_coverage(template_coverage,
             else:
                 task.run(template_cut, science.clone(), sources.copy(deep=True))
         _run_and_check_coverage(template_coverage=0.09)
-        _run_and_check_coverage(template_coverage=0.19)
+        _run_and_check_coverage(template_coverage=0.15)
         _run_and_check_coverage(template_coverage=0.7)
 
     def test_clear_template_mask(self):
@@ -1052,12 +1052,17 @@ def test_incomplete_template_coverage(self):
 
         science_height = science.getBBox().getDimensions().getY()
 
-        def _run_and_check_coverage(template_coverage):
+        def _run_and_check_coverage(template_coverage,
+                                    requiredTemplateFraction=0.1,
+                                    minTemplateFractionForExpectedSuccess=0.2):
             template_cut = template.clone()
             template_height = int(science_height*template_coverage + border)
             template_cut.image.array[:, template_height:] = 0
             template_cut.mask.array[:, template_height:] = template_cut.mask.getPlaneBitMask('NO_DATA')
-            task = self._setup_subtraction()
+            task = self._setup_subtraction(
+                requiredTemplateFraction=requiredTemplateFraction,
+                minTemplateFractionForExpectedSuccess=minTemplateFractionForExpectedSuccess
+            )
             if template_coverage < task.config.requiredTemplateFraction:
                 doRaise = True
             elif template_coverage < task.config.minTemplateFractionForExpectedSuccess:
@@ -1070,7 +1075,7 @@ def _run_and_check_coverage(template_coverage):
             else:
                 task.run(template_cut, science.clone(), sources.copy(deep=True))
         _run_and_check_coverage(template_coverage=0.09)
-        _run_and_check_coverage(template_coverage=0.19)
+        _run_and_check_coverage(template_coverage=0.15)
         _run_and_check_coverage(template_coverage=.7)
 
     def test_clear_template_mask(self):

From be2d8d73389b102a126c771ac0b12c44ff4fc849 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 18:20:02 -0700
Subject: [PATCH 18/24] Combine all fixup commits from review

---
 python/lsst/ip/diffim/detectAndMeasure.py    | 14 +++++-----
 python/lsst/ip/diffim/makeKernelBasisList.py |  6 ++--
 python/lsst/ip/diffim/subtractImages.py      | 27 ++++++------------
 python/lsst/ip/diffim/utils.py               |  5 ++--
 tests/test_detectAndMeasure.py               | 29 +++++++++++---------
 tests/test_utils.py                          |  5 ++++
 6 files changed, 42 insertions(+), 44 deletions(-)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 4515e72b..0009b613 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -355,8 +355,9 @@ class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
     )
     badMaskPlanes = lsst.pex.config.ListField(
         dtype=str,
-        doc="Detections whose footprint peak lies on a pixel with any of these "
-            "mask planes set will be rejected before measurement.",
+        doc="Detections whose footprint peak lies on a pixel with any of these"
+            " mask planes set will be rejected before measurement."
+            " Any missing mask planes will be silently ignored.",
         default=("NO_DATA", "BAD", "SAT", "EDGE"),
     )
     raiseOnBadSubtractionRatio = pexConfig.Field(
@@ -967,7 +968,7 @@ def _rejectBadMaskedDetections(self, initialDiaSources, mask):
         Parameters
         ----------
         initialDiaSources : `lsst.afw.table.SourceCatalog`
-            The catalog of detected sources, before measurement.
+            The catalog of detected peaks.
         mask : `lsst.afw.image.Mask`
             Mask of the image used for detection.
 
@@ -987,8 +988,7 @@ def _rejectBadMaskedDetections(self, initialDiaSources, mask):
             return initialDiaSources
 
         badBitMask = mask.getPlaneBitMask(presentPlanes)
-        x0 = mask.getX0()
-        y0 = mask.getY0()
+        x0, y0 = mask.getXY0()
 
         selector = np.ones(len(initialDiaSources), dtype=bool)
         for i, source in enumerate(initialDiaSources):
@@ -1004,7 +1004,7 @@ def _rejectBadMaskedDetections(self, initialDiaSources, mask):
         if nRejected > 0:
             self.log.info("Rejected %d detections on pixels with bad mask "
                           "planes %s before measurement.",
-                          nRejected, list(presentPlanes))
+                          nRejected, presentPlanes)
         return initialDiaSources[selector].copy(deep=True)
 
     def _removeBadSources(self, diaSources):
@@ -1412,7 +1412,7 @@ class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
         name="{fakesType}{coaddName}Diff_scoreTempExp",
     )
     scoreMeasuredExposure = pipeBase.connectionTypes.Output(
-        doc="Maximum likelihood image for detection.",
+        doc="Score image after backgrond subtraction and detection.",
         dimensions=("instrument", "visit", "detector"),
         storageClass="ExposureF",
         name="{fakesType}{coaddName}Diff_scoreExp",
diff --git a/python/lsst/ip/diffim/makeKernelBasisList.py b/python/lsst/ip/diffim/makeKernelBasisList.py
index a686b85c..e212ec49 100644
--- a/python/lsst/ip/diffim/makeKernelBasisList.py
+++ b/python/lsst/ip/diffim/makeKernelBasisList.py
@@ -293,10 +293,8 @@ def _calculateBasisSigmas(referenceSigma, targetSigma, basisMinSigma, basisGauss
 
     # The smallest basis Gaussian sigma is kernelSigma/basisGaussBeta to center
     # the geometric basis spacing around kernelSigma, but bounded below by
-    # basisMinSigma. Using max() here (rather than a kernelSigma-vs-threshold
-    # branch) keeps the basis continuous as kernelSigma crosses
-    # basisMinSigma*basisGaussBeta, avoiding a discontinuous jump in the AL
-    # basis that produced unstable kernels.
+    # basisMinSigma. Keep the larger of the two so that we don't have a
+    # discontinuous jump in the AL basis and unstable kernel solutions.
     if basisNGauss > 1:
         basisSigmaGauss = [max(basisMinSigma, kernelSigma/basisGaussBeta), ]
     else:
diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 2a4d3291..18b15b76 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -1232,11 +1232,6 @@ def run(self, template, science, sources, visitSummary=None):
         """
         self._prepareInputs(template, science, visitSummary=visitSummary)
 
-        # Reflect the PSF kernel in both axes so the preconvolution implements
-        # a matched-filter correlation rather than a convolution, which
-        # centers the score on the true source location even for asymmetric
-        # PSFs. The same reflected kernel is reused as ``preConvKernel`` in
-        # the decorrelation step so that both stages are self-consistent.
         convolutionKernel = self._makePreconvolutionKernel(science.psf)
         matchedScience = self._convolveExposure(science, convolutionKernel, self.convolutionControl,
                                                 interpolateBadMaskPlanes=True)
@@ -1265,11 +1260,11 @@ def run(self, template, science, sources, visitSummary=None):
 
     @staticmethod
     def _flagScoreEdge(exposure, innerBBox):
-        """Set the ``EDGE`` mask bit on pixels outside ``innerBBox``.
+        """Set the EDGE mask bit on pixels outside a known-valid region.
 
         Parameters
         ----------
-        exposure : `~lsst.afw.image.Mask`
+        exposure : `~lsst.afw.image.Exposure`
             Exposure that will be modified in place. Must have
             an ``EDGE`` mask plane.
         innerBBox : `~lsst.geom.Box2I`
@@ -1300,7 +1295,7 @@ def _makePreconvolutionKernel(psf):
         the image with the PSF, so peaks in the output align with the PSF's
         centroid — even for asymmetric PSFs. The kernel is evaluated at the
         PSF's average position and returned as a constant
-        `~lsst.afw.math.FixedKernel`.
+        `~lsst.afw.math.Kernel`.
 
         Parameters
         ----------
@@ -1309,25 +1304,21 @@ def _makePreconvolutionKernel(psf):
 
         Returns
         -------
-        kernel : `~lsst.afw.math.FixedKernel`
+        kernel : `~lsst.afw.math.Kernel`
             The PSF reflected about both axes, normalized to sum to one.
 
         Raises
         ------
-        RuntimeError
+        ValueError
             Raised if the PSF kernel has an even size along either axis.
-            The reflection ``[::-1, ::-1]`` only preserves the peak pixel
-            (and therefore the kernel center) for odd-sized kernels; an
-            even-sized kernel would be silently off-center by one pixel
-            after reflection and produce a misaligned matched filter.
+            It's not possible to center an even-sized kernel.
         """
         avgPos = psf.getAveragePosition()
         localKernel = psf.getLocalKernel(avgPos)
         dims = localKernel.getDimensions()
-        width, height = dims.getX(), dims.getY()
-        if width % 2 == 0 or height % 2 == 0:
-            raise RuntimeError(
-                f"Preconvolution requires an odd-sized PSF kernel, got {width}x{height}. "
+        if dims.x % 2 == 0 or dims.y % 2 == 0:
+            raise ValueError(
+                f"Preconvolution requires an odd-sized PSF kernel, got {dims.x}x{dims.y}. "
             )
         kimg = lsst.afw.image.ImageD(dims)
         localKernel.computeImage(kimg, True)  # normalize to unit sum
diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py
index debed678..5d902cdd 100644
--- a/python/lsst/ip/diffim/utils.py
+++ b/python/lsst/ip/diffim/utils.py
@@ -511,14 +511,15 @@ def setSourceFootprints(sources, kernelSize):
     sources : `lsst.afw.table.SourceCatalog`
         The source catalog to add footprints to. Modified in place.
     kernelSize : `int`
-        Width (in pixels) of the PSF-matching kernel. The footprint side is
+        Width (in pixels) of the PSF-matching kernel. The footprint size is
         ``2*kernelSize + 1`` so that the kernel can be convolved across the
         source without falling off the stamp.
 
     Returns
     -------
     sources : `lsst.afw.table.SourceCatalog`
-        The same (modified) source catalog, returned for convenience.
+        The input `sources` catalog that was modified in place,
+        returned for convenience.
     """
     size = 2*kernelSize + 1
     for source in sources:
diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index ab32f4da..487091c2 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -447,12 +447,14 @@ def test_mask_cosmic_rays(self):
         # Add CR-like shape and check that CR is detected
         # Pick two locations on top of sources, since that is what is likely to
         # be missed in the first stage of CR rejection.
-        difference.image.array[crX0:crX0+1, crY0:crY0+5] += 1234  # Arbitrary but large flux for the CRs
-        difference.image.array[crX1:crX1+5, crY1:crY1+1] += 1234
+        crMaskSetInput = np.zeros(difference.image.array.shape, bool)
+        crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
+        crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
+        difference.image.array[crMaskSetInput] += 1234
         output = detectionTask.run(science, matchedTemplate, difference)
         crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
-        self.assertTrue(np.all(crMaskSet[crX0:crX0+1, crY0:crY0+5]))
-        self.assertTrue(np.all(crMaskSet[crX1:crX1+5, crY1:crY1+1]))
+        self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
+        self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
 
     def test_missing_mask_planes(self):
         """Check that detection runs with missing mask planes.
@@ -1249,7 +1251,7 @@ def _detection_wrapper(positive=True):
             difference.maskedImage -= matchedTemplate.maskedImage
             score = subtractTask._convolveExposure(difference, scienceKernel,
                                                    subtractTask.convolutionControl)
-            # Keep a copy of the score image before detection so we can verify
+            # Record the score image metric before detection so we can verify
             # that the background was actually subtracted from it.
             originalScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
             originalDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
@@ -1260,9 +1262,9 @@ def _detection_wrapper(positive=True):
             # The difference image should have been modified in place by the
             # background subtraction, so its median should shift noticeably.
             subtractedScoreMedian = np.median(score.image.array[np.isfinite(score.image.array)])
-            self.assertNotAlmostEqual(originalScoreMedian, subtractedScoreMedian, places=3)
+            self.assertLess(np.abs(subtractedScoreMedian), np.abs(originalScoreMedian))
             subtractedDiffimMedian = np.median(difference.image.array[np.isfinite(difference.image.array)])
-            self.assertNotAlmostEqual(originalDiffimMedian, subtractedDiffimMedian, places=3)
+            self.assertLess(np.abs(subtractedDiffimMedian), np.abs(originalDiffimMedian))
 
             refIds = []
             scale = 1. if positive else -1.
@@ -1299,8 +1301,7 @@ def test_mask_cosmic_rays(self):
         difference.maskedImage -= matchedTemplate.maskedImage
         score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
         output = detectionTask.run(science, matchedTemplate, difference, score, sources)
-        crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
-        self.assertTrue(np.all(crMaskSet == 0))
+        self.assertFalse(np.any((output.subtractedMeasuredExposure.mask.array & crMask) > 0))
 
         crX0 = round(sources.getX()[0] - science.getX0())
         crY0 = round(sources.getY()[0] - science.getY0())
@@ -1309,13 +1310,15 @@ def test_mask_cosmic_rays(self):
         # Inject CR-like shapes into the difference image. CR detection runs
         # on the difference, and the mask should propagate to the measured
         # exposure returned by the task.
-        difference.image.array[crX0:crX0+1, crY0:crY0+5] += 1234
-        difference.image.array[crX1:crX1+5, crY1:crY1+1] += 1234
+        crMaskSetInput = np.zeros(difference.image.array.shape, bool)
+        crMaskSetInput[crX0:crX0+1, crY0:crY0+5] = True
+        crMaskSetInput[crX1:crX1+5, crY1:crY1+1] = True
+        difference.image.array[crMaskSetInput] += 1234
         score = subtractTask._convolveExposure(difference, scienceKernel, subtractTask.convolutionControl)
         output = detectionTask.run(science, matchedTemplate, difference, score, sources)
         crMaskSet = (output.subtractedMeasuredExposure.mask.array & crMask) > 0
-        self.assertTrue(np.all(crMaskSet[crX0:crX0+1, crY0:crY0+5]))
-        self.assertTrue(np.all(crMaskSet[crX1:crX1+5, crY1:crY1+1]))
+        self.assertFalse(np.any(crMaskSet[~crMaskSetInput]))
+        self.assertTrue(np.all(crMaskSet[crMaskSetInput]))
 
 
 class TestNegativePeaks(lsst.utils.tests.TestCase):
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 74a8457d..34a54b7d 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -27,6 +27,7 @@
 import lsst.utils.tests
 
 from test_detectAndMeasure import makeVisitInfo, MockResponse
+from utils import makeTestImage
 
 
 class UtilsTest(lsst.utils.tests.TestCase):
@@ -43,3 +44,7 @@ def test_populate_sattle_raises(self):
         with mock.patch('requests.put', return_value=response):
             with self.assertRaises(requests.exceptions.HTTPError):
                 populate_sattle_visit_cache(visit_info)
+
+    def test_raise_on_even_kernel(self):
+        with self.assertRaises(ValueError):
+            makeTestImage(kernelSize=32)

From 4cbfe53e143956e290d89c5f219d9b19761949f1 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 18:31:18 -0700
Subject: [PATCH 19/24] Clean up diffim unit tests

Consolidate setting common values, but there should be no behavior changes.
---
 tests/test_detectAndMeasure.py | 129 +++++++++++++++++-----------
 tests/test_subtractTask.py     | 152 ++++++++++++++++++---------------
 2 files changed, 162 insertions(+), 119 deletions(-)

diff --git a/tests/test_detectAndMeasure.py b/tests/test_detectAndMeasure.py
index 487091c2..42e2097c 100644
--- a/tests/test_detectAndMeasure.py
+++ b/tests/test_detectAndMeasure.py
@@ -168,7 +168,8 @@ def test_detection_xy0(self):
         noiseLevel = 1.
         staticSeed = 1
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         science.getInfo().setVisitInfo(makeVisitInfo())
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
@@ -204,7 +205,8 @@ def test_raise_bad_psf(self):
         noiseLevel = 1.
         staticSeed = 1
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
         difference = science.clone()
@@ -237,7 +239,8 @@ def test_measurements_finite(self):
         transientSeed = 6
         xSize = 256
         ySize = 256
-        kwargs = {"psfSize": 2.4, "x0": 0, "y0": 0,
+        psfSize = 2.4
+        kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
                   "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
                                          nSrc=1, **kwargs)
@@ -285,7 +288,8 @@ def test_remove_unphysical(self):
         staticSeed = 1
         xSize = 256
         ySize = 256
-        kwargs = {"psfSize": 2.4, "xSize": xSize, "ySize": ySize}
+        psfSize = 2.4
+        kwargs = {"psfSize": psfSize, "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
                                          nSrc=1, **kwargs)
         science.getInfo().setVisitInfo(makeVisitInfo())
@@ -331,19 +335,22 @@ def test_detect_transients(self):
         noiseLevel = 1.
         staticSeed = 1
         transientSeed = 6
+        nTransients = 10
+        transientFlux = 1000
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
         scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
 
         # Configure the detection Task
         detectionTask = self._setup_detection(doMerge=False, doSkySources=True)
         kwargs["seed"] = transientSeed
-        kwargs["nSrc"] = 10
-        kwargs["fluxLevel"] = 1000
+        kwargs["nSrc"] = nTransients
+        kwargs["fluxLevel"] = transientFlux
 
         # Run detection and check the results
-        def _detection_wrapper(positive=True):
+        def _run_and_check_detections(positive=True):
             transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
             science = scienceBase.clone()
             if positive:
@@ -361,8 +368,8 @@ def _detection_wrapper(positive=True):
             scale = 1. if positive else -1.
             for diaSource in output.diaSources:
                 self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
-        _detection_wrapper(positive=True)
-        _detection_wrapper(positive=False)
+        _run_and_check_detections(positive=True)
+        _run_and_check_detections(positive=False)
 
     def test_detect_transients_with_background(self):
         """Run detection on a difference image containing transients and a background.
@@ -371,12 +378,15 @@ def test_detect_transients_with_background(self):
         noiseLevel = 1.
         staticSeed = 1
         transientSeed = 6
+        nTransients = 10
+        transientFlux = 1000
         fluxLevel = 500
         xSize = 512
         ySize = 512
         x0 = 123
         y0 = 456
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
                   "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0}
         params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
 
@@ -389,11 +399,11 @@ def test_detect_transients_with_background(self):
         # Configure the detection Task
         detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
         kwargs["seed"] = transientSeed
-        kwargs["nSrc"] = 10
-        kwargs["fluxLevel"] = 1000
+        kwargs["nSrc"] = nTransients
+        kwargs["fluxLevel"] = transientFlux
 
         # Run detection and check the results
-        def _detection_wrapper(positive=True):
+        def _run_and_check_detections(positive=True):
             transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
             science = scienceBase.clone()
             if positive:
@@ -411,8 +421,8 @@ def _detection_wrapper(positive=True):
             scale = 1. if positive else -1.
             for transient in transientSources:
                 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
-        _detection_wrapper(positive=True)
-        _detection_wrapper(positive=False)
+        _run_and_check_detections(positive=True)
+        _run_and_check_detections(positive=False)
 
     def test_mask_cosmic_rays(self):
         """Run detection on a difference image containing a cosmic ray.
@@ -423,7 +433,9 @@ def test_mask_cosmic_rays(self):
         fluxLevel = 500
         xSize = 400
         ySize = 400
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "xSize": xSize, "ySize": ySize}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
+                  "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
         crMask = science.mask.getPlaneBitMask("CR")
@@ -462,7 +474,8 @@ def test_missing_mask_planes(self):
         # Set up the simulated images
         noiseLevel = 1.
         fluxLevel = 500
-        kwargs = {"psfSize": 2.4, "fluxLevel": fluxLevel, "addMaskPlanes": []}
+        psfSize = 2.4
+        kwargs = {"psfSize": psfSize, "fluxLevel": fluxLevel, "addMaskPlanes": []}
         # Use different seeds for the science and template so every source is a diaSource
         science, sources = makeTestImage(seed=5, noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         science.getInfo().setVisitInfo(makeVisitInfo())
@@ -488,11 +501,12 @@ def test_detect_dipoles(self):
         xSize = 300
         ySize = 300
         kernelSize = 31
+        psfSize = 2.4
         # Avoid placing sources near the edge for this test, so that we can
         # easily check that the correct number of sources are detected.
         templateBorderSize = kernelSize//2
         dipoleFlag = "ip_diffim_DipoleFit_classification"
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
                   "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
                   "xSize": xSize, "ySize": ySize}
         dipoleFlag = "ip_diffim_DipoleFit_classification"
@@ -535,7 +549,8 @@ def test_sky_sources(self):
         transientFluxLevel = 1000.
         transientFluxRange = 1.5
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         science.getInfo().setVisitInfo(makeVisitInfo())
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
@@ -580,7 +595,8 @@ def test_exclude_mask_detections(self):
         transientSeed = 6
         fluxLevel = 500
         radius = 2
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         science.getInfo().setVisitInfo(makeVisitInfo())
         detector = DetectorWrapper(numAmps=1).detector
@@ -601,7 +617,7 @@ def test_exclude_mask_detections(self):
         kwargs["fluxLevel"] = 1000
 
         # Run detection and check the results
-        def _detection_wrapper(setFlags=True):
+        def _run_and_check_detections(setFlags=True):
             transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
             difference = science.clone()
             difference.maskedImage -= matchedTemplate.maskedImage
@@ -626,8 +642,8 @@ def _detection_wrapper(setFlags=True):
                             self._check_diaSource(transientSources, diaSource, refIds=refIds)
                     else:
                         self._check_diaSource(transientSources, diaSource, refIds=refIds)
-        _detection_wrapper(setFlags=False)
-        _detection_wrapper(setFlags=True)
+        _run_and_check_detections(setFlags=False)
+        _run_and_check_detections(setFlags=True)
 
     def test_fake_mask_plane_propagation(self):
         """Test that we have the mask planes related to fakes in diffim images.
@@ -721,7 +737,9 @@ def test_mask_streaks(self):
         fluxLevel = 500
         xSize = 400
         ySize = 400
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "xSize": xSize, "ySize": ySize}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
+                  "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
 
@@ -779,7 +797,8 @@ def _setup_sattle_tests(self):
         noiseLevel = 1.
         staticSeed = 1
         fluxLevel = 500
-        shared_kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+        psfSize = 2.4
+        shared_kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
                          "x0": 12345, "y0": 67890}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
                                          **shared_kwargs)
@@ -906,11 +925,12 @@ def test_detection_xy0(self):
         staticSeed = 1
         fluxLevel = 500
         kernelSize = 31
+        psfSize = 2.4
         # Buffer source positions away from the score image's EDGE strip.
         # Preconvolution adds ~kernelSize//2 pixels of EDGE on top of the
         # EDGE already set by detection smoothing on the science image.
         templateBorderSize = kernelSize//2
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "x0": 12345, "y0": 67890,
                   "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
@@ -954,7 +974,8 @@ def test_measurements_finite(self):
         transientSeed = 6
         xSize = 256
         ySize = 256
-        kwargs = {"psfSize": 2.4, "x0": 0, "y0": 0,
+        psfSize = 2.4
+        kwargs = {"psfSize": psfSize, "x0": 0, "y0": 0,
                   "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(seed=staticSeed, noiseLevel=noiseLevel, noiseSeed=6,
                                          nSrc=1, **kwargs)
@@ -996,8 +1017,11 @@ def test_detect_transients(self):
         noiseLevel = 1.
         staticSeed = 1
         transientSeed = 6
+        nTransients = 10
+        transientFlux = 1000
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
         scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
         scienceKernel = scienceBase.psf.getKernel()
@@ -1006,11 +1030,11 @@ def test_detect_transients(self):
         # Configure the detection Task
         detectionTask = self._setup_detection(doMerge=False)
         kwargs["seed"] = transientSeed
-        kwargs["nSrc"] = 10
-        kwargs["fluxLevel"] = 1000
+        kwargs["nSrc"] = nTransients
+        kwargs["fluxLevel"] = transientFlux
 
         # Run detection and check the results
-        def _detection_wrapper(positive=True):
+        def _run_and_check_detections(positive=True):
             """Simulate positive or negative transients and run detection.
 
             Parameters
@@ -1040,8 +1064,8 @@ def _detection_wrapper(positive=True):
             for diaSource, goodSrcFlag in zip(output.diaSources, goodSrcFlags):
                 if goodSrcFlag:
                     self._check_diaSource(transientSources, diaSource, refIds=refIds, scale=scale)
-        _detection_wrapper(positive=True)
-        _detection_wrapper(positive=False)
+        _run_and_check_detections(positive=True)
+        _run_and_check_detections(positive=False)
 
     def test_detect_dipoles(self):
         """Run detection on a difference image containing dipoles.
@@ -1056,11 +1080,12 @@ def test_detect_dipoles(self):
         xSize = 300
         ySize = 300
         kernelSize = 31
+        psfSize = 2.4
         # Avoid placing sources near the edge for this test, so that we can
         # easily check that the correct number of sources are detected.
         templateBorderSize = kernelSize//2
         dipoleFlag = "ip_diffim_DipoleFit_classification"
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel, "fluxRange": fluxRange,
                   "nSrc": nSources, "templateBorderSize": templateBorderSize, "kernelSize": kernelSize,
                   "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
@@ -1101,7 +1126,8 @@ def test_sky_sources(self):
         transientFluxLevel = 1000.
         transientFluxRange = 1.5
         fluxLevel = 500
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
         transients, transientSources = makeTestImage(seed=transientSeed, psfSize=2.4,
@@ -1149,7 +1175,8 @@ def test_exclude_mask_detections(self):
         transientSeed = 6
         fluxLevel = 500
         radius = 2
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel}
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         science.getInfo().setVisitInfo(makeVisitInfo())
         detector = DetectorWrapper(numAmps=1).detector
@@ -1172,7 +1199,7 @@ def test_exclude_mask_detections(self):
         kwargs["fluxLevel"] = 1000
 
         # Run detection and check the results
-        def _detection_wrapper(setFlags=True):
+        def _run_and_check_detections(setFlags=True):
             transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
             difference = science.clone()
             difference.maskedImage -= matchedTemplate.maskedImage
@@ -1199,8 +1226,8 @@ def _detection_wrapper(setFlags=True):
                             self._check_diaSource(transientSources, diaSource, refIds=refIds)
                     else:
                         self._check_diaSource(transientSources, diaSource, refIds=refIds)
-        _detection_wrapper(setFlags=False)
-        _detection_wrapper(setFlags=True)
+        _run_and_check_detections(setFlags=False)
+        _run_and_check_detections(setFlags=True)
 
     def test_detect_transients_with_background(self):
         """Background measured on the difference image should be subtracted
@@ -1210,24 +1237,27 @@ def test_detect_transients_with_background(self):
         noiseLevel = 1.
         staticSeed = 1
         transientSeed = 6
+        nTransients = 10
+        transientFlux = 1000
         fluxLevel = 500
         xSize = 512
         ySize = 512
         x0 = 123
         y0 = 456
         kernelSize = 31
+        psfSize = 2.4
         templateBorderSize = kernelSize//2
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
                   "xSize": xSize, "ySize": ySize, "x0": x0, "y0": y0,
                   "kernelSize": kernelSize, "templateBorderSize": templateBorderSize}
-        params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
+        chebyshev_params = [2.2, 2.1, 2.0, 1.2, 1.1, 1.0]
 
         # The background model must cover the grown image bbox, otherwise it
         # would be extrapolated outside its declared domain.
         bbox2D = lsst.geom.Box2D(
             lsst.geom.Point2D(x0 - templateBorderSize, y0 - templateBorderSize),
             lsst.geom.Extent2D(xSize + 2*templateBorderSize, ySize + 2*templateBorderSize))
-        background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
+        background_model = afwMath.Chebyshev1Function2D(chebyshev_params, bbox2D)
         scienceBase, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6,
                                              background=background_model, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
@@ -1237,10 +1267,10 @@ def test_detect_transients_with_background(self):
         # Configure the detection Task
         detectionTask = self._setup_detection(doMerge=False, doSubtractBackground=True)
         kwargs["seed"] = transientSeed
-        kwargs["nSrc"] = 10
-        kwargs["fluxLevel"] = 1000
+        kwargs["nSrc"] = nTransients
+        kwargs["fluxLevel"] = transientFlux
 
-        def _detection_wrapper(positive=True):
+        def _run_and_check_detections(positive=True):
             transients, transientSources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=8, **kwargs)
             science = scienceBase.clone()
             if positive:
@@ -1270,8 +1300,8 @@ def _detection_wrapper(positive=True):
             scale = 1. if positive else -1.
             for transient in transientSources:
                 self._check_diaSource(output.diaSources, transient, refIds=refIds, scale=scale)
-        _detection_wrapper(positive=True)
-        _detection_wrapper(positive=False)
+        _run_and_check_detections(positive=True)
+        _run_and_check_detections(positive=False)
 
     def test_mask_cosmic_rays(self):
         """Cosmic rays detected on the difference image should propagate
@@ -1283,7 +1313,8 @@ def test_mask_cosmic_rays(self):
         fluxLevel = 500
         xSize = 400
         ySize = 400
-        kwargs = {"seed": staticSeed, "psfSize": 2.4, "fluxLevel": fluxLevel,
+        psfSize = 2.4
+        kwargs = {"seed": staticSeed, "psfSize": psfSize, "fluxLevel": fluxLevel,
                   "xSize": xSize, "ySize": ySize}
         science, sources = makeTestImage(noiseLevel=noiseLevel, noiseSeed=6, **kwargs)
         matchedTemplate, _ = makeTestImage(noiseLevel=noiseLevel/4, noiseSeed=7, **kwargs)
diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py
index 728ae362..954d0273 100644
--- a/tests/test_subtractTask.py
+++ b/tests/test_subtractTask.py
@@ -40,6 +40,10 @@
 
 
 class AlardLuptonSubtractTestBase:
+    goodPsfSize = 2.0
+    midPsfSize = 2.4
+    badPsfSize = 2.8
+
     def _setup_subtraction(self, fluxField="truth_instFlux", errField="truth_instFluxErr", **kwargs):
         """Setup and configure the image subtraction PipelineTask.
 
@@ -93,8 +97,9 @@ def test_mismatched_template(self):
         """
         xSize = 200
         ySize = 200
-        science, sources = makeTestImage(psfSize=2.4, xSize=xSize + 20, ySize=ySize + 20)
-        template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
+        science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
+                                    doApplyCalibration=True)
         task = self._setup_subtraction()
         with self.assertRaises(AssertionError):
             task.run(template, science, sources)
@@ -105,10 +110,10 @@ def test_mismatched_filter(self):
         """
         xSize = 200
         ySize = 200
-        science, sources = makeTestImage(psfSize=2.4, xSize=xSize + 20, ySize=ySize + 20,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20,
                                          band="g", physicalFilter="g noCamera")
-        template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True,
-                                    band="not-g", physicalFilter="not-g noCamera")
+        template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
+                                    doApplyCalibration=True, band="not-g", physicalFilter="not-g noCamera")
         task = self._setup_subtraction()
         with self.assertRaises(AssertionError):
             task.run(template, science, sources)
@@ -119,10 +124,10 @@ def test_incomplete_template_coverage(self):
         xSize = 400
         ySize = 400
         nSources = 80
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6, nSrc=nSources,
-                                         xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7, nSrc=nSources,
-                                    templateBorderSize=border, doApplyCalibration=True,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
+                                         nSrc=nSources, xSize=xSize, ySize=ySize)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
+                                    nSrc=nSources, templateBorderSize=border, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
 
         science_height = science.getBBox().getDimensions().getY()
@@ -155,8 +160,8 @@ def _run_and_check_coverage(template_coverage,
 
     def test_clear_template_mask(self):
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
         task = self._setup_subtraction(mode="convolveTemplate")
@@ -196,8 +201,8 @@ def test_equal_images(self):
         with the same size psf in the template and science.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction()
         output = task.run(template, science, sources)
@@ -220,8 +225,9 @@ def test_equal_images_missing_mask_planes(self):
         mask planes.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6, addMaskPlanes=[])
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
+                                         addMaskPlanes=[])
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True, addMaskPlanes=[])
         task = self._setup_subtraction()
         output = task.run(template, science, sources)
@@ -243,8 +249,8 @@ def test_psf_size(self):
         fwhmExposureBuffer and fwhmExposureGrid parameters are set.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
 
         schema = afwTable.ExposureTable.makeMinimalSchema()
@@ -304,8 +310,8 @@ def test_auto_convolveTemplate(self):
         the template psf is the smaller.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction(mode="convolveTemplate")
         output = task.run(template.clone(), science.clone(), sources)
@@ -319,8 +325,8 @@ def test_auto_convolveScience(self):
         the science psf is the smaller.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction(mode="convolveScience")
         output = task.run(template.clone(), science.clone(), sources)
@@ -337,8 +343,9 @@ def test_science_better(self):
         statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
 
         def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
-            science, sources = makeTestImage(psfSize=2.0, noiseLevel=scienceNoiseLevel, noiseSeed=6)
-            template, _ = makeTestImage(psfSize=2.3, noiseLevel=templateNoiseLevel, noiseSeed=7,
+            science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel,
+                                             noiseSeed=6)
+            template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                         templateBorderSize=20, doApplyCalibration=True)
             task = self._setup_subtraction(mode="convolveScience")
             output = task.run(template, science, sources)
@@ -372,8 +379,9 @@ def test_template_better(self):
         statsCtrlDetect = makeStats(badMaskPlanes=("EDGE", "BAD", "NO_DATA"))
 
         def _run_and_check_images(statsCtrl, statsCtrlDetect, scienceNoiseLevel, templateNoiseLevel):
-            science, sources = makeTestImage(psfSize=2.4, noiseLevel=scienceNoiseLevel, noiseSeed=6)
-            template, _ = makeTestImage(psfSize=2.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
+            science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel,
+                                             noiseSeed=6)
+            template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                         templateBorderSize=20, doApplyCalibration=True)
             task = self._setup_subtraction()
             output = task.run(template, science, sources)
@@ -408,9 +416,9 @@ def test_symmetry(self):
         noiseLevel = 1.
         # Don't include a border for the template, in order to make the results
         #  comparable when we swap which image is treated as the "science" image.
-        science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel,
+        science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel,
                                          noiseSeed=6, templateBorderSize=0, doApplyCalibration=True)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel,
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel,
                                     noiseSeed=7, templateBorderSize=0, doApplyCalibration=True)
         task = self._setup_subtraction(mode='auto')
 
@@ -438,8 +446,9 @@ def test_few_sources(self):
         """
         xSize = 256
         ySize = 256
-        science, sources = makeTestImage(psfSize=2.4, nSrc=10, xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, nSrc=10, xSize=xSize, ySize=ySize, doApplyCalibration=True)
+        science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize,
+                                    doApplyCalibration=True)
         task = self._setup_subtraction()
         sources = sources[0:1]
         with self.assertRaises(InsufficientKernelSourcesError):
@@ -451,8 +460,8 @@ def test_kernel_source_selector(self):
         xSize = 256
         ySize = 256
         nSourcesSimulated = 20
-        sciencePsfSize = 2.4
-        templatePsfSize = 2.0
+        sciencePsfSize = self.midPsfSize
+        templatePsfSize = self.goodPsfSize
         science, sources = makeTestImage(psfSize=sciencePsfSize, nSrc=nSourcesSimulated,
                                          xSize=xSize, ySize=ySize)
         template, _ = makeTestImage(psfSize=templatePsfSize, nSrc=nSourcesSimulated,
@@ -553,7 +562,7 @@ def test_background_subtraction(self):
         ySize = 512
         x0 = 123
         y0 = 456
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20,
                                     xSize=xSize, ySize=ySize, x0=x0, y0=y0,
                                     doApplyCalibration=True)
@@ -561,7 +570,7 @@ def test_background_subtraction(self):
 
         bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
         background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
-        science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
                                          background=background_model,
                                          xSize=xSize, ySize=ySize, x0=x0, y0=y0)
         # Don't use ``self._setup_subtraction()`` here.
@@ -640,8 +649,8 @@ def _run_and_check_images(science, template, sources, statsCtrl,
             varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
             self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
 
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=scienceNoiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         # Verify that the variance plane of the difference image is correct
         #  when the template and science variance planes are correct
@@ -709,8 +718,8 @@ def _run_and_check_images(science, template, sources, statsCtrl,
             varMean = computeRobustStatistics(output.difference.variance, output.difference.mask, statsCtrl)
             self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
 
-        science, sources = makeTestImage(psfSize=2.0, noiseLevel=scienceNoiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=templateNoiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         # Verify that the variance plane of the difference image is correct
         #  when the template and science variance planes are correct
@@ -743,12 +752,12 @@ def test_exposure_properties_convolve_template(self):
         noiseLevel = 1.
         seed = 37
         rng = np.random.RandomState(seed)
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
         psf = science.psf
         psfAvgPos = psf.getAveragePosition()
         psfSize = getPsfFwhm(science.psf)
         psfImg = psf.computeKernelImage(psfAvgPos)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
 
         # Generate a random aperture correction map
@@ -790,8 +799,8 @@ def test_exposure_properties_convolve_science(self):
         noiseLevel = 1.
         seed = 37
         rng = np.random.RandomState(seed)
-        science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         psf = template.psf
         psfAvgPos = psf.getAveragePosition()
@@ -856,13 +865,14 @@ def test_fake_mask_plane_propagation(self):
         """
         xSize = 200
         ySize = 200
-        science, sources = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize)
+        science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize)
         science_fake_img, science_fake_sources = makeTestImage(
-            psfSize=2.4, xSize=xSize, ySize=ySize, seed=7, nSrc=2, noiseLevel=0.25, fluxRange=1
+            psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=7, nSrc=2, noiseLevel=0.25, fluxRange=1
         )
-        template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
+                                    doApplyCalibration=True)
         tmplt_fake_img, tmplt_fake_sources = makeTestImage(
-            psfSize=2.4, xSize=xSize, ySize=ySize, seed=9, nSrc=2, noiseLevel=0.25, fluxRange=1
+            psfSize=self.midPsfSize, xSize=xSize, ySize=ySize, seed=9, nSrc=2, noiseLevel=0.25, fluxRange=1
         )
         # created fakes and added them to the images
         science.image.array += science_fake_img.image.array
@@ -911,8 +921,8 @@ def test_metadata_metrics(self):
         that the difference image limiting magnitude is calculated correctly,
         both with a "good" and "bad" seeing template.
         """
-        science, sources = makeTestImage(psfSize=2.8, noiseLevel=1)
-        template_good, _ = makeTestImage(psfSize=2.4, doApplyCalibration=True, noiseLevel=0.25,
+        science, sources = makeTestImage(psfSize=self.badPsfSize, noiseLevel=1)
+        template_good, _ = makeTestImage(psfSize=self.midPsfSize, doApplyCalibration=True, noiseLevel=0.25,
                                          templateBorderSize=20)
         template_bad, _ = makeTestImage(psfSize=9.5, doApplyCalibration=True, noiseLevel=0.25,
                                         templateBorderSize=20)
@@ -1003,8 +1013,9 @@ def test_mismatched_template(self):
         """
         xSize = 200
         ySize = 200
-        science, sources = makeTestImage(psfSize=2.4, xSize=xSize + 20, ySize=ySize + 20)
-        template, _ = makeTestImage(psfSize=2.4, xSize=xSize, ySize=ySize, doApplyCalibration=True)
+        science, sources = makeTestImage(psfSize=self.midPsfSize, xSize=xSize + 20, ySize=ySize + 20)
+        template, _ = makeTestImage(psfSize=self.midPsfSize, xSize=xSize, ySize=ySize,
+                                    doApplyCalibration=True)
         task = self._setup_subtraction()
         with self.assertRaises(AssertionError):
             task.run(template, science, sources)
@@ -1016,9 +1027,9 @@ def test_equal_images(self):
         noiseLevel = 1.
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
         task = self._setup_subtraction()
@@ -1044,9 +1055,9 @@ def test_incomplete_template_coverage(self):
         border = 20
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=border, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
 
@@ -1082,9 +1093,9 @@ def test_clear_template_mask(self):
         noiseLevel = 1.
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
         diffimEmptyMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]
@@ -1127,13 +1138,13 @@ def test_agnostic_template_psf(self):
         noiseLevel = .3
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel,
                                          noiseSeed=6, templateBorderSize=0,
                                          xSize=xSize, ySize=ySize)
-        template1, _ = makeTestImage(psfSize=2.8, noiseLevel=noiseLevel,
+        template1, _ = makeTestImage(psfSize=self.badPsfSize, noiseLevel=noiseLevel,
                                      noiseSeed=7, doApplyCalibration=True,
                                      xSize=xSize, ySize=ySize)
-        template2, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel,
+        template2, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel,
                                      noiseSeed=8, doApplyCalibration=True,
                                      xSize=xSize, ySize=ySize)
         task = self._setup_subtraction()
@@ -1164,8 +1175,9 @@ def test_few_sources(self):
         """
         xSize = 256
         ySize = 256
-        science, sources = makeTestImage(psfSize=2.4, nSrc=10, xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, nSrc=10, xSize=xSize, ySize=ySize, doApplyCalibration=True)
+        science, sources = makeTestImage(psfSize=self.midPsfSize, nSrc=10, xSize=xSize, ySize=ySize)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, nSrc=10, xSize=xSize, ySize=ySize,
+                                    doApplyCalibration=True)
         task = self._setup_subtraction()
         sources = sources[0:1]
         with self.assertRaises(InsufficientKernelSourcesError):
@@ -1180,7 +1192,7 @@ def test_background_subtraction(self):
         ySize = 512
         x0 = 123
         y0 = 456
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20,
                                     xSize=xSize, ySize=ySize, x0=x0, y0=y0,
                                     doApplyCalibration=True)
@@ -1188,7 +1200,7 @@ def test_background_subtraction(self):
 
         bbox2D = lsst.geom.Box2D(lsst.geom.Point2D(x0, y0), lsst.geom.Extent2D(xSize, ySize))
         background_model = afwMath.Chebyshev1Function2D(params, bbox2D)
-        science, sources = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
                                          background=background_model,
                                          xSize=xSize, ySize=ySize, x0=x0, y0=y0)
         # Don't use ``self._setup_subtraction()`` here.
@@ -1271,9 +1283,9 @@ def _run_and_check_images(science, template, sources, statsCtrl,
                                               statsCtrl)
             self.assertFloatsAlmostEqual(varMean, scienceNoise + templateNoise, rtol=0.1)
 
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=scienceNoiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=scienceNoiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=templateNoiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=templateNoiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
         # Verify that the variance plane of the Score image is correct
@@ -1307,13 +1319,13 @@ def test_exposure_properties(self):
         noiseLevel = 1.
         xSize = 400
         ySize = 400
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6,
                                          xSize=xSize, ySize=ySize)
         psf = science.psf
         psfAvgPos = psf.getAveragePosition()
         psfSize = getPsfFwhm(science.psf)
         psfImg = psf.computeKernelImage(psfAvgPos)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True,
                                     xSize=xSize, ySize=ySize)
 
@@ -1349,8 +1361,8 @@ def test_runSimplifiedTaskWithExistingKernel(self):
         `AlardLuptonSubtractTask` if it uses the AL kernel.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         alTask = AlardLuptonSubtractTest._setup_subtraction(AlardLuptonSubtractTest())
         task = self._setup_subtraction(useExistingKernel=True)
@@ -1366,8 +1378,8 @@ def test_runSimplifiedTaskWithSourceDetection(self):
         reasonable output.
         """
         noiseLevel = 1.
-        science, sources = makeTestImage(psfSize=2.4, noiseLevel=noiseLevel, noiseSeed=6)
-        template, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=7,
+        science, sources = makeTestImage(psfSize=self.midPsfSize, noiseLevel=noiseLevel, noiseSeed=6)
+        template, _ = makeTestImage(psfSize=self.goodPsfSize, noiseLevel=noiseLevel, noiseSeed=7,
                                     templateBorderSize=20, doApplyCalibration=True)
         task = self._setup_subtraction(useExistingKernel=False,
                                        fluxField="base_PsfFlux_instFlux",

From b1ec8e8fdb01cc979326f6400114f3d287268787 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 18:35:13 -0700
Subject: [PATCH 20/24] Pass keywords by name to `kernel.computeImage`

---
 python/lsst/ip/diffim/psfMatch.py       | 11 ++++-------
 python/lsst/ip/diffim/subtractImages.py |  4 ++--
 2 files changed, 6 insertions(+), 9 deletions(-)

diff --git a/python/lsst/ip/diffim/psfMatch.py b/python/lsst/ip/diffim/psfMatch.py
index 7c46f9e4..fb4677df 100644
--- a/python/lsst/ip/diffim/psfMatch.py
+++ b/python/lsst/ip/diffim/psfMatch.py
@@ -595,12 +595,9 @@ def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
         """
         # What is the final kernel sum
         kImage = afwImage.ImageD(spatialKernel.getDimensions())
-        # Evaluate the kernel at the cell-set center, not at the world
-        # origin (the default (0, 0) lies outside the cellset bbox for any
-        # non-origin region and produces a meaningless kernel sum).
-        regionBBox = kernelCellSet.getBBox()
-        xcen, ycen = regionBBox.getCenter()
-        kSum = spatialKernel.computeImage(kImage, False, xcen, ycen)
+        # Evaluate the kernel at the cell-set center, not at the world origin.
+        xcen, ycen = kernelCellSet.getBBox().getCenter()
+        kSum = spatialKernel.computeImage(kImage, doNormalize=False, x=xcen, y=ycen)
         self.log.info("Final spatial kernel sum %.3f", kSum)
 
         # Look at how well conditioned the matrix is
@@ -778,7 +775,7 @@ def _createPcaBasis(self, kernelCellSet, nStarPerCell, ps):
         for j in range(nToUse):
             # Check for NaNs?
             kimage = afwImage.ImageD(pcaBasisList[j].getDimensions())
-            pcaBasisList[j].computeImage(kimage, False)
+            pcaBasisList[j].computeImage(kimage, doNormalize=False)
             if not (True in np.isnan(kimage.array)):
                 trimBasisList.append(pcaBasisList[j])
 
diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 18b15b76..92784b59 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -709,7 +709,7 @@ def runConvolveScience(self, template, science, psfMatchingKernel, backgroundMod
 
         kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions())
         xcen, ycen = bbox.getCenter()
-        norm = psfMatchingKernel.computeImage(kernelImage, False, xcen, ycen)
+        norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False, x=xcen, y=ycen)
 
         matchedScience = self._convolveExposure(science, psfMatchingKernel,
                                                 self.convolutionControl,
@@ -1321,7 +1321,7 @@ def _makePreconvolutionKernel(psf):
                 f"Preconvolution requires an odd-sized PSF kernel, got {dims.x}x{dims.y}. "
             )
         kimg = lsst.afw.image.ImageD(dims)
-        localKernel.computeImage(kimg, True)  # normalize to unit sum
+        localKernel.computeImage(kimg, doNormalize=True)  # normalize to unit sum
         # Reflect about the kernel center. PSF kernels are odd-sized,
         # so ``[::-1, ::-1]`` places the peak at the same pixel.
         kimg.array[...] = kimg.array[::-1, ::-1]

From 48c9c1d7d017c56d5861faeedb983b21fa79b7a4 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 18:42:23 -0700
Subject: [PATCH 21/24] Move duplicated PSF width checks in makeKernel to
 common method

Also remove the redundant calculation of the kernel source list in subtractImages - makeCandidateList was effectively called twice.
---
 python/lsst/ip/diffim/makeKernel.py     | 112 +++++++++++++-----------
 python/lsst/ip/diffim/subtractImages.py |  15 +---
 2 files changed, 65 insertions(+), 62 deletions(-)

diff --git a/python/lsst/ip/diffim/makeKernel.py b/python/lsst/ip/diffim/makeKernel.py
index ad0a1d2d..e647ca23 100644
--- a/python/lsst/ip/diffim/makeKernel.py
+++ b/python/lsst/ip/diffim/makeKernel.py
@@ -122,6 +122,9 @@ def run(self, template, science, kernelSources, preconvolved=False,
             Typically the output of `MakeKernelTask.selectKernelSources`.
         preconvolved : `bool`, optional
             Was the science image convolved with its own PSF?
+        templateFwhmPix, scienceFwhmPix : `float`, optional
+            FWHM of the template or science PSF, in pixels.
+            Will be recalculated if not specified.
 
         Returns
         -------
@@ -133,33 +136,9 @@ def run(self, template, science, kernelSources, preconvolved=False,
                 Spatially varying background-matching function.
         """
         kernelCellSet = self._buildCellSet(template.maskedImage, science.maskedImage, kernelSources)
-        if (scienceFwhmPix is None) or (templateFwhmPix is None):
-            # Calling getPsfFwhm on template.psf fails on some rare occasions when
-            # the template has no input exposures at the average position of the
-            # stars. So we try getPsfFwhm first on template, and if that fails we
-            # evaluate the PSF on a grid specified by fwhmExposure* fields.
-            # To keep consistent definitions for PSF size on the template and
-            # science images, we use the same method for both.
-            try:
-                templateFwhmPix = getPsfFwhm(template.psf)
-                scienceFwhmPix = getPsfFwhm(science.psf)
-            except (InvalidParameterError, RangeError):
-                self.log.debug("Unable to evaluate PSF at the average position. "
-                               "Evaluting PSF on a grid of points."
-                               )
-                templateFwhmPix = evaluateMeanPsfFwhm(template,
-                                                      fwhmExposureBuffer=self.config.fwhmExposureBuffer,
-                                                      fwhmExposureGrid=self.config.fwhmExposureGrid
-                                                      )
-                scienceFwhmPix = evaluateMeanPsfFwhm(science,
-                                                     fwhmExposureBuffer=self.config.fwhmExposureBuffer,
-                                                     fwhmExposureGrid=self.config.fwhmExposureGrid
-                                                     )
-            # Apply the Gaussian approximation FWHM boost only when
-            # the size was auto-computed. An explicit ``scienceFwhmPix`` is
-            # assumed already to reflect the preconvolved effective width.
-            if preconvolved:
-                scienceFwhmPix *= np.sqrt(2)
+        templateFwhmPix, scienceFwhmPix = self._checkPsfWidths(template, science,
+                                                               templateFwhmPix, scienceFwhmPix,
+                                                               preconvolved=False)
         basisList = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix,
                                              metadata=self.metadata)
         spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
@@ -182,36 +161,18 @@ def selectKernelSources(self, template, science, candidateList=None, preconvolve
             Sources to check as possible kernel candidates.
         preconvolved : `bool`, optional
             Was the science image convolved with its own PSF?
+        templateFwhmPix, scienceFwhmPix : `float`, optional
+            FWHM of the template or science PSF, in pixels.
+            Will be recalculated if not specified.
 
         Returns
         -------
         kernelSources : `lsst.afw.table.SourceCatalog`
             Kernel candidates with appropriate sized footprints.
         """
-        if (scienceFwhmPix is None) or (templateFwhmPix is None):
-            # Calling getPsfFwhm on template.psf fails on some rare occasions when
-            # the template has no input exposures at the average position of the
-            # stars. So we try getPsfFwhm first on template, and if that fails we
-            # evaluate the PSF on a grid specified by fwhmExposure* fields.
-            # To keep consistent definitions for PSF size on the template and
-            # science images, we use the same method for both.
-            try:
-                templateFwhmPix = getPsfFwhm(template.psf)
-                scienceFwhmPix = getPsfFwhm(science.psf)
-            except (InvalidParameterError, RangeError):
-                self.log.debug("Unable to evaluate PSF at the average position. "
-                               "Evaluting PSF on a grid of points."
-                               )
-                templateFwhmPix = evaluateMeanPsfFwhm(template,
-                                                      fwhmExposureBuffer=self.config.fwhmExposureBuffer,
-                                                      fwhmExposureGrid=self.config.fwhmExposureGrid
-                                                      )
-                scienceFwhmPix = evaluateMeanPsfFwhm(science,
-                                                     fwhmExposureBuffer=self.config.fwhmExposureBuffer,
-                                                     fwhmExposureGrid=self.config.fwhmExposureGrid
-                                                     )
-            if preconvolved:
-                scienceFwhmPix *= np.sqrt(2)
+        templateFwhmPix, scienceFwhmPix = self._checkPsfWidths(template, science,
+                                                               templateFwhmPix, scienceFwhmPix,
+                                                               preconvolved=False)
         kernelSize = self.makeKernelBasisList(templateFwhmPix, scienceFwhmPix)[0].getWidth()
         kernelSources = self.makeCandidateList(template, science, kernelSize,
                                                candidateList=candidateList,
@@ -454,3 +415,52 @@ def _adaptCellSize(self, candidateList):
             New dimensions to use for the kernel.
         """
         return self.kConfig.sizeCellX, self.kConfig.sizeCellY
+
+    def _checkPsfWidths(self, template, science, templateFwhmPix, scienceFwhmPix, preconvolved=False):
+        """Check the science and template FWHM.
+
+        Parameters
+        ----------
+        template : `lsst.afw.image.Exposure`
+            Exposure that will be convolved.
+        science : `lsst.afw.image.Exposure`
+            The exposure that will be matched.
+        templateFwhmPix, scienceFwhmPix : `float`, optional
+            FWHM of the template or science PSF, in pixels.
+            Will be recalculated if not specified.
+        preconvolved : `bool`, optional
+            Was the science image convolved with its own PSF?
+
+        Returns
+        -------
+        templateFwhmPix, scienceFwhmPix : `float`
+            The desired PSF FWHM
+        """
+        if (scienceFwhmPix is None) or (templateFwhmPix is None):
+            # Calling getPsfFwhm on template.psf fails on some rare occasions when
+            # the template has no input exposures at the average position of the
+            # stars. So we try getPsfFwhm first on template, and if that fails we
+            # evaluate the PSF on a grid specified by fwhmExposure* fields.
+            # To keep consistent definitions for PSF size on the template and
+            # science images, we use the same method for both.
+            try:
+                templateFwhmPix = getPsfFwhm(template.psf)
+                scienceFwhmPix = getPsfFwhm(science.psf)
+            except (InvalidParameterError, RangeError):
+                self.log.debug("Unable to evaluate PSF at the average position. "
+                               "Evaluting PSF on a grid of points."
+                               )
+                templateFwhmPix = evaluateMeanPsfFwhm(template,
+                                                      fwhmExposureBuffer=self.config.fwhmExposureBuffer,
+                                                      fwhmExposureGrid=self.config.fwhmExposureGrid
+                                                      )
+                scienceFwhmPix = evaluateMeanPsfFwhm(science,
+                                                     fwhmExposureBuffer=self.config.fwhmExposureBuffer,
+                                                     fwhmExposureGrid=self.config.fwhmExposureGrid
+                                                     )
+            # Apply the Gaussian approximation FWHM boost only when
+            # the size was auto-computed. An explicit ``scienceFwhmPix`` is
+            # assumed already to reflect the preconvolved effective width.
+            if preconvolved:
+                scienceFwhmPix *= np.sqrt(2)
+        return templateFwhmPix, scienceFwhmPix
diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 92784b59..56920a3e 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -613,17 +613,10 @@ def runKernelSourceDetection(self, template, science):
                 PSF-matching kernel.
         """
         kernelSize = self.makeKernel.makeKernelBasisList(
-            self.templatePsfSize, self.sciencePsfSize)[0].getWidth()
-        candidateList = self.makeKernel.makeCandidateList(template, science, kernelSize,
-                                                          candidateList=None,
-                                                          sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize)
-        sources = self.makeKernel.selectKernelSources(template, science,
-                                                      candidateList=candidateList,
-                                                      preconvolved=False,
-                                                      templateFwhmPix=self.templatePsfSize,
-                                                      scienceFwhmPix=self.sciencePsfSize)
-
-        # return sources
+            self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
+        sources = self.makeKernel.makeCandidateList(template, science, kernelSize,
+                                                    candidateList=None,
+                                                    sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize)
         return self._sourceSelector(template, science, sources, fallback=True)
 
     def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None):

From 964707e2d7aeb806dd047ff239c6136255add3da Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 18:48:31 -0700
Subject: [PATCH 22/24] Test a range of psf sizes in equal images unit test

---
 tests/test_subtractTask.py | 64 ++++++++++++++++++++------------------
 1 file changed, 33 insertions(+), 31 deletions(-)

diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py
index 954d0273..fa30f736 100644
--- a/tests/test_subtractTask.py
+++ b/tests/test_subtractTask.py
@@ -475,6 +475,7 @@ def _run_and_check_sources(sourcesIn, maxKernelSources=1000, minKernelSources=3)
                                            )
             task.templatePsfSize = templatePsfSize
             task.sciencePsfSize = sciencePsfSize
+            task.matchedPsfSize = sciencePsfSize
             # Verify that source flags are not set in the input catalog
             # Note that this will use the last flag in the list for the rest of
             #  the test.
@@ -516,37 +517,38 @@ def test_order_equal_images(self):
         noiseLevel = .1
         seed1 = 6
         seed2 = 7
-        science1, sources1 = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=seed1,
-                                           clearEdgeMask=True)
-        template1, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=seed2,
-                                     templateBorderSize=0, doApplyCalibration=True,
-                                     clearEdgeMask=True)
-        task1 = self._setup_subtraction(mode="convolveTemplate")
-        results_convolveTemplate = task1.run(template1, science1, sources1)
-
-        science2, sources2 = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=seed1,
-                                           clearEdgeMask=True)
-        template2, _ = makeTestImage(psfSize=2.0, noiseLevel=noiseLevel, noiseSeed=seed2,
-                                     templateBorderSize=0, doApplyCalibration=True,
-                                     clearEdgeMask=True)
-        task2 = self._setup_subtraction(mode="convolveScience")
-        results_convolveScience = task2.run(template2, science2, sources2)
-        bbox = results_convolveTemplate.difference.getBBox().clippedTo(
-            results_convolveScience.difference.getBBox())
-        diff1 = science1.maskedImage.clone()[bbox]
-        diff1 -= template1.maskedImage[bbox]
-        diff2 = science2.maskedImage.clone()[bbox]
-        diff2 -= template2.maskedImage[bbox]
-        self.assertFloatsAlmostEqual(results_convolveTemplate.difference[bbox].image.array,
-                                     diff1.image.array,
-                                     atol=noiseLevel*5.)
-        self.assertFloatsAlmostEqual(results_convolveScience.difference[bbox].image.array,
-                                     diff2.image.array,
-                                     atol=noiseLevel*5.)
-        diffErr = noiseLevel*2
-        self.assertMaskedImagesAlmostEqual(results_convolveTemplate.difference[bbox].maskedImage,
-                                           results_convolveScience.difference[bbox].maskedImage,
-                                           atol=diffErr*5.)
+        for psfSize in [self.midPsfSize, self.goodPsfSize, self.badPsfSize]:
+            science1, sources1 = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed1,
+                                               clearEdgeMask=True)
+            template1, _ = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed2,
+                                         templateBorderSize=0, doApplyCalibration=True,
+                                         clearEdgeMask=True)
+            task1 = self._setup_subtraction(mode="convolveTemplate")
+            results_convolveTemplate = task1.run(template1, science1, sources1)
+
+            science2, sources2 = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed1,
+                                               clearEdgeMask=True)
+            template2, _ = makeTestImage(psfSize=psfSize, noiseLevel=noiseLevel, noiseSeed=seed2,
+                                         templateBorderSize=0, doApplyCalibration=True,
+                                         clearEdgeMask=True)
+            task2 = self._setup_subtraction(mode="convolveScience")
+            results_convolveScience = task2.run(template2, science2, sources2)
+            bbox = results_convolveTemplate.difference.getBBox().clippedTo(
+                results_convolveScience.difference.getBBox())
+            diff1 = science1.maskedImage.clone()[bbox]
+            diff1 -= template1.maskedImage[bbox]
+            diff2 = science2.maskedImage.clone()[bbox]
+            diff2 -= template2.maskedImage[bbox]
+            self.assertFloatsAlmostEqual(results_convolveTemplate.difference[bbox].image.array,
+                                         diff1.image.array,
+                                         atol=noiseLevel*5.)
+            self.assertFloatsAlmostEqual(results_convolveScience.difference[bbox].image.array,
+                                         diff2.image.array,
+                                         atol=noiseLevel*5.)
+            diffErr = noiseLevel*2
+            self.assertMaskedImagesAlmostEqual(results_convolveTemplate.difference[bbox].maskedImage,
+                                               results_convolveScience.difference[bbox].maskedImage,
+                                               atol=diffErr*5.)
 
     def test_background_subtraction(self):
         """Check that we can recover the background,

From 3157ca6a4c7c83a2e79a574e16da49fb20a7b742 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 19:01:02 -0700
Subject: [PATCH 23/24] Refactor psfSize calculations to make sure it is done
 only once

---
 python/lsst/ip/diffim/subtractImages.py | 39 +++++++++++++------------
 1 file changed, 21 insertions(+), 18 deletions(-)

diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 56920a3e..9263e5c8 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -451,9 +451,11 @@ def run(self, template, science, sources, visitSummary=None):
                 Sources from the input catalog that were used to construct the
                 PSF-matching kernel.
         """
+        self.usePreconvolution = False
         self._prepareInputs(template, science, visitSummary=visitSummary)
 
         convolveTemplate = self.chooseConvolutionMethod(template, science)
+        self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
 
         kernelResult = self.runMakeKernel(template, science, sources=sources,
                                           convolveTemplate=convolveTemplate,
@@ -505,6 +507,8 @@ def chooseConvolutionMethod(self, template, science):
         RuntimeError
             If an unsupported convolution mode is supplied.
         """
+        if self.usePreconvolution:
+            raise RuntimeError("Choosing a convolution method is incompatible with preconvolution!")
         if self.config.mode == "auto":
             convolveTemplate = _shapeTest(template,
                                           science,
@@ -529,7 +533,7 @@ def chooseConvolutionMethod(self, template, science):
         return convolveTemplate
 
     def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False):
-        """Construct the PSF-matching kernel.
+        """Construct the PSF-matching kernel. Not used for preconvolution.
 
         Parameters
         ----------
@@ -561,7 +565,8 @@ def runMakeKernel(self, template, science, sources=None, convolveTemplate=True,
                 Sources from the input catalog that were used to construct the
                 PSF-matching kernel.
         """
-
+        if self.usePreconvolution:
+            raise RuntimeError("Incorrect matching kernel calculation configured.")
         if convolveTemplate:
             reference = template
             target = science
@@ -898,8 +903,7 @@ def _convolveExposure(self, exposure, kernel, convolutionControl,
         else:
             return convolvedExposure[bbox]
 
-    def _sourceSelector(self, template, science, sources, fallback=False, preconvolved=False,
-                        scienceFwhmPix=None):
+    def _sourceSelector(self, template, science, sources, fallback=False, preconvolved=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 checking the science and template
@@ -921,11 +925,6 @@ def _sourceSelector(self, template, science, sources, fallback=False, preconvolv
         preconvolved : `bool`, optional
             If set, exclude a wider buffer around the edge of the image to
             account for an extra convolution.
-        scienceFwhmPix : `float`, optional
-            Effective FWHM of the science image in pixels to use for sizing
-            the kernel basis and source footprints. If `None`, the
-            width is derived from ``self.sciencePsfSize`` (boosted by
-            ``sqrt(2)`` in the ``preconvolved`` case).
 
         Returns
         -------
@@ -943,9 +942,10 @@ def _sourceSelector(self, template, science, sources, fallback=False, preconvolv
             selected = self.fallbackSourceSelector.selectSources(sources).selected
         else:
             selected = self.sourceSelector.selectSources(sources).selected
-        if scienceFwhmPix is None:
-            scienceFwhmPix = self.sciencePsfSize*np.sqrt(2) if preconvolved else self.sciencePsfSize
-        kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, scienceFwhmPix)[0].getWidth()
+        # It is OK to use just self.matchedPsfSize for the science PSF here,
+        # since we are just using it to calculate the size of the matching
+        # kernel.
+        kSize = self.makeKernel.makeKernelBasisList(self.templatePsfSize, self.matchedPsfSize)[0].getWidth()
         selectSources = sources[selected].copy(deep=True)
         # Set the footprints, to be used in `makeKernel` and `checkMask`.
         kernelSources = setSourceFootprints(selectSources, kernelSize=kSize)
@@ -1223,19 +1223,20 @@ def run(self, template, science, sources, visitSummary=None):
                 Final kernel used to PSF-match the template to the science
                 image.
         """
+        self.usePreconvolution = True
         self._prepareInputs(template, science, visitSummary=visitSummary)
+        self.matchedPsfSize = self.sciencePsfSize*np.sqrt(2)
 
         convolutionKernel = self._makePreconvolutionKernel(science.psf)
         matchedScience = self._convolveExposure(science, convolutionKernel, self.convolutionControl,
                                                 interpolateBadMaskPlanes=True)
         self.metadata["convolvedExposure"] = "Preconvolution"
 
-        self.preconvolvedSciencePsfSize = self.sciencePsfSize*np.sqrt(2)
-        self.log.info("Preconvolved science PSF FWHM: %f pixels", self.preconvolvedSciencePsfSize)
-        self.metadata["preconvolvedSciencePsfSize"] = self.preconvolvedSciencePsfSize
+        self.matchedPsfSize = self.sciencePsfSize*np.sqrt(2)
+        self.log.info("Preconvolved science PSF FWHM: %f pixels", self.matchedPsfSize)
+        self.metadata["preconvolvedSciencePsfSize"] = self.matchedPsfSize
         try:
-            kernelSources = self._sourceSelector(template, matchedScience, sources, preconvolved=True,
-                                                 scienceFwhmPix=self.preconvolvedSciencePsfSize)
+            kernelSources = self._sourceSelector(template, matchedScience, sources, preconvolved=True)
             subtractResults = self.runPreconvolve(template, science, matchedScience,
                                                   kernelSources, convolutionKernel)
 
@@ -1367,7 +1368,7 @@ def runPreconvolve(self, template, science, matchedScience, kernelSources, preCo
         kernelResult = self.makeKernel.run(template[innerBBox], matchedScience[innerBBox], kernelSources,
                                            preconvolved=True,
                                            templateFwhmPix=self.templatePsfSize,
-                                           scienceFwhmPix=self.preconvolvedSciencePsfSize)
+                                           scienceFwhmPix=self.matchedPsfSize)
 
         matchedTemplate = self._convolveExposure(template, kernelResult.psfMatchingKernel,
                                                  self.convolutionControl,
@@ -1572,9 +1573,11 @@ def run(self, template, science, visitSummary=None, inputPsfMatchingKernel=None)
             Raised if fraction of good pixels, defined as not having NO_DATA
             set, is less then the configured requiredTemplateFraction
         """
+        self.usePreconvolution = False
         self._prepareInputs(template, science, visitSummary=visitSummary)
 
         convolveTemplate = self.chooseConvolutionMethod(template, science)
+        self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize
 
         if self.config.useExistingKernel:
             psfMatchingKernel = inputPsfMatchingKernel

From 703b5776012e2c273f93e65a253b4b1416f95992 Mon Sep 17 00:00:00 2001
From: Ian Sullivan <sullii@uw.edu>
Date: Sun, 3 May 2026 19:42:40 -0700
Subject: [PATCH 24/24] Refactor CR detection and background subtraction

Use common code for both standard and Score versions.
---
 python/lsst/ip/diffim/detectAndMeasure.py | 203 +++++++++++++++-------
 1 file changed, 139 insertions(+), 64 deletions(-)

diff --git a/python/lsst/ip/diffim/detectAndMeasure.py b/python/lsst/ip/diffim/detectAndMeasure.py
index 0009b613..1cd9bdda 100644
--- a/python/lsst/ip/diffim/detectAndMeasure.py
+++ b/python/lsst/ip/diffim/detectAndMeasure.py
@@ -664,20 +664,16 @@ def run(self, science, matchedTemplate, difference, kernelSources=None,
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
-        if self.config.doSubtractBackground:
-            # Run background subtraction before clearing the mask planes
-            detectionExposure = difference.clone()
-            background = self.subtractInitialBackground.run(detectionExposure).background
-        else:
-            detectionExposure = difference
-            background = afwMath.BackgroundList()
+        # Run background subtraction (if configured) before clearing the mask
+        # planes.
+        detectionExposure, background = self._subtractInitialBackground(differenceExposure=difference)
 
         self._prepareInputs(detectionExposure)
 
-        if self.config.doFindCosmicRays and not self.config.doSubtractBackground:
-            # Detect and interpolate over any remaining cosmic rays
-            # This only needs to run once, right before the final detection.
-            self.findAndMaskCosmicRays(detectionExposure)
+        # Detect and interpolate over any remaining cosmic rays. When
+        # background subtraction is enabled, this is deferred to the second
+        # pass so it runs only once on a properly background-subtracted image.
+        self._findInitialCosmicRays(differenceExposure=difference)
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't
@@ -691,23 +687,12 @@ def run(self, science, matchedTemplate, difference, kernelSources=None,
         )
 
         if self.config.doSubtractBackground:
-            # Run background subtraction again after detecting peaks
-            # but before measurement
-            # First update the mask using the detection image
-            difference.setMask(detectionExposure.mask)
-            background = self.subtractFinalBackground.run(difference).background
-
-            if self.config.doFindCosmicRays:
-                # Detect and interpolate over any remaining cosmic rays
-                self.findAndMaskCosmicRays(difference)
-
-            # Re-run detection to get final footprints
-            table = afwTable.SourceTable.make(self.schema)
-            results = self.detection.run(
-                table=table,
-                exposure=difference,
+            # Refit the background using the first-pass detection mask, then
+            # find cosmic rays and re-run detection on the difference.
+            results, background = self._subtractFinalBackgroundAndRedetect(
+                differenceExposure=difference,
+                firstPassDetectionExposure=detectionExposure,
                 doSmooth=True,
-                background=background
             )
         measurementResults.differenceBackground = background
 
@@ -755,6 +740,110 @@ def _prepareInputs(self, difference):
                 mask.addMaskPlane(mp)
         mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)
 
+    def _subtractInitialBackground(self, *, differenceExposure, scoreExposure=None):
+        """Optionally fit and subtract an initial background.
+
+        Parameters
+        ----------
+        differenceExposure : `lsst.afw.image.ExposureF`
+            Template-subtracted difference image. Used only as the input to
+            background fitting; never modified.
+        scoreExposure : `lsst.afw.image.ExposureF`, optional
+            Maximum-likelihood score image, when detection is to be run on
+            the score rather than on the difference.
+
+        Returns
+        -------
+        detectionExposure : `lsst.afw.image.ExposureF`
+            Exposure to use for first-pass detection.
+        background : `lsst.afw.math.BackgroundList`
+            The fit background, or an empty list when background subtraction
+            is disabled.
+        """
+        if not self.config.doSubtractBackground:
+            return (scoreExposure if scoreExposure is not None else differenceExposure,
+                    afwMath.BackgroundList())
+        if scoreExposure is None:
+            detectionExposure = differenceExposure.clone()
+            background = self.subtractInitialBackground.run(detectionExposure).background
+        else:
+            background = self.subtractInitialBackground.run(differenceExposure.clone()).background
+            detectionExposure = scoreExposure.clone()
+            detectionExposure.image -= background.getImage()
+        return detectionExposure, background
+
+    def _findInitialCosmicRays(self, *, differenceExposure, scoreExposure=None):
+        """Find cosmic rays before the first detection pass, if requested.
+
+        Cosmic-ray finding is only performed in the first pass when background
+        subtraction is disabled; with ``doSubtractBackground`` set it is
+        deferred to `_subtractFinalBackgroundAndRedetect` so that it runs on
+        the final background-subtracted image.
+
+        Parameters
+        ----------
+        differenceExposure : `lsst.afw.image.ExposureF`
+            Difference exposure on which cosmic rays are detected; its mask
+            is updated in place.
+        scoreExposure : `lsst.afw.image.ExposureF`, optional
+            Maximum-likelihood score image, used as the detection exposure
+            in the score-image variant. When supplied, the new CR mask bits
+            are OR'd onto its mask.
+        """
+        if not self.config.doFindCosmicRays or self.config.doSubtractBackground:
+            return
+        self.findAndMaskCosmicRays(differenceExposure)
+        if scoreExposure is not None:
+            scoreExposure.mask.array |= differenceExposure.mask.array
+
+    def _subtractFinalBackgroundAndRedetect(self, *, differenceExposure,
+                                            firstPassDetectionExposure,
+                                            scoreExposure=None, doSmooth=True):
+        """Refit the background using the first-pass detection mask, then
+        re-find cosmic rays (if requested) and run final detection.
+
+        Parameters
+        ----------
+        differenceExposure : `lsst.afw.image.ExposureF`
+            Difference exposure used to fit the final background and find
+            cosmic rays. Modified in place.
+        firstPassDetectionExposure : `lsst.afw.image.ExposureF`
+            Exposure that was used for the first detection pass; its mask
+            (carrying first-pass detection bits) is copied onto
+            ``differenceExposure`` before refitting.
+        scoreExposure : `lsst.afw.image.ExposureF`, optional
+            Maximum-likelihood score image. If supplied, the fitted background
+            is also subtracted from it, and it is used for detection.
+        doSmooth : `bool`, optional
+            Whether the detection task should smooth the input before
+            thresholding. Should be `True` unless the score image is used.
+
+        Returns
+        -------
+        detectResults : `lsst.pipe.base.Struct`
+            Positive and negative peaks from detection, to be passed to source
+            measurement.
+        background : `lsst.afw.math.BackgroundList`
+            The refit background.
+        """
+        differenceExposure.setMask(firstPassDetectionExposure.mask)
+        background = self.subtractFinalBackground.run(differenceExposure).background
+        if scoreExposure is not None:
+            scoreExposure.image -= background.getImage()
+        if self.config.doFindCosmicRays:
+            self.findAndMaskCosmicRays(differenceExposure)
+            if scoreExposure is not None:
+                scoreExposure.mask.array |= differenceExposure.mask.array
+        detectionExposure = scoreExposure if scoreExposure is not None else differenceExposure
+        table = afwTable.SourceTable.make(self.schema)
+        detectResults = self.detection.run(
+            table=table,
+            exposure=detectionExposure,
+            doSmooth=doSmooth,
+            background=background,
+        )
+        return detectResults, background
+
     def processResults(self, science, matchedTemplate, difference, sources, idFactory,
                        kernelSources=None, positives=None, negatives=None, measurementResults=None):
         """Measure and process the results of source detection.
@@ -1480,31 +1569,26 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
         if idFactory is None:
             idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
 
-        if self.config.doSubtractBackground:
-            # Background is measured on the difference image, and then
-            # subtracted from both the difference and score images.
-            # The preconvolution kernel is normalized to 1, so the
-            # same background level applies to both.
-            background = self.subtractInitialBackground.run(difference.clone()).background
-            detectionScoreExposure = scoreExposure.clone()
-            detectionScoreExposure.image -= background.getImage()
-        else:
-            detectionScoreExposure = scoreExposure
-            background = afwMath.BackgroundList()
+        # Background is measured on the difference image, and then subtracted
+        # from both the difference and score images. The preconvolution kernel
+        # is normalized to 1, so the same background level applies to both.
+        detectionScoreExposure, background = self._subtractInitialBackground(
+            differenceExposure=difference, scoreExposure=scoreExposure,
+        )
 
         self._prepareInputs(detectionScoreExposure)
 
-        if self.config.doFindCosmicRays and not self.config.doSubtractBackground:
-            # Detect cosmic rays on the difference image and propagate the mask
-            # to the score image.
-            self.findAndMaskCosmicRays(difference)
-            scoreExposure.mask.array |= difference.mask.array
+        # Cosmic rays are detected on the difference image and the resulting
+        # mask is propagated to the score image used for detection.
+        self._findInitialCosmicRays(
+            differenceExposure=difference, scoreExposure=scoreExposure,
+        )
 
         # Don't use the idFactory until after deblend+merge, so that we aren't
         # generating ids that just get thrown away (footprint merge doesn't
         # know about past ids).
         table = afwTable.SourceTable.make(self.schema)
-        results = self.detection.run(
+        detectResults = self.detection.run(
             table=table,
             exposure=detectionScoreExposure,
             doSmooth=False,
@@ -1513,22 +1597,13 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
 
         if self.config.doSubtractBackground:
             # Refine the background with detected peaks flagged in the mask.
-            # The refinement is fit on the difference and also applied
-            # to the score image.
-            difference.setMask(detectionScoreExposure.mask)
-            background = self.subtractFinalBackground.run(difference).background
-            scoreExposure.image -= background.getImage()
-
-            if self.config.doFindCosmicRays:
-                self.findAndMaskCosmicRays(difference)
-                scoreExposure.mask.array |= difference.mask.array
-
-            table = afwTable.SourceTable.make(self.schema)
-            results = self.detection.run(
-                table=table,
-                exposure=scoreExposure,
+            # The refinement is fit on the difference and also applied to the
+            # score image, where the second-pass detection runs.
+            detectResults, background = self._subtractFinalBackgroundAndRedetect(
+                differenceExposure=difference,
+                firstPassDetectionExposure=detectionScoreExposure,
+                scoreExposure=scoreExposure,
                 doSmooth=False,
-                background=background,
             )
         measurementResults.differenceBackground = background
         measurementResults.scoreMeasuredExposure = scoreExposure
@@ -1538,8 +1613,8 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
 
         if self.config.doDeblend:
             sources, positives, negatives = self._deblend(difference,
-                                                          results.positive,
-                                                          results.negative)
+                                                          detectResults.positive,
+                                                          detectResults.negative)
 
             self.processResults(science, matchedTemplate, difference,
                                 sources, idFactory, kernelSources,
@@ -1549,11 +1624,11 @@ def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources
 
         else:
             positives = afwTable.SourceCatalog(self.schema)
-            results.positive.makeSources(positives)
+            detectResults.positive.makeSources(positives)
             negatives = afwTable.SourceCatalog(self.schema)
-            results.negative.makeSources(negatives)
+            detectResults.negative.makeSources(negatives)
             self.processResults(science, matchedTemplate, difference,
-                                results.sources, idFactory, kernelSources,
+                                detectResults.sources, idFactory, kernelSources,
                                 positives=positives,
                                 negatives=negatives,
                                 measurementResults=measurementResults)