L2ノルム3

脳波と異なり、脳磁図では大脳皮質内の錐体細胞層だけが信号源であるとされています。 そこで、いっそのこと大脳皮質だけに双極子があると仮定してL2ノルムを解いてみます。
格子点？の設定は、Neuromag付属のSegLabを使います。 詳細は 皮髄境界mesh作成 を参考にしてください。
FTPを使ってPCへ転送後、Fiff2MatFileのmeshボタンを押し、MATファイル化します。

これで準備OKです。
stacksize(63000000);
loadmatfile(xgetfile('*_mesh.mat'));
表示してみます。
scf();set(gcf(),'visible','off');
param3d(NodesPts(1,:),NodesPts(2,:),NodesPts(3,:));
h=gca();h.children.mark_mode='on';
h.isoview='on';h.box='off';
set(gcf(),'visible','on');


因みに格子点数は43154個でした。
では実際の格子点を皮髄境界に限定したL2ノルムを作成してみます。
loadmatfile('.../xxx.mat');//SEF
getf('.../CalcLFM.sce');
getf('.../CalcSSP2.sce');
Ch=ConversionCoordination(MEG2HEAD,ChInform);
Center=[0,0,40]/1000;
R=CalcTangentialVectors(NodesPts(1:3,:)',Center);
wm=1;// gradiometer : magnetometer

次にRの容量が大きいので分割し、ファイルに保存しながら磁場導出行列を計算していきます。
nq=size(NodesPts,2);RP=R.P;RV2=R.V2;RV3=R.V3;R.V1=[];
n=1:10000;R.P=RP(n,:);R.V2=RV2(n,:);R.V3=RV3(n,:);
L1=real(CalcVectorViewLFM(R,Ch,Center,wm));
save('L1.dat',L1);clear('L1');
n=n+10000;R.P=RP(n,:);R.V2=RV2(n,:);R.V3=RV3(n,:);
L2=real(CalcVectorViewLFM(R,Ch,Center,wm));
save('L2.dat',L2);clear('L2');
n=n+10000;R.P=RP(n,:);R.V2=RV2(n,:);R.V3=RV3(n,:);
L3=real(CalcVectorViewLFM(R,Ch,Center,wm));
save('L3.dat',L3);clear('L3');
n=n+10000;R.P=RP(n,:);R.V2=RV2(n,:);R.V3=RV3(n,:);
L4=real(CalcVectorViewLFM(R,Ch,Center,wm));
save('L4.dat',L4);clear('L4');
n=40001:nq;R.P=RP(n,:);R.V2=RV2(n,:);R.V3=RV3(n,:);
L5=real(CalcVectorViewLFM(R,Ch,Center,wm));
save('L5.dat',L5);clear('L5');

保存したファイルを読み出しながら磁場導出行列を連結させます。
chn=size(Ch.P,2);
LFM=zeros(size(NodesPts,2)*2,chn);
load('L1.dat');n=1:10000;LFM([n,n+nq],:)=L1';clear('L1');
load('L2.dat');n=n+10000;LFM([n,n+nq],:)=L2';clear('L2');
load('L3.dat');n=n+10000;LFM([n,n+nq],:)=L3';clear('L3');
load('L4.dat');n=n+10000;LFM([n,n+nq],:)=L4';clear('L4');
load('L5.dat');n=40001:nq;LFM([n,n+nq],:)=L5';clear('L5');
save('L5.dat',LFM);
LL=zeros(chn,chn);
n=1:(chn/2);
L1=LFM(:,n);L2=LFM(:,n+chn/2);clear('LFM');
L1=L1';
LL(n,n+chn/2)=L1*L2;
L1=L1';L2=L2';
LL(n+chn/2,n)=L2*L1;
L2=L1';
LL(n,n)=L2*L1;
clear('L1');clear('L2');
load('L5.mat');
L1=LFM(:,n+chn/2);clear('LFM');
L2=L1';
LL(n+chn/2,n+chn/2)=L2*L1;
clear('L1');clear('L2');
load('L5.dat');L1=LFM(1:nq,:);L2=LFM(nq+1:$,:);clear('LFM');
相当面倒なことをしていますが、200MBのデータの行列の計算を分割しているためです。 尚本来作成したL1〜L5.datファイルを消すコマンドは mdelete()なのですが、壊れていて使えません。 explorerか何かで消してください。
磁場導出行列の特異値分解はメモリが足りず、
[U,S,V]=svd(LFM',0);
でも解決しないので諦めました。 そこで以下のように微小な単位行列を組み込むことで逆行列が求まるようにします。

LLI=inv(LL+mean(LL)/1000*eye(LL));
L1=L1*LLI;L2=L2*LLI;
clear('LL');clear('LLI');
// MEG data, noise data, variance.
meg=CalcSSP2(ave_0001(:,1:306));
meg=meg';time=ave_0001(:,$);

22msec、86msecの電流源推定を行います。
t=min(find(time>=22));
//t=min(find(time>=86));
q=((L1*meg(:,t)).^2+(L2*meg(:,t)).^2).^0.5;
q=abs(q)'*1e+9;//nA
scf();set(gcf(),'visible','off');
ncolor=32;
set(gcf(),'color_map',jetcolormap(ncolor));
subplot(121);
P=NodesPts(1:3,:);
param3d(P(1,:),P(2,:),P(3,:));
h=gca();h.box='off';
h.children.mark_mode='on';h.children.mark_size=1;
h.children.foreground=1;
for nn=1:ncolor-1;...
PP=P;n=1:nq;n(q < nn/ncolor*4)=[];PP=PP(:,n);...
if size(PP,2)==0;break;end; param3d(PP(1,:),PP(2,:),PP(3,:));...
h.children(1).mark_mode='on';h.children(1).mark_size=round(nn/32*5+3);...
h.children(1).foreground=nn+1;...
end;
h.rotation_angles=[70,200];
h.axes_visible='off';h.isoview='on';
subplot(122);Matplot((ncolor:-1:1)');
hh=gca();hh.axes_visible='off';hh.box='off';
hh.tight_limits='on';hh.axes_bounds=[0.75,0.2,0.05,0.6];
h.axes_bounds=[0,0,0.9,1];h.margins=[0,0,0,0];
set(gcf(),'visible','on');


青が0nA、赤が4nAです。 上記のγの値を大きくすれば電流は小さくなるようです。 特異値分解の際の次数に相当するような気がします。