平滑化2

Matlabにはsmooth3という3次元平滑化を関数Mファイルがあります。 spm_smoothto8bit.mと比較しました。
spm_smoothto8bitはgaussフィルター(正規分布)のFWHM(半値幅)を引数としています。 正規分布は以下の式です。

μとσは平均値と標準偏差です。μがゼロ、σが1のときの正規分布は数のようになります。
str='-exp(-(x.^2)/2)/sqrt(2*pi)'; % 平均0、標準偏差の正規分布×-1
normpdf01=inline(str,'x'); % 正規分布の定義
x=-5:0.1:5;
y=-normpdf01(x); % statictics toolboxがあればy=-normpdf(x,0,1);と等価
figure;set(gcf,'color',[1,1,1]);plot(x,y);

正規分布は上に凸で、平均値のとき最大値を取ることがわかります。
正規分布最大値の半分のときのxを求めます。

検算します。平均値ゼロ、標準偏差1のときのxを求めます。
[xmax,ymax]=fminsearch(str,5) %-1をかけて最小値を求める。
xmax =
0
ymax =
-0.3989
x1=fzero(strcat(str,sprintf('-%0.8f/2',ymax)),0); %ピーク値の半分の幅を求める。
x1=abs(x1-xmax)
x1 =
1.1774
sqrt(log(4))
ans =
1.1774
一致しました。

が正規分布のFWHMになります。標準偏差は1で固定なんらどうか?

次に
doc smooth3
とすると、フィルターがgaussianの場合、[3,3,3]のサイズの標準偏差が0.65であると記載されています。 標準偏差を1固定と仮定すると1サイズあたり、
0.65/3*log(256)^0.5
ans =
0.5102
のfwhmに相当します。
[filename,pathname]=uigetfile('*.img');
fid=fopen(strcat(pathname,filename),'rb');
V=fread(fid,inf,'int16');
fclose(fid);
V=reshape(int16(V),256,256,86); % doulbeよりint16の方が高速
tic;V1=smooth3(V,'gaussian',3);toc % fwhm 1.5
tic;V2=smooth3(V,'gaussian',5);toc % fwhm 2.5
tic;V3=smooth3(V,'gaussian',9);toc % fwhm 4.5
K1=[rot90(squeeze(V(:,128,:))),rot90(squeeze(int16(V1(:,128,:))));...
rot90(squeeze(int16(V2(:,128,:)))),rot90(squeeze(int16(V3(:,128,:))))];
figure;imagesc(K1,[0,384]);colormap(gray);daspect([2,240/256,240/256]);

・・・何か勘違いしているのかもしれませんが、標準偏差も変化しているのかぼかしの効果がはっきりしません。 また、V1,V2,V3の計算時間は4.7970、20.078、120.8900秒でした。
そこでfilterにboxを使ってみました。
tic;V1=smooth3(V,'box',3);toc % simple 3D-convolution?
tic;V2=smooth3(V,'box',5);toc % simple 3D-convolution?
tic;V3=smooth3(V,'box',9);toc % simple 3D-convolution?
K1=[rot90(squeeze(V(:,128,:))),rot90(squeeze(int16(V1(:,128,:))));...
rot90(squeeze(int16(V2(:,128,:)))),rot90(squeeze(int16(V3(:,128,:))))];
figure;imagesc(K1,[0,384]);colormap(gray);daspect([2,240/256,240/256]);

なんかいい感じです。V1,V2,V3の計算時間は4.8120、20.1090、118.8600秒でした。
spm_smoothto8bit.mではどうなるでしょうか?
tic;W1=spm_smoothto8bit(VF,1); % fwhm 1
tic;W2=spm_smoothto8bit(VF,2); % fwhm 2
tic;W3=spm_smoothto8bit(VF,4); % fwhm 4
K2=[rot90(squeeze(V(:,128,:))),rot90(squeeze(int16(W1.dat(:,128,:))));...
rot90(squeeze(int16(W2.dat(:,128,:)))),rot90(squeeze(int16(W3.dat(:,128,:))))];
figure;imagesc(K2,[0,384]);colormap(gray);daspect([2,240/256,240/256]);

W1,W2,W3の計算時間は4.4220、4.6250、4.8900秒でした。
smooth3の3番目の引数と、spm_smoothto8bitの2番目の引数FWHMの換算はどうも違うようです。 演算速度はspm_smoothto8bitの方が圧倒的に速いので、SPMではこちらを採用しているのだと思います。