DM-43490: Put correct size to PixelGrid for PSF estimation

python/lsst/meas/extensions/piff/piffPsfDeterminer.py

@@ -80,6 +80,10 @@ class PiffPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass):
         "e.g. 0.5 is equal to 2x oversampling",
         default=1,
     )
+    modelSize = pexConfig.Field[int](
+        doc="Internal model size for PIFF (typically odd, but not enforced)",
+        default=25,
+    )
     outlierNSigma = pexConfig.Field[float](
         doc="n sigma for chisq outlier rejection",
         default=4.0
@@ -137,7 +141,30 @@ def setDefaults(self):
         # i) aperture flux with 12 pixel radius can be compared to PSF flux.
         # ii) fake sources injected to match the 12 pixel aperture flux get
         #     measured correctly
+        # stampSize should also be at least sqrt(2)*modelSize/samplingSize so
+        # that rotated images have support in the model
+
         self.stampSize = 25
+        # Resize the stamp to accommodate the model, but only if necessary.
+        if self.useCoordinates == "sky":
+            self.stampSize = max(25, 2*int(0.5*self.modelSize*np.sqrt(2)/self.samplingSize) + 1)
+
+    def validate(self):
+        super().validate()
+
+        if (stamp_size := self.stampSize) is not None:
+            model_size = self.modelSize
+            sampling_size = self.samplingSize
+            if self.useCoordinates == 'sky':
+                min_stamp_size = int(np.sqrt(2) * model_size / sampling_size)
+            else:
+                min_stamp_size = int(model_size / sampling_size)
+
+            if stamp_size < min_stamp_size:
+                msg = (f"PIFF model size of {model_size} is too large for stamp size {stamp_size}. "
+                       f"Set stampSize >= {min_stamp_size}"
+                       )
+                raise pexConfig.FieldValidationError(self.__class__.modelSize, self, msg)
 
 
 def getGoodPixels(maskedImage, zeroWeightMaskBits):
@@ -331,12 +358,14 @@ def determinePsf(
             stampSize = psfCandidateList[0].getWidth()
 
         scale = exposure.getWcs().getPixelScale().asArcseconds()
+
         match self.config.useCoordinates:
             case 'field':
                 detector = exposure.getDetector()
                 pix_to_field = detector.getTransform(PIXELS, FIELD_ANGLE)
                 gswcs = UVWcsWrapper(pix_to_field)
                 pointing = None
+
             case 'sky':
                 gswcs = CelestialWcsWrapper(exposure.getWcs())
                 skyOrigin = exposure.getWcs().getSkyOrigin()
@@ -346,6 +375,7 @@ def determinePsf(
                     ra*galsim.degrees,
                     dec*galsim.degrees
                 )
+
             case 'pixel':
                 gswcs = galsim.PixelScale(scale)
                 pointing = None
@@ -385,7 +415,7 @@ def determinePsf(
             'model': {
                 'type': 'PixelGrid',
                 'scale': scale * self.config.samplingSize,
-                'size': stampSize,
+                'size': self.config.modelSize,
                 'interp': self.config.interpolant
             },
             'interp': {

tests/test_psf.py

@@ -38,6 +38,7 @@
 from lsst.meas.base import SingleFrameMeasurementTask
 from lsst.meas.extensions.piff.piffPsfDeterminer import PiffPsfDeterminerConfig, PiffPsfDeterminerTask
 from lsst.meas.extensions.piff.piffPsfDeterminer import _validateGalsimInterpolant
+from packaging import version
 
 
 def psfVal(ix, iy, x, y, sigma1, sigma2, b):
@@ -55,6 +56,41 @@ def psfVal(ix, iy, x, y, sigma1, sigma2, b):
             + b*np.exp(-0.5*(u**2 + (v*ab)**2)/sigma2**2))/(1 + b)
 
 
+def make_wcs(angle_degrees=None):
+    """Make a simple SkyWcs that is rotated around the origin.
+
+    Parameters
+    ----------
+    angle_degrees : `float`, optional
+        The angle to rotate the WCS by, in degrees.
+
+    Returns
+    -------
+    wcs : `~lsst.afw.geom.SkyWcs`
+        The WCS object.
+    """
+    cdMatrix = np.array([
+        [1.0, 0.0],
+        [0.0, 1.0]
+    ]) * 0.2 / 3600
+
+    if angle_degrees is not None:
+        angle_radians = np.radians(angle_degrees)
+        cosang = np.cos(angle_radians)
+        sinang = np.sin(angle_radians)
+        rot = np.array([
+            [cosang, -sinang],
+            [sinang, cosang]
+        ])
+        cdMatrix = np.dot(cdMatrix, rot)
+
+    return afwGeom.makeSkyWcs(
+        crpix=geom.PointD(0, 0),
+        crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
+        cdMatrix=cdMatrix,
+    )
+
+
 class SpatialModelPsfTestCase(lsst.utils.tests.TestCase):
     """A test case for SpatialModelPsf"""
 
@@ -101,11 +137,7 @@ def setUp(self):
         self.exposure = afwImage.makeExposure(self.mi)
         self.exposure.setPsf(measAlg.DoubleGaussianPsf(self.ksize, self.ksize,
                                                        1.5*sigma1, 1, 0.1))
-        cdMatrix = np.array([1.0, 0.0, 0.0, 1.0]) * 0.2/3600
-        cdMatrix.shape = (2, 2)
-        wcs = afwGeom.makeSkyWcs(crpix=geom.PointD(0, 0),
-                                 crval=geom.SpherePoint(0.0, 0.0, geom.degrees),
-                                 cdMatrix=cdMatrix)
+        wcs = make_wcs()
         self.exposure.setWcs(wcs)
 
         #
@@ -193,6 +225,8 @@ def setUp(self):
     def setupDeterminer(
         self,
         stampSize=None,
+        kernelSize=None,
+        modelSize=25,
         debugStarData=False,
         useCoordinates='pixel',
         downsample=False,
@@ -204,6 +238,19 @@ def setupDeterminer(
         ----------
         stampSize : `int`, optional
             Set ``config.stampSize`` to this, if not None.
+        kernelSize : `int`, optional
+            Cutout size for the PSF candidates. This is unused if ``stampSize``
+            if provided and its value is used for cutout size instead.
+        modelSize : `int`, optional
+            Internal model size for PIFF.
+        debugStarData : `bool`, optional
+            Include star images used for fitting in PSF model object?
+        useCoordinates : `str`, optional
+            Spatial coordinates to regress against for PSF modelling.
+        downsample : `bool`, optional
+            Whether to downsample the PSF candidates before modelling?
+        withlog : `bool`, optional
+            Should Piff produce chatty log messages?
         """
         starSelectorClass = measAlg.sourceSelectorRegistry["objectSize"]
         starSelectorConfig = starSelectorClass.ConfigClass()
@@ -220,14 +267,18 @@ def setupDeterminer(
         self.starSelector = starSelectorClass(config=starSelectorConfig)
 
         makePsfCandidatesConfig = measAlg.MakePsfCandidatesTask.ConfigClass()
+        if kernelSize:
+            makePsfCandidatesConfig.kernelSize = kernelSize
         if stampSize is not None:
             makePsfCandidatesConfig.kernelSize = stampSize
+
         self.makePsfCandidates = measAlg.MakePsfCandidatesTask(config=makePsfCandidatesConfig)
 
         psfDeterminerConfig = PiffPsfDeterminerConfig()
         psfDeterminerConfig.spatialOrder = 1
-        if stampSize is not None:
-            psfDeterminerConfig.stampSize = stampSize
+        psfDeterminerConfig.stampSize = stampSize
+        psfDeterminerConfig.modelSize = modelSize
+
         psfDeterminerConfig.debugStarData = debugStarData
         psfDeterminerConfig.useCoordinates = useCoordinates
         if downsample:
@@ -237,7 +288,7 @@ def setupDeterminer(
 
         self.psfDeterminer = PiffPsfDeterminerTask(psfDeterminerConfig)
 
-    def subtractStars(self, exposure, catalog, chi_lim=-1):
+    def subtractStars(self, exposure, catalog, chi_lim=-1.):
         """Subtract the exposure's PSF from all the sources in catalog"""
         mi, psf = exposure.getMaskedImage(), exposure.getPsf()
 
@@ -307,7 +358,7 @@ def checkPiffDeterminer(self, **kwargs):
             chiLim = 7.0
         else:
             numAvail = len(psfCandidateList)
-            chiLim = 6.1
+            chiLim = 6.4
 
         self.assertEqual(metadata['numAvailStars'], numAvail)
         self.assertEqual(sum(self.catalog['use_psf']), metadata['numGoodStars'])
@@ -385,10 +436,6 @@ def testPiffDeterminer_debugStarData(self):
         """Test Piff with debugStarData=True."""
         self.checkPiffDeterminer(debugStarData=True)
 
-    def testPiffDeterminer_skyCoords(self):
-        """Test Piff sky coords."""
-        self.checkPiffDeterminer(useCoordinates='sky')
-
     def testPiffDeterminer_downsample(self):
         """Test Piff determiner with downsampling."""
         self.checkPiffDeterminer(downsample=True)
@@ -397,6 +444,47 @@ def testPiffDeterminer_withlog(self):
         """Test Piff determiner with chatty logs."""
         self.checkPiffDeterminer(withlog=True)
 
+    def testPiffDeterminer_stampSize26(self):
+        """Test Piff with a psf stampSize of 26."""
+        with self.assertRaises(ValueError):
+            self.checkPiffDeterminer(stampSize=26)
+
+    def testPiffDeterminer_modelSize26(self):
+        """Test Piff with a psf stampSize of 26."""
+        with self.assertRaises(ValueError):
+            self.checkPiffDeterminer(modelSize=26, stampSize=25)
+
+    def testPiffDeterminer_skyCoords(self):
+        """Test Piff sky coords."""
+
+        self.checkPiffDeterminer(useCoordinates='sky')
+
+    @lsst.utils.tests.methodParameters(angle_degrees=[0, 35, 77, 135])
+    def testPiffDeterminer_skyCoords_with_rotation(self, angle_degrees):
+        """Test Piff sky coords with rotation."""
+
+        wcs = make_wcs(angle_degrees=angle_degrees)
+        self.exposure.setWcs(wcs)
+        self.checkPiffDeterminer(useCoordinates='sky', kernelSize=35)
+
+    # TODO: Merge this case with the above in DM-44467.
+    @unittest.skipUnless(version.parse(galsim.version) >= version.parse("2.5.2"),
+                         reason="Requires GalSim >= 2.5.2",
+                         )
+    def testPiffDeterminer_skyCoords_rot45(self):
+        """Test Piff sky coords."""
+
+        wcs = make_wcs(angle_degrees=45)
+        self.exposure.setWcs(wcs)
+        self.checkPiffDeterminer(useCoordinates='sky', kernelSize=35)
+
+    def testPiffDeterminer_skyCoords_failure(self, angle_degrees=135):
+        """Test that using small PSF candidates with sky coordinates fails."""
+        wcs = make_wcs(angle_degrees=angle_degrees)
+        self.exposure.setWcs(wcs)
+        with self.assertRaises(ValueError):
+            self.checkPiffDeterminer(useCoordinates='sky', stampSize=15)
+
 
 class PiffConfigTestCase(lsst.utils.tests.TestCase):
     """A test case to check for valid Piff config"""