Performing a phase correlation with fft in R

Asked 3/6, 2015 at 21:1 Answered 9/6, 2015 at 4:53

Solved r image-processing fft cross-correlation

I am trying to implement a 2d phase correlation algorithm in R using a recipe from Wikipedia (http://en.wikipedia.org/wiki/Phase_correlation) in order to track the movement between 2 images. These images (frames) were captured with a camera shaking in the wind and the ultimate goal is to remove the shake in these and subsequent frames. The two example images and the R code are below:

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)

The vector shift, to my understanding, should contain information on the transitive shift (dx and dy) that best corrects these two images. However the shift variable gives dx=1 and dy=1, which I assume indicates no shift in either the x or y direction. This occurs for subsequent frames where there are visible shifts or several pixels in both the x and y direction.

Do any of y'all see an error in my code/formulas? Or do I need to try something fancier like filtering the images first before I do a phase correlation?

Cheers gals and guys!

Filly answered 3/6, 2015 at 21:1 Comment(1)

Hi I just came across this while researching for image stitching using phase correlation. May I know what's the reason behind using R and not simply Matlab or Python (if licensing costs is the main reason behind not using Matlab) – Candiscandle 13/9, 2017 at 14:58

The code looks correct from what I know about phase correlation. If I understand what you want correctly, you are trying to use phase correlation to determine the offset between two images given that their homographies are nothing more than horizontal and vertical offsets. The fact that you're only getting the shift to be at the origin is most likely due to your images lacking sufficient high frequency information in order to properly determine a good shift.

Try these two images instead (these were from the Wikipedia article you referenced, but I extracted them out and saved them as individual images):

Image #1 Image #2

When I run these two images with your R code, I get this for my phase correlation map. Bear in mind that your images were actually saved as .png, so I had to change the library to library(png) and I used readPNG instead of readTIFF. Keep that in mind when you try and run your code with the above example images:

enter image description here

Also, the location of where the maximum peak occurred was:

> shift
     row col
[1,] 132 153

This tells us that the image shifted over by 132 rows and 153 columns. Take note that this is with respect to the centre of the image. If you want to determine the actual offset, you'll need to subtract this by half the rows for the vertical coordinate and half the columns for the horizontal coordinate.

Therefore, the code works totally fine... it's just that your images lack sufficient high frequency information for the phase correlation to work. What correlation is trying to do in this case is that we're trying to find "similar" variations between each image. If there are a lot of variations between each image and are very similar, then phase correlation will work well. However, if we don't have that much variation, then phase correlation won't work.

Why is that the case? The basis behind phase correlation is that we assume that the image is corrupted with Gaussian white noise, and so if we correlate white noise with itself (from one image to another) it will give a very nice high peak at where the offset or the shift is and almost zero everywhere. Because of the fact that your images lack a lot of high frequency information and the fact that the images are clean, then phase correlation actually won't work. Therefore, what some people actually suggest is to pre-whiten your image so that the image contains white noise so that you can get the nice peak at where the offset should be that we're talking about.

However, just to make sure that you eliminate any false maximums, it is a good idea to also smoothen the cross-correlation matrix in the frequency domain (r in your R code) so that there is a high probability that there will only be one true maximum. Using a Gaussian filter in the frequency / FFT domain should work fine.

In any case, I don't see much variation in your images and so something to take away from this is that you gotta make sure your image has a lot of high frequency information for this to work!

Atropine answered 4/6, 2015 at 2:59 Comment(7)

@MikeWise - Thanks :)! I don't really post in the R tag, but this question was too tempting to pass up. – Atropine 4/6, 2015 at 3:14

I used to do a lot of FFT stuff years ago, and I looked at this last night, but I had no idea after reading it and not enough energy to try and refresh all that old knowledge. But it looks to me like you nailed it. – Incorporeity 4/6, 2015 at 3:17

Cheers for the advice, @Atropine ! I really appreciate your thorough answer. I am going to work on adding some white noise into the two images above and see what happens. Ill post any progress I make as I make it! – Filly 4/6, 2015 at 20:36

@rayreng I was able to stabilize the two frames I posted as well as the 300 subsequent with your suggestions/answer. So thanks! I added white noise to the images, which helped but I also found it useful to apply a Gaussian smoothing filter to the normalized cross correlation matrix (the matrix r in the code contained within my original question). That eliminated all unwanted maximums within the r matrix. Here is a link showing the stabilized video next to the shaky original. Thanks again for your help! (youtube.com/watch?v=irDFk2kbKaE) – Filly 5/6, 2015 at 15:25

@AlexMiller - Ah yes. I forgot to mention that smoothing will also help spread the peak of where the shift is detected... I should probably add that in just to be complete but that's truly awesome. Glad you got it to work! BTW, if you don't mind... would you be able to add the code that puts in these enhancements for self containment? I'm not that proficient in R and I'd like to see how you ended up doing all of that :) – Atropine 5/6, 2015 at 15:26

Nice one! The wiki article also mentions that you can extend this method to also extract possible differences in rotation or scaling by first converting the images to log-polar coordinates. Would you know how to do that by any chance? Article is B. S Reddy and B. N. Chatterji, “An FFT-based technique for translation, rotation, and scale-invariant image registration”, IEEE Transactions on Image Processing 5, no. 8 (1996): 1266–1271. – Delaine 24/12, 2015 at 9:50

@rayryeng-ReinstateMonica Hey, I'm trying to do a Python version of this code but I'm not getting the same result, can you look for me, please? pastebin.com/DG3Ed6Um – Culliton 2/4, 2020 at 23:33

Below is qualitative description of the routine followed by the R code used to efficiently and robustly find the phase correlation between the two images posted in the original question. Thanks to @rayreng for the advice and pointing me in the right direction!

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

The vector shift indicates that the image was shifted in the x positive direction and in negative y direction. In other words the second image (f2) was moved roughly to the upper right. The values of the vector shift will vary slightly with each attempt due to the noise introduced into the second image along with the smoothing operator from the Gaussian filter on the r matrix. I have noticed that a phase correlation similar to the one outlined above works better on larger images/matrices. To see the results of the above algorithm, visit the stabilized video at https://www.youtube.com/watch?v=irDFk2kbKaE.

Filly answered 9/6, 2015 at 4:53 Comment(1)

Hot tags

Godot Unity Godot Help Programming Godot 4.X GUI GDScript 3D 2D Physics CSharp Godot 3.X VR XR Projects C++

Recommended topics

Hot tags