Demo iahotelling

Illustrate the Hotelling Transform

Script

Image creation

A binary image with an ellipse

1. import numpy as np
2. import numpy.linalg as LA
3. import ia636 as ia
4. 
5. f = ia.iagaussian((100,100), array([[45,50]]).T, [[40*40,10*10],[10*10,20*20]])
6. f = f > (f.max() * np.exp(-1./2))
7. adshow(f, title='Ellipse')

Ellipse

Feature vector

The coordinates of each 1 pixel are the features: cols and rows coordinates.

1. X = np.nonzero(f)
2. X = np.array(X).T
3. n,m = X.shape
4. print 'Number of samples (n):', n
5. print 'Number of attributes (m):', m
6. print 'X=\n',
7. print X
Number of samples (n): 2471
Number of attributes (m): 2
X=
[[ 6 44]
 [ 6 45]
 [ 6 46]
 ..., 
 [84 54]
 [84 55]
 [84 56]]

Mean vector

The mean value is the centroid.

1. mu = X.mean(axis=0)
2. print 'Centroid = ',mu
Centroid =  [ 45.  50.]

Computing the covariance matrix

The covariance matrix is given by:

It can be computed by the numpy function cov or directly from the equation:

1. C = np.cov(X-mu, rowvar=False)
2. C1 = ((X-mu).T).dot(X-mu)/(n-1)
3. print 'Covariance matrix computed usint np.cov\n',C.round(2)
4. print 'Covariance matrix computed from the definition\n',C1.round(2)
Covariance matrix computed usint np.cov
[[ 398.01   24.4 ]
 [  24.4    98.72]]
Covariance matrix computed from the definition
[[ 398.01   24.4 ]
 [  24.4    98.72]]

Computing eigenvalues and eigenvectors

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

01. [e,V] = LA.eigh(C)
02. print 'e:',e.round(2)
03. print 'V:\n',V.round(2)
04. index = np.argsort(e)[::-1]
05. V = V[:,index]
06. e = e[index]
07. print '\nEigenvalues vector\n'
08. print e.round(2)
09. print '\nEigenvectors matrix\n'
10. print V.round(2)
11. print '\nInverse is transpose\n',
12. print LA.inv(V).round(2)
e: [  96.74  399.98]
V:
[[ 0.08 -1.  ]
 [-1.   -0.08]]

Eigenvalues vector

[ 399.98   96.74]

Eigenvectors matrix

[[-1.    0.08]
 [-0.08 -1.  ]]

Inverse is transpose
[[-1.   -0.08]
 [ 0.08 -1.  ]]

Measure the angle of inclination

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

1. #sigma = np.sqrt(e)
2. theta = arctan(V[0,0]/V[1,0])*180./np.pi
3. print 'The inclination angle is', theta.round(3),' degrees'
The inclination angle is 85.37  degrees

Plot the eigenvectors and centroid

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

1. adshow(f, title='Ellipse')

Ellipse

Compute the Hotelling transform

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

01. Xnew = (X-mu).dot(V)
02. Xnew_mu = Xnew.mean(axis=0)
03. print 'Mean value of the transform\n'
04. print Xnew_mu.round(2)
05. C_Xnew = np.cov(Xnew,rowvar=False)
06. print '\nCovariance matrix of the transform\n'
07. print C_Xnew.round(2)
08. C_Xnew = np.where(C_Xnew <1e-10, 0, C_Xnew)
09. print '\nStandard deviation matrix of the transformed data\n'
10. print np.sqrt(C_Xnew).round(2)
Mean value of the transform

[ 0.  0.]

Covariance matrix of the transform

[[ 399.98    0.  ]
 [   0.     96.74]]

Standard deviation matrix of the transformed data

[[ 20.     0.  ]
 [  0.     9.84]]

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.

The black dots in the transformed image is due to the direct geometric transform.

1. Xnew = np.rint(Xnew + mu).astype(int)
2. g = np.zeros(f.shape, np.bool)
3. g[Xnew[:,1],Xnew[:,0]] = True
4. adshow(g, title='Ellipse w/ straightened axes')

Ellipse w/ straightened axes

Image read and display

The RGB color image is read and displayed

1. f = adread('boat.ppm')
2. adshow(f,title='Original RGB image')

