空間周波数フィルタリング

function imageFiltered = SpatialFilter(image, cutoff, pass)

% image: square array
% cutoff: cycles per image
% pass: 1) high-pass. 2) low-pass.


%% FFT
dc = mean(image(:));
image = image – dc;
imageFFT = fft2(image);

%% Power spectrum
powerS = log10(abs(fftshift(imageFFT)));
figure;
imagesc(powerS, [min(powerS(:)) max(powerS(:))]);
colormap(gray(256));
axis image;

%% Spatial Filter
[m, n] = size(image);
[x, y] = meshgrid(1:m, 1:n);
if pass == 1
filter = ((x-((m+1)/2)).^2 + (y-((n+1)/2)).^2) > cutoff^2;
else
filter = ((x-((m+1)/2)).^2 + (y-((n+1)/2)).^2) < (cutoff+1)^2;
end
filter = ifftshift(filter);

%% Filtering and inverse FFT
imageFilteredFFT = imageFFT.*filter;
imageFiltered = ifft2(imageFilteredFFT, ‘symmetric’);
imageFiltered = imageFiltered + dc;


関数は以上です。Psychtoolbox 関数は使用していません。フィルタをかける画像が必要なので、ここでは cameraman.tif を読み込んで、変数の型を double に変更します。
image = imread(‘cameraman.tif’);
image = cast(image, ‘double’);
imagesc(image, [0 255]);
colormap(gray);
axis image;

imageFiltered = SpatialFilter(image, 30, 1); を実行すると、配列変数 image に対して 30 cycles per image 未満の空間周波数をカットするハイパスフィルタをかけます(30 cycles はカットされません)。また、SpatialFilter() はパワースペクトルを表示します。

出力されたハイパス画像を表示するには以下を実行します。
imagesc(imageFiltered, [min(imageFiltered(:)) max(imageFiltered(:))]);
colormap(gray(256));
axis image;

元画像は256階調のグレースケールですが、フィルタをかけた画像の最小値は −31.7 、最大値は 301 になっていることに注意してください。実験刺激をフィルタリングする場合、フィルタされた画像の値が出力範囲(例えば0~255)におさまっている必要があります。

imageFiltered = SpatialFilter(image, 30, 2); を実行すると、配列変数 image に対して 30 cycles per image を超える空間周波数をカットするローパスフィルタをかけます(30 cycles はカットされません)。

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s