diff --git a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py index 4670968..239b402 100644 --- a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py +++ b/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': { diff --git a/tests/test_psf.py b/tests/test_psf.py index 096023b..332f904 100644 --- a/tests/test_psf.py +++ b/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"""