# Demo iaphasecorr

namespace: ia636 iaphasecorrdemo

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

This function is required by the process of estimating the angle of rotation. 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 numpy import *
2 from ia636 import *
3
4 def iatopolar(f,domain,ang_range=2*pi):
5     f = array(f)
6     m,n = f.shape
7     dm,dn = domain
8     Ry,Rx = floor(array(f.shape)/2)
9
10     b = min(Ry,Rx)/dm
11     a = ang_range/dn
12
13     y,x = indices(domain)
14
15     XI = Rx + (b*y)*cos(a*x)
16     YI = Ry + (b*y)*sin(a*x)
17
18     g = iainterpollin(f, array([YI.ravel(), XI.ravel()]))
19     g.shape = domain
20     return g
```

## Testing iatopolar() function

```1 #import scipy
2
3 f = ianormalize(iagaussian((151,151), [[75],[75]], [[800,0],[0,800]]), [0,200])
4
5 adshow(f, "Image is cartesian coordinates")
6 g = iatopolar(f,(150,300))
7 adshow(g, "Image is 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 = iatopolar(iadftview(F),(F.shape[0]/2,360), pi)
5     pH = iatopolar(iadftview(H),(H.shape[0]/2,360), pi)
6     return iaphasecorr(pF, pH)
```

The function can be applied as follows.

``` 1 f = adread("cameraman.pgm")
2 t = zeros(array(f.shape)+200)
3 t[100:f.shape[0]+100,100:f.shape[1]+100] = f
4 f = t
5
6 t1 = array([
7             [1,0,-f.shape[0]/2.],
8             [0,1,-f.shape[1]/2.],
9             [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 def iatrphasecorr2d(f,h):
2     rg = iarotphasecorr2d(f,h)
3     peak = unravel_index(argmax(rg), rg.shape)
4     ang = (float(peak[1])/rg.shape[1])*180
5     h_rot = scipy.ndimage.interpolation.rotate(h, -ang, reshape=False)
6     g = iaphasecorr(f,h_rot)
7     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.

``` 1 t3 = array([
2             [1,0,50],
3             [0,1,32],
4             [0,0,1]]);
5
6 T = dot(t3,T)
7 h = iaffine(f,T,0)
8 h.shape = f.shape
9
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(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)))
```