非公式広島大学病院脳磁図室ホームページMRIcro&SPM篇

Affine変換

spm_affreg.mを2回実行して得られた返り値が、Affine変換の変換行列となっています。引数はspm_smoothto8bit.mで平滑処理後に信号強度を微調整をされたものです。
変数を追ってみます。 Affine変換の最適化の条件、変数aflagはspm_normalise.m内の値を使うこととし、重み付けはなしとします。
[filename,pathname]=uigetfile('*sn.mat');
load(strcat(pathname,filename));
VF1=spm_smoothto8bit(VF,8);
VG1=VG; % 既に平滑化済み
VF1.pinfo(1:2,:)=VF1.pinfo(1:2,:)/spm_global(VF1); %平均値/8以下のvoxelを無視して補正
VG1.pinfo(1:2,:)=VG1.pinfo(1:2,:)/spm_global(VG1);
aflags=struct('sep',max(flags.smoref,flags.smosrc), 'regtype',flags.regtype,...
'WG',[],'WF',[],'globnorm',0);
M0=eye(4);
[M1,scal1]=spm_affreg(VG1,VF1,aflags,M0); % 5mm以内を近似?
aflags.sep=4; %
[M2,scal2]=spm_affreg(VG1,VF1,aflags,M1); % 8/2mm以内を近似?
原点(0,0,0)にある一辺200mmの立方体がどのように変形するか見てみます。
[x,y,z]=ndgrid(-1:2:1,-1:2:1,-1:2:1);
V=[x(:),y(:),z(:)]'*10; % 一辺20の立方体の座標
W=zeros(3,4,6);
W(:,:,1)=V(:,[1,2,4,3]);
W(:,:,2)=V(:,[3,4,8,7]);
W(:,:,3)=V(:,[7,8,6,5]);
W(:,:,4)=V(:,[5,6,2,1]);
W(:,:,5)=V(:,[5,1,3,7]);
W(:,:,6)=V(:,[2,6,8,4]);
figure;set(gcf,'color',[1,1,1],'renderer','zbuffer');
W=[reshape(W,3,4*6);ones(1,4*6)];
W0=M0*W; % 0回目のAffine変換
W0=reshape(W0(1:3,:),3,4,6);
h0=fill3(squeeze(W0(1,:,:)),squeeze(W0(2,:,:)),squeeze(W0(3,:,:)),zeros(4,6));
set(h0,'FaceColor','none','EdgeColor',[0,0,1],'LineWidth',2);hold on;
W1=M1*W; % 1回目のAffine変換
W1=reshape(W1(1:3,:),3,4,6);
h1=fill3(squeeze(W1(1,:,:)),squeeze(W1(2,:,:)),squeeze(W1(3,:,:)),zeros(4,6));
set(h1,'FaceColor','none','EdgeColor',[0,0.7,0],'LineWidth',2);
W2=M2*W; % 2回目のAffine変換
W2=reshape(W2(1:3,:),3,4,6);
h2=fill3(squeeze(W2(1,:,:)),squeeze(W2(2,:,:)),squeeze(W2(3,:,:)),zeros(4,6));
set(h2,'FaceColor','none','EdgeColor',[1,0,0],'LineWidth',2);
daspect([1,1,1]);axis tight;

青→緑→赤の順でAffine変換が最適化されたことを示します。青→緑がおおまかに位置合わせを行う'coarse Affine registration'で、緑→赤が微調整を行う'fine Affine registration'です。緑→赤はほとんど変化していません。

Affine変換後の画像を見るにはspm_slice_vol.dllを使います。spm_slice_vol.dllは
変換後のz=0のx-y平面画像=spm_slice_vol(三次元配列,Affine変換の逆行列,[x,y],補間形式)
なる構文を取ります。x,yは[0,0]から[x,y]の平面という意味です。
VV=spm_read_vols(VF);;
K=eye(4);K(1:3,4)=[120,120,0]';
W=spm_slice_vol(VV,inv(K*VF.mat),[240,240],1);
figure;imagesc(rot90(W),[0,384]);colormap(gray);daspect([1,1,1]);

冠状断にします。
CR=eye(4);CR(2:3,2:3)=rot90(CR(2:3,2:3));
W=spm_slice_vol(VV,inv(K*CR*VF.mat),[240,240],1);
imagesc(rot90(W),[0,384]);daspect([1,1,1]);

Affine変換後の水平断・冠状断・矢状断を表示します。
SG=eye(4);SG(1:3,1:3)=rot90(SG(1:3,1:3));
M2=VG.mat*inv(VF.mat*Affine);
W=zeros(480,480);
W(1:240,1:240)=rot90(spm_slice_vol(VV,inv(K*M2*VF.mat),[240,240],1));
W(1:240,241:end)=rot90(spm_slice_vol(VV,inv(K*CR*M2*VF.mat),[240,240],1));
W(241:end,1:240)=rot90(rot90(spm_slice_vol(VV,inv(K*SG*M2*VF.mat),[240,240],1)));
imagesc(W,[0,384]);daspect([1,1,1]);

Analyse formatにして、MRIcroで脳表を見てみます。
VVV=zeros(240,240,86*2+1);
KK=eye(4);
for z=-86:86;...
K(3,4)=z;...
VVV(:,:,z+87)=spm_slice_vol(VV,inv(KK*K*M2*VF.mat),[240,240],1);...
end;
VVV=int16(VVV);
[filename,pathname]=uiputfile('*.img');
fid=fopen(strcat(pathname,filename,'.img'),'wb');
fwrite(fid,VVV,'int16');
fclose(fid);
MRIcroを立ち上げ、Headerを下記の条件とし、*.hdrの形で保存します。

次に、*.imgを読み込み、skull stripボタンを押して脳表画像を作成しました。

・・・Affine変換だけだと標準脳ぽくないです。