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 all 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
44 changes: 41 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,26 @@ 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}'
)
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}')
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.

Why this? I don't think there is any restriction that either of these must be odd. Odd model size is common, but not really necessary. Most PSF profiles have a lot of flux in the central 1 pixel, so odd helps model that better. But there could certainly be cases where even works better. (I did test both odd and even for DES, and the odd sizes in fact worked better in those tests -- I haven't switched to even sizes since then.)
And I usually use even stampSize actually. The default in Piff is 32, which is much larger than our modelSize in DES, even without any rotator.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I agree that there is no reason in principle.

In practice it is nice to have odd for some applications. For example, odd means the model can be centered on that center pixel, and then convolutions with the PSF image will not cause a shift.

Copy link
Member

@arunkannawadi arunkannawadi 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.

Hmm, stampSize is enforced to be odd upstream as well, so even if we drop that requirement here, you can be guaranteed that stampSize will be odd. Dropping that requirement upstream will have to get discussed within the Science Pipelines team but I can go with whatever is the consensus is on modelSize. But based on Mike's comment, may be we should not make it a requirement just yet.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

We want it for shear

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.

In practice it is nice to have odd for some applications. For example, odd means the model can be centered on that center pixel, and then convolutions with the PSF image will not cause a shift.

This consideration is orthogonal to the reasons for having the model use an odd size. You can draw it onto an odd-sized image (and center it) regardless of what the underlying model uses.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

OK.

I thought that re-gridding the model from a center-near-pixel-boundary to center-at-pixel-center would introduce larger interpolation errors. If that's not the case then I withdraw my assertion

Copy link
Member

Choose a reason for hiding this comment

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

I'll be moving around the code a bit after merging, so I will drop the requirement then.

Copy link
Contributor

Choose a reason for hiding this comment

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

I still think people should tend to use odd modelSize in most cases for other reasons. Mostly that typical PSFs have a hot spot right at the center, which an odd model size is better at describing. But making it a requirement feels a little heavy handed. E.g. modeling very out of focus stars, which don't have that central peak, might work better with an even model size for some reason.


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)
86 changes: 75 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,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"""
Expand All @@ -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

Expand Down
Loading