сравнивая ДПФ и БПФ для неравномерно расположенных точек u, v - PullRequest
Новая
       96

сравнивая ДПФ и БПФ для неравномерно расположенных точек u, v

0 голосов
/ 05 февраля 2020

Я хотел бы вычислить действительные и мнимые части V (u, v) из I (x, y) (2D-изображение с известным масштабом пикселей)

enter image description here

где точки u, v выбраны неравномерно (координаты u, v даны в единицах длин волн или обратных радианах, а x, y в единицах радианов). Стандартный подход заключается в выполнении дискретного преобразования Фурье, но это быстро становится очень медленным для многих точек u, v.

Я видел где-то альтернативный подход к ускорению вещей.

i) БПФ для изображения I (x, y).

ii) интерполяция I_fft в сетку Фурье-пространства.

iii) прогнозирование V (u, v) ) из интерполированного I_fft.

Ниже приведен пример сценария, демонстрирующий два метода. Я считаю, что второй метод значительно отличается от первого метода (ДПФ; «правда»). Что-то не так с моим расчетом? Полагаю, это неудивительно, но я удивлен разницей.

Если предположить, что в моих расчетах нет ничего плохого, есть ли лучшее ускорение? Я начал изучать программное обеспечение для неоднородного быстрого преобразования Фурье (NUFFT) (например, pynufft на GitHub "https://github.com/jyhmiinlin/pynufft"), но мне не удалось найти работающий пример для этого. Кто-нибудь использовал это раньше для этой цели? Есть ли другие альтернативные подходы, которые я могу использовать?

Scipt: 

import numpy as np
import matplotlib.pyplot as plt
from astropy import units
from scipy import fftpack, interpolate

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)
uv_distance = np.sqrt(u**2.0 + v**2.0)

# Generate visibilities using dFT (method 1)
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

# Generate visibilities using FFT (method 2)
image_FFT = fftpack.fftshift(fftpack.fft2(fftpack.fftshift(image)))
u_grid = fftpack.fftshift(np.fft.fftfreq(n_pixels, d=pixel_scale * units.arcsec.to(units.rad)))
v_grid = fftpack.fftshift(np.fft.fftfreq(n_pixels, d=pixel_scale * units.arcsec.to(units.rad)))
spliner = interpolate.RectBivariateSpline(u_grid, v_grid, image_FFT.real, kx=1, ky=1)
splinei = interpolate.RectBivariateSpline(u_grid, v_grid, image_FFT.imag, kx=1, ky=1)
visibilities_real__FFT = spliner.ev(u, v)
visibilities_imag__FFT = splinei.ev(u, v)

# NOTE : Figure
figure, axes = plt.subplots(
    nrows=1,
    ncols=2,
    figsize=(10,5)
)
im0 = axes[0].imshow(
    image_FFT.real,
    extent=[np.min(u_grid)*10**-3.0, np.max(u_grid)*10**-3.0, np.min(v_grid)*10**-3.0, np.max(v_grid)*10**-3.0],
    cmap="jet")
plt.colorbar(im0, ax=axes[0])
im1 = axes[1].imshow(
    image_FFT.imag,
    extent=[np.min(u_grid)*10**-3.0, np.max(u_grid)*10**-3.0, np.min(v_grid)*10**-3.0, np.max(v_grid)*10**-3.0],
    cmap="jet")
plt.colorbar(im1, ax=axes[1])
axes[0].plot(u*10**-3.0, v*10**-3.0, linestyle="None", marker=".", color="black")
axes[1].plot(u*10**-3.0, v*10**-3.0, linestyle="None", marker=".", color="black")
axes[0].set_xlabel(r"u (k$\lambda$)", fontsize=15)
axes[1].set_xlabel(r"u (k$\lambda$)", fontsize=15)
axes[0].set_ylabel(r"v (k$\lambda$)", fontsize=15)
axes[1].set_yticks([])
axes[0].set_title("Real", fontsize=15)
axes[1].set_title("Imag", fontsize=15)
plt.subplots_adjust(wspace=0.0)

plt.figure()
plt.scatter(
    visibilities_real__dFT,
    visibilities_imag__dFT,
    marker="o",
    alpha=0.85,
    label="dFT"
)
plt.scatter(
    visibilities_real__FFT,
    visibilities_imag__FFT,
    marker="o",
    alpha=0.85,
    label="pynufft"
)
for i in range(N):
    plt.plot(
        [visibilities_real__dFT[i], visibilities_real__FFT[i]],
        [visibilities_imag__dFT[i], visibilities_imag__FFT[i]],
        color="black"
    )
plt.xlabel("Real", fontsize=15)
plt.ylabel("Imag", fontsize=15)
plt.legend()
plt.show()
...