Neuromagのphantomでいろいろ計算

Neuromagのphantomを使ってGOFを計算してみます。

Contents

Sarvas式の作成

functionからreturnまでの%を削除してreturnまでをコピーしてhns_Sarvas2.mという名前で保存してください。 qは[x,y,z,qx,qy,qz]の電流双極子です。 SはNs個×[x,y,z,zx,zy,zz]のセンサー位置と単位ベクトル化した法線方向です。 pは[x,y,z]にある仮想球の中心座標です。

% function Bz=hns_Sarvas2(q,S,p)
% % Sarvas' rule
% % q = 1x [x,y,z,qx,qy,qz]      current dipole
% % S = Ns x [x,y,z,zx,zy,zz] sensor position & direction
% % p = [x,y,z]               sphere centroid
% % written by Akia Hashzume
% Ns=size(S,1);
% R=S(:,1:3)-ones(Ns,1)*p;
% R0=ones(Ns,1)*(q(1:3)-p);
% A=R-R0;
%
% r =sqrt(sum(R.^2,2));
% a =sqrt(sum(A.^2,2));
% ar=sum(A.*R,2);
% F=(a.*(r.*a+ar))*ones(1,3);
% n=a+2*r+ar./a;
% nn=(a.^2)./r+a+n;
% n=n*ones(1,3);
% nn=nn*ones(1,3);
% nablaF=nn.*R-n.*R0;
% Q=ones(Ns,1)*q(4:6);
%
% QR0=[Q(:,2).*R0(:,3)-Q(:,3).*R0(:,2),...
%      Q(:,3).*R0(:,1)-Q(:,1).*R0(:,3),...
%      Q(:,1).*R0(:,2)-Q(:,2).*R0(:,1)];
% QR0R=sum(QR0.*R,2)*ones(1,3);
% B=QR0-(QR0R.*nablaF)./F;
% B=B./F*(1e-7)*(1e+15);% fT
% Bz=S(:,4).*B(:,1)+S(:,5).*B(:,2)+S(:,6).*B(:,3);
% return;

データの読み込み

ファイルを http://meg.aalip.jp/matlab/data/dipole1.mat から入手し,作業ディレクトリに保存してください。 F.dipoleは19.9795msの時、頭座標上の仮想球中心[0,0,0]とし、 Xfitで単一電流双極子で推定した結果です。 F.dipole(:,1)は全センサーを、 F.dipole(:,2)は平面型グラジオメーターのみを、 F.dipole(:,3)はマグネトメーターのみを使って推定した結果です。 F.vv_gra,F.vv_mag,F.vv_loutは波形やトポグラフィを描くときのセンサーの二次元座標です。 graは平面型グラジオメーターのみのとき、magはマグネトメーターのみのとき、loutは全306chのときです。

load dipole1.mat;
meg=F.ave_0001(:,1:306);
time=F.ave_0001(:,end);
SSP=eye(306);
for n=1:10 % SSPの計算
    EigenVector=eval(sprintf('F.PCA_%04.0f',n));
    SSP=SSP-EigenVector*EigenVector';
end;
megSSP=meg*SSP;
mag=3:3:306;
gra=setdiff(1:306,mag);
figure('color',[1,1,1]);
subplot(221);plot(time,meg(:,gra),'b');axis tight;title('gradiometer');
subplot(222);plot(time,meg(:,mag),'g');axis tight;title('magnetometer');
subplot(223);plot(time,megSSP(:,gra),'b');axis tight;title('gradiometer + SSP');
subplot(224);plot(time,megSSP(:,mag),'g');axis tight;title('magnetometer + SSP');

センサー座標と法線方向の計算

% まずセンサー座標から頭座標へ変換
MTX=F.MEG2HEAD(:,1:3);
MT0=F.MEG2HEAD(:,4);
Ch=MTX*F.ChInform(1:3,1:306)+MT0*ones(1,306);
X=MTX*F.ChInform(4:6,1:306);
Y=MTX*F.ChInform(7:9,1:306);
Z=MTX*F.ChInform(10:12,1:306);

