Skip to content

DM-53465: Fix PSF normalization #171

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

Merged
arunkannawadi merged 3 commits into from
Dec 2, 2025
Merged

DM-53465: Fix PSF normalization #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
12dc6b8
Improve the log message when PSF is not normalized
arunkannawadi Dec 2, 2025
fca02cb
Calculate the padding and zero-pad PSF image
arunkannawadi Dec 2, 2025
386051a
Add a unit test checking PSF normalization
arunkannawadi Dec 2, 2025
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
28 changes: 24 additions & 4 deletions python/lsst/drp/tasks/assemble_cell_coadd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -323,6 +323,20 @@ def __init__(self, *args, **kwargs):
self.makeSubtask("scale_zero_point")

self.psf_warper = afwMath.Warper.fromConfig(self.config.psf_warper)
if (warping_kernel_name := self.config.psf_warper.warpingKernelName.lower()).startswith("lanczos"):
psf_padding = 2 * int(warping_kernel_name.lstrip("lanczos")) - 1
self.log.debug(
"Padding PSF image by %d pixels since the warping kernel is %s.",
psf_padding,
self.config.psf_warper.warpingKernelName,
)
else:
psf_padding = 10
Copy link
Contributor

@taranu taranu Dec 2, 2025

Choose a reason for hiding this comment

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

The only other non-lanczos option is bilinear; why 10 for that and only 5 for lanczos3?

Copy link
Member Author

@arunkannawadi arunkannawadi Dec 2, 2025

Choose a reason for hiding this comment

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

5 for Lanczos3 is well-motivated because that is the support of the kernel and verified that making it 6 or 7 made no difference to the output. I didn't want to make any strong assumptions about the support of non-Lanczos kernels and took a more conservative approach of 10. I also didn't want to introduce it as a config option.

self.log.info(
"Padding PSF image by %d pixels since the warping kernel is not Lanczos.",
psf_padding,
)
self.psf_padding = psf_padding

def runQuantum(self, butlerQC, inputRefs, outputRefs):
# Docstring inherited.
Expand Down Expand Up @@ -790,8 +804,13 @@ def run(
undistorted_psf_im.getDimensions(),
), "PSF image does not share the coordinates of the 'calexp'"

# Convert the PSF image from Image to MaskedImage.
undistorted_psf_maskedImage = afwImage.MaskedImageD(image=undistorted_psf_im)
# Convert the PSF image from Image to MaskedImage and
# zero-pad the image.
undistorted_psf_bbox = undistorted_psf_im.getBBox()
undistorted_psf_maskedImage = afwImage.MaskedImageD(
undistorted_psf_bbox.dilatedBy(self.psf_padding)
)
undistorted_psf_maskedImage.image[undistorted_psf_bbox].array[:, :] = undistorted_psf_im.array
# TODO: In DM-43585, use the variance plane value from noise.
undistorted_psf_maskedImage.variance += 1.0 # Set variance to 1

Expand All @@ -809,11 +828,12 @@ def run(
psf_stacker = psf_stacker_gc[cellInfo.index]
psf_stacker.add_masked_image(warped_psf_maskedImage, weight=weight)

if warped_psf_maskedImage.image.array.sum() < 0.995:
if not (0.995 < (psf_normalization := warped_psf_maskedImage.image.array.sum()) < 1.005):
Copy link
Contributor

@taranu taranu Dec 2, 2025

Choose a reason for hiding this comment

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

I think at some point this needs to actually do something rather than warn, like try to increase the padding and then fall back to renormalizing by the sum (which should only happen if the PSF cut off, i.e. >0 at the edge).

Copy link
Member Author

@arunkannawadi arunkannawadi Dec 2, 2025

Choose a reason for hiding this comment

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

Maybe, but for now, I just want to monitor this to make sure this never happens. Piff guarantees that the model it returns is normalized and lanczos resampling shouldn't introduce deviations that are larger than 0.005. If the PSF normalization is off by more than 0.5%, it means something with the assumptions made is incorrect and I'd like to know that before padding is auto-adjusted.

self.log.warning(
"PSF for %s in %s lost more than 0.5 per cent of flux",
"PSF image for %s in %s is not normalized to 1.0, but instead %f",
warp_input.dataId,
cellInfo.index,
psf_normalization,
)

if (ap_corr_map := warp.getInfo().getApCorrMap()) is not None:
Expand Down
7 changes: 7 additions & 0 deletions tests/test_assemble_cell_coadd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,13 @@ def test_inputs_sorted(self):
for _, singleCellCoadd in self.result.multipleCellCoadd.cells.items():
self.checkSortOrder(singleCellCoadd.inputs)

def test_psf_normalization(self):
"""Check that the sum of PSF images is close to 1."""
self.runTask()
for cellId, singleCellCoadd in self.result.multipleCellCoadd.cells.items():
with self.subTest(x=repr(cellId.x), y=repr(cellId.y)):
self.assertFloatsAlmostEqual(singleCellCoadd.psf_image.array.sum(), 1.0, rtol=None, atol=1e-7)


class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
pass
Expand Down