diff --git a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py index 1c6f530..aefb427 100644 --- a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py +++ b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py @@ -78,6 +78,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", + default=25, + ) outlierNSigma = pexConfig.Field[float]( doc="n sigma for chisq outlier rejection", default=4.0 @@ -123,11 +127,17 @@ class PiffPsfDeterminerConfig(BasePsfDeterminerTask.ConfigClass): def setDefaults(self): super().setDefaults() - # stampSize should be at least 25 so that + # stampSize for PIFF means the size of input stamps extracted for each + # object. + # + # stampSize should be at least modelSize*sqrt(2)/samplingSize + # to ensure that rotated images have support in the model + # also we want to ensure # 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 - self.stampSize = 25 + # Default modelSize is 25, and 35 is about sqrt(2) * 25 + self.stampSize = 35 def getGoodPixels(maskedImage, zeroWeightMaskBits): @@ -291,6 +301,7 @@ def determinePsf( psfCellSet : `None` Unused by this PsfDeterminer. """ + if self.config.stampSize: stampSize = self.config.stampSize if stampSize > psfCandidateList[0].getWidth(): @@ -300,13 +311,17 @@ def determinePsf( self.log.debug("stampSize not set. Using candidate size.") stampSize = psfCandidateList[0].getWidth() + _check_stamp_size(self.config) + 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() @@ -316,6 +331,7 @@ def determinePsf( ra*galsim.degrees, dec*galsim.degrees ) + case 'pixel': gswcs = galsim.PixelScale(scale) pointing = None @@ -355,7 +371,7 @@ def determinePsf( 'model': { 'type': 'PixelGrid', 'scale': scale * self.config.samplingSize, - 'size': stampSize, + 'size': self.config.modelSize, 'interp': self.config.interpolant }, 'interp': { @@ -404,4 +420,26 @@ def determinePsf( return PiffPsf(drawSize, drawSize, piffResult), None +def _check_stamp_size(config): + coord_sys = config.useCoordinates + stamp_size = config.stampSize + model_size = config.modelSize + sampling_size = config.samplingSize + if coord_sys == '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: + raise ValueError( + 'for {coord_sys} coords and modelSize {model_size} ' + f'we require stampSize >= ' + f'{min_stamp_size}, got {stamp_size}' + ) + if stamp_size % 2 == 0: + raise ValueError(f'stampSize should be odd, got {stamp_size}') + if model_size % 2 == 0: + raise ValueError(f'modelSize should be odd, got {model_size}') + + measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask) diff --git a/tests/test_psf.py b/tests/test_psf.py index fdad1e1..8414a0a 100644 --- a/tests/test_psf.py +++ b/tests/test_psf.py @@ -92,11 +92,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 = _get_wcs() self.exposure.setWcs(wcs) # @@ -183,7 +179,8 @@ def setUp(self): def setupDeterminer( self, - stampSize=None, + stampSize=35, + modelSize=25, debugStarData=False, useCoordinates='pixel' ): @@ -209,14 +206,15 @@ def setupDeterminer( self.starSelector = starSelectorClass(config=starSelectorConfig) makePsfCandidatesConfig = measAlg.MakePsfCandidatesTask.ConfigClass() - if stampSize is not None: - makePsfCandidatesConfig.kernelSize = stampSize + 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 @@ -284,7 +282,8 @@ def checkPiffDeterminer(self, **kwargs): self.assertNotIn('image', psf._piffResult.stars[0].data.__dict__) # Test how well we can subtract the PSF model - self.subtractStars(self.exposure, self.catalog, chi_lim=6.1) + # note chi_lim increased from 6.1 to 6.5 + self.subtractStars(self.exposure, self.catalog, chi_lim=6.5) # Test bboxes for point in [ @@ -345,8 +344,50 @@ def testPiffDeterminer_debugStarData(self): def testPiffDeterminer_skyCoords(self): """Test Piff sky coords.""" + + self.checkPiffDeterminer(useCoordinates='sky') + + def testPiffDeterminer_skyCoords_rot35(self): + """Test Piff sky coords.""" + + wcs = _get_wcs(angle_degrees=35) + self.exposure.setWcs(wcs) + self.checkPiffDeterminer(useCoordinates='sky') + + def testPiffDeterminer_skyCoords_rot77(self): + """Test Piff sky coords.""" + + # causes a relatively bad fit + wcs = _get_wcs(angle_degrees=77) + self.exposure.setWcs(wcs) self.checkPiffDeterminer(useCoordinates='sky') + def testPiffDeterminer_skyCoords_rot135(self): + """Test Piff sky coords.""" + + wcs = _get_wcs(angle_degrees=135) + self.exposure.setWcs(wcs) + self.checkPiffDeterminer(useCoordinates='sky') + + # 45 degrees causes singular matrix in PIFF + # interestingly 135 does not + def testPiffDeterminer_skyCoords_rot45(self): + """Test Piff sky coords.""" + + wcs = _get_wcs(angle_degrees=45) + self.exposure.setWcs(wcs) + self.checkPiffDeterminer(useCoordinates='sky') + + 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) + class PiffConfigTestCase(lsst.utils.tests.TestCase): """A test case to check for valid Piff config""" @@ -371,6 +412,29 @@ def testValidateGalsimInterpolant(self): self.assertTrue(eval(f"galsim.{interp}")) +def _get_wcs(angle_degrees=None): + 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 TestMemory(lsst.utils.tests.MemoryTestCase): pass