diff --git a/python/lsst/pipe/tasks/calibrateImage.py b/python/lsst/pipe/tasks/calibrateImage.py
index 1c7b9aada..427545364 100644
--- a/python/lsst/pipe/tasks/calibrateImage.py
+++ b/python/lsst/pipe/tasks/calibrateImage.py
@@ -319,7 +319,9 @@ def setDefaults(self):
         self.photometry.match.sourceSelection.retarget(sourceSelector.NullSourceSelectorTask)
 
         # All sources should be good for PSF summary statistics.
+        # TODO: These should both be changed to calib_psf_used with DM-41640.
         self.compute_summary_stats.starSelection = "calib_photometry_used"
+        self.compute_summary_stats.starSelector.flags.good = ["calib_photometry_used"]
 
 
 class CalibrateImageTask(pipeBase.PipelineTask):
diff --git a/python/lsst/pipe/tasks/computeExposureSummaryStats.py b/python/lsst/pipe/tasks/computeExposureSummaryStats.py
index cb2d01c32..8c0e1e8ca 100644
--- a/python/lsst/pipe/tasks/computeExposureSummaryStats.py
+++ b/python/lsst/pipe/tasks/computeExposureSummaryStats.py
@@ -34,6 +34,7 @@
 import lsst.afw.math as afwMath
 import lsst.afw.image as afwImage
 import lsst.geom as geom
+from lsst.meas.algorithms import ScienceSourceSelectorTask
 from lsst.utils.timer import timeMethod
 
 
@@ -55,9 +56,13 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config):
         default=("NO_DATA", "SUSPECT"),
     )
     starSelection = pexConfig.Field(
-        doc="Field to select sources to be used in the PSF statistics computation.",
+        doc="Field to select full list of sources used for PSF modeling.",
         dtype=str,
-        default="calib_psf_used"
+        default="calib_psf_used",
+    )
+    starSelector = pexConfig.ConfigurableField(
+        target=ScienceSourceSelectorTask,
+        doc="Selection of sources to compute PSF star statistics.",
     )
     starShape = pexConfig.Field(
         doc="Base name of columns to use for the source shape in the PSF statistics computation.",
@@ -88,6 +93,27 @@ class ComputeExposureSummaryStatsConfig(pexConfig.Config):
         default=("BAD", "CR", "EDGE", "INTRP", "NO_DATA", "SAT", "SUSPECT"),
     )
 
+    def setDefaults(self):
+        super().setDefaults()
+
+        self.starSelector.setDefaults()
+        self.starSelector.doFlags = True
+        self.starSelector.doSignalToNoise = True
+        self.starSelector.doUnresolved = False
+        self.starSelector.doIsolated = False
+        self.starSelector.doRequireFiniteRaDec = False
+        self.starSelector.doRequirePrimary = False
+
+        self.starSelector.signalToNoise.minimum = 50.0
+        self.starSelector.signalToNoise.maximum = 1000.0
+
+        self.starSelector.flags.bad = ["slot_Shape_flag", "slot_PsfFlux_flag"]
+        # Select stars used for PSF modeling.
+        self.starSelector.flags.good = ["calib_psf_used"]
+
+        self.starSelector.signalToNoise.fluxField = "slot_PsfFlux_instFlux"
+        self.starSelector.signalToNoise.errField = "slot_PsfFlux_instFluxErr"
+
 
 class ComputeExposureSummaryStatsTask(pipeBase.Task):
     """Task to compute exposure summary statistics.
@@ -129,6 +155,11 @@ class ComputeExposureSummaryStatsTask(pipeBase.Task):
     ConfigClass = ComputeExposureSummaryStatsConfig
     _DefaultName = "computeExposureSummaryStats"
 
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self.makeSubtask("starSelector")
+
     @timeMethod
     def run(self, exposure, sources, background):
         """Measure exposure statistics from the exposure, sources, and
@@ -243,10 +274,14 @@ def update_psf_stats(self, summary, psf, bbox, sources=None, image_mask=None, so
             # good sources).
             return
 
-        psf_mask = sources[self.config.starSelection] & (~sources[self.config.starShape + '_flag'])
-        nPsfStar = psf_mask.sum()
+        # Count the total number of psf stars used (prior to stats selection).
+        nPsfStar = sources[self.config.starSelection].sum()
+        summary.nPsfStar = int(nPsfStar)
+
+        psf_mask = self.starSelector.run(sources).selected
+        nPsfStarsUsedInStats = psf_mask.sum()
 
-        if nPsfStar == 0:
+        if nPsfStarsUsedInStats == 0:
             # No stars to measure statistics, so we must return the defaults
             # of 0 stars and NaN values.
             return
@@ -263,12 +298,15 @@ def update_psf_stats(self, summary, psf, bbox, sources=None, image_mask=None, so
         psfYY = psf_cat[self.config.psfShape + '_yy']
         psfXY = psf_cat[self.config.psfShape + '_xy']
 
-        starSize = (starXX*starYY - starXY**2.)**0.25
+        # Use the trace radius for the star size.
+        starSize = np.sqrt(starXX/2. + starYY/2.)
+
         starE1 = (starXX - starYY)/(starXX + starYY)
         starE2 = 2*starXY/(starXX + starYY)
         starSizeMedian = np.median(starSize)
 
-        psfSize = (psfXX*psfYY - psfXY**2)**0.25
+        # Use the trace radius for the psf size.
+        psfSize = np.sqrt(psfXX/2. + psfYY/2.)
         psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
         psfE2 = 2*psfXY/(psfXX + psfYY)
 
@@ -279,9 +317,8 @@ def update_psf_stats(self, summary, psf, bbox, sources=None, image_mask=None, so
 
         psfStarDeltaSizeMedian = np.median(starSize - psfSize)
         psfStarDeltaSizeScatter = sigmaMad(starSize - psfSize, scale='normal')
-        psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian**2.
+        psfStarScaledDeltaSizeScatter = psfStarDeltaSizeScatter/starSizeMedian
 
-        summary.nPsfStar = int(nPsfStar)
         summary.psfStarDeltaE1Median = float(psfStarDeltaE1Median)
         summary.psfStarDeltaE2Median = float(psfStarDeltaE2Median)
         summary.psfStarDeltaE1Scatter = float(psfStarDeltaE1Scatter)
diff --git a/python/lsst/pipe/tasks/selectImages.py b/python/lsst/pipe/tasks/selectImages.py
index 3b9fb246d..f07a6b45d 100644
--- a/python/lsst/pipe/tasks/selectImages.py
+++ b/python/lsst/pipe/tasks/selectImages.py
@@ -258,7 +258,7 @@ class PsfWcsSelectImagesConfig(pipeBase.PipelineTaskConfig,
     maxScaledSizeScatter = pexConfig.Field(
         doc="Maximum scatter in the size residuals, scaled by the median size",
         dtype=float,
-        default=0.009,
+        default=0.019,
         optional=True,
     )
     maxPsfTraceRadiusDelta = pexConfig.Field(