How do you get a properly normalized rendering of an undersampled ImagePSF?

[rknop]

My situation is I have a PSF that is undersampled -- for example, a gaussian with 1σ width of 0.42 detector pixels, and a total flux of 1.

I create a ImagePSF that is oversampled by a factor of (say) 5.

If I run this through a photutils.psf.PSFPhotometry fitter, it all works well.

However, if I evaluate the ImagePSF (i.e. use either the __call__ method or the evaluate method), the sum of the returned image stamp is not 1.0, but something larger. Because the fit flux from PSFPhotometry comes out right, whatever the fitter is doing exactly compensates for (or, more likely, implicitly assumes the same interpolation) the non-normalization that results on the image from evaluating the PSF.

If I do this what a PSF that's not undersampled (e.g. if the Gaussian 1σ width is 3 detector pixels), then the evaluated image does come out summing to 1.0.

My question is: is there a method I can call to evaluate PSF that is undersampled on the image but stored oversampled, and get back an image clip that is normalized right (in my case, sums to 1, assuming a big enough clip size)?

Here is example code that shows what I'm talking about>

import numpy as np

import astropy.table
import photutils.psf


# You probably want stampsize = 2 * 5 * 2√(2ln2) * max(sigmax,sigmay)
def make_gaussian_psf( oversamp=3, sigmax=1.2, sigmay=None, stampsize=None, x0=511, y0=511 ):
    if oversamp % 2 != 1:
        raise ValueError( "Please use an odd oversampling factor." )

    sigmay = sigmax if sigmay is None else sigmay

    # Default to a stamp size of 2 * 5 FHWMs (so "radius" 5 FWHM)
    if stampsize is None:
        stampsize = int( np.ceil( 10 * 2.35482 * max( sigmax, sigmay ) ) )
        stampsize += 1 if stampsize % 2 ==0 else 0
    if stampsize % 2 != 1:
        raise ValueError( "stampsize must be odd" )

    oversampsize = oversamp * stampsize
    ctr = oversampsize // 2
    xvals = np.arange( -ctr, ctr+1 )
    yvals = np.arange( -ctr, ctr+1 )
    ovsigmax = oversamp * sigmax
    ovsigmay = oversamp * sigmay
    data = np.exp( -( xvals[np.newaxis,:]**2 / ( 2. * ovsigmax**2 ) +
                      yvals[:,np.newaxis]**2 / ( 2. * ovsigmay**2 ) ) )
    # photutils expects an oversampled PSF to be normalized like this:
    data /= data.sum() / ( oversamp**2 )

    psf = photutils.psf.ImagePSF( data, flux=1, x_0=x0, y_0=y0, oversampling=oversamp )

    return psf


