Phase correlation with polar coordinates

Definition: Using phase correlation the angle can be easily found out.

Steps:
  • Both images must be the same size.
  • FFT from the images.
  • Polar transform.
  • FFT from the polar spectrum.
  • Phase coorelation from the polar spectrum. (angle)
  • Inverse FFT
 1 def ictransformacaopolar(f,m,n):
 2 
 3     m,n = f.shape
 4     Rx = floor(m/2)
 5     Ry = floor(n/2)
 6 
 7     f[Rx,:] = (f[Rx-1,:]+f[Rx+1,:])/2
 8     f[:,Ry] = (f[:,Ry-1]+f[:,Ry+1])/2
 9 
10     b = min(Rx,Ry)/m
11     a = pi/n
12 
13     xl,yl = indices((m,n))
14 
15     XI = Rx + (b*yl)*cos(a*xl)
16     YI = Ry + (b*yl)*sin(a*xl)
17 
18     g = iainterpolclosest(f, (XI,YI))
19 
20     return g
 1 from ia636 import *
 2 from numpy import *
 3 
 4 def phase_correlation_p(f,g): #both images have to be the same
 5     F = fft.fftn(f.astype(float))
 6     G = fft.fftn(f.astype(float))
 7 
 8     pF = ictransformacaopolar(F,F.shape[0],F.shape[1])
 9     pG = ictransformacaopolar(G,G.shape[0],G.shape[1])
10 
11     pFF = fft.fftn(pF)
12     pGG = fft.fftn(pG)
13 
14     crossproduct = (pFF*pGG.conj())/(abs(pFF*pGG.conj()))
15     peak = fft.ifftn(crossproduct).real
16 
17     #p, pos = max(ravel(peak)), argmax(ravel(peak))
18     #(row, col) = iaind2sub(peak.shape, pos)
19 
20     return peak #=angle
 1 from numpy import *
 2 
 3 f = mmreadgray("courseIA369O1S2011/jan_14_0/fprint.jpg")
 4 adshow(f, title ='reference image')
 5 T = array([[1,0,0],[0,1,0],[0,0,1]], 'd')
 6 angle = math.radians(5)
 7 Rotation = array([[cos(angle),-sin(angle),0],[sin(angle),cos(angle),0],[0,0,1]],'d')
 8 T = dot(T,Rotation)
 9 g = iaffine(f,T)
10 adshow(g, title ='sensing image')
11 
12 #Polar coordinates
13 angle = phase_correlation_p(f,g)
14 adshow(ianormalize(angle), title ='phase correlation')
15 adshow(ianormalize(angle[0:10]), title ='phase correlation')
16 adshow(ianormalize(angle[::40,::40]), title ='phase correlation')
17 #print "peaks:", angle
18 maxvalue = unravel_index(argmax(angle), angle.shape)
19 ans = mod((float(maxvalue[1])/angle.shape[1])*2*pi,pi)
20 print "angle:", ans
angle: 0.0

reference image

sensing image

phase correlation

phase correlation

phase correlation