Retry warped PSF evaluation if transform fails

for initial rough PSF size estimate

Author: yalsayyad

python/lsst/ip/diffim/modelPsfMatch.py

@@ -45,6 +45,55 @@
 sigma2fwhm = 2.*np.sqrt(2.*np.log(2.))
 
 
+def _safeEstimateRoughShape(psf, position, jitter=2.0, maxTries=5, log=None):
+
+    """
+    Try to compute the shape of `psf` at `position`.
+
+    If it fails, perturb the position until success or maxTries.
+
+    Parameters
+    ----------
+    psf : `lsst.afw.detection.Psf`
+        The PSF object.
+    position : `lsst.geom.Point2D`
+        The nominal position at which to evaluate the PSF shape.
+    jitter : `float`
+        Size of random offset in pixels to try when retrying.
+    maxTries : `int`
+        Maximum number of attempts (including the original).
+    log : `lsst.utils.logging.LsstLogAdapter` optional
+        logger to use for logging retries.
+
+    Returns
+    -------
+    shape : lsst.afw.geom.ellipses.Quadrupole
+        The shape of the PSF at the (possibly perturbed) position.
+    """
+    try:
+        return psf.computeShape(position)
+    except pexExceptions.Exception as err:
+        if log:
+            log.info(f"safeEstimateRoughShape: initial evaluation of PSF failed with message: {str(err)}."
+                     " Retrying nearby")
+        firstErr = err
+
+    # random offsets must be deterministic
+    offsets = np.random.default_rng(seed=40).uniform(-jitter, jitter, size=(maxTries - 1, 2))
+    candidates = [position + geom.Extent2D(dx, dy) for dx, dy in offsets]
+    for i, candidate in enumerate(candidates):
+        try:
+            shape = psf.computeShape(candidate)
+            if log:
+                log.info(f"safeEstimateRoughShape succeeded on try {i + 2} at position {candidate}")
+            return shape
+        except pexExceptions.Exception:
+            continue
+
+    # If nothing worked, raise AlgorithmError from original position
+    raise PsfComputeShapeError(position) from firstErr
+
+
 def nextOddInteger(x):
     nextInt = int(np.ceil(x))
     return nextInt + 1 if nextInt%2 == 0 else nextInt
@@ -175,15 +224,12 @@ def run(self, exposure, referencePsfModel, kernelSum=1.0):
         # exposure's bounding box in DM-32756.
         sciAvgPos = exposure.getPsf().getAveragePosition()
         modelAvgPos = referencePsfModel.getAveragePosition()
-        try:
-            fwhmScience = exposure.getPsf().computeShape(sciAvgPos).getDeterminantRadius()*sigma2fwhm
-        except pexExceptions.RangeError as err:
-            raise WarpedPsfTransformTooBigError(
-                f"Unable to compute the FWHM of the science Psf at {sciAvgPos}"
-                "due to an unexpectedly large transform."
-            ) from err
-        except pexExceptions.Exception as err:
-            raise PsfComputeShapeError(sciAvgPos) from err
+
+        # TODO If DM-52537 eliminates the need for _safeEstimateRoughShape,
+        # the remove it, and replace next line with a direct call to
+        # computeShape in a try/except block.
+        fwhmScience = _safeEstimateRoughShape(exposure.getPsf(), sciAvgPos,
+                                              log=self.log).getDeterminantRadius()*sigma2fwhm
         fwhmModel = referencePsfModel.computeShape(modelAvgPos).getDeterminantRadius()*sigma2fwhm
 
         basisList = makeKernelBasisList(self.kConfig, fwhmScience, fwhmModel, metadata=self.metadata)