def make_fake_image( skynoise=10., starflux=100000., nx=1024, ny=1024, x=511, y=511,
                     oversamp=3, sigmax=1.2, sigmay=None, stampsize=None ):
    rng = np.random.default_rng()
    image = rng.normal( scale=skynoise, size=(ny, nx) )

    xc = int( np.floor( x + 0.5 ) )
    yc = int( np.floor( y + 0.5 ) )

    psf = make_gaussian_psf( oversamp=oversamp, sigmax=sigmax, sigmay=sigmay, stampsize=stampsize, x0=x, y0=y )
    xmin = -( stampsize // 2 ) + xc
    xmax = xmin + stampsize
    ymin = -( stampsize // 2 ) + yc
    ymax = ymin + stampsize
    # WORRY : did I get the x, y order right?
    xvals, yvals = np.meshgrid( np.arange( xmin, xmax ), np.arange( ymin, ymax ) )

    # Not adding poisson noise to our PSF; deal
    # Using the astropy FittableModel interface
    image[ ymin:ymax, xmin:xmax ] += psf( xvals, yvals ) * starflux
    # Directly calling ImagePSF.evaluate yields exactly the same result (no surprise)
    # image[ ymin:ymax, xmin:xmax ] += psf.evaluate( xvals, yvals, starflux, x, y )
    
    return image, psf


# **********************************************************************
# Actual experiments

# Have a plenty-sampled PSF with a sigma of 1.2, stampsize 29, oversampling by 3)
psf = make_gaussian_psf( oversamp=3, sigmax=1.2, stampsize=29, x0=511, y0=511 )
xvals = np.arange( 511 - (29//2), 511 + (29//2) + 1 )
yvals = np.arange( 511 - (29//2), 511 + (29//2) + 1 )
xvals, yvals = np.meshgrid( xvals, yvals )
print( f"1.2-sigma PSF clip sums to {psf( xvals, yvals ).sum()}" )

# Have a PSF with a 1-pixel FWHM (so σ = 1/(2√(2ln2)) = 0.42 ), stampsize 11, oversamped 5x
psf = make_gaussian_psf( oversamp=5, sigmax=0.42, stampsize=11, x0=511, y0=511 )
xvals = np.arange( 511 - (29//2), 511 + (29//2) + 1 )
yvals = np.arange( 511 - (29//2), 511 + (29//2) + 1 )
xvals, yvals = np.meshgrid( xvals, yvals )
print( f"0.42-sigma PSF clip sums to {psf( xvals, yvals ).sum()}" )


# Now let's try some PSF photometry and see how it does
# We are always going to use sky noise 10, and just not bother
#  with star poisson noise for now (because we're naughty)

# First with the happy big star
image, psf = make_fake_image( oversamp=3, sigmax=1.2, stampsize=29 )
noiseim = np.full_like( image, 10. )
fit_shape = ( 15, 15 )
photor = photutils.psf.PSFPhotometry( psf, fit_shape )
phot = photor( image, error=noiseim,
               init_params=astropy.table.Table( { 'x_init': [510.5], 'y_init': [513.], 'flux_init': [95000] } ) )
print( f"Fit flux of 100000-flux star for sigma = 1.2: {phot['flux_fit'][0]}; "
       f"(x,y)=({phot['x_fit'][0]}, {phot['y_fit'][0]}" )

# Now with a little star
image, psf = make_fake_image( oversamp=5, sigmax=0.42, stampsize=11 )
fit_shape = ( 7, 7 )
photor = photutils.psf.PSFPhotometry( psf, fit_shape )
phot = photor( image, error=noiseim,
               init_params=astropy.table.Table( { 'x_init': [511.2], 'y_init': [510.9], 'flux_init': [120000] } ) )
print( f"Fit flux of 100000-flux star for sigma = 0.42: {phot['flux_fit'][0]}; "
       f"(x,y)=({phot['x_fit'][0]}, {phot['y_fit'][0]}" )


# RESULTS OF A RUN:
#    root@655734bb3457:/snappl/experimentation# python play_with_photutils.py
#    1.2-sigma PSF clip sums to 1.0000000000018092
#    0.42-sigma PSF clip sums to 1.1267689222010988
#    Fit flux of 100000-flux star for sigma = 1.2: 99833.32184334665; (x,y)=(511.00130567398907, 510.99895105013127
#    Fit flux of 100000-flux star for sigma = 0.42: 100005.90585339173; (x,y)=(511.0002593061244, 511.0001099684552
#
# CONCLUSION:
#
# For PSFs whose FWHM is less than ~2 pixels on the original
# image, the stamps you get from the __call__ method of a photutils
# ImagePSF are not properly normalize.  (Interestingly, they are too big
# by the same factor that lanczos4 interpolation gives....)  However,
# photutils PSFPhotometry works in such a way that it takes this into
# account (...or just doesn't use __call__) and produces reliable results.
#
# Which is all well and good if you just want to use PSFPhotometry, but
# if you need thumbnails for other purposes (e.g. scene modelling),
# that is very sad.

[rknop]

As a note : simply dividing the returned clip by its sum is not good enough, because I need to be able to work in cases where a non trivial amount of the PSF's energy is outside the bounds of the clip.

[larrybradley]

Oversampled ImagePSF models should be normalized such that the sum of the array values is oversampling**2 (so 25 in your case). Is that what you are seeing?

[rknop]

If you look at the example code, I divide the image that I use to create the ImagePSF by oversampling**2, so that it will work properly with PSFPhotometry.

The issue isn't the image that's inside the ImagePSF object; it's the issue of a rendered PSF at the image scale, i.e. what you get by calling the ImagePSF object's __call__ method. With the same oversampling factor, a PSF with σ=3 pixels gets a rendered clip that sums to 1, but a PSF with σ=0.42 pixels gets a rendered clip that sums to 1.13.

[rknop]

And, I don't think that's a bug -- since the fitting methods come out with the right flux. So, whatever astropy's fitting is doing, it's assuming that __call__ will work this way.

What I want is a way to get an image clip out that is properly normalized.

[larrybradley]

This will require a deeper dive, but I have a suspicion that this is because when the model is evaluated on a discrete grid, the resulting values are simply samples from the underlying image (ref: #1262). That is how astropy modeling works. If I shift your example by 0.5 pixel in x, I get a sum of 0.9968356805968168. In principle, the sum should not change simply by shifting the model.

psf = make_gaussian_psf( oversamp=5, sigmax=0.42, stampsize=11, x0=511.5, y0=511 )
xvals = np.arange(511 - (29//2), 511 + (29//2) + 1)
yvals = np.arange(511 - (29//2), 511 + (29//2) + 1)
xvals, yvals = np.meshgrid(xvals, yvals)
psf(xvals, yvals).sum()

PSFPhotometry isn't doing anything special. I think the results are self consistent there because it evaluates the model the same way as the simulated image was created.

If this is indeed the cause, the solution would be to make a new "integrating" ImagePSF model class, but it's a bit tricky to implement (and will be slower).

[rknop]

Yeah, OK. I was hoping that this had been thought through before. I was also hoping that somebody knew math well enough to know if there was an interpolation method that would do a good job of faking the integration.