From e051711f916958445cb3d4adbc6cf782c4aaaf26 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Thu, 4 Apr 2024 13:57:08 -0400 Subject: [PATCH 1/5] add modelSize to piff config The modelSize is the internal PIFF model grid size. This gets sent to the piff config rather than stampSize for the "size" entry. modelSize defaults to 25. Update the stampSize to default to sqrt(2) * 25 ~ 35 to accomodate rotations Add a check that stamp size is big enough --- .../meas/extensions/piff/piffPsfDeterminer.py | 40 +++++++++++++++++-- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py index 1c6f530..47b9c1d 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,22 @@ 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}' + ) + + measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask) From 7fe21e99d6a2aa0fe46aef1c8ff2b869230e2b3a Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Thu, 4 Apr 2024 13:58:40 -0400 Subject: [PATCH 2/5] update tests to include rotation This is an important check when working in sky coords Oddly at 45 degrees piff is giving a singlular matrix error and fails to fit. We should track this down before accepting this PR. 135 works fine so I wonder if there is something sepcial about 45 for this test --- tests/test_psf.py | 76 ++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 65 insertions(+), 11 deletions(-) diff --git a/tests/test_psf.py b/tests/test_psf.py index fdad1e1..064df35 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,40 @@ 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') + class PiffConfigTestCase(lsst.utils.tests.TestCase): """A test case to check for valid Piff config""" @@ -371,6 +402,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 From 6d161106140849297fbaf6c379d1c76c4683a4f6 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 5 Apr 2024 14:20:28 -0400 Subject: [PATCH 3/5] use property list to convert to galsim wcs Will work for FITS wcs --- .../meas/extensions/piff/piffPsfDeterminer.py | 19 +++++++++---------- tests/test_psf.py | 12 ++++++------ 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py index 47b9c1d..7c0c266 100644 --- a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py +++ b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py @@ -31,7 +31,7 @@ import lsst.meas.algorithms as measAlg from lsst.meas.algorithms.psfDeterminer import BasePsfDeterminerTask from .piffPsf import PiffPsf -from .wcs_wrapper import CelestialWcsWrapper, UVWcsWrapper +from .wcs_wrapper import UVWcsWrapper def _validateGalsimInterpolant(name: str) -> bool: @@ -270,6 +270,12 @@ def _computeWeightAlternative(maskedImage, maxSNR): return weightArr +def _make_gs_wcs(dm_wcs): + plist = dm_wcs.getFitsMetadata() + wcs = galsim.FitsWCS(header=plist) + return wcs + + class PiffPsfDeterminerTask(BasePsfDeterminerTask): """A measurePsfTask PSF estimator using Piff as the implementation. """ @@ -323,15 +329,8 @@ def determinePsf( pointing = None case 'sky': - gswcs = CelestialWcsWrapper(exposure.getWcs()) - skyOrigin = exposure.getWcs().getSkyOrigin() - ra = skyOrigin.getLongitude().asDegrees() - dec = skyOrigin.getLatitude().asDegrees() - pointing = galsim.CelestialCoord( - ra*galsim.degrees, - dec*galsim.degrees - ) - + gswcs = _make_gs_wcs(exposure.getWcs()) + pointing = gswcs.toWorld(gswcs.origin) case 'pixel': gswcs = galsim.PixelScale(scale) pointing = None diff --git a/tests/test_psf.py b/tests/test_psf.py index 064df35..bd10306 100644 --- a/tests/test_psf.py +++ b/tests/test_psf.py @@ -371,12 +371,12 @@ def testPiffDeterminer_skyCoords_rot135(self): # 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_skyCoords_rot45(self): + """Test Piff sky coords.""" + + wcs = _get_wcs(angle_degrees=45) + self.exposure.setWcs(wcs) + self.checkPiffDeterminer(useCoordinates='sky') class PiffConfigTestCase(lsst.utils.tests.TestCase): From cb3fe089b0c8208cc9ea0a3fea029d9125682b27 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Fri, 5 Apr 2024 14:44:12 -0400 Subject: [PATCH 4/5] revert to using CelestialWcsWrapper We think its an issue with the underlying galsim code rather than the wrapper --- .../meas/extensions/piff/piffPsfDeterminer.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py index 7c0c266..47b9c1d 100644 --- a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py +++ b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py @@ -31,7 +31,7 @@ import lsst.meas.algorithms as measAlg from lsst.meas.algorithms.psfDeterminer import BasePsfDeterminerTask from .piffPsf import PiffPsf -from .wcs_wrapper import UVWcsWrapper +from .wcs_wrapper import CelestialWcsWrapper, UVWcsWrapper def _validateGalsimInterpolant(name: str) -> bool: @@ -270,12 +270,6 @@ def _computeWeightAlternative(maskedImage, maxSNR): return weightArr -def _make_gs_wcs(dm_wcs): - plist = dm_wcs.getFitsMetadata() - wcs = galsim.FitsWCS(header=plist) - return wcs - - class PiffPsfDeterminerTask(BasePsfDeterminerTask): """A measurePsfTask PSF estimator using Piff as the implementation. """ @@ -329,8 +323,15 @@ def determinePsf( pointing = None case 'sky': - gswcs = _make_gs_wcs(exposure.getWcs()) - pointing = gswcs.toWorld(gswcs.origin) + gswcs = CelestialWcsWrapper(exposure.getWcs()) + skyOrigin = exposure.getWcs().getSkyOrigin() + ra = skyOrigin.getLongitude().asDegrees() + dec = skyOrigin.getLatitude().asDegrees() + pointing = galsim.CelestialCoord( + ra*galsim.degrees, + dec*galsim.degrees + ) + case 'pixel': gswcs = galsim.PixelScale(scale) pointing = None From 5529420d5f7ef2326464a47860d335264b9381b2 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Thu, 11 Apr 2024 11:23:18 -0400 Subject: [PATCH 5/5] add check that stampSize/modelSize are odd --- python/lsst/meas/extensions/piff/piffPsfDeterminer.py | 4 ++++ tests/test_psf.py | 10 ++++++++++ 2 files changed, 14 insertions(+) diff --git a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py index 47b9c1d..aefb427 100644 --- a/python/lsst/meas/extensions/piff/piffPsfDeterminer.py +++ b/python/lsst/meas/extensions/piff/piffPsfDeterminer.py @@ -436,6 +436,10 @@ def _check_stamp_size(config): 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 bd10306..8414a0a 100644 --- a/tests/test_psf.py +++ b/tests/test_psf.py @@ -378,6 +378,16 @@ def testPiffDeterminer_skyCoords_rot45(self): 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"""