% pop_icathresh() - main menu for choosing threshold for component
%                      rejection in EEGLAB. 
%
% Usage:
%   >> [OUTEEG rej] = pop_icathresh(INEEG, threshval, rejmethod, 
%                                      rejvalue, interact);
%
% Inputs:
%   INEEG      - input dataset
%   threshval  - values of thresholds for each of the 3 statistical
%                measures. Default is [] and the program uses the value
%                in the dataset.
%   rejmethod  - either 'percent', 'dataset' or 'current'. 'percent'
%                will reject a given percentage of components with
%                the highest value in one or several statistical 
%                measure. 'dataset' will use an other dataset for
%                calibration. 'current' will use the current dataset
%                for calibration. Default is 'current'.    
%   rejvalue   - percentage if rejmethod is 'percent', dataset number
%                if rejmethod is 'dataset' (no input if 'current'). If it
%                is a percentage, '25' for instance means that the 25%  
%                independent components with the highest values of one  
%                statistical measure are  rejected. Note that, one can 
%                also enter one percentage per statistical value (such as 
%                25 20 30).
%   interact   - interactive windows or just rejection
%
% Inputs:
%   OUTEEG     - output dataset with updated thresholds 
%   rej        - new rejection array   
%
% Graphic interface:
%   The graphic interface is divided into 3 parts. On the top, the 
% experimenter chooses the basis for the component rejection (see 
% rejmethod input). On the middle panel, the user can tune thresholds
% manually and plot the distribution of statistical values. Each time 
% that setting of one of these two top panels is modified, the curve
% on the third panel are redrawn.   
%   The third panel is composed of 3 graphs, one for each statistical 
% measure. The blue curve on each graph indicates the accuracy of the
% measure to detect artifactual components and non-artifactual 
% components of a dataset. Note that 'artifactual components are 
% defined with the rejection methods fields on the top (if you use 
% the percentage methods, these artifactual components may not actually
% be artifactual). For accurate rejection, one should first reject
% artifactual component manually and then set up the threshold to
% approximate this manual rejection. The green curve indicate the
% rejection when all measure are considered. Points on the blue and 
% on the green curves indicate the position of the current threshold on
% the parametric curve. Not that for the green cumulative curve, this
% point is at the same location for all graphs.
%   How to tune the threshold? To tune the threshold, one should try to 
% maximize the detection of non-artifactual components (because it is
% preferable to miss some artifactual components than to classify as
% artifactual non-artifact components). 
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eeglab()

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program 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 2 of the License, or
% (at your option) any later version.
%
% This program 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 this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%PROBLEM: when the % is set and we change manually yhe threshold,
% the software comes back to the percentage rejection

% 01-25-02 reformated help & license -ad 

function [EEG, rej, com] = pop_icathresh( EEG, threshval, rejmethod, rejvalue, interact); 

com = [];
rej = EEG.reject.gcompreject;
if nargin < 1
	help pop_icathresh;
	return;	
end;

if nargin < 2
	threshval = [];
end;
if nargin < 3
	rejmethod = 'current';
end;
if nargin < 4
	rejvalue = 25;
end;
if nargin < 5
	interact = 1;
end;
if ~isempty(threshval)
	EEG.reject.threshentropy = threshval(1);
	EEG.reject.threshkurtact = threshval(2);
	EEG.reject.threshkurtdist = threshval(3);
end;

tagmenu = 'pop_icathresh';

if ~isempty( findobj('tag', tagmenu))
	error('cannot open two identical windows, close the first one first');
end;

% the current rejection will be stored in userdata of the figure
% --------------------------------------------------------------
gcf = figure('visible', 'off', 'numbertitle', 'off', 'name', 'Choose thresholds', 'tag', tagmenu, 'userdata', rej);
pos = get(gca, 'position');
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)]./100;
axis off;

% definition of callbacks
% -----------------------

cb_cancel  = [  'userdat = get(gcbo, ''userdata'');' ... % restore thresholds
				'EEG.reject.threshentropy  = userdat{1};' ...
				'EEG.reject.threshkurtact  = userdat{2};' ...
				'EEG.reject.threshkurtdist = userdat{3};' ...
				'clear userdat; close(gcbf);' ];

