Italiano - English

Implementazione di Gaussian Smoothing ricorsivo

Il filtro gaussiano è molto utilizzato per il preprocessing di segnali. In particolare è un filtro passa basso utile a rimuovere il rumore di alta frequenza (Smoothing). In questa pagina viene illustrata una implementazione ricorsiva 2D secondo il metodo di T.Youg, Lucas J, van Vliet

Filtro Gaussiano 2D Ricorsivo

L'algoritmo si basa sulla trasformata Z di una approssimazione della trasformata di Fourier della funzione gaussiana continua. La trasformata produce due equazioni alle differenze che applicate al segnale in ingresso calcolano il segnale filtrato.

Gli esperimenti eseguiti mostrano che la tecnica pruduce risultati molto vicini all'algoritmo di convoluzione standard tramite FFT. Dal punto di vista delle prestazioni l'algoritmo risulta piu' lento della convoluzione tramite FFT tuttavia la sua implementazione è decisamente semplice, al contrario di una FFT ottimizzata.

Di seguito viene illustrata una implementazione con Matlab(r). I risultati vengono poi confrontati con una convoluzione con FTT

Copyright (c) 2010 PkLab.net

Indice

Un po di teoria sul filtro gaussiano ricorsivo

L'equazione per la gaussiana vale:

La gaussiana nelle 2 dimensioni vale

Dato che la gaussiana è una funzione separabile la sua versione 1d puo' essere applicata ripetutamente a tutte le dimensioni del segnale da filtrare.

Data una immagine I, per ottenere l'immagine filtrata Ig applicare il filtro Gaussiano 1d a cascata sulle due dimensioni dell'immagine

L'algoritmo ricorsivo approssima la funzione esponenziale:

con la seguente funzione polinomiale:

Questa approssimazione introduce un errore epsilon facilmente calcolabile:

clear;
T = -5:0.01:5;
% Pure exponential Gaussian
Gt = exp(-T.^2/2) / sqrt(2*pi);

% Polinomial Gaussian
a0 = 2.490895;  a2 = 1.466003;
a4 = -0.024393; a6 = 0.178257;
At = 1./(a0 + a2*T.^2 + a4*T.^4 + a6*T.^6);

% Approsimation error rate relative to maximun G(t=0)
Et = 100*(Gt-At)*sqrt(2*pi);

% Locate max and min error
maxe=find(Et==max(Et));
mine=find(Et==min(Et));

iptsetpref('ImshowBorder','loose');
iptsetpref('ImshowInitialMagnification','fit');
iptsetpref('ImshowAxesVisible','on');
win = 340;
h = 1;

%plot results
h=figure(h);
plot(T,Gt,'b');
hold on;
plot(T,At,'m');
plot(T,Et,'r')

stem(T(maxe),Et(maxe),'--');text(T(maxe)+.2,Et(maxe),'Max \epsilon [%]');
stem(T(mine),Et(mine),'--');text(T(mine)+.2,Et(mine),'Min \epsilon [%]');
legend('g(t) \it Pure Gaussian ','a(t) \it Approsimate Gaussian','\epsilon(t) \it Relative Error %');
grid;
hold off;

Come si nota la funzione polinomiale approssimata è molto vicina alla funzione esponenziale. L'errore di approssimazione e' circa 0.6% rapportato al valore massimo per l'esponenziale. Relativamente all'andamento dell'errore, i valori negativi corrispondono ad una leggera amplificazione della funzione approssimata rispetto alla funzione esponenziale. Allo stesso modo i valori positivi indicano una leggera riduzione.

Particolarmente interessante e' l'abbinamento amplificazione sulla cuspide seguita e preceduta da una da riduzione. Ciò produce un aumento del contrasto in corrispondenza degli edge in una immagine

Algoritmo per il filtro gaussiano ricorsivo

Dalla approssimazione polinomiale sopra descritta, con opportuni passaggi e trasformate si ottengono le due equazioni alle differenze:

Equazione alle differenze in avanti

