function [norm] = electrodenormalize(cfg);

% ELECTRODENORMALIZE rotates and translates electrode positions to standard
% or averaged standard electrode positions.
%
% Use as
%   [elec] = electrodenormalize(cfg)
% The configuration can contain the following options
%   cfg.channel        = cell-array, see CHANNELSELECTION
%   cfg.casesensitive  = 'yes' (default) or 'no' determines whether string comparisons
%                        of electrode labels are case sensitive or not
%
% Different methods are implemented. You can apply a spatial deformation
% (i.e. 'warp') that minimizes the distance between the electrodes and the
% averaged standard. Alternatively, you can apply a rigid body realignment
% based on three fiducial locations.
%
%   cfg.method         = 'warp' or 'realignfiducial'
%   cfg.fiducial       = cell-array with the name of three fiducials used
%                        for realigning, default is {'nasion', 'lpa', 'rpa'}
%
% The electrode set of interest that will be normalized is specified as
%   cfg.elecfile       = string with filename, or alternatively
%   cfg.elec           = structure with electrode definition
%
% The electrode set(s) that serve as the standard can be specified either
% as electrode structures (i.e. when they are already read in memory) or as
% electrode files.
%
%   cfg.template       = single electrode set that serves as standard
%   cfg.template{1..N} = list of electrode sets that are averaged into the standard

% Copyright (C) 2005, Robert Oostenveld
%
% $Log: electrodenormalize.m,v $
% Revision 1.5  2005/05/17 17:50:37  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.4  2005/03/21 15:49:43  roboos
% added cfg.casesensitive for string comparison of electrode labels
% added cfg.feedback and cfg.plot3d option for debugging
% changed output: now ALL electrodes of the input are rerurned, after applying the specified transformation
% fixed small bug in feedback regarding distarnce prior/after realignfiducials)
% added support for various warping strategies, a.o. traditional, rigidbody, nonlin1-5, etc.
%
% Revision 1.3  2005/03/16 09:18:56  roboos
% fixed bug in fprintf feedback, instead of giving mean squared distance it should give mean distance before and after normalization
%
% Revision 1.2  2005/01/18 12:04:39  roboos
% improved error handling of missing fiducials
% added other default fiducials
% changed debugging output
%
% Revision 1.1  2005/01/17 14:56:06  roboos
% new implementation
%

% set the defaults
if ~isfield(cfg, 'method'),   cfg.method = 'unknown';    end
if ~isfield(cfg, 'channel'),  cfg.channel = 'all';       end
if ~isfield(cfg, 'feedback'), cfg.feedback = 'no';       end
if ~isfield(cfg, 'plot3d'),   cfg.plot3d = 'no';         end
if ~isfield(cfg, 'casesensitive'), cfg.casesensitive = 'yes'; end

% this is a common mistake which can be accepted
if strcmp(cfg.method, 'realignfiducials')
  cfg.method = 'realignfiducial';
end

if strcmp(cfg.feedback, 'yes')
  % use the global fb field to tell the warping toolbox to print feedback
  global fb
  fb = 1;
else
  global fb
  fb = 0;
end

% get the template electrode definitions
if ~iscell(cfg.template)
  cfg.template = {cfg.template};
end
Ntemplate = length(cfg.template);
for i=1:Ntemplate
  if isstruct(cfg.template{i})
    template(i) = cfg.template{i};
  else
    template(i) = read_fcdc_elec(cfg.template{i});
  end
end

% get the electrode definition that should be warped
if isfield(cfg, 'elec')
  elec = cfg.elec;
else
  elec = read_fcdc_elec(cfg.elecfile);
end

% remember the original electrode locations and labels
orig = elec;

% convert all labels to lower case for string comparisons
if strcmp(cfg.casesensitive, 'no')
  for i=1:length(elec.label)
    elec.label{i} = lower(elec.label{i});
  end
  for j=1:length(template)
    for i=1:length(template(j).label)
      template(j).label{i} = lower(template(j).label{i});
    end
  end
end

if strcmp(cfg.method, 'warp')
  % rename the default warp to one of the method recognized by the warping toolbox
  cfg.method = 'traditional';
end

