numpy 频域图像旋转和缩放?

kfgdxczn  于 2023-01-26  发布在  其他
关注(0)|答案(4)|浏览(186)

我正在编写一些代码,利用相位相关性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
lymnna71

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]旋转。我不能完全解释这里到底发生了什么,但你也得到了我以前得到的相同效果。为了在傅里叶域中旋转原始图像,你必须首先将2D fftShift应用到原始图像M,然后计算FFT,旋转,IFFT,然后应用ifftShift。这样,图像的旋转中心和傅里叶域的中心就同步了。
AFAI还记得我们在两个独立的数组中旋转实部和虚部,然后合并它们。我们还测试了复数上的各种插值算法,没有太大的影响:)。它在我们的包pytom中。
然而,这可能是超级少,但与两个额外的移位不是真的很快,除非你指定一些时髦的数组索引算法。

mznpcxlj

mznpcxlj2#

那么,旋转和缩放图像的结果是旋转和缩放(与反比例)傅立叶变换。
还要注意,旋转和缩放在像素数量上都是线性的,而FFT是O(wlogwh*logh),所以实际上最终并不那么昂贵。

hrysbysz

hrysbysz3#

我意识到这有点晚了,但我只是想在回顾移位不变性的基础知识时回答这个问题。问题是你在旋转之前扩展了傅里叶空间(比如,考虑到混叠)。看看旋转图像的FT:轴向尖峰(别名)出现在边缘,在那里它们不在傅立叶旋转的IFT中。
你应该先旋转然后再处理混叠。因为你要考虑混叠(以周期=像素数循环傅立叶空间),然后通过旋转来消除这种影响,所以你会导致混叠出现在最终的图像中。本质上,你是在分散傅立叶混叠,从而将图像空间混叠拉到一起。
对于90度旋转,旋转工作平滑,因为没有混叠;k空间的角完美地匹配。

agxfikkp

agxfikkp4#

我敢肯定OP早就忘记了这个问题,但我最近花了整个下午试图回答这个问题,并认为我会分享我的结果,以防帮助其他人!
我发现有效的方法是结合原来的帖子和前面的答案,具体来说,需要执行两轮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)

我们可以看到,结果几乎与在真实的空间中应用旋转相同,只是在对象的边缘有一些值得注意的例外:

相关问题