Equazione alle differenze all'indietro

E' introdotto il parametro q funzione del parametro sigma della gaussiana. I parametri B,b0,b1,b2,b3 sono funzione di q. L'algoritmo è molto semplice

Nota:L'algoritmo è definito solo per sigma>=0.5

Inquanto approssimazione, l'algoritmo è soggetto ad un errore che decresce con l'aumentare del parametro sigma scelto. Per sigma=5 vale circa 0.5E-02. A tal proposito è bene ricordare che tutti gli algoritmi per il filtro gaussiano sono soggetti ad errore di approssimazione, normalmente dipendente dal parametro sigma

Gaussian Filters Accuracy

Per maggiori dettagli consultare [1]

Costruzione immagine in ingresso

sigma = 2.6;
step=49;
thin=1;
num = 5;
step=step+thin;
dim = (num*step)+thin;
Ig = ones(dim,dim);
for g=1:step:dim
  Ig(g:g+thin,:) = 0;
  Ig(:,g:g+thin) = 0;
end
I = im2double(Ig); % conversione immagine in double
h=figure(h+1);set(h,'Position',[100 100 win win]);
imshow(Ig),title('Input Image');

La conversione in double viene eseguita per applicare il metodo FFT

Calcolo del parametro q

Il parametro q è funzione di sigma ed è definito come segue:

S=0.5:0.01:8;
Q = S;
for n=1:length(S)
  if (S(n) >= 2.5)
    Q(n) = 0.98711 * S(n) - 0.96330;
  elseif (S(n) >= 0.5)
    Q(n) = 3.97156 - 4.14554 * sqrt(1 - 0.26891 * S(n));
  end;
end
h=figure(h+1);set(h,'Position',[100 100 win win]);
plot(S,Q);grid;xlabel('sigma');ylabel('q');title('q = f (sigma)');


parametro q in funzione di sigma

disp(sprintf('Sigma = %3.2f',sigma));

if (sigma >= 2.5)
  q = 0.98711 * sigma - 0.96330;
elseif (sigma >= 0.5)
  q = 3.97156 - 4.14554 * sqrt(1 - 0.26891 * sigma);
else
  error('Gaussian sigma must be greater or equal than 0.5');
end;
Sigma = 2.60

Calcolare i coefficienti del filtro gaussiano

Il parametro q ed i coefficienti sono stati derminati con la trasformata Z

b0 = 1.57825 + (2.44413 * q) + (1.4281 * q^2) + (0.422205 * q^3);
b1 = (2.44413 * q)  + (2.85619 * q^2) + (1.26661 * q^3);
b2 = -(1.4281 * q^2) - (1.26661 * q^3);
b3 = 0.422205 * q^3;
B  = 1 - ( ( b1 + b2 + b3 )/ b0 );

Per migliorare l'efficienza del calcolo, i coefficienti vengono scalati per il fattore b0

B1 = b1/b0;
B2 = b2/b0;
B3 = b3/b0;

Risoluzione dei bordi

Ogni convoluzione pone il problema della gestione dei bordi. Per ogni pixel, il calcolo del filtro ricorsvivo richiede la presenza di almeno 3 pixel precedenti e 3 pixel successivi. Questo implica che le prime 3 colonne (o righe) non possono essere calcolate. Per risolvere questo problema l'immagine viene allargata di una quantità di pixel sufficiente per applicare l'algoritmo, ovvero 3 pixel.

Se nell'intorno dei bordi, immagine non è interessante ai fini delle elaborazioni successive(come spesso accade), l'allargamento si può evitare, considerato che per la gaussiana ricorsiva, riguarda 3 pixel, indipendentemente dal valore di sigma

Nel caso di convoluzione con kernel gaussiano standard, i pixel interessati sono normalmente più numerosi e comunque in numero dipendente da sigma. Ad esempio per un sigma=2.6 il kernel per la convoluzione gaussiana ha dimensioni 9x9 In questi casi è consigliabile allargare l'immagine della quantità pari alla metà del kernel.

