Add option to have InterpolatedImage de-pixelize the image

[rmjarvis]

One of the persistent issues with using PSF models such as those from PSFEx or Piff is that they generally model the "effective PSF", including the pixel convolution, rather than the underlying true PSF. This means that any objects that use this PSF object have to be rendered with method='no_pixel', rather than 'auto', or especially 'phot'. Using photon shooting would implicitly convolve a second pixel, so lead to significant biases in the resulting images.

This PR adds the capability of fitting for a "depixelized" image to input into InterpolatedImage, which would represent the underlying profile, without the pixel convolution. The expected user API for this is simply to add the parameter depixelize=True to the InterpolatedImage constructor.

Then, to use a Piff model, for instance, with photon shooting, one can draw it normally, and then model it as InterpolatedImage(piff_image, depixelize=True). Then convolve this by the galaxy profile and you can render the combination using method='phot'.

It's a big matrix equation, so it scales badly for large images, but for the kinds of images we usually have for PSFs, it takes around 0.1-0.2 seconds, which doesn't seem to bad for many use cases. Especially if you're in a testing regime where you can treat the PSF as constant for a reasonable number of galaxies.

[rmjarvis]

Yes, that's a very good point about the caching possibility. I'll have to think about how best to implement the cache, but that definitely seems worth doing.

[rmjarvis]

Hi Matt, I think I have the cache working. Thanks for the idea!

Some timings:

• 32x32:
  • First time: 0.013 s
  • Second time: 0.0007 s
• 48x48:
  • First time: 0.13 s
  • Second time: 0.0019 s
• 51x51:
  • First time: 0.18 s
  • Second time: 0.0034 s
• 64x64:
  • First time: 0.70 s
  • Second time: 0.0045 s
• 80x80:
  • First time: 2.5 s
  • Second time: 0.010 s
• 96x96:
  • First time: 7.3 s
  • Second time: 0.021 s
• 128x128:
  • First time: 42 s
  • Second time: 0.051 s

So this is a huge improvement. Even 128x128 is pretty acceptable after the initial pass. The first pass time looks close to the expected scaling of O(Npix^3), and subsequent ones are only O(Npix^2).

[arunkannawadi]

This seems like a useful feature to have. Thanks Mike and Matt!

[rmjarvis]

Thanks Matt. I'll let this sit another day or two before merging in case anyone else wants to have a look.

[jmeyers314]

Super cool. My only comment is that I think it may be helpful if the InterpolatedImage docstring had more guidance about setting x_interpolant and maxk when using this feature (similar to comments you made in the test suite). I see some moderately significant swings depending on these variables:

import numpy as np
import galsim
import matplotlib.pyplot as plt

def tryit(fwhm=0.6, maxk=50, interp=galsim.Lanczos(11)):    
    psf = galsim.Kolmogorov(fwhm=fwhm).shift(0.01, 0.04).shear(g1=0.1, g2=-0.02)
    psfimg_real = psf.drawImage(scale=0.2, nx=64, ny=64, dtype=float, method='real_space')
    psfimg_pixel = psf.drawImage(scale=0.2, nx=64, ny=64, dtype=float)
    psfimg_nopixel = psf.drawImage(scale=0.2, nx=64, ny=64, method='no_pixel', dtype=float)
    # normalize
    psfimg_pixel /= np.sum(psfimg_pixel.array)
    psfimg_nopixel /= np.sum(psfimg_nopixel.array)
    # reconstruct psf from psfimg_pixel
    ii = galsim.InterpolatedImage(psfimg_pixel, depixelize=True, _force_maxk=maxk, x_interpolant=interp)
    # draw, implicitly convolving with pixel
    iiimg = ii.drawImage(scale=0.2, nx=64, ny=64, dtype=float)
    return psfimg_real, psfimg_pixel, psfimg_nopixel, iiimg

fig, axes = plt.subplots(ncols=2, nrows=1, figsize=(8, 4))
fwhms = np.linspace(0.3, 2.0, 50)
for interp, maxk in [
    (galsim.Quintic(), 0.0), 
    (galsim.Quintic(), 50),
    (galsim.Lanczos(5), 0.0), 
    (galsim.Lanczos(5), 50), 
    (galsim.Lanczos(11), 50)
]:
    max_resids = []
    T_resids = []
    for fwhm in fwhms:
        psfimg_real, psfimg_pixel, psfimg_nopixel, iiimg = tryit(fwhm, interp=interp, maxk=maxk)        
        mom_psf = galsim.hsm.FindAdaptiveMom(psfimg_pixel)
        mom_ii = galsim.hsm.FindAdaptiveMom(iiimg)
        max_resids.append(np.max(np.abs(iiimg.array-psfimg_pixel.array)))
        T_resids.append((mom_ii.moments_sigma - mom_psf.moments_sigma)/mom_psf.moments_sigma)
    axes[0].plot(fwhms, max_resids, label=str(interp)+','+str(maxk))
    axes[1].plot(fwhms, T_resids, label=str(interp)+','+str(maxk))
axes[0].legend()
axes[1].legend()
axes[0].set_xlabel("FWHM (arcsec)")
axes[1].set_xlabel("FWHM (arcsec)")
axes[0].set_ylabel("max(|psf_eff - psf_reconstructed|)")
axes[1].set_ylabel("dT/T")
axes[0].set_yscale("log")
axes[1].set_yscale("log")
plt.tight_layout()
plt.show()

gives me

[rmjarvis]

@jmeyers314 I added "We recommend using a Lanczos interpolant with this option for best results. (Higher order tends to work better here.)" to the discussion of the depixelize option. I think the maxk issue I identified is independent of depixelize and is instead a problem intrinsic to how we calculate maxk for InterpolatedImages. So instead of adding something there, I added documentation for the _force_stepk and _force_maxk parameters:

        _force_stepk:       Override the normal stepk calculation (using gsparams.folding_threshold)
                            and force stepk to the given value. [default: 0]
        _force_maxk:        Override the normal maxk calculation (using gsparams.maxk_threshold)
                            and force maxk to the given value.  This option in particular can help
                            reduce FFT artifacts in a manner that is currently unobtainable by
                            lowering maxk_threshold. [default: 0]