eeg_interp_sph_spline.m

Matlab下的EEG处理程序库

function [Ej,U] = eeg_interp_sph_spline(Z,E,Ipoints)
% EEG_INTERP_SPH_SPLINE - Spherical Spline Interpolation of Potential
%
% Useage:   [Ej,U] = eeg_interp_sph_spline(voltage,Ei,Ipoints)
%
% where:    'voltage' is an EEG/ERP measurement at time t from
%           electrode positions Ei = [X Y Z] on a scalp surface.  All
%           input arrays are the same size, assumed (Nelec x 1).
%           The origin of Ei = [X Y Z] is assumed (0,0,0).
%
%           Ej => [x y z] points on a hemisphere (z > 0)
%           U  => spherical spline interpolated voltages at Ej
%
%           Ej = is generated by 'elec_sphere_points' with 
%                Rpoints=Ipoints and Epoints=16.
%
% Notes:    This function calculates the spherical spline of 
%           Perrin et al (1989).  Electroenceph. & Clin. 
%             Neurophysiology, 72: 184-187. Corrigenda (1990),
%             Electroenceph. & Clin. Neurophysiology, 76: 565.
%             (see comments in the .m file for details).

% $Revision: 1.2 $ $Date: 2003/03/02 03:20:43 $

% Licence:  GNU GPL, no implied or express warranties
% History:  08/01  [e-mail], with
%                  mathematical and matlab advice from
%                  Dr. Murk Bottema (Flinders University of SA)
%
%           02/02  Needs testing & verification!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

msg = sprintf('...EEG_INTERP_SPH_SPLINE: Needs further testing & verification!\n');
error(msg);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check for correct size & orientation of E & Z
[e1,e2] = size(E);
[v1,v2] = size(Z);
if e1 < e2, E = E.'; [e1,e2] = size(E); end
if v1 < v2, Z = Z.'; [v1,v2] = size(Z); end
if ~and(isequal(e1,v1),and(isequal(e2,3),isequal(v2,1)))
    error('ERROR: Ei must be Nx3 & Voltage must be Nx1 matrices');
end
n = e1; % The number of electrodes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computations described by Perrin et al. (1989)
%
% The value at any point (E) on a sphere is calculated by a spherical
% spline interpolation function (U) of the scalp potentials (z(i)) at
% all of i electrodes (i=1:n).  The spherical projection of any electrode
% i is denoted E(i) and we have (Eq. 1):
% 
% U(E) = Co + ( for i=1:n, sum = (sum + (C(i) * g(cos(E,E(i)))) )
% 
% The terms of this function are explained as follows.
% 
% The array C(i) are the spherical spline coefficients.  These are
% solutions of the simultaneous equations (Eq. 2):
% 
% a) GC + TCo = Z
% b) T'C = 0
% 
% with:

T =  ones(size(Z));   % T' = ones(1,n)
                      % C' = c(i), for i=1:n (unknown)
                      % Z' = z(i), for i=1:n (already input)

%           G = g(x) = g( cos(Ei,Ej) )
%
% where:    cos(Ei,Ej) is the cosine of the angle between the
%           spherical projections of electrodes Ei and Ej &
%           the function g(x) is the sum of the series (Eq. 3),
%
% g(x) = 1/4pi * for n=1:inf, sum = sum + ( ( (2*n+1)/(n^m * (n+1)^m) ) * Pn(x) );
%
% where m is a constant > 1 and Pn(x) is the nth degree Legendre
% polynomial.  Perrin et al. (1989) evaluated m=1:6 and recommend m=4,
% for which the first 7 terms of Pn(x) are sufficient to obtain
% a precision of 10^-6 for g(x) (ie, the above for loop is n=1:7).
% Perrin et al. (1989) recommend tabulation of g(x) for
% x = linspace(-1,1,2000) to be used as a lookup given actual
% values for cos(E(i),E(j)).

% first obtain spherical projections of Ei (assuming not already input)
fprintf('...EEG_INTERP_SPH_SPLINE: Projecting electrodes to sphere.\n');
[r,x,y,z] = elec_proj_sph(E(:,1),E(:,2),E(:,3));
Ei = [ x y z ]; clear x y z;

% Calculate the cosines, if Ei is Nx3 & Ej is Mx3, COS is NxM
fprintf('...EEG_INTERP_SPH_SPLINE: Calculating cosines.\n');
COS = cosines(Ei,Ei);

% Calc g(x) size NxM
fprintf('...EEG_INTERP_SPH_SPLINE: Calculating g(x).\n');
G = g(COS); clear COS;


% OK, now we have all the components to solve Eq. 2
% a) GC + TCo = Z
% b) T'C = 0

	% First add ones to last row & column of Gx
	GT = [G T];
	GT(n+1,:) = [T;1]';
	% Append 0 to Z (V becomes size(n+1,1))
	V = [Z;0];
	
	% solve for C + TCo = inv(G)*Z (or C = G\V; see "help slash")    
    C = GT\V; clear GT V;
    
	%A\B is the matrix division of A into B, which is roughly the
	%same as INV(A)*B, except it is computed by Gaussian elimination.
	%B/A is the matrix division of A into B, which is roughly the
	%same as B*INV(A)
    
	% Calculate inverse of Gx
	% Ginv = inv(GT);
	% EMSE NOTE: The inversion of G is sensitive to errors because 
	% electrode locations are not really on a sphere. Therefore, 
	% the solution should be regularized.  In EMSE v4.2, the 
	% regularization involves Singluar Value Decomposition:
	% G = UWV'
	% inv(G) approx = V W- U'
	%
	% W- = [w-(i,j)], w-(i,j) { 1/w(ij) for w(i,j) > cutoff
	%                           0       otherwise
	%
	% cutoff is determined by a smoothing parameter:
	%
	% cutoff = k * 10 ^ -lambda
	%
	% where k = { 10^-5 for laplacian
	%             10^-7 for voltage
	% and lambda is a the smoothing parameter.
	% Smaller values of lambda (<0), fewer terms are
	% truncated and the solution is less smooth.
	%
	%[U,W,S] = svd(GT);  %GT = U*W*S';
	%Winv = inv(W);
	%cutoff = 10^-7 * 10^0;
	%Winv(find(Winv<cutoff)) = 0;
	%GTinv = S * Winv * U.';

	Co = C(1);
	Ci = C(2:end);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Obtain interpolated potentials at E (Eq.1, Perrin et al., 1989)
% U(E) = c(0) + ( for i=1:n, sum = (sum + (c(i) * g(cos(E,E(i)))) )
% U(E) = c(0) + sum( Ci * g(x) )

% Generate spherical points for interpolation, Ej
fprintf('...EEG_INTERP_SPH_SPLINE: Generating spherical interpolation points.\n');
[x,y,z] = elec_sphere_points(16,Ipoints,r);
Ej      = [x y z]; clear x y z;

fprintf('...EEG_INTERP_SPH_SPLINE: Spherical Interpolation Progress:\n');
rows = 80;  % progress indicator

for e = 1:length(Ej)
    
    B = Ej(e,:);
    Cos = cosines(Ei,B); % Nx1
    
    Gx  = g(Cos);  % Nx1
    
    CiGx = Ci .* Gx'; % array multiplication!
    
    U(e,1) = Co + sum(CiGx);
    
    % progress indicator
    fprintf('.'); if isequal(e,rows), fprintf('\n'); rows = rows + 80; end
end
fprintf('\n');
return







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function [Gx] = g(X)
	% Solve Eq. 3 Perrin et al. (1989)
	%
	% g(x) = 1/4pi * (for n=1:inf, sum = sum + ( ( (2*n+1)/(n^m * (n+1)^m) ) * Pn(x) ) );
	%
	% where m is a constant > 1 and Pn(x) is the nth degree Legendre
	% polynomial.  Perrin et al. (1989) evaluated m=1:6 and recommend m=4,
	% for which the first 7 terms of Pn(x) are sufficient to obtain
	% a precision of 10^-6 for g(x) (ie, the above for loop is n=1:7).
    
    m = 4;    N = 7;
    
  	% Perrin et al. (1989) recommend tabulation of g(x) for
	% x = linspace(-1,1,2000) to be used as a lookup given actual
	% values for cos(Ei,Ej).
    if isempty(X),
        msg = sprintf('...Cosine matrix empty for g(x).\n');
        error(msg);
        %msg = sprintf('...Warning, cosine matrix empty, generating x.\n');
        %x = linspace(-1,1,2000);
    end
    
    P = legendre(N,X);
    %P = LEGENDRE(N,X) computes the associated Legendre functions 
    %of degree N and order m = 0, 1, ..., N, evaluated for each element
    %of X.
    %In general, P has one more dimension than X.
    %Each element P(m+1,i,j,k,...) contains the associated Legendre
    %function of degree N and order m evaluated at X(i,j,k,...).
    
    for n=1:N, Series(n) = (2*n + 1) / (n^m * (n+1)^m); end
    
    switch ndims(P)
    case 2,
        P = P(2:N+1,:);
        tmp = ones(N,size(X,1));
        for n=1:N, tmp(n,:) = Series(n) * P(n,:); end
    case 3,
        P = P(2:N+1,:,:);
        tmp = ones(N,size(X,1),size(X,2));
        for n=1:N, tmp(n,:,:) = Series(n) * P(n,:,:); end
    otherwise
    end
    %fprintf('size P = '); size(P)  % debug
    
    Gx = (1/(4*pi)) * squeeze(sum(tmp));
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Cos] = cosines(A,B)
    
    for     a = 1:size(A,1),  Aa = A(a,:); A_len = sqrt( sum(Aa.^2));
        for b = 1:size(B,1),  Bb = B(b,:); B_len = sqrt( sum(Bb.^2));

            if( Aa == Bb ) Cos(a,b) = 0;
            else           Cos(a,b) = dot(Aa,Bb) / (A_len * B_len);
            end
        end
    end
return