平面型グラジオメータの等磁場線図

平面型グラジオメータがある場合、そのままでは等磁場線図を書くことができません。 Neuromagは一旦7cmの球面上の格子点を設定して最小ノルム法で電流解を求め、その電流解がDewar面上のマグネトメータ上で鎖交する磁場から等磁場線図を作成しています。

仮想球の設定
仮想球の作成
導出磁場行列の計算
磁場計算行列の逆行列＝空間フィルタ（最小L2ノルム）の作成
最小L2ノルムの検算
平面型グラジオメータのデータを等磁場線図に

仮想球の設定

% functionからreturnまでの%を削除してreturnまでをコピーしてhns_MakeSphere.mという名前で保存してください。

% function R=hns_MakeSphere(centroid,radius,seg)
% if nargin==0;centroid=[0,0,0.04];end;
% if nargin<2;radius=0.05:0.01:0.07;end;
% if nargin<3;seg=15;end;
% ang=[];span=pi/seg;
% Phi=[-0.5,1]*pi/2;
% phi=linspace(Phi(1),Phi(2),fix((Phi(2)-Phi(1))/span));
% phi(end)=[];
%
% for n=1:length(phi);
%     N=round(cos(phi(n))*2*pi/span);
%     theta=linspace(0,pi*2,N+1)+n*pi/seg*2;
%     ang=[ang,[theta(1:end-1);phi(n)*ones(1,N)]];
% end;
%
% theta=ang(1,:);phi=ang(2,:);
% P=[cos(theta).*cos(phi);sin(theta).*cos(phi);sin(phi)];
% V1=[cos(theta).*sin(phi);sin(theta).*sin(phi);-cos(phi)];
% V2=[P(2,:).*V1(3,:)-P(3,:).*V1(2,:);...
%     P(3,:).*V1(1,:)-P(1,:).*V1(3,:);...
%     P(1,:).*V1(2,:)-P(2,:).*V1(1,:)];
%
% p=[];v1=[];v2=[];
% for n=1:length(radius);% 重殻構造用
%     theta=ang(1,:)+n/2*pi;
%     phi=ang(2,:);
%     P=[cos(theta).*cos(phi);sin(theta).*cos(phi);sin(phi)];
%     V1=[cos(theta).*sin(phi);sin(theta).*sin(phi);-cos(phi)];
%     V2=[P(2,:).*V1(3,:)-P(3,:).*V1(2,:);...
%         P(3,:).*V1(1,:)-P(1,:).*V1(3,:);...
%         P(1,:).*V1(2,:)-P(2,:).*V1(1,:)];
%     頭頂部
%     P(:,end+1)=[0,0,1]';
%     V1(:,end+1)=[1,0,0]';
%     V2(:,end+1)=[0,1,0]';
%
%     p =[p,P*radius(n)];
%     v1=[v1,V1];
%     v2=[v2,V2];
% end;
%
% P=p;V1=v1;V2=v2;
% P=P+centroid'*ones(1,size(P,2));
% Q1=[P;V1];% 1Aとする
% Q2=[P;V2];
% R.P=P';
% R.Q1=Q1';
% R.Q2=Q2';
% return;

仮想球の作成

仮想球を作成します。球の中心座標を[0,0,0.04]m、半径を0.07m、緯度方向の分割数を15とします。

centroid=[0,0,0.04];
radius=0.07;
seg=40;
Sp=hns_MakeSphere(centroid,radius,seg);
P=Sp.P;
Q1=Sp.Q1;
Q2=Sp.Q2;
figure('color',[1,1,1]);
plot3(P(:,1),P(:,2),P(:,3),'ro');hold on;
h1=quiver3(Q1(:,1),Q1(:,2),Q1(:,3),Q1(:,4),Q1(:,5),Q1(:,6));
h2=quiver3(Q2(:,1),Q2(:,2),Q2(:,3),Q2(:,4),Q2(:,5),Q2(:,6));
set([h1(:);h2(:)],'color',[0,0,1]);view([0,90]);
axis tight;daspect([1,1,1]);rotate3d on;

導出磁場行列の計算

仮想球上の格子点とセンサーの位置関係から導出磁場行列を求めます。

centroid=[0,0,0.04];
Q=[Q1;Q2];
nq=size(Q,1);
% Dewar
R=hns_GetCoils(0);
nsns=size(R.sensor,1);% radial magnetometer
L=zeros(nq,nsns);
for n=1:nsns
    L(:,n)=hns_Sarvas(Q,R.sensor(n,:),centroid);
end;
Ldwmag=L';% センサー数 × 格子点数x2　に

% Neuromag
R=hns_GetCoils(1);
nsns=size(R.sensorpla1,1);% radial planar gradiometer
L=zeros(nq,nsns);
for n=1:nsns
    L(:,n)=hns_Sarvas(Q,R.sensorpla1(n,:),centroid)...
          -hns_Sarvas(Q,R.sensorpla2(n,:),centroid);
end;
Lvvpla=1/1.68*L';% センサー数 × 格子点数x2　に
nsns=size(R.sensormag1,1);% radial magnetometer
L=zeros(nq,nsns);
for n=1:nsns
    L(:,n)=hns_Sarvas(Q,R.sensormag1(n,:),centroid)...
          +hns_Sarvas(Q,R.sensormag2(n,:),centroid)...
          +hns_Sarvas(Q,R.sensormag3(n,:),centroid)...
          +hns_Sarvas(Q,R.sensormag4(n,:),centroid);
end;
Lvvmag=1/4*L';% センサー数 × 格子点数x2　に

% PQ1400R
R=hns_GetCoils(2);
nsns=size(R.sensorgra1,1);% radial axial gradiometer
L=zeros(nq,nsns);
for n=1:nsns
    L(:,n)=hns_Sarvas(Q,R.sensorgra1(n,:),centroid)...
          -hns_Sarvas(Q,R.sensorgra2(n,:),centroid);
end;
Lpqgra=L';% センサー数 × 格子点数x2　に
nsns=size(R.sensorpla1,1);% tangential planar gradiometer
L=zeros(nq,nsns);
for n=1:nsns
    L(:,n)=hns_Sarvas(Q,R.sensorpla1(n,:),centroid)...
          -hns_Sarvas(Q,R.sensorpla2(n,:),centroid);
end;
Lpqpla=L';% センサー数 × 格子点数x2　に

磁場計算行列の逆行列＝空間フィルタ（最小L2ノルム）の作成

格子点上に電流源推定をする最小ノルム法を持ちいた空間フィルタを作成します。逆行列の近似はいろいろあるのですが、できるだけランクを落としたくないのと、 Matlabの計算精度に任せた方がいいと思うので安直にpinv関数を使うことにします。

MNdwmag=pinv(Ldwmag);
MNvvpla=pinv(Lvvpla);
MNvvmag=pinv(Lvvmag);
MNpqgra=pinv(Lpqgra);
MNpqpla=pinv(Lpqpla);

最小L2ノルムの検算

[0.07,0,0.04]mのところに[0,1,1]nAの電流が流れた時に各脳磁計のコイルで検出された磁場から上で作成した空間フィルタを使って電流源推定します。結果は既知電流の両側に渦巻くよう電流源が推定されます。電流の総和が最小となるというフィルタなので、タイトルに電流量の総和[nA]を書いてみました。 1nAが数倍に大きくなって推定されてしまいました。

centroid=[0,0,0.04];
q=[[0.07,0,0.04],[0,1,1]*1e-9];
R=hns_GetCoils(0);% Dewar
B=hns_Sarvas2(q,R.sensor,centroid);
qs1=MNdwmag*B;
R=hns_GetCoils(1);% Neuromag System
Bpla=hns_Sarvas2(q,R.sensorpla1,centroid)...
    -hns_Sarvas2(q,R.sensorpla2,centroid);
Bpla=Bpla/1.68;
qs2=MNvvpla*Bpla;
Bmag=hns_Sarvas2(q,R.sensormag1,centroid)...
    +hns_Sarvas2(q,R.sensormag2,centroid)...;
    +hns_Sarvas2(q,R.sensormag3,centroid)...
    +hns_Sarvas2(q,R.sensormag4,centroid);
Bmag=Bmag/4;
qs3=MNvvmag*Bmag;
R=hns_GetCoils(2);% PQ1400R
Bgra=hns_Sarvas2(q,R.sensorgra1,centroid)...
    -hns_Sarvas2(q,R.sensorgra2,centroid);
qs5=MNpqgra*Bgra;
Bpla=hns_Sarvas2(q,R.sensorpla1,centroid)...
    -hns_Sarvas2(q,R.sensorpla2,centroid);
