ft_rejectcomponent.m

function [data] = ft_rejectcomponent(cfg, comp, data)

% FT_REJECTCOMPONENT backprojects an ICA (or similar) decomposition to the
% channel level after removing the independent components that contain
% the artifacts. This function does not automatically detect the artifact
% components, you will have to do that yourself.
%
% Use as
%    [data] = ft_rejectcomponent(cfg, comp)
% or as
%    [data] = ft_rejectcomponent(cfg, comp, data)
%
% where the input comp is the result of FT_COMPONENTANALYSIS. The output
% data will have the same format as the output of FT_PREPROCESSING.
%
% An optional input argument data can be provided. In that case
% componentanalysis will do a subspace projection of the input data
% onto the space which is spanned by the topographies in the unmixing
% matrix in comp, after removal of the artifact components.  Please use
% this option of including data as input, if you wish to use the output
% data.grad in further computation, for example for leadfield computation.
%
% The configuration structure can contain
%   cfg.component  = list of components to remove, e.g. [1 4 7] or see FT_CHANNELSELECTION
%   cfg.demean     = 'no' or 'yes', whether to demean the input data (default = 'yes')
%   cfg.updatesens = 'no' or 'yes' (default = 'yes')
%
% 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.
%
% See also FT_COMPONENTANALYSIS, FT_PREPROCESSING

% Copyright (C) 2005-2014, 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 comp data
ft_preamble provenance comp data
ft_preamble trackconfig

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

% set the defaults
cfg.component       = ft_getopt(cfg, 'component',  []);
cfg.demean          = ft_getopt(cfg, 'demean',     'yes');
cfg.feedback        = ft_getopt(cfg, 'feedback',   'text');
cfg.updatesens      = ft_getopt(cfg, 'updatesens', 'yes');

% the data can be passed as input arguments or can be read from disk
hascomp = exist('comp', 'var');
hasdata = exist('data', 'var');

if hascomp && hasdata
  % check if the input data is valid for this function
  istlck  = ft_datatype(data, 'timelock');  % this will be temporary converted into raw
  data    = ft_checkdata(data, 'datatype', 'raw');
  comp    = ft_checkdata(comp, 'datatype', 'comp');
  label   = data.label;
  nchans  = length(data.label);
  ncomps  = length(comp.label);
elseif hascomp
  % check if the input data is valid for this function
  istlck  = ft_datatype(comp, 'timelock');  % this will be temporary converted into raw
  comp    = ft_checkdata(comp, 'datatype', 'raw+comp');
  label   = comp.topolabel;
  ncomps  = length(comp.label);
else
  ft_error('incorrect number of input arguments');
end

% cfg.component can be indicated by number or by label
cfg.component = ft_channelselection(cfg.component, comp.label);
reject = match_str(comp.label, cfg.component);

if isempty(reject)
  ft_warning('no components were selected for rejection');
end

if min(reject)<1
  ft_error('you cannot remove components that are not present in the data');
end

if max(reject)>ncomps
  ft_error('you cannot remove components that are not present in the data');
end

if hasdata && strcmp(cfg.demean, 'yes')
  % optionally perform baseline correction on each trial
  ft_info('baseline correcting data \n');
  for trial=1:numel(data.trial)
    data.trial{trial} = ft_preproc_baselinecorrect(data.trial{trial});
  end
end

% set the rejected component amplitudes to zero
ft_info('removing %d components\n', length(reject));
ft_info('keeping %d components\n',  ncomps-length(reject));

% create a projection matrix by subtracting the subspace spanned by the
% topographies of the to-be-removed components from identity
[seldat, selcomp] = match_str(label, comp.topolabel);

if hasdata && length(seldat)~=length(label)
  ft_warning('the subspace projection is not guaranteed to be correct for non-orthogonal components');
end

