function varargin = zscore_gui(action,varargin)
%ZSCORE_GUI - : management of zscore visualization
% function varargin = zscore_gui(action,varargin)

%<autobegin> ---------------------- 26-May-2004 11:34:45 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: GUI and Related
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_color_scheme.m
%   toolbox\bst_imaging_display.m
%   toolbox\bst_layout.m
%   toolbox\bst_message_window.m
%   toolbox\dataplot_cb.m
%   toolbox\findclosest.m
%   toolbox\inorcol.m
%   toolbox\patch_swell.m
%   toolbox\save_fieldnames.m
%   toolbox\togglebuttoncolor.m
%   toolbox\vertices_connectivity.m
%   toolbox\zscore_gui.m  NOTE: Routine calls itself explicitly
%
% Application data and their calls in this file:
%   'TileType'
%   
%   setappdata(zscore_fig,'TileType','T')
%   
%
% Figure Files opened by this function:
%   'zscore_gui.fig'
%
%   Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
%   <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure zscore_gui.fig
%   uicontrol:edit:Baseline "" uses Callback for zscore_gui DefineBaseline
%   uicontrol:edit:CorrectionFactor "" uses Callback for zscore_gui ZScoreThreshold
%   uicontrol:edit:ZScoreThreshold "5" uses Callback for zscore_gui ZScoreThreshold
%   uicontrol:pushbutton:AnalyzeCorticalMap "Spatiotemp. Segment" uses Callback for
%     zscore_gui AnalyzeCorticalMap
%   uicontrol:pushbutton:AverageMaps "Avrg" uses Callback for zscore_gui AverageMaps
%   uicontrol:pushbutton:pushbutton55 "Load" uses Callback for zscore_gui LoadBaseline
%   uicontrol:pushbutton:pushbutton56 "Fetch files" uses Callback for zscore_gui ZScoreBatch
%   uicontrol:pushbutton:Qpushbutton7 "Close" uses Callback for close(gcbf)
%   uicontrol:pushbutton:STHistogram "Histo" uses Callback for zscore_gui STHistogram
%   uicontrol:togglebutton:IntSTH "Integ." uses Callback for zscore_gui IntegSTHistogram
%   uicontrol:togglebutton:ZScore "Z-score" uses Callback for togglebuttoncolor ; zscore_gui ZScore
%   uicontrol:togglebutton:ZScoreThresholdApply "Apply" uses Callback for
%     togglebuttoncolor; zscore_gui ZScoreThreshold
%
% At Check-in: $Author: Mosher $  $Revision: 12 $  $Date: 5/26/04 10:02a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 26-May-2004 11:34:45 -----------------------



global data

ImageGridTime = data.Display.Time(1):data.Display.SamplingRate:data.Display.Time(end);

switch(action)
	
case 'create'
	zscore_fig = openfig('zscore_gui.fig','reuse');
	setappdata(zscore_fig,'TileType','T')
	bst_color_scheme(zscore_fig)
	bst_layout('align',zscore_fig,2,2,4)
	set(zscore_fig,'visible','on')
	
	%----------------------------------------------------------------------------
	
case 'DefineBaseline'
	hZScore_fig = guihandles(gcbf);
	set(hZScore_fig.ZScore,'Value',1)
	dataplot_cb('ToggleButtonColor',hZScore_fig.ZScore)
	
	if ~strcmp(get(hZScore_fig.Baseline,'String'),'File')
		set(hZScore_fig.Baseline,'Userdata',0) % FLAG: Baseline defined on current file
	end
	
	togglebuttoncolor
	zscore_gui('ZScore',1) % Force ZScore computation
	
	%----------------------------------------------------------------------------
	
case 'LoadBaseline' % Define Baseline from a different result file
	cd(Users.STUDIES)
	[resfile, respath] = uigetfile('*result*.mat','Please Select a Result File Containing the Activation Baseline');
	if resfile == 0, return, end
	Zscore_fig = findobj(0,'Tag','tesselation_select');
	hZScore_fig = guihandles(Zscore_fig);
	set(hZScore_fig.LoadBaseline,'userdata',fullfile(respath, resfile));
	set(hZScore_fig.Baseline,'String','File');
	set(hZScore_fig.Baseline,'Userdata',1) % FLAG: Baseline defined on a pre-loaded file
	
	dataplot_cb DefineBaseline
	
	%----------------------------------------------------------------------------
	