drawgraphs  = [ 'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'if get( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'') == 1, ' ... % test if percentage
				'   perc = str2num(get( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''String''));' ... 
				'   if length(perc < 2), perc = [perc perc perc]; end;' ...
				'   perc = round((100-perc) * length( EEG.stats.compenta) / 100);' ... % convert the percentage
				'   method = zeros( size( EEG.stats.compenta ) );' ...
				'   [tmprej tmpindex] = sort(EEG.stats.compenta(:));' ...
				'   method( tmpindex(perc(1)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''entstring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold 
				'   [tmprej tmpindex] = sort(EEG.stats.compkurta(:));' ...
				'   method( tmpindex(perc(2)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''kurtstring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold 
				'   [tmprej tmpindex] = sort(EEG.stats.compkurtdist(:));' ...
				'   method( tmpindex(perc(3)+1:end) ) = 1;' ...
				'   set( findobj(''parent'', fig, ''tag'', ''kurtdiststring''), ''string'', num2str(tmprej(perc(1)+1)));' ... % update threshold 
				'   allvalues = [EEG.stats.compenta(:) EEG.stats.compkurta(:) EEG.stats.compkurtdist(:) ];' ...	
				'   clear perc tmprej tmpindex;' ...
                'end;' ...
				'if get( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'') == 1, ' ... % test if other dataset
				'   di = str2num(get( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''String''));' ... 
				'   if isempty( di ), clear fig; return; end;' ...
				'   method    = ALLEEG(di).reject.gcompreject'';' ... 
				'   allvalues = [ALLEEG(di).stats.compenta(:) ALLEEG(di).stats.compkurta(:) ALLEEG(di).stats.compkurtdist(:) ];' ...	
				'   clear di;' ...
                'end;' ...
				'if get( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'') == 1, ' ... % test if current dataset
				'   method = EEG.reject.gcompreject'';' ... 
				'   allvalues = [EEG.stats.compenta(:) EEG.stats.compkurta(:) EEG.stats.compkurtdist(:) ];' ...	
                'end;' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphent'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [0 EEG.reject.threshkurtact EEG.reject.threshkurtdist ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshentropy, ''Entropy of activity'',''Artifact detection (%)'', ''Non-artifact detection (%)'');' ...
				'set(gca, ''tag'', ''graphent'');' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphkurt'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [EEG.reject.threshentropy 0 EEG.reject.threshkurtdist ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshkurtact, ''Kurtosis of activity'', ''Artifact detection (%)'', '''');' ...
				'set(gca, ''tag'', ''graphkurt'');' ...
				'axes( findobj( ''parent'', fig, ''tag'', ''graphkurtdist'')); cla;' ...
                '[tmp1 tmp2 tmp3] = drawproba( method, allvalues, [EEG.reject.threshentropy EEG.reject.threshkurtact 0 ], ' ...
				'   { ''&'' ''|'' }, EEG.reject.threshkurtdist, ''Kurtosis of topography'', ''Artifact detection (%)'', '''');' ...
				'set(gca, ''tag'', ''graphkurtdist'');' ...
				'clear method allvalues fig;' ...
				];
                                    
cb_percent = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''on'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''off''); clear fig;' drawgraphs ];
				
cb_other   = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''off'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''on''); clear fig;' drawgraphs ];
				
cb_current = [  'if isempty( gcbf ), fig = gcf; else fig = gcbf; end;' ... % determine the figure handler
                'set( findobj(''parent'', fig, ''tag'', ''Icurrent''), ''value'', 1);'...
				'set( findobj(''parent'', fig, ''tag'', ''Iother''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercent''), ''value'', 0);' ...
				'set( findobj(''parent'', fig, ''tag'', ''Ipercenttext''), ''enable'', ''off'');' ...
				'set( findobj(''parent'', fig, ''tag'', ''Iothertext''), ''enable'', ''off''); clear fig;' drawgraphs ];

