function [scd] = ft_scalpcurrentdensity(cfg, data)

% FT_SCALPCURRENTDENSITY computes an estimate of the SCD using the
% second-order derivative (the surface Laplacian) of the EEG potential
% distribution
%
% The relation between the surface Laplacian and the SCD is explained
% in more detail on http://tinyurl.com/ptovowl.
%
% Use as
%   [data] = ft_scalpcurrentdensity(cfg, data)
% or
%   [timelock] = ft_scalpcurrentdensity(cfg, timelock)
% where the input data is obtained from FT_PREPROCESSING or from
% FT_TIMELOCKANALYSIS. The output data has the same format as the input
% and can be used in combination with most other FieldTrip functions
% such as FT_FREQANALYSIS or FT_TOPOPLOTER.
%
% The configuration should contain
%   cfg.method       = 'finite' for finite-difference method or
%                      'spline' for spherical spline method
%                      'hjorth' for Hjorth approximation method
%   cfg.elec         = structure with electrode positions or filename, see FT_READ_SENS
%   cfg.trials       = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.feedback     = string, 'no', 'text', 'textbar', 'gui' (default = 'text')
%
% The finite method require the following
%   cfg.conductivity = conductivity of the scalp (default = 0.33 S/m)
%
% The spline and finite method require the following
%   cfg.conductivity = conductivity of the scalp (default = 0.33 S/m)
%   cfg.lambda       = regularization parameter (default = 1e-05)
%   cfg.order        = order of the splines (default = 4)
%   cfg.degree       = degree of legendre polynomials (default for
%                       <=32 electrodes  = 9,
%                       <=64 electrodes  = 14,
%                       <=128 electrodes = 20,
%                       else             = 32
%
% The hjorth method requires the following
%   cfg.neighbours   = neighbourhood structure, see FT_PREPARE_NEIGHBOURS
%
% For the spline method you can specify the following
%   cfg.badchannel      = cell-array, see FT_CHANNELSELECTION for details (default = [])
%
% Note that the scalp conductivity, electrode dimensions and the potential
% all have to be expressed in the same SI units, otherwise the units of
% the SCD values are not scaled correctly. The spatial distribution still
% will be correct.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% The 'finite' method implements
%   TF Oostendorp, A van Oosterom; The surface Laplacian of the potential:
%   theory and application. IEEE Trans Biomed Eng, 43(4): 394-405, 1996.
%   G Huiskamp; Difference formulas for the surface Laplacian on a
%   triangulated sphere. Journal of Computational Physics, 2(95): 477-496,
%   1991.
%
% The 'spline' method implements
%   F. Perrin, J. Pernier, O. Bertrand, and J. F. Echallier.
%   Spherical splines for scalp potential and curernt density mapping.
%   Electroencephalogr Clin Neurophysiol, 72:184-187, 1989
% including their corrections in
%   F. Perrin, J. Pernier, O. Bertrand, and J. F. Echallier.
%   Corrigenda: EEG 02274, Electroencephalography and Clinical
%   Neurophysiology 76:565.
%
% The 'hjorth' method implements
%   B. Hjort; An on-line transformation of EEG scalp potentials into
%   orthogonal source derivation. Electroencephalography and Clinical
%   Neurophysiology 39:526-530, 1975.
%
% See also FT_PREPROCESSING, FT_TIMELOCKANALYSIS, FT_FREQNALYSIS, FT_TOPOPLOTER.

% Copyright (C) 2004-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'trial'}); % prevent accidental typos, see issue 1729

% set the defaults
cfg.method          = ft_getopt(cfg, 'method',       'spline');
cfg.conductivity    = ft_getopt(cfg, 'conductivity', 0.33); % in S/m
cfg.trials          = ft_getopt(cfg, 'trials',       'all', 1);
cfg.feedback        = ft_getopt(cfg, 'feedback',     'text');
cfg.badchannel      = ft_getopt(cfg, 'badchannel',     {});

switch cfg.method
  case 'hjorth'
    cfg = ft_checkconfig(cfg, 'required', {'neighbours'});
  case 'spline'
    cfg.lambda  = ft_getopt(cfg, 'lambda', 1e-5);
    cfg.order   = ft_getopt(cfg, 'order', 4);
    cfg.degree  = ft_getopt(cfg, 'degree', []);

    if isempty(cfg.degree) % determines degree of Legendre polynomials bases on number of electrodes
      nchan = numel(data.label);
      if nchan<=32
        cfg.degree = 9;
      elseif nchan<=64
        cfg.degree = 14;
      elseif nchan<=128
        cfg.degree = 20;
      else
        cfg.degree = 32;
      end
    end
  otherwise
    cfg = ft_checkconfig(cfg); % perform a simple consistency check
end

% store original datatype
dtype = ft_datatype(data);

% check if the input data is valid for this function
data = ft_checkdata(data, 'datatype', 'raw', 'feedback', 'yes', 'ismeg', []);

% get the electrode positions
tmpcfg = cfg;
tmpcfg.senstype = 'EEG';
elec = ft_fetch_sens(tmpcfg, data);