tic   % start performace counter

useborder = true; %
if (useborder)
  border=3;
  G = padarray(I,[border border],'replicate','both');
else
  G = I;
end

[R C] = size(G);
W = G; %preallocate matrix

Applico equazione alle differenze in avanti sulle colonne

Una immagine è un segnale a 2 dimensioni(righe e colonne). Utilizzando la proprietà della gaussiana di essere separabile, l'algoritmo viene applicato a cascata prima lungo una direzione e sul risultato viene riapplicato lungo l'altra direzione.

Ogni riga dell'immagine viene considerata come un singolo vettore i cui elmenti sono i pixel delle colonne per la data riga. Dunque la riga è costante mentre la colonna e variabile.

ovvero

NOTA: nel codice seguente l'immagine la variabile G contiene l'immagine I precedentemente allargata per risolvere i bordi.

for r = 1:R
  W(r,1) = B * G(r,1);
  W(r,2) = B * G(r,2) + B1 * W(r,1);
  W(r,3) = B * G(r,3) + B1 * W(r,2) + B2 * W(r,1);
  for c = 4:C
    W(r,c) = B * G(r,c) + ...
      B1 * W(r,c-1) + B2 * W(r,c-2) + B3 * W(r,c-3);
  end
end

Applico equazione alle differenze all'indietro sulle colonne

Equazione all'indietro viene applicata dall'ultimo elemento verso il primo

ovvero

for r = 1:R
  G(r,C)   = B * W(r,C);
  G(r,C-1) = B * W(r,C-1) + B1 * G(r,C);
  G(r,C-2) = B * W(r,C-2) + B1 * G(r,C-1) + B2 * G(r,C);

  % dall'ultima colonna verso la prima
  for c = C-3:-1:1
    G(r,c) = B * W(r,c) + ...
      B1 * G(r,c+1) + B2 * G(r,c+2) + B3 * G(r,c+3);
  end;

end;

A questo punto G contiene l'immagine filtrata lungo la direzione delle colonne

rtime = toc;
h=figure(h+1);set(h,'Position',[h*100 100 win win]);
imshow(G);title('Gaussian applicato alle colonne');
tic

Applico equazione alle differenze in avanti sulle righe

Se al risultato G applichiamo i due passi precedenti utilizzando la riga come variabile e la colonna come costante otterremo l'immagine filtrata nelle due direzioni. L'equazione in avanti sulle righe diventa la seguente

NOTA: nel codice seguente la variabile G contiene l'immagine allargata e filtrata lungo le colonne

for c=1:C
  W(1,c) = B * G(1,c);
  W(2,c) = B * G(2,c) + B1 * W(1,c);
  W(3,c) = B * G(3,c) + B1 * W(2,c) + B2 * W(1,c);
  for r=4:R
    W(r,c) = B * G(r,c) + ...
      B1 * W(r-1,c) + B2 * W(r-2,c) + B3 * W(r-3,c);
  end
end

Applico equazione alle differenze all'indietro sulle righe

Allo stesso modo del passo precedente l'equazione all'indietro diventa

for c = 1:C
  % Backward equation is applied from last row down to first row
  G(R,c)   = B * W(R,c);
  G(R-1,c) = B * W(R-1,c) + B1 * G(R,c);
  G(R-2,c) = B * W(R-2,c) + B1 * G(R-1,c) + B2 * G(R,c);
  for r = R-3:-1:1
    G(r,c) = B * W(r,c) + ...
      B1 * G(r+1,c) + B2 * G(r+2,c) + B3 * G(r+3,c);
  end;
end;

Ripristino dimensione immagine

Ora la variabile G contiene immagine filtrata lungo entrambe le direzioni, a meno degli eventuali bordi aggiunti per applicare l'agoritmo

I bordi aggiunti in precedenza vengono rimossi cosi' G contiene il risultato

if (useborder)
  G = G (border + 1:R-border, border + 1:C-border);
end

rtime = rtime+toc;  % stop performace counter

