ali*_*i_m 7 math numpy transformation fft image-processing
我正在编写一些代码来恢复测试图像相对于使用相位相关的模板的旋转、缩放和平移,a la Reddy & Chatterji 1996。我采用原始测试图像的 FFT 以找到比例因子和旋转角度,但我需要旋转和缩放的测试图像的 FFT以获得平移。
现在我可以在空间域中应用旋转和缩放,然后进行 FFT,但这似乎有点低效 - 是否可以直接在频域中获得旋转/缩放图像的傅立叶系数?
编辑 1: 好的,我按照 user1816548 的建议进行了尝试。对于 90 度的倍数,我可以得到模糊合理的旋转,尽管图像的极性发生了奇怪的变化。不是 90o 倍数的角度给了我非常滑稽的结果。
编辑 2: 我对图像应用了零填充,并且在旋转 FFT 时将其边缘包裹起来。我很确定我正在围绕 FFT 的 DC 分量旋转,但是对于不是 90o 倍数的角度,我仍然得到奇怪的结果。
示例输出:
可执行的 Numpy/Scipy 代码:
import numpy as np
from scipy.misc import lena
from scipy.ndimage.interpolation import rotate,zoom
from scipy.fftpack import fft2,ifft2,fftshift,ifftshift
from matplotlib.pyplot import subplots,cm
def testFourierRotation(angle):
M = lena()
newshape = [2*dim for dim in M.shape]
M = procrustes(M,newshape)
# rotate, then take the FFT
rM = rotate(M,angle,reshape=False)
FrM = fftshift(fft2(rM))
# take the FFT, then rotate
FM = fftshift(fft2(M))
rFM = rotatecomplex(FM,angle,reshape=False)
IrFM = ifft2(ifftshift(rFM))
fig,[[ax1,ax2,ax3],[ax4,ax5,ax6]] = subplots(2,3)
ax1.imshow(M,interpolation='nearest',cmap=cm.gray)
ax1.set_title('Original')
ax2.imshow(rM,interpolation='nearest',cmap=cm.gray)
ax2.set_title('Rotated in spatial domain')
ax3.imshow(abs(IrFM),interpolation='nearest',cmap=cm.gray)
ax3.set_title('Rotated in Fourier domain')
ax4.imshow(np.log(abs(FM)),interpolation='nearest',cmap=cm.gray)
ax4.set_title('FFT')
ax5.imshow(np.log(abs(FrM)),interpolation='nearest',cmap=cm.gray)
ax5.set_title('FFT of spatially rotated image')
ax6.imshow(np.log(abs(rFM)),interpolation='nearest',cmap=cm.gray)
ax6.set_title('Rotated FFT')
fig.tight_layout()
pass
def rotatecomplex(a,angle,reshape=True):
r = rotate(a.real,angle,reshape=reshape,mode='wrap')
i = rotate(a.imag,angle,reshape=reshape,mode='wrap')
return r+1j*i
def procrustes(a,target,padval=0):
b = np.ones(target,a.dtype)*padval
aind = [slice(None,None)]*a.ndim
bind = [slice(None,None)]*a.ndim
for dd in xrange(a.ndim):
if a.shape[dd] > target[dd]:
diff = (a.shape[dd]-target[dd])/2.
aind[dd] = slice(np.floor(diff),a.shape[dd]-np.ceil(diff))
elif a.shape[dd] < target[dd]:
diff = (target[dd]-a.shape[dd])/2.
bind[dd] = slice(np.floor(diff),target[dd]-np.ceil(diff))
b[bind] = a[aind]
return b
Run Code Online (Sandbox Code Playgroud)
我不确定这是否已经解决,但我相信我有解决您的问题的解决方案,即您在第三张图中观察到的效果:
您观察到的这种奇怪的效果是由于您实际计算 FFT 的来源。本质上,FFT 从阵列的第一个像素开始M[0][0]
。但是,您定义了围绕 的旋转M[size/2+1,size/2+1]
,这是正确的方法,但错误 :) 。傅立叶域已计算出M[0][0]
!如果你现在在傅立叶域中旋转,你是在旋转M[0][0]
而不是围绕M[size/2+1,size/2+1]
。我无法完全解释这里到底发生了什么,但你也得到了和我以前一样的效果。为了在傅立叶域中旋转原始图像,您必须首先将 2DfftShift
应用于原始图像 M,然后计算 FFT、旋转、IFFT,然后应用ifftShift
。这样您的图像旋转中心和傅立叶域的中心就会同步。
AFAI 记得我们还在两个单独的数组中旋转实部和虚部,然后将它们合并。我们还在复数上测试了各种插值算法,但效果不大:)。它在我们的包pytom 中。
然而,这可能会非常少,但是两个额外的移位并不是很快,除非你指定了一些时髦的数组索引算法。
那么,旋转和缩放的图像会产生旋转和缩放(具有逆比例)傅立叶变换。
另请注意,旋转和缩放在像素数量上都是线性的,而 FFT 是 O(w*logw*h*logh),因此最终实际上并没有那么昂贵。