% select channels and trials of interest
tmpcfg = keepfields(cfg, {'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
tmpcfg.channel = elec.label;
data = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

Ntrials = numel(data.trial);

if isempty(cfg.badchannel)
  % check if the first sample of the first trial contains NaNs; if so treat it as a bad channel
  cfg.badchannel = ft_channelselection(find(isnan(data.trial{1}(:,1))), data.label);
end

% match the order of the data channels with the channel positions, order them according to the data
[datindx, elecindx] = match_str(data.label, elec.label);
[goodindx, tmp]     = match_str(data.label, setdiff(data.label, cfg.badchannel, 'stable'));

if ~isempty(cfg.badchannel)
  ft_info('detected channel %s as bad\n', cfg.badchannel{:});
  tmpcfg         = [];
  tmpcfg.channel = data.label(goodindx);
  data           = ft_selectdata(tmpcfg, data);
end

allchanpos  = elec.chanpos(elecindx,:); % the position of all channels, ordered according to the data
goodchanpos = allchanpos(goodindx,:);   % the position of good channels

% compute SCD for each trial
if strcmp(cfg.method, 'spline')
  fprintf('Checking spherical fit... ');
  [c, r] = fitsphere(allchanpos);
  d = allchanpos - repmat(c, size(allchanpos,1), 1);
  d = sqrt(sum(d.^2, 2));
  d = mean(abs(d) / r);
  if abs(d-1) > 0.1
    ft_warning('bad spherical fit (residual: %.2f%%). The interpolation will be inaccurate.', 100*(d-1));
  elseif abs(d-1) < 0.01
    fprintf('perfect spherical fit (residual: %.1f%%)\n', 100*(d-1));
  else
    fprintf('good spherical fit (residual: %.1f%%)\n', 100*(d-1));
  end
  % Builds the spatial filter only once.
  fprintf('Calculating the filter to build the SCD.\n');
  [WVo, WLo] = sphsplint(goodchanpos, allchanpos, cfg.order, cfg.degree, cfg.lambda);
  % Creates a montage to apply the spatial filter.
  montage.tra      = WLo;
  montage.labelold = elec.label(elecindx(goodindx));
  montage.labelnew = elec.label(elecindx);
  % Applies the montage to both the data and electrode definition
  scd  = ft_apply_montage(data, montage);
  elec = ft_apply_montage(elec, montage);

elseif strcmp(cfg.method, 'finite')
  if ~isempty(cfg.badchannel)
    ft_error('the method "%s" does not support the specification of bad channels', cfg.method);
  end
  % the finite difference approach requires a triangulation
  prj = elproj(allchanpos);
  tri = delaunay(prj(:,1), prj(:,2));
  % the new electrode montage only needs to be computed once for all trials
  montage.tra = lapcal(allchanpos, tri);
  montage.labelold = data.label;
  montage.labelnew = data.label;
  % apply the montage to the data, also update the electrode definition
  scd  = ft_apply_montage(data, montage);
  elec = ft_apply_montage(elec, montage);

elseif strcmp(cfg.method, 'hjorth')
  if ~isempty(cfg.badchannel)
    ft_error('the method "%s" does not support the specification of bad channels', cfg.method);
  end
  % convert the neighbourhood structure into a montage
  labelnew = {};
  labelold = {};
  for i=1:length(cfg.neighbours)
    labelnew = cat(2, labelnew, cfg.neighbours(i).label);
    labelold = cat(2, labelold, cfg.neighbours(i).neighblabel(:)');
  end
  labelold = cat(2, labelnew, labelold);
  labelold = unique(labelold);
  tra = zeros(length(labelnew), length(labelold));
  for i=1:length(cfg.neighbours)
    thischan   = match_str(labelold, cfg.neighbours(i).label);
    thisneighb = match_str(labelold, cfg.neighbours(i).neighblabel);
    tra(i, thischan) = 1;
    tra(i, thisneighb) = -1/length(thisneighb);
  end
  % combine it in a montage
  montage.tra = tra;
  montage.labelold = labelold;
  montage.labelnew = labelnew;
  % apply the montage to the data, also update the electrode definition
  scd  = ft_apply_montage(data, montage);
  elec = ft_apply_montage(elec, montage);

else
  ft_error('unknown method "%s"', cfg.method);
end

if strcmp(cfg.method, 'spline') || strcmp(cfg.method, 'finite')
  % correct the units
  ft_warning('trying to correct the units, assuming uV and mm');
  for trlop=1:Ntrials
    % The surface laplacian is proportional to potential divided by squared distance which means that, if
    % - input potential is in uV, which is 10^6 too large
    % - units of electrode positions are in mm, which is 10^3 too large
    % these two cancel out against each other. Hence the computed laplacian
    % is in SI units (MKS).
    scd.trial{trlop} = cfg.conductivity * -1 * scd.trial{trlop};
  end
  fprintf('output surface laplacian is in V/m^2\n');
else
  fprintf('output Hjorth filtered potential is in uV\n');
end

% Adds the electrode definition to the data.
scd.elec = elec;

% convert back to input type if necessary
switch dtype
  case 'timelock'
    scd = ft_checkdata(scd, 'datatype', 'timelock');
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous data

% rename the output variable to accommodate the savevar postamble
data = scd;

ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data