Gaussian Smooting with Convolution (Use FFT)

Qui viene riportato a scopo dimostrativo e comparativo lo smoothing gaussiano operato con la procedura standard di convoluzione

Nota:Matlab utilizzaza FFT per la convoluzione

% Calcolo della dimensione del kernel per la convoluzione
cutAt = .0001;
t = 1:30;
ssq = sigma^2;
width = find( exp(-(t.^2)/(2*ssq))/(sigma*sqrt(2*pi)) >cutAt,1,'last');
if isempty(width)
  width = 1;
end

Compilo il kernel con i valori della gaussiana 1d

t = (-width:width);
gau = exp(-(t.*t)/(2*ssq))/(sigma*sqrt(2*pi));

Convoluzione con la gaussiana 1d , prima sulle righe con gau poi sulle colonne utilizzando il kernel trasposto gau'

tic
Cx = imfilter(I,gau, 'conv', 'replicate');
Cxy = imfilter(Cx, gau' , 'conv', 'replicate');
ftime=toc;

Ripristino scala dell'immagine

La gaussiana introduce un fattore di scala pari al suo valore massimo, che dato da x=0 nel caso 1d e x=0 & y=0 nel caso 2d

Nel caso dell'algoritmo ricorsivo, il parametro B porta in scala il segnale filtrato. La convoluzione come sopra operata non richiede la scalatura dell'output

Risultati

Di seguito vengono riportati i risultati comparati con l'algoritmo di riferimento che fa uso della convoluzione (FFT Gaussian).

% seleziono un profilo
profile = round(2.5*step);

h=figure(h+1);set(h,'Position',[h*100 100 win win]);
imshow(Cxy);title('FFT Gaussian');
% mostro il profilo selezionato sull'immagine
hold on;plot([0 C],[profile profile],'-.');hold off

h=figure(h+1);set(h,'Position',[h*100 100 win win]);
imshow(G);title('Recursive Gaussian Image');
% mostro il profilo selezionato sull'immagine
hold on;plot([0 C],[profile profile],'m-.');hold off

% Immagine differenza tra le due immagini filtrate
DeltaG = (Cxy - G);
h=figure(h+1);set(h,'Position',[h*100 100 win win]);
imshow(DeltaG,[]);title('Immagine Delta (FFT - Recursive)');

% Istogramma differenza
h=figure(h+1);set(h,'Position',[h*100 100 win win]);
hist(DeltaG(:),1000);
title('Difference distribution (FFT - Recursive)');grid;


Come era atteso esiste una differenza di risultato tra i due algoritmi. Se si considera il risultato ottenuto tramite la FFT come gold standard, allora la differenza tra le due immagini filtrate rappresenta l immagine delta dell'algoritmo ricorsivo (l'immagine è opportunamente scalata nei valori)

Osservando l immagine delta si nota il reticolo chiaro. Questi sono gli effetti della approssimazione polinomiale sopra descritti, infatti seguono l'andamento dell'errore epsilon.

L'istogramma dell immagine delta fornisce una indicazione visuale della distribuzione di frequenza della differenza stessa. Si nota che la maggior parte della differenza si concentra intorno al valore medio vicino allo 0

disp(sprintf('Mean error = %d',mean2(DeltaG)));
disp(sprintf('Var error  = %d',var(DeltaG(:))));
Mean error = 5.923444e-004
Var error  = 3.221474e-005

Per osservare a fondo le implicazione di questa differenza è possibile osservare le derivate di un profilo (esempio lungo la linea tratteggiata) che attraversa le due immagini filtrate

%derivata lungo il profilo selezionato
dg = diff(G(profile,:));
dc = diff(Cxy(profile,:));

%plot delle derivate
h = figure(h+1);set(h,'Position',[h*100 100 win win]);
plot(dc);
hold on;
plot(dg,'m');
seledge = 2*step;
select =seledge-20:seledge+20;
mind = 1.1*min(dg(select));
maxd = 1.1*max(dg(select));

