在 R 中使用 fft 执行相位相关

Question

我正在尝试使用维基百科 (http://en.wikipedia.org/wiki/Phase_correlation) 中的方法在 R 中实现二维相位相关算法，以跟踪 2 个图像之间的运动。这些图像（帧）是用在风中摇晃的相机拍摄的，最终目标是消除这些帧和后续帧中的抖动。两个示例图像和 R 代码如下：

## we will need the tiff library 
library(tiff)

## read in the tiff files 
f1=as.matrix(readTIFF('f1.tiff',native=TRUE))
f2=as.matrix(readTIFF('f2.tiff',native=TRUE))

## take the fft of the first  frame
F1 <- fft(f1)
## take the Conjugate fft of the second frame
F2.c <- Conj(fft(f2))

## calculate the cross power spectrum according to the wiki article
R <- (F1*F2.c)/abs(F1*F2.c)
## take the inverse fft of R
r <- fft(R,inv=TRUE)/length(R)
## because the zero valued imaginary numbers are not needed
r <- Re(r)

## show the normalized cross-correlation
image(r)

## find the max in the cross correlation matrix, or the phase shift -
## between the two images
shift <- which(r==max(r),arr.ind=TRUE)

根据我的理解，矢量 shift 应该包含有关最能校正这两个图像的传递位移（dx 和 dy）的信息。然而，shift 变量给出 dx=1 和 dy=1，我假设这表示在 x 或 y 方向上没有移动。对于在 x 和 y 方向上都有可见偏移或多个像素的后续帧，会发生这种情况。

你们有没有看到我的 code/formulas 中有错误？或者我是否需要先尝试一些更高级的东西，比如在进行相位相关之前先过滤图像？

小伙伴们干杯！

Answer 1

根据我对相位相关的了解，代码看起来是正确的。如果我理解你想要的是正确的，那么你正在尝试使用相位相关来确定两个图像之间的偏移，因为它们的单应性只不过是水平和垂直偏移。您只能在原点移动的事实很可能是由于您的图像缺乏足够的高频信息以正确确定良好的移动。

试试这两张图片（这些来自您引用的维基百科文章，但我将它们提取出来并保存为单独的图片）：

当我运行这两个图像与你的 R 代码时，我得到这个作为我的相位相关图。请记住，您的图像实际上保存为 .png，因此我不得不将库更改为 library(png)，并且我使用 readPNG 而不是 readTIFF。当您尝试运行您的代码与上述示例图像时请记住这一点：

此外，最大峰值出现的位置是：

> shift
     row col
[1,] 132 153

这告诉我们图像移动了 132 行和 153 列。请注意，这是 相对于图像的中心 。如果要确定实际偏移量，则需要将其减去垂直坐标行的一半和水平坐标列的一半。

因此，代码工作得很好......只是你的图像缺乏足够的高频信息来使相位相关工作。在这种情况下，相关性试图做的是我们试图找到每个图像之间的 "similar" 变化。如果每个图像之间有很多变化并且非常相似，那么相位相关将很有效。但是，如果我们没有那么大的变化，那么相位相关将不起作用。

为什么会这样？相位相关背后的基础是我们假设图像被高斯白噪声破坏，因此如果我们将白噪声与其自身相关联（从一个图像到另一个图像），它将在偏移或偏移处给出一个非常好的高峰几乎是零。由于您的图像缺少大量高频信息并且图像很干净，因此相位相关实际上不起作用。因此，有些人实际上建议 pre-whiten 你的图像，这样图像就包含白噪声，这样你就可以在我们应该偏移的地方获得漂亮的峰值重新谈论。

但是，为了确保消除任何错误的最大值，最好也平滑频域中的 cross-correlation 矩阵（r 在您的 R 代码中）以便很可能只有一个真正的最大值。在频域/FFT 域中使用高斯滤波器应该可以正常工作。

在任何情况下，我都没有看到您的图像有太大变化，因此需要注意的是，您必须确保您的图像具有大量高频信息才能正常工作！

Answer 2

下面是例程的定性描述，后面是 R 代码，用于有效和稳健地找到原始问题中发布的两个图像之间的相位相关性。感谢@rayreng 的建议并为我指明了正确的方向！

读入两张图片
为第二张图片添加噪点
使用 fft 将两个图像转换为频谱
使用乘法对转换后的图像进行卷积
return 到具有逆 fft 的空间域。这是您的归一化互相关
重新排列归一化互相关矩阵，使零频位于矩阵的中间（类似于matlab中的fftshift）
使用二维高斯分布平滑归一化互相关矩阵
通过识别平滑归一化相关矩阵的最大索引值来确定偏移

请注意，此函数还使用自定义二维高斯生成器（见下文）和类似于 matlabs fftshift 的自定义函数。

 ### R CODE ###
 rm(list=ls())
 library(tiff)
 ## read in the tiff images 
 f1 <- readTIFF('f1.tiff',native=TRUE)
 f1 <- matrix(f1,ncol=ncol(f1),nrow=nrow(f1)) 


 ## take the fft of f1 
 F1 <- fft(f1)

 ## what is the subsequent frame?
 f2 <-readTIFF('f2.tiff',native=TRUE)
 f2 <- matrix(f2,ncol=ncol(f2),nrow=nrow(f2))

 ## create a vector of random noise the same length as f2
 noise.b <- runif(length(f2),min(range(f2)),max(range(f2)))
 ## add the noise to the f2
 f2 <- noise.b+f2   

## take the fft and conjugate of the f2
F2.c <- Conj(fft(f2))

## calculate the cross-power spectrum
R <- (F1*F2.c)/abs(F1*F2.c)
## calculate the normalized cross-correlation with fft inverse
r <- fft(R,inv=TRUE)/length(R)
## rearrange the r matrix so that zero frequency component is in the -
## middle of the matrix.  This is similar to the function - 
## fftshift in matlab 

fftshift <- function(x) {
if(class(x)=='matrix') {
    rd2 <- floor(nrow(x)/2)
    cd2 <- floor(ncol(x)/2)

    ## Identify the first, second, third, and fourth quadrants 
    q1 <- x[1:rd2,1:cd2]
    q2 <- x[1:rd2,(cd2+1):ncol(x)]
    q3 <- x[(rd2+1):nrow(x),(cd2+1):ncol(x)]
    q4 <- x[(rd2+1):nrow(x),1:cd2]

    ## rearrange the quadrants 
    centered.t <- rbind(q4,q1)
    centered.b <- rbind(q3,q2)
    centered <- cbind(centered.b,centered.t)

    return(Re(centered))             
}
if(class(x)!='matrix') {
    print('sorry, this class of input x is not supported yet')
    }
}

## now use the defined function fftshift on the matrix r
r <- fftshift(r)
r <- Re(r)

## try and smooth the matrix r to find the peak!
## first build a 2d gaussian matrix  
sig = 5 ## define a sigma 

## determine the rounded half dimensions of the r matrix 
x.half.dim <- floor(ncol(r)/2)
y.half.dim <- floor(nrow(r)/2)

## create x and y vectors that correspond to the indexed row
## and column values of the r matrix.  
xs <- rep(-x.half.dim:x.half.dim,ncol(r))
ys <- rep(-y.half.dim:y.half.dim,each=nrow(r))

## apply the gaussian blur formula 
## (http://en.wikipedia.org/wiki/Gaussian_blur)
## to every x and y indexed value
gaus <- 1/(2*pi*sig^2) * exp(-(xs^2 + ys^2)/(2*sig^2))
gaus <- matrix(gaus,ncol=ncol(r),nrow=nrow(r),byrow=FALSE)

## now convolve the gaus with r in the frequency domain
r.filt <- fft((fft(r)*fft(gaus)),inv=TRUE)/length(r)
r.filt <- fftshift(Re(r.filt))

## find dx and dy with the peak in r    
min.err <- which(r.filt==max(r.filt),arr.ind=TRUE)

## how did the image move from the previous? 
shift <- (dim(f1)+3)/2-min.err

矢量偏移表示图像在x 正方向和y 负方向偏移。换句话说，第二张图片（f2）被粗略地移动到了右上角。由于第二张图像中引入的噪声以及来自 r 矩阵上的高斯滤波器的平滑算子，每次尝试时矢量偏移的值都会略有不同。我注意到类似于上面概述的相位相关在较大的 images/matrices 上效果更好。要查看上述算法的结果，请访问 https://www.youtube.com/watch?v=irDFk2kbKaE 处的稳定视频。

在 R 中使用 fft 执行相位相关

Performing a phase correlation with fft in R

r

image-processing

fft

cross-correlation