在R中使用FFT进行相位相关性计算

Question

在R中使用FFT进行相位相关性计算

3

我正在尝试使用维基百科上的一份配方（http://en.wikipedia.org/wiki/Phase_correlation）在R中实现2D相位相关算法，以便跟踪两个图像之间的移动。这些图像（帧）是由于摄像机在风中晃动而捕获的，最终目标是消除这些和随后的帧中的抖动。下面是两个示例图像和R代码：

frame 1 frame 2

## 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)

< p >据我理解，向量shift应包含有关最佳矫正这两幅图像的平移位移信息（dx和dy）。但是，变量shift给出dx=1和dy=1，我认为这表示x或y方向都没有移动。这在后续帧中发生，在这些帧中存在可见的x和y方向上的多个像素移动。

您是否看到我的代码/公式中的任何错误？还是我需要尝试一些更高级的东西，比如在进行相位相关之前先过滤图像？

祝大家好运！

- Alex Witsil

嗨，我在研究使用相位相关进行图像拼接时偶然发现了这个。请问为什么要使用R而不是简单地使用Matlab或Python（如果许可证成本是不使用Matlab的主要原因）？ - bFig8

2个回答

4

以下是用于有效且稳健地查找原始问题中发布的两个图像之间的相位相关性的R代码所遵循的常规定性描述。感谢@rayreng的建议和指引我正确的方向！

Read in both images
Add noise to the second image
Transform both images into the frequency spectrum using fft
convolve transformed images using a multiplication
return to the spatial domain with a inverse fft. This is your normalized cross correlation
Rearrange the normalized cross correlation matrix so that the zero frequency is in the middle of the matrix (similar to fftshift in matlab)
Use a 2d Gaussian distribution to smooth the normalized cross correlated matrix
Determine the shift by identifying the maximum indexed value of the smoothed normalized correlated matrix

Be aware that this function also uses a custom 2d Gaussian generator (see below) and a custom function similar to 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。

- Alex Witsil

1

不错！维基百科的文章还提到，您可以通过先将图像转换为对数极坐标来扩展此方法以提取旋转或缩放的可能差异。您知道如何做到这一点吗？文章是B. S Reddy和B. N. Chatterji，“用于平移、旋转和尺度不变图像配准的基于FFT的技术”，IEEE图像处理交易5，第8期（1996年）：1266-1271。 - Tom Wenseleers

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- rayryeng · Accepted Answer

根据我对相位相关的了解，代码看起来是正确的。如果我正确理解你的意思，你试图使用相位相关来确定两个图像之间的偏移量，给定它们的单应性只是水平和垂直偏移。您只得到从原点开始的移位，这很可能是由于您的图像缺乏足够的高频信息，以便正确确定良好的移位。

请尝试使用以下两个图像（这些来自您引用的维基百科文章，但我将它们提取出来并将它们保存为单独的图像）：

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

此外，最大峰值发生的位置是：

> shift
     row col
[1,] 132 153

这告诉我们图像向下移动了132行和向右移动了153列。请注意，这是相对于图像中心的位置。如果您想确定实际偏移量，您需要从垂直坐标减去半行数，并从水平坐标减去半列数。

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

为什么会这样？相位相关背后的基础是我们假设图像受到高斯白噪声的污染，因此，如果我们将白噪声与自身（从一个图像到另一个图像）相关联，它将在偏移或移位处给出非常好的高峰，而几乎在其他地方都是零。由于您的图像缺乏许多高频信息并且图像干净，因此相位相关实际上不起作用。因此，一些人建议预处理图像，使图像包含白噪声，以便您可以在所讨论的偏移位置获得漂亮的峰值。

但是，为了确保消除任何虚假的最大值，在频率域中（您R代码中的r）平滑交叉相关矩阵也是一个好主意，以便只有一个真正的最大值的概率很高。在频率/ FFT域中使用高斯过滤器应该可以正常工作。

无论如何，我没有看到您的图像有太多变化，因此需要注意的是，您必须确保您的图像具有很多高频信息才能使此方法起作用！