Lista de Alunos | Ex.13 | Ex.15

Ex14

  1. Implementar a correlação de fase 2D para translação E rotação. A implementação do Pedro pode servir como ponto de partida, mas é necessário torná-la mais genérica.
  2. Ler e entender cada passo do algoritmo apresentado no artigo: Volume Registration using the 3-D pseudo-polar Fourier transform

Programas Pedro

ICInterpolacao2

 1 def icinterpolacao2(f, XI, YI):
 2 
 3     m,n = XI.shape
 4 
 5     g = zeros((m,n))
 6     for i in arange(m):
 7         for j in arange(n):
 8 
 9             if XI[i,j]==floor(XI[i,j]):
10                 px1 = 1
11                 px2 = 0
12             else:
13                 px1 = XI[i,j]-floor(XI[i,j])
14                 px2 = ceil(XI[i,j])-XI[i,j]
15 
16             if YI[i,j]==floor(YI[i,j]):
17                 py1 = 1
18                 py2 = 0
19             else:
20                 py1 = YI[i,j]-floor(YI[i,j])
21                 py2 = ceil(YI[i,j])-YI[i,j]
22 
23             p1 = (px1*py1)*f[floor(YI[i,j]),floor(XI[i,j])]
24             p2 = (px1*py2)*f[ceil(YI[i,j]),floor(XI[i,j])]
25             p3 = (px2*py1)*f[floor(YI[i,j]),ceil(XI[i,j])]
26             p4 = (px2*py2)*f[ceil(YI[i,j]),ceil(XI[i,j])]
27             g[i,j] = round(p1+p2+p3+p4)
28 
29     return g

ICTransformaçãoPolar

 1 def ictransformacaopolar(f, MP, NP):
 2 
 3     m,n = f.shape
 4     R = floor(m/2)
 5 
 6     f[R,:] = (f[R-1,:]+f[R+1,:])/2
 7     f[:,R] = (f[:,R-1]+f[:,R+1])/2
 8 
 9     b = R/MP
10     a = pi/NP
11 
12     xl,yl = meshgrid(arange(NP),arange(MP))
13 
14     XI = R + (b*yl)*cos(a*xl)
15     YI = R + (b*yl)*sin(a*xl)
16 
17     g = icinterpolacao2(f, XI, YI)
18     g = g[50:200,:]
19 
20     return g, a

Execução do programa:

 1 im = mmreadgray('db1_a/1_1.tif')
 2 f = 255*ones((389,389))
 3 f[0:374,0:388] = im
 4 mmshow(f)
 5 
 6 MP=500
 7 NP=500
 8 
 9 g,a = ictransformacaopolar(f, MP, NP)
10 mmshow(g)
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

ICPicos

 1 def icpicos(z):
 2 
 3     from ia636 import *
 4 
 5     lin,col = z.shape
 6 
 7     x1 = z-take(z,(arange(col)-1)%col,1)
 8     x2 = z-take(z,(arange(col)+1)%col,1)
 9     x3 = z-take(z,(arange(lin)-1)%lin,0)
10     x4 = z-take(z,(arange(lin)+1)%lin,0)
11 
12     mask = (x1>0)&(x2>0)&(x3>0)&(x4>0)
13 
14     picos = z[mask]
15     picos = sort(picos)
16 
17     peaks = zeros((5,2))
18 
19     for i in arange(5):
20         teste = ones(z.shape)
21         teste2 = cumsum(teste)-1
22         teste2.shape = z.shape
23         teste3 = (z == picos[i])
24         ind = teste2[teste3]
25         y0,x0 = iaind2sub(z.shape, ind)
26         peaks[i,0] = floor(x0)
27         peaks[i,1] = floor(y0)
28 
29     return peaks
<string>:1: SyntaxWarning: import * only allowed at module level

Execução do programa:

1 z = array([[1,0,0,0,0],[0,2,0,0,0],[0,0,3,0,0],[0,0,0,4,0],[0,0,0,0,5],])
2 peaks = icpicos(z)
3 print peaks
[[ 0.  0.]
 [ 1.  1.]
 [ 2.  2.]
 [ 3.  3.]
 [ 4.  4.]]

ICCorrelaçãoDeFase1Nova

Realiza a correlação de fase entre as imagens do espectro de Fourier de cada uma das imagens para determinar os 5 ângulos candidatos.

 1 def iccorrelacaodefase1nova(f1, f2):
 2 
 3     from numpy.fft import *
 4     from ia636 import *
 5 
 6     m1,n1 = f1.shape
 7     m2,n2 = f2.shape
 8 
 9     m = max(m1,m2) - max(m1,m2)%2 + 1;
10     n = max(n1,n2) - max(n1,n2)%2 + 1;
11 
12     ff1 = zeros((m,n))
13     ff2 = zeros((m,n))
14     ff1[0:m1,0:n1] = f1
15     ff2[0:m2,0:n2] = f2
16 
17     FF1 = fft2(ff1)
18     FF2 = fft2(ff2)
19 
20     FF = FF1*conj(FF2)
21     NCS = FF/abs(FF)
22 
23     PC = real(ifft2(NCS))
24     PC = iapconv(PC,[[0.16,0.4,0.16],[0.4,1,0.4],[0.16,0.4,0.16]])
25 
26     peaks = icpicos(PC)
27     for i in arange(5):
28         x0 = peaks[i,0]
29         y0 = peaks[i,1]
30         x0 = -(x0-1)
31         y0 = -(y0-1)
32         if x0<-(n/2):
33             x0 = x0 + n
34         elif x0>n/2:
35             x0 = x0 - n
36         if y0<-(m/2):
37             y0 = y0 + m
38         elif y0>(m/2):
39             y0 = y0 - m
40         peaks[i,0] = x0
41         peaks[i,1] = y0
42 
43     return peaks
<string>:1: SyntaxWarning: import * only allowed at module level
<string>:1: SyntaxWarning: import * only allowed at module level

