Introdução
O filtro dito como homomórfico é baseado nos princípios de iluminância e reflectância - conceitos já vistos antes na disciplina - para a realização da filtragem. Matematicamente falando, podemos mostrar a seguinte relação: f(x,y) = i(x,y)r(x,y) Do qual, o "i" (iluminância) representa as variações espaciais lentas, ditas de baixas frequências, e o "r"(reflectância) representa as variações espaciais rápidas, ditas de altas frequências.
A iluminância representa a quantidade de luz incidente sobre a imagem/pixel. Já a reflectância demonstra quanto dessa luz é refletida e depende da superficie pela qual a luz incide.
O filtro homomórfico promete corrigir a má iluminação da cena. Para que isso seja bem articulado, devemos atuar separadamente nas componentes de iluminância e reflectância, para que possamos atenuar as frequencias baixas (i) e manter as frequências altas (r).
Código no OpenCV
O código implementado pode ser visto a seguir:
Listagem 1. equalize.cpp
#include <iostream>
#include <opencv2/opencv.hpp>
#include <math.h>
#define RADIUS 1
using namespace cv;
using namespace std;
// troca os quadrantes da imagem da DFT
void deslocaDFT(Mat& image ){
Mat tmp, A, B, C, D;
// se a imagem tiver tamanho impar, recorta a regiao para
// evitar cópias de tamanho desigual
image = image(Rect(0, 0, image.cols & -2, image.rows & -2));
int cx = image.cols/2;
int cy = image.rows/2;
// reorganiza os quadrantes da transformada
// A B -> D C
// C D B A
A = image(Rect(0, 0, cx, cy));
B = image(Rect(cx, 0, cx, cy));
C = image(Rect(0, cy, cx, cy));
D = image(Rect(cx, cy, cx, cy));
// A <-> D
A.copyTo(tmp); D.copyTo(A); tmp.copyTo(D);
// C <-> B
C.copyTo(tmp); B.copyTo(C); tmp.copyTo(B);
}
int main(int argc, char** argv){
Mat imaginaryInput, complexImage, multsp;
Mat padded, filter, mag;
Mat image, imagegray, tmp;
Mat_<float> realInput, zeros;
vector<Mat> planos;
// habilita/desabilita ruido
int noise=0;
// frequencia do ruido
int freq=10;
// ganho inicial do ruido
float gain=1;
// valor do ruido
float mean;
// guarda tecla capturada
char key;
// valores ideais dos tamanhos da imagem
// para calculo da DFT
int dft_M, dft_N;
// abre a câmera default
// cap.open(0);
// if(!cap.isOpened())
// return -1;
// captura uma imagem para recuperar as
// informacoes de gravação
//cap >> image;
image = imread(argv[1],CV_LOAD_IMAGE_GRAYSCALE);
// identifica os tamanhos otimos para
// calculo do FFT
dft_M = getOptimalDFTSize(image.rows);
dft_N = getOptimalDFTSize(image.cols);
// realiza o padding da imagem
copyMakeBorder(image, padded, 0,
dft_M - image.rows, 0,
dft_N - image.cols,
BORDER_CONSTANT, Scalar::all(0));
// parte imaginaria da matriz complexa (preenchida com zeros)
zeros = Mat_<float>::zeros(padded.size());
// prepara a matriz complexa para ser preenchida
complexImage = Mat(padded.size(), CV_32FC2, Scalar(0));
// a função de transferência (filtro frequencial) deve ter o
// mesmo tamanho e tipo da matriz complexa
filter = complexImage.clone();
// cria uma matriz temporária para criar as componentes real
// e imaginaria do filtro ideal
tmp = Mat(dft_M, dft_N, CV_32F);
// prepara o filtro passa-baixas ideal
for(int i=0; i<dft_M; i++){
for(int j=0; j<dft_N; j++){
if((i-dft_M/2)*(i-dft_M/2)+(j-dft_N/2)*(j-dft_N/2) < RADIUS*RADIUS){
tmp.at<float> (i,j) = 1.0;
}
}
}
float GAMMA_H = 1;
float GAMMA_L = 0.25;
float c = 1;
float D, D0;
D0 = RADIUS*RADIUS;
// prepara o filtro homomórfico
for(int i=0; i<dft_M; i++){
for(int j=0; j<dft_N; j++){
D = (i-dft_M/2)*(i-dft_M/2)+(j-dft_N/2)*(j-dft_N/2);
tmp.at<float> (i,j) = (GAMMA_H - GAMMA_L)*(1.0-exp(-1.0*c*((D*D)/((D0*D0))))) + GAMMA_L;
}
}
// cria a matriz com as componentes do filtro e junta
// ambas em uma matriz multicanal complexa
Mat comps[]= {tmp, tmp};
merge(comps, 2, filter);
imagegray = imread(argv[1],CV_LOAD_IMAGE_GRAYSCALE);
imshow("original", imagegray);
// realiza o padding da imagem
copyMakeBorder(imagegray, padded, 0,
dft_M - image.rows, 0,
dft_N - image.cols,
BORDER_CONSTANT, Scalar::all(0));
// limpa o array de matrizes que vao compor a
// imagem complexa
planos.clear();
// cria a compoente real
realInput = Mat_<float>(padded);
// insere as duas componentes no array de matrizes
planos.push_back(realInput);
planos.push_back(zeros);
// combina o array de matrizes em uma unica
// componente complexa
merge(planos, complexImage);
// calcula o dft
dft(complexImage, complexImage);
// realiza a troca de quadrantes
deslocaDFT(complexImage);
// aplica o filtro frequencial
mulSpectrums(complexImage,filter,complexImage,0);
// limpa o array de planos
planos.clear();
// separa as partes real e imaginaria para modifica-las
split(complexImage, planos);
// usa o valor medio do espectro para dosar o ruido
mean = abs(planos[0].at<float> (dft_M/2,dft_N/2));
// insere ruido coerente, se habilitado
if(noise){
// F(u,v) recebe ganho proporcional a F(0,0)
planos[0].at<float>(dft_M/2 +freq, dft_N/2 +freq) +=
gain*mean;
planos[1].at<float>(dft_M/2 +freq, dft_N/2 +freq) +=
gain*mean;
// F*(-u,-v) = F(u,v)
planos[0].at<float>(dft_M/2 -freq, dft_N/2 -freq) =
planos[0].at<float>(dft_M/2 +freq, dft_N/2 +freq);
planos[1].at<float>(dft_M/2 -freq, dft_N/2 -freq) =
-planos[1].at<float>(dft_M/2 +freq, dft_N/2 +freq);
}
// recompoe os planos em uma unica matriz complexa
merge(planos, complexImage);
// troca novamente os quadrantes
deslocaDFT(complexImage);
cout << complexImage.size().height << endl;
// calcula a DFT inversa
idft(complexImage, complexImage);
// limpa o array de planos
planos.clear();
// separa as partes real e imaginaria da
// imagem filtrada
split(complexImage, planos);
// normaliza a parte real para exibicao
normalize(planos[0], planos[0], 0, 1, CV_MINMAX);
imshow("filtrada", planos[0]);
for(;;){
key = (char) waitKey(10);
if( key == 27 ) break; // esc pressed!
switch(key){
// aumenta a frequencia do ruido
case 'q':
freq=freq+1;
if(freq > dft_M/2-1)
freq = dft_M/2-1;
break;
// diminui a frequencia do ruido
case 'a':
freq=freq-1;
if(freq < 1)
freq = 1;
break;
// amplifica o ruido
case 'x':
gain += 0.1;
break;
// atenua o ruido
case 'z':
gain -= 0.1;
if(gain < 0)
gain=0;
break;
// insere/remove ruido
case 'e':
noise=!noise;
break;
}
}
imwrite("dark_out.jpeg", planos[0]);
return 0;
}Resultados
Podemos ver o resultado da aplicação do filtro a seguir:
Figura 1. Entrada do fitro homomórfico
Figura 2. Saída do filtro homomórfico.