%%%% ロボットビジョン(5T) sample05_filter.m  2015/5/28
clear;
Iin = imread('image01.bmp');
I = double(Iin);

Idft = fftshift(fft2(I));   % 2次元FFT後，DCが中心に来るようにシフト

% 入力画像のスペクトル
Mag_Idft = abs(Idft); 
[n,m] = size(Mag_Idft);
for j = 1 : n
    for i = 1 : m
        if Mag_Idft(j,i) > 50000;  % 表示のため，適当な値で切る
            Mag_Idft(j,i) = 50000;
        end;
    end;
end;

% 円形ローパスフィルタ-----------------------------------------
Ilow = Idft;
for j = 1 : n
    for i = 1 : m
        r = sqrt((j-m/2)^2 + (i-n/2)^2);
        if r > 30  %円形ローパスフィルタのサイズ（半径） 
            Ilow(j,i) = 0;
        end;
    end;
end;
Ilow_idft = ifft2(fftshift(Ilow));  % 2次元逆FFT
% 円形ローパスフィルタ後の画像のスペクトル
Mag_Ilow = abs(Ilow); 
for j = 1 : n
    for i = 1 : m
        if Mag_Ilow(j,i) > 50000;  % 表示のため，適当な値で切る
            Mag_Ilow(j,i) = 50000;
        end;
    end;
end;


% 円形ハイパスフィルタ-----------------------------------------
Ihi = Idft;
for j = 1 : n
    for i = 1 : m
        r = sqrt((j-m/2)^2 + (i-n/2)^2);
        if r < 40 %円形ハイパスフィルタのサイズ（半径）
            Ihi(j,i) = 0;
        end;
    end;
end;
Ihi_idft = ifft2(fftshift(Ihi));  % 2次元逆FFT
% 円形ハイパスフィルタ後の画像のスペクトル
Mag_Ihi = abs(Ihi); 
for j = 1 : n
    for i = 1 : m
        if Mag_Ihi(j,i) > 50000;  % 表示のため，適当な値で切る
            Mag_Ihi(j,i) = 50000;
        end;
    end;
end;


% ガウス型ローパスフィルタ-----------------------------------------
Igl = Idft;
pi = 3.1415;
sigma = 20;   % ガウス関数の標準偏差
for j = 1 : n
    for i = 1 : m
        gaussian(j,i) = exp(-((i-m/2)^2+(j-n/2)^2)/(2*sigma^2));
        Igl(j,i) = Igl(j,i) * gaussian(j,i);
    end;
end;
Igl_idft = ifft2(fftshift(Igl));  % 2次元逆FFT
% ガウス型ローパスフィルタ後の画像のスペクトル
Mag_Igl = abs(Igl); 
for j = 1 : n
    for i = 1 : m
        if Mag_Igl(j,i) > 500;  % 表示のため，適当な値で切る
            Mag_Igl(j,i) = 500;
        end;
    end;
end;

figure(1);
subplot(2,4,1); %画像を2x4並べて表示
imagesc(I); colormap(gray);
subplot(2,4,2);
imagesc(Mag_Idft); colormap(gray);

subplot(2,4,3);
imagesc(abs(Ilow_idft)); colormap(gray);
subplot(2,4,4);
imagesc(Mag_Ilow); colormap(gray);

subplot(2,4,5);
imagesc(abs(Ihi_idft)); colormap(gray);
subplot(2,4,6);
imagesc(Mag_Ihi); colormap(gray);

subplot(2,4,7);
imagesc(abs(Igl_idft)); colormap(gray);
subplot(2,4,8);
imagesc(Mag_Igl); colormap(gray);

figure(2); % ガウス関数表示
mesh(gaussian);

%%%% EOF