Dewar上の磁場分布2
Scilabを使って、Neuromag社のSource Modellingの等磁場線図もどきを作りたいと思います。
センサー信号→球表面メッシュ点の電流源→Dewarメッシュ点の磁場
という順序で計算していきます。
球表面やDewarのメッシュ点情報は
Dewar.mat にあります。
まずセンサー情報をdevice座標からhead座標に座標変換します。
stacksize(50000000);
loadmatfile('.../Dewar.mat');
loadmatfile('.../xxx.mat');
getf('.../Sarvas.sce');
getf('.../CalcSSP2.sce');
// conversion of coordination system from device to head
MTX=MEG2HEAD(1:3,1:3);
Ch.P=MTX*ChInform(1:3,1:306)+MEG2HEAD(1:3,4)*ones(1,306);
Ch.X=MTX*ChInform(4:6,1:306);
Ch.Y=MTX*ChInform(7:9,1:306);
Ch.Z=MTX*ChInform(10:12,1:306);
次に、球のメッシュ点、法線方向の座標変換後、
接平面方向の単位ベクトルを計算します。
球面上の点を
とおくと、(x,y,z)における接平面のΦ方向のベクトルは
となります。
// spherical nodes definition
radius=70/1000;
Center=[0,0,40]/1000;
Sp.P=MTX*Sphere.nodes*radius+Center'*ones(1,size(Sphere.nodes,2));
Sp.V=MTX*Sphere.nodes;
// calculation of tangential vectors at sphere nodes
SinPhi=Sp.V(3,:);
N=find(SinPhi>-0.6);
// Bottom nodes on sphere are not shown in Source Modelling
//phi=asin(Sp.V(3,:));
CosPhi=(1-SinPhi.^2).^0.5;
check=find(CosPhi==0);
CosPhi(check)=1;
CosTheta=Sp.V(1,:)./CosPhi;
SinTheta=Sp.V(2,:)./CosPhi;
Sp.V1=[-CosTheta.*SinPhi;-SinTheta.*SinPhi;CosPhi];
Sp.V1(:,check)=[1,0,0]'*ones(1,length(check));
Sp.V1=Sp.V1./(ones(3,1)*sum(Sp.V1.^2,1).^0.5);
Sp.V2=[Sp.V(2,:).*Sp.V1(3,:)-Sp.V(3,:).*Sp.V1(2,:);...
Sp.V(3,:).*Sp.V1(1,:)-Sp.V(1,:).*Sp.V1(3,:);...
Sp.V(1,:).*Sp.V1(2,:)-Sp.V(2,:).*Sp.V1(1,:)];
Sp.V2=Sp.V2./(ones(3,1)*sum(Sp.V2.^2,1).^0.5);
Sp.V1=Sp.V1(:,N);
Sp.V2=Sp.V2(:,N);
Sp.V =Sp.V(:,N);
Sp.P =Sp.P(:,N);
球の接平面方向のベクトルとセンサー面の法線方向のベクトルを表示します。
scf();set(gcf(),'visible','off');
param3d(Sp.P(1,:),Sp.P(2,:),Sp.P(3,:));
h=gca();h.box='off';xset('use color',0);
h.children.mark_mode='on';h.children.mark_size=3;
P=[Sp.P;Sp.P+0.01*Sp.V1];
for x=1:size(P,2);param3d(P([1,4],x),P([2,5],x),P([3,6],x));end;
P=[Sp.P;Sp.P+0.01*Sp.V2];
for x=1:size(P,2);param3d(P([1,4],x),P([2,5],x),P([3,6],x));end;
param3d(Ch.P(1,:),Ch.P(2,:),Ch.P(3,:));
h.children(1).mark_mode='on';h.children(1).mark_size=3;
P=[Ch.P;Ch.P+0.01*Ch.Z];
for x=1:size(P,2);param3d(P([1,4],x),P([2,5],x),P([3,6],x));end;
h.box='off';h.axes_visible='off';h.isoview='on';
set(gcf(),'visible','on');
次に接平面ベクトル方向に1Aが流れたとして磁場導出行列を求めます。
因みに各チャンネルにおける信号の強さは以下の式で求めます。
ここで
bkはk番目のチャンネルの磁場信号、
Nkはk番目のチャンネルの計算点総数、
wkpはk番目のチャンネルのp番目の計算点の重み付け、
B (rkp )はk番目のチャンネルのp番目の計算点の磁場ベクトル、
nkpはk番目のチャンネルのp番目の計算点の法線ベクトル、
です。具体的には下記の表を用います。
coil | Nk | position r/mm | wkp |
---|
magnetometer | 4 | ±12.9,±12.9,0.3 | 1/4 |
---|
gradiometer | 2 | ±8.4, 0,0.3 | ±1/1.68cm |
---|
N=length(N);
LFM=zeros(N*2,306);
Q1=[Sp.P;Sp.V1]';
Q2=[Sp.P;Sp.V2]';
C1=3:3:306;
C2=1:3:306;
C3=2:3:306;
P1=Ch.P(:,C2)+8.4/1000*Ch.X(:,C2)+0.3/1000*Ch.Z(:,C2);
P2=Ch.P(:,C2)-8.4/1000*Ch.X(:,C2)+0.3/1000*Ch.Z(:,C2);
P3=Ch.P(:,C3)+8.4/1000*Ch.X(:,C3)+0.3/1000*Ch.Z(:,C3);
P4=Ch.P(:,C3)-8.4/1000*Ch.X(:,C3)+0.3/1000*Ch.Z(:,C3);
P5=Ch.P(:,C1)+12.9/1000*Ch.X(:,C1)+12.9/1000*Ch.Y(:,C1)+0.3/1000*Ch.Z(:,C1);
P6=Ch.P(:,C1)+12.9/1000*Ch.X(:,C1)-12.9/1000*Ch.Y(:,C1)+0.3/1000*Ch.Z(:,C1);
P7=Ch.P(:,C1)-12.9/1000*Ch.X(:,C1)+12.9/1000*Ch.Y(:,C1)+0.3/1000*Ch.Z(:,C1);
P8=Ch.P(:,C1)-12.9/1000*Ch.X(:,C1)-12.9/1000*Ch.Y(:,C1)+0.3/1000*Ch.Z(:,C1);
P=[P1,P2,P3,P4,P5,P6,P7,P8];
Z=Ch.Z;
Z=[Z(:,C2),Z(:,C2),Z(:,C3),Z(:,C3),Z(:,C1),Z(:,C1),Z(:,C1),Z(:,C1)];
wd=1/1.68; // weight of gradiometers
wm=1; // weight ratio of magnetometer/gradiometer
for x=1:N;...
LQ1=sum(Sarvas(Q1(x,:),P',Center).*Z',2)';...
LQ2=sum(Sarvas(Q2(x,:),P',Center).*Z',2)';...
// gradiometers
L12=wd*(LQ1(:,1:102)-LQ1(:,103:204));...
L22=wd*(LQ2(:,1:102)-LQ2(:,103:204));...
L13=wd*(LQ1(:,205:306)-LQ1(:,307:408));...
L23=wd*(LQ2(:,205:306)-LQ2(:,307:408));...
// magnetometers
L11=wm*(LQ1(:,409:510)+LQ1(:,511:612)+LQ1(:,613:714)+LQ1(:,715:816))/4;...
L21=wm*(LQ2(:,409:510)+LQ2(:,511:612)+LQ2(:,613:714)+LQ2(:,715:816))/4;...
LFM(x,:) =[L11,L12,L13];...
LFM(x+N,:)=[L21,L22,L23];...
end;
LFM=LFM'*1e+15;
ここで以下の式を考えます。
bcoilは各チャンネルの磁場信号、
qsphereは球表面メッシュ点上の電流の大きさ、
LcoilはSarvas式に基づいて計算した、
球表面上の接平面方向の直交する2つの単位ベクトル方向に1A流れたとして計算した磁場導出行列です。
という制約を設け、Lagrangeの未定係数法を用いると、球面上の電流は
で求まります。解の安定性を図る?ため、以下のように特異値分解を行い、特異値の上位n個までを用いることとします。
Unは行列Uの1〜n列だけの行列、
Λnは行列Λの1〜n行だけの行列です。
これを元の式に代入して球面上の電流を求めると以下のようになります。
こうして解かないと行列Lcoilがfull rankにならず、
逆行列が求まりません。
nn=120;
[U,S,V]=svd(LFM);
Ssvd=diag(S);Ssvd((nn+1):$)=0;
Ssvd=[diag(Ssvd),zeros(size(S,1),size(S,2)-size(S,1))];
LL=Ssvd(1:nn,:)*V';
これで準備はできました。
次に脳磁図信号を球面メッシュ点上の電流で表現してみます。
meg=CalcSSP2(ave_0001(:,1:306))';
time=ave_0001(:,$);
t=min(find(time>22));
MEG=[wm*meg(3:3:306,t);meg(1:3:306,t);meg(2:3:306,t)];
Q=LL'*inv(LL*LL')*(U(:,1:nn)'*MEG);
q1=Q(1:N)';q2=Q((N+1):$)';
q=[q1.*Sp.V1(1,:)+q2.*Sp.V2(1,:);...
q1.*Sp.V1(2,:)+q2.*Sp.V2(2,:);...
q1.*Sp.V1(3,:)+q2.*Sp.V2(3,:)];
scf();set(gcf(),'visible','off');
param3d(Sp.P(1,:),Sp.P(2,:),Sp.P(3,:));
h=gca();h.box='off';h.axes_visible='off';
h.children.mark_mode='on';h.children.mark_size=3;
P=[Sp.P;Sp.P+q*3e+6];
for x=1:size(P,2);param3d(P([1,4],x),P([2,5],x),P([3,6],x));end;
h.children(:).foreground=5;h.isoview='on';
set(gcf(),'visible','on');
真上および真横から見た球面上の電流です。
最後に球面メッシュ点上の電流からDewar表面上の磁場を計算します。そこで
から、まず磁場導出行列Ldewarを求める必要があります。
// conversion of coordination system from device to head
VV.P=MEG2HEAD(1:3,1:3)*VectorView.nodes;
VV.P=VV.P+MEG2HEAD(1:3,4)*ones(1,size(VectorView.nodes,2));
// calculation of radial vectors at mesh points on the dewar
T=VectorView.triangle;
VP=zeros(VV.P);
for x=1:size(VV.P,2);
k=(T(1,:)-x).*(T(2,:)-x).*(T(3,:)-x);...
k(k~=0)=1;k=(1:length(k)).*k;...
k(k==0)=[];TT=T;TT(:,k)=[];...
for y=1:size(TT,2);...
if TT(2,y)==x;TT(:,y)=TT([2,3,1],y);...
elseif TT(3,y)==x;TT(:,y)=TT([3,1,2],y);...
end;end;...
V1=VV.P(:,TT(2,:))-VV.P(:,TT(1,:));...
V2=VV.P(:,TT(3,:))-VV.P(:,TT(1,:));...
VP(:,x)=sum([V1(2,:).*V2(3,:)-V1(3,:).*V2(2,:);...
V1(3,:).*V2(1,:)-V1(1,:).*V2(3,:);...
V1(1,:).*V2(2,:)-V1(2,:).*V2(1,:)],2);end;
VV.V=VP./(ones(3,1)*sum(VP.^2,1).^0.5);
ベクトルの外積で重み付けをして各メッシュ点における法線ベクトルを求めています。
Ld=zeros(N*2,size(VV.P,2));
for x=1:N;...
Ld(x,:)= sum(Sarvas(Q1(x,:),VV.P',Center).*VV.V',2)';...
Ld(x+N,:)=sum(Sarvas(Q2(x,:),VV.P',Center).*VV.V',2)';...
end;
磁場導出行列が求まりました。
Bd=Ld'*Q*1e+15;
scf();set(gcf(),'visible','off');
T=VectorView.triangle;
T=T($:-1:1,:);
ntri=size(T,2);
X=matrix(VV.P(1,T),3,ntri);
Y=matrix(VV.P(2,T),3,ntri);
Z=matrix(VV.P(3,T),3,ntri);
B=matrix(Bd(T),3,ntri);
scale=32;scale2=scale*2;
Bz=B/200*scale+scale;Bz(Bz<1)=1;Bz(Bz>scale2)=scale2;
set(gcf(),'color_map',jetcolormap(scale2));
plot3d(X,Y,Z);ha=gca();f=ha.children.data;
TL=tlist(['3d','x','y','z','color'],f.x,f.y,f.z,Bz);
ha.children.data=TL;
ha.children.color_mode=-1;
ha.children.color_flag=3;
ha.isoview='on';ha.box='off';ha.margins=[0,0,0,0];
ha.axes_visible='off';ha.rotation_angles=[0,-90];
set(gcf(),'visible','on');
ha.rotation_angles=[30,180];
ホントは等磁場線図にしたいのですが、曲面の等高線の描き方がわからないので、ここまでです。