空間周波数の操作３（周波数解析のグラフ）

function spectrum3

close all;

clear all;

global MYFILENAME; % 別ファイルでも使うのでグローバル変数

imgFileName = 'grayRabbit.jpg'; % 画像のファイル名

im1 = imread(imgFileName); % 画像のデータを読み込む。

% 画像のファイル名を分解

[pathstr, MYFILENAME, ext] = fileparts(imgFileName);

%% 必要に応じて正規化

im2 = double(im1);

imMean=mean(im2(:)); % 平均輝度

imSD=std(im2(:)); % 輝度の標準偏差

im2=(im2-imMean)./imSD; % 平均を０、標準偏差を１に。

im2 = im2.*25; % 標準偏差の変更（25は適当）

im2 = im2+128; % 平均値の変更（128は256階調の半分）

% 正規化しないときは次の行をコメントアウト

%im2 = im1;

% 正規化した画像をjpgファイルとして出力。

imwrite(uint8(im2), [MYFILENAME '-norm.jpg'], 'jpeg');

%% 座標の設定

imSize = size(im1);

halfImgSize = imSize(2)/2; % 画像が正方形、つまり、imSize(1) = imSize(2)を前提としている。

[X, Y] = meshgrid(-halfImgSize:(halfImgSize-1), halfImgSize:-1:-(halfImgSize-1)); % 行がY方向、列がX方向であることに注意。Y軸は上が正であることに注意。

[TH, R] = cart2pol(X,Y); % X,Yを極座標へ変換

%% 出力結果を表示するfigureの設定

ScreenSize = get(0, 'ScreenSize');

figure('Position', [ScreenSize(3)./4, ScreenSize(4)./4, 1000, 800]); % 画像の位置

%% 正規化した画像を表示

subplot(2,3,1);

ShowGrayScale(im2);

%% 画像データにハニング窓（四角形）をかける    

absXY = max(abs(X), abs(Y)); 

% absXYは0からhalfImgSizeを取り得る。

% 画像の中心部では、absXYは0で、hanningの値は1となる。

% absXYが最大(halfImgSize）のとき、hanningの値は0となる。

hanning1 = 0.5 - 0.5 * cos(2*pi*(halfImgSize/2+absXY./2)/halfImgSize);

im3 = double(im2) .* hanning1;

subplot(2,3,2);

ShowGrayScale(im3);

%% 画像データにハニング窓（円形）をかける

R2 = min(R, halfImgSize); % halfImgSizeを超えない半径のデータを作る。

% R2は0からhalfImgSizeを取り得る。

% 画像の中心部では、R2は0で、hanningの値は1となる。

% R2が最大(halfImgSize）のとき、hanningの値は0となる。

hanning2 = 0.5 - 0.5 * cos(2*pi*(halfImgSize/2+R2./2)/halfImgSize);

im4 = double(im2) .* hanning2;

subplot(2,3,3);

ShowGrayScale(im4);

%% im2のFFT

Fim2 = fftshift(fft2(im2));

% スペクトル画像

subplot(2,3,4);

ShowGrayScale(log(abs(Fim2)));

%ShowGrayScale(log(Fim2.*conj(Fim2))); % こういう方法もある。

%% im3のFFT

Fim3 = fftshift(fft2(im3));

% スペクトル画像

subplot(2,3,5);

ShowGrayScale(log(abs(Fim3)));

%% im4のFFT

Fim4 = fftshift(fft2(im4));

% スペクトル画像

subplot(2,3,6);

ShowGrayScale(log(abs(Fim4)));

% 周波数解析

ShowSpectrum(R, Fim4);

function ShowGrayScale(imageData)

% 画像をグレースケールで表示する。

imagesc(imageData);

axis('image');

axis('off');

colormap(gray);

function ShowSpectrum(Rdata, Fdata)

% スペクトル解析のグラフ（横軸は周波数（半径）cycle/image）を描く

% Rdataは、画像中心からの半径のデータ

% Fdataは、画像データをフーリエ変換した結果（RdataとFdataのサイズは同じ。ただしFdataは複素数）

global MYFILENAME;

% 中心からの半径の最大値

Rmax = floor(size(Rdata, 1)/2);

Intensity = zeros(1,Rmax); % 半径ごとのスペクトルの合計値を保管

outputData = zeros(Rmax, 2); % CSVファイルに出力するデータ

for n = 1:Rmax;

    mask = (Rdata > n-1 & Rdata <= n); % ある半径の1ピクセル以内のデータだけを取り出す。

    % （パワー）スペクトルの計算

    spect = mask.*abs(Fdata);

    %spect = mask.*(abs(Fdata).^2); % パワースペクトル

    %spect = mask.*(Fdata.*conj(Fdata)); % パワースペクトル（上の式と結果は同じ）

    ind = find(mask); % indはmask=1の場所を表す列ベクトル

    Intensity(n) = sum(sum(spect)) / size(ind, 1); % 正規化（spectの数は半径ごとに異なるため）

    outputData(n, :) = [n Intensity(n)];

end

figure;

loglog(1:Rmax,Intensity,'o'); % 対数軸でのグラフ

xlabel('Spatial Frequency (cycle per image)'), ylabel('Fourier Power');

title(MYFILENAME);

% すべてのパワースペクトルのデータを出力

dlmwrite([MYFILENAME '-ps1.csv'], outputData, 'precision', 20); % 有効桁数に注意