Wstęp do analizy obrazu: Różnice pomiędzy wersjami

Z Brain-wiki
(Utworzono nową stronę "Laboratorium_EEG/ Wstęp do analizy obrazu")
 
 
(Nie pokazano 13 wersji utworzonych przez 2 użytkowników)
Linia 1: Linia 1:
 
[[Laboratorium_EEG]]/ Wstęp do analizy obrazu
 
[[Laboratorium_EEG]]/ Wstęp do analizy obrazu
 +
 +
=Częstości w obrazie i dwuwymiarowa transformata Fouriera=
 +
O analizie obrazów za pomocą Pythona można poczytać [[Pracownia_EEG/analiza_obrazu|tu]].
 +
 +
Obrazy możemy analizować tak jakby były to sygnały 2D.
 +
Aby uświadomić sobie związki między jednowymiarowym sygnałem a obrazem — sygnałem 2D — zrobimy następujące ćwiczenie:
 +
 +
Generujemy w jednym wymiarze funkcję sinus o 10 okresach na 100 punktów:
 +
<source lang = matlab>
 +
X=100; %pikseli
 +
Fs=1/X; %czestosc probkowania co jeden piksel
 +
dx=1;
 +
x=0:dx:X-dx;
 +
fx=10/X; % czestosc sinusa 10 okresów na X
 +
syg=sin(2*pi*x*fx+pi/3);
 +
figure(1)
 +
plot(x,syg)
 +
</source>
 +
 +
Ten obrazek jest znajomy. Teraz przekształcimy go w obraz:
 +
<source lang = matlab>
 +
figure(2)
 +
% teraz zrobimy z niego obraz 2D
 +
Y=X/2;
 +
SYG=zeros(Y,X);
 +
 +
for i=1:Y
 +
    SYG(i,:)=syg;
 +
end
 +
 +
 +
subplot(2,2,1)
 +
imagesc(SYG)
 +
set(gca,'Ydir','normal');
 +
colormap(gray)
 +
</source>
 +
 +
Te zmieniające intensywność prążki to właśnie obrazek przedstawiający funkcję sinus w kierunku X.
 +
Zobaczmy jak wygląda jego transformata Fouriera:
 +
<source lang = matlab>
 +
subplot(2,2,2)
 +
S=fft2(SYG);
 +
SK_X=-X/2:X/2;
 +
SK_Y=-X/2:X/2;
 +
imagesc(SK_X,SK_Y,abs(fftshift(S)))
 +
set(gca,'Ydir','normal');
 +
colormap(gray)
 +
</source>
 +
 +
Oraz ten sam sygnał z dodaną składową stałą i jego transformata:
 +
<source lang = matlab>
 +
subplot(2,2,3)
 +
 +
imagesc(SYG+1)
 +
set(gca,'Ydir','normal');
 +
colormap(gray)
 +
 +
subplot(2,2,4)
 +
S=fft2(SYG+1);
 +
imagesc(SK_X,SK_Y,abs(fftshift(S)))
 +
set(gca,'Ydir','normal');
 +
colormap(gray)
 +
pause
 +
 +
</source>
 +
 +
=Korelacja i splot dla obrazu=
 +
Przyjmijmy konwencję:
 +
* obraz = duża macierz
 +
* jądro splotu/korelacji = mniejsza macierz
 +
 +
załóżmy, że mamy obraz A:
 +
<source lang = matlab>
 +
A = [17  24  1  8  15
 +
    23  5  7  14  16
 +
      4  6  13  20  22
 +
    10  12  19  21  3
 +
    11  18  25  2  9];
 +
</source>
 +
 +
i jądro korelacji 
 +
 +
<source lang = matlab>
 +
h = [8  1  6
 +
    3  5  7
 +
    4  9  2];
 +
</source>
 +
 +
wówczas operacja korelacji dla elementu A(2,4) jest następująca:
 +
# kładziemy macierz h elementem środkowym (5) na elemencie A(2,4) (czyli 14)
 +
# wymnażamy przez siebie wszystkie pokrywające się elementy i sumujemy iloczyny
 +
 +
 +
<source lang = matlab>
 +
A(2,4)=1*8 + 8*1 +15*6 +
 +
  7*3 + 14*5 + 16*7+
 +
  13*4 + 20*9 +22*2
 +
</source>
 +
 +
W matlabie mamy funkcję <tt>filter2(h,A)</tt>, która wykonuje powyższą operację- dwuwymiarową korelację.
 +
 +