qs6=MNpqpla*Bpla;
V1=Sp.Q1(:,4:6)';
V2=Sp.Q2(:,4:6)';
for n=1:5;
    switch n
        case 1;qs=qs1;t='Dewar: Magnetometer';
        case 2;qs=qs2;t='Neuromag System: Planar Gradiometer';
        case 3;qs=qs3;t='Neuromag System: Magnetometer';
        case 4;qs=qs5;t='PQ1400R: Axial Gradiometer';
        case 5;qs=qs6;t='PQ1400R: Planar Gradiometer';
    end;
    figure('color',[1,1,1]);
    nq=size(qs,1)/2;
    q1=qs(1:nq)*ones(1,3);
    q2=qs(nq+(1:nq))*ones(1,3);
    q1=q1*1e-9;
    q2=q2*1e-9;
    Q=[P,q1.*V1'+q2.*V2'];
    plot3(Q(:,1),Q(:,2),Q(:,3),'.','color',[1,1,1]*0.8);hold on;
    qq(4:6)=q(4:6)*3e+7;% 見えやすいように矢印の大きさ補正
    h=quiver3(q(1),q(2),q(3),qq(4),qq(5),qq(6),'color',[0,0,1]);
    set(h,'linewidth',2);
    h=quiver3(Q(:,1),Q(:,2),Q(:,3),Q(:,4),Q(:,5),Q(:,6),'color',[1,0,0]);
    set(h,'linewidth',2);
    daspect([1,1,1]);axis tight;view([90,0]);
    rotate3d on;
    st=[t,sprintf(' %0.2f nA',sum(abs(qs))*1e+9)];
    title(st);
end;

平面型グラジオメータのデータを等磁場線図に

最小L2ノルムの空間フィルタを使って各脳磁計のデータをDewar表面上のマグネトメータに変換する行列を作成します。 7cmの球面上の電流→各脳磁計の各コイルのデータ→変換行列→Dewar上のデータに変換して等磁場線図を書いてみます。結果は・・・頭頂部が少し変ですが、まあ良しとします。

Tvvpla2dwmag=Ldwmag*MNvvpla;
Tvvmag2dwmag=Ldwmag*MNvvmag;
Tpqgra2dwmag=Ldwmag*MNpqgra;
Tpqpla2dwmag=Ldwmag*MNpqpla;
q=[[0.07,0,0.04],[0,1,1]*1e-9];
centroid=[0,0,0.04];
R=hns_GetCoils(0);
Bdwmag=hns_Sarvas2(q,R.sensor,centroid);
seg=100;
px=R.ChPos(:,1);
py=R.ChPos(:,2);
xx=linspace(min(px),max(px),seg);
yy=linspace(min(py),max(py),seg);
[xx,yy]=meshgrid(xx,yy);
LevelStep=10;

R=hns_GetCoils(1);
Bvvpla=1/1.68*(hns_Sarvas2(q,R.sensorpla1,centroid)...
              -hns_Sarvas2(q,R.sensorpla2,centroid));
Bvvmag=1/4*(hns_Sarvas2(q,R.sensormag1,centroid)...
           +hns_Sarvas2(q,R.sensormag2,centroid)...
           +hns_Sarvas2(q,R.sensormag3,centroid)...
           +hns_Sarvas2(q,R.sensormag4,centroid));
R=hns_GetCoils(2);
Bpqgra=hns_Sarvas2(q,R.sensorgra1,centroid)...
      -hns_Sarvas2(q,R.sensorgra2,centroid);
Bpqpla=hns_Sarvas2(q,R.sensorpla1,centroid)...
      -hns_Sarvas2(q,R.sensorpla2,centroid);

for n=1:5
    switch n;
        case 1;B=Bdwmag;t='Dewar';
        case 2;B=Tvvpla2dwmag*Bvvpla;t='Neuromag System Planar Gradioemter';
        case 3;B=Tvvmag2dwmag*Bvvmag;t='Neuromag System Magnetometer';
        case 4;B=Tpqgra2dwmag*Bpqgra;t='PQ1400R Axial Gradiometer';
        case 5;B=Tpqpla2dwmag*Bpqpla;t='PQ1400R Planar Gradiometer';
    end;
    figure('color',[1,1,1]);
    [x,y,z]=griddata(px,py,B,xx,yy);
    LevelList=fix(min(B)/LevelStep)*LevelStep:LevelStep:-LevelStep;
    plot(px,py,'.','color',[1,1,1]*0.8);hold on;
    contour(x,y,z,'LineColor',[0,0,1],'LevelList',LevelList);
    contour(x,y,z,'LineColor',[1,1,1]*0.1,'LevelList',0);
    LevelList=LevelStep:LevelStep:fix(max(B)/LevelStep)*LevelStep;
    contour(x,y,z,'LineColor',[1,0,0],'LevelList',LevelList);
    t=[t,sprintf(' %0.1f ~ %0.1f',min(B),max(B))];
    axis tight;axis off;title(t);
end;