以下代码在通过傅里叶相移移动图像时创建伪影:
def phase_shift(fimage, dx, dy):
# Shift the phase of the fourier transform of an image
dims = fimage.shape
x, y = np.meshgrid(np.arange(-dims[1] / 2, dims[1] / 2), np.arange(-dims[0] / 2, dims[0] / 2))
kx = -1j * 2 * np.pi * x / dims[1]
ky = -1j * 2 * np.pi * y / dims[0]
shifted_fimage = fimage * np.exp(-(kx * dx + ky * dy))
return shifted_fimage
实际移动图像并获取移动图像的用法:
def translate_by_phase_shift(image, dx, dy):
# Get the fourier transform
fimage = np.fft.fftshift(np.fft.fftn(image))
# Phase shift
shifted_fimage = phase_shift(fimage, dx, dy)
# Inverse transform -> translated image
shifted_image = np.real(np.fft.ifftn(np.fft.ifftshift(shifted_fimage)))
return shifted_image
cv2.warpAffine()
使用相同的班次。
[更新]建议使用的注释之一
scipy.ndimage.fourier.fourier_shift()
fourier_shifted_image = fourier_shift(np.fft.fftn(image), shift)
shifted_image = np.fft.ifftn(fourier_shifted_image)
并绘制了真实的部分(
shifted_image.real
)
phase_shift()
[更新]现在我们排除了我的phase_shift()函数,这里有一个可复制的代码,前提是您从这里下载图像阵列:
https://www.dropbox.com/s/dmbv56xfqkv8qqz/image.npy?dl=0
import os
import numpy as np
import matplotlib
matplotlib.use('TKAgg')
import matplotlib.pyplot as plt
from scipy.ndimage.fourier import fourier_shift
# Load the image (update path according to your case)
image = np.load(os.path.expanduser('~/DDS/46P_Wirtanen/image.npy'))
# Shift vector
shift = np.array([-3.75, -7.5 ])
# Phase-shift
fourier_shifted_image = fourier_shift(np.fft.fftn(image), shift)
shifted_image = np.fft.ifftn(fourier_shifted_image)
interp_method = 'hanning'
zoomfov = [1525, 1750, 1010, 1225]
vmin = np.percentile(image, 0.1)
vmax = np.percentile(image, 99.8)
fig, ax = plt.subplots(1,2, figsize=(14, 6), sharex=True,sharey=True)
ax[0].imshow(image, origin='lower', cmap='gray', vmin=vmin, vmax=vmax, interpolation=interp_method)
ax[0].set_title('Original image')
ax[1].imshow(shifted_image.real, origin='lower', cmap='gray', vmin=vmin, vmax=vmax, interpolation=interp_method)
ax[1].set_title('with scipy.ndimage.fourier.fourier_shift()')
plt.axis(zoomfov)
plt.tight_layout()
plt.show()
输出如下所示:
[更新]
根据Cris的回复,我使用opencv中的其他插值方法对强度进行对数缩放,得出了类似的结论:该伪影确实也存在于图像中的Lanczos标志中
要达到此目的的代码:
# Compare interpolation methods
import cv2
# Fourier phase shift.
fourier_shifted = fourier_shift(np.fft.fftn(image), shift)
fourier_shifted_image = np.fft.ifftn(fourier_shifted).real
# Use opencv
Mtrans = np.float32([[1,0,shift[1]],[0,1, shift[0]]])
shifted_image_cubic = cv2.warpAffine(image, Mtrans, image.shape[::-1], flags=cv2.INTER_CUBIC)
shifted_image_lanczos = cv2.warpAffine(image, Mtrans, image.shape[::-1], flags=cv2.INTER_LANCZOS4)
zoomfov = [1525, 1750, 1010, 1225]
pmin = 2
pmax = 99.999
fig, ax = plt.subplots(1,3, figsize=(19, 7), sharex=True,sharey=True)
ax[0].imshow(fourier_shifted_image, origin='lower', cmap='gray',
vmin=np.percentile(fourier_shifted_image, pmin), vmax=np.percentile(fourier_shifted_image, pmax),
interpolation=interp_method, norm=LogNorm())
add_rectangle(zoomfov, ax[0])
ax[0].set_title('shifted with Fourier phase shift')
ax[1].imshow(shifted_image_cubic, origin='lower', cmap='gray',
vmin=np.percentile(shifted_image_cubic, pmin), vmax=np.percentile(shifted_image_cubic, pmax),
interpolation=interp_method, norm=LogNorm())
add_rectangle(zoomfov, ax[1])
ax[1].set_title('with cv2.warpAffine(...,flags=cv2.INTER_CUBIC)')
ax[2].imshow(shifted_image_lanczos, origin='lower', cmap='gray',
vmin=np.percentile(shifted_image_lanczos, pmin), vmax=np.percentile(shifted_image_lanczos, pmax),
interpolation=interp_method, norm=LogNorm())
#ax[2].imshow(shifted_image.real, origin='lower', cmap='gray', vmin=np.percentile(Llights_prep[frame], pmin), vmax=np.percentile(Llights_prep[frame], pmax), interpolation=interp_method)
add_rectangle(zoomfov, ax[2])
ax[2].set_title('with cv2.warpAffine(...,flags=cv2.INTER_LANCZOS4) ')
plt.axis(zoomfov)
plt.tight_layout()
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=None)
plt.show()
为了回答克里斯的问题,我们的业余成像系统(直径为130毫米的成像系统)确实无法避免采样不足的恒星。我天真地应用了相同的算法,而不是我在专业的、更大的仪器上使用的算法,而这些仪器并没有显示出这个问题。