if strcmp(cfg.plot3d, 'yes')
  % create an empty figure
  figure
  axis equal
  axis vis3d
  hold on
  xlabel('x')
  ylabel('y')
  zlabel('z')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if any(strcmp(cfg.method, {'rigidbody', 'globalrescale', 'traditional', 'nonlin1', 'nonlin2', 'nonlin3', 'nonlin4', 'nonlin5'}))
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % determine electrode selection and overlapping subset for warping
  cfg.channel = channelselection(cfg.channel, elec.label);
  for i=1:Ntemplate
    cfg.channel = channelselection(cfg.channel, template(i).label);
  end

  % make subselection of electrodes
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label = elec.label(datsel);
  elec.pnt   = elec.pnt(datsel,:);
  for i=1:Ntemplate
    [cfgsel, datsel] = match_str(cfg.channel, template(i).label);
    template(i).label = template(i).label(datsel);
    template(i).pnt   = template(i).pnt(datsel,:);
  end

  % compute the average of the template electrode positions
  all = [];
  for i=1:Ntemplate
    all = cat(3, all, template(i).pnt);
  end
  avg    = mean(all,3);
  stderr = std(all, [], 3);

  fprintf('warping electrodes to template... '); % the newline comes later
  [norm.pnt, norm.m] = warp_pnt(elec.pnt, avg, cfg.method);
  norm.label = elec.label;

  dpre  = mean(sqrt(sum((avg - elec.pnt).^2, 2)));
  dpost = mean(sqrt(sum((avg - norm.pnt).^2, 2)));
  fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);

  if strcmp(cfg.plot3d, 'yes')
    % plot all electrodes before warping
    my_plot3(elec.pnt, 'r.');
    my_plot3(elec.pnt(1,:), 'r*');
    my_plot3(elec.pnt(2,:), 'r*');
    my_plot3(elec.pnt(3,:), 'r*');
    my_text3(elec.pnt(1,:), elec.label{1}, 'color', 'r');
    my_text3(elec.pnt(2,:), elec.label{2}, 'color', 'r');
    my_text3(elec.pnt(3,:), elec.label{3}, 'color', 'r');

    % plot all electrodes after warping
    my_plot3(norm.pnt, 'm.');
    my_plot3(norm.pnt(1,:), 'm*');
    my_plot3(norm.pnt(2,:), 'm*');
    my_plot3(norm.pnt(3,:), 'm*');
    my_text3(norm.pnt(1,:), norm.label{1}, 'color', 'm');
    my_text3(norm.pnt(2,:), norm.label{2}, 'color', 'm');
    my_text3(norm.pnt(3,:), norm.label{3}, 'color', 'm');

    % plot the template electrode locations
    my_plot3(avg,      'b.');
    my_plot3(avg(1,:), 'b*');
    my_plot3(avg(2,:), 'b*');
    my_plot3(avg(3,:), 'b*');
    my_text3(avg(1,:), norm.label{1}, 'color', 'b');
    my_text3(avg(2,:), norm.label{2}, 'color', 'b');
    my_text3(avg(3,:), norm.label{3}, 'color', 'b');

    % plot lines connecting the input/warped electrode locations with the template locations
    my_line3(elec.pnt, avg, 'color', 'r');
    my_line3(norm.pnt, avg, 'color', 'm');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(cfg.method, 'realignfiducial')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % try to determine the fiducials automatically if not specified
  option1 = {'nasion' 'left' 'right'};
  option2 = {'nasion' 'lpa' 'rpa'};
  option3 = {'nz' 'lpa' 'rpa'};
  if ~isfield(cfg, 'fiducial')
    if length(match_str(elec.label, option1))==3
      cfg.fiducial = option1;
    elseif length(match_str(elec.label, option2))==3
      cfg.fiducial = option2;
    elseif length(match_str(elec.label, option3))==3
      cfg.fiducial = option3;
    else
      error('could not determine three fiducials, please specify cfg.fiducial')
    end
  end
  fprintf('using fiducials {''%s'', ''%s'', ''%s''}\n', cfg.fiducial{1}, cfg.fiducial{2}, cfg.fiducial{3});

  % determine electrode selection
  cfg.channel = channelselection(cfg.channel, elec.label);
  [cfgsel, datsel] = match_str(cfg.channel, elec.label);
  elec.label = elec.label(datsel);
  elec.pnt   = elec.pnt(datsel,:);

  if length(cfg.fiducial)~=3
    error('you must specify three fiducials');
  end

  % do case-insensitive search for fiducial locations
  nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
  if length(nas_indx)~=1 || length(lpa_indx)~=1 || length(rpa_indx)~=1
    error('not all fiducials were found in the electrode set');
  end
  elec_nas = elec.pnt(nas_indx,:);
  elec_lpa = elec.pnt(lpa_indx,:);
  elec_rpa = elec.pnt(rpa_indx,:);

  % find the matching fiducials in the template and average them
  templ_nas = [];
  templ_lpa = [];
  templ_rpa = [];
  for i=1:Ntemplate
    nas_indx = match_str(lower(template(i).label), lower(cfg.fiducial{1}));
    lpa_indx = match_str(lower(template(i).label), lower(cfg.fiducial{2}));
    rpa_indx = match_str(lower(template(i).label), lower(cfg.fiducial{3}));
    if length(nas_indx)~=1 || length(lpa_indx)~=1 || length(rpa_indx)~=1
      error(sprintf('not all fiducials were found in template %d', i));
    end
    templ_nas(end+1,:) = template(i).pnt(nas_indx,:);
    templ_lpa(end+1,:) = template(i).pnt(lpa_indx,:);
    templ_rpa(end+1,:) = template(i).pnt(rpa_indx,:);
  end
  templ_nas = mean(templ_nas,1);
  templ_lpa = mean(templ_lpa,1);
  templ_rpa = mean(templ_rpa,1);

  % realign both to a common coordinate system
  elec2common  = headcoordinates(elec_nas, elec_lpa, elec_rpa);
  templ2common = headcoordinates(templ_nas, templ_lpa, templ_rpa);

  % compute the combined transform and realign the electrodes to the template
  norm       = [];
  norm.m     = elec2common * inv(templ2common);
  norm.pnt   = warp3d(norm.m, elec.pnt, 'homogeneous');
  norm.label = elec.label;

  nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
  dpre  = mean(sqrt(sum((elec.pnt([nas_indx lpa_indx rpa_indx],:) - [templ_nas; templ_lpa; templ_rpa]).^2, 2)));
  nas_indx = match_str(lower(norm.label), lower(cfg.fiducial{1}));
  lpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{2}));
  rpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{3}));
  dpost = mean(sqrt(sum((norm.pnt([nas_indx lpa_indx rpa_indx],:) - [templ_nas; templ_lpa; templ_rpa]).^2, 2)));
  fprintf('mean distance between fiducials prior to realignment %f, after realignment %f\n', dpre, dpost);

  if strcmp(cfg.plot3d, 'yes')
    % plot the first three electrodes before transformation
    my_plot3(elec.pnt(1,:), 'r*');
    my_plot3(elec.pnt(2,:), 'r*');
    my_plot3(elec.pnt(3,:), 'r*');
    my_text3(elec.pnt(1,:), elec.label{1}, 'color', 'r');
    my_text3(elec.pnt(2,:), elec.label{2}, 'color', 'r');
    my_text3(elec.pnt(3,:), elec.label{3}, 'color', 'r');

    % plot the template fiducials
    my_plot3(templ_nas, 'b*');
    my_plot3(templ_lpa, 'b*');
    my_plot3(templ_rpa, 'b*');
    my_text3(templ_nas, ' nas', 'color', 'b');
    my_text3(templ_lpa, ' lpa', 'color', 'b');
    my_text3(templ_rpa, ' rpa', 'color', 'b');

    % plot all electrodes after transformation
    my_plot3(norm.pnt, 'm.');
    my_plot3(norm.pnt(1,:), 'm*');
    my_plot3(norm.pnt(2,:), 'm*');
    my_plot3(norm.pnt(3,:), 'm*');
    my_text3(norm.pnt(1,:), norm.label{1}, 'color', 'm');
    my_text3(norm.pnt(2,:), norm.label{2}, 'color', 'm');
    my_text3(norm.pnt(3,:), norm.label{3}, 'color', 'm');
  end

