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

Question

我正在尝试使用来自维基百科 (http://en.wikipedia.org/wiki/Phase_correlation) 的配方在 R 中实现 2d 相位相关算法，以跟踪 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）的信息。然而，移位变量给出 dx=1 和 dy=1，我假设这表明在 x 或 y 方向上都没有移位。对于在 x 和 y 方向上都有可见偏移或多个像素的后续帧，会发生这种情况。

你们中有人在我的代码/公式中看到错误吗？还是在进行相位相关之前，我是否需要先尝试一些更高级的方法，例如过滤图像？

男女朋友们干杯！

Answer 1

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

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

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

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

> shift
     row col
[1,] 132 153

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

因此，代码完全可以正常工作...只是您的图像缺少足够的高频信息来使相位相关起作用。在这种情况下，相关性试图做的是我们试图找到每个图像之间的“相似”变化。如果每个图像之间有很多变化并且非常相似，那么相位相关就会很好地工作。但是，如果我们没有那么多变化，那么相位相关将不起作用。

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

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

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

Answer 2

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

1. 阅读两张图片
2. 为第二张图片添加噪点
3. 使用 fft 将两个图像转换为频谱
4. 使用乘法对变换后的图像进行卷积
5. 以逆 fft 返回到空间域。这是您的标准化互相关
6. 重新排列归一化的互相关矩阵，使零频率在矩阵的中间（类似于matlab中的fftshift）
7. 使用二维高斯分布平滑归一化互相关矩阵
8. 通过确定平滑归一化相关矩阵的最大索引值来确定偏移

9. 请注意，此函数还使用自定义二维高斯生成器（见下文）和类似于 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 矩阵上的高斯滤波器的平滑算子，矢量偏移的值会随着每次尝试而略有不同。我注意到类似于上面概述的相位相关在较大的图像/矩阵上效果更好。要查看上述算法的结果，请访问https://www.youtube.com/watch?v=irDFk2kbKaE 的稳定视频。