Exercício 09 - Filtragem no domínio da frequência

Autor: Wagner Machado do Amaral
Data: 18/05/2009

Enunciado: ia636-2009:exercicio9

Índice

  • iacos - iacos com seus parâmetros calculados no domínio da frequência
  • filtros - Filtragem no domínio da frequência

iacos com seus parâmetros calculados no domínio da frequência

A função iacos pode ser implementada com seus parâmetros calculados no domínio da frequência, conforme abaixo:

 1 m = 200. # colunas
 2 n = 100. # linhas
 3 T = 200. # frequencia da onda
 4 
 5 F = zeros([n,m],dtype=complex) # espectro com amplitude zeros em todas as frequências
 6 
 7 # atribui 1 para cossenoides horizontais de frequencia
 8 # determinada pelo parametro freq no sentido negativo
 9 v = 0
10 u = round(m/T)
11 F[v, u] = complex( 1, 0)
12 F[v, m-u] = F[v, u]
13 
14 print 'Simetrico? ', iaisdftsym(F)
15 
16 f = abs(iaidft(F))
17 
18 mmshow(ianormalize(f,[0,256]))
Simetrico?  True
Warning: downcasting image from double to uint16 (may lose precision)

 1 m = 200 # colunas
 2 n = 100 # linhas
 3 freq = 2 # frequencia da onda
 4 
 5 F = zeros([n,m],dtype=complex) # espectro com amplitude zeros em todas as frequências
 6 
 7 # atribui 1 para cossenoides horizontais de frequencia
 8 # determinada pelo parametro freq no sentido negativo
 9 v = (n/2)
10 u = (m/2 - freq)
11 F[v, u] = complex( 1, 0)
12 
13 # atribui 1 para cossenoides horizontais de frequencia
14 # determinada pelo parametro freq no sentido positivo
15 v = (n/2)
16 u = (m/2 + freq)
17 F[v, u] = complex( 1, 0)
18 
19 mmshow(ianormalize(abs(F),[0,256]),title="Espectro gerado manualmente.")
20 mmshow(ianormalize(abs(iadft(F)),[0,256]),title="Inversa da Fourier aplicada no espectro.")
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Espectro gerado manualmente.

Inversa da Fourier aplicada no espectro.

Teoricamente, se aplicarmos uma rotação no espectro o resultado da transformada inversa de Fourier também será rotacionado.

Sendo P(u,v) um ponto no espectro, o ponto rotacionado P(u,v)' pode ser determinado por:

Dessa forma, u' e v' podem ser determinados por:

 1 m = 256 # colunas
 2 n = 256 # linhas
 3 freq = 2 # frequencia da onda
 4 t = pi/4 # angulo theta de rotação da senoide
 5 
 6 F = zeros([n,m],dtype=complex)
 7 
 8 v = (n/2)
 9 u = (m/2 - freq)
10 uR = u*cos(t)+v*-sin(t)
11 vR = u*sin(t)+v*cos(t)
12 F[ uR , vR ] = complex( 1, 0);
13 
14 v = (n/2)
15 u = (m/2 + freq)
16 uR = u*cos(t)+v*-sin(t)
17 vR = u*sin(t)+v*cos(t)
18 F[ uR , vR ] = complex( 1, 0);
19 
20 mmshow(ianormalize(abs(F),[0,256]),title="Espectro gerado manualmente.")
21 mmshow(ianormalize(abs(iadft(F)),[0,256]),title="Inversa da Fourier aplicada no espectro.")
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Espectro gerado manualmente.

Inversa da Fourier aplicada no espectro.

Por algum motivo o espectro não está como esperado.

Porém a Inversa da Fourier aplicada no espectro gerou uma onda senoidal rotacionada conforme esperado.

Filtragem no domínio da frequência

Segue abaixo a implementação de algumas funções utilizadas para realizar Filtragem no domínio da frequência:

Função recostruir_imagem

Aplica filtro no domínio da frequencia de uma imagem.

1 def recostruir_imagem(f,H):
2   F = iadft(f) #transformada de fourier da imagem f
3   FH = F * iaifftshift(H) #aplica filtro no espectro
4   g=abs(iaidft(FH)) #reconstroe imagem com espectro alterado pelo filtro
5 
6   iashow(ianormalize(H,[0,255]),title='Filtro')
7   iashow(iadftview(FH),title='Filtro combinado com espectro')
8 
9   return g

Função filtro_circulo

Cria filtro em formato de circulo.

Sendo:

Passa baixa:

