% POP_EEGTHRESH - reject artifacts by detecting outlier values.  This has 
%                   long been a standard method for selecting data to reject.
%                   Applied either for electrode data or component activations.
% Usage:
%   >> pop_eegthresh( INEEG, typerej); % pop-up interactive window
%   >> [EEG Indexes] = pop_eegthresh( INEEG, typerej, elec_comp, lowthresh, ...
%                upthresh, starttime, endtime, superpose, reject);
%
% Graphic interface:
%   "Electrode|Component indices(s)" - [edit box] indices of the electrode(s) or 
%                 component(s) to take into consideration. Same as the 'elec_comp'
%                 parameter from the command line.
%   "Minimum rejection threshold(s)" - [edit box] lower threshold limit(s) 
%                 (in uV|std. dev.). Sets command line parameter 'lowthresh'.
%   "Maximum rejection threshold(s)" - [edit box] upper threshold limit(s) 
%                  (in uV|std. dev.). Sets command line parameter 'upthresh'.
%   "Start time limit(s)" - [edit box] starting time limit(s) (in seconds). 
%                 Sets command line parameter 'starttime'.
%   "End time limit(s)" - [edit box] ending time limit(s) (in seconds). 
%                 Sets command line parameter 'endtime'.
%   "Display previous rejection marks: " - [Checkbox]. Sets the command line
%                 input option 'eegplotplotallrej'.
%   "Reject marked trials: " - [Checkbox]  Sets the command line
%                 input option 'eegplotreject'.
%
% Inputs:
%   INEEG      - input EEG dataset
%   typerej    - type of rejection (0 = independent components; 1 = raw
%              data). Default is 1. For independent components, before
%              thresholding the activations are normalized (to have std. dev. 1).
%   elec_comp  - [e1 e2 ...] electrode|component numbers to take 
%              into consideration for rejection
%   lowthresh  - lower threshold limit (in uV|std. dev. For components, the 
%              threshold(s) are in std. dev.). Can be an array if more than one 
%              electrode|component number is given in elec_comp (above). 
%              If fewer values than the number of electrodes|components, the 
%              last value is used for the remaining electrodes|components. 
%   upthresh   - upper threshold limit (in uV|std dev) (see lowthresh above)
%   starttime  - rejection window start time(s) in seconds (see lowthresh above)
%   endtime    - rejection window end time(s) in seconds (see lowthresh)
%   superpose  - [0|1] 0=do not superpose rejection markings on previous
%              rejection marks stored in the dataset: 1=show both current and
%              previously marked rejections using different colors. {Default: 0}.
%   reject     - [1|0] 0=do not actually reject the marked trials (but store the 
%              marks: 1=immediately reject marked trials. {Default: 1}.
% Outputs:
%   Indexes    - index of rejected trials
%     When EEGPLOT is called, modifications are applied to the current 
%     dataset at the end of the call to EEGPLOT when the user presses 
%     the 'Reject' button.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: EEGTHRESH, EEGLAB, EEGPLOT, POP_REJEPOCH 

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad 
% 03-07-02 added srate argument to eegplot call -ad

function [EEG, Irej, com] = pop_eegthresh( EEG, icacomp, elecrange, negthresh, posthresh, ...
   						starttime, endtime, superpose, reject, topcommand)

Irej = [];
com = '';
if nargin < 1
   help pop_eegthresh;
   return;
end
if nargin < 2
   icacomp = 1;
end

if icacomp == 0
	if isempty( EEG(1).icasphere )
		disp('Error: you must run ICA first'); return;
	end
end
if exist('reject') ~= 1
    reject = 1;
end

if nargin < 3
    
    % which set to save
    % -----------------
    promptstr = { fastif(icacomp,'Electrode (indices(s), Ex: 2 4 5):'   , 'Component (indices, Ex: 2 6:8 10):'), ...
                  fastif(icacomp,'Minimum rejection threshold(s) (uV, Ex:-20 -10 -15):', 'Minimum rejection threshold(s) (std. dev,  Ex: -3 -2.5 -2):'), ...
                  fastif(icacomp,'Maximum rejection threshold(s) (uV, Ex: 20 10 15):'  , 'Maximum rejection threshold(s) (std. dev, Ex: 2 2 2.5):'), ...
                  'Start time limit(s) (seconds, Ex -0.1 0.3):', ...
                  'End time limit(s) (seconds, Ex 0.2):', ...
                  'Display previous rejection marks', ...
                  'Reject marked trial(s)' };
    inistr = { fastif(icacomp, ['1:' int2str(EEG(1).nbchan)], '1:5'), ...
               fastif(icacomp, '-10', '-20'),  ...
               fastif(icacomp, '10', '20'), ...
               num2str(EEG(1).xmin), ...
               num2str(EEG(1).xmax), ...
               '0', ...
               '0' };
    
    g1 = [1 0.1 0.75];
    g2 = [1 0.22 0.85];
    geometry = {g1 g1 g1 g1 g1 1 g2 g2};
    uilist = {...
              { 'Style', 'text', 'string', promptstr{1}} {} { 'Style','edit'      ,'string' ,inistr{1}  'tag' 'cpnum'}...
              { 'Style', 'text', 'string', promptstr{2}} {} { 'Style','edit'      ,'string' ,inistr{2} 'tag' 'lowlim' }...
              { 'Style', 'text', 'string', promptstr{3}} {} { 'Style','edit'      ,'string' ,inistr{3}  'tag' 'highlim'}...
              { 'Style', 'text', 'string', promptstr{4}} {} { 'Style','edit'      ,'string' ,inistr{4}  'tag' 'starttime'}...
              { 'Style', 'text', 'string', promptstr{5}} {} { 'Style','edit'      ,'string' ,inistr{5}  'tag' 'endtime'}...
              {}...
              { 'Style', 'text', 'string', promptstr{6}} {} { 'Style','checkbox'  ,'string'  ,' ' 'value' str2double(inistr{6}) 'tag','rejmarks' }...
              { 'Style', 'text', 'string', promptstr{7}} {} { 'Style','checkbox'  ,'string'  ,' ' 'value' str2double(inistr{7}) 'tag' 'rejtrials'} ...
               };
    figname = fastif(icacomp == 0, 'Rejection abnormal comp. values -- pop_eegthresh()','Rejection abnormal elec. values -- pop_eegthresh()');
    result = inputgui( geometry,uilist,'pophelp(''pop_eegthresh'');', figname);
    
    size_result = size(result);
    if size_result(1) == 0
        return;
    end
    elecrange    = result{1};
    negthresh    = result{2};
    posthresh    = result{3};
    starttime    = result{4};
    endtime      = result{5};
    superpose    = result{6};
    reject       = result{7};
end

if ischar(elecrange) % convert arguments if they are in text format
    calldisp = 1;
    elecrange = eval( [ '[' elecrange ']' ]  );
    negthresh = eval( [ '[' negthresh ']' ]  );
    posthresh = eval( [ '[' posthresh ']' ]  );
    if ischar(starttime)
        starttime = eval( [ '[' starttime ']' ]  );
    end
    if ischar(endtime)
        endtime   = eval( [ '[' endtime ']' ]  );
    end
else
    calldisp = 0;
end

% process multiple datasets
% -------------------------
if length(EEG) > 1
    if nargin < 2
        [ EEG, com ] = eeg_eval( 'pop_eegthresh', EEG, 'warning', 'on', 'params', { icacomp, elecrange, negthresh, posthresh, starttime, endtime, superpose, reject } );
    else
        [ EEG, com ] = eeg_eval( 'pop_eegthresh', EEG, 'params', { icacomp, elecrange, negthresh, posthresh, starttime, endtime, superpose, reject } );
    end
    Irej = [];
    return;
end

if any(starttime < EEG.xmin) 
 fprintf('Warning : starttime inferior to minimum time, adjusted\n'); 
	starttime(find(starttime < EEG.xmin)) = EEG.xmin; 
end
if any(endtime   > EEG.xmax) 
	fprintf('Warning : endtime superior to maximum time, adjusted\n'); 
	endtime(find(endtime > EEG.xmax)) = EEG.xmax;
end
if isempty(elecrange)
    elecrange = 1:EEG.nbchan;
end

if icacomp == 1
	[Itmp Irej NS Erejtmp] = eegthresh( EEG.data, EEG.pnts, elecrange, negthresh, posthresh, [EEG.xmin EEG.xmax], starttime, endtime);
    tmpelecIout = zeros(EEG.nbchan, EEG.trials);
    tmpelecIout(elecrange,Irej) = Erejtmp;
else
    icaacttmp = eeg_getdatact(EEG, 'component', elecrange);
	[Itmp Irej NS Erejtmp] = eegthresh( icaacttmp, EEG.pnts, 1:length(elecrange), negthresh, posthresh, [EEG.xmin EEG.xmax], starttime, endtime);
    tmpelecIout = zeros(size(EEG.icaweights,1), EEG.trials);
    tmpelecIout(elecrange,Irej) = Erejtmp;
end

fprintf('%d channel selected\n', size(elecrange(:), 1));
fprintf('%d/%d trials marked for rejection\n', length(Irej), EEG.trials);
tmprejectelec = zeros( 1, EEG.trials);
tmprejectelec(Irej) = 1;

rej  = tmprejectelec;
rejE = tmpelecIout;
if calldisp
    if icacomp == 1
        macrorej = 'EEG.reject.rejthresh';
        macrorejE = 'EEG.reject.rejthreshE';
    else
        macrorej  = 'EEG.reject.icarejthresh';
        macrorejE = 'EEG.reject.icarejthreshE';
    end
	
	colrej = EEG.reject.rejthreshcol;
	eeg_rejmacro; % script macro for generating command and old rejection arrays
	     
    if icacomp == 1
        eegplot( EEG.data(elecrange,:,:), 'srate', EEG.srate, 'limits', [EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:}); 
    else
        eegplot( icaacttmp, 'srate', EEG.srate, 'limits', [EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:}); 
    end
else 
    if reject == 1
        EEG = pop_rejepoch(EEG, rej, 0);
    end
end
if ~isempty(rej)
    if icacomp	== 1
        EEG.reject.rejthresh  = rej;
        EEG.reject.rejthreshE = rejE;
    else
        EEG.reject.icarejthresh  = rej;
        EEG.reject.icarejthreshE = rejE;
    end
end

%com = sprintf('Indexes = pop_eegthresh( %s, %d, [%s], [%s], [%s], [%s], [%s], %d, %d);', ...
%   inputname(1), icacomp, num2str(elecrange),  num2str(negthresh), ...
%   num2str(posthresh), num2str(starttime ) , num2str(endtime), superpose, reject ); 
com = [ com sprintf('EEG = pop_eegthresh(EEG,%s);', ...
		vararg2str({icacomp,elecrange,negthresh,posthresh,starttime,endtime,superpose,0})) ]; % reject is always set to 0 because trials are rejected in eegplot
if nargin < 3
	Irej = com;
end

% reject artifacts in a sequential fashion to save memory (ICA ONLY)
% -------------------------------------------------------
function [Irej, Erej] = thresh( data, elecrange, timerange, negthresh, posthresh, starttime, endtime);
    Irej    = [];
    Erej    = zeros(size(data,1), size(data,2));
    for index = 1:length(elecrange)
       tmpica = data(index,:,:);
       tmpica = reshape(tmpica, 1, size(data,2)*size(data,3));
       
       % perform the rejection
       % ---------------------	
	   tmpica = (tmpica-mean(tmpica,2)*ones(1,size(tmpica,2)))./ (std(tmpica,0,2)*ones(1,size(tmpica,2)));
	   [I1, Itmprej, NS, Etmprej] = eegthresh( tmpica, size(data,2), 1, negthresh, posthresh, timerange, starttime, endtime);
 	   Irej = union_bc(Irej, Itmprej);
 	   Erej(elecrange(index),Itmprej) = Etmprej;
	end