diff --git a/python/lsst/ip/diffim/subtractImages.py b/python/lsst/ip/diffim/subtractImages.py
index 84f2933c..3482e378 100644
--- a/python/lsst/ip/diffim/subtractImages.py
+++ b/python/lsst/ip/diffim/subtractImages.py
@@ -193,12 +193,33 @@ class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config):
         doc="Minimum signal to noise ratio of detected sources "
         "to use for calculating the PSF matching kernel."
     )
+    detectionThresholdMax = lsst.pex.config.Field(
+        dtype=float,
+        default=500,
+        doc="Maximum signal to noise ratio of detected sources "
+        "to use for calculating the PSF matching kernel."
+    )
+    maxKernelSources = lsst.pex.config.Field(
+        dtype=int,
+        default=1000,
+        doc="Maximum number of sources to use for calculating the PSF matching kernel."
+        "Set to -1 to disable."
+    )
+    minKernelSources = lsst.pex.config.Field(
+        dtype=int,
+        default=3,
+        doc="Minimum number of sources needed for calculating the PSF matching kernel."
+    )
     badSourceFlags = lsst.pex.config.ListField(
         dtype=str,
         doc="Flags that, if set, the associated source should not "
         "be used to determine the PSF matching kernel.",
         default=("sky_source", "slot_Centroid_flag",
-                 "slot_ApFlux_flag", "slot_PsfFlux_flag", ),
+                 "slot_ApFlux_flag", "slot_PsfFlux_flag",
+                 "base_PixelFlags_flag_interpolated",
+                 "base_PixelFlags_flag_saturated",
+                 "base_PixelFlags_flag_bad",
+                 ),
     )
     excludeMaskPlanes = lsst.pex.config.ListField(
         dtype=str,
@@ -416,7 +437,9 @@ def run(self, template, science, sources, finalizedPsfApCorrCatalog=None,
             raise RuntimeError("Cannot handle AlardLuptonSubtract mode: %s", self.config.mode)
 
         try:
-            selectSources = self._sourceSelector(sources, science.mask)
+            sourceMask = science.mask.clone()
+            sourceMask.array |= template[science.getBBox()].mask.array
+            selectSources = self._sourceSelector(sources, sourceMask)
             if convolveTemplate:
                 self.metadata.add("convolvedExposure", "Template")
                 subtractResults = self.runConvolveTemplate(template, science, selectSources)
@@ -752,20 +775,31 @@ def _sourceSelector(self, sources, mask):
                 flags *= ~sources[flag]
             except Exception as e:
                 self.log.warning("Could not apply source flag: %s", e)
-        sToNFlag = (sources.getPsfInstFlux()/sources.getPsfInstFluxErr()) > self.config.detectionThreshold
+        signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr()
+        sToNFlag = signalToNoise > self.config.detectionThreshold
         flags *= sToNFlag
+        sToNFlagMax = signalToNoise < self.config.detectionThresholdMax
+        flags *= sToNFlagMax
         flags *= self._checkMask(mask, sources, self.config.excludeMaskPlanes)
-        selectSources = sources[flags]
+        selectSources = sources[flags].copy(deep=True)
+        if (len(selectSources) > self.config.maxKernelSources) & (self.config.maxKernelSources > 0):
+            signalToNoise = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
+            indices = np.argsort(signalToNoise)
+            indices = indices[-self.config.maxKernelSources:]
+            flags = np.zeros(len(selectSources), dtype=bool)
+            flags[indices] = True
+            selectSources = selectSources[flags].copy(deep=True)
+
         self.log.info("%i/%i=%.1f%% of sources selected for PSF matching from the input catalog",
                       len(selectSources), len(sources), 100*len(selectSources)/len(sources))
-        if len(selectSources) < self.config.makeKernel.nStarPerCell:
+        if len(selectSources) < self.config.minKernelSources:
             self.log.error("Too few sources to calculate the PSF matching kernel: "
                            "%i selected but %i needed for the calculation.",
-                           len(selectSources), self.config.makeKernel.nStarPerCell)
+                           len(selectSources), self.config.minKernelSources)
             raise RuntimeError("Cannot compute PSF matching kernel: too few sources selected.")
         self.metadata.add("nPsfSources", len(selectSources))
 
-        return selectSources.copy(deep=True)
+        return selectSources
 
     @staticmethod
     def _checkMask(mask, sources, excludeMaskPlanes):
diff --git a/tests/test_subtractTask.py b/tests/test_subtractTask.py
index ea83c84c..1e3eeefb 100644
--- a/tests/test_subtractTask.py
+++ b/tests/test_subtractTask.py
@@ -402,6 +402,52 @@ def test_few_sources(self):
                                     "Cannot compute PSF matching kernel: too few sources selected."):
             task.run(template, science, sources)
 
+    def test_kernel_source_selector(self):
+        """Test with only 1 source, to check that we get a useful error.
+        """
+        xSize = 256
+        ySize = 256
+        nSourcesSimulated = 20
+        science, sources = makeTestImage(psfSize=2.4, nSrc=nSourcesSimulated,
+                                         xSize=xSize, ySize=ySize)
+        template, _ = makeTestImage(psfSize=2.0, nSrc=nSourcesSimulated,
+                                    xSize=xSize, ySize=ySize, doApplyCalibration=True)
+        badSourceFlag = "slot_Centroid_flag"
+
+        def _run_and_check_sources(sourcesIn, maxKernelSources=1000, minKernelSources=3):
+            sources = sourcesIn.copy(deep=True)
+            # Verify that source flags are not set in the input catalog
+            self.assertEqual(np.sum(sources[badSourceFlag]), 0)
+            config = subtractImages.AlardLuptonSubtractTask.ConfigClass()
+            config.badSourceFlags = [badSourceFlag, ]
+            config.maxKernelSources = maxKernelSources
+            config.minKernelSources = minKernelSources
+
+            task = subtractImages.AlardLuptonSubtractTask(config=config)
+            nSources = len(sources)
+            # Flag a third of the sources
+            sources[0:: 3][badSourceFlag] = True
+            nBadSources = np.sum(sources[badSourceFlag])
+            if maxKernelSources > 0:
+                nGoodSources = np.minimum(nSources - nBadSources, maxKernelSources)
+            else:
+                nGoodSources = nSources - nBadSources
+
+            signalToNoise = sources.getPsfInstFlux()/sources.getPsfInstFluxErr()
+            signalToNoise = signalToNoise[~sources[badSourceFlag]]
+            signalToNoise.sort()
+            selectSources = task._sourceSelector(sources, science.mask)
+            self.assertEqual(nGoodSources, len(selectSources))
+            signalToNoiseOut = selectSources.getPsfInstFlux()/selectSources.getPsfInstFluxErr()
+            signalToNoiseOut.sort()
+            self.assertFloatsAlmostEqual(signalToNoise[-nGoodSources:], signalToNoiseOut)
+
+        _run_and_check_sources(sources)
+        _run_and_check_sources(sources, maxKernelSources=len(sources)//3)
+        _run_and_check_sources(sources, maxKernelSources=-1)
+        with self.assertRaises(RuntimeError):
+            _run_and_check_sources(sources, minKernelSources=1000)
+
     def test_order_equal_images(self):
         """Verify that the result is the same regardless of convolution mode
         if the images are equivalent.