DM-32411: Make WarpedPsf delegate to source PSF's computeImage.

[TallJimbo]

No description provided.

[jmeyers314]

I've been playing around with the following script to see how this PR is working:

import lsst.afw.detection as afwDetection
import lsst.afw.image as afwImage
import lsst.afw.math as afwMath
import lsst.afw.geom as afwGeom
import lsst.geom as geom
import numpy as np
from lsst.meas.algorithms import WarpedPsf


class PyGaussianPsf(afwDetection.Psf):
    # Like afwDetection.GaussianPsf, but handles computeImage exactly instead of
    # via interpolation.  This is a subminimal implementation.  It works for the
    # tests here but isn't fully functional as a Psf class.

    def __init__(self, width, height, sigma):
        afwDetection.Psf.__init__(self, isFixed=False)
        self.dimensions = geom.Extent2I(width, height)
        self.sigma = sigma

    def _doComputeKernelImage(self, position=None, color=None):
        bbox = self.computeBBox(position, color)
        img = afwImage.Image(bbox, dtype=np.float64)
        x, y = np.ogrid[bbox.minY:bbox.maxY+1, bbox.minX:bbox.maxX+1]
        rsqr = x**2 + y**2
        img.array[:] = np.exp(-0.5*rsqr/self.sigma**2)
        img.array /= np.sum(img.array)
        return img

    def _doComputeImage(self, position=None, color=None):
        bbox = self.computeBBox(position, color)
        img = afwImage.Image(bbox, dtype=np.float64)
        y, x = np.ogrid[float(bbox.minY):bbox.maxY+1, bbox.minX:bbox.maxX+1]
        x -= (position.x - np.floor(position.x+0.5))
        y -= (position.y - np.floor(position.y+0.5))
        rsqr = x**2 + y**2
        img.array[:] = np.exp(-0.5*rsqr/self.sigma**2)
        img.array /= np.sum(img.array)
        img.setXY0(geom.Point2I(
            img.getX0() + np.floor(position.x+0.5),
            img.getY0() + np.floor(position.y+0.5)
        ))
        return img

    def _doComputeBBox(self, position=None, color=None):
        dims = self.dimensions
        return geom.Box2I(geom.Point2I(-dims/2), dims)

    def _doComputeShape(self, position=None, color=None):
        return afwGeom.ellipses.Quadrupole(self.sigma**2, self.sigma**2, 0.0)

psf = PyGaussianPsf(15, 15, 2.0)
# psf = afwDetection.GaussianPsf(15, 15, 2.0)

transform = afwGeom.makeTransform(geom.AffineTransform([0.1, 0.1]))
warpingControl = afwMath.WarpingControl("lanczos3")
wPsf = WarpedPsf(psf, transform, warpingControl)

trueImg = psf.computeImage(geom.Point2D(0.0, 0.5))
warpedImg = wPsf.computeImage(geom.Point2D(0.0, 0.5))

print(warpedImg.array.shape)

One thing I've noticed is that the WarpedPsf that emerges with this PR has a bbox 1 pixel larger than the pre-PR version. Aside from that, I haven't noticed any particular improvement yet, but maybe I'm not properly oriented.

[jmeyers314]

Thinking about this some more... for a given WarpedPsf evaluation, there are essentially 3 pixel grids in play. There's the real grid where the undistorted Psf draws with computeImage, there's the hypothetical grid where the undistorted Psf draws with computeKernelImage, and there's the real grid where WarpedPsf draws with computeImage. For a pure shift, that means it must be possible to arrange for the Psf.computeKernelImage grid to exactly overlap the WarpedPsf.computeImage grid, and provided the Psf is analytic, say, no interpolation would be required. The above snippet features a purely analytic Psf though, and I can't seem to find a setup that elides the interpolation.

[TallJimbo]

At some level, the motivation for this change is not so much to avoid interpolation as to ensure the interpolation of the real image is the same as the interpolation of the PSF model. So even if we had a Psf interface for native evaluation with some transform applied, we wouldn't necessarily want to use it in coadd PSF (which is really all WarpedPsf is for, at least at present).

the WarpedPsf that emerges with this PR has a bbox 1 pixel larger than the pre-PR version

Is this just for computeImage, or for computeKernelImage, too? I think we want to guarantee that the latter has odd size (though I don't care exactly what it is; that's just a question of rounding, I think).

[arunkannawadi]

Is this just for computeImage, or for computeKernelImage, too? I think we want to guarantee that the latter has odd size (though I don't care exactly what it is; that's just a question of rounding, I think).

For both computeImage and computeKernelImage, the bbox is coming out to be of shape (16, 16).

[TallJimbo]

Ok, thanks for the confirmation. I'll go get a fix on the branch ASAP.

[arunkannawadi]

Josh's script now returns (17, 17) whereas the main branch returns (15, 15). I edited geom::ceil in your recent commit to geom::floor and it returned (15, 15). Do you see any other failures if we used geom::floor?

[TallJimbo]

I don't anticipate anything concrete going wrong; I was leaning on the side of growing too much rather than growing too little, but there are effectively two places where we round up in that routine, so it's also not unreasonable to change that one to round down. Please go ahead and commit that change if it seems to make sizes stay more consistent.