daniel_chiu_Ex9

Autor: daniel_chiu
Data: 20/05/2009

Item 1

A imagem a ser utilizada e seu espectro encontram-se abaixo.

 1 IMG_FP=mmreadgray('db1_a/1_1.tif')
 2 DFT_FP=iadft(IMG_FP)
 3 
 4 iashow(IMG_FP,title='(1)Impressão digital')
 5 iashow(iadftview(DFT_FP),title='(2)DFT da figura (1)')
 6 
 7 IMG_BC=mmreadgray('Attachments/ia636-2009/exercicio9/barcode.tif')
 8 DFT_BC=iadft(IMG_BC)
 9 iashow(IMG_BC,title='(3)Código de barras')
10 iashow(iadftview(DFT_BC),title='(4)DFT da figura (3)')

(1)Impressão digital

(2)DFT da figura (1)

(3)Código de barras

(4)DFT da figura (3)

Filtro passa-baixa ideal

Um filtro passa-baixa ideal é um filtro onde todas as frequências dentro de uma figura como um círculo, quadrado ou retângulo passam sem atenuação e as frequências de fora são filtradas. Se a figura de um círculo é utilizada como filtro, a frequência de corte é dado pelo raio do círculo.

Onde, é dado por:

Para realizar a filtragem, deve-se deslocar a transformada da imagem ao centro de seu espectro. Para isso, a função iafftshift será utilizada. A partir daí, multiplica-se a transformada deslocada pelo filtro, que deve ser do mesmo tamanho da imagem.

Neste exemplo, a frequência de corte é 100.

1 P,Q=DFT_FP.shape
2 LOWPASS_AUX=iacircle([P,Q],100,[P/2,Q/2])
3 LOWPASS=ianormalize(LOWPASS_AUX,[0,1])
4 iashow(LOWPASS,title='(5)Filtro passa-baixa')
5 
6 DFT_LP_FP=DFT_FP*iaifftshift(LOWPASS)
7 iashow(iadftview(DFT_LP_FP), title='(6)Filtro multiplicado pelo espectro (2)')
8 IMG_LP_FP=abs(iaidft(DFT_LP_FP)).astype(uint16)
9 iashow(IMG_LP_FP,title='(7)Imagem da impressão filtrada')
Warning: downcasting image from double to uint16 (may lose precision)

(5)Filtro passa-baixa

(6)Filtro multiplicado pelo espectro (2)

(7)Imagem da impressão filtrada

1 P,Q=DFT_BC.shape
2 LOWPASS_AUX=iacircle([P,Q],100,[P/2,Q/2])
3 LOWPASS=ianormalize(LOWPASS_AUX,[0,1])
4 iashow(LOWPASS,title='(8)Filtro passa-baixa')
5 
6 DFT_LP_BC=DFT_BC*iaifftshift(LOWPASS)
7 iashow(iadftview(DFT_LP_BC), title='(9)Filtro multiplicado pelo espectro (4)')
8 IMG_LP_BC=abs(iaidft(DFT_LP_BC)).astype(uint16)
9 iashow(IMG_LP_BC,title='(10)Imagem das barras filtrada')
Warning: downcasting image from double to uint16 (may lose precision)

(8)Filtro passa-baixa

(9)Filtro multiplicado pelo espectro (4)

(10)Imagem das barras filtrada

Filtro passa-alta ideal

De maneira análoga ao filtro passa-baixa, o filtro passa-alta deixa passar as frequências do lado de fora do círculo e atenua as que estão do lado de dentro. A equação do filtro circular ideal é:

Onde, é dado por:

1 P,Q=DFT_FP.shape
2 HIGHPASS_AUX=iacircle([P,Q],100,[P/2,Q/2])
3 HIGHPASS=1-ianormalize(HIGHPASS_AUX,[0,1])
4 iashow(HIGHPASS,title='(14)Filtro passa-alta')
5 
6 DFT_HP_FP=DFT_FP*iaifftshift(HIGHPASS)
7 iashow(iadftview(DFT_HP_FP), title='(15)Filtro multiplicado pelo espectro (4)')
8 IMG_HP_FP=abs(iaidft(DFT_HP_FP)).astype(uint16)
9 iashow(IMG_HP_FP,title='(16)Imagem da impressão filtrada')
Warning: downcasting image from double to uint16 (may lose precision)

