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

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

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

# 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