在频域中旋转和缩放图像？

Question

我正在编写一些代码，使用相位相关法a la Reddy & Chatterji 1996来恢复测试图像相对于模板的旋转、缩放和平移。我对原始测试图像进行FFT以找到比例因子和旋转角度，但我需要对旋转和缩放后的测试图像进行FFT才能获得平移信息。

现在我可以在空间域中应用旋转和缩放，然后进行FFT，但这似乎有点低效 - 是否可能直接在频率域中获得旋转/缩放后的图像的傅里叶系数？

编辑1： 好的，我按照user1816548的建议尝试了一下。对于90度的倍数角度，我可以得到看起来比较合理的旋转结果，但是图像的极性会发生奇怪的变化。对于不是90度的倍数角度，我的结果非常离谱。

编辑2： 我对图像进行了零填充，并在旋转时将FFT的边缘包装起来。我非常确定我正在围绕FFT的DC分量进行旋转，但是对于不是90度的倍数角度，我仍然得到奇怪的结果。

示例输出：

---

可执行的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

- ali_m

---

我也对这个很感兴趣。如果你查看空间旋转图像的FFT，你会看到从水平轴上的400处向垂直中心线延伸的尖峰。但是这些线条在旋转后的FFT中并不存在。另外，你只绘制了振幅，相位图呢？你能把它们也发出来吗？谢谢。 - Emily L.

Answer 1

我不确定这个问题是否已经解决，但我相信我有关于你第三张图中观察到的效果的解决方案：

你观察到的奇怪效果是由于你实际计算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]旋转。我无法完全解释这里发生了什么，但你得到了我曾经得到的相同效果。为了在傅里叶域中旋转原始图像，你必须首先对原始图像M应用二维fftShift，然后计算FFT，旋转，IFFT，然后应用ifftShift。这样你的图像旋转中心和傅里叶域的中心就会同步。

据我记得，我们还在分别旋转实部和虚部的两个数组，然后将它们合并。我们还对复数进行了各种插值算法的测试，但效果不大。这在我们的包pytom中。

然而，除非你指定一些花哨的数组索引算术，否则这可能会很慢，而且需要两次额外的移位。

- El Dude

Answer 2

我相信提问者早已忘记了这个问题，但我最近花了整个下午来回答这个问题，并想在这里分享我的结果，以帮助其他人！

我发现有效的方法是将原始帖子和之前的回答结合起来。具体来说，需要在实空间和傅里叶空间中分别执行两轮fftshift/ifftshift。

完整的过程如下图所示： 复制此代码的方法如下：

im = plt.imread(fpath) # this example has shape (1000,1000)

# Pad array such that image is twice the original size
# to avoid Fourier wraparound artefact
impad = np.pad(im,(500,500))

# Apply fftshift in real space before fft2 (!)
fft = np.fft.fftshift(impad)

# Do the regular fft --> fftshift --> rotate --> ifftshift --> ifft2
fft = np.fft.fft2(fft)
fft = np.fft.fftshift(fft)
fft = ndimage.rotate(fft,30.7,reshape=False)
fft = np.fft.ifftshift(fft)
out = np.fft.ifft2(fft)

# Apply ifftshift in real space after ifft2 (!)
out = np.fft.ifftshift(out)

# Crop back to original size
out=out[500:-500,500:-500]

# Take absolute value
im_rotated = abs(out)

我们可以看到，使用傅里叶变换旋转的结果几乎与在实际空间中应用旋转相同，在物体的边缘有一些明显的异常情况：

- George Lewis

Answer 3

旋转和缩放图像会导致旋转和缩放（具有反比例缩放）的傅里叶变换。

此外，请注意，旋转和缩放在像素数量上都是线性的，而FFT的时间复杂度是O(w*logw*h*logh)，因此最终实际上并不那么昂贵。

- tjltjl

---

嗨，谢谢你的回复。我可能最终会在空间域中进行转换，但无论如何，我都想弄清楚在傅里叶域中如何工作 - 请参见我的编辑。 - ali_m

确保您的旋转以正确的位置为中心-根据您的约定，它是中心或左上角。基本上，图像的最亮点是您应该围绕旋转的点。您可以通过取一些相当规则的图像并转换图像及其FT来检查它们之间的关系。 - tjltjl

你是在DCT周围精确地旋转还是只是大致旋转？因为数组长度是偶数（例如8），如果不围绕DCT像素的中心旋转，则DCT周围的旋转将会出错。我认为你得到的结果只是围绕DCT像素的角落旋转所得到的结果... - tjltjl

好的，我尝试围绕DCT像素的确切中心旋转（fftshift将其放置在偶数和奇数大小的维度上的ceil((dim-1)/2.)）。即使在旋转后FFT中最亮的像素仍保持在同一位置，但在进行IFFT时仍会出现奇怪的“万花筒”效果。 - ali_m

我得等有空的时候再好好看看这个... 这变得越来越有趣了 :) - tjltjl

Answer 4

我意识到现在回答有些晚了，但是我想在复习移不变性的基础知识时在这里回答这个问题。问题在于你在旋转之前扩展了傅里叶空间（也就是考虑混淆），看一下旋转图像的傅里叶变换：轴向尖峰（别名）出现在边缘，在傅里叶旋转的反变换中则没有。

你应该先旋转，然后再处理混淆。因为你正在考虑混淆（在像素数量等于周期的情况下循环你的傅里叶空间），然后通过旋转抛弃了那个努力，导致别名出现在最终图像中。本质上，你正在扩散傅里叶别名，因此将图像空间别名聚集在一起。

90度旋转的情况下旋转平稳无误，因为没有混淆；k空间的角落完美匹配。

- Chris L.