else
  error('unknown method');
end

% apply the spatial transformation to all electrodes, and replace the
% electrode labels by their case-sensitive original values
if any(strcmp(cfg.method, {'rigidbody', 'globalrescale', 'traditional', 'nonlin1', 'nonlin2', 'nonlin3', 'nonlin4', 'nonlin5'}))
  norm.pnt   = warp3d(norm.m, orig.pnt, cfg.method);
else
  norm.pnt   = warp3d(norm.m, orig.pnt, 'homogenous');
end
norm.label = orig.label;

% add version information to the configuration
try
  % get the full name of the function
  cfg.version.name = mfilename('fullpath');
catch
  % required for compatibility with Matlab versions prior to release 13 (6.5)
  [st, i] = dbstack;
  cfg.version.name = st(i);
end
cfg.version.id = '$Id: electrodenormalize.m,v 1.5 2005/05/17 17:50:37 roboos Exp $';

% remember the configuration
norm.cfg = cfg;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some simple SUBFUNCTIONs that facilitate 3D plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function h = my_plot3(xyz, varargin)
h = plot3(xyz(:,1), xyz(:,2), xyz(:,3), varargin{:});
function h = my_text3(xyz, varargin)
h = text(xyz(:,1), xyz(:,2), xyz(:,3), varargin{:});
function my_line3(xyzB, xyzE, varargin)
for i=1:size(xyzB,1)
  line([xyzB(i,1) xyzE(i,1)], [xyzB(i,2) xyzE(i,2)], [xyzB(i,3) xyzE(i,3)], varargin{:})
end
