% GetCSD - Compute Current Source Density (CSD) estimates using the spherical spline
%          surface Laplacian algorithm suggested by Perrin et al. (1989, 1990)
%
% (published in appendix of Kayser J, Tenke CE, Clin Neurophysiol 2006;117(2):348-368)
%
% Usage: [CE] = GetCSD(EEGmont, ERPdata);
%
% Implementation of algorithms described by Perrin, Pernier, Bertrand, and
% Echallier in Electroenceph Clin Neurophysiol 1989;72(2):184-187, and 
% Corrigenda EEG 02274 in Electroenceph Clin Neurophysiol 1990;76:565.
% 
% Input parameters:
%   EEGmont = name of EEG montage spherical coordinates file in ASCII format,
%             with rows consisting of seven columns: electrode label, two 
%             spherical angles (theta, phi), and Cartesian coordinates x,y,z
%   ERPdata = name of ERP data file in ASCII format stored as 
%             electrodes-by-samples (rows-by-columns) matrix
%
% Output parameter:
%        CE = CSD estimates as electrodes-by-samples matrix
%
% Copyright (C) 2003 by Jürgen Kayser (Email: kayserj@pi.cpmc.columbia.edu)
% GNU General Public License (http://www.gnu.org/licenses/gpl.txt)
% Updated: $Date: 2005/02/11 14:00:00 $ $Author: jk $
%
function [CE] = GetCSD(EEGmont, ERPdata)
ERP = load(ERPdata);                     % load electrodes-by-samples matrix
[Elab,Theta,Phi,X,Y,Z] = ...             % read EEG montage
   textread(EEGmont,'%s%f%f%f%f%f');
ThetaRad = (2 * pi * Theta) / 360;       % convert Theta and Phi to radians ...
PhiRad = (2 * pi * Phi) / 360;           % ... and Cartesian coordinates ...
[X,Y,Z] = sph2cart(ThetaRad,PhiRad,1.0); % ... for optimal resolution
nElec = length(Elab);                    % determine size of EEG montage
EF(nElec,nElec) = 0;                     % initialize interelectrode matrix ...
for i = 1:nElec; for j = 1:nElec;        % ... and compute all cosine distances
  EF(i,j) = 1 - ( ( (X(i) - X(j))^2 + ...
    (Y(i) - Y(j))^2 + (Z(i) - Z(j))^2 ) / 2 );
end; end;
m = 4; N = 50;                           % set m constant and N iterations
G(nElec,nElec) = 0; H(nElec,nElec) = 0;  % claim memory for G- and H-matrices
for i = 1:nElec; for j = 1:nElec;
  P = zeros(N);                          % compute Legendre polynomial
  for n = 1:N;
    p = legendre(n,EF(i,j));
    P(n) = p(1);
  end;
  g = 0.0; h = 0.0;                      % compute h- and g-functions
  for n = 1:N;
    g = g + ( (( 2.0*n+1.0) * P(n)) / ((n*n+n)^m    ) );
    h = h + ( ((-2.0*n-1.0) * P(n)) / ((n*n+n)^(m-1)) );
  end;
  G(i,j) =  g / 4.0 / pi;                % finalize cell of G-matrix
  H(i,j) = -h / 4.0 / pi;                % finalize cell of H-matrix
end; end;
CE = CSD(ERP, G, H);                     % compute CSD for ERP data