% CSD - Current Source Density (CSD) transformation based on spherical spline
%       surface Laplacian as suggested by Perrin et al. (1989, 1990)
%
% (published in appendix of Kayser J, Tenke CE, Clin Neurophysiol 2006;117(2):348-368)
%
% Usage: [X, Y] = CSD(Data, G, H, lambda, head);
%
% 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:
%   Data = surface potential electrodes-by-samples matrix
%      G = g-function electrodes-by-electrodes matrix
%      H = h-function electrodes-by-electrodes matrix
% lambda = smoothing constant lambda (default = 1.0e-5)
%   head = head radius (default = no value for unit sphere [µV/m²])
%          specify a value [cm] to rescale CSD data to smaller units [µV/cm²]
%          (e.g., use 10.0 to scale to more realistic head size)
%
% Output parameter:
%      X = current source density (CSD) transform electrodes-by-samples matrix
%      Y = spherical spline surface potential (SP) interpolation electrodes-
%          by-samples matrix (only if requested)
%
% 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 $
%        - code compression and comments 
% Updated: $Date: 2007/02/07 11:30:00 $ $Author: jk $
%        - recommented rescaling (unit sphere [µV/m²] to realistic head size [µV/cm²])
%   Fixed: $Date: 2009/05/16 11:55:00 $ $Author: jk $
%        - memory claim for output matrices used inappropriate G and H dimensions
%   Added: $Date: 2009/05/21 10:52:00 $ $Author: jk $
%        - error checking of input matrix dimensions
%
function [X, Y] = CSD(Data, G, H, lambda, head)
[nElec,nPnts] = size(Data);            % get data matrix dimensions
if ~(size(G,1) == size(G,2)) | ...     % matrix dimension error checking: 
   ~(size(H,1) == size(H,2)) | ...     % G and H matrix must be nElec-by-nElec
   ~(size(G,1) == nElec) | ...         
   ~(size(H,1) == nElec)
   X = NaN; Y = NaN;                   % default to invalid output 
   help CSD
   disp(sprintf(strcat('*** Error: G (%d-by-%d) and H (%d-by-%d) matrix', ...
                       ' dimensions must match rows (%d) of data matrix'), ...
                       size(G),size(H),nElec));
   return 
end   
mu = mean(Data);                       % get grand mean
Z = (Data - repmat(mu,nElec,1));       % compute average reference
Y = Data;  X = Y;                      % claim memory for output matrices
if nargin < 5; head = 1.0; end;        % initialize scaling variable [µV/m²]
head = head * head;                    % or rescale data to head sphere [µV/cm²]
if nargin < 4; lambda = 1.0e-5; end;   % initialize smoothing constant
for e = 1:size(G,1);                   % add smoothing constant to diagonale
  G(e,e) = G(e,e) + lambda; 
end; 
Gi = inv(G);                           % compute G inverse
for i = 1:size(Gi,1);                  % compute sums for each row
  TC(i) = sum(Gi(i,:));
end;
sgi = sum(TC);                         % compute sum total
for p = 1:nPnts
  Cp = Gi * Z(:,p);                    % compute preliminary C vector
  c0 = sum(Cp) / sgi;                  % common constant across electrodes
  C = Cp - (c0 * TC');                 % compute final C vector
  for e = 1:nElec;                     % compute all CSDs ...
    X(e,p) = sum(C .* H(e,:)') / head; % ... and scale to head size
  end;
  if nargout > 1; for e = 1:nElec;     % if requested ...
    Y(e,p) = c0 + sum(C .* G(e,:)');   % ... compute all SPs
  end; end;
end;

% function [G,H] = GetGH (M, m)
% 
% Modification of GetCSD.m to obtain the G- and H-matrices, which are 
% invariant for a given EEG montage, which are then used as input with
% the EEG/ERP data to the CSD.m routine
%
% Usage: [G,H] = GetGH (M, m)
% 
% Input parameter:
%         M = a MatLab structure for an EEG montage consisting of <n> sites
%             as returned by the function ExtractEEGMontage.m;
%             M consists of:  <n> cell string array 'lab'
%                             <n> double 'theta'
%                             <n> double 'phi'
%                             <n * 2> double matrix 'xy'
%                             ('xy' values are not used in this function)
%         m = positive integer constant that affects spherical spline flexibility
%             values of m can be:    2 - most flexible interpolation
%                                    3 - medium flexibility
%                                    4 - more rigid splines (default)
%                                5..10 - increasingly rigid splines  
%
% Output parameters:
%         G = G electrodes-by-electrodes matrix (SP spline)
%         H = H electrodes-by-electrodes matrix (CSD spline)
%
% Notes: 1) GetGH.m differs from GetCSD.m by removing the first three lines
%           in the GetCSD.m code because spherical angles are defined through
%           <M>, and by removing the last line handling the EEG/ERP data
%        2) In contrast to GetCSD.m, the <m> constant may be optionally
%           specified as an input parameter (default: m = 4)
%        3) Using the G and H matrices with data not conforming to the original
%           EEG montage will not produce valid CSDs
%        4) The orientation of Theta and Phi is defined in the appendix of
%           Kayser J, Tenke CE, Clin Neurophysiol 2006;117(2):348-368), and
%           is not identical to other spherical notifications using the
%           same spherical angle names (e.g., as in EEGlab)
%
% -----------------------------------------------------------------------
% 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 six 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 $
% -----------------------------------------------------------------------
% Modified: $Date: 2009/05/12 19:12:00 $ $Author: jk $
%        - eliminate <ERPdata> input argument and direct use of CSD.m
%        - return G and H transformation matrices
%    Added: $Date: 2009/05/14 17:20:00 $ $Author: jk $
%        - input argument <m> to specify constant for spline flexibility
%    Added: $Date: 2009/05/22 11:28:00 $ $Author: jk $
%        - progress bar to indicate status of iteration 
%  Changed: $Date: 2010/04/11 19:46:00 $ $Author: jk $
%        - allow more rigid values for spline constant m [2..10] 
%
function [G,H] = GetGH (M, m)
ThetaRad = (2 * pi * M.theta) / 360;     % convert Theta and Phi to radians ...
PhiRad = (2 * pi * M.phi) / 360;         % ... and Cartesian coordinates ...
[X,Y,Z] = sph2cart(ThetaRad,PhiRad,1.0); % ... for optimal resolution
nElec = length(M.lab);                   % 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;
if nargin < 2
   m = 4;                                % set m constant default
end
if ~ismember(m,[2:10])                    % verify m constant
%   disp(sprintf('Error: Invalid m = %d [use 2, 3, or 4]',[m]));
   disp(sprintf('Error: Invalid m = %d [use an integer between 2 and 10]',[m]));
   G = NaN;
   H = NaN;
   return
end
disp(sprintf('Spline flexibility:  m = %d',[m])); 
N = 50;                                  % set N iterations
G(nElec,nElec) = 0; H(nElec,nElec) = 0;  % claim memory for G- and H-matrices
fprintf('%d iterations for %d sites [',N,nElec); % intialize progress bar
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
  if j == 1; fprintf('*'); end;          % show progress
  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;
disp(']');                               % finalize progress bar