lionisEX13

Autor: lionis
Data: 02/07/2009

Exercício 1

1. Corrigir a função de geração da gaussiana que está na toolbox do Adessowiki

 1 from numpy import *
 2 from numpy.linalg import *
 3 from ia636 import iaind2sub
 4 
 5 def iagaussian2(s, mu, sigma):
 6     s, mu, sigma = asarray(s), asarray(mu), asarray(sigma)
 7 
 8     if (s.size == 1):
 9         print "Caso 1D"
10         x = arange(s)
11         d = x - mu
12         k = (1. * d) / sigma
13         k *= k
14     else:
15         print "Caso ",len(s),"D"
16 
17         num_elem = prod(s)
18         num_dim = len(s)
19 
20         x = zeros((num_elem, num_dim))
21         x[:,num_dim-1] = arange(num_elem)
22         #print "x="; print x
23 
24         for i in range(1,num_dim):
25             aux=x[:,num_dim-i]
26             x[:,num_dim-i-1],x[:,num_dim-i]=iaind2sub([s[num_dim-i-1],s[num_dim-i]],aux)
27         x=x.astype(uint16)
28 
29         d = x - resize(reshape(mu,(mu.size,1)), (num_elem,mu.size))
30 
31         if sigma.size == 1:
32             tmp = 1./sigma
33         else:
34             tmp = inv(sigma)
35 
36         k = dot(d*d, tmp)
37         k = sum(k,axis=1)
38 
39     g = zeros(s, float)
40     aux = ravel(g)
41     if sigma.size == 1:
42         tmp = sigma
43     else:
44         tmp = det(sigma)
45     aux[:] = 1./(sqrt(2*pi)*tmp) * exp(-1./2 * k)
46     return g

Caso 1D

1 f = iagaussian2(100, 50, 10)
2 g = ianormalize(f, [0,1])
3 mmplot([[g]])
Caso 1D

Caso n-D

1 f = iagaussian2([100,80], [50,30], [[20*20,0],[0,30*30]])
2 g = ianormalize(f, [0,255]).astype(uint16)
3 iashow(g, title='Gaussiana 2D')
Caso  2 D

Gaussiana 2D

Exercício 2

2. Corrigir transformada de Hoteling da toolbox do Adessowiki

Corrigir a demonstração da transformada de Hoteling da toolbox do Adessowiki e complementar com a reconstrução da imagem colorida do barquinho, retirando a última componente da imagem aplicada a Trans. de Hoteling e verifica se o erro médio quadrático é de fato o valor do menor autovetor.

Image creation

1 f = iagaussian2([100,100], [45,60], [[1,0],[0,1]]) > 0.0000001
2 
3 g = ianormalize(f, [0,255]).astype(uint16)
4 iashow(f, title='Gaussiana 2D')
Caso  2 D

Gaussiana 2D

Feature vector

 1 from ia636 import iaind2sub
 2 import numpy.oldnumeric.mlab as MLab
 3 
 4 #f = iagaussian2([8,6], [4,2], [[1,0],[0,1]]) > 0.0000001
 5 
 6 non_zero = asarray(nonzero(ravel(f)))
 7 rows, cols = iaind2sub(f.shape, non_zero)
 8 x = concatenate((rows, cols), 0)
 9 mx = x.mean(axis=1)
10 
11 print "Centroid=",mx
Centroid= [ 45.  60.]

Computing eigenvalues and eigenvectors

1 Cx = cov(x)
2 print "Covariance Matrix="; print Cx.round()
3 [aval, avec] = eig(Cx)
4 aux = argsort(aval)[::-1]
5 print "aux="; print aux
6 aval = diag(take(aval, aux))
7 print "Eigenvalues="; print aval.round()
8 avec = take(avec, aux, axis=1)
9 print "Eigenvectors="; print avec.round()
Covariance Matrix=
[[ 8.  0.]
 [ 0.  8.]]
aux=
[1 0]
Eigenvalues=
[[ 8.  0.]
 [ 0.  8.]]
Eigenvectors=
[[ 0.  1.]
 [ 1.  0.]]

Measure the angle of inclination

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 ",theta," degrees"
angle is  90.0  degrees

Plot the eigenvectors and centroid

 1 iashow(f, title="Elipse original")
 2 
 3 non_zero = asarray(nonzero(ravel(f)))
 4 rows, cols = iaind2sub(f.shape, non_zero)
 5 mmplot([[rows, 100-cols]])
 6 
 7 # esboça os autovetores na imagem
 8 x1, y1 = arange(mx[0], 100-mx[0]+a1*vec1[0], 0.1), arange(mx[1], 100-mx[1]-a1*vec1[1], -0.1)
 9 mmplot([[x1, 100-y1]]) # largest in green
10 
11 x2, y2 = arange(mx[0], mx[0]-a2*vec2[0], 0.1), arange(100-mx[1], mx[1]+a2*vec2[1],0.042)
12 
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 16, in <module>
TypeError: 'numpy.ndarray' object is not callable

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

Elipse original

Compute the Hotelling transform

 1 xt = x-mx
 2 print "xt="; print xt
 3 print transpose(xt)
 4 
 5 y = transpose(dot(avec, transpose(x-mx)))
 6 my = y.mean(axis=0)
 7 print my
 8 Cy = cov(y)
 9 print Cy
10 print floor(0.5 + sqrt(Cy))
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 1, in <module>
ValueError: shape mismatch: objects cannot be broadcast to a single shape

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

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 iashow(f)
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 2, in <module>
  File "xs_runner_python.py", line 138, in mmshow
    else:                im = mmgshow(f)
TbxError: mm_gshow: can only process 1 slice gray-scale 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 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 = MLab.mean(x)
 3 print mx
 4 Cx = MLab.cov(x)
 5 print Cx
 6 [aval, avec] = MLab.eig(Cx)
 7 aux = argsort(aval)[::-1]
 8 aval = MLab.diag(take(aval, aux))
 9 print aval
10 avec = take(avec, aux)
11 print 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 ]

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 y = transpose(dot(avec, transpose(x-mx)))
2 my = MLab.mean(y)
3 print my
4 Cy = MLab.cov(y)
5 print Cy
6 Cy = Cy * (1 - (-1e-10 < Cy < 1e-10))
7 print floor(0.5 + sqrt(Cy))
-5.02408318147e-16
4310.55813695
66.0

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

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 1, in <module>
IndexError: invalid index

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