diff --git a/python/lsst/ip/diffim/imageDecorrelation.py b/python/lsst/ip/diffim/imageDecorrelation.py
index 13fe1748..604c0981 100644
--- a/python/lsst/ip/diffim/imageDecorrelation.py
+++ b/python/lsst/ip/diffim/imageDecorrelation.py
@@ -28,10 +28,12 @@
 import lsst.meas.algorithms as measAlg
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
+from lsst.pex.exceptions import InvalidParameterError
 from lsst.utils.timer import timeMethod
 
 from .imageMapReduce import (ImageMapReduceConfig, ImageMapReduceTask,
                              ImageMapper)
+from .utils import computeAveragePsf
 
 __all__ = ("DecorrelateALKernelTask", "DecorrelateALKernelConfig",
            "DecorrelateALKernelMapper", "DecorrelateALKernelMapReduceConfig",
@@ -224,6 +226,9 @@ def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatching
             variance = scienceExposure.variance.array
             # Variance plane of the convolved image, before convolution.
             targetVariance = templateExposure.variance.array
+            # If the template is convolved, the PSF of the difference image is
+            # that of the science image
+            psfImg = scienceExposure.getPsf().computeKernelImage(geom.Point2D(xcen, ycen))
         else:
             # We convolved the science image
             self.log.info("Decorrelation after science image convolution")
@@ -233,6 +238,15 @@ def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatching
             variance = templateExposure.variance.array
             # Variance plane of the convolved image, before convolution.
             targetVariance = scienceExposure.variance.array
+            # If the science image is convolved, the PSF of the difference image
+            # is that of the template image, and will be a CoaddPsf which might
+            # not be defined for the entire image.
+            # Try the simple calculation first, and fall back on a slower method
+            # if the PSF of the template is not defined in the center.
+            try:
+                psfImg = templateExposure.getPsf().computeKernelImage(geom.Point2D(xcen, ycen))
+            except InvalidParameterError:
+                psfImg = computeAveragePsf(templateExposure, psfExposureBuffer=0.05, psfExposureGrid=100)
 
         # The maximal correction value converges to sqrt(targetVarianceMean/varianceMean).
         # Correction divergence warning if the correction exceeds 4 orders of magnitude.
@@ -247,7 +261,6 @@ def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatching
 
         kArr = kimg.array
         diffimShape = subtractedExposure.image.array.shape
-        psfImg = subtractedExposure.getPsf().computeKernelImage(geom.Point2D(xcen, ycen))
         psfShape = psfImg.array.shape
 
         if preConvMode:
diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py
index 6d3c3cf1..c5681312 100644
--- a/python/lsst/ip/diffim/utils.py
+++ b/python/lsst/ip/diffim/utils.py
@@ -1116,7 +1116,7 @@ def getPsfFwhm(psf, average=True, position=None):
 
 def evaluateMeanPsfFwhm(exposure: afwImage.Exposure,
                         fwhmExposureBuffer: float, fwhmExposureGrid: int) -> float:
-    """Get the median PSF FWHM by evaluating it on a grid within an exposure.
+    """Get the mean PSF FWHM by evaluating it on a grid within an exposure.
 
     Parameters
     ----------
@@ -1142,7 +1142,8 @@ def evaluateMeanPsfFwhm(exposure: afwImage.Exposure,
 
     See Also
     --------
-    getPsfFwhm
+    `getPsfFwhm`
+    `computeAveragePsf`
     """
 
     psf = exposure.psf
@@ -1173,6 +1174,71 @@ def evaluateMeanPsfFwhm(exposure: afwImage.Exposure,
     return np.nanmean(width)
 
 
+def computeAveragePsf(exposure: afwImage.Exposure,
+                      psfExposureBuffer: float, psfExposureGrid: int) -> afwImage.ImageD:
+    """Get the average PSF by evaluating it on a grid within an exposure.
+
+    Parameters
+    ----------
+    exposure : `~lsst.afw.image.Exposure`
+        The exposure for which the average PSF is to be computed.
+        The exposure must contain a `psf` attribute.
+    psfExposureBuffer : `float`
+        Fractional buffer margin to be left out of all sides of the image
+        during the construction of the grid to compute average PSF in an
+        exposure.
+    psfExposureGrid : `int`
+        Grid size to compute the average PSF in an exposure.
+
+    Returns
+    -------
+    psfImage : `~lsst.afw.image.Image`
+        The average PSF across the exposure.
+
+    Raises
+    ------
+    ValueError
+        Raised if the PSF cannot be computed at any of the grid points.
+
+    See Also
+    --------
+    `evaluateMeanPsfFwhm`
+    """
+
+    psf = exposure.psf
+
+    bbox = exposure.getBBox()
+    xmax, ymax = bbox.getMax()
+    xmin, ymin = bbox.getMin()
+
+    xbuffer = psfExposureBuffer*(xmax-xmin)
+    ybuffer = psfExposureBuffer*(ymax-ymin)
+
+    nImg = 0
+    psfArray = None
+    for (x, y) in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, psfExposureGrid),
+                                    np.linspace(ymin+ybuffer, ymax-ybuffer, psfExposureGrid)
+                                    ):
+        pos = geom.Point2D(x, y)
+        try:
+            singleImage = psf.computeKernelImage(pos)
+        except InvalidParameterError:
+            _LOG.debug("Unable to compute PSF image at position (%f, %f).", x, y)
+            continue
+
+        if psfArray is None:
+            psfArray = singleImage.array
+        else:
+            psfArray += singleImage.array
+        nImg += 1
+
+    if psfArray is None:
+        raise ValueError("Unable to compute PSF image at any position on the exposure.")
+
+    psfImage = afwImage.ImageD(psfArray/nImg)
+    return psfImage
+
+
 def detectTestSources(exposure):
     """Minimal source detection wrapper suitable for unit tests.