[Master Index]
[Index for Toolbox]
eeg_sph
(Toolbox/eeg_sph.m in BrainStorm 2.0 (Alpha))
Function Synopsis
G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
Help Text
EEG_SPH - Calculate the electric potential , spherical head, arbitrary orientation
function G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
function G = eeg_sph(L,Channel,Param,Order);
L is 3 x nL, each column a source location
Channel is the channel structure, same for Param
Order is
-1 current dipole
0 focal(magnetic) dipole % NOT SUPPORTED
1 1st order multipole % NOT SUPPORTED
Verbose : toggle Verbose mode
Cross-Reference Information
This function calls
- dlegpoly C:\BrainStorm_2001\Toolbox\dlegpoly.m
- dotprod C:\BrainStorm_2001\Toolbox\dotprod.m
- good_channel C:\BrainStorm_2001\Toolbox\good_channel.m
- rownorm C:\BrainStorm_2001\Toolbox\rownorm.m
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\eeg_sph.m
function G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
%EEG_SPH - Calculate the electric potential , spherical head, arbitrary orientation
% function G = eeg_sph(L,Channel,Param,Order,Verbose,varargin);
% function G = eeg_sph(L,Channel,Param,Order);
% L is 3 x nL, each column a source location
% Channel is the channel structure, same for Param
% Order is
% -1 current dipole
% 0 focal(magnetic) dipole % NOT SUPPORTED
% 1 1st order multipole % NOT SUPPORTED
% Verbose : toggle Verbose mode
%<autobegin> ---------------------- 14-Jun-2004 17:10:26 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Forward Modeling
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\dlegpoly.m
% toolbox\dotprod.m
% toolbox\good_channel.m
% toolbox\rownorm.m
%
% Subfunctions in this file, in order of occurrence in file:
% G = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in,lam_berg_in)
%
% At Check-in: $Author: Mosher $ $Revision: 19 $ $Date: 6/14/04 3:37p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-2004
%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
% University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
% Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
% CNRS, Hopital de la Salpetriere, Paris, France
%
% See BrainStorm website at http://neuroimage.usc.edu for further information.
%
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 14-Jun-2004 17:10:26 -----------------------
% /---Script Authors--------------------------------------\
% | |
% | *** John Ermer, Ph.D. |
% | Signal % Image Processing Institute |
% | University of Southern California |
% | Los Angeles, CA, USA |
% | |
% | *** John C. Mosher, Ph.D. |
% | Biophysics Group |
% | |
% | *** Sylvain Baillet Ph.D. |
% | Cognitive Neuroscience & Brain Imaging Laboratory |
% | CNRS UPR640 - LENA |
% | Hopital de la Salpetriere, Paris, France |
% | sylvain.baillet@chups.jussieu.fr |
% | |
% \-------------------------------------------------------/
%
% Date of creation: October, 25 1999
%
% Script History -----------------------------------------------------------------------------------------------------------
% SB 19-Nov-2002 : Edited Header
% Updated management of EEG reference
% JCM 20-Nov-2002 : Fixed headers to have only one autocomments block
% SB 09-Mar-2004 : Added verbose mode
% --------------------------------------------------------------------------------------------------------------------------
if nargin < 5
Verbose = 1; % Default
end
%-----------------------------
nmax = 80;
%-----------------------------
% EEG Channels
EEGndx = good_channel(Channel,[],'EEG');
% Reference Channel
REFndx = good_channel(Channel,[],'EEG REF');
% EEG Channel locations
% Add EEG reference at the end of Channel and Param structures
Re = [Channel(EEGndx).Loc,Channel(REFndx).Loc]'; % Electrode location array
if length(Param) ~= length([EEGndx,REFndx])
Param(REFndx) = Param(EEGndx(1));
end
Param = Param([EEGndx,REFndx]);
Rq = L';
center = [Param.Center]';
Re = Re - center;
Rq = Rq - repmat(center(1,:),size(Rq,1),1); % Back to origin [0 0 0] for the sensors and the sources
clear tmp
switch(Param(1).EEGType)
case 'EEG_SINGLE'
R = Param(1).Radii(end);
sigma = Param(1).Conductivity(end);
otherwise
R = Param(1).Radii;
sigma = Param(1).Conductivity;
end
switch(Param(1).EEGType)
case 'EEG_BERG'
method = 2;
mu_berg_in = Param(1).Berg.mu;
lam_berg_in = Param(1).Berg.lam;
case 'EEG_3SHELL'
method = 1;
mu_berg_in = [];
lam_berg_in = [];
otherwise
method = 1;
mu_berg_in = [];
lam_berg_in = [];
end
if 1 % Projection of the EEG sensors on the sphere
[theta phi Re_sph] = cart2sph(Re(:,1),Re(:,2),Re(:,3));
Re_sph = R(end)*ones(size(Re_sph));
[Re(:,1) Re(:,2) Re(:,3)] = sph2cart(theta,phi,Re_sph);
end
Gtmp = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in',lam_berg_in');
if isempty(REFndx) % Average Reference
G = NaN * zeros(length(EEGndx),size(Gtmp,2));
%G(EEGndx,:) = Gtmp - repmat(mean(Gtmp),size(Gtmp,1),1);
G = Gtmp - repmat(mean(Gtmp),size(Gtmp,1),1);
clear Gtmp
else % Specific electrode as reference
% Remove its lead field from all others
% Note: reference leadfield is at the end of original G matrix (from gain_sph)
G = NaN * zeros(length(EEGndx),size(Gtmp,2));
%G(EEGndx,:) = Gtmp(1:end-1,:) - repmat(Gtmp(end,:),size(Gtmp,1)-1,1);
%G = G(EEGndx,:);
G = Gtmp(1:end-1,:) - repmat(Gtmp(end,:),size(Gtmp,1)-1,1);
% % Reference lead field is zero
% G(REFndx,:) = zeros(1,size(G,2));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Subfunctions -------------------------------------------------------------------------------
%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function G = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in,lam_berg_in)
%GAINP_SPH6X EEG Multilayer Spherical Forward Model
% function G = gainp_sph6x(Rq,Re,R,sigma,nmax,method,mu_berg_in,lam_berg_in)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% EEG MULTILAYER SPHERICAL FORWARD MODEL (gainp_sph6x.m)
%
% This function computes the voltage potential forward gain matrix for an array of
% EEG electrodes on the outermost layer of a single/multilayer conductive sphere.
% Each region of the multilayer sphere is assumed to be concentric with
% isontropic conductivity. EEG sensors are assumed to be located on the surface
% of the outermost sphere.
%
% Calculation of the electric potiential is performed using either (user-specified)
% of the following methods (Ref: Z. Zhang "A fast method to compute surface
% potentials generated by dipoles within multilayer anisotropic spheres"
% (Phys. Med. Biol. 40, pp335-349,1995)
%
% 1) Closed Form Solution (Single Shell Case Only). See formulas (1H,1H')
%
% 2) Series Expansion using Legendre Polynomials. See formulas (1I,2I,3I and 4I)
%
% 3) Series Approximiation of a Multilayer Sphere as three dipoles in a
% single shell using "Berg/Sherg" parameter approximation.
% See formulas (1i',5i" and 6i)
%
% Dipole generator(s) are assumed to be interior to the innermost "core" layer. For those
% dipoles external to the sphere, the dipole "image" is computed and used determine the
% gain function. The exception to this is the Legendre Method where all dipoles MUST be
% interior to the innermost "core" layer.
%
% INPUTS (Required):
% Rq : dipole location(in meters) P x 3
% Re : 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 scalar
% If not specified, a default value based on outermost
% dipole magnitude is computed. (Note: This parameter
% is ignored when Berg Approximation is commanded)
% method : Method used for computing forward potential
% 1=Legendre Series Approx; 2=Berg Parameter Approx
% (Note: Default and all other values invoke Legendre
% Series Approx. Exception is single-shell case where
% closed form solution is always used) scalar
% mu_berg_in: User specified initial value for Berg eccentricity
% factors (Required if Berg Method is commanded) 3 x 1
% lam_berg_in: User specified initial value for Berg magnitude
% factors (Required if Berg Method is commanded) 3 x 1
%
% WHERE: M=# of sensors; P=# of dipoles; NL = # of sphere layers
%
% OUTPUTS:
% G : EEG forward model gain matrix M x (3*P)
%
% External Functions and Files:
% dlegpoly.m; rownorm.m; dotprod.m: USC/LANL MEG/EEG Toolbox
% zhang_fit.m: External Function used to fit Berg Parameters (Zhang Eq# 5i")
%
% - John Ermer 6/3/99
% - 8/9/99: Modified to use EM image for dipoles external to brain
% (Applies to single-shell and Berg Methods only) (John Ermer)
% - 10/31/99: Optimized Processing Associated with Dipoles falling outside sphere
% (Applies to single-shell and Berg Methods only) (John Ermer)
% - 01/26/00: Corrected minor dimension error which caused program to fault when
% external dipoles and multiple sensors were present. (John Ermer)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% THIS PART CHECKS INPUT PARAMETERS FOR THEIR DIMENSION AND VALIDITY %%%
%
NL = length(R); % # of concentric sphere layers
P = size(Rq,1);
M = size(Re,1);
%
if R(1)~= min(R)
error('Head radii must be specified from innermost to outmost layer!!! ')
end
%
if size(Rq,2) ~= 3
error('Dipole location must have three columns!!!')
end
%
if nargin < 6, % Check # of input terms to see if method is specified
method = 1; % Default Method = Legendre Series Expansion
else
if (method>2)|(method<0)
method = 0; % Default Method = Legendre Series Expansion
end
end
%
%%% This part pre-initializes parameters used in future calculations %%%
%
G = zeros(M,3*P); % Pre-Allocate Gain Matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This part computes the potential for a dipole contained within a single-layer
%%% homogeneous sphere using closed form formula
%%% The EEG single-shell solution uses the vector form (which avoids computationally
%%% expensive intrinsic functions) described by Mosher et al ("EEG and MEG: Forward
%%% Solutions for inverse problems" IEEE BME Trans March 1999)
%
if NL == 1 % Single Shell Case (Closed Form Solution)
%
Re_mag = repmat(R(NL),P,M); %(PxM)
Re_mag_sq = repmat(R(NL)*R(NL),P,M); %(PxM)
Rq_mag = rownorm(Rq); %(Px1)
%
Rq1 = Rq; %(Px1)
Rq1_mag = Rq_mag; %(Px1)
Rq1_mag_sq = Rq_mag.*Rq_mag; %(Px1)
Re_dot_Rq1 = Rq1*Re'; %(PxM)
%
const = 4.0*pi*sigma(NL);
term = 1./(const*Rq1_mag_sq); %(Px1)
%
%%% This part checks for the presence of Berg dipoles which are external to
%%% the sphere. For those dipoles external to the sphere, the dipole parameters
%%% are replaced with the electrical image (internal to sphere) of the dipole
%
nx = find(Rq1_mag > R(NL));
%
if nx>0
Rq1_temp = Rq1(nx,:);
Rq1(nx,:) = R(NL)*R(NL)*Rq1_temp./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,3);
Rq1_mag(nx,1) = rownorm(Rq1(nx,:));
Rq1_mag_sq(nx,1) = Rq1_mag(nx,1).*Rq1_mag(nx,1);
Re_dot_Rq1(nx,:) = R(NL)*R(NL)*Re_dot_Rq1(nx,:)./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,M);
term(nx,:) = (R(NL)./rownorm(Rq1_temp)).*term(nx,:);
% was : term(nx,:) = (R(NL)/rownorm(Rq1_temp))'.*term(nx,:);
end
%
%%% Calculation of Forward Gain Matrix Contribution due to K-th Berg Dipole
%
Rq1_mag = repmat(Rq1_mag,1,M); %(PxM)
Rq1_mag_sq = repmat(Rq1_mag_sq,1,M); %(PxM)
term = repmat(term,1,M); %(PxM)
%
d_mag = reshape( rownorm(reshape(repmat(Re,1,P)',3,P*M)' ...
-repmat(Rq1,M,1)) ,P,M); %(PxM)
d_mag_cub = d_mag.*d_mag.*d_mag; %(PxM)
F_scalar = d_mag.*(Re_mag.*d_mag+Re_mag_sq-Re_dot_Rq1); %(PxM)
c1 = term.*(2*( (Re_dot_Rq1-Rq1_mag_sq)./d_mag_cub) ...
+ 1./d_mag - 1./Re_mag); %(PxM)
c2 = term.*((2./d_mag_cub) + (d_mag+Re_mag)./(Re_mag.*F_scalar));
%
G = G + reshape(repmat((c1 - c2.*Re_dot_Rq1)',3,1),M,3*P) ...
.*repmat(reshape(Rq1',1,3*P),M,1) ...
+ reshape(repmat((c2.*Rq1_mag_sq)',3,1),M,3*P) ...
.*repmat(Re,1,P);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%% This part computes the potential for a dipole contained within a multi-layer
%%% isontropic sphere using Legendre Polynomial Expansion (Zhang Eqs 1H, 1H')
%%% (Code based on gainp_sph.m by CCH, Aug/20/1995)
%
elseif (NL>1)&(method==1) % Multi-shell Case Using Legendre Expansion
%
Rq_mag = rownorm(Rq);
%
if ~all( Rq_mag < R(1)+eps ) % check if dipoles within the brain
errordlg('Legendre method assumes all dipoles(s) inside brain layer - please modify dipole locations OR use the "Berg" approach')
return
end
%
% compute weights fn. fn depends only on the radii and cdv
Rq_mag = rownorm(Rq);
Re_mag = R(NL); % Radius of outermost layer (Sensor distance from origin
%
%
if nargin < 5, % check # of inputs to see if nmax was specified
nmax = fix(10/(1-max(Rq_mag)/Re_mag)); % Default for # Legendre Series Terms
end
%
Ren = Re/Re_mag;
Rqn = Rq./[Rq_mag,Rq_mag,Rq_mag];
for k = 1:NL-1
s(k) = sigma(k)/sigma(k+1);
end
a = Re_mag./R;
ainv = R/Re_mag;
sm1 = s-1;
twonp1 = 2*[1:nmax]+1;
twonp1 = twonp1(:);
f = zeros(nmax,1);
%
for n = 1:nmax
np1 = n+1;
Mc = eye(2);
for k = 2:NL-1,
Mc = Mc*[n+np1*s(k), np1*sm1(k)*a(k)^twonp1(n);...
n*sm1(k)*ainv(k)^twonp1(n) , np1+n*s(k)];
end;
Mc(2,:) = [n*sm1(1)*ainv(1)^twonp1(n) , np1+n*s(1)]*Mc;
Mc = Mc/(twonp1(n))^(NL-1);
f(n) = n/(n*Mc(2,2)+np1*Mc(2,1));
end;
%
onevec = ones(M,1);
wtemp = ((twonp1./[1:nmax]').*f)/(4*pi*sigma(NL)*R(NL)^2);
n = [1:nmax]';
nm1 = n-1;
for i = 1:P, % loop over all dipoles
rqn = [Rqn(i,1)*onevec,Rqn(i,2)*onevec,Rqn(i,3)*onevec];
cosgamma = dotprod(rqn,Ren);
rqn = rqn(1,:);
[Pl,dP] = dlegpoly(nmax,cosgamma); % evaluate legendre poly and its derivative
ratio = (Rq_mag(i)/Re_mag).^nm1;
z = Ren- cosgamma*rqn;
w = wtemp.*ratio;
Gterm1 = Pl'*(w.*n);
Gterm2 = dP'*w;
G(:,3*i-2:3*i) = Gterm1*rqn + [z(:,1).*Gterm2,z(:,2).*Gterm2,z(:,3).*Gterm2];
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%----------------------------------------------------------------------------------%%
%%% This part computes the potential for a dipole contained within a multi-layer
%%% isontropic sphere using Berg Parameter Approximation (Zhang Eqs 1i',5i" and 6i)
%%% The EEG single-shell solution uses the vector form (which avoids computationally
%%% expensive intrinsic functions) described by Mosher et al ("EEG and MEG: Forward
%%% Solutions for inverse problems" IEEE BME Trans March 1999)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
elseif (NL>1)&(method==2)
%
if (~exist('mu_berg_in')|~exist('lam_berg_in'))
error('Berg Parameters have not been specified!!!')
elseif size(mu_berg_in)~=size(lam_berg_in)
error('Berg Parameters are of Unequal Lengths!!!')
else
J = length(mu_berg_in);
mu_berg = mu_berg_in;
lam_berg = lam_berg_in;
end
%
Re_mag = repmat(R(NL),P,M); %(PxM)
Re_mag_sq = repmat(R(NL)*R(NL),P,M); %(PxM)
Rq_mag = rownorm(Rq); %(Px1)
Rq_mag_sq = Rq_mag.*Rq_mag; %(Px1)
Re_dot_Rq = Rq*Re'; %(PxM)
%
for k=1:J
%
Rq1 = mu_berg(k)*Rq; %(Px3)
Rq1_mag = mu_berg(k)*Rq_mag; %(Px1)
Rq1_mag_sq = (mu_berg(k)*mu_berg(k))*Rq_mag_sq; %(Px1)
Re_dot_Rq1 = mu_berg(k)*Re_dot_Rq; %(PxM)
%
const = 4.0*pi*sigma(NL);
const1 = const/lam_berg(k);
term = 1./(const1*Rq1_mag_sq); %(PxM)
%
%%% This part checks for the presence of Berg dipoles which are external to
%%% the sphere. For those dipoles external to the sphere, the dipole parameters
%%% are replaced with the electrical image (internal to sphere) of the dipole
%
nx = find(Rq1_mag > R(NL));
%
if nx>0
Rq1_temp = Rq1(nx,:);
Rq1(nx,:) = R(NL)*R(NL)*Rq1_temp./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,3);
Rq1_mag(nx,1) = rownorm(Rq1(nx,:));
Rq1_mag_sq(nx,1) = Rq1_mag(nx,1).*Rq1_mag(nx,1);
Re_dot_Rq1(nx,:) = R(NL)*R(NL)*Re_dot_Rq1(nx,:)./repmat((rownorm(Rq1_temp).*rownorm(Rq1_temp)),1,M);
term(nx,:) = (R(NL)/rownorm(Rq1_temp))'.*term(nx,:);
end
%
%%% Calculation of Forward Gain Matrix Contribution due to K-th Berg Dipole
%
Rq1_mag = repmat(Rq1_mag,1,M); %(PxM)
Rq1_mag_sq = repmat(Rq1_mag_sq,1,M); %(PxM)
term = repmat(term,1,M); %(PxM)
%
d_mag = reshape( rownorm(reshape(repmat(Re,1,P)',3,P*M)' ...
-repmat(Rq1,M,1)) ,P,M); %(PxM)
d_mag_cub = d_mag.*d_mag.*d_mag; %(PxM)
%
F_scalar = d_mag.*(Re_mag.*d_mag+Re_mag_sq-Re_dot_Rq1); %(PxM)
%
c1 = term.*(2*( (Re_dot_Rq1-Rq1_mag_sq)./d_mag_cub) ...
+ 1./d_mag - 1./Re_mag); %(PxM)
c2 = term.*((2./d_mag_cub) + (d_mag+Re_mag)./(Re_mag.*F_scalar));
%
G = G + reshape(repmat((c1 - c2.*Re_dot_Rq1)',3,1),M,3*P) ...
.*repmat(reshape(Rq1',1,3*P),M,1) ...
+ reshape(repmat((c2.*Rq1_mag_sq)',3,1),M,3*P) ...
.*repmat(Re,1,P);
end
%
%---------------------------------------------------------------------------------------
%
end % End of Check for Forward Model Calculation Method
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON