eeg_lap_sph_spline.m

Matlab下的EEG处理程序库

function [S,Li] = eeg_interp_sph_spline(V,X,Y,Z,Npoints)
% EEG_INTERP_SPH_SPLINE - Spherical Spline Laplacian of Potential
%
% Useage:   [S,Li] = eeg_lap_sph_spline(voltage,X,Y,Z,Npoints)
%
% where:    'voltage' is an EEG/ERP measurement at time t from
%           electrode positions (X,Y,Z) on a scalp surface.  All
%           input arrays are the same size, assumed (Nelec x 1).
%           The origin of X,Y,Z is assumed (0,0,0).
%
%           S  => the x,y,z points on a hemisphere
%           Li => spherical spline Laplacian of voltages
%
%           S is generated by 'elec_sphere_points' with 
%           Rpoints=Npoints and Epoints=16.
%
% Notes:    This function calculates the spherical spline Laplacian 
%           of Perrin et al (1989).  Electroenceph. & Clin. 
%             Neurophysiology, 72: 184-187.
%

% $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)
%
%           08/01  Needs testing & verification!
%                  With large electrode arrays, the result should be
%                  regularized (not sure why).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Check for correct size & orientation of X,Y,Z,V
    [x1,x2] = size(X);
    [y1,y2] = size(Y);
    [z1,z2] = size(Z); 
    [v1,v2] = size(V); 
    if ~and(isequal(x1,y1,z1,v1),isequal(x2,y2,z2,v2))
        error('ERROR: all X, Y, Z, V must be size (Nx1)'); 
    end
    if x1 < x2, X = X'; [x1,x2] = size(X); end
    if y1 < y2, Y = Y'; [y1,y2] = size(Y); end
    if z1 < z2, Z = Z'; [z1,z2] = size(Z); end
    if v1 < v2, V = V'; [v1,v2] = size(V); end
    if ~(isequal(x2,y2,z2,v2,1))
        error('ERROR: all X, Y, Z, V must be size (Nx1)');
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Computations derived from Perrin et al. (1989)

    % G * C = V,  where:
    %
    % G  (NxN)    spherical spline
    % C  (Nx1)    spline coefficients
    % V  (Nx1)    voltage potentials at (X,Y,Z) electrode locations
    %
    % solve for C = V * Sp' (C = V\Sp; (see "help slash"))

    % First calc G, where G is a function of the cosine of 
    % the angle (theta) between electrode point vectors
    A   = [ X Y Z ];
    COS = cosines(A,A);
    G   = spheric_spline(COS);
    C   = spline_coefficients(G,V); % spline coefficients (Co,C1,...,Cn)
    Co  = C(1);
    Ci  = C(2:end);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Obtain interpolated potentials at S (eq.1, Perrin et al., 1989)
    % V(s) = Co + sum( Ci * g(x) )

    % get spherical electrode radius
    [r,x,y,z] = elec_fit_sphere(X,Y,Z,0,0,0,0);
    % Generate spherical points for interpolation
    [x,y,z]   = elec_sphere_points(16,Npoints,r);
    S         = [x y z];
    
    fprintf('...Spherical Interpolation Progress:\n');
    rows = 80;  % progress indicator
    
    for p = 1:length(S)
        
        B = S(p,:);
        Cos = cosines(A,B);         %(1xN)
        
        Gx  = spheric_spline(Cos);  %(1xN)
        
        CiGx = Ci .* Gx';           %(1xN)
        
        Vi(p,1) = Co + sum(CiGx);
        
        fprintf('.');  % progress indicator
        if( p == rows ) fprintf('\n'); rows = rows + 80; end
    end
    
return





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve eq. 3 Perrin et al. (1989)
% g(COS) = 1/4pi * sum[n=1:inf] (( (2*n+1)/(n^m * (n+1)^m) ) * Pn(COS));

function [Gx] = spheric_spline(Cosine)
    
    m = 4;    N = 7;    % gives accuracy of 10^-6
    
    P = legendre(N,Cosine);
    %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,...).
    
    ndim = ndims(P);
    switch ndim
    case 2, P = P(2:N+1,:);
    case 3, P = P(2:N+1,:,:);
    case 4, P = P(2:N+1,:,:,:);
    otherwise
    end
    
    k = (1/4 * pi);
    
    for n = 1:(N),  Series(n,1) = (2*n + 1) / (n^(m-1) * (n+1)^(m-1));  end
    
    if min(size(Cosine)) == 1,    Gx      = k * ( Series' * P );
    else
        for i = 1:length(Cosine), Gx(i,:) = k * ( Series' * P(:,:,i) );
        end
    end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve eq. 2 Perrin et al. (1989)
function [C] = spline_coefficients(Gx,V)
    
    % add ones to first row & column of Gx
    tmp =  ones(size(V));       Gx   = [tmp Gx];
    tmp = [ones(size(V))' 1]';  Gx   = [tmp Gx']';
    
    Gx(1,1) = 0;    % according to Murk
    %Gx(1,1) = 1;    % according to EMSE, Greenblatt
    
    CoV = [0 V']';
    
    C = Gx\CoV;
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