Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add modelSize to piff config #23

Merged
merged 5 commits into from
Apr 12, 2024
Merged

Add modelSize to piff config #23

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e051711
add modelSize to piff config
esheldon Apr 4, 2024
7fe21e9
update tests to include rotation
esheldon Apr 4, 2024
6d16110
use property list to convert to galsim wcs
esheldon Apr 5, 2024
cb3fe08
revert to using CelestialWcsWrapper
esheldon Apr 5, 2024
5529420
add check that stampSize/modelSize are odd
esheldon Apr 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 37 additions & 3 deletions python/lsst/meas/extensions/piff/piffPsfDeterminer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
esheldon marked this conversation as resolved.
Show resolved Hide resolved
)
outlierNSigma = pexConfig.Field[float](
doc="n sigma for chisq outlier rejection",
default=4.0
Expand Down Expand Up @@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that in the default processing we don't use sky coordinates, I'm inclined to keep this as 25 in the default here. For your testing purposes, you should be able to override this value.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rmjarvis Is there any disadvantage to having modelSize == stampSize?

Copy link
Contributor

@rmjarvis rmjarvis Apr 12, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. You definitely don't want stampSize == modelSize. It should be larger for 3 reasons.

  1. The Lanczos interpolation extends the active modeling space out from the modelSize by interpolant.xrange(). I think your default is Lanczos(11), so it will reach out by 11 pixels on all sides. That alone suggests you want stampSize = modelSize + 22.
  2. Even if your centroids are absolutely perfect, there is a sub-pixel variation of half a pixel in some direction. So you would want to add one more pixel of data around the edges to make sure you are properly constraining all the pixels around the edge. That adds another +2 to the stamp size. Maybe a little more if you don't trust your centroids to be very sub-pixel accurate.
  3. Finally, the rotation. This multiplies that calculation by up to sqrt(2). (1/max(|sin(theta)|,|cos(theta)|) if you want to calculate it explicitly for a given WCS.)

And honestly, there is barely any down side to using a larger than needed stamp size. Piff will just ignore any pixels that are not being used to constrain the model. So I'd recommend the default being much larger in fact. Use 80x80 or something.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The default I put in was 35 which was designed to account for your part 3. in the list, since the default modelSize was set to 25. It seems like 1. and 2. imply we should go even larger

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The Lanczos interpolation extends the active modeling space out from the modelSize by interpolant.xrange()

Mike, could you confirm that this is needed even if the modeling is happening in pixel coordinates as opposed to sky coordinates?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"Needed"? No. Not technically. But it will produce inaccurate PSF models near the edges if you don't.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, I remember the compute scales as cubic power of the size, but does it as modelSize**3 or stampSize**3? I hope the former.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

modelSize**6 actually.



def getGoodPixels(maskedImage, zeroWeightMaskBits):
Expand Down Expand Up @@ -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():
Expand All @@ -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()
Expand All @@ -316,6 +331,7 @@ def determinePsf(
ra*galsim.degrees,
dec*galsim.degrees
)

case 'pixel':
gswcs = galsim.PixelScale(scale)
pointing = None
Expand Down Expand Up @@ -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': {
Expand Down Expand Up @@ -404,4 +420,22 @@ def determinePsf(
return PiffPsf(drawSize, drawSize, piffResult), None


def _check_stamp_size(config):
arunkannawadi marked this conversation as resolved.
Show resolved Hide resolved
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}'
)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note to self: this is a check only on the config, so should be implemented in a validate method.


measAlg.psfDeterminerRegistry.register("piff", PiffPsfDeterminerTask)
76 changes: 65 additions & 11 deletions tests/test_psf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

#
Expand Down Expand Up @@ -183,7 +179,8 @@ def setUp(self):

def setupDeterminer(
self,
stampSize=None,
stampSize=35,
modelSize=25,
debugStarData=False,
useCoordinates='pixel'
):
Expand All @@ -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
arunkannawadi marked this conversation as resolved.
Show resolved Hide resolved

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

Expand Down Expand Up @@ -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 [
Expand Down Expand Up @@ -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"""
Expand All @@ -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

Expand Down