平滑化
spm_smoothto8bit.mが平滑化を行う関数Mファイルです。
どの程度平滑化されているかMatlabで検討してみます。
まずoriginalのデータを読みます。SPMのnormaliseを実行したときに保存される*sn.matを読みます。
[filename,pathname]=uigetfile('*sn.mat'); % sn.matファイルを選択
load(strcat(pathname,filename)); % *sn.matファイルを読み込む
データが読み込まれました。但し、ここには実際の三次元配列データはありません。
関数spm_read_volsを使って三次元配列を取得します。
VF.dat=int16(spm_read_vols(VF)); % 三次元配列を取得 三次元2byte整数配列化
figure;colormap(gray);
subplot(121);X=rot90(squeeze(VF.dat(round(VF.dim(1)/2),:,:))); % 画像を90度回転
imagesc(X,[0,384]);daspect(abs([VF.mat(3,3),VF.mat(2,2),VF.mat(1,1)])); % 縦横比補正
subplot(122);X=rot90(squeeze(VF.dat(:,round(VF.dim(2)/2),:)));
imagesc(X,[0,384]);daspect(abs([VF.mat(3,3),VF.mat(2,2),VF.mat(1,1)]));
h=findobj(gcf,'type','axes');set(h,'visible','off'); %目盛りを消去
set(h(2),'position',[0,0,0.5,1]);set(h(1),'position',[0.5,0,0.5,1]);
次にspm_smoothto8bit.mを使います。
help spm_smoothto8bit
とすると使い方が表示されます。書式は
平滑後データ=spm_smoothto8bit(平滑前データ,GaussフィルターのFWHM)
だそうです。
FWHMはfull-width at half-maximumの略です。日本語では半値幅と呼ばれ、
ピークの高さの半分における幅を示すものだそうです。
FWHMの初期値は8mmで、この値はflags.smosrcに保存されています。
clear X h;
VF=rmfield(VF,'dat'); %このフィールドがあるとerrorがでます
V0=spm_smoothto8bit(VF,flags.smosrc);
V0という変数ができました。
V0.datなる正の整数1byteからなる３次元配列とV0.pinfが新たな変更点です。
V0.datを見てみます。
figure;imagesc(rot90(squeeze(V0.dat(:,128,:))),[0,1.5*256]);
colormap(gray);daspect([2,240/256,240/256]);
平滑化されています。半値幅をいろいろ変えてみます。
V1=spm_smoothto8bit(VF,1);
V2=spm_smoothto8bit(VF,2);
V4=spm_smoothto8bit(VF,4);
V8=V0;
V16=spm_smoothto8bit(VF,16);
V32=spm_smoothto8bit(VF,32);
X=[rot90(squeeze(V1.dat(:,:,50))),rot90(squeeze(V2.dat(:,:,50))),...
rot90(squeeze(V4.dat(:,:,50)));rot90(squeeze(V8.dat(:,:,50))),...
rot90(squeeze(V16.dat(:,:,50))),rot90(squeeze(V32.dat(:,:,50)))];
imagesc(X,[0,1.5*256]);daspect([240/256,240/256,2]);
1,2,4,8,16,32mmの半値幅です
。半値幅が大きくなるに従い、ぼかしが強くなっているのがわかります。
spm_smmothto8bit.mは三次元配列の値をスライス毎に変化させています。
変化させた情報はpinfoに書かれています。
被験者脳の三次元配列と、半値幅1mmのGaussフィルター処理後の三次元配列で比較してみます。
V=spm_read_vols(VF); % 被験者脳データを再読み込み
X=double(squeeze(V(:,128,:))); %被験者脳
Y=double(squeeze(V1.dat(:,128,:)));% Gaussフィルター処理後
figure;imagesc([rot90(X),rot90(Y)],[0,384]);
colormap(gray);daspect([2,240/256,240/256]);
スライス毎に縞々模様が見えます。pinfoを使って復元します。
Z=Y*diag(V1.pinfo(1,:))+ones(size(Y,1),1)*V1.pinfo(2,:);
figure;imagesc([rot90(X),rot90(Z)],[0,384]);
colormap(gray);daspect([2,240/256,240/256]);
ほぼ復元されています。Gaussフィルター前後の差分をみてみます。
ZZ=X-Z;
figure;imagesc([rot90(ZZ+192),rot90(ZZ*5+192)],[0,384]);
colormap(gray);daspect([2,240/256,240/256]);
右図は差分を５倍にして表示したものです。
尚、spm_normalise.mの中ではこのpinfoは平均値の1/8以下のデータは切り捨てた値で再補正されています。
spm_global.mなるもので、割り算してありますが、
x=rand(3,5,7);
spm_global(x)
と
x=x(:);
mean(x(x>mean(x)/8)
は等価です。
半値幅とpinfoの関係を見てみます。
K1=[V1.pinfo(1,:);V2.pinfo(1,:);V4.pinfo(1,:);...
V8.pinfo(1,:);V16.pinfo(1,:);V32.pinfo(1,:)];
K2=[V1.pinfo(2,:);V2.pinfo(2,:);V4.pinfo(2,:);...
V8.pinfo(2,:);V16.pinfo(2,:);V32.pinfo(2,:)];
figure;set(gcf,'color',[1,1,1]);
subplot(121);plot(K1');title('scale');xlim([1,86]);
subplot(122);plot(K2');title('offset');xlim([1,86]);
半値幅が小さいほど、スライス毎の縞々模様が目立つことがわかります。
たぶんspm_smoothto8bit.mは三次元配列に任意の半値幅の正規分布でx方向、y方向、z方向にぼかしをいれる関数なのだと思います。