Passa alta:

1 def filtro_circulo(size,D0,tipo_filtro):
2   N, M = size[0], size[1]
3   (u, v) = iameshgrid(range(M), range(N))
4   D = sqrt( (u-M/2)**2 + (v-N/2)**2 )
5   if( tipo_filtro == 'passabaixa' ):
6     H = D<=D0
7   elif( tipo_filtro == 'passaalta' ):
8     H = D>D0
9   return H

Função filtro_butterworth

Cria filtro butterworth.

Sendo:

Passa baixa:

Passa alta:

 1 def filtro_butterworth(size,D0,n,tipo_filtro):
 2   N, M = size[0], size[1]
 3   (u, v) = iameshgrid(range(M), range(N))
 4   D = sqrt((u-M/2)**2 + (v-N/2)**2)
 5 
 6   if( tipo_filtro == 'passabaixa' ):
 7     H = 1 / ( 1 + (D/D0)**2*n )
 8   elif( tipo_filtro == 'passaalta' ):
 9     H = 1 / ( 1 + (D0/D)**2*n )
10   return H

Função filtro_passa_faixa

Cria filtro passa faixa em formato de circulo.

1 def filtro_passa_faixa(size,frequencia_corte):
2   C_maior = filtro_circulo(size,frequencia_corte[1],'passabaixa')
3   C_menor = filtro_circulo(size,frequencia_corte[0],'passaalta')
4   H = C_maior * C_menor
5   return H

Função filtro_retangulo

Cria filtro em formato de retangulo.

 1 def filtro_retangulo(size,frequencia_corte,tipo_filtro):
 2   N, M = size[0], size[1]
 3   if ( frequencia_corte*2 > M ):
 4     rect_width = M
 5   else:
 6     rect_width = frequencia_corte*2
 7 
 8   if ( frequencia_corte*2 > N ):
 9     rect_height = N
10   else:
11     rect_height = frequencia_corte*2
12 
13   if( tipo_filtro == 'passabaixa' ):
14     H = iarectangle([N,M], [rect_height,rect_width], [N/2.,M/2.])
15   elif( tipo_filtro == 'passaalta' ):
16     H = -iarectangle([N,M], [rect_height,rect_width], [N/2.,M/2.])
17 
18   return H

Testando as funções de filtro no dominio da frequencia

1 imgname = 'db1_a/1_1.tif'
2 f =adreadgray(imgname)
3 iashow(f,title='Imagem original: '+imgname)
4 F = iadft(f)
5 iashow(iadftview(F),title='Espectro da Transformada de Fourier')

Imagem original: db1_a/1_1.tif

Espectro da Transformada de Fourier

1 frequencia_corte = 50
2 tipo_filtro = 'passabaixa'
3 H = filtro_circulo(f.shape,frequencia_corte,tipo_filtro)
4 g = recostruir_imagem(f,H)
5 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro

1 frequencia_corte = 100
2 tipo_filtro = 'passaalta'
3 H = filtro_circulo(f.shape,frequencia_corte,tipo_filtro)
4 g = recostruir_imagem(f,H)
5 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro

1 frequencia_corte = 50
2 tipo_filtro = 'passabaixa'
3 H = filtro_retangulo(f.shape,frequencia_corte,tipo_filtro)
4 g = recostruir_imagem(f,H)
5 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro

1 frequencia_corte = 50
2 tipo_filtro = 'passabaixa'
3 n=1
4 H = filtro_butterworth(f.shape,frequencia_corte,n,tipo_filtro)
5 g = recostruir_imagem(f,H)
6 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro

1 frequencia_corte = 50
2 tipo_filtro = 'passabaixa'
3 n=100
4 H = filtro_butterworth(f.shape,frequencia_corte,n,tipo_filtro)
5 g = recostruir_imagem(f,H)
6 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro

1 frequencia_corte = 50
2 tipo_filtro = 'passaalta'
3 n=1
4 H = filtro_butterworth(f.shape,frequencia_corte,n,tipo_filtro)
5 g = recostruir_imagem(f,H)
6 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro

1 frequencia_corte = [20,50]
2 H = filtro_passa_faixa(f.shape,frequencia_corte)
3 g = recostruir_imagem(f,H)
4 iashow( g , title='Imagem reconstruída pelo novo espectro' )
Warning: downcasting image from double to uint16 (may lose precision)
Warning: downcasting image from double to uint16 (may lose precision)

Filtro

Filtro combinado com espectro

Imagem reconstruída pelo novo espectro