(14)Filtro passa-alta

(15)Filtro multiplicado pelo espectro (4)

(16)Imagem da impressão filtrada

1 P,Q=DFT_BC.shape
2 HIGHPASS_AUX=iacircle([P,Q],100,[P/2,Q/2])
3 HIGHPASS=1-ianormalize(HIGHPASS_AUX,[0,1])
4 iashow(HIGHPASS,title='(11)Filtro passa-alta')
5 
6 DFT_HP_BC=DFT_BC*iaifftshift(HIGHPASS)
7 iashow(iadftview(DFT_HP_BC), title='(12)Filtro multiplicado pelo espectro (2)')
8 IMG_HP_BC=abs(iaidft(DFT_HP_BC)).astype(uint16)
9 iashow(IMG_HP_BC,title='(13)Imagem da impressão filtrada')
Warning: downcasting image from double to uint16 (may lose precision)

(11)Filtro passa-alta

(12)Filtro multiplicado pelo espectro (2)

(13)Imagem da impressão filtrada

Item 2

No espectro da imagem original, os pontos estão distribuídos de maneira simétrica, como é possível ver na Tabela 1 impressa abaixo. Cada par de pontos simétrico no domínio da frquência corresponde a uma imagem cossenoidal. Dessa maneira, a imagem abaixo pode ser descrita através de uma composição de imagens cossenodais. Abaixo, podemos perceber isso ao reconstituir a imagem com algumas dessas componentes.

 1 def dwsc_dftdecompose():
 2 
 3     # Criação da imagem
 4     f=hstack(((ones((128,32))*200),ones((128,64))*50))
 5     f=hstack((f,ones((128,32))*200))
 6     f=f.astype(uint16)
 7 
 8     iashow(iapad(f),title='Imagem original')
 9 
