Box2DKernel & even oversampling

convolution_demo.py

import numpy as np
import scipy.ndimage
import scipy.interpolate
import scipy.special
from astropy.convolution import convolve, Box2DKernel


# some function we can oversample
def make_fake_image(oversample):
    pts = np.linspace(-10, 10, 10 * oversample + 1)
    xx, yy = np.meshgrid(pts, pts)
    zz = np.hypot(xx, yy)
    im = (2 * scipy.special.j1(zz) / (zz + (zz == 0)))**2
    im[zz == 0] = 1
    return xx, yy, im


def convolve_astropy(im, sz):
    return convolve(im, Box2DKernel(sz))


def convolve_scipy(im, sz):
    k2 = np.ones(sz * sz).reshape(sz, sz) / (sz * sz)
    c2 = scipy.ndimage.convolve(im, k2)
    if (sz % 2) == 0:
        c2 = scipy.ndimage.shift(c2, (0.5, 0.5))
        # scipy.ndimage.convolution by even kernel moves the center of the
        # array by a half pixel; we take this back out.
    return c2
    

# what's a truth value?  Let's define the 13x oversampled image to be truth.
# then for comparison we have to sample it at the appropriate places

def demo():
    truex, truey, trueimage = make_fake_image(13)
    trueimage = convolve_scipy(trueimage, 13)
    truei = scipy.interpolate.RegularGridInterpolator(
        (truey[:, 0], truex[0]), trueimage, method='cubic')
    for oversample in range(1, 14):
        x0, y0, im = make_fake_image(oversample)
        ca = convolve_astropy(im, oversample)
        cs = convolve_scipy(im, oversample)
        ct = truei(np.array([x0, y0]).T)
        sigma = np.std(ca - ct)
        print('%03d %8.6f %8.6f' % (
            oversample, np.std(ca - ct), np.std(cs - ct)))