Bardzo podobna jest operacja dwuwymiarowego splotu. Jedyna różnica jest taka, że przed operacją mnożenia i dodawania jądro jest obracane o 180 stopni.
 +
 +
<source lang = matlab>
 +
>> rot90(rot90(h))
 +
 +
ans =
 +
 +
    2    9    4
 +
    7    5    3
 +
    6    1    8
 +
</source>
 +
Operacja splotu jest zaimplementowana jako <tt>conv2</tt>
 +
 +
Zarówno operacja splotu jak i korelacji jest kiepsko określona na brzegach obrazu, tam gdzie jądro wystaje poza brzeg.
 +
W obu przypadkach wykonywane jest domyślnie dopełnianie obrazu zerami.
 +
Efekt ten widać w poniższym przykładzie jako czarne krawędzie.
 +
 +
<source lang = matlab>
 +
I = single(imread('coins.png'));
 +
h = ones(5,5) / 25;
 +
I2 = filter2(h,I);
 +
I3 = conv2(h,I);
 +
subplot(2,2,1)
 +
imagesc(I), title('Pierwotny obraz');
 +
subplot(2,2,2)
 +
imagesc(I2), title('Filtrowany obraz')
 +
subplot(2,2,3)
 +
imagesc(I3), title('Splatany obraz')
 +
colormap('gray')
 +
</source>
 +
 +
=Dwuwymiarowe filtry FIR=
 +
Są naturalnym rozszerzeniem jednowymiarowych filtrów FIR i dają się łatwo zapisać w postaci macierzowej.
 +
 +
==Metoda transformacji częstości==
 +
 +
Projektujemy filtr jednowymiarowy i transformujemy go na 2D:
 +
 +
<source lang = matlab>
 +
b = remez(10,[0 0.4 0.6 1],[1 1 0 0]);
 +
h = ftrans2(b);
 +
[H,w] = freqz(b,1,64,'whole');
 +
colormap(jet(64))
 +
plot(w/pi-1,fftshift(abs(H)))
 +
figure, freqz2(h,[32 32])
 +
</source>
 +
 +
Zprojektujmy i zastosujmy filtr górnoprzepustowy do obrazka:
 +
<source lang = matlab>
 +
b = fir1(14,0.25,'high');
 +
h = ftrans2(b);
 +
[H,w] = freqz(b,1,64,'whole');
 +
colormap(jet(64))
 +
subplot(2,2,1)
 +
plot(w/pi-1,fftshift(abs(H)))
 +
subplot(2,2,2), freqz2(h,[32 32])
 +
 +
load spine
 +
 +
subplot(2,2,3)
 +
imagesc(X)
 +
colormap bone
 +
 +
 +
subplot(2,2,4)
 +
imagesc(filter2(h,X))
 +
 +
</source>
 +
 +
==Kilka konkretnych filtrów==
 +
A teraz kilka powszechnie stosowanych filtrów. Proszę zbadać funkcję odpowiedzi impulsowej i działanie filtra na nasz przykładowy obrazek.
 +
;Filtr uśredniający:
 +
<source lang = matlab>
 +
1/9 1/9 1/9
 +
1/9 1/9 1/9
 +
1/9 1/9 1/9
 +
</source>
 +
;Filtr wyostrzający:
 +
<source lang = matlab>
 +
0 -1 0
 +
-1 5 -1
 +
0 -1 0
 +
</source>
 +
 +
;Filtr uwypuklający:
 +
<source lang = matlab>
 +
-1 -1 -1 -1 0
 +
-1 -1 -1 0 1
 +
-1 -1 0 1 1
 +
-1 0 1 1 1
 +
0 1 1 1 1
 +
</source>
 +
 +
=Filtry nieliniowe=
 +
Dotychczas rozważane filtry obliczały lokalną korelację jądra filtru z kolejnymi fragmentami obrazu, czyli wykonywały tam operację liniową. Łatwo sobie wyobrazić, że w kolejnych &bdquo;okienkach&rdquo; filtra będziemy wykonywać dowolną inną operację np. medianę i jej wartość wpisywać jako element wynikowego obrazka.
 +
 +
Wczytajmy przykładowy obrazek:
 +
<source lang = matlab>
 +
I = imread('eight.tif');
 +
subplot(2,2,1), imshow(I)
 +
title('oryginalny obraz')
 +
</source>
 +
 +
dodajemy szum:
 +
 +
<source lang = matlab>
 +
