Komputerowe Przetwarzanie Obrazów

Szybka Transformata Fouriera

1. Generowanie sygnałów 1D o różnych częstotliwościach oraz dodawanie szumu.

t = 0:0.001:2;
x1 = sin( 3*pi*2*t );
plot(t,x1)
title('Czestotliwosc rowna 3 Hz')
xlabel('Czas w sekundach')

pause
x2 = sin( 5*pi*2*t );
figure
plot(t,x2)
title('Czestotliwosc rowna 5 Hz')
xlabel('Czas w sekundach')

pause
x = x1+x2;

figure
plot(t,x)
title('Polaczone sygnaly o czestotliwosciach 3 oraz 5 Hz')
xlabel('Czas w sekundach')

pause

x = x + 0.8*(1 - 2*rand(size(x)));
plot(t,x)
title('Polaczone sygnaly z dodanym szumem')
xlabel('Czas w sekundach')

pause
close all

2. Analiza częstotliwości występujących w sygnale

t = 0:0.001:2;
x1 = sin( 25*pi*2*t );
plot(t,x1)
title('Czestotliwosc rowna 25 Hz')
xlabel('Czas w sekundach')
pause

x2 = sin( 40*pi*2*t );
figure
plot(t,x2)
title('Czestotliwosc rowna 40 Hz')

xlabel('Czas w sekundach')

pause
x = x1+x2;

figure
plot(t,x)
title('Polaczone sygnaly o czestotliwosciach 25 oraz 40 Hz')
xlabel('Czas w sekundach')

pause

Y = fft(x,512);
Y2 = Y.* conj(Y) / 512;

%tylko pierwsza polowa wspolczynników jest znacząca
f = 1000*(0:256)/512;
plot(f,Y2(1:257))
title('Frequency content of x')
xlabel('frequency (Hz)')

pause
close all

3. Jednowymiarowa transformata Fouriera i odwrotna transformata Fouriera

len = 512
t = 0:len;
x = sin(0.08*t)+sin(0.03*t);
plot(t,x);
title('Oryginalny sygnal')
xlabel('Czas')

pause

y=fft(x)
Y2 = y.* conj(y);
figure
plot(Y2(1:len/2))
title('Wspolczynniki rozwiniecia za pomoca FFT')
xlabel('Czestotliwosc w Hz')

pause
x2 = ifft(y);

figure
plot(t,x2)
title('Sygnal odzyskany po zastosowaniu odwrotnej transformaty Fouriera')
xlabel('Czas')

pause
close all

4. Usuwanie szumów z sygnału 1D za pomocą FFT (Fast Fourier Transform). Zwróć uwagę na

zniekształcenia na lewym i prawym skraju odfiltrowanego sygnału.

len = 512
t = 0:len;
x = sin(0.08*t)+sin(0.03*t) + sin(0.2*t);
plot(t,x);
title('Oryginalny sygnal')
xlabel('Czas w sekundach')

pause
x = x + 2*(-1 + 2*rand(size(x)));
figure
plot(t,x);
title('Oryginalny sygnal z dodanym szumem')
xlabel('Czas')
pause

y=fft(x)

Y2 = y.* conj(y);
figure
plot(Y2(1:len/2))
title('Wspolczynniki FFT')
xlabel('Czestotliwosc [Hz]')

pause
max(Y2)
min(Y2)
mask = Y2 > 0.2*max(Y2);

Y2 = Y2 .* mask;
y = y .* mask;
figure
plot(Y2(1:len/2))
title('Wspolczynniki FFT po usunieciu zbyt slabych czestotliwosci')
xlabel('Czestotliwosc [Hz]')

pause
x2 = ifft(y);
figure
plot(t,x2)
title('Odzyskany sygnal')
xlabel('Czas')
pause; close all;

5. Dwuwymiarowa FFT oraz IFFT (Inverse Fast Fourier Transform). Zwróć uwagę, że rysunek

transformaty pokazuje większą energię dla dużych horyzontalnych częstotliwości niż dla
dużych wertykalnych częstotliwości. Wynika to z faktu, iż na oryginalnym rysunku szybsze
zmiany następują w przekroju poziomym (tzn. jest większa częstotliwość zmian), dlatego, że
prostokąt jest węższy w poziomie.

size = [32 32]
%size = [512 512]
x = zeros(size);
x(size(1)/4:3*size(1)/4, 3*size(2)/8 : 5*size(2)/8) = ones(size(1)/2+1, (5/8 - 3/8)*size(2)+1);
imshow(x,'notruesize')
title('Obraz oryginalny')

pause
F = fft2(x);
F2 = log(abs(F));
figure
imshow(F2,[-1 5],'notruesize'); colormap(jet); colorbar
title('FFT')
%size(F2)

pause
F = fft2(x,256,256);
F2 = log(abs(F));
figure
imshow(F2,[-1 5],'notruesize'); colormap(jet); colorbar
title('FFT z uzupelnieniem zerami do rozmiaru 256x256')
%size(F2)

pause
%przesuniecie ukladu wspolrzednych tak, aby czestotliwosc 0 byla w srodku
F2 = fftshift(F);
F2 = log(abs(F2));
figure
imshow(F2,[-1 5],'notruesize'); colormap(jet); colorbar
title('FFT i przesuniecie ukladu wspolrzednych tak, aby czestotliwosc 0 byla w srodku')

pause
F = fft2(x);
x2 = ifft2(F);
figure
imshow(real(x2),'notruesize')
title('Odwrotna FFT')

pause
close all

6. Przykłady działania FFT na różnych obrazach