% gradiometer
Sp1=[Ch(:,gra)+0.0084*X(:,gra)+0.0003*Z(:,gra);Z(:,gra)]';
Sp2=[Ch(:,gra)-0.0084*X(:,gra)+0.0003*Z(:,gra);Z(:,gra)]';
% magnetometer
Sm1=[Ch(:,mag)+0.0129*X(:,mag)+0.0129*Y(:,mag)+0.0003*Z(:,mag);Z(:,mag)]';
Sm2=[Ch(:,mag)+0.0129*X(:,mag)-0.0129*Y(:,mag)+0.0003*Z(:,mag);Z(:,mag)]';
Sm3=[Ch(:,mag)-0.0129*X(:,mag)+0.0129*Y(:,mag)+0.0003*Z(:,mag);Z(:,mag)]';
Sm4=[Ch(:,mag)-0.0129*X(:,mag)-0.0129*Y(:,mag)+0.0003*Z(:,mag);Z(:,mag)]';

磁場の計算

青が実測値、緑が計算値です。 magnetometerのGOFはよくありません。 SSP使わないほうが、GOFはいいみたいです・・・

q=[F.dipole(3:5,2);F.dipole(7:9,2)*1e-9]';p=[0,0,0];
Bp=1/1.68*(hns_Sarvas2(q,Sp1,p)-hns_Sarvas2(q,Sp2,p));
q=[F.dipole(3:5,3);F.dipole(7:9,3)*1e-9]';p=[0,0,0];
Bm=1/4*(hns_Sarvas2(q,Sm1,p)+hns_Sarvas2(q,Sm2,p)+hns_Sarvas2(q,Sm3,p)+hns_Sarvas2(q,Sm4,p));
[dummy,t]=min(abs(time-F.dipole(1)));
grach=1:204;magch=1:102;
figure('color',[1,1,1]);

BP=meg(t,gra)';    gof=1-(BP-Bp)'*(BP-Bp)/(BP'*BP);
BPs=megSSP(t,gra)';gofs=1-(BPs-Bp)'*(BPs-Bp)/(BPs'*BPs);

subplot(221);
plot(grach,[BP,Bp]');axis tight;
title(sprintf('GRA GOF %0.2f',gof*100));
subplot(222);
plot(grach,[BPs,Bp]');axis tight;
title(sprintf('GRA+SSP GOF %0.2f',gofs*100));

BM=meg(t,mag)';    gof=1-(BM-Bm)'*(BM-Bm)/(BM'*BM);
BMs=megSSP(t,mag)';gofs=1-(BMs-Bm)'*(BMs-Bm)/(BMs'*BMs);

subplot(223);
plot(magch,[BM,Bm]');axis tight;
title(sprintf('MAG GOF %0.2f',gof*100));
subplot(224);
plot(magch,[BMs,Bm]');axis tight;
title(sprintf('MAG+SSP GOF %0.2f',gofs*100));

topography化

P=F.vv_mag(:,2:3);seg=100;
[X,Y]=meshgrid(linspace(min(P(:,1)),max(P(:,1)),seg),linspace(min(P(:,2)),max(P(:,2)),seg));

figure('color',[1,1,1]);
for n=1:6
    subplot(2,3,n);
    switch n
        case 1;B=sqrt(BP(1:2:end).^2+(BP(2:2:end).^2));k=0;
        case 2;B=sqrt(BPs(1:2:end).^2+(BPs(2:2:end).^2));k=0;
        case 3;B=sqrt(Bp(1:2:end).^2+(Bp(2:2:end).^2));k=0;
        case 4;B=BM;k=-1;
        case 5;B=BMs;k=-1;
        case 6;B=Bm;k=-1;
    end;
    Z=griddata(P(:,1),P(:,2),B,X,Y);
    surf(Z);set(gca,'clim',[k,1]*max(abs(B)));
    axis tight;axis off;view([0,90]);
    colorbar;
    switch n
        case 1;title('GRA measured');
        case 2;title('GRA+SSP measured');
        case 3;title('GRA simulated');
        case 4;title('MAG measured');
        case 5;title('MAG+SSP meadured');
        case 6;title('MAG simulated');
    end;
end;
set(findobj('type','surface'),'linestyle','none');

Published with MATLAB® 7.14