1. Objetivo
Vimos como tratar imagens no domínio espacial. Agora devemos aprender a tratar imagens no domínio da frequência, nos aprofundando no que é esse domínio, suas diferenças com o domínio espacial e as vantages de utilizar filtros de frequência.
2. Fourier e o Domínio da Frequência
Usualmente, estamos acostumados a trabalhar com o domínio do espaço e tempo. É o mais natural para nós seres humanos, pois é o que estamos habituados em nosso cotidiano. Porém o matemático e físico francês Jean-Baptiste Joseph Fourier nos introduziu suas Série e Transformada, que nos permitem representar sinais no domínio da frequência. Isso significa que podemos ver como os sinais se repetem. A série de Fourier nos permite representar qualquer sinal periódico como uma soma de senóides, e a transformada nos permite levar essas funções para o domínio da frequência. A transformada de Fourier é dada por:
\$F(u,v) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)*e^{-j2\pi(ux + vy)} dx dy\$
E a transformada inversa é dada por:
\$f(x,y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(u,v)*e^{j2\pi(ux + vy)} du dv\$
Quando lidamos com computadores, devemos utilizar a transformada discreta de Fourier, pois os sinais são amostrados pelos equipamentos eletrônicos de forma discreta.
\$F(u,v) = \frac{1}{MN} \sum_{0}^{M-1} \sum_{0}^{N-1} f(x,y)*e^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})}\$
A transformada pode ser algo extremamente custoso computacionalmente, devido à grande quantidade de exponeciais. Porém existe o algoritmo da Transformada Rápida de Fourier, que funciona bem com imagens que possuem dimensões múltiplas de 2, 3 ou 5. A fft tem um ganho significativo de desempenho em relação à transformada dft, e por isso iremos utilizá-la. Para isso, aumentaremos o tamanho da imagem com uma técnica chamada de padding.
| O padding consiste em copiar a imagem para uma imagem maior, deixando espaços vazios (0) nas bordas. Isso não altera a transformada da imagem. |
Alguns componentes da imagem podem ser separados em componentes de alta frequência e de baixa frequência. Os sinais de baixa frequência, quando tratamos de imagens, representam as cores e forma gerais da imagem. Já as componentes de alta frequência representam os detalhes da imagem. Isso será mais bem compreendido no tópico seguinte.
2.1. Trocando Regiões
Para facilitar nossos cálculos, devemos deslocar o espectro de fourier para o centro da figura. Isso não altera a transformada devido à sua propriedade de periodicidade. Para isso, devemos multiplicar todos os seus elementos por um termo. Utilizando as propriedades da transformada de Fourier, esse cálculo é substituído por uma troca nas regiões imagem.
\$\mathfrak{F}(f(x,y)(-1)^(x+y)) = F(u-M/2,v-N/2)\$
3. Visualização do Espectro de Frequência e Filtros
Iremos aplicar a transformada de Fourier à figura 2. O seu espectro de amplitude é mostrado logo em seguida, na figura 3. O código utilizado pode ser baixado aqui ou pode ser visto de forma semelhante na documentação do OpenCV.
Para utilizarmos bem o domínio da frequência, precisamos entender os tipos de filtros que são utilizados. Os filtros passa-alta deixam passar as frequências mais altas (figura 4). Os filtros passa-baixa deixam passar as frequências mais baixas (figura 5). Temos também os filtros passa-faixa, nos quais apenas as componentes em uma determinada faixa de frequência são mantidas.
4. O Filtro Homomórfico
Cada pixel de uma imagem pode ser separado em componentes distintas de iluminância e reflectância. A iluminância (\$i(x,y)\$) representa a quantidade de luz que incide sobre o pixel, e apresenta variações lentas (baixa frequência). Já a reflectância (\$i(x,y)\$) indica quanto dessa luz incidente é refletida, dependendo do material e apresentando variações mais rápidas (alta frequência). Desse modo, cada pixel pode ser descrito por:
\$f(x,y) = i(x,y)r(x,y)\$
O filtro homomórfico tem o objetivo de corrigir a má iluminação de uma cena. Para isso, devemos atuar separadamente nas componentes de iluminância e reflectância, e assim atenuar a iluminância. Como \$\mathfrak{F}(i(x,y)r(x,y)) ne \mathfrak{F}(i(x,y))\mathfrak{F}(i(x,y))\$, devemos executar as seguintes operações para a aplicação do filtro:
\$z(x,y) = ln(i(x,y)r(x,y)) = ln(i(x,y)) + ln(r(x,y))\$
\$S(u,v) = \mathfrak{F}(z(x,y))*H(u,v)\$
\$s(x,y) = \mathfrak{F}^-1(S(u,v))\$
E a imagem filtrada, \$g(x,y)\$ é dada por:
\$g(x,y) = e^(s(x,y))\$
A componente \$H(u,v)\$ é o filtro homomórfico na frequência. Ele pode ser descrito por:
\$H(u,v) = (\gamma_H - \gamma_L)(1-e^(-c((D^2(u,v))/D_0^2))) + \gamma_L\$
5. O Código
A implementação do filtro Homomófico é mostrada abaixo, na listagem homomorfico.cpp
#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/imgproc.hpp>
using namespace cv;
using namespace std;
void deslocaDFT(Mat& image ){
Mat tmp, A, B, C, D;
image = image(Rect(0, 0, image.cols & -2, image.rows & -2));
int cx = image.cols/2;
int cy = image.rows/2;
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.copyTo(tmp); D.copyTo(A); tmp.copyTo(D);
C.copyTo(tmp); B.copyTo(C); tmp.copyTo(B);
}
int gammaL_slider = 2, gammaH_slider = 20, sharpC_slider = 1, cutoff_slider = 5;
const int gammaL_max = 10, gammaH_max = 50, sharpC_max = 100, cutoff_max = 200;
int gammaL, gammaH, sharpC, cutoff;
Mat im, imFiltered,padded;
int dft_M, dft_N;
Mat homomorphicFilter(double gl, double gh, double c, double d0){
Mat filter = Mat(padded.size(), CV_32FC2, Scalar(0));
Mat tmp = Mat(dft_M, dft_N, CV_32F);
for(int i=0; i<dft_M; i++){
for(int j=0; j<dft_N; j++){
tmp.at<float> (i,j) = (gh - gl)*(1 - exp(-c*(( (i-dft_M/2)*(i-dft_M/2) + (j-dft_N/2)*(j-dft_N/2) ) / (d0*d0) ))) + gl;
}
}
Mat comps[]= {tmp,tmp};
imshow("Filter", tmp);
merge(comps, 2, filter);
// normalize(filter,filter,0,1,CV_MINMAX);
return filter;
}
void applyFilter(void){
vector<Mat> planos; planos.clear();
Mat zeros = Mat_<float>::zeros(padded.size());
Mat realInput = Mat_<float>(padded);
Mat complex;
realInput += Scalar::all(1);
log(realInput,realInput);
//normalize(realInput, realInput, 0, 1, CV_MINMAX);
//imshow("logimage",realInput);
planos.push_back(realInput);
planos.push_back(zeros);
merge(planos, complex);
dft(complex, complex);
deslocaDFT(complex);
resize(complex,complex,padded.size());
normalize(complex,complex,0,1,CV_MINMAX);
Mat filter = homomorphicFilter(gammaL,gammaH,sharpC,cutoff);
mulSpectrums(complex,filter,complex,0);
deslocaDFT(complex);
idft(complex, complex);
//normalize(complex, complex, 0, 1, CV_MINMAX);
planos.clear();
split(complex, planos);
exp(planos[0],planos[0]);
//planos[0] -= Scalar::all(1);
normalize(planos[0], planos[0], 0, 1, CV_MINMAX);
imFiltered = planos[0].clone();
}
void on_trackbar(int, void*){
gammaL = (double) gammaL_slider/10;
gammaH = (double) gammaH_slider/10;
sharpC = (double) sharpC_slider;
cutoff = (double) cutoff_slider;
applyFilter();
imshow("Homomorphic",imFiltered);
}
int main(int argc, char* argv[]){
im = imread(argv[1],CV_LOAD_IMAGE_GRAYSCALE);
namedWindow("Homomorphic", WINDOW_NORMAL);
namedWindow("original",WINDOW_NORMAL);
namedWindow("Filter",WINDOW_NORMAL);
imshow("original",im);
if(!im.data){
cout<<"Nao abriu a imagem!\n";
return -1;
}
dft_M = getOptimalDFTSize(im.rows);
dft_N = getOptimalDFTSize(im.cols);
copyMakeBorder(im, padded, 0, dft_M - im.rows, 0, dft_N - im.cols, BORDER_CONSTANT, Scalar::all(0));
imFiltered = padded.clone();
cout<<"original: "<<im.rows<<'x'<<im.cols<<endl;
cout<<"padded: "<<padded.rows<<'x'<<padded.cols<<endl;
char TrackbarName[50];
sprintf( TrackbarName, "Gamma L x %d", gammaL_max );
createTrackbar( TrackbarName, "Homomorphic", &gammaL_slider, gammaL_max, on_trackbar);
sprintf( TrackbarName, "Gamma H x %d", gammaH_max );
createTrackbar( TrackbarName, "Homomorphic", &gammaH_slider, gammaH_max, on_trackbar);
sprintf( TrackbarName, "C x %d", sharpC_max );
createTrackbar( TrackbarName, "Homomorphic", &sharpC_slider, sharpC_max, on_trackbar);
sprintf( TrackbarName, "Cutoff Frequency x %d", cutoff_max );
createTrackbar( TrackbarName, "Homomorphic", &cutoff_slider, cutoff_max, on_trackbar);
//on_trackbar(0,0);
waitKey(0);
return 0;
}5.1. Trocando Regiões da Imagem
A primeira função implementada é a de trocar regiões da imagem. Como mencionado anteriormente, para que os nossos filtros funcionem corretamente, devemos deslocar a transformada de Fourier da função. Isso é feito de modo semelhante ao exercício de trocar regiões (esse!), com a função deslocaDFT().
void deslocaDFT(Mat& image ){
Mat tmp, A, B, C, D;
image = image(Rect(0, 0, image.cols & -2, image.rows & -2));
int cx = image.cols/2;
int cy = image.rows/2;
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.copyTo(tmp); D.copyTo(A); tmp.copyTo(D);
C.copyTo(tmp); B.copyTo(C); tmp.copyTo(B);
}
A linha image = image(Rect(0, 0, image.cols & -2, image.rows & -2)); implementa um operador lógico bit a bit (&), e faz com que se a imagem tiver um número ímpar de linhas ou colunas, a última linha/coluna é excluída. Com isso temos regiões de tamanhos idênticos.
|
5.2. Declarando as Variáveis Globais
Declaramos algumas variáveis globais, que devem ser modificadas pela função das trackbar e pelo main(). São elas:
-
int gammaL_slider = 2, gammaH_slider = 20, sharpC_slider = 1, cutoff_slider = 5;: Variáveis das Trackbars; -
const int gammaL_max = 10, gammaH_max = 50, sharpC_max = 100, cutoff_max = 200;: Valores máximos das Trackbars; -
int gammaL, gammaH, sharpC, cutoff;: Parâmetros do filtro (\$\gamma_L,\gamma_H,c,D_0\$); -
Mat im, imFiltered,padded;: Matrizes (imagem original, filtrada e com padding); -
int dft_M, dft_N;: Tamanhos ideais para a realização da Transformada;
5.3. Implementação do Filtro
Criamos uma imagem de zeros com o tamanho ideal para a transformada. A equação do filtro homomófirco é implementada através dos dois loops encadeados em cada um dos pixels de uma imagem temporária. Criamos um vetor de duas posições com essa imagem temporária repetida, através da função merge(). O vetor deve ser igual nas duas posĩções pois o filtro é aplicado nas duas camadas da imagem complexa, através da função mulSpectrums(), que será vista posteriormente. A função retorna uma matriz que é o filtro.
Mat homomorphicFilter(double gl, double gh, double c, double d0){
Mat filter = Mat(padded.size(), CV_32FC2, Scalar(0));
Mat tmp = Mat(dft_M, dft_N, CV_32F);
for(int i=0; i<dft_M; i++){
for(int j=0; j<dft_N; j++){
tmp.at<float> (i,j) = (gh - gl)*(1 - exp(-c*(( (i-dft_M/2)*(i-dft_M/2) + (j-dft_N/2)*(j-dft_N/2) ) / (d0*d0) ))) + gl;
}
}
Mat comps[]= {tmp,tmp};
imshow("Filter", tmp);
merge(comps, 2, filter);
// normalize(filter,filter,0,1,CV_MINMAX);
return filter;
}5.4. main()
No main() do programa, recebemos como argumento a imagem, abrimos e realizamos o padding. Para sabermos as dimensões ideais para a transformada, utilizamos a função getOptimalDFTSize(). Já para realizar o padding, utilizamos a função copyMakeBorder(), cujos parâmetros são explicados na documentação do OpenCV.
Em seguida, criamos as Trackbars de modo semelhante ao exemplo anterior do Tilt-Shift.
int main(int argc, char* argv[]){
im = imread(argv[1],CV_LOAD_IMAGE_GRAYSCALE);
namedWindow("Homomorphic", WINDOW_NORMAL);
namedWindow("original",WINDOW_NORMAL);
namedWindow("Filter",WINDOW_NORMAL);
imshow("original",im);
if(!im.data){
cout<<"Nao abriu a imagem!\n";
return -1;
}
dft_M = getOptimalDFTSize(im.rows);
dft_N = getOptimalDFTSize(im.cols);
copyMakeBorder(im, padded, 0, dft_M - im.rows, 0, dft_N - im.cols, BORDER_CONSTANT, Scalar::all(0));
imFiltered = padded.clone();
cout<<"original: "<<im.rows<<'x'<<im.cols<<endl;
cout<<"padded: "<<padded.rows<<'x'<<padded.cols<<endl;
char TrackbarName[50];
sprintf( TrackbarName, "Gamma L x %d", gammaL_max );
createTrackbar( TrackbarName, "Homomorphic", &gammaL_slider, gammaL_max, on_trackbar);
sprintf( TrackbarName, "Gamma H x %d", gammaH_max );
createTrackbar( TrackbarName, "Homomorphic", &gammaH_slider, gammaH_max, on_trackbar);
sprintf( TrackbarName, "C x %d", sharpC_max );
createTrackbar( TrackbarName, "Homomorphic", &sharpC_slider, sharpC_max, on_trackbar);
sprintf( TrackbarName, "Cutoff Frequency x %d", cutoff_max );
createTrackbar( TrackbarName, "Homomorphic", &cutoff_slider, cutoff_max, on_trackbar);
on_trackbar(0,0);
waitKey(0);
return 0;
}5.5. Utilizando as Trackbars
O método que é chamado quando modificamos os valores das Trackbars é o void on_trackbar(int, void*). Esse método converte os valores das Trackbars para pontos flutuantes e chama a função applyFilter(), que aplica o filtro homomórfico. Por fim, a imagem é mostrada.
void on_trackbar(int, void*){
gammaL = (double) gammaL_slider/10;
gammaH = (double) gammaH_slider/10;
sharpC = (double) sharpC_slider;
cutoff = (double) cutoff_slider;
applyFilter();
imshow("Homomorphic",imFiltered);
}5.6. Aplicando o Filtro
O filtro é aplicado através da função applyFilter(). Primeiramente, criamos a parte real e imaginária da imagem. A parte real será a própria imagem (já com o padding) e a parte imaginária será uma matriz de zeros. Aplicamos o \$ln()\$ na parte real, como descrito anteriormente, e realizamos a transformada, com a função dft().
| Antes de aplicarmos o ln, devemos somar 1 à todos os elementos da função para evitar \$ln(0)\$ |
Uma vez transformada, deslocamos as regiões para ajustar a transformada. Calculamos o filtro e utilizamos a função mulSpectrums() para aplicarmos o filtro no domínio da frequência. Deslocamos novamente as regiões e convertemos novamente para o domínio espacial, utilizando a transformada inversa (idft()).
A função resize() é utilizada pois a função deslocaDFT() pode remover uma linha ou coluna, e o filtro e a imagem devem ter as mesmas dimensões para que a função mulSpectrums()`funcione corretamente.
|
Por fim, fazemos a exponencial da imagem, e a normalizamos, para a correta exibição na tela. Devemos pegar somente a parte real da imagem, e por isso utilizamos a função split().
void applyFilter(void){
vector<Mat> planos; planos.clear();
Mat zeros = Mat_<float>::zeros(padded.size());
Mat realInput = Mat_<float>(padded);
Mat complex;
realInput += Scalar::all(1);
log(realInput,realInput);
//normalize(realInput, realInput, 0, 1, CV_MINMAX);
//imshow("logimage",realInput);
planos.push_back(realInput);
planos.push_back(zeros);
merge(planos, complex);
dft(complex, complex);
deslocaDFT(complex);
resize(complex,complex,padded.size());
normalize(complex,complex,0,1,CV_MINMAX);
Mat filter = homomorphicFilter(gammaL,gammaH,sharpC,cutoff);
mulSpectrums(complex,filter,complex,0);
deslocaDFT(complex);
idft(complex, complex);
//normalize(complex, complex, 0, 1, CV_MINMAX);
planos.clear();
split(complex, planos);
exp(planos[0],planos[0]);
planos[0] -= Scalar::all(1);
normalize(planos[0], planos[0], 0, 1, CV_MINMAX);
imFiltered = planos[0].clone();
}6. A Interface do Programa
Assim que o programa é aberto, são mostradas 3 janelas:
-
A primeira mostra a imagem filtrada, juntamente com as 4 trackbars de controle;
-
A segunda mostra a imagem original
-
A terceira mostra o filtro
7. Corrigindo a Luminosidade de Imagens
A seguir teremos dois exemplos da aplicação do filtro. O primeiro é um exemplo cĺássico da imagem do túnel, que é mostrado no site do Matlab. A correção de luminosidade é mostrada logo em seguida.
A segunda cena foi preparada pelo autor do texto utilizando um quarto escuro e um led. A foto original é mostrada, assim como o resultado da filtragem logo em seguida.
