SPMの関数

備忘録です。
関数機能
spm_defaults.m初期値設定
spm_dicom_headers.mDICOM header取得
spm_matrix.mAffine行列作成
spm_normalise.m標準脳化
spm_read_vols.m三次元配列取得
spm_sample_vol.dll三次元補間
spm_segment.m脳切り出し
spm_slice_vol.dllAffine変換後のslice画像取得
spm_smoothto8bit.m平滑化
spm_vol.mAnalyseまたはMNICファイルのheader取得
spm_write_sn.m標準脳化後データ出力
spm_write_vol.mAnalyse形式で保存

spm_defaults.m

書式 spm_defaults
SPMの初期値を設定する

実行例
spm_defaults;

spm_dicom_headers.m

書式 hdr=spm_dicom_headers(P) Matlab Image Processing Toolboxのdicominfo関数と同じ機能
中身を見るにはhdr{1}などとする

例 DICOMのheader取得
[filename1,pathname1]=uigetfile('*'); % DICOM file 1
[filename2,pathname2]=uigetfile('*'); % DICOM file 2
P(1,:)=strcat(pathname1,filename1);
P(2,:)=strcat(pathname2,filename2);
hdr=spm_dicom_headers(P);

spm_matrix.m

書式 A=spm_matrix(P) 以下の順で座標変換行列を作成する
軸向き変化→拡大縮小→Z軸回転→Y軸回転→X軸回転→平行移動

例 Affine変換行列の作成例
P(1:3)=[7,11,13];%x,y,zを7,11,13平行移動
P(4:6)=[30,45,70]/180*pi;%x,y,z軸を30,45,70度回転
P(7:9)=[240/256,240/256,2];% x,y,zのvoxwlの大きさは240/256,240/256,2mm
A=spm_marix(P);

spm_normalise.m

書式 prm=spm_normalise(VG,VF,matname,VWG,VWF,flags) *sn.mat内の構造体変数VFやVGでも可
VWF,VWGがない場合は''とする
平滑化→Affine変換→非線形変換という手順を経る

例 SPMの初期設定値で標準脳化
spm_defaults; % 初期値設定
[filename,pathname]=uigetfile('*.mnc'); % 標準脳ファイル名
VG=strcat(pathname,filename);
[filename,pathname]=uigetfile('*.img'); % 被験者脳ファイル名
VF=strcat(pathname,filename);
matname=strcat(pathname,'abc.mat'); % 保存ファイル名
prm=spm_normalise(VG,VF,matname,'','',defaults.normalise.estimate);

spm_read_vols.m

書式 VO=spm_read_vols(V,mask) VF.datを含まないこと

例 被験者解剖脳の読み込み
[filename,pathname]=uigetfile('*sn.mat');
load(strcat(pathname,filename)); % *sn.mat読み込み
VO=spm_spm_read_vols(VF);
低レベルのデータ読み込み
[filename,pathname]=uigetfile('*.img'); % *.imgファイル選択
fid=fopen(strcat(pathname,filename),'rb'); % *.imgを開く
V=fread(fid,inf,'int16'); % 全データを2byte整数と仮定して読み込み
fcloes(fid); % *.imgを閉じる
V=int16(V); % 2byte整数に変換
V=reshape(256,256,86); % 三次元配列化

spm_sample_vol.dll

書式 [VO,X,Y,Z]=spm_sample_vol(V,x,y,z,hold) ndgridかmeshgridが必要
reshapeで配列を組みなおす必要がある
interp3では大規模三次元配列の補間は不可能

例 5x7x11の配列を10x14x22の配列に線形補間
V=rand(5,7,11);
x=linspace(1,size(V,1),size(V,1)*2);
y=linspace(1,size(V,2),size(V,2)*2);
z=linspace(1,size(V,3),size(V,3)*2);
[x,y,z]=ndgrid(x,y,z);
VO=spm_sample_vol(V,x,y,z,1);
VO=reshape(VO,size(x));

spm_segment.m

書式 VO=spm_segment(PF,PG,flags); 解剖画像をAffine変換し標準脳化した後、灰白質・白質・髄液に分画する。
それぞれのtemplateは.../spm2/apriori/*.mnc
非線形変換は行われない。

例 灰白質・白質・髄液に分画
[filename,pathname]=uigetfile('*.mnc'); %標準脳template
PG=strcat(pathname,filename);
[filename,pathname]=uigetfile('*.img'); %解剖画像
PF=strcat(pathname,filename);
spm_defaults;
VO=spm_segment(PF,PG,defaults.segment.estimate);

spm_slice_vol.dll

書式 X=spm_slice_vol(V,A,dim,hold) Affine変換後、0スライス目の[0,0]〜[dim(1),dim(2)]の画像が出力される。
Aを調整してスライス枚数を変化させる必要がある。
冠状断、矢状断を得るにもAを変化させる

例 5x7x11の配列を10x14x22の配列に線形補間
[filename,pathname]=uigetfile('*sn.mat');
load(strcat(pathname,filename));
VF.dat=spm_read_vols(VF);
K=eye(4);K(1:3,1:3)=eye(3)*2; % voxelを1mm単位に
K(3,4)=-VF.dim(3)*2/2; % Z軸中央部分の水平断面
XYimg=spm_slice_vol(VF.dat,Affine*inv(K),VG.dim([1,2])*2,1);
K=eye(4);K(1:3,1:3)=eye(3)*2; % 1mm単位に
K(2:3,2:3)=rot90(K(2:3,2:3)); % YとZを変換→冠状断
K(3,4)=-VG.dim(2); %X軸中央に
YZimg=spm_slice_vol(VF.dat,Affine*inv(K),VG.dim([1,3])*2,1);
K=eye(4);K(1:3,1:3)=eye(3)*2; % 1mm単位に
K(1:3,1:3)=rot90(K(1:3,1:3)); % XとZを変換→矢状断
K(3,4)=-VG.dim(1); % X軸中央に
ZYimg=spm_slice_vol(VF.dat,Affine*inv(K),VG.dim([3,2])*2,1);

spm_smoothto8bit.m

書式 VO=spm_smoothto8bit(V,fwhm) VF.datを含まないこと
fwhmのdefaultはflags.smosrcの8

例 8mmのGaussフィルターで平滑化
[filename,pathname]=uigetfile('*sn.mat');
load(strcat(pathname,filename)); % *sn.mat読み込み
VO=spm_smoothto8bit(VF,flags.smosrc);

spm_vol.m

書式 V=spm_vol(P) 三次元データは含まれない
analyse formatを読むspm_vol_ana.mと minc formatを読むspm_vol_minc.mが本体
更にspm_vol_minc..はspm_read_hdr.mが、spm_vol_minc.はspm_read_netcdf.mが本体。
三次元データの読み込みはspm_read_volsを使う

例 *.imgのデータ読み込み
[filename,pathname]=uigetfile('*.img');
V=spm_vol(strcat(pathname,filename));
V.dat=spm_read_vols(V); % 三次元データ取得

例 *.imgのデータ読み込み
[filename,pathname]=uigetfile('*.img');
V=spm_vol_ana(strcat(pathname,filename));
V.dat=spm_read_vols(V); % 三次元データ取得

spm_write_sn.m

書式 msk=spm_write_sn(V,prm,flags,'mask') 機能画像の三次元配列の中心と解剖画像の三次元配列の中心を一致させておくこと
非線形変換の情報prm.Trがないと線形変換のみ実行される

例 標準脳化後データの取得
spm_defaults; % 初期値設定
[filename,pathname]=uigetfile('*.mnc'); % 標準脳ファイル名
VG=strcat(pathname,filename);
[filename,pathname]=uigetfile('*.img'); % 被験者脳ファイル名
VF=strcat(pathname,filename);
matname=strcat(pathname,'abc.mat'); % 保存ファイル名
prm=spm_normalise(VF,VG,matname,'','',defaults.normalise.estimate);
msk=spm_write_sn(VF,prm);

spm_write_vol.m

書式 VO=spm_write_vol(VI,dat) VI.fnameと同じ名前+拡張子.matの変換行列も保存される。
三次元配列の輝度変化に注意
dim(1)はMRIcroのX size×-1
警告 Cant get default Analyze orientation - assuming flippedとでるが無視します。
fwriteを使って保存し、MRIcroでheaderを書き足すほうがわかりやすい

任意の三次元配列をAnalyse形式にする
[filename,pathname]=uiputfile('*.img');
V.fname=strcat(pathname,filename);
V.dat=int16(rand(256,256,86)*2^8); % 256x256x86の2byte整数と仮定
V.dim=[size(V.dat),4]; % 4は2byte整数の意味
V.pinfo=[1,0,0]'; % 各スライスの信号強度は三次元配列の値をそのまま使う
V.mat=diag([-0.9375,0.9375,2,1]); % voxelの大きさは0.9375x0.9375x2mm
VO=spm_write_vol(V,V.dat); % 警告が出るが無視します

例 *.imgを*_new.imgという名前で保存(fwrite)
[filename,pathname]=uigetfile('*.img');
V=spm_vol(strcat(pathname,filename));
V.dat=spm_read_vols(V);
fid=fopen(V.fname,'wb');
fwrite(fid,V.dat,'int16');
fclose(fid); % 後でMRIcroのheaderを書き足すこと!