EEG potencial for multilayer model P = vb_bem_eeg_linear_grk(Peeg,Vs,Fs,V,J) OUTPUTS: P : EEG Lead field matrix ( M x NP) P : 多層モデルのEEG センサポテンシャル INPUTS (Required): V : dipole location(in meters) NP x 3 J : dipole moment NP x 3 Vs : vertex coordinate for boundary surfaces [Nvertex x 3] Fs : patch index for boundary surfaces [Npatch x 3] Peeg : transform matrix to map surface potencial to EEG [Nsensor x Nvertex] Peeg : 境界頂点係数行列 Vs : 境界三角面の座標 Fs : 境界三角面頂点インデックス V : 電流双極子座標 (m) J : 電流双極子モーメント 2004-02-06 Taku Yoshioka 2004-12-26 M. Sato modified Copyright (C) 2011, ATR All Rights Reserved. License : New BSD License(see VBMEG_LICENSE.txt)
0001 function P = vb_bem_eeg_linear_grk(Peeg,Vs,Fs,V,J) 0002 % EEG potencial for multilayer model 0003 % P = vb_bem_eeg_linear_grk(Peeg,Vs,Fs,V,J) 0004 % OUTPUTS: 0005 % P : EEG Lead field matrix ( M x NP) 0006 % P : 多層モデルのEEG センサポテンシャル 0007 % 0008 % INPUTS (Required): 0009 % V : dipole location(in meters) NP x 3 0010 % J : dipole moment NP x 3 0011 % Vs : vertex coordinate for boundary surfaces [Nvertex x 3] 0012 % Fs : patch index for boundary surfaces [Npatch x 3] 0013 % Peeg : transform matrix to map surface potencial to EEG [Nsensor x Nvertex] 0014 % 0015 % Peeg : 境界頂点係数行列 0016 % Vs : 境界三角面の座標 0017 % Fs : 境界三角面頂点インデックス 0018 % V : 電流双極子座標 (m) 0019 % J : 電流双極子モーメント 0020 % 0021 % 2004-02-06 Taku Yoshioka 0022 % 2004-12-26 M. Sato modified 0023 % 0024 % Copyright (C) 2011, ATR All Rights Reserved. 0025 % License : New BSD License(see VBMEG_LICENSE.txt) 0026 0027 % P0 : 境界頂点上のポテンシャル (Nvertex,Ndipole) 0028 0029 Nskip = 100; 0030 NSKIP = 1000; 0031 0032 % 変数の次元 0033 NJ = size(V,1); % 電流双極子 0034 NF = size(Fs,1); % 境界三角面 0035 NV = size(Vs,1); % 境界頂点 0036 NS = size(Peeg,1); % EEG センサ 0037 0038 % ポテンシャル 0039 P0 = zeros(NV,1); 0040 P = zeros(NS,NJ); % 境界面上のポテンシャル 0041 0042 % 三角面の頂点 triangle vertex 0043 x1 = Vs(Fs(:,1),:); 0044 x2 = Vs(Fs(:,2),:); 0045 x3 = Vs(Fs(:,3),:); 0046 0047 % 三角面の面積 0048 s0 = vb_cross2(x2-x1,x3-x1); 0049 S = sqrt(sum(s0.^2,2)); 0050 0051 %%%% 観測点でのポテンシャル計算 %%% 0052 0053 % 三角面の内点 0054 ip1 = (4*x1+x2+x3)/6; 0055 ip2 = (x1+4*x2+x3)/6; 0056 ip3 = (x1+x2+4*x3)/6; 0057 0058 for n=1:NJ, 0059 % n-番目の双極子 0060 x0 = V(n,:); 0061 J0 = J(n,:); 0062 0063 % 三角面内点上の双極子ポテンシャル 0064 p1 = vb_dipole_pot(J0,x0,ip1); 0065 p2 = vb_dipole_pot(J0,x0,ip2); 0066 p3 = vb_dipole_pot(J0,x0,ip3); 0067 0068 P0 = zeros(NV,1); 0069 0070 for i = 1:NF 0071 ix1 = Fs(i,1); 0072 ix2 = Fs(i,2); 0073 ix3 = Fs(i,3); 0074 P0(ix1) = P0(ix1)+S(i)/36*(4*p1(i)+p2(i)+p3(i)); 0075 P0(ix2) = P0(ix2)+S(i)/36*(p1(i)+4*p2(i)+p3(i)); 0076 P0(ix3) = P0(ix3)+S(i)/36*(p1(i)+p2(i)+4*p3(i)); 0077 end 0078 0079 % 境界面ポテンシャル 0080 P(:,n) = Peeg*P0; 0081 0082 if mod(n,Nskip)==0, fprintf('-'); end; 0083 if mod(n,NSKIP)==0, fprintf(' %3.0f %% done\n', 100*n/NJ); end; 0084 end;