sn.matの内容

標準脳化後に作成された*sn.matをみてみます。 因みに*sn.matの内容はspm_normalise.mの返り値です。
ところどころ間違っているかもしれません。

[filename,pathname]=uigetfile('*sn.mat'); % *sn.matを選択
load(strcat(pathname,filename)); % *sn.matを開く
help spm_normalise
help spm_affreg
help spm_vol
 などから得られる説明を参考にしました。
*sn.matには、 VGは選択した標準脳のheader情報、 VFは被験者の脳のheader情報、 flagsは各種parameterの値、 AffineはVGからVFへの線形変換行列、 Trはたぶん非線形変換に関する行列が含まれています。

VFとVGは構造体になっています。
fnameファイル名
dimx,y,zの三次元配列の要素数とデータ型
mat三次元配列を実空間座標に変換する行列
pinfo(1,:)(:,:,z)平面上の信号のscale
pinfo(2,:)(:,:,z)平面上の信号のoffset
pinfo(3,:)(:,:,z)平面へのoffset?
n?
descrip撮影時間やファイル形式など
privatefnameのheader情報（構造体の入れ子）
dat実際の三次元配列（削除されている）
pinfoは(:,:,z)平面の信号を
new_val=rand(256,256)*50; % 256×256からなる平面上の信号強度と仮定
z=50; % 50番目のスライスと仮定
original_val=new_val*pinfo(1,z)+pinfo(2,z)
 に変換したことを示すそうです。

datは削除されています。読み込むにはspm_read_volsを使います。
[VF.dat,xyz]=spm_read_vols(VF);
figure;imagesc(rot90(squeeze(VF.dat(round(VF.dim(1)/2),:,:))),[0,384]);
daspect(abs([VF.mat(3,3),VF.mat(2,2),VF.mat(1,1)]));colormap(gray);
xyzはVF.datのx,y,z座標からなる行列です。
clear xyz; % to save memory
VG.dat=spm_read_vols(VG);
figure;imagesc(rot90(squeeze(VG.dat(round(VG.dim(1)/2),:,:))),[0,1]);% 最大強度は1
daspect(abs([VG.mat(3,3),VG.mat(2,2),VG.mat(1,1)]));colormap(gray);
flagも構造体になっています。
smosrc被験脳の平滑化の半値幅8
smoref標準脳の平滑化の半値幅0
regtype規則化の方式
'none','rigid','subj','mni'がある 		mni
weight重み付け（病変脳用？）''
cutoff離散cosine変換の切捨て?30
nits非線形変換の繰り返し回数16
reg規則化の変換量?0.1
wtsrc被験者脳の重み付けの何か？0
graphicsspm_normalise_disp.mの実行1
versionversion'@(#)spm_normalise.m 2.8 03/03/04'
dateSPM実行日Matlabのdate関数
spmから起動した場合には、初期値は
spm_defaults
 で定義されたものが使われます。

被験者脳の標準脳への変換行列は変数Affineが該当します。すなわち

ここで、

は被験者脳および標準脳の三次元配列を、

は標準脳の三次元配列を実空間の座標に変換する行列（VG.mat）、 被験者脳の実空間座標を標準脳の実空間座標に変換する行列（SPM実行時にLinear {affine} componentとして表示される行列）、 被験者脳の三次元配列を実空間の座標に変換する行列（VF.mat）です。
変数Affineは

です。逆行列になっています。
線形変換だけでよかったら、脳磁図の電流双極子の座標を標準脳の座標に変換できそうです。
M=VG.mat*inv(VF.mat*Affine);
dip=[70,0,0]/1000; % MRI座標上で(70,0,0)mmのところに双極子があると仮定
new_dip=M*[dip,1]';
new_dip(end)=[];
Trは非線形変換に関するもののようです。 たぶん、三次元離散cosine変換後のx,y,z方向のvoxelの数＋移動方向のベクトルからなる四次元配列になっています。 cuttoffは、spm_default.mでは25に、cutoffが未定の場合spm_normaliseでは自動的に30になっています。 spm_normalise.mをみてみると、三次元離散cosine変換後のx,y,z方向の配列の大きさは
tmp = sqrt(sum(VG(1).mat(1:3,1:3).^2));
k = max(round((VG(1).dim(1:3).*tmp)/cutoff),[1 1 1]);
 で定義されていました。 コードを見る限りはcutoffはvoxelの大きさを意味するように見えます。
図にしてみます。
nx=1:size(Tr,1); % x方向のvoxel数
ny=1:size(Tr,2); % y方向のvoxel数
nz=1:size(Tr,3); % z方向のvoxel数
vx=squeeze(Tr(:,:,:,1)); % 各voxelのx方向の移動量
vy=squeeze(Tr(:,:,:,2)); % 各voxelのy方向の移動量
vz=squeeze(Tr(:,:,:,3)); % 各voxelのz方向の移動量
[nx,ny,nz]=ndgrid(nx,ny,nz);
figure;set(gcf,'color',[1,1,1]);
quiver3(nx,ny,nz,vx,vy,vz,2) % scaleは2
set(findobj('type','line'),'linewidth',2); %線の太さを2に
hold on; % 上書き許可
plot3(nx(:),ny(:),nz(:),'b.'); % 起点を青点で表示
daspect([1,1,1]); % 縦横高さを1:1:1にする
view([-80,40]); % 視点を変える
原点(0,0,0)近くで移動量が大きいことがわかります。 MRIのk-spaceみたいに原点近くに大部分の信号が集まっているのだと思います。 Matlabのsignal processing toolboxの関数idct使って逆三次元離散コサイン変換を試みてみます。
for x=1:7;for y=1:9;...
vx(x,y,:)=idct(squeeze(vx(x,y,:)));...
vy(x,y,:)=idct(squeeze(vy(x,y,:)));...
vz(x,y,:)=idct(squeeze(vz(x,y,:)));...
end;end;
for y=1:9;for z=1:7;...
vx(:,y,z)=idct(squeeze(vx(:,y,z)));...
vy(:,y,z)=idct(squeeze(vy(:,y,z)));...
vz(:,y,z)=idct(squeeze(vz(:,y,z)));...
end;end;
for x=1:7;for z=1:7;...
vx(x,:,z)=idct(squeeze(vx(x,:,z)));...
vy(x,:,z)=idct(squeeze(vy(x,:,z)));...
vz(x,:,z)=idct(squeeze(vz(x,:,z)));...
end;end;
figure;set(gcf,'color',[1,1,1]);
quiver3(nx,ny,nz,vx,vy,vz,2) % scaleは2
set(findobj('type','line'),'linewidth',2); %線の太さを2に
hold on; % 上書き許可
plot3(nx(:),ny(:),nz(:),'b.'); % 起点を青点で表示
daspect([1,1,1]); % 縦横高さを1:1:1にする
view([-80,40]); % 視点を変える
たぶん25mmボクセル毎の移動量を意味しているものと思われます。 標準脳から被験者脳への移動量か、被験者脳から標準脳への移動量かはよくわかりません。