我正在编写一些代码,利用相位相关性a la Reddy & Chatterji 1996恢复测试图像相对于模板的旋转、缩放和平移。我对原始测试图像进行FFT,以找到缩放因子和旋转Angular ,但我需要对 * 旋转和缩放 * 测试图像进行FFT,以获得平移。
现在,我可以在空间域中应用旋转和缩放,然后进行FFT,但这似乎有点低效-是否可以在频域中 * 直接 * 获得旋转/缩放图像的傅立叶系数?
**编辑1:**好吧,我按照用户1816548的建议做了一个尝试。对于90度的倍数,我可以得到看起来模糊的旋转,尽管图像极性会有奇怪的变化。而不是90度的倍数给予我得到的结果相当滑稽。
**编辑2:**我对图像进行了零填充,并且在旋转时包裹了FFT的边缘。我非常确定我是围绕FFT的直流分量旋转的,但对于不是90 °倍数的Angular ,我仍然会得到奇怪的结果。
示例输出:
可执行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
4条答案
按热度按时间lymnna711#
我不确定这个问题是否已经解决,但我相信我已经解决了你关于第三张图中观察到的效应的问题:
你观察到的这种奇怪的效果是由于你实际计算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中。
然而,这可能是超级少,但与两个额外的移位不是真的很快,除非你指定一些时髦的数组索引算法。
mznpcxlj2#
那么,旋转和缩放图像的结果是旋转和缩放(与反比例)傅立叶变换。
还要注意,旋转和缩放在像素数量上都是线性的,而FFT是O(wlogwh*logh),所以实际上最终并不那么昂贵。
hrysbysz3#
我意识到这有点晚了,但我只是想在回顾移位不变性的基础知识时回答这个问题。问题是你在旋转之前扩展了傅里叶空间(比如,考虑到混叠)。看看旋转图像的FT:轴向尖峰(别名)出现在边缘,在那里它们不在傅立叶旋转的IFT中。
你应该先旋转然后再处理混叠。因为你要考虑混叠(以周期=像素数循环傅立叶空间),然后通过旋转来消除这种影响,所以你会导致混叠出现在最终的图像中。本质上,你是在分散傅立叶混叠,从而将图像空间混叠拉到一起。
对于90度旋转,旋转工作平滑,因为没有混叠;k空间的角完美地匹配。
agxfikkp4#
我敢肯定OP早就忘记了这个问题,但我最近花了整个下午试图回答这个问题,并认为我会分享我的结果,以防帮助其他人!
我发现有效的方法是结合原来的帖子和前面的答案,具体来说,需要执行两轮fftshift/ifftshift;一个在真实的空间,另一个在傅立叶空间。
整个过程如下图所示:
复制此操作的代码是:
我们可以看到,结果几乎与在真实的空间中应用旋转相同,只是在对象的边缘有一些值得注意的例外: