EEG potential for dipoles
in multilayer spherical forward model
Legendre expansion calculation
eeg = vb_eeg_legendre(Xq, Q, Xe, R, sigma, Nmax)
INPUTS (Required):
Xq : dipole location(in meters) NP x 3
Q : dipole direction NP x 3
Xe : EEG sensors(in meters) on the scalp M x 3
R : radii(in meters) of sphere from
innermost to outermost NL x 1
sigma: conductivity from innermost to outermost NL x 1
INPUTS (Optional):
Nmax : # of terms used in Truncated Legendre Series
M = # of sensors
NP = # of dipoles
NL = # of sphere
OUTPUTS:
eeg : EEG Lead field matrix ( NP x M )
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL.46, 245-259, 1999
John C. Mosher, Richard M. Leahy and Paul S. Lewis
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function eeg = vb_eeg_legendre(Xq, Q, Xe, R, sigma, Nmax) 0002 % EEG potential for dipoles 0003 % in multilayer spherical forward model 0004 % Legendre expansion calculation 0005 % eeg = vb_eeg_legendre(Xq, Q, Xe, R, sigma, Nmax) 0006 % INPUTS (Required): 0007 % Xq : dipole location(in meters) NP x 3 0008 % Q : dipole direction NP x 3 0009 % Xe : EEG sensors(in meters) on the scalp M x 3 0010 % R : radii(in meters) of sphere from 0011 % innermost to outermost NL x 1 0012 % sigma: conductivity from innermost to outermost NL x 1 0013 % 0014 % INPUTS (Optional): 0015 % Nmax : # of terms used in Truncated Legendre Series 0016 % 0017 % M = # of sensors 0018 % NP = # of dipoles 0019 % NL = # of sphere 0020 % 0021 % OUTPUTS: 0022 % eeg : EEG Lead field matrix ( NP x M ) 0023 % 0024 % IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL.46, 245-259, 1999 0025 % John C. Mosher, Richard M. Leahy and Paul S. Lewis 0026 % 0027 % Copyright (C) 2011, ATR All Rights Reserved. 0028 % License : New BSD License(see VBMEG_LICENSE.txt) 0029 0030 M = size(Xe,1); % # of sensors 0031 NP = size(Xq,1); % # of dipoles 0032 NL = length(R); % # of sphere 0033 0034 % Radius of dipoles 0035 Rq = sqrt(sum(Xq.^2,2)); % NP x 1 0036 0037 % Radius of outermost layer (Sensor distance from origin) 0038 Rmax = R(NL); 0039 0040 if R(1)~= min(R) 0041 fprintf('Head radii must be R(1) < ... < R(NL)\n') 0042 return 0043 end 0044 0045 if ~all( Rq < R(1) ) 0046 % check if dipoles within the innermost layer 0047 fprintf('All dipoles must be inside innermost layer\n') 0048 return 0049 end 0050 0051 % EEG potential 0052 eeg = zeros(NP,M); 0053 0054 if nargin < 6, 0055 % Default # of Legendre expansion 0056 Nmax = min( fix(10/log(Rmax/max(Rq))), 200); 0057 end 0058 0059 % Normalized position (unit vector) 0060 Re = sqrt(sum(Xe.^2, 2)); 0061 Xe = Xe./Re(:,ones(1,3)); % EEG sensor ( M x 3 ) 0062 Xq = Xq./Rq(:,ones(1,3)); % dipole ( NP x 3 ) 0063 0064 % Weight factor of legendre expansion 0065 f = vb_multi_shell_param(R, sigma, Nmax); % Eq.(16) & (17) 0066 0067 npow = 2*[1:Nmax]+1; 0068 0069 % Coefficient of legendre expansion in Eq.(15) 0070 fn = ( (npow./[1:Nmax]).*f' )/(4*pi*sigma(NL)*R(NL)^2); % ( 1 x Nmax ) 0071 0072 n = [1:Nmax]; 0073 nm1 = n-1; 0074 0075 % loop for dipoles 0076 for i = 1:NP, 0077 Xqn = Xq(i,:)'; 0078 Qq = Q(i,:); 0079 0080 % angle between dipole moment and position vector 0081 cosalpha = Qq*Xqn; 0082 0083 % angle between dipole position and sensor position vectors 0084 cosgamma = Xe*Xqn; % ( M x 1 ) 0085 0086 % evaluate legendre polynomial and its derivative 0087 [Pn, dP] = vb_legendre_grad(Nmax,cosgamma); % ( Nmax x M ) 0088 0089 % power of radius 0090 ratio = (Rq(i)/Rmax).^nm1; 0091 0092 z = Xe*Qq' - cosgamma*cosalpha; % ( M x 1 ) 0093 0094 w = fn.*ratio; % (1 x Nmax) 0095 0096 eeg1 = (w.*n) * Pn; % (1 x M) 0097 eeg2 = w * dP; % (1 x M) 0098 0099 eeg(i,:) = eeg1.*cosalpha + eeg2.*z'; 0100 end 0101 0102 %eeg = eeg';