Wstęp do analizy obrazu
Laboratorium_EEG/ Wstęp do analizy obrazu
Spis treści
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:
- 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
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.
