computing the FT for non-uniform u,v coordinates

[Sketos]

I would like to compute the real and imaginary parts of V(u,v) from I(x,y) (a 2D image with a known pixel scale) for non-uniformly spaced u,v coordinates, where,

V(u,v) = sum_i^I sum_j^J I(x_i,y_j) exp(-2 * pi * i * (x_i * u + y_j * v))

where the u,v data are in units of wavelengths (or inverse radians) and x,y are in units of radians.

Is there a way to perform this calculation with nufft or does it only works for when x, y are non-uniformly space?

Below is an example script of the discrete FT (which becomes very slow for a large number of u,v points) that I want to compare against nufft.

import numpy as np
import matplotlib.pyplot as plt
from astropy import units

def z(x, y, sigma_x, sigma_y):
    return 1.0 / (2.0 * np.pi * sigma_x * sigma_y) * np.exp(-(x**2.0 / (2.0 * sigma_x**2) + y**2.0 / (2.0 * sigma_y**2.0)))


size = 5.0 # The angular size of the image in arcsec.
n_pixels = int(2**6.0)

# compute the pixel scale
pixel_scale = size / n_pixels

# generate the x,y regular grid
x_rad = np.linspace(
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)
y_rad = np.linspace(
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)
x_rad, y_rad = np.meshgrid(
    x_rad,
    y_rad
)

# make an image (example 1)
sigma_x = 0.5 # in units of arcsec
sigma_y = 0.5 # in units of arcsec
image = z(
    x=x_rad,
    y=y_rad,
    sigma_x=sigma_x * units.arcsec.to(units.rad),
    sigma_y=sigma_y * units.arcsec.to(units.rad)
)

# Generate the non-uniform u,v coordinates (in units of inverse radians or wavelengths)
N = 100
u = np.random.uniform(
    -1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    +1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    N)
v = np.random.uniform(
    -1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    +1.0 / 2.0 / (pixel_scale * units.arcsec.to(units.rad)) / 2.0,
    N)

# Generate visibilities using dFT
visibilities = np.zeros(
    shape=(u.shape[-1]),
    dtype="complex"
)
image_flatten = np.ndarray.flatten(image)
x_rad_flatten = np.ndarray.flatten(x_rad)
y_rad_flatten = np.ndarray.flatten(y_rad)
for j in range(u.shape[-1]):
    for i in range(image_flatten.shape[-1]):
        visibilities[j] += image_flatten[i] * np.exp(-2j * np.pi * (x_rad_flatten[i] * u[j] + y_rad_flatten[i] * v[j]))
visibilities_real__dFT = visibilities.real
visibilities_imag__dFT = visibilities.imag

plt.figure()
plt.scatter(
    visibilities_real__dFT,
    visibilities_imag__dFT,
    marker="o",
    alpha=0.85,
    label="dFT"
)
plt.xlabel("Real", fontsize=15)
plt.ylabel("Imag", fontsize=15)
plt.legend()
plt.show()

[ahbarnett]

Yes: use the 2d type-2 nufft.
PS you'll find FINUFFT faster
https://github.com/flatironinstitute/finufft