Demo iahotelling

Illustrate the Hotelling Transform

Script

Image creation

A binary image with an ellipsis (demo)

1 #from ia636 import iagaussian
2 
3 f = gigaussian([100,100], [45,50], [[20,5],[5,10]]) > 0.0000001
4 iashow(f)
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 3, in <module>
  File "<string>", line 26, in gigaussian
NameError: global name 'iaind2sub' is not defined

------------------------------------------------------------

Feature vector

The coordinates of each 1 pixel are the features: cols and rows coordinates. The mean value is the centroid.

 1 #from ia636 import iaind2sub
 2 
 3 #import numpy.oldnumeric.mlab as MLab
 4 
 5 (row, col) = iaind2sub(f.shape, nonzero(ravel(f)))
 6 rows,cols = row[0,:],col[0,:]
 7 
 8 x = zeros([rows.shape[0],2])
 9 x[:,0],x[:,1] = cols,rows
10 
11 mx = mean(x,0)
12 print mx
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 5, in <module>
NameError: name 'iaind2sub' is not defined

------------------------------------------------------------

Computing eigenvalues and eigenvectors

The eigenvalues and eigenvectors are computed from the covariance. The eigenvalues are sorted in decrescent order

1 Cx = cov(transpose(x))
2 print"\n covariance: \n",Cx
3 [aval, avec] = eig(Cx)
4 aux = argsort(aval)[::-1]
5 aval = diag(take(aval, aux))
6 print"\n Eigenvalues: \n",aval
7 avec = take(avec, aux, axis=1)
8 print"\n Eigenvectors: \n",avec
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 1, in <module>
NameError: name 'x' is not defined

------------------------------------------------------------

Measure the angle of inclination

The direction of the eigenvector of the largest eigenvalue gives the inclination of the elongated figure

1 a1 = sqrt(aval[0,0])
2 a2 = sqrt(aval[1,1])
3 vec1 = avec[:,0]
4 vec2 = avec[:,1]
5 theta = arctan(1.*vec1[1]/vec1[0])*180/pi
6 print 'angle is %3f degrees' % (theta)
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 1, in <module>
NameError: name 'aval' is not defined

------------------------------------------------------------

Plot the eigenvectors and centroid

The eigenvectors are placed at the centroid and scaled by the square root of its correspondent eigenvalue

 1 #from ia636 import iacontour
 2 
 3 iashow(f)
 4 x_,y_ = iaind2sub(f.shape, nonzero(ravel(iacontour(f))))
 5 x_,y_ = x_[0,:],y_[0,:]
 6 mmplot([[y_, 100-x_]])
 7 
 8 # esboça os autovetores na imagem
 9 x1, y1 = arange(mx[0], 100-mx[0]+a1*vec1[0], -0.1), arange(mx[1], mx[1]-a1*vec1[1],0.235)
10 mmplot([[x1, 100-y1]]) # largest in green
11 
12 x2, y2 = arange(mx[0], mx[0]-a2*vec2[0], 0.1), arange(mx[1], mx[1]+a2*vec2[1],0.042)
13 mmplot([[x2, 100-y2]]) # smaller in blue
14 
15 mmplot([[[mx[0], 100-mx[1]]]])# centroid in magenta
16 #g('set data style points')
17 #g('set xrange [0:100]')
18 #g('set yrange [0:100]')
19 # g.plot(d0, d1, d2, d3)
20 # showfig(mx)
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 3, in <module>
NameError: name 'f' is not defined

------------------------------------------------------------

Compute the Hotelling transform

The Hotelling transform, also called Karhunen-Loeve (K-L) transform, or the method of principal components, is computed below

1 mx = resize(mx,[x.shape[0],2])
2 y = transpose(dot(avec, transpose(x-mx)))
3 my = mean(y)
4 print my
5 Cy = cov(y)
6 print "\n",Cy
7 print "\n",floor(0.5 + sqrt(Cy))
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 1, in <module>
NameError: name 'mx' is not defined

------------------------------------------------------------

Display the transformed data

The centroid of the transformed data is zero (0,0). To visualize it as an image, the features are translated by the centroid of the original data, so that only the rotation effect of the Hotelling transform is visualized.

1 #from ia636 import iasub2ind
2 
3 ytrans = floor(0.5 + y + resize(mx, (x.shape[0], 2)))
4 g = zeros(f.shape)
5 i = iasub2ind(f.shape, ytrans[:,1], ytrans[:,0])
6 put(g, i, 1)
7 iashow(g)
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 3, in <module>
NameError: name 'y' is not defined

------------------------------------------------------------

Image read and display

The RGB color image is read and displayed

1 f = iaread('boat.ppm')
2 adshow(f)

Extracting and display RGB components

The color components are stored in the third dimension of the image array

1 r = f[0,:,:]
2 g = f[1,:,:]
3 b = f[2,:,:]
4 iashow(r)
5 iashow(g)
6 iashow(b)

Feature vector: R, G and B values

The features are the red, green and blue components. The mean vector is the average color in the image. The eigenvalues and eigenvectors are computed. The dimension of the covariance matrix is 3x3 as there are 3 features in use

 1 x = 1.*concatenate((ravel(r)[:,newaxis], ravel(g)[:,newaxis], ravel(b)[:,newaxis]), 1)
 2 mx = mean(x,0)
 3 print mx
 4 Cx = cov(x,rowvar=0)
 5 print"\n", Cx
 6 [aval, avec] = eig(Cx)
 7 aux = argsort(aval)[::-1]
 8 aval = diag(take(aval, aux))
 9 print"\n", aval
10 avec = take(avec, aux,1)
11 print"\n", avec
[ 101.26770732   94.9191999    87.41316573]

[[ 3497.95108557  3273.20413327  3086.73663437]
 [ 3273.20413327  3239.67139072  3103.11181452]
 [ 3086.73663437  3103.11181452  3032.39588874]]

[[ 9572.66965395     0.             0.        ]
 [    0.           180.2290775      0.        ]
 [    0.             0.            17.11963358]]

[[-0.59515464 -0.76785348  0.2370485 ]
 [-0.58010086  0.20637147 -0.78796815]
 [-0.55612404  0.60647494  0.5682554 ]]

Hotelling transform

The K-L transform is computed. The mean vector is zero and the covariance matrix is decorrelated. We can see the values of the standard deviation of the first, second and third components

1 mx1 = tile(mx,(x.shape[0],1))
2 
3 y = transpose(dot(avec, transpose(x-mx1)))
4 my = mean(y,0)
5 print "\n",my
6 Cy = cov(y,rowvar=0)
7 print "\n",Cy
8 Cy = Cy * (1 - all(Cy) < 1e-10)
9 print "\n",floor(0.5 + sqrt(Cy))
[ -5.02408318e-16   1.09747058e-13  -1.29345410e-13]

[[ 4310.55813695  4236.27947486 -2084.89696396]
 [ 4236.27947486  4226.86030824 -2135.00580178]
 [-2084.89696396 -2135.00580178  1232.59991984]]

[[ 66.  65.  NaN]
 [ 65.  65.  NaN]
 [ NaN  NaN  35.]]

Displaying each component of the transformed image

The transformed features are put back in three different images with the g1 with the first component (larger variance) and g3 with the smaller variance

1 g1 = y[:,0]
2 g2 = y[:,1]
3 g3 = y[:,2]
4 g1 = reshape(g1, r.shape)
5 g2 = reshape(g2, r.shape)
6 g3 = reshape(g3, r.shape)
7 iashow(g1)
8 iashow(g2)
9 iashow(g3)
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)