>
%
% 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:
%
% $$\mathbf{if}(\sigma \geq 2.5) \;\mathbf{then} \;q = 0.98711 \sigma - 0.96330$$
%
% $$\mathbf{if}(\sigma \geq 0.5) \;\mathbf{then} \;q = 3.97156 - 4.14554 \sqrt{1 - 0.26891 \sigma}$$
%
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)');
%%
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;
%% 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.
%
% $$W[r,c] = B\,I[r,c] + (b_{1}w[r,c-1] + b_{2}w[r,c-2] +
% b_{3}w[r,c-3])/b0$$
%
% ovvero
%
% $$W[r,c] = B\,I[r,c] + B_{1}w[r,c-1] + B_{2}w[r,c-2] + B_{3}w[r,c-3]$$
%
% *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
%
% $$G[r,c] = B\,W[r,c] + (b_{1}G[r,c+1] + b_{2}G[r,c+2] + b_{3}
% G[r,c+3])/b0$$
%
% ovvero
%
% $$G[r,c] = B\,W[r,c] + B_{1}G[r,c+1] + B_{2}G[r,c+2] + B_{3} G[r,c+3]$$
%
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
%
% $$W[r,c] = B\,G[r,c] + B_{1}W[r-1,c] + B_{2}W[r-2,c] + B_{3}W[r-3,c]$$
%
% *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
%
% $$G[r,c] = B\,W[r,c] + B_{1}G[r+1,c] + B_{2}G[r+1,c] + B_{3} G[r+3,c]$$
%
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
%
% $$G(x)=\frac{1}{\sigma \, \sqrt[]{2 \pi }}e^{-\frac{x^2}{2\sigma^2}}
% \Rightarrow G(x=0)=\frac{1}{\sigma \, \sqrt[]{2 \pi }}$$
%
% $$G(x,y)=\frac{1}{2\pi \sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}
% \Rightarrow G(x=0,y=0)=\frac{1}{2 \pi \sigma^2}$$
%
% 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(:))));
%%
% 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
%
% # L'implementazione dell'algoritmo ricorsivo qui realizzata può essere
% ottimizzata
% # L'implementazione dell'algoritmo FFT implementata da Matlab è altamente
% ottimizzata
% # 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;
%% Riferimenti
% # Young I.T. and L.J. Van Vliet, _Recursive Implementation of the Gaussian
% Filter._ Signal Processing, 1995. 44(2): p. 139-151
%
% # Young I.T., Gerbrands J.J and L.J. Van Vliet, _Image Processing
% Fundamentals 2.3_: p 57
% # LIBROW, _Gaussian filter, or Gaussian blur_
%
% # Image Processing Learning Resources, _Gaussian Smoothing_
%
% # WikiPedia, _Gaussian blur_
%