10     # Cálculo da transformada de Fourier
11     F=iadft(f)
12     E=iadftview(F)
13     print 'Tabela 1\n',(F[0,:]).round(3)
14     iashow(E,title='Espectro da imagem')
15     mmplot([[E[E.shape[0]/2,:]]],ptitle='Grafico do espectro')
16 
17     # Decomposição da imagem em imagens cossenoidais
18     # Aqui se mostra a contribuição de algumas frequências (2^i-1 até 64)
19     F_recon=zeros_like(F)
20     for i in range(7):
21         F_decom=zeros_like(F)
22         F_decom[0,2**i-1]=F[0,2**i-1]
23         if i>0:
24             F_decom[0,128-2**i+1]=F[0,128-2**i+1]
25         F_recon+=F_decom
26         f_decom=(iaidft(F_decom).real).astype(uint16)
27         f_recon=(iaidft(F_recon).real).astype(uint16)
28         iashow(iapad(f_decom),title='Imagem componente frequência '+str(2**i-1))
29         iashow(iapad(f_recon),title='Imagem reconstruída')
30 
31 dwsc_dftdecompose()
Tabela 1
[  2.04800000e+06    +0.j   7.82121290e+05+19200.j  -0.00000000e+00    -0.j
  -2.60288049e+05-19200.j  -0.00000000e+00    +0.j   1.55669487e+05+19200.j
  -0.00000000e+00    +0.j  -1.10652326e+05-19200.j  -0.00000000e+00    -0.j
   8.55014830e+04+19200.j  -0.00000000e+00    +0.j  -6.93798850e+04-19200.j
  -0.00000000e+00    -0.j   5.81192300e+04+19200.j  -0.00000000e+00    -0.j
  -4.97741280e+04-19200.j  -0.00000000e+00    +0.j   4.33145060e+04+19200.j
   0.00000000e+00    -0.j  -3.81438490e+04-19200.j  -0.00000000e+00    +0.j
   3.38927400e+04+19200.j  -0.00000000e+00    -0.j  -3.03201130e+04-19200.j
   0.00000000e+00    +0.j   2.72619050e+04+19200.j  -0.00000000e+00    +0.j
  -2.46025260e+04-19200.j  -0.00000000e+00    +0.j   2.22581360e+04+19200.j
   0.00000000e+00    -0.j  -2.01663910e+04-19200.j   0.00000000e+00    +0.j
   1.82799200e+04+19200.j  -0.00000000e+00    -0.j  -1.65620340e+04-19200.j
  -0.00000000e+00    -0.j   1.49838270e+04+19200.j  -0.00000000e+00    -0.j
  -1.35221660e+04-19200.j   0.00000000e+00    +0.j   1.21582660e+04+19200.j
   0.00000000e+00    +0.j  -1.08766660e+04-19200.j  -0.00000000e+00    -0.j
   9.66446800e+03+19200.j   0.00000000e+00    +0.j  -8.51077500e+03-19200.j
  -0.00000000e+00    -0.j   7.40625700e+03+19200.j  -0.00000000e+00    +0.j
  -6.34282300e+03-19200.j   0.00000000e+00    -0.j   5.31335600e+03+19200.j
  -0.00000000e+00    -0.j  -4.31150400e+03-19200.j  -0.00000000e+00    -0.j
   3.33151600e+03+19200.j   0.00000000e+00    -0.j  -2.36809400e+03-19200.j
  -0.00000000e+00    +0.j   1.41627700e+03+19200.j  -0.00000000e+00    -0.j
  -4.71334000e+02-19200.j  -0.00000000e+00    +0.j  -4.71334000e+02+19200.j
  -0.00000000e+00    +0.j   1.41627700e+03-19200.j   0.00000000e+00    +0.j
  -2.36809400e+03+19200.j   0.00000000e+00    -0.j   3.33151600e+03-19200.j
   0.00000000e+00    -0.j  -4.31150400e+03+19200.j  -0.00000000e+00    +0.j
   5.31335600e+03-19200.j  -0.00000000e+00    -0.j  -6.34282300e+03+19200.j
  -0.00000000e+00    -0.j   7.40625700e+03-19200.j  -0.00000000e+00    +0.j
  -8.51077500e+03+19200.j   0.00000000e+00    -0.j   9.66446800e+03-19200.j
  -0.00000000e+00    +0.j  -1.08766660e+04+19200.j  -0.00000000e+00    -0.j
   1.21582660e+04-19200.j   0.00000000e+00    -0.j  -1.35221660e+04+19200.j
   0.00000000e+00    +0.j   1.49838270e+04-19200.j   0.00000000e+00    +0.j
  -1.65620340e+04+19200.j  -0.00000000e+00    -0.j   1.82799200e+04-19200.j
   0.00000000e+00    +0.j  -2.01663910e+04+19200.j   0.00000000e+00    -0.j
   2.22581360e+04-19200.j   0.00000000e+00    +0.j  -2.46025260e+04+19200.j
   0.00000000e+00    -0.j   2.72619050e+04-19200.j   0.00000000e+00    +0.j
  -3.03201130e+04+19200.j  -0.00000000e+00    -0.j   3.38927400e+04-19200.j
   0.00000000e+00    -0.j  -3.81438490e+04+19200.j  -0.00000000e+00    -0.j
   4.33145060e+04-19200.j   0.00000000e+00    -0.j  -4.97741280e+04+19200.j
   0.00000000e+00    -0.j   5.81192300e+04-19200.j   0.00000000e+00    +0.j
  -6.93798850e+04+19200.j   0.00000000e+00    -0.j   8.55014830e+04-19200.j
   0.00000000e+00    +0.j  -1.10652326e+05+19200.j   0.00000000e+00    -0.j
   1.55669487e+05-19200.j   0.00000000e+00    -0.j  -2.60288049e+05+19200.j
   0.00000000e+00    -0.j   7.82121290e+05-19200.j]

Imagem original

Espectro da imagem

Imagem componente frequência 0

Imagem reconstruída

Imagem componente frequência 1

Imagem reconstruída

Imagem componente frequência 3

Imagem reconstruída

Imagem componente frequência 7

Imagem reconstruída

Imagem componente frequência 15

Imagem reconstruída

Imagem componente frequência 31

