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ć 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.