Demo iaphasecorr

namespace:ia636
page: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")

Image is cartesian coordinates

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

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

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

Contributions

  • André Luis da Costa, 1st semester 2011