Demo Phase Correlation

Illustrate using Phase Correlation to estimate rotation and translation between images.

Description

In this lesson we explain how to use Phase Correlation to estimate the angle of rotation and the translation between 2D images.

Converting an image from Cartesian to Polar coordinates

It converts a plane from coordinates to , with and . Notice that the domain in polar coordinates must be informed explicitaly and will influence in the angle resolution.

1. from ia636 import *
2. 
3. f = ianormalize(iagaussian((151,151), [[75],[75]], [[800,0],[0,800]]), [0,200]).astype(uint8)
4. 
5. adshow(iaisolines(f,10,3), "Image in cartesian coordinates")
6. g = iapolar(f,(150,300),pi)
7. adshow(iaisolines(g.astype(int),10,3), "Image in polar coordinates")
8. #adshow(g, "Image in polar coordinates")

Image in cartesian coordinates

Image in polar coordinates

Estimating the angle of rotation

The following function will be used to estimate the angle of rotation between 2D images.

1. def iarotphasecorr2d(f,h):
2.     F = fft.fftn(f)
3.     H = fft.fftn(h)
4.     pF = iapolar(iadftview(F),(F.shape[0]/2,360),pi)
5.     pH = iapolar(iadftview(H),(H.shape[0]/2,360),pi)
6.     return iaphasecorr(pF, pH)

The function can be applied as follows.

01. f = adread("cameraman.pgm")
02. t = zeros(array(f.shape)+200)
03. t[100:f.shape[0]+100,100:f.shape[1]+100] = f
04. f = t
05. 
06. t1 = array([
07.             [1,0,-f.shape[0]/2.],
08.             [0,1,-f.shape[1]/2.],
09.             [0,0,1]]);
10. 
11. t2 = array([
12.             [1,0,f.shape[0]/2.],
13.             [0,1,f.shape[1]/2.],
14.             [0,0,1]]);
15. 
16. theta = radians(30)
17. r1 = array([
18.         [cos(theta),-sin(theta),0],
19.         [sin(theta),cos(theta),0],
20.         [0,0,1]]);
21. 
22. T = dot(t2,dot(r1,t1))
23. h = iaffine(f,T,0)
24. h.shape = f.shape
25. 
26. adshow(f, "Original image")
27. adshow(h, "Image rotated by 30°")
28. 
29. # Execute the function with the image and its rotated sample
30. pg = iarotphasecorr2d(f,h)
31. pg = ianormalize(pg)
32. 
33. # Get the point of higher value
34. peak = unravel_index(argmax(pg), pg.shape)
35. 
36. # Calculate the angle
37. ang = (float(peak[1])/pg.shape[1])*180
38. 
39. adshow(pg, "Correlation map in polar coordinates - Peak at %s, corresponds to angle %f°"%(str(peak),ang))

Original image

Image rotated by 30°

Correlation map in polar coordinates - Peak at (0, 60), corresponds to angle 30.000000°

Estimating the angle of rotation and the translation

Now we will compute the angle of rotation and the translation. The function below first find the angle of rotation; after that, it rotate the image and find the translation. Two phase correlation maps are returned: one for the translation and other for rotation.

1. import scipy
2. 
3. def iatrphasecorr2d(f,h):
4.     rg = iarotphasecorr2d(f,h)
5.     peak = unravel_index(argmax(rg), rg.shape)
6.     ang = (float(peak[1])/rg.shape[1])*180
7.     h_rot = scipy.ndimage.interpolation.rotate(h, -ang, reshape=False)
8.     g = iaphasecorr(f,h_rot)
9.     return g, rg

The following code find the angle of rotation and the translation. Then, the original image is obtained from the rotated and translated sample image.

01. t3 = array([
02.             [1,0,50],
03.             [0,1,32],
04.             [0,0,1]]);
05. 
06. T = dot(t3,T)
07. h = iaffine(f,T,0)
08. h.shape = f.shape
09. 
10. adshow(f, "Original image")
11. adshow(h, "Image rotated by 30° and translated by (50,32)")
12. 
13. g, rg = iatrphasecorr2d(f,h)
14. g = ianormalize(g)
15. rg = ianormalize(rg)
16. 
17. trans_peak = unravel_index(argmax(g), g.shape)
18. rot_peak = unravel_index(argmax(rg), rg.shape)
19. ang = (float(rot_peak[1])/rg.shape[1])*180
20. trans = (array(h.shape)-array(trans_peak))
21. 
22. adshow(g, "Translation correlation map - Peak %s, \n corresponds to translation %s"%(str(trans_peak), str(tuple(trans))))
23. adshow(ianormalize(rg), "Rotation correlation map - Peak %s, corresponds to angle %f°"%(str(rot_peak),ang))
24. 
25. t4 = array([
26.             [1,0,-trans[0]],
27.             [0,1,-trans[1]],
28.             [0,0,1]]);
29. 
30. 
31. theta1 = radians(-ang)
32. r2 = array([
33.         [cos(theta1),-sin(theta1),0],
34.         [sin(theta1),cos(theta1),0],
35.         [0,0,1]]);
36. 
37. T1 = dot(t4,dot(t2,dot(r2,t1)))
38. f1 = iaffine(h,T1,0)
39. f1.shape = h.shape
40. adshow(f1, "Sample image rotated and translated by %f° and %s, respectively"%(-ang,tuple(-trans)))

Original image

Image rotated by 30° and translated by (50,32)

Translation correlation map - Peak (397, 454), corresponds to translation (59, 2)

Rotation correlation map - Peak (0, 60), corresponds to angle 30.000000°

Sample image rotated and translated by -30.000000° and (-59, -2), respectively

See Also

References

  1. B. Srinivasa Reddy and B. N. Chatterji. An FFT-Based Technique for Translation, Rotation, and Scale-Invariant Image Registration. IEEE Trans. on Image Processing. vol 5. n. 8, 1996. PDF at IEEE

Contributions

  • André Luis da Costa, 1st semester 2011