x = imread('fig1.bmp');
x = rgb2gray(x);
%x = imread('portret.jpg');
x = double(x)/255;

imshow(x,'notruesize')
title('Obraz oryginalny')

pause
F = fft2(x);
F2 = log(1 + abs(F));
minn = min(min(F2))
maxx = max(max(F2))
figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('FFT')

pause
F = fft2(x,512,512);
F2 = log(1 + abs(F));
figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('FFT z uzupelnieniem zerami do rozmiaru 512x512')

pause
%przesuniecie ukladu wspolrzednych tak, aby czestotliwosc 0 byla w srodku
F2 = fftshift(F);
F2 = log(1+abs(F2));

figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('FFT i przesuniecie ukladu wspolrzednych tak, aby czestotliwosc 0 byla w srodku')

pause
F = fft2(x);
x2 = ifft2(F);
figure
imshow(real(x2),'notruesize')
title('Odwrotna FFT')

pause
close all

7. Operacje w przestrzeni częstotliwości – usuwanie częstotliwości zbyt małych, zbyt dużych, itd.

x = imread('fig1.bmp');
x = rgb2gray(x);

%x = imread('portret.jpg');
x = double(x)/255;

imshow(x,'notruesize')
title('Obraz oryginalny')

pause
F = fft2(x);
F2 = log(1 + abs(F));
minn = min(min(F2))
maxx = max(max(F2))
figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('FFT')

pause
F = fft2(x,512,512);
F2 = log(1 + abs(F));
figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('FFT z uzupelnieniem zerami do rozmiaru 512x512')

pause
%przesuniecie ukladu wspolrzednych tak, aby czestotliwosc 0 byla w srodku
F2 = fftshift(F);
F2 = log(1+abs(F2));

figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('FFT i przesuniecie ukladu wspolrzednych tak, aby czestotliwosc 0 byla w srodku')

pause
F = fft2(x);
F2 = fftshift(F);
F2 = log(1 + abs(F2));
minn = min(min(F2))
maxx = max(max(F2))
mask = F2 >= 0.07*maxx;
%mask2 = F2 <= 0.9*maxx;
%mask = mask & mask2;
F = F .* mask;
F2 = F2 .* mask;

figure
imshow(F2,[minn maxx],'notruesize'); colormap(jet); colorbar
title('Usuniecie zbyt malych czestotliwosci')

pause
x2 = ifft2(F);

figure
imshow(real(x2),'notruesize')
title('Odwrotna FFT')

pause
close all

8. Konwolucja za pomocą FFT.

clear;
x = imread('portret.jpg');
s = size(x)
x = imread('portret.jpg');
x = double(x)/255;
imshow(x,'notruesize')
title('Obraz oryginalny')

pause
mask = [-1,-2,-1; 0, 0, 0; 1, 2, 1];
mask = rot90(mask,2);
mask(s(1),s(2)) = 0; % Zero-pad mask to be 8-by-8;

%KONWOLUCJA
x2 = ifft2(fft2(x).*fft2(mask));
x2 = real(x2); % Remove imaginary part caused by roundoff error

figure
imshow(x2,'notruesize');
title('po konwolucji');

pause

mask = [-1,-2,-1; 0, 0, 0; 1, 2, 1];
x3 = filter2(mask, x);
%x3 = mat2gray(x3);
figure
imshow(x3,'notruesize');
title('po filtrowaniu');

pause
close all

9. TEMPLATE MATCHING: Korelacja za pomocą FFT - „feature detection” na przykładzie

lokalizacji literki „a”.

WERSJA 1:

bw = imread('text.tif');
a=bw(59:71,81:91); %Extract one of the letters "a" from the image.
%bw = rot90(bw);

imshow(bw);
figure, imshow(a);
figure, imshow(rot90(a,2));

s = size(bw)
s2 = size(a)

pause
C = real(ifft2(fft2(bw) .* fft2(rot90(a,2),256,256)));
figure, imshow(C,[])%Display, scaling data to appropriate range.
maxx = max(C(:)) %Find max pixel value in C.

%thresh = 45; %Use a threshold that's a little less than max.
thresh = 0.9*maxx; %Use a threshold that's a little less than max.
figure, imshow(C > thresh)%Display showing pixels over threshold.

pause
[r c] = find(C > thresh);
Z = zeros(size(bw));
for i =1:size(r)
disp(sprintf('%d %d',r(i),c(i)));
Z(r(i):r(i)+s2(1)-1, c(i):c(i)+s2(2)-1) = a;
end

figure
imshow(Z,[]);
pause
close all

WERSJA 2:

bw = imread('text.tif');
a=bw(59:71,81:91); %Extract one of the letters "a" from the image.
bw = rot90(bw);
imshow(bw);
figure, imshow(a);
figure, imshow(rot90(a,3));

s = size(bw)
s2 = size(a)

pause
C = real(ifft2(fft2(bw) .* fft2(rot90(a,3),256,256)));
figure, imshow(C,[])%Display, scaling data to appropriate range.
maxx = max(C(:)) %Find max pixel value in C.

%thresh = 45; %Use a threshold that's a little less than max.
thresh = 0.9*maxx; %Use a threshold that's a little less than max.
figure, imshow(C > thresh)%Display showing pixels over threshold.

pause
[r c] = find(C > thresh);
Z = zeros(size(bw));
for i =1:size(r)
disp(sprintf('%d %d',r(i),c(i)));
Z(r(i):r(i)+s2(2)-1, c(i):c(i)+s2(1)-1) = rot90(a,1);
end

figure
imshow(Z,[]);
pause
close all