Execução do programa:

1 f1 = mmreadgray('db1_a/1_1.tif')
2 f2 = mmreadgray('db1_a/1_2.tif')
3 peaks = iccorrelacaodefase1nova(f1, f2)
4 print peaks
[[  52.   79.]
 [  48.  167.]
 [ 162.    9.]
 [  57.  -96.]
 [ -40.  -36.]]

ICCorrelaçãoDeFase2Nova

Realiza a correlação de fase entre as imagens das impressões digitais, para cada um dos 5 ângulos candidatos.

 1 def iccorrelacaodefase2nova(f1, f2):
 2 
 3     from numpy.fft import *
 4     from ia636 import *
 5 
 6     m1,n1 = f1.shape
 7     m2,n2 = f2.shape
 8 
 9     m = max(m1,m2) - max(m1,m2)%2 + 1;
10     n = max(n1,n2) - max(n1,n2)%2 + 1;
11 
12     ff1 = 255*ones((m,n),dtype=uint8)
13     ff2 = 255*ones((m,n),dtype=uint8)
14     ff1[0:m1,0:n1] = f1
15     ff2[0:m1,0:n1] = f2
16 
17     FF1 = fft2(ff1)
18     FF2 = fft2(ff2)
19 
20     FF = FF1*conj(FF2)
21     NCS = FF/abs(FF)
22 
23     PC = real(ifft2(NCS))
24     mmshow(ianormalize(PC,[0,255]))
25     PC = iapconv(PC,[[0.16,0.4,0.16],[0.4,1,0.4],[0.16,0.4,0.16]])
26     mmshow(ianormalize(PC,[0,255]))
27     PC2 = sort(ravel(PC))
28 
29     pcmax = PC2[-1]
30 
31     return pcmax
<string>:1: SyntaxWarning: import * only allowed at module level
<string>:1: SyntaxWarning: import * only allowed at module level

Execução do programa:

1 f1 = mmreadgray('db1_a/1_1.tif')
2 f2 = mmreadgray('db1_a/1_2.tif')
3 pcmax = iccorrelacaodefase2nova(f1, f2)
4 print pcmax
Warning: downcasting image from int32 to uint16 (may lose precision)
Warning: downcasting image from int32 to uint16 (may lose precision)
0.144156443338

ICProjeto

Realiza a comparação entre duas impressões digitais e retorna a correlação de fase entre elas.

 1 from ia636 import *
 2 
 3 def icprojeto(nome1,nome2):
 4 
 5     from numpy.fft import *
 6 
 7     im1 = mmreadgray(nome1)
 8     im2 = mmreadgray(nome2)
 9     m1,n1 = im1.shape
10     m2,n2 = im2.shape
11 
12     dim = max(max(m1,n1),max(m2,n2));
13     if dim%2==0:
14         dim += 1
15 
16     im1mod = 255*ones((dim,dim),dtype=uint8)
17     im2mod = 255*ones((dim,dim),dtype=uint8)
18     im1mod[0:m1,0:n1] = im1
19     im2mod[0:m2,0:n2] = im2
20 
21     IM1 = abs(fftshift(fft2(im1mod)))
22     IM2 = abs(fftshift(fft2(im2mod)))
23 
24     mmshow(im1mod)
25     mmshow(im2mod)
26     mmshow(ianormalize(log(abs(IM1)+1),[0,255]))
27     mmshow(ianormalize(log(abs(IM2)+1),[0,255]))
28 
29     IM1t,a = ictransformacaopolar(IM1,500,500)
30     IM2t,a = ictransformacaopolar(IM2,500,500)
31     picos = iccorrelacaodefase1nova(IM1t,IM2t)
32 
33     mmshow(ianormalize(IM1t,[0,255]))
34     mmshow(ianormalize(IM2t,[0,255]))
35     print picos
36 
37     mc = zeros((1,5))
38 
39     maxc = 0
40     return maxc
<string>:3: SyntaxWarning: import * only allowed at module level

Execução do programa:

1 nome1 = 'db1_a/1_1.tif'
2 nome2 = 'db1_a/1_2.tif'
3 
4 maxc = icprojeto(nome1,nome2)
5 print maxc
ERROR execute

*** Exception while evaluating code:
Sandbox thread timeout [30 s]

Parte restante do programa:

for i in arange(5):

x1 = picos[i,0]

theta = (-x1*a)%(2*pi)

theta1 = floor(2*theta/pi)

theta2 = theta%(pi/2)

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

im2r = icrotacao(im2,theta1);

im2r = icgeorigidbranco(im2r,theta2);

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

mc[i] = iccorrelacaodefase2nova(im1,im2r)

maxc = mc.max

Ver o exercício do Fernando ia636-2009:CRoNuXs-EX4 de transformação geométrica.