La convoluzione

La convoluzione rappresenta l'operatore fondamentale utilizzato per
filtrare linearmente dei segnali al fine di evidenziare delle frequenze o
attenuarne altre.

Contents

Introduzione alla convoluzione

Come noto un filtro è un sistema che trasforma un segnale in un altro, esso stesso è definito tramite il segnale di risposta all'impulso che si vuole ottenere. In ambito digitale il segnale impulso non è la delta di Dirac ma più semplicemente un singolo valore 1 in un segnale contenente soli valori nulli. Si verifichi l'implementazione della convoluzione qui.

% Definiamo l'impulso discreto:
I = zeros(10,10);
I(1,1) = 1;

% Creiamo una risposta all'impulso desiderata:
g = fspecial('gaussian',9,1);
g = g * 10/ sum(g(:)); % Solo per visualizzazione.

% Verifichiamo che l'impulso si comporti come identità:
res = convolution(g,I);
figure;
subplot(131); surf(I);   axis equal; title('Impulso');
subplot(132); surf(g);   axis equal; title('Filtro Gaussiano');
subplot(133); surf(res); axis equal; title('Risultato della convoluzione');

Il guadagno di un filtro è dato dalla somma dei suoi coefficienti:

img = double(imread('cameraman.tif'))/256;
g = fspecial('gaussian',9,1);
g1 = g / sum(g(:));
g2 = g1 * 2;
figure;
subplot(221); imshow(img); title('Immagine originale');
subplot(222); surf(g1); title('Filtro gaussiano');
subplot(223); imshow(convolution(g1,img)); title('Guadagno 1');
subplot(224); imshow(convolution(g2,img)); title('Guadagno 2');

Filtri reali e ideali

% Inizializzazione:
img = double(imread('cameraman.tif'))/256;
g = fspecial('gaussian',256,1);

% Lavoriamo nella frequenza:
IMG = fftshift(fft2(img));
G = fftshift(fft2(g));
IMGO = IMG.*G;
imgo = abs(fftshift(ifft2(IMGO)));

% Vediamo il risultato:
figure;
subplot(221); surf(g); title('Filtro Gaussiano');
subplot(222); surf(abs(G)); title('Filtro Gaussiano in frequenza');
subplot(223); imshow(img); title('Immagine');
subplot(224); imshow(imgo); title('Immagine filtrata');

Siamo riusciti a ritenere le basse frequenze con un'approssimazione di un
filtro reale, se provassimo a forzare il filtro reale?
% Lavoriamo nella frequenza:
Gmag = abs(G) > 0.5;
Gphase = angle(G);
G = Gmag.*exp(1i*Gphase);
IMGO = IMG.*G;
imgo = abs(fftshift(ifft2(IMGO)));

% Vediamo il risultato:
figure;
subplot(121); imshow(abs(G)); title('Filtro ideale in frequenza');
subplot(122); imshow(imgo); title('Immagine filtrata');

Filtri di convoluzione separabili

Un filtro lineare generico di NxN coefficienti può essere convoluto con
un'immagine MxM tramite l'algoritmo vanilla con complessità:

$$O(N^2(N+M)^2)$$

Alcuni filtri però sono separabili, ovvero possono essere rappresentati come convoluzione di filtri di dimensionalità inferiore. Ad esempio:

$$e^{-\frac{1}{2}\frac{x^2+y^2}{\sigma^2}} = e^{-\frac{1}{2}\frac{x^2}{\sigma^2}}e^{-\frac{1}{2}\frac{y^2}{\sigma^2}}$$

quindi nel dominio delle frequenze moltiplicare per un filtro gaussiano bivariato può essere ridotto a due moltiplicazioni per due filtri gaussiani monovariati, quindi nel dominio spaziale:

$$img \star g_{x,y} = img \star g_x \star g_y$$

la cui complessità computazionale è:

$$O(2N(N+M)^2)$$

% Scomponiamo il filtro di prima:
gx = g1(ceil(size(g1,1)/2),:); gx = gx / sum(gx);
gy = g1(:,ceil(size(g1,1)/2)); gy = gy / sum(gy);

% Applichiamolo in entrambe i modi:
imgo1 = convolution(g1,img);
imgo2 = convolution(gy,convolution(gx,img));

% Confrontiamo i risultati:
fprintf('Errore: %g\n',max(abs(imgo1(:)-imgo2(:))));

Errore: 1.33227e-15

Filtri Gaussiani e DoG

Un filtro gaussiano rappresenta un'approssimazione di un filtro
passa-basso ideale, nel dominio delle frequenze permette di ritenere
frequenze al di sotto di una determinata frequenza di soglia legata alla
deviazione standard come segue:

$$\phi = \frac{1}{2\pi\sigma}$$

Se si è interessati a una determinata banda di frequenza è possibile
comporre i risultati dell'applicazione di due differenti filtri, ad
esempio se si intendono mantenere le frequenze tra 0.08 (sigma = 2) e
0.16 (sigma = 1) è possibile applicare due filtri e calcolarne la
differenza:

$$img \star g_{\sigma=1} - img \star g_{\sigma=2} = img \star (g_{\sigma=1} - g_{\sigma=2})$$

% Preparazione dei filtri:
g1 = fspecial('gaussian',13,1);
g2 = fspecial('gaussian',13,2);

% Applicazione:
imgo = filter2(g1-g2,img);

% Vediamo il risultato:
figure;
subplot(121); imshow(img); title('Immagine originale');
subplot(122); imshow(imgo*2); title('Filtro passa-banda'); % *2 per visualizzazione.

Il risultato mostra solamente le componenti in una determinata banda di
frquenza, ovvero i bordi dell'immagine appartenenti a quella banda.
Un'immagine può essere scomposta in componenti a frequenze differenti
evidenziando valor medio e bordi a differenti frequenze. DoG sta per
Difference of Gaussians e permette di ricavare per ogni pixel dei
coefficienti che ne descrivano in frequenza il comportamento.
Avendo scelto in questo esempio 2 frequenze abbiamo 3 componenti:
% Componente a frequenza più bassa:
imglow = filter2(g2,img);

% La banda di frequenza selezionata prima:
imgmid = imgo;

% Il resto a frequenza più alta:
imghigh = img - imglow - imgmid;

% Vediamoli:
figure;
subplot(221); imshow(imglow);
subplot(222); imshow(imgmid);
subplot(223); imshow(imghigh);
subplot(224); imshow(imglow+imgmid+imgmid);


Published with MATLAB® R2012b