Original RGB image

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. adshow(r,title='Red component')
5. adshow(g,title='Green component')
6. adshow(b,title='Blue component')

Red component

Green component

Blue component

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

01. X = np.array([r.ravel(),g.ravel(),b.ravel()]).T
02. print X.shape
03. V,e,mu = ia.iapca(X)
04. print 'Mean value of each channel (R,G,B)\n'
05. print mu.round(2)
06. Cx = np.cov(X,rowvar=False)
07. print '\nCovariance matrix\n'
08. print Cx.round(2)
09. print '\nEigenvalues matrix\n'
10. print e.round(2)
11. print '\nEigenvectors matrix\n'
12. print V.round(2)
(65792, 3)
Mean value of each channel (R,G,B)

[ 101.27   94.92   87.41]

Covariance matrix

[[ 3497.95  3273.2   3086.74]
 [ 3273.2   3239.67  3103.11]
 [ 3086.74  3103.11  3032.4 ]]

Eigenvalues matrix

[ 9572.67   180.23    17.12]

Eigenvectors matrix

[[-0.6  -0.77  0.24]
 [-0.58  0.21 -0.79]
 [-0.56  0.61  0.57]]

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

01. y = (X-mu).dot(V)
02. my = mean(y,axis=0)
03. print 'Mean value of the transform\n'
04. print my.round(3)
05. Cy = np.cov(y,rowvar=False)
06. print '\nCovariance matrix of the transform\n'
07. print Cy.round(2)
08. Cy = np.where(Cy <1e-10, 0, Cy)
09. print '\nStandard deviation matrix of the transformed data\n'
10. print np.sqrt(Cy).round(2)
Mean value of the transform

[ 0.  0.  0.]

Covariance matrix of the transform

[[ 9572.67    -0.      -0.  ]
 [   -0.     180.23    -0.  ]
 [   -0.      -0.      17.12]]

Standard deviation matrix of the transformed data

[[ 97.84   0.     0.  ]
 [  0.    13.42   0.  ]
 [  0.     0.     4.14]]

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].reshape(r.shape)
2. g2 = y[:,1].reshape(r.shape)
3. g3 = y[:,2].reshape(r.shape)
4. adshow(255-ia.ianormalize(g1),title='Variance:%f' % (Cy[0,0],))
5. adshow(ia.ianormalize(g2),title='Variance:%f' % (Cy[1,1],))
6. adshow(ia.ianormalize(g3),title='Variance:%f' % (Cy[2,2],))

Variance:9572.669654

Variance:180.229077

Variance:17.119634

Full reconstruction

One can reconstruct the image by multiplying the transformed image by the transpose of eigenvector matrix and adding the mean value of the image.

1. y1=y.copy()
2. X1 = y1.dot(V.T) + mu
3. f1=np.zeros_like(f)
4. for i in range(3):
5.     f1[i,:,:]=reshape(X1[:,i],r.shape)
6. print 'Mean square error = ',np.mean((f1-f)**2).round(2)
7. print 'Max square error = ',np.max((f1-f)**2).round(2)
8. adshow(f,title='original')
9. adshow(f1,title='Recomposed image')
Mean square error =  0.08
Max square error =  1

original

Recomposed image

Partial reconstruction

The component corresponding to the smallest variance is erased from the transformed image. It is noticeable that the mean square error is the smallest eigenvalue.

01. y1=y.copy()
02. y1[:,2]=0
03. X1 = y1.dot(V.T) + mu
04. print 'iaimginfo(X1):', ia.iaimginfo(X1)
05. X1n = ia.ianormalize(X1).astype('uint8')
06. f1=np.zeros_like(f)
07. for i in range(3):
08.     f1[i,:,:]=reshape(X1n[:,i],r.shape)
09. print 'Mean square error = ',np.mean((1.*f1-f)**2).round(2)
10. print 'Max square error = ',np.max((1.*f1-f)**2).round(2)
11. adshow(f,title='original')
12. adshow(f1,title='Recomposed image')
13. adshow(ia.ianormalize(np.abs(1*f1-f)),'normalized absolute error')
iaimginfo(X1): <type 'numpy.ndarray'> (65792, 3) float64 -5.213513 217.375115
Mean square error =  445.43
Max square error =  1764.0

original

Recomposed image

normalized absolute error

References

See Also