Imagem reconstruída

Imagem componente frequência 63

Imagem reconstruída

Item 3 - Construção de uma função de filtragem

Baseando-se nas funções dos filtros abaixo, será construída uma função python que filtra a imagem, utilizando um filtro circular ideal, Butterworth ou Gaussiano.

Sendo: e a frequência de corte, as equações dos filtros são:

Filtro circular ideal:

  • Passa-baixa
  • Passa-alta

Filtro circular Butterworth de ordem :

  • Passa-baixa
  • Passa-alta

Filtro Gaussiano

  • Passa-baixa
  • Passa-alta

A implementação encontra-se abaixo. Clique em Change Parameters para mudar os parâmetros de entrada.

 1 def dwsc_circfilter(image,highpass_lowpass,cutoff_freq,type,Butterworth_order=1,show_filter=0,show_spectrum=0):
 2 
 3     N,M=image.shape
 4 
 5     u,v=iameshgrid(arange(M),arange(N))
 6     D0=cutoff_freq
 7     D=zeros((N,M))
 8 
 9     D=sqrt((u-M/2)**2+(v-N/2)**2)
10 
11     if type==0: #ideal
12         caption='Filtro ideal'
13         if highpass_lowpass==0: #passa-baixa
14             caption+=' passa-baixa'
15             H=(D<=D0)
16         else: #passa-alta
17             caption+=' passa-alta'
18             H=(D>D0)
19         H=H.astype(float)
20 
21 
22     elif type==1: #Butterworth
23         caption='Filtro Butterworth'
24         if highpass_lowpass==0: #passa-baixa
25             caption+=' passa-baixa'
26             H=1.0/(1+(D/D0)**(2*n))
27         else: #passa-alta
28             caption+=' passa-alta'
29             H=1.0/(1+(D0/D)**(2*n))
30 
31 
32     elif type==2: #Gaussiano
33         caption='Filtro Gaussiano'
34         if highpass_lowpass==0: #passa-baixa
35             caption+=' passa-baixa'
36             H=exp(-D**2/(2*D0**2))
37         else: #passa-alta
38             caption+=' passa-alta'
39             H=1-exp(-D**2/(2*D0**2))
40 
41     if show_filter!=0:
42         adshow(iapad(ianormalize(H,[0,255])),title=caption)
43 
44     F=iadft(image)
45     G=F*iaifftshift(H)
46 
47     if show_spectrum!=0:
48         adshow(iadftview(F),title='Espectro da imagem original')
49         adshow(iadftview(G),title='Espectro da imagem filtrada')
50 
51     filtered_image=abs(iaidft(G)).astype(uint16)
52 
53     return filtered_image
54 
55 parser = OptionParser()
56 parser.add_option("--imagem", type='string',\
57 default='Attachments/ia636-2009/exercicio9/barcode.tif',\
58 help='copie e cole o link de uma das imagens da lista em http://parati.dca.fee.unicamp.br/adesso/wiki/ia636-2009/images/view/')
59 parser.add_option("--tipo_corte", type='int', default=0, help='Passa-baixa=0, Passa-alta=1')
60 parser.add_option("--freq_corte", type='int', default=100, help='Frequência de corte')
61 parser.add_option("--tipo_filtro", type='int', default=0, help='Ideal=0, Butterworth=1, Gaussiano=2')
62 parser.add_option("--ordem_butterworth", type='int', default=2, help='Ordem do filtro Butterworth')
63 parser.add_option("--mostrar_filtros", type='int',default=1, help='Mostra filtros e espectros')
64 
65 opt, arg = parser.parse_args()
66 
67 f=mmreadgray(opt.imagem)
68 hi_lo_pass=opt.tipo_corte
69 cutoff_freq=opt.freq_corte
70 type=opt.tipo_filtro
71 n=opt.ordem_butterworth
72 show_filter=opt.mostrar_filtros
73 
74 g=dwsc_circfilter(f,hi_lo_pass,cutoff_freq,type,n,show_filter,show_filter)
75 
76 adshow(f, title='Imagem original')
77 adshow(g, title='Imagem filtrada')

Imagem original

Imagem filtrada