J = imnoise(I,'salt & pepper',0.02);
 +
subplot(2,2,2), imshow(J)
 +
</source>
 +
 +
Filtrujemy obrazek filtrem uśredniającym
 +
 +
<source lang = matlab>
 +
K = filter2(fspecial('average',3),J)/255;
 +
subplot(2,2,3), imshow(K)
 +
</source>
 +
 +
A teraz medianowym
 +
 +
<source lang = matlab>
 +
L = medfilt2(J,[3 3]);
 +
subplot(2,2,4), imshow(L)
 +
</source>
 +
=Ćwiczenia=
 +
* Proszę pobrać i wczytać do Matlaba [http://www.fuw.edu.pl/~jarekz/SYGNALY/TF/lena.mat zdjęcie].
 +
** Analogicznie jak w przykładzie z kręgosłupem, proszę zaprojektować filtr dolnoprzepustowy, oraz filtr górnoprzepustowy i obejrzeć uzyskiwane obrazy i ich widma.
 +
** Filtr pasmowo-przepustowy można zaprojektować za pomocą funkcji <tt>remez</tt>.
 +
 +
 +
* Proszę pobrać i wczytać do Matlaba [http://www.fuw.edu.pl/~jarekz/SYGNALY/TF/tfr.mat obraz]. Analogicznie jak w przykładzie z kręgosłupem, proszę zaprojektować filtr dolnoprzepustowy, który usunie z obrazu pięć struktur o wysokich częstościach przestrzennych.

Aktualna wersja na dzień 09:31, 22 mar 2016

Laboratorium_EEG/ Wstęp do analizy obrazu

Częstości w obrazie i dwuwymiarowa transformata Fouriera

O analizie obrazów za pomocą Pythona można poczytać tu.

Obrazy możemy analizować tak jakby były to sygnały 2D. Aby uświadomić sobie związki między jednowymiarowym sygnałem a obrazem — sygnałem 2D — zrobimy następujące ćwiczenie:

Generujemy w jednym wymiarze funkcję sinus o 10 okresach na 100 punktów: 

X=100; %pikseli
Fs=1/X; %czestosc probkowania co jeden piksel
dx=1;
x=0:dx:X-dx;
fx=10/X; % czestosc sinusa 10 okresów na X
syg=sin(2*pi*x*fx+pi/3);
figure(1)
plot(x,syg)

Ten obrazek jest znajomy. Teraz przekształcimy go w obraz: 

figure(2)
% teraz zrobimy z niego obraz 2D
Y=X/2;
SYG=zeros(Y,X);

for i=1:Y
    SYG(i,:)=syg;
end


subplot(2,2,1)
imagesc(SYG)
set(gca,'Ydir','normal');
colormap(gray)

Te zmieniające intensywność prążki to właśnie obrazek przedstawiający funkcję sinus w kierunku X. Zobaczmy jak wygląda jego transformata Fouriera: 

subplot(2,2,2)
S=fft2(SYG);
SK_X=-X/2:X/2;
SK_Y=-X/2:X/2;
imagesc(SK_X,SK_Y,abs(fftshift(S)))
set(gca,'Ydir','normal');
colormap(gray)

Oraz ten sam sygnał z dodaną składową stałą i jego transformata: 

subplot(2,2,3)

imagesc(SYG+1)
set(gca,'Ydir','normal');
colormap(gray)

subplot(2,2,4)
S=fft2(SYG+1);
imagesc(SK_X,SK_Y,abs(fftshift(S)))
set(gca,'Ydir','normal');
colormap(gray)
pause

Korelacja i splot dla obrazu

Przyjmijmy konwencję:

  • obraz = duża macierz
  • jądro splotu/korelacji = mniejsza macierz

załóżmy, że mamy obraz A: 

A = [17  24   1   8  15
     23   5   7  14  16
      4   6  13  20  22
     10  12  19  21   3
     11  18  25   2   9];

i jądro korelacji 

h = [8   1   6
     3   5   7
     4   9   2];

wówczas operacja korelacji dla elementu A(2,4) jest następująca:

  1. kładziemy macierz h elementem środkowym (5) na elemencie A(2,4) (czyli 14)
  2. wymnażamy przez siebie wszystkie pokrywające się elementy i sumujemy iloczyny


A(2,4)=1*8 + 8*1 +15*6 +
	   7*3 + 14*5 + 16*7+
	   13*4 + 20*9 +22*2

W matlabie mamy funkcję filter2(h,A), która wykonuje powyższą operację- dwuwymiarową korelację.

Bardzo podobna jest operacja dwuwymiarowego splotu. Jedyna różnica jest taka, że przed operacją mnożenia i dodawania jądro jest obracane o 180 stopni. 

>> rot90(rot90(h))

ans =

     2     9     4
     7     5     3
     6     1     8

Operacja splotu jest zaimplementowana jako conv2

Zarówno operacja splotu jak i korelacji jest kiepsko określona na brzegach obrazu, tam gdzie jądro wystaje poza brzeg. W obu przypadkach wykonywane jest domyślnie dopełnianie obrazu zerami. Efekt ten widać w poniższym przykładzie jako czarne krawędzie. 

I = single(imread('coins.png'));
h = ones(5,5) / 25;
I2 = filter2(h,I);
I3 = conv2(h,I);
subplot(2,2,1)
imagesc(I), title('Pierwotny obraz');
subplot(2,2,2)
imagesc(I2), title('Filtrowany obraz')
subplot(2,2,3)
imagesc(I3), title('Splatany obraz')
colormap('gray')

Dwuwymiarowe filtry FIR

Są naturalnym rozszerzeniem jednowymiarowych filtrów FIR i dają się łatwo zapisać w postaci macierzowej.

Metoda transformacji częstości

Projektujemy filtr jednowymiarowy i transformujemy go na 2D: 

b = remez(10,[0 0.4 0.6 1],[1 1 0 0]);
h = ftrans2(b);
[H,w] = freqz(b,1,64,'whole');
colormap(jet(64))
plot(w/pi-1,fftshift(abs(H)))
figure, freqz2(h,[32 32])

Zprojektujmy i zastosujmy filtr górnoprzepustowy do obrazka: 

b = fir1(14,0.25,'high');
h = ftrans2(b);
[H,w] = freqz(b,1,64,'whole');
colormap(jet(64))
subplot(2,2,1)
plot(w/pi-1,fftshift(abs(H)))
subplot(2,2,2), freqz2(h,[32 32])

load spine 

subplot(2,2,3)
imagesc(X)
colormap bone


subplot(2,2,4)
imagesc(filter2(h,X))

Kilka konkretnych filtrów

A teraz kilka powszechnie stosowanych filtrów. Proszę zbadać funkcję odpowiedzi impulsowej i działanie filtra na nasz przykładowy obrazek.

Filtr uśredniający
1/9	1/9	1/9
1/9	1/9	1/9
1/9	1/9	1/9
Filtr wyostrzający
0	-1	0
-1	5	-1
0	-1	0
Filtr uwypuklający
-1	-1	-1	-1	0
-1	-1	-1	0	1
-1	-1	0	1	1
-1	0	1	1	1
0	1	1	1	1

Filtry nieliniowe

Dotychczas rozważane filtry obliczały lokalną korelację jądra filtru z kolejnymi fragmentami obrazu, czyli wykonywały tam operację liniową. Łatwo sobie wyobrazić, że w kolejnych „okienkach” filtra będziemy wykonywać dowolną inną operację np. medianę i jej wartość wpisywać jako element wynikowego obrazka.

Wczytajmy przykładowy obrazek: 

I = imread('eight.tif');
subplot(2,2,1), imshow(I)
title('oryginalny obraz')

dodajemy szum: 

J = imnoise(I,'salt & pepper',0.02);
subplot(2,2,2), imshow(J)

Filtrujemy obrazek filtrem uśredniającym 

K = filter2(fspecial('average',3),J)/255;
subplot(2,2,3), imshow(K)

A teraz medianowym 

 L = medfilt2(J,[3 3]);
subplot(2,2,4), imshow(L)

Ćwiczenia

  • Proszę pobrać i wczytać do Matlaba zdjęcie.
    • Analogicznie jak w przykładzie z kręgosłupem, proszę zaprojektować filtr dolnoprzepustowy, oraz filtr górnoprzepustowy i obejrzeć uzyskiwane obrazy i ich widma.
    • Filtr pasmowo-przepustowy można zaprojektować za pomocą funkcji remez.


  • Proszę pobrać i wczytać do Matlaba obraz. Analogicznie jak w przykładzie z kręgosłupem, proszę zaprojektować filtr dolnoprzepustowy, który usunie z obrazu pięć struktur o wysokich częstościach przestrzennych.
Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Wstęp_do_analizy_obrazu&oldid=4986”