cb_entropy  = [ 'EEG.reject.threshentropy = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];
cb_kurtact  = [ 'EEG.reject.threshkurtact = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];
cb_kurtdist  = [ 'EEG.reject.threshkurtdist = str2num(get(gcbo, ''string''));' ...
				drawgraphs ];

if interact
	cb_calrej = [ 'ButtonName=questdlg2( ''This will erase previous projections'', ''Confirmation'', ''CANCEL'', ''OK'', ''OK'');' ]
else
	cb_calrej = [ 'ButtonName= ''OK''' ];	
end;
cb_calrej = [ cb_calrej ...
    		  'switch ButtonName,' ...
       	 	  '   case ''OK'',' ... 
			  '    rej1 = find(EEG.stats.compenta > EEG.reject.threshentropy);' ...
			  '    rej2 = find(EEG.stats.compkurta > EEG.reject.threshkurtact);' ...
			  '    rej3 = find(EEG.stats.compkurtdist > EEG.reject.threshkurtdist);' ...
			  '    EEG.reject.gcompreject = (rej1 & rej2) | rej3;' ...
			  '    clear rej1 rej2 rej3;' ...
			  'end; clear ButtonName;' ]; 

% default value for rejection methods
% -----------------------------------
rejvalother = ''; 
rejvalpercent = '25'; 
switch rejmethod
	case 'percent', rejvalpercent = num2str( rejvalue);
	case 'dataset', rejvalother = num2str( rejvalue); 
end;
    				
% -----------------------------------------------------
allh = supergui(gcf, { 	[1] ...
					[2 1] [2 1] [2 1] ...
					[1] ... 
					[1] ... 
					[1.5 0.5 1] [1.5 0.5 1] [1.5 0.5 1] ...
					[1] [1] [1] [1] [1] [1] [1] [1] [1] [1] ...
					[1 1 1 1 1] ...
					}, [], ...
	{ 'Style', 'text', 'string', 'Calibration method', 'FontSize', 13, 'fontweight', 'bold' }, ...
	...
	{ 'style', 'checkbox', 'String', '% of artifactual components (can also put one % per rejection)', 'tag', 'Ipercent', 'value', 1, 'callback', cb_percent}, ...
	{ 'style', 'edit', 'String', rejvalpercent, 'tag', 'Ipercenttext', 'callback', drawgraphs }, ...
	...
	{ 'style', 'checkbox', 'String', 'Specific dataset (enter dataset number)', 'tag', 'Iother', 'value', 0, 'callback', cb_other}, ...
	{ 'style', 'edit', 'String', rejvalother, 'tag', 'Iothertext', 'enable', 'off', 'callback', drawgraphs }, ...
	...
	{ 'style', 'checkbox', 'String', 'Current dataset', 'tag', 'Icurrent', 'value', 0, 'callback', cb_current}, ...
	{ }, ...
	...
	{ }, ...
	...
	{ 'Style', 'text', 'string', 'Threshold values', 'FontSize', 13, 'fontweight', 'bold' }, ...
	...
	{ 'style', 'text', 'String', 'Entropy of activity threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshentropy), 'tag', 'entstring', 'callback',  cb_entropy }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compenta,20); title(''Entropy of activity'');' }, ...
	...
	{ 'style', 'text', 'String', 'Kurtosis of activity threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshkurtact), 'tag', 'kurtstring', 'callback', cb_kurtact }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compkurta,200); title(''Kurtosis of activity'');' },  ...
	...
	{ 'style', 'text', 'String', 'Kurtosis of topography threshold' }, ...
	{ 'style', 'edit', 'String', num2str(EEG.reject.threshkurtdist), 'tag', 'kurtdiststring', 'callback', cb_kurtdist }, ...
	{ 'style', 'pushbutton', 'String', 'Histogram', 'callback', 'figure; hist(EEG.stats.compkurtdist,20); title(''Kurtosis of topography'');' }, ...
	...
	{ }, { }, { }, { }, { }, { }, { }, { }, { }, { }, ...
	...
	{ 'style', 'pushbutton', 'String', 'Cancel', 'callback', cb_cancel, 'userdata', { EEG.reject.threshentropy EEG.reject.threshkurtact EEG.reject.threshkurtdist }}, ...
	{ 'style', 'pushbutton', 'String', 'Auto thresh'  , 'callback', '', 'enable', 'off' }, ...
	{ 'style', 'pushbutton', 'String', 'Help'  , 'callback', 'pophelp(''pop_icathresh'');' }, ...
	{ 'style', 'pushbutton', 'String', 'Accept thresh.' , 'callback', 'close(gcbf);' }, ...
	{ 'style', 'pushbutton', 'String', 'Calc. rejection' , 'callback', 'close(gcbf);' } ...
	...
	);

h = axes('position', [0  15 30 30].*s+q, 'tag', 'graphent');
h = axes('position', [35 15 30 30].*s+q, 'tag', 'graphkurt');
h = axes('position', [70 15 30 30].*s+q, 'tag', 'graphkurtdist');

rejmethod
switch rejmethod
	case 'dataset', eval([ 'gcbf = [];' cb_other]);
	case 'current', eval([ 'gcbf = [];' cb_current]);
	otherwise, eval([ 'gcbf = [];' cb_percent]);
end;

	
return;	 