if hasdata
  mixing   = comp.topo(selcomp,:);
  unmixing = comp.unmixing(:,selcomp);
  
  % I am not sure about this, but it gives comparable results to the ~hasdata case
  % when comp contains non-orthogonal (=ica) topographies, and contains a complete decomposition
  
  montage     = [];
  montage.tra = eye(length(selcomp)) - mixing(:, reject)*unmixing(reject, :);
  % we are going from data to components, and back again
  montage.labelold = comp.topolabel(selcomp);
  montage.labelnew = comp.topolabel(selcomp);
  
  keepunused = 'yes'; % keep the original data which are not present in the mixing provided
  bname = 'reject';
  
else
  mixing = comp.topo(selcomp, :);
  mixing(:, reject) = 0;
  
  montage     = [];
  montage.tra = mixing;
  % we are going from components to data
  montage.labelold = comp.label;
  montage.labelnew = comp.topolabel(selcomp);
  
  keepunused = 'no'; % don't need to keep the original rejected components
  bname = 'invcomp';
  
  % create the initial data structure, remove all component details
  data = keepfields(comp, {'trial', 'time', 'label', 'fsample', 'grad', 'elec', 'opto', 'trialinfo', 'sampleinfo'});
end % if hasdata

% apply the linear projection to the data
data = ft_apply_montage(data, montage, 'keepunused', keepunused, 'feedback', cfg.feedback, 'showcallinfo', cfg.showcallinfo);

sensfield = cell(0,1);
if isfield(data, 'grad')
  sensfield{end+1} = 'grad';
end
if isfield(data, 'elec')
  sensfield{end+1} = 'elec';
end
if isfield(data, 'opto')
  sensfield{end+1} = 'opto';
end

% apply the linear projection also to the sensor description
if ~isempty(sensfield)
  if  strcmp(cfg.updatesens, 'yes')

    for m = 1:numel(sensfield)
      ft_info('also applying the backprojection matrix to the %s structure\n', sensfield{m});
      
      % the balance field is needed to keep the sequence of linear projections
      if ~isfield(data.(sensfield{m}), 'balance')
        data.(sensfield{m}).balance.current = 'none';
      end
      
      % keepunused = 'yes' is required to get back e.g. reference or otherwise
      % unused sensors in the sensor description. the unused components need to
      % be removed in a second step
      sens = ft_apply_montage(data.(sensfield{m}), montage, 'keepunused', 'yes', 'balancename', bname, 'feedback', cfg.feedback);
      
      % remove the unused channels from the grad/elec/opto
      [junk, remove]    = match_str(comp.label, sens.label);
      sens.tra(remove,:) = [];
      sens.label(remove) = [];
      sens.chanpos(remove,:) = [];
      if isfield(sens, 'chanori')
        sens.chanori(remove,:) = [];
      end
      
      % there could have been sequential subspace projections, so the
      % invcomp-field may have been renamed into invcompX. If this it the case,
      % take the one with the highest suffix
      invcompfield = bname;
      if  ~isfield(sens.balance, invcompfield)
        for k = 10:-1:1
          if isfield(sens.balance, [bname num2str(k)])
            invcompfield = [invcompfield num2str(k)];
            break;
          end

        end
      end
      
      % remove the unused components from the balancing
      [junk, remove]    = match_str(comp.label, sens.balance.(invcompfield).labelnew);
      sens.balance.(invcompfield).tra(remove, :)   = [];
      sens.balance.(invcompfield).labelnew(remove) = [];
      data.(sensfield{m})  = sens;
    end
    
  else
    for m = 1:numel(sensfield)
      ft_info('not applying the backprojection matrix to the %s structure\n', sensfield{m});
      % simply copy it over
      comp.(sensfield{m}) = data.(sensfield{m});
    end
  end
end % if sensfield

if istlck
  % convert the raw structure back into a timelock structure
  data = ft_checkdata(data, 'datatype', 'timelock');
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble trackconfig
ft_postamble previous   comp data
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data