% seleziono area per lo zoom
rectangle('Position',[2*step-20,mind,40,maxd-mind],...
          'Curvature',[1,1],...
         'LineWidth',1,'LineStyle','--')

hold off;
grid;title('Derivata lungo profilo selezionato');
legend('FFT Gaussian','Recursive Gaussian');
xlabel(sprintf('Riga %3.0d',profile));

% Mostro zoom della derivata nel zona selezionata
h=figure(h+1);
set(h,'Position',[h*100 100 win win]);

plot(select,dc(select));
hold on;

% mostro edge dell'immagine originale
rectangle('Position',[seledge,mind,thin,maxd-mind],'FaceColor','black');

plot(select,dc(select));      %derivata della gaussiana FFT
plot(select,dg(select),'m');  %derivata della gaussiana ricorsiva

hold off;
grid;title('Zoom derivata su edge');
%legend('FFT Gaussian','Recursive Gaussian');
xlabel(sprintf('Riga %3.0d',profile));

Come si nota, in corrispondenza degli edge, la derivata del per l'algoritmo ricorsivo è leggermente maggiore rispetto al metodo FFT. Anche questo risultato era atteso dopo l'analisi dell'errore epsilon

Velocità

Riguardo alle prestazioni dell'algoritmo in termini di velocità si misura una decisa inversione rispetto a quanto dichiarato in [1].

FFT exec time = 25.094 ms 
       Recursive exec time = 42.957 ms 
 Recursive/FFT Speed ratio = 1.71 
 

I tempi di esecuzione su Pentium M 1.86Gz WinXp Sp.3

In particolare l'algoritmo ricorsivo risulta essere circa 2 volte piu' lento di Gaussian FFT, mentre in [1] è riportato essere circa 3 volte piu' veloce.

A tal proposito va considerato che

  1. L'implementazione dell'algoritmo ricorsivo qui realizzata può essere ottimizzata
  2. L'implementazione dell'algoritmo FFT implementata da Matlab è altamente ottimizzata
  3. L'articolo [1] risale al 1995 e si riferisce ad hardware e tecniche ormai obsolete.
% disp(sprintf('            FFT exec time = %3.3f ms',ftime*1000));
% disp(sprintf('      Recursive exec time = %3.3f ms',rtime*1000));
% disp(sprintf('Recursive/FFT Speed ratio = %3.2f',rtime/ftime));

Esempio con immagini reali

CallGaussian;
Mean error = 3.272933e-003
 Var error = 2.617660e-004

Riferimenti

  1. Young I.T. and L.J. Van Vliet, Recursive Implementation of the Gaussian Filter. Signal Processing, 1995. 44(2): p. 139-151 http://www.ph.tn.tudelft.nl/~lucas/publications/1995/SP95TYLV/SP95TYLV.html
  2. Young I.T., Gerbrands J.J and L.J. Van Vliet, Image Processing Fundamentals 2.3: p 57
  3. LIBROW, Gaussian filter, or Gaussian blur http://www.librow.com/articles/article-9
  4. Image Processing Learning Resources, Gaussian Smoothing http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm
  5. WikiPedia, Gaussian blur http://en.wikipedia.org/wiki/Gaussian_blur
Vota questa pagina:

0 Commenti:

Lascia il tuo commento:

Note:
  • La tua email non è obligatoria e non sarà visibile in alcun modo
  • Si prega di inviare solo commenti relativi a questa pagina
  • Commenti inappropriati o offensivi saranno modificati o eliminati
  • Codici HTML non sono consentiti. Prego usare i BB code:
    [b]bold[/b], [u]underline[/u], [i]italic[/i], [code]code[/code]
Il codice, le illustrazioni e gli esempi riportati in questa pagina sono solo a scopo illustrativo. L'autore non prende alcuna responsabilità per il loro utilizzo da parte dell'utente finale.
Questo materiale è di proprietà di Pk Lab ed è utilizzabile liberamente a condizione di citarne la fonte.