case 'ZScore' % Switches from Z-Score to Absolute current density mapping 
	
	ZScore_fig = gcbf;
	hZScore_fig = guihandles(ZScore_fig);
	
	switch get(hZScore_fig.ZScore,'Value')
	case 1 % Switch to ZScore map 
		
		% Are we switching from ZScore to Absolute mapping or vice-versa ?
		if get(hZScore_fig.Baseline,'Userdata') == 1  % Baseline is defined from an existing result file, load it
			data.Results = load(get(hZScore_fig.LoadBaseline,'Userdata'));
			FlagBaselineFile = 1;
		else % else use current loaded results
			FlagBaselineFile = 0;
		end
		
		if isempty(data.Results), return, end % No results were loaded
		
		test = 0; % By default, consider not computing a new ZScore - check if a previously computation is available - 
		% if test is set to 1 anytime below, this will conduct to a new ZScore computation either because the baseline has changed 
		% or the computation is conducted on a new set of data.
		if nargin == 1
			varargin{1} = 0;
		end
		
		if (~isfield(data.Results,'ZScore') | varargin{1} == 1) % No ZScore was defined beforehand or computation is forced
			data.Results.ZScore.Baseline = [-Inf Inf]; % Dummy Score
			test = 1; % Force ZScore computation
		end
		
		%Get baseline time limits
		if FlagBaselineFile == 1
			set(hZScore_fig.Baseline,'String',...
				[num2str(ImageGridTime(1)*1000,'%3.1f'),' ',num2str(ImageGridTime(end)*1000,'%3.1f')]);
			
			baseline = [ImageGridTime(1) ImageGridTime(end)];
		else
			baseline = str2num(get(hZScore_fig.Baseline,'String'))/1000;
		end
		
		if length(baseline) < 2
			errordlg('Please enter 2 time values to define the baseline interval')
			return
		end
		
		if isfield(data.Results.ZScore,'FlagBaselineFile') % Was this ZScore computed from a stored baseline ?
			if data.Results.ZScore.FlagBaselineFile == 0 % No - check new baseline limits 
				% Get exact time samples closest to these values
				if round(1000*baseline(1)) < round(1000*ImageGridTime(1)) | round(1000*baseline(2)) > round(1000*ImageGridTime(end)) % Accuracy set to the ms
					errordlg('Please set baseline time extremas within the available time window for this result file')
					return    
				end
				
				[tmp, tmin] = min(abs(ImageGridTime-baseline(1)));
				[tmp, tmax] = min(abs(ImageGridTime-baseline(2)));
				
				baseline_new = [tmin tmax];
				
				test =  (sum(baseline_new == ...
					data.Results.ZScore.Baseline) ~= 2); % New baseline requested - recompute Z-Score
				
			elseif test ~= 1 % Yes - Baseline was stored in a file - do not recompute the ZScore
				
				[tmp, tmin] = min(abs(ImageGridTime-baseline(1)));
				[tmp, tmax] = min(abs(ImageGridTime-baseline(2)));
				
				test = 0;
				
			elseif varargin{1} == 1 % Force computation
				[tmp, tmin] = min(abs(ImageGridTime-baseline(1)));
				[tmp, tmax] = min(abs(ImageGridTime-baseline(2)));
				
				test = 1;
				
			end
			
		else
			[tmp, tmin] = min(abs(ImageGridTime-baseline(1)));
			[tmp, tmax] = min(abs(ImageGridTime-baseline(2)));
			
			test = 1;
		end
		
		if 1%test
			
			set(ZScore_fig,'Pointer','watch')
			
			bst_message_window('wrap',...
				{'Computing a new Z-Score map from this baseline. . .',...
					sprintf('Baseline is %d time samples large (out of %d samples available)',...
					tmax-tmin+1, length(ImageGridTime)),...
				})
			
			% Store original results patter in memory
			if 1%(~isfield(data.Results,'ImageGridAmpOrig') & FlagBaselineFile == 0) | varargin{1} == 1
				data.Results.ImageGridAmpOrig = data.Results.ImageGridAmp;
			end
			
			if FlagBaselineFile == 0 % remove the zero here !!!
				data.Results.ImageGridAmp = data.Results.ImageGridAmpOrig;
			end
			
			%data.Results.ImageGridAmp = abs(data.Results.ImageGridAmp);
			
			
			% Compute ZScores on that baseline
			data.Results.ZScore.Baseline = [tmin tmax];
			data.Results.ZScore.BaselineTime = [baseline(1) baseline(2)];
			
			data.Results.ImageGridAmp = double(data.Results.ImageGridAmp);
			
			% Average activity and its standard deviation over the baseline, for each source
			data.Results.ZScore.ImageGridZ.mean = ...
				mean((data.Results.ImageGridAmp(:,tmin:tmax)),2);
			data.Results.ZScore.ImageGridZ.std = ...
				std((data.Results.ImageGridAmp(:,tmin:tmax)),0,2);
			data.Results.ZScore.ImageGridZ.std = ...
				data.Results.ZScore.ImageGridZ.std + max(data.Results.ZScore.ImageGridZ.std)*eps; % Avoid devide-by-zero issues
			
			if ~isempty(get(hZScore_fig.CorrectionFactor,'String'))
				data.Results.ZScore.ImageGridZ.std = str2num(get(hZScore_fig.CorrectionFactor,'String'))*data.Results.ZScore.ImageGridZ.std;
			end
			
			% Compute Z-Scores
			if FlagBaselineFile == 1 % Now load current results
				ZScore = data.Results.ZScore; 
				data.Results = get(hZScore_fig.ResultFiles,'Userdata');
				data.Results.ZScore = ZScore; clear ZScore
				if ~isfield(data.Results,'ImageGridAmpOrig') 
					data.Results.ImageGridAmpOrig = [];
				end
				if isempty(data.Results.ImageGridAmpOrig)
					data.Results.ImageGridAmpOrig = data.Results.ImageGridAmp; % Store original current density map 
				end
				data.Results.ImageGridAmp = data.Results.ImageGridAmpOrig;
				data.Results.ZScore.BaselineFile = get(hZScore_fig.LoadBaseline,'Userdata');
				data.Results.ZScore.FlagBaselineFile = 1;
			else
				data.Results.ZScore.FlagBaselineFile = 0;
				data.Results.ZScore.BaselineFile = '';
			end
			
			% 			bst_message_window('Updating Result File...'), drawnow
			% 			% 			data.Results.ImageGridAmp = single(data.Results.ImageGridAmpOrig);
			% 			% 			tmp = data.Results.ImageGridAmpOrig;
			% 			% 			data.Results.ImageGridAmpOrig = [];
			% 			% 			save_fieldnames(data.Results,get(hZScore_fig.Refresh,'Userdata'));
			% 			ZScore = data.Results.ZScore;
			% 			Users = get_user_directory;
			% 			save(fullfile(Users.STUDIES,data.Results.FileName{1}),'ZScore','-append');
			% 			clear ZScore
			% 			%			data.Results.ImageGridAmpOrig = tmp; clear tmp
			%			data.Results.ImageGridAmp = double(data.Results.ImageGridAmp);
			
			data.Results.ImageGridAmp = (data.Results.ImageGridAmp)-repmat(data.Results.ZScore.ImageGridZ.mean,1,size(data.Results.ImageGridAmp,2));
			iStd = spdiags(1./data.Results.ZScore.ImageGridZ.std, 0, length(data.Results.ZScore.ImageGridZ.std), length(data.Results.ZScore.ImageGridZ.std)) ;
			
			%data.Results.ZScore.ImageGridZ.Amp = data.Results.ImageGridAmp;
			
			bst_message_window('Computing a new Z-Score map from this baseline...-> DONE'), drawnow
			
		else % Compute ZScore from saved mean and std maps
			
			data.Results.ImageGridAmpOrig = double(data.Results.ImageGridAmp);
			%data.Results.ImageGridAmp= data.Results.ZScore.ImageGridZ.Amp ;
			
			data.Results.ImageGridAmp = (data.Results.ImageGridAmp)-repmat(data.Results.ZScore.ImageGridZ.mean,1,size(data.Results.ImageGridAmp,2));
			iStd = spdiags(1./data.Results.ZScore.ImageGridZ.std, 0, length(data.Results.ZScore.ImageGridZ.std), length(data.Results.ZScore.ImageGridZ.std)) ;
			
		end
		
		%data.Results.ImageGridAmp = rand(size(iStd*data.Results.ImageGridAmp));
		data.Results.ImageGridAmp = single(abs(iStd*data.Results.ImageGridAmp));
		
		set(ZScore_fig,'Pointer','arrow')
		
		
	case 0 % Switch to absolute current density mapping
		
		bst_message_window('Previous Z-Score map on same baseline is used'), drawnow
		data.Results.ImageGridAmp = data.Results.ImageGridAmpOrig;
		set(hZScore_fig.ZScoreThresholdApply,'Value',0)
		dataplot_cb('ToggleButtonColor',hZScore_fig.ZScoreThresholdApply)
	end
	
	% Update display
	global hDataplot
	ctime = get(hDataplot.TimeSlider,'Value');
	tindx = findclosest(ctime,[data.Display.Time(1):data.Display.SamplingRate:data.Display.Time(end)]');
	
	
	if ishandle(data.Results.Display.Handles.hPatch) % surface handle for result display is valid
		data.Results.Display.OPTIONS.Time = tindx;
		[tmp,OPTIONS] = bst_imaging_display(data.Results.ImageGridAmp,data.Results.Display.OPTIONS); clear tmp
	end
	
	%----------------------------------------------------------------------------
case 'ZScoreThreshold'
	%Thresholded ZScore map in terms of multiples of sigmas in the baseline
	
	Zscore_fig = findobj(0,'Tag','tesselation_select');
	hZScore_fig = guihandles(Zscore_fig);
	TessWin  = findobj(0,'tag','tessellation_window'); 
	%if isempty(TessWin), return, end
	%hTessWin = guihandles(TessWin);
	
	%TessWinUserData = get(TessWin,'Userdata');
	
	% Are we switching from ZScore to Absolute mapping or vice-versa ?
	data.Results = get(hZScore_fig.ResultFiles,'Userdata');
	
	switch get(hZScore_fig.ZScoreThresholdApply,'Value')
	case 1 % Switch to ZScore map and threshold it
		if isempty(data.Results), return, end % No results were loaded
		
		if ~isfield(data.Results,'ZScore') % No ZScore was defined beforehand
			dataplot_cb ZScore % Compute the ZScore map
			data.Results = get(hZScore_fig.ResultFiles,'Userdata');
		end
		
		set(hZScore_fig.ZScore,'Value',1);
		dataplot_cb('ToggleButtonColor',hZScore_fig.ZScore)
		
		ZThres = str2num(get(hZScore_fig.ZScoreThreshold,'String'));
		if isempty(ZThres)
			set(hZScore_fig.ZScoreThreshold,'String',2); % Apply Default
			ZThres = str2num(get(hZScore_fig.ZScoreThreshold,'String'));
		end
		
		data.Results.ZScore.ImageGridZ.ZThres = ZThres;
		
		set(hZScore_fig.ZScoreThreshold,'Userdata',ZThres / max(data.Results.ImageGridAmp(:)))
		
	case 0
		% Don't apply thresholded map
		set(hZScore_fig.ZScoreThreshold,'Userdata',0)
		
	end
	
	
	%data.Results.ImageGridAmp(iZeroed) = eps ; %Thresholded Map
	
	% Update display
	set(hZScore_fig.ResultFiles,'Userdata',data.Results);
	
	dataplot_cb ScaleColormap
	%dataplot_cb tesselation_select_done
	%----------------------------------------------------------------------------
case 'IntegSTHistogram' % Integrate Saptio-temporal histogram overtime
	Zscore_fig = findobj(0,'Tag','tesselation_select');
	hZScore_fig = guihandles(Zscore_fig);
	
	avResults = get(hZScore_fig.ResultFiles,'Userdata');
	if ~isfield(avResults,'HistoSrc'), return, end % This is not an histogram result set
	
	% Integrate over time
	figure, imagesc(avResults.ImageGridTime, 1:size(avResults.ImageGridAmp,1),avResults.ImageGridAmp)
	colorbar
	avResults.ImageGridAmp = max(avResults.ImageGridAmp,[],2) *ones(1,size(avResults.ImageGridAmp,2));
	
	
	set(hZScore_fig.ResultFiles,'Userdata',avResults);
	
	
	%----------------------------------------------------------------------------
	
case 'STHistogram' % Spatio-temporal histogram of activations (in terms of ZScore) on a series of result files
	
	Zscore_fig = findobj(0,'Tag','tesselation_select');
	hZScore_fig = guihandles(Zscore_fig);
	
	% Look for data.Results file in the current study folder
	% Get the path
	[path,file,ext] = fileparts(Users.CurrentData.StudyFile);
	% Dir list and look for result files
	cd(fullfile(Users.STUDIES,path))
	ResFiles = dir('*result*.mat');
	for k = 1:length(ResFiles)
		ResFiles(k).name = strrep(ResFiles(k).name,'.mat','');
	end
	[ResSelect,ok] = listdlg('Liststring',{ResFiles.name},'Selectionmode','multiple','Name','Please select one or sevral result files',...
		'listsize',[400 400]);
	
	if ok == 0 
		return
	end
	
	hw = waitbar(0,['Computing the spatio-temporal histogram of activation maps through a set of ',int2str(length(ResSelect)),' result files...']);
	nDataFiles = 0;
	
	% Get Value for Z threshold
	ZThres = str2num(get(hZScore_fig.ZScoreThreshold,'String'));
	if isempty(ZThres)
		set(hZScore_fig.ZScoreThreshold,'String',2); % Apply Default
		ZThres = str2num(get(hZScore_fig.ZScoreThreshold,'String'));
	end
	
	
	for k = 1:length(ResSelect)
		dataplot_cb('LoadResultFile',ResFiles(ResSelect(k)).name)
		ResFiles(ResSelect(k)).name
		data.Results = get(hZScore_fig.ResultFiles,'Userdata');
		
		if k == 1
			avResults = data.Results;
			avResults.ImageGridAmp = 0* avResults.ImageGridAmp; % ImageGridAmp becomes a hit-count table - 
			% cell (i,j) will indicate how many times source i had a ZScore > ZThres at time j
			
			avResults.HistoSrc = cell(length(ResSelect));
		end
		
		if isfield(data.Results,'ZScore')
			
			%             % Adapt length to the length of shortest result file
			%             cropMin = 0;
			%             cropMax = 0;
			%             if data.Results.Time(1) > avResults.Time(1)       
			%                 cropMin = data.Results.Time(1)-avResults.Time(1);
			%                 avResults.ImageGridAmp = avResults.ImageGridAmp(:,cropMin+1:end);
			%                 avResults.ImageGridTime = avResults.ImageGridTime(:,cropMin+1:end);
			%                 avResults.Fsynth = avResults.Fsynth(:,cropMin+1:end);
			%                 avResults.ZScore.ImageGridZ.Amp = avResults.ZScore.ImageGridZ.Amp(:,cropMin+1:end);
			%                 
			%             end
			%             if data.Results.Time(1) < avResults.Time(1)       
			%                 cropMin = -data.Results.Time(1)+avResults.Time(1);
			%                 %                 data.Results.ImageGridAmp = data.Results.ImageGridAmp(:,cropMin+1:end);
			%                 data.Results.ZScore.ImageGridZ.Amp = data.Results.ZScore.ImageGridZ.Amp(:,cropMin+1:end);
			%                 data.Results.Time = data.Results.Time(cropMin+1:end);    
			%             end
			%             if data.Results.Time(end) > avResults.Time(end)       
			%                 cropMax = data.Results.Time(end)-avResults.Time(end);
			%                 %                 data.Results.ImageGridAmp = data.Results.ImageGridAmp(:,1:end-cropMax);
			%                 data.Results.ZScore.ImageGridZ.Amp = data.Results.ZScore.ImageGridZ.Amp(:,1:end-cropMax);
			%                 data.Results.Time = data.Results.Time(1:end-cropMax);    
			%             end
			%             if data.Results.Time(end) < avResults.Time(end)       
			%                 cropMax = -data.Results.Time(end)+avResults.Time(end);
			%                 avResults.ImageGridAmp = avResults.ImageGridAmp(:,1:end-cropMax);
			%                 avResults.ImageGridTime = avResults.ImageGridTime(:,1:end-cropMax);
			%                 avResults.Fsynth = avResults.Fsynth(:,1:end-cropMax);
			%                 avResults.ZScore.ImageGridZ.Amp = avResults.ZScore.ImageGridZ.Amp(:,1:end-cropMax);
			%             end
			%             
			avResults.Time = data.Results.Time;
			
			data.Results.ZScore.ImageGridZ.Amp = (data.Results.ImageGridAmp)-repmat(data.Results.ZScore.ImageGridZ.mean,1,size(data.Results.ImageGridAmp,2));
			iStd = spdiags(1./data.Results.ZScore.ImageGridZ.std, 0, length(data.Results.ZScore.ImageGridZ.std), length(data.Results.ZScore.ImageGridZ.std)) ;
			data.Results.ZScore.ImageGridZ.Amp = abs(iStd*data.Results.ZScore.ImageGridZ.Amp);
			
			avResults.HistoSrc{k} = sparse(data.Results.ZScore.ImageGridZ.Amp > ZThres);
			% avResults.HistoSrc{k} stores the source indices and time where sources had greater ZScore that ZThres
			avResults.ImageGridAmp = avResults.ImageGridAmp + avResults.HistoSrc{k} ;
			
			if k == 1
				prev = figure;
			end
			
			nDataFiles = nDataFiles + 1;    
			k
			figure(prev)
			plot(avResults.ImageGridTime, sum(avResults.HistoSrc{k},1))
			
			
		else % No ZScore available - skip file
			
			bst_message_window([ResFiles(ResSelect(k)).name,' does not contain any ZScore map - skip file'])
			
		end
		
		waitbar(k/length(ResSelect))
		
	end
	
	avResults.ImageGridAmp = 100*avResults.ImageGridAmp/nDataFiles; % Scale to 100% hits
	
	cd(Users.STUDIES);
	[path,file,ext] = fileparts(Users.CurrentData.StudyFile);
	cd(path)
	c = clock;
	ResFiletmp = strrep(ResFiles(ResSelect(end)).name,Users.STUDIES,'');
	ResFiletmp = strrep(ResFiletmp,'.mat','');
	Iunder = findstr(ResFiletmp,'_');
	DataFile= [ResFiletmp(1:Iunder(end)-1)];
	
	newname = ([DataFile sprintf('_HISTO_%02.0f%02.0f',c(4:5)) ext]);
	i = 0;
	while(exist(newname,'file')),
		i = i+1; % subtract another minute
		c(5) = mod(c(5) - 1,60);
		newname = ([DataFile sprintf('_HISTO_%02.0f%02.0f',MethodeCode,c(4),c(5)) ext]);
	end
	
	save_fieldnames(avResults,newname)
	
	delete(hw)
	
	msgbox(['Average activation maps saved in ', newname])
	set(hZScore_fig.ResultFiles,'Userdata',avResults)
	
	% Refresh results list
	dataplot_cb mesh_rendering
	ResFiles = get(hZScore_fig.ResultFiles,'String');
	iFile = find(strcmp(ResFiles,strrep(newname,'.mat','')));
	set(hZScore_fig.ResultFiles,'Value',iFile)
	
	
	
	%----------------------------------------------------------------------------
	
case 'AverageMaps' % Average activations across multiple result files
	
	Zscore_fig = findobj(0,'Tag','tesselation_select');
	hZScore_fig = guihandles(Zscore_fig);
	
	% Look for data.Results file in the current study folder
	% Get the path
	[path,file,ext] = fileparts(Users.CurrentData.StudyFile);
	% Dir list and look for result files
	cd(fullfile(Users.STUDIES,path))
	ResFiles = dir('*result*.mat');
	for k = 1:length(ResFiles)
		ResFiles(k).name = strrep(ResFiles(k).name,'.mat','');
	end
	[ResSelect,ok] = listdlg('Liststring',{ResFiles.name},'Selectionmode','multiple','Name','Please select one or sevral result files',...
		'listsize',[400 400]);
	
	if ok == 0 
		return
	end
	
	% Now compute ZScore for each of the selected files
	hw = waitbar(0,['Averaging activation maps through a set of ',int2str(length(ResSelect)),' result files...']);
	for k = 1:length(ResSelect)
		dataplot_cb('LoadResultFile',ResFiles(ResSelect(k)).name)
		ResFiles(ResSelect(k)).name
		data.Results = get(hZScore_fig.ResultFiles,'Userdata');
		
		if k == 1
			avResults = data.Results;
			avResults.ImageGridAmp = avResults.ImageGridAmp/length(ResSelect);    
		end
		
		% !!!!!!Suppose all result files were computed on the same time window - discard what's commented below
		
		%         % Adapt length to the length of shortes result file
		%         cropMin = 0;
		%         cropMax = 0;
		%         if data.Results.Time(1) > avResults.Time(1)       
		%             cropMin = data.Results.Time(1)-avResults.Time(1);
		%             avResults.ImageGridAmp = avResults.ImageGridAmp(:,cropMin+1:end);
		%             avResults.ImageGridTime = avResults.ImageGridTime(:,cropMin+1:end);
		%             avResults.Time = avResults.Time(:,cropMin+1:end);
		%             avResults.Fsynth = avResults.Fsynth(:,cropMin+1:end);
		%     
		%             if isfield(data.Results,'ZScore')
		%                 avResults.ZScore.ImageGridZ.Amp = avResults.ZScore.ImageGridZ.Amp(:,cropMin+1:end);
		%             end 
		%         end
		%         if data.Results.Time(1) < avResults.Time(1)       
		%             cropMin = -data.Results.Time(1)+avResults.Time(1);
		%             data.Results.ImageGridAmp = data.Results.ImageGridAmp(:,cropMin+1:end);
		%             ImageGridTime = ImageGridTime(:,cropMin+1:end);
		%             data.Results.Time = data.Results.Time(:,cropMin+1:end);
		%             if isfield(data.Results,'ZScore')
		%                 data.Results.ZScore.ImageGridZ.Amp = data.Results.ZScore.ImageGridZ.Amp(:,cropMin+1:end);
		%             end
		%         end
		%         if data.Results.Time(end) > avResults.Time(end)       
		%             cropMax = data.Results.Time(end)-avResults.Time(end);
		%             data.Results.ImageGridAmp = data.Results.ImageGridAmp(:,1:end-cropMax);
		%             ImageGridTime = ImageGridTime(:,1:end-cropMax);
		%             data.Results.Time = data.Results.Time(:,1:end-cropMax);
		%             if isfield(data.Results,'ZScore')
		%                 data.Results.ZScore.ImageGridZ.Amp = data.Results.ZScore.ImageGridZ.Amp(:,1:end-cropMax);
		%             end
		%         end
		%         if data.Results.Time(end) < avResults.Time(end)       
		%             cropMax = -data.Results.Time(end)+avResults.Time(end);
		%             avResults.ImageGridAmp = avResults.ImageGridAmp(:,1:end-cropMax);
		%             avResults.ImageGridTime = avResults.ImageGridTime(:,1:end-cropMax);
		%             avResults.Fsynth = avResults.Fsynth(:,1:end-cropMax);
		%             avResults.Time = avResults.Time(:,1:end-cropMax);
		%             if isfield(data.Results,'ZScore')
		%                 avResults.ZScore.ImageGridZ.Amp = avResults.ZScore.ImageGridZ.Amp(:,1:end-cropMax);
		%             end
		%         end 
		
		if k == 1
			prev = figure;
			if isfield(data.Results,'ZScore')
				data.Results.ZScore.ImageGridZ.Amp = (data.Results.ImageGridAmp)-repmat(data.Results.ZScore.ImageGridZ.mean,1,size(data.Results.ImageGridAmp,2));
				iStd = spdiags(1./data.Results.ZScore.ImageGridZ.std, 0, length(data.Results.ZScore.ImageGridZ.std), length(data.Results.ZScore.ImageGridZ.std)) ;
				data.Results.ZScore.ImageGridZ.Amp = abs(iStd*data.Results.ZScore.ImageGridZ.Amp);
			end
		end
		
		if k >1
			avResults.ImageGridAmp = avResults.ImageGridAmp + (data.Results.ImageGridAmp)/length(ResSelect);
			if isfield(data.Results,'ZScore')
				avResults.ZScore.ImageGridZ.Amp = avResults.ZScore.ImageGridZ.Amp + (data.Results.ZScore.ImageGridZ.Amp)/length(ResSelect);
			end
		end
		figure(prev)
		plot(ImageGridTime,mean(abs(data.Results.ImageGridAmp),1),'b'), hold on
		plot(avResults.ImageGridTime,mean(abs(avResults.ImageGridAmp),1),'r'), hold off
		
		waitbar(k/length(ResSelect))
		
	end
	
	%avResults.ImageGridAmp = avResults.ImageGridAmp/length(ResSelect);
	%     if isfield(avResults,'ZScore')
	%         avResults.ZScore.ImageGridZ.Amp = avResults.ZScore.ImageGridZ.Amp/length(ResSelect);
	%     end
	
	cd(Users.STUDIES);
	[path,file,ext] = fileparts(Users.CurrentData.StudyFile);
	cd(path)
	c = clock;
	ResFiletmp = strrep(ResFiles(ResSelect(end)).name,Users.STUDIES,'');
	ResFiletmp = strrep(ResFiletmp,'.mat','');
	Iunder = findstr(ResFiletmp,'_');
	DataFile= [ResFiletmp(1:Iunder(end-1)-1)];
	
	newname = ([DataFile sprintf('_AVRresults_%02.0f%02.0f',c(4:5)) ext]);
	i = 0;
	while(exist(newname,'file')),
		i = i+1; % subtract another minute
		c(5) = mod(c(5) - 1,60);
		newname = ([DataFile sprintf('_AVRresults_%02.0f%02.0f',MethodeCode,c(4),c(5)) ext]);
	end
	
	save_fieldnames(avResults,newname)
	
	delete(hw)
	
	msgbox(['Average activation maps saved in ', newname])
	set(hZScore_fig.ResultFiles,'Userdata',avResults)
	
	% Refresh results list
	dataplot_cb mesh_rendering
	ResFiles = get(hZScore_fig.ResultFiles,'String');
	iFile = find(strcmp(ResFiles,strrep(newname,'.mat','')));
	set(hZScore_fig.ResultFiles,'Value',iFile)
	
	%----------------------------------------------------------------------------
case 'ZScoreBatch' % Compute ZScore on a set of Result files
	
	% Look for data.Results file in the current study folder
	% Get the path
	[path,file,ext] = fileparts(Users.CurrentData.StudyFile);
	% Dir list and look for result files
	cd(fullfile(Users.STUDIES,path))
	ResFiles = dir('*result*.mat');
	for k = 1:length(ResFiles)
		ResFiles(k).name = strrep(ResFiles(k).name,'.mat','');
	end
	[ResSelect,ok] = listdlg('Liststring',{ResFiles.name},'Selectionmode','multiple','Name','Please select one or sevral result files',...
		'listsize',[400 400]);
	
	if ok == 0 
		return
	end
	
	% Now compute ZScore for each of the selected files
	hw = waitbar(0,['Computing ZScore for the selected set of ',int2str(length(ResSelect)),' result files...']);
	for k = 1:length(ResSelect)
		dataplot_cb('LoadResultFile',ResFiles(ResSelect(k)).name)
		ResFiles(ResSelect(k)).name
		Zscore_fig = findobj(0,'Tag','tesselation_select');
		hZScore_fig = guihandles(Zscore_fig);
		set(hZScore_fig.ZScore,'Value',1)
		dataplot_cb('ToggleButtonColor',hZScore_fig.ZScore)
		dataplot_cb('ZScore',1) % Force the computation of the ZScore for each of the file (ie force test to 1 in dataplot_cb ZScore)
		waitbar(k/length(ResSelect))
	end
	set(hZScore_fig.Baseline,'Userdata',0) 
	
	delete(hw)
	

	%----------------------------------------------------------------------------
	
	
case 'AnalyzeCorticalMap' % Spatio-temporal segmentation of the cortical current density maps
	
	Zscore_fig = findobj(0,'Tag','tesselation_select');
	hTessSelect = guihandles(Zscore_fig);
	available = findobj(Zscore_fig,'Tag','available');
	removed =  findobj(Zscore_fig,'Tag','removed');
	availableID = get(removed,'Value');
	IDs = get(removed,'String');
	if isempty(IDs), return, end
	
	% Identification of the active mesh surfaces
	Visu = get(DATAPLOT,'Userdata');
	
	try
		load(Visu.Tesselation,'Comment')
	catch
		cd(Users.SUBJECTS)
		load(Visu.Tesselation,'Comment')
	end
	
	if ~iscell(IDs)
		imesh = find(strcmp(IDs,Comment));
	end
	
	nsurf = imesh(availableID); % Selected surface
	
	TessWin  = findobj(0,'tag','tessellation_window'); 
	hTessWin = guihandles(TessWin);
	TessWinUserData = get(TessWin,'Userdata');
	
	if ~isfield(TessWinUserData,'VertConn') % Vertex connectivity not available here - load or compute
		nsurfaces = length(Comment);
		
		% ----------------- Load current surface parameters - 
		load(Visu.Tesselation,'Faces','Vertices')
		FV.faces = Faces{nsurf}; clear Faces
		FV.vertices = Vertices{nsurf}'; clear Vertices
		% ----------------- Load current surface parameters - DONE
		
		load(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile),'VertConn'); % Is the vertex connectivity available ?
		
		if isempty(VertConn) % Vertex connectivity was not defined before
			VertConn = cell(nsurfaces,1);
		else
			if exist(VertConn,'file')
				load(VertConn)
			else
				VertConn = cell(nsurfaces,1);
			end
		end
		
		if isempty(VertConn{nsurf}) % Compute the vertex connectivity
			[pathname, filename] = fileparts(Visu.Tesselation);    
			VertConn{nsurf} = vertices_connectivity(FV);
			clear FV
			VertConnFile = [filename,'_vertconn.mat'];
			
			% Save this file in the current subject folder
			[path,name,ext] = fileparts(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile));
			if ~isempty(path)
				cd(path)
			else
				cd(Users.SUBJECTS)
			end
			
			save(VertConnFile,'VertConn'); 
			
			VertConn = VertConnFile;
			save(Users.CurrentData.SubjectFile,'VertConn','-append') 
		else % Load VC from the file
			load(Users.CurrentData.SubjectFile,'VertConn');
		end
		load(VertConn);
		VertConn = VertConn{nsurf};
		TessWinUserData.VertConn = VertConn; clear VertConn;    
	end
	
	
	% Get Z threshold
	ZThres = str2num(get(hTessSelect.ZScoreThreshold,'String'))  ;
	
	% Vertices which maximum Z score is above ZThres
	data.Results = get(hTessSelect.ResultFiles,'Userdata');
	ImageGridAmp =  data.Results.ImageGridAmp; clear data.Results
	
	% Get baseline info
	data.Results = get(hTessSelect.ResultFiles,'userdata');
	AnalyzeWin = [30 90]/1000; % in sec
	[tmp, tmin] = min(abs(ImageGridTime-AnalyzeWin(1)));
	[tmp, tmax] = min(abs(ImageGridTime-AnalyzeWin(2)));
	clear data.Results
	
	%Baseline = data.Results.ZScore.Baseline; clear data.Results
	%Time = Baseline(end)+1:size(ImageGridAmp,2); %setdiff([1:size(ImageGridAmp,2)],[Baseline(1):Baseline(end)]);
	
	Time = [tmin:tmax];
	
	MaxVerts = max(ImageGridAmp(:,Time)');
	iThres = find(MaxVerts > ZThres);
	% Threshold in terms of correlation
	CorrThres = .9; 
	
	Scout{1} = [];
	iScout = 1;
	imax = 1;
	
	% Sort MaxVerts
	[s, iMaxVertsSort] = sort(MaxVerts); clear s
	
	CONNEX = 1; % if 1; enforce connexity of each clusters 
	
	
	while ~isempty(imax)
		%%%% Segmentation process
		if iScout == 1 % We are starting from scratch
			% Initialization
			% Start region growing from the vertex with maximum ZScore
			imax = iMaxVertsSort(end);
		else
			tmp = setdiff(iMaxVertsSort,[Scout{:}]); 
			[mm, imax] = max(MaxVerts(tmp)); clear mm
			imax = tmp(imax);
		end
		
		% Does imax checks iThres
		if isempty(find(iThres == imax))
			break
		end
		
		
		if CONNEX == 1
			% Look for his neighbors
			voiz = patch_swell(imax,TessWinUserData.VertConn);
			% Select those that check the Zthres condition
			voiz = intersect(iThres,voiz);
			
		else % Just look for all the points above ZThres
			voiz = setdiff(iThres,imax);
		end
		
		% Remove those previously selected
		voiz = setdiff(voiz,[Scout{:}]);
		
		% Compute correlation coefficients of their time series
		%CorrCoeff = abs(corrcoef(ImageGridAmp([imax,voiz],Time)'));
		if ~isempty(voiz)
			CorrCoeff = abs(ImageGridAmp(imax,Time)/norm(ImageGridAmp(imax,Time)))...
				*abs(ImageGridAmp(voiz,Time)'*inorcol(ImageGridAmp(voiz,Time)'));
		else
			CorrCoeff = CorrThres/10;
		end
		
		% Apply Correlation Threshold
		%icorr = find(CorrCoeff(1,2:end)>CorrThres);
		icorr = find(CorrCoeff>CorrThres);
		
		Scout{iScout} = [imax,voiz(icorr)];
		
		while ~isempty(icorr)
			
			if CONNEX == 1
				voiz = patch_swell(voiz(icorr),TessWinUserData.VertConn);
				% Select those that check the Zthres condition
				voiz = intersect(iThres,voiz);
			else
				voiz = setdiff(iThres,imax);
			end
			
			% Remove those previously selected
			voiz = setdiff(voiz,[Scout{:}]);
			
			% Compute correlation coefficients of their time series
			%CorrCoeff = abs(corrcoef(ImageGridAmp([imax,voiz],Time)'));
			if ~isempty(voiz)
				CorrCoeff = abs(ImageGridAmp(imax,Time)/norm(ImageGridAmp(imax,Time)))...
					*abs(ImageGridAmp(voiz,Time)'*inorcol(ImageGridAmp(voiz,Time)'));
				
				CorrCoeffBaseLine = abs(ImageGridAmp(imax,125:500)/norm(ImageGridAmp(imax,125:500)))...
					*abs(ImageGridAmp(voiz,125:500)'*inorcol(ImageGridAmp(voiz,125:500)'));      
			else
				
				CorrCoeff = CorrThres/10;
			end
			
			% Apply Correlation Threshold
			%icorr = find(CorrCoeff(1,2:end)>CorrThres);
			icorr = find(CorrCoeff>CorrThres);
			
			Scout{iScout} = [Scout{iScout},voiz(icorr)];
			
		end
		
		if length([Scout{:}]) ~= length(unique([Scout{:}]))
			keyboard
		end
		
		length(iThres) - length([Scout{:}])
		
		CorticalScouts.CorticalProbePatches(iScout) = Scout(iScout); 
		CorticalScouts.CorticalMarkersLabels{iScout} = int2str(iScout);
		CorticalScouts.CorticalSpots(iScout) = imax;
		if iScout == 1
			load(Visu.Tesselation,'Vertices')
		end
		CorticalScouts.CorticalMarkers(iScout,:) = Vertices{nsurf}(:,Scout{iScout}(1));
		
		iScout = iScout+1;
		bst_message_window('')    
		
		
	end
	
	newScout = 1;
	
	for iScout = 1:length(Scout)
		if length(Scout{iScout}) == 1 % Single point scout: do we keep it ?
			% Check if its neighbors with activation < ZThres have time series that correlate with his with corrcoeff > CorrThres
			voiz = patch_swell(Scout{iScout},TessWinUserData.VertConn);
			
			if ~isempty(voiz)
				CorrCoeff = abs(ImageGridAmp(imax,Time)/norm(ImageGridAmp(imax,Time)))...
					*abs(ImageGridAmp(voiz,Time)'*inorcol(ImageGridAmp(voiz,Time)'));
			else
				CorrCoeff = CorrThres/10;
			end
			
			icorr = find(CorrCoeff>CorrThres);
			
			if isempty(icorr)
				bst_message_window('Remove spurious Scout')
				%               Scout{iScout} = [];
				%               iScout = iScout - 1;
			else
				newCorticalScouts.CorticalProbePatches(newScout) = Scout(iScout); 
				newCorticalScouts.CorticalMarkersLabels{newScout} = int2str(newScout);
				newCorticalScouts.CorticalSpots(newScout) = Scout{iScout}(1);
				if newScout == 1
					load(Visu.Tesselation,'Vertices')
				end
				newCorticalScouts.CorticalMarkers(newScout,:) = Vertices{nsurf}(:,Scout{iScout}(1));
				newScout = newScout + 1;
			end
			
		else
			newCorticalScouts.CorticalProbePatches(newScout) = Scout(iScout); 
			newCorticalScouts.CorticalMarkersLabels{newScout} = int2str(newScout);
			newCorticalScouts.CorticalSpots(newScout) = Scout{iScout}(1);
			if newScout == 1
				load(Visu.Tesselation,'Vertices')
			end
			newCorticalScouts.CorticalMarkers(newScout,:) = Vertices{nsurf}(:,Scout{iScout}(1));
			newScout = newScout + 1;
		end
		
		clear CorticalScouts
		CorticalScouts = newCorticalScouts;
		
	end
	
	
	CorticalScouts.VertConn = TessWinUserData.VertConn;
	CorticalScouts.CorticalProbeDepth = ones(1,length(Scout));
	
	cd(Users.STUDIES)
	
	length(CorticalScouts.CorticalProbePatches)
	
	save test CorticalScouts
	%----------------------------------------------------------------------------
	
	
end
