All-zero row and column in the PSF model

[kirxkirx]

Strangely, the model PSF constructed with EPSFBuilder always has an all-zero row and an all-zero column next to the edges of the model PSF image. Here is how it looks like (note the black row on the bottom and the black column on the right, while the corresponding top row and left column have a range of pixels values):

I was experimenting with extracting PSF from smallish cutouts which made it easier to notice (with a large cutout the PSF would fall all the way down to zero anyhow). First noticed on real data, I realized that the all-zero bands appear also when reconstructing PSF from a simulated image with Gaussian stars. Here is the code to reproduce the above image (based on test_epsf.py):

import numpy as np
import matplotlib.pyplot as plt

from scipy.spatial import cKDTree
from astropy.table import Table
from astropy.visualization import simple_norm
from photutils.datasets import make_gaussian_prf_sources_image
from astropy.nddata import NDData
from photutils.psf import extract_stars, EPSFBuilder

# Generate a test image with Gaussian stars
shape = (750, 750)
nstars = 100
rng = np.random.default_rng(0)
xx = rng.uniform(low=0, high=shape[1], size=nstars)
yy = rng.uniform(low=0, high=shape[0], size=nstars)
# enforce a minimum separation
min_dist = 25
coords = [(yy[0], xx[0])]
for xxi, yyi in zip(xx, yy):
    newcoord = [yyi, xxi]
    dist, _ = cKDTree([newcoord]).query(coords, 1)
    if np.min(dist) > min_dist:
        coords.append(newcoord)
        yy, xx = np.transpose(coords)
        zz = rng.uniform(low=0, high=200000, size=len(xx))
# define a table of model parameters
stddev = 2.0
sources = Table()
sources['amplitude'] = zz
sources['x_0'] = xx
sources['y_0'] = yy
sources['sigma'] = np.zeros(len(xx)) + stddev
sources['theta'] = 0.0
data = make_gaussian_prf_sources_image(shape, sources)
nddata = NDData(data)
init_stars = Table()
init_stars['x'] = xx.astype(int)
init_stars['y'] = yy.astype(int)

# Display the generated test image
norm = simple_norm(data, 'sqrt', percent=99)
fig, ax = plt.subplots(figsize=(5,5))
ax.imshow(data, norm=norm, cmap='gray')
ax.set_title('Simulated image')
plt.show()

# Extract PSF from the simulated image
size = 11
oversampling = 4
stars = extract_stars(nddata, init_stars, size=size)
epsf_builder = EPSFBuilder(oversampling=oversampling, maxiters=15,progress_bar=False, norm_radius=25,recentering_maxiters=15)
epsf, fitted_stars = epsf_builder(stars)

# Display oversmpled PSF
norm = simple_norm(epsf.data, 'sqrt', percent=99)
fig, ax = plt.subplots(figsize=(5,5))
ax.imshow(epsf.data, norm=norm, cmap='gray')
ax.set_title('Model PSF')
ax.set_xlabel('X [pixels]')
ax.set_ylabel('Y [pixels]')
plt.show()
print(epsf.data)

# Display downsampled PSF
model_psf = epsf.data[::oversampling, ::oversampling]
norm = simple_norm(model_psf, 'sqrt', percent=99)
fig, ax = plt.subplots(figsize=(5,5))
ax.imshow(model_psf, norm=norm, cmap='gray')
ax.set_title('Model (downsampled) PSF')
ax.set_xlabel('X [pixels]')
ax.set_ylabel('Y [pixels]')
plt.show()
print(model_psf)



# Print version info
import photutils
print('\n photutils version: ')
print(photutils.__version__)

import astropy
print('\n astropy version: ')
print(astropy.__version__)

print('\n numpy version: ')
print(np.__version__)

The final lines of this code's output showing the numerical values of the image and the photutils/astropy/numpy versions on my machine:

[[5.34884925e-05 1.65266053e-04 4.30301723e-04 8.82556579e-04
  1.43279909e-03 1.84049074e-03 1.86156378e-03 1.48370997e-03
  9.32943599e-04 4.62742539e-04 1.81343626e-04 0.00000000e+00]
 [1.68164780e-04 2.79156978e-04 7.46585480e-04 1.55276556e-03
  2.52816107e-03 3.22377239e-03 3.21844399e-03 2.51643878e-03
  1.54051939e-03 7.38353698e-04 2.77416584e-04 0.00000000e+00]
 [4.36908762e-04 7.36746093e-04 1.97514726e-03 4.11116497e-03
  6.69872498e-03 8.54524163e-03 8.53375181e-03 6.67176669e-03
  4.08329663e-03 1.95625675e-03 7.34668659e-04 0.00000000e+00]
 [8.93128915e-04 1.53563363e-03 4.11520741e-03 8.56607166e-03
  1.39573624e-02 1.78037681e-02 1.77801834e-02 1.39012092e-02
  8.50809205e-03 4.07604573e-03 1.53059626e-03 0.00000000e+00]
 [1.42770481e-03 2.50528209e-03 6.71211603e-03 1.39719249e-02
  2.27663441e-02 2.90409250e-02 2.90018436e-02 2.26737553e-02
  1.38769815e-02 6.64813126e-03 2.49638838e-03 0.00000000e+00]
 [1.78331757e-03 3.19887284e-03 8.56953853e-03 1.78386666e-02
  2.90680397e-02 3.70811897e-02 3.70314987e-02 2.89512996e-02
  1.77185499e-02 8.48860383e-03 3.18725844e-03 0.00000000e+00]
 [1.74836072e-03 3.19563659e-03 8.56413013e-03 1.78274289e-02
  2.90507594e-02 3.70604466e-02 3.70118045e-02 2.89360352e-02
  1.77092919e-02 8.48420113e-03 3.18555858e-03 0.00000000e+00]
 [1.35292246e-03 2.49890613e-03 6.70027289e-03 1.39474242e-02
  2.27284615e-02 2.89947247e-02 2.89560716e-02 2.26380541e-02
  1.38548243e-02 6.63768414e-03 2.49240475e-03 0.00000000e+00]
 [8.17550154e-04 1.52971704e-03 4.10366591e-03 8.54234493e-03
  1.39200735e-02 1.77573885e-02 1.77337367e-02 1.38642689e-02
  8.48506640e-03 4.06507649e-03 1.52651843e-03 0.00000000e+00]
 [3.89808658e-04 7.33293932e-04 1.96742461e-03 4.09552423e-03
  6.67346989e-03 8.51299927e-03 8.50158080e-03 6.64650177e-03
  4.06767512e-03 1.94876800e-03 7.31876961e-04 0.00000000e+00]
 [1.47290729e-04 2.74493398e-04 7.38852204e-04 1.53800270e-03
  2.50600915e-03 3.19664294e-03 3.19222093e-03 2.49545774e-03
  1.52711179e-03 7.31517171e-04 2.74846223e-04 0.00000000e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
  0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00]]

 photutils version: 
1.6.0

 astropy version: 
5.1

 numpy version: 
1.23.3

The oversampled PSF also has one all-zero row and one all-zero column:

What am I doing wrong? Or maybe the all-zero row and column are supposed to be there for a reason? I'm concerned because I'm using the PSF model for image subtraction. I'm trying to track down the origins of image subtraction imperfections.