Matlab画像処理「周波数フィルタ」（2次元FFT）

初投稿です。対戦よろしくお願いします。

画像処理には様々な手法がありますが、この記事で紹介するのは周波数フィルタ処理です。

そもそも知らねーよ！って人はこちらの記事が参考になるかと思います。
https://www.visco-tech.com/newspaper/column/detail20/

matlabで実装できたので、紹介します。

今回実装したのはローパスフィルタです。すなわち高周波ノイズ（点のように狭い領域でのノイズ）を除去します。

出力
左がフィルタ前、右がフィルタ後（伝わりづらい…）

3つの関数にまとめましたので、
%%%ここから下をコピペする%%%以下をコピペすればすぐ使えると思います。

関数の説明

①function [pic2] = FFTfilter(pic1,Gain)
②function[pic2] = OddEsc(pic1)
③function [pic2] =PicConv(pic1)

①は名前の通りFFT計算をする箇所です。
計算の都合上、画像のサイズが奇数だった場合ピクセルを追加して偶数にしています。その処理を関数②に投げています。
周波数領域のトリミングをする際に、画像の象限を反対にすると都合が良いので、③でその処理をしています。

参考：
西住工房様（ https://algorithm.joho.info/image-processing/fourier-transform/ ）

便利なところ
カラー、白黒に関係なく変換できます。
FFTGain（0~1）の調整（0に近いほどフィルタが強くなります。）が出来ます。

欠点
ハイパスフィルタ、バンドパスフィルタの能力はないです。
もし追加するのであれば

if i < sz1/2-a || i > sz1/2+a || j < sz2/2-b || j > sz2/2+b

の部分が周波数トリミングに対応しているのでここの条件を変えてください。

画像サイズが奇数だったときに追加する処理を行うので右端（あるいは下端）の1列に輝度ゼロの領域が生まれることがあります。

以下サンプルコード

%%%画像の読み込み%%%
test = imread("C:\hogehoge\test.jpg");

%%%フィルタ処理(0.1はフィルタ透過率)%%%
FFTGain = 0.1 
test2 = FFTfilter(test,FFTGain);

%%%出力%%%
subplot(1,2,1)
imshow(test)
subplot(1,2,2)
imshow(test2)
imwrite(test2,"C:\hogehoge\test2.jpg")


%%%ここから下をコピペする%%%
function [pic2] = FFTfilter(pic1,Gain)
pic1 = OddEsc(pic1);

sz1 = size(pic1,1);
sz2 = size(pic1,2);
sz3 = size(pic1,3);

ftest=fft2(pic1);
ftest=PicConv(ftest);
a = ceil(sz1*Gain);
b = ceil(sz2*Gain);

for k =1:sz3
    for j=1:sz2
        for i=1:sz1
            if i < sz1/2-a || i > sz1/2+a || j < sz2/2-b || j > sz2/2+b
            ftest(i,j,k) = 0;
            end
        end
    end
end

ftest = PicConv(ftest);
pic2 = real(ifft2(ftest))/255;


end

function[pic2] = OddEsc(pic1)
sz1 = size(pic1,1);
sz2 = size(pic1,2);
sz3 = size(pic1,3);

if rem(sz1,2)==1
    sz1 = sz1+1;
    pic1(sz1,sz2,sz3) = 0;
end

if rem(sz2,2)==1
    sz2 = sz2+1;
    pic1(sz1,sz2,sz3) = 0;
end

pic2 = pic1;

end

function [pic2] =PicConv(pic1)
a = size(pic1,1);
b = size(pic1,2);
c = size(pic1,3);
pic2 = zeros(a,b,c);

pic2(1:a/2,1:b/2,:) = pic1(a/2+1:a,b/2+1:b,:);
pic2(a/2+1:a,1:b/2,:) = pic1(1:a/2,b/2+1:b,:);
pic2(1:a/2,b/2+1:b,:) = pic1(a/2+1:a,1:b/2,:);
pic2(a/2+1:a,b/2+1:b,:) = pic1(1:a/2,1:b/2,:);

end

干渉計の干渉縞のノイズ除去などに使ってみてください。それでは。