function varargout = manage_corticalscouts(action,varargin)
%MANAGE_CORTICALSCOUTS - Application M-file for manage_corticalscouts.fig
% function varargout = manage_corticalscouts(action,varargin)

%<autobegin> ---------------------- 26-May-2004 11:30:55 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\brainstorm_filt.m
%   toolbox\bst_color_scheme.m
%   toolbox\bst_layout.m
%   toolbox\bst_message_window.m
%   toolbox\bst_wavedata_display.m
%   toolbox\curvature_cortex.m
%   toolbox\dataplot_cb.m
%   toolbox\findclosest.m
%   toolbox\inorcol.m
%   toolbox\mutincomp.m
%   toolbox\norcol.m
%   toolbox\norlig.m
%   toolbox\patch_swell.m
%   toolbox\smooth_cortex.m
%   toolbox\togglebuttoncolor.m
%   toolbox\tri_area.m
%   toolbox\vertices_connectivity.m
%   toolbox\view_surface.m
%
% Application data and their calls in this file:
%   'FVCData'
%   'Scouts'
%   'TileType'
%   'TriArea'
%   
%   setappdata(ParamFig,'TileType','Tile')
%   setappdata(data.Results.Display.Handles.Figures.TimeSeries,'TileType','T')
%   setappdata(data.Results.Display.Handles.hPatch,'FVCData',...
%   setappdata(data.Results.Display.Handles.hPatch,'Scouts',Scouts)
%   setappdata(data.Results.Display.Handles.hPatch,'Scouts',Scouts);
%   setappdata(data.Results.Display.Handles.hPatch,'TriArea',TriArea)
%   setappdata(fActivation,'TileType','T')
%   
%   FVCData = getappdata(data.Results.Display.Handles.hPatch,'FVCData');
%   Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
%   TriArea = getappdata(data.Results.Display.Handles.hPatch,'TriArea');
%
% Figure Files opened by this function:
%   'manage_corticalscouts.fig'
%
%   Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
%   <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure manage_corticalscouts.fig
%   uicontrol:checkbox:AbsoluteCurrent "Absolute values" uses Callback for
%     manage_corticalscouts mutincomp_meanmax
%   uicontrol:pushbutton:DeleteCorticalSpot "Delete" uses Callback for
%     manage_corticalscouts DeleteCorticalSpot
%   uicontrol:pushbutton:Hpushbutton35 "Activate" uses Callback for
%     manage_corticalscouts CorticalSpotActivity
%   uicontrol:pushbutton:LoadCorticalSpot "Load" uses Callback for
%     manage_corticalscouts LoadCorticalSpot
%   uicontrol:pushbutton:pushbutton38 "+" uses Callback for manage_corticalscouts GrowCorticalArea
%   uicontrol:pushbutton:pushbutton39 "-" uses Callback for manage_corticalscouts ReduceCorticalArea
%   uicontrol:pushbutton:pushbutton41 "MRI" uses Callback for
%     manage_corticalscouts ImportScoutsFromMRI
%   uicontrol:pushbutton:pushbutton42 "Threshold" uses Callback for
%      manage_corticalscouts CorticalScoutFromThres
%   uicontrol:pushbutton:pushbutton43 "Corr. with sensor" uses Callback for
%     manage_corticalscouts CorrelationMapFromSensor
%   uicontrol:pushbutton:Qpushbutton7 "Close" uses Callback for close(gcbf)
%   uicontrol:pushbutton:RenameCorticalSpot "Rename" uses Callback for
%     manage_corticalscouts RenameCorticalSpot
%   uicontrol:pushbutton:SaveCorticalSpot "Save" uses Callback for
%     manage_corticalscouts SaveCorticalSpot
%   uicontrol:pushbutton:SelectCorticalSpot "Probe" uses Callback for
%     manage_corticalscouts SelectCorticalSpot
%   uicontrol:radiobutton:AllCorticalArea "All" uses Callback for
%     manage_corticalscouts mutincomp_meanmax
%   uicontrol:radiobutton:MaxCorticalArea "Max" uses Callback for
%     manage_corticalscouts mutincomp_meanmax
%   uicontrol:radiobutton:MeanCorticalArea "Mean" uses Callback for
%     manage_corticalscouts mutincomp_meanmax
%   uicontrol:togglebutton:pushbutton40 "Show selected" uses Callback for
%     manage_corticalscouts ShowClusters
%
% At Check-in: $Author: Mosher $  $Revision: 20 $  $Date: 5/26/04 9:59a $
%
% 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:30:55 -----------------------

% /---Script Author--------------------------------------\
% |                                                      |
% | *** Sylvain Baillet, Ph.D.                           |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   | 
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% |                                                      |
% \------------------------------------------------------/

DefaultScoutColor = [0 .3 0]; % Default color for selected cortical scout area
DefaultLabelColor = [1 1 1]; % default text color for scout labels
    
global data % Get informations from the main Viewer window

switch(action)
    
case 'create' % Create display-parameter window 
    
    ParamFig = openfig('manage_corticalscouts.fig','reuse');
    setappdata(ParamFig,'TileType','Tile')
    bst_color_scheme(ParamFig);
    bst_layout('align',ParamFig,2,2,4);
    set(ParamFig,'Visible','on')
    
    % Get previous scouts in memory
    if isappdata(data.Results.Display.Handles.hPatch,'Scouts')
        Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
        nScouts = length(Scouts);
        hScoutsManagerFig = guihandles(ParamFig);
        set(hScoutsManagerFig.CorticalSpotList,'String',{Scouts(:).Label},...
            'Max',nScouts,'Value',1)    
    end


    %----------------------------------------------------------------------------
    
case 'SelectCorticalSpot'
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    
    figure(data.Results.Display.Handles.Figures.Surfaces)
    hold on
    
    if nargin == 1
        ginput(1);
    else
        data.Results.Display.Handles.hPatch = (varargin{2});
        figure(data.Results.Display.Handles.Surfaces)
        hold on
    end
    
    scurrent = plot3(0,0,0,'o','parent',gca,'Markersize',5,'Markerfacecolor','r');
    set(scurrent,'visible','off')
    
    vertices = get(data.Results.Display.Handles.hPatch,'vertices');
    % Search for intersection with cortical nodes which 
    % are considered as visible only, ie FVCData ~- NaN
    vert_isnan = get(data.Results.Display.Handles.hPatch,'FaceVertexCdata');
    if size(vert_isnan,2) == 1
        vert_isnan = find(~isnan(vert_isnan));
    else
        vert_isnan = sum(vert_isnan ,2);
        vert_isnan = find(~isnan(vert_isnan));
    end
    
    vertices = vertices(vert_isnan,:);
    
    x = vertices(:,1)';
    y = vertices(:,2)';
    z = vertices(:,3)';		
    clear vertices
    
    if nargin == 1 % Vertex selection from a manual probe
        
        % Get surface vertex which distance to ray entering the viewing box is minimum
        
        C = get(gca,'Currentpoint');		
        
        % Unit vector on traced ray accross viewing box from entering point
        u = (C(1,:)-C(2,:))'/norm(C(1,:)-C(2,:));
        
        inorm=  1:length(x);
        tempo1 = [C(1,1)-x(inorm);C(1,2)-y(inorm);C(1,3)-z(inorm)];
        DIR = u * ones(1,length(inorm));
        tempo = cross(tempo1,DIR);
        clear DIR
        % Compute distance of all vertices to ray
        distdroite = norcol(tempo);
        distboite = norcol(tempo1);
        distdroite = distdroite/max(distdroite);
        distboite = distboite/max(distboite);
        [minn,im] = min(.04*distboite+0.96*distdroite);				
        clear tempo tempo1
        
    else % Point location is passed as an argument 
        inorm=  1:length(x);
        im = varargin{1};
    end
    
    set(scurrent,'Xdata',x(inorm(im)),'Ydata',y(inorm(im)),'Zdata',z(inorm(im))) 
    drawnow
    
    set(scurrent,'visible','on')
    if isempty(im)
        warndlg('Please reselect your location','Location error')	
    end
    
    % Store cortical scout information as appdata of corresponding patch object
    Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
    iScout = length(Scouts) + 1;
    
    Scouts(iScout).Vertex = vert_isnan(inorm(im)); % Vertex index for current scout
    Scouts(iScout).Handle = scurrent; % Handle to plot of scout location
    
    htext = text(1.05*x(inorm(im)),1.05*y(inorm(im)),1.05*z(inorm(im)),int2str(length(Scouts)));
    set(htext,'FontWeight','normal','color','g','Fontname','helvetica','Fontsize',10,'FontUnits','Point')
    Scouts(iScout).HandleLabels = htext; % Handle of text object for scout label 
    Scouts(iScout).Label = get(htext,'String');
    
    setappdata(data.Results.Display.Handles.hPatch,'Scouts',Scouts);
    
    nScouts = length(Scouts);
    
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    set(hScoutsManagerFig.CorticalSpotList,'String',{Scouts(:).Label},...
        'Max',nScouts,'Value',1)    
    
    rotate3d on
    
    
    %----------------------------------------------------------------------------
case 'DeleteCorticalSpot'
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    
    %figure(data.Results.Display.Handles.Figures.Surfaces)
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    SelectedScout = get(hScoutsManagerFig.CorticalSpotList,'Value');
    
    % Update label information
    Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
    RemainingScouts = setdiff([1:length(Scouts)],SelectedScout);
    
    % Delete graphical objects 
    if ~isempty(Scouts)
        delete(Scouts(SelectedScout).Handle)
        delete(Scouts(SelectedScout).HandleLabels)
        Scouts =  Scouts(RemainingScouts);
        set(hScoutsManagerFig.CorticalSpotList,'String',{Scouts(RemainingScouts).Label},'Value',1,'Max',length(Scouts));
    else
        set(hScoutsManagerFig.CorticalSpotList,'String',{''},'Value',1,'Max',1);
    end
    setappdata(data.Results.Display.Handles.hPatch,'Scouts',Scouts);
    
    %----------------------------------------------------------------------------
case 'RenameCorticalSpot'
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    
    figure(data.Results.Display.Handles.Figures.Surfaces)
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    SelectedScout = get(hScoutsManagerFig.CorticalSpotList,'Value');
    SelectedScout = SelectedScout(1); % only one at a time please
    
    Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
    
    Title = '';
    Prompt = sprintf('Please enter a new label for %s',Scouts(SelectedScout).Label);
    tmp = inputdlg(Prompt,Title);
    if isempty(tmp), return, end
    Scouts(SelectedScout).Label = tmp{:};
    set(hScoutsManagerFig.CorticalSpotList,'String',{Scouts(:).Label})
    
    figure(data.Results.Display.Handles.Figures.Surfaces)
    findobj(0,'tag','tessellation_window'); 
    set(Scouts(SelectedScout).HandleLabels,'String',Scouts(SelectedScout).Label)
    setappdata(data.Results.Display.Handles.hPatch,'Scouts',Scouts);
    
    %----------------------------------------------------------------------------
    
case {'GrowCorticalArea','ReduceCorticalArea'} % Grow/shrink selected scout area 
    % by including connected successive rings of connected vertices
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    if ~isfield(data.Results.Display.OPTIONS,'VertConn') % Vertex connectivity not available here - load or compute
        
        
        % Is vertex connectivity available ?
        load(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile),'VertConn'); 
        if ~exist('VertConn','var') % Vertex connectivity not computed yet
            VertConn = '';
        end
        
        VertConnFile = VertConn; % Filename containing connectivity information
        [pathname, filename] = fileparts(fullfile(Users.SUBJECTS,data.Results.OPTIONS.GridLoc));    
        
        cd(pathname)
        if ~exist(VertConnFile,'file')
            
            % How many surfaces in tessellation file ?
            % should be only one but older BrainStorm files may contain several
            load(fullfile(Users.SUBJECTS,data.Results.OPTIONS.GridLoc),'Comment');
            nsurfaces = length(Comment); % Number of surfaces in tessellation file
            VertConn = cell(nsurfaces,1); % Create cell array where to store vertex connectivity values
            
            VertConn{data.Results.OPTIONS.iGrid(1)} = vertices_connectivity(fv); % compute vertex connectivity 
            [pathname, filename] = fileparts(fullfile(Users.SUBJECTS,data.Results.OPTIONS.GridLoc));    
            cd(pathname)
            if ~exist(VertConnFile,'file')
                VertConnFile = [filename,'_vertconn.mat']; % Where to save it
            end
            
            % Save 
            save(fullfile(pathname,VertConnFile),'VertConn'); 
            % Update tessellation file information
            tmp = VertConn; % CBB VertConn should be stored tessellation file not in separate file (SB, 09-Oct-2003)
            VertConn = VertConnFile;
            save(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile),'VertConn','-append'); 
            data.Results.Display.OPTIONS.VertConn = tmp; clear tmp nsurfaces
            
        else % load existing vertex connectivity map
            
            load(fullfile(pathname,VertConnFile))
            data.Results.Display.OPTIONS.VertConn = VertConn; clear VertConn
            
        end
        
    end
    
    % Now grow a patch around the selected probe by adding a ring neighbors
    
    SelectedScout = get(hScoutsManagerFig.CorticalSpotList,'Value');
    Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
    
    
    switch(action)
    case 'GrowCorticalArea'
        Scouts(SelectedScout).Vertex = [Scouts(SelectedScout).Vertex,...
                patch_swell([Scouts(SelectedScout).Vertex],data.Results.Display.OPTIONS.VertConn)];
    case 'ReduceCorticalArea'
        % Reduce more gradually the area of the cortical patch to be monitored
        if length(Scouts(SelectedScout).Vertex)>1
            Scouts(SelectedScout).Vertex = [Scouts(SelectedScout).Vertex(1:end-1)];
        end
    end
    
    % Now show clusters
    % Keep track of the original activation vertex colordata
    if ~isappdata(data.Results.Display.Handles.hPatch,'FVCData')
        setappdata(data.Results.Display.Handles.hPatch,'FVCData',...
            get(data.Results.Display.Handles.hPatch,'FaceVertexCData')); 
    end
    FVCData = getappdata(data.Results.Display.Handles.hPatch,'FVCData');
    
    % RGB Values
    Scouts(SelectedScout).Color = DefaultScoutColor; % CBB | User may assign/alter scout patch color ?
    FVCData(Scouts(SelectedScout).Vertex,:) = repmat(Scouts(SelectedScout).Color,length(Scouts(SelectedScout).Vertex),1);
    
    set(data.Results.Display.Handles.hPatch,'FaceVertexCData',FVCData);
    
    % Compute corresponding surface area
    if ~isappdata(data.Results.Display.Handles.hPatch,'TriArea')
        TriArea = 1e4*tri_area(get(data.Results.Display.Handles.hPatch,'Faces'), get(data.Results.Display.Handles.hPatch,'Vertices'));
        setappdata(data.Results.Display.Handles.hPatch,'TriArea',TriArea)
    else
        TriArea = getappdata(data.Results.Display.Handles.hPatch,'TriArea');
    end
    
    % CBB make the following area computation faster... 
    kk = 0;
    for k = Scouts(SelectedScout).Vertex
        kk = kk+1;
        [i,j] = find(get(data.Results.Display.Handles.hPatch,'Faces')-k == 0);
        tri{kk}= i';
    end
    if isempty(tri{1})
        tria = [];
    else
        tria = unique([tri{:}]);
    end
    o_area = sum(TriArea(tria))/3; %in cm2
    
    set(hScoutsManagerFig.ScoutArea,'String',sprintf('%3.1f cm2',o_area))
    
    setappdata(data.Results.Display.Handles.hPatch,'Scouts',Scouts)    
    
    %----------------------------------------------------------------------------
    
case 'ShowClusters' % Show all cortical scouts/clusters at once
    
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    togglebuttoncolor(gcbo)
    
    Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
    
    % Which clusters were selected ?
    SelectedScout = get(hScoutsManagerFig.CorticalSpotList,'Value');
    
    % Original current mapping
    FVCData = getappdata(data.Results.Display.Handles.hPatch,'FVCData');
    
    if get(gcbo,'Value') % Show selected clusters
        FVCData([Scouts(SelectedScout).Vertex],:) = repmat(DefaultScoutColor,length([Scouts(SelectedScout).Vertex]),1);
    end
    
    set(data.Results.Display.Handles.hPatch,'FaceVertexCData',FVCData);
    
    %----------------------------------------------------------------------------
    
case 'SaveCorticalSpot' % Save cortical scout locations 
    % Cortical scouts need to be indexed with both the corresponding cortical surface name 
    % and the current data set name 
    
    data.Results.Display.Handles.hPatch  = findobj(0,'tag','tessellation_window'); 
    CorticalScouts = get(data.Results.Display.Handles.hPatch,'Userdata');
    
    CorticalScouts.HandleLabels = get(CorticalScouts.HandleLabels,'String');
    CorticalScouts.Handle = [get(CorticalScouts.Handle,'Xdata');...
            get(CorticalScouts.Handle,'Ydata');...
            get(CorticalScouts.Handle,'Zdata')];
    if iscell( CorticalScouts.Handle)
        CorticalScouts.Handle = [CorticalScouts.Handle{:}];
    end
    
    CorticalScouts.Handle = reshape(CorticalScouts.Handle,length(CorticalScouts.HandleLabels),3);
    % CorticalScouts.Handle is now Number of Scouts x 3
    
    ScoutsManagerFig = findobj(0,'Tag','tesselation_select');
    hSelect = guihandles(ScoutsManagerFig);
    
    % What's the name of the current cortical envelope ?
    ActiveTess = get(hSelect.removed,'String'); % Find the active Cortical surface
    iCortex = get(hSelect.removed,'Value'); 
    if iscell(ActiveTess)
        ActiveTess= ActiveTess{iCortex};
    end
    
    if isfield(CorticalScouts,'CorticalProbePatches') % Process members of cortical patches to get their x,y,z's
        Visu = get(DATAPLOT,'Userdata');
        try
            load(Visu.Tesselation,'Comment','Vertices')
        catch
            cd(Users.SUBJECTS)
            load(Visu.Tesselation,'Comment','Vertices')
        end
        % What's the current surface ?
        imesh = find(strcmp(Comment,ActiveTess));
        Vertices = Vertices{imesh};
        CorticalScouts.CorticalProbePatchesXYZ = cell(size(CorticalScouts.CorticalProbePatches));
        for k=1:length(CorticalScouts.CorticalProbePatches)   % for each area
            CorticalScouts.CorticalProbePatchesXYZ{k} = ...
                Vertices(:,CorticalScouts.CorticalProbePatches{k});
        end
        
    end
    
    % What's the name of the current data set ?
    ResultFiles = get(hSelect.ResultFiles,'String');
    ResultFile = ResultFiles{get(hSelect.ResultFiles,'Value')};
    
    cd(Users.STUDIES)
    dirr = fileparts(Users.CurrentData.StudyFile);
    cd(dirr);
    
    [filename, pathname] = uiputfile('CorticalScout.mat', 'Save Cortical Scout in...');
    save(fullfile(pathname,filename), 'ResultFile','ActiveTess','CorticalScouts','-mat');
    
    %----------------------------------------------------------------------------
    
case 'LoadCorticalSpot' % Load cortical scout locations 
    data.Results.Display.Handles.hPatch  = findobj(0,'tag','tessellation_window'); 
    
    if nargin == 1
        
        CorticalScouts = get(data.Results.Display.Handles.hPatch,'Userdata');
        if ~isempty(CorticalScouts) % some exist - remove
            if isfield(CorticalScouts,'CorticalMarkers')
                delete(CorticalScouts.Handle)
                delete(CorticalScouts.HandleLabels)
                clear CorticalScouts
            end
        end
        
        cd(Users.STUDIES)
        dirr = fileparts(Users.CurrentData.StudyFile);
        cd(dirr);
        
        [filename, pathname] = uigetfile('CorticalScout.mat', 'Save Cortical Scout in...');
        if filename == 0
            return
        end
        load(fullfile(pathname,filename), 'ResultFile','ActiveTess','CorticalScouts','-mat');
    else
        load(varargin{1}, 'ResultFile','ActiveTess','CorticalScouts','-mat');
    end
    
    figure(data.Results.Display.Handles.hPatch), hold on
    MainAxes = findobj(data.Results.Display.Handles.hPatch,'Tag','MainAxes');
    vertices = get(findobj(MainAxes,'type','patch'),'vertices');
    x = vertices(:,1)';
    y = vertices(:,2)';
    z = vertices(:,3)';		
    clear vertices
    
    CorticalScouts.Handle = zeros(length(CorticalScouts.Vertex),1);
    for k = 1:length(CorticalScouts.Vertex)
        scurrent = plot3(0,0,0,'o','parent',MainAxes,'Markersize',5,'Markerfacecolor','r');
        set(scurrent,'visible','off')
        set(scurrent,'Xdata',x(CorticalScouts.Vertex(k)),...
            'Ydata',y(CorticalScouts.Vertex(k)),...
            'Zdata', z(CorticalScouts.Vertex(k)))
        set(scurrent,'visible','on')
        CorticalScouts.Handle(k) = scurrent;
        
        r_text =1.1;
        
        if iscell(CorticalScouts.HandleLabels)
            htext = text(r_text*x(CorticalScouts.Vertex(k)),...
                r_text*y(CorticalScouts.Vertex(k)),...
                r_text*z(CorticalScouts.Vertex(k)),...
                CorticalScouts.HandleLabels{k});
        else
            htext = text(r_text*x(CorticalScouts.Vertex(k)),...
                r_text*y(CorticalScouts.Vertex(k)),...
                r_text*z(CorticalScouts.Vertex(k)),...
                CorticalScouts.HandleLabels);
            
        end
        
        
        set(htext,'FontWeight','normal','color','g','Fontname','helvetica',...
            'Fontsize',10,'FontUnits','Point')
        
        Label_tmp(k) = htext;
    end
    
    ScoutsManagerFig = findobj(0,'Tag','tesselation_select');
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    set(hScoutsManagerFig.CorticalSpotList,'String',CorticalScouts.HandleLabels,...
        'Value',1,'Max', length(CorticalScouts.HandleLabels))
    
    CorticalScouts.HandleLabels = Label_tmp; clear Label_tmp;
    
    set(data.Results.Display.Handles.hPatch,'Userdata',CorticalScouts);
    
    figure(data.Results.Display.Handles.hPatch)
    rotate3d on
    
    %----------------------------------------------------------------------------
    
case 'mutincomp_meanmax' % Mutually icompatible checkboxes
    handles = guihandles(gcbf);
    
    mutincomp([handles.MeanCorticalArea,handles.MaxCorticalArea,handles.AllCorticalArea])
    h = findobj([handles.MeanCorticalArea,handles.MaxCorticalArea,handles.AllCorticalArea],'Value',1); 
    if length(h)>1
        set(h(2:end),'Value',0) 
    end   
    
    %----------------------------------------------------------------------------
    
case 'CorticalSpotActivity'
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    % Get scout information
    SelectedScout = get(hScoutsManagerFig.CorticalSpotList,'Value'); % Which scout was selected ?
    Scouts = getappdata(data.Results.Display.Handles.hPatch,'Scouts');
    
    % What kind of display is wanted here ?
    MEAN = get(hScoutsManagerFig.MeanCorticalArea,'Value');
    MAX = get(hScoutsManagerFig.MaxCorticalArea,'Value');
    ALL = get(hScoutsManagerFig.AllCorticalArea,'Value');
    
    % -> Prepare cell array containing time series to display
    for k = 1:length(SelectedScout) % For each cortical survey area
        if MEAN == 1% Average of the patch activity at each time instant
            if get(hScoutsManagerFig.AbsoluteCurrent,'Value')
                CorticalActivity{k} = mean(abs(double(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:))),1);
            else
                CorticalActivity{k} = sign(mean(double(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:))))...
                    .* mean(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:),1) ;
            end
        elseif MAX==1 % Strongest
            if get(hScoutsManagerFig.AbsoluteCurrent,'Value')
                CorticalActivity{k} = max(abs(double(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:))),[],1);
            else
                CorticalActivity{k} = max(double(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:)),[],1);
            end
            
        else % Display every source activity in the patch
            % All sources need to be displayed
            if get(hScoutsManagerFig.AbsoluteCurrent,'Value')
                CorticalActivity{k} = abs(double(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:)));
            else
                CorticalActivity{k} = double(data.Results.ImageGridAmp(Scouts(SelectedScout(k)).Vertex,:));
            end
            % Detect Maxima in amplitudes (needed for proper display in the offset plot)
            maxCorticalActivity{k} = sign(min(CorticalActivity{k})).*max(abs(CorticalActivity{k}));
        end
    end
    % <- Prepare cell array containing time series to display
    
    % Now display these time series
    if ~isfield(data.Results.Display.Handles.Figures,'TimeSeries')
        fActivation = figure;
        data.Results.Display.Handles.Figures.TimeSeries = fActivation;
        setappdata(fActivation,'TileType','T')
        bst_color_scheme(fActivation)
        bst_layout('align',fActivation,2, 2, 3)
        bst_layout('align',data.Results.Display.Handles.Figures.Surfaces,2, 2, 1)
    else
        figure(data.Results.Display.Handles.Figures.TimeSeries)
        clf
        setappdata(data.Results.Display.Handles.Figures.TimeSeries,'TileType','T')
        bst_color_scheme(data.Results.Display.Handles.Figures.TimeSeries)
        bst_layout('align',data.Results.Display.Handles.Figures.TimeSeries,2, 2, 3)
        bst_layout('align',data.Results.Display.Handles.Figures.Surfaces,2, 2, 1)
    
    end
    
    % Display set-up
    tail = {'MEAN','MAX','ALL'};
    set(gcf,'Tag','source_time_series','Name',['Activation time series: ', tail{find([MEAN,MAX,ALL])}])
    Color = get(gca,'ColorOrder'); % Default Colormapping for line plots
    
    ColumnLayout = 1; % Flag
    OPTIONS.nColumnLayout = 1; % Number of columns
    OverlayAllChannels = 0;
    nDispTypes = 1;
    nSets = length(SelectedScout);
    OPTIONS.Time = linspace(data.Display.Time(1),data.Display.Time(end),data.Display.NumberOfSamples);
    OPTIONS.CurrentTime = dataplot_cb('GET','current_time');
    
    % Initialize column plot structure
    tmp = deal(struct('lines',{cell(1,nSets)},'axes',zeros(OPTIONS.nColumnLayout,1),'ytick',[]));
    PlotHandles.PlotColumn = tmp;
    clear tmp
    
    % Display parameter for proper layout of time series
    OPTIONS.Fmax = [];
    for iSet = 1:nSets
        OPTIONS.Fmax(iSet) = max(max(abs(CorticalActivity{iSet})));
    end
    
    OPTIONS.WaveLabel = {Scouts.Label};
    OPTIONS.DataType = 'source';
    OPTIONS.Modality = 1;
    
    OPTIONS.LineColor = num2cell(rand(length(CorticalActivity),3),2);

    switch(get(hScoutsManagerFig.popupScoutPlotType,'Value'))
    case 1 % Overlay all timeseries
    case 2 % Column plot
        OPTIONS.DisplayTypes = [1 0 0];
        if 1,%size(cat(1,CorticalActivity{:}),1)/length(CorticalActivity) > 1 % Several sources per scouts
            OPTIONS.Modality = 1:length(CorticalActivity);
            OPTIONS.SelectedChannels{1} = 1:size(CorticalActivity{1},1);
            for k = 2:length(CorticalActivity)
                OPTIONS.SelectedChannels{k} = size(CorticalActivity{k-1},1)+1:size(CorticalActivity{k-1},1)+size(CorticalActivity{k},1);
            end
        end
        CorticalActivity = {cat(1,CorticalActivity{:})};    
    end
    OPTIONS.DisplayTypes = [1 0 0];
    [data.Results.Display.Scouts.Handles,OPTIONS] = ...
        bst_wavedata_display(CorticalActivity,OPTIONS);

    
%----------------------------------------------------------------------------

case 'ImportScoutsFromMRI' % Create a scout file from points selected in the MR
    
    if ~ishandle(data.Results.Display.Handles.Figures.Surfaces)
        return
    end
    ScoutsManagerFig = gcbf;
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    cd(Users.STUDIES);
    [filename, pathname] = uigetfile('*.txt', 'Select an MRI Point File');
    if filename == 0, return, end
    CorticalScouts = struct('MarkersLabels','','CorticalSpots',[],'CorticalMarkers',[]);
    
    [Number,Labels,X,Y,Z,Xmri,Ymri,Zmri] = ...
        textread(fullfile(pathname,filename),'%d %s %f %f %f %f %f %f',-1);
    clear Number Xmri Ymri Zmri
    
    Vertices = get(data.Results.Display.Handles.hPatch,'Vertices');
    X = X/1000;
    Y = Y/1000;
    Z = Z/1000;
    
    for k = 1:length(X) % Find closest vertex in the tessellation
        [mm Scouts.Vertex(k)] = min(norcol(...
            [Vertices(1,:)-X(k); ...
                Vertices(2,:)-Y(k); ...
                Vertices(3,:)- Z(k)])); 
        %Scouts.Handle(k,:) = Vertices(:,Scouts.Vertex(k))';
    end
    %Scouts.Handle = Vertices(:,Scouts.Vertex)';
    clear Vertices
    
    % Create Cortical Scout File
    cd(Users.STUDIES);
    [filename, pathname] = uiputfile(['CorticalScoutFromMRI.mat'], 'Save Cortical Scout in...');
    save(fullfile(pathname,filename), 'ResultFile','ActiveTess','CorticalScouts','-mat');
    
    % Visualization
    dataplot_cb('LoadCorticalSpot',fullfile(pathname,filename))
    
    %----------------------------------------------------------------------------
    
case 'CorrelationWithDepthElectrodes' % Compute time-series correlation between recorded signals from
    % depth electrodes and corresponding locations in imaging map
    
    data.Results.Display.Handles.hPatch  = findobj(0,'tag','tessellation_window'); 
    %CorticalScouts = get(data.Results.Display.Handles.hPatch,'Userdata');
    ScoutsManagerFig = findobj(0,'Tag','tesselation_select');
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    
    % Load specific cortical scout file (depth electrodes)
    %DepthScouts = load('I:\Data\Studies\Epilepsie\02May02\DepthElectrodes_Scouts.mat');
    DepthScouts = load('C:\Users\Silvin_E\BrainStorm2001\Studies\Epilepsie\02May02\DepthElectrodes_Scouts.mat');
    Visu = get(DATAPLOT,'Userdata');
    % Identification of the active mesh surfaces
    try
        load(Visu.Tesselation,'Comment')
    catch
        cd(Users.SUBJECTS)
        load(Visu.Tesselation,'Comment')
    end
    available = findobj(ScoutsManagerFig,'Tag','available');
    removed =  findobj(ScoutsManagerFig,'Tag','removed');
    availableID = get(removed,'Value');
    IDs = get(removed,'String');
    if isempty(IDs), return, end
    
    if ~iscell(IDs)
        imesh = find(strcmp(IDs,Comment));
    end
    
    nsurf = imesh(availableID); % Selected surface
    load(Visu.Tesselation,'Faces','Vertices')
    Faces = Faces{nsurf};
    Vertices = Vertices{nsurf};
    % Compute tri area
    A = tri_area(Faces,100*Vertices');
    
    
    % First ask 
    %  on which time window to compute correlation
    Results = get(hScoutsManagerFig.ResultFiles,'Userdata');
    Prompt = {...
            'Time window for correlation computation (in ms):',...
        };
    
    Answer = inputdlg(Prompt,'',1,{sprintf('%3.2f %3.2f',Results.ImageGridTime(1)*1000, Results.ImageGridTime(end)*1000')});
    if isempty(Answer{1})
        return
    end
    
    CorrTimeWin = str2num(Answer{1})/1000; % Time window for correlation computation
    clear Answer
    
    
    % Now for each depth electrode - compute correlation between recorded signal and estimated source time courses 
    for k = 1:length(DepthScouts.CorticalScouts.Vertex)
        % Find sensor in the data's Channel structure
        DepthScouts.CorticalScouts.HandleLabels{k}
        
        CorrRefChan = find(strcmpi(DepthScouts.CorticalScouts.HandleLabels{k},{Visu.Channel.Name}));
        if isempty(CorrRefChan), return, end % Channel name not found
        CorrRefChanName = Visu.Channel(CorrRefChan).Name;
        F = Visu.Data.F(CorrRefChan,:); 
        % Find time indices for correlation computation
        % Trim out source and sensor time series
        Results = get(hScoutsManagerFig.ResultFiles,'Userdata');
        F = F(Results.Time);
        TimeInd = findclosest(CorrTimeWin,Results.ImageGridTime');
        
        load(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile),'VertConn'); % Is the vertex connectivity available ?
        load(fullfile(Users.SUBJECTS,fileparts(Users.CurrentData.SubjectFile),VertConn))
        % Now grow a patch around the selected probe and selected neighbors
        patchh = ...
            [DepthScouts.CorticalScouts.Vertex(k),patch_swell(DepthScouts.CorticalScouts.Vertex(k),VertConn{nsurf})];
        patchh = [patchh,patch_swell(patchh,VertConn{nsurf})];
        %patchh = [patchh,patch_swell(patchh,VertConn{nsurf})];
        %patchh = [patchh,patch_swell(patchh,VertConn{nsurf})];
        kkk = 0;
        for kk=patchh % For each element of the patch
            kkk = kkk+1;
            [i,j] = find(Faces-kk == 0); % Find all triangles for which current parch element is a vertex
            if ~isempty(i)
                %tri{kkk} = i';
                VertArea(kkk)= sum(A(i))/3; % Each vertex contributes to 1/3 of the each triangle surface it belongs to
            else
                kkk = kkk-1;
            end
        end
        
        %tri = unique([tri{:}]); % 
        %area(k) = sum(A(tri))/3;% in cm2
        
        area(k) = sum(VertArea); % Total area covered by path k
        
        clear tri kk kkk VertArea
        % take average activity within patch
        %         [maxx ind] = max(abs(Results.ImageGridAmp(patch,TimeInd(1):TimeInd(end))));
        %         Time = TimeInd(1):TimeInd(end);
        %         for time = 1:length(Time)
        %             tmp(time) = maxx(time).*sign(Results.ImageGridAmp(patch(ind(time)),Time(time)));
        %         end
        %         Results.ImageGridAmp = tmp;
        
        %Results.ImageGridAmp = mean((Results.ImageGridAmp(patch,TimeInd(1):TimeInd(end))));
        %Results.ImageGridAmp = Results.ImageGridAmp - mean(Results.ImageGridAmp);
        
        [Results.ImageGridAmp,filtwts] = brainstorm_filt(Results.ImageGridAmp(patchh,:),1250,0,70);
        [F,filtwts] = brainstorm_filt(F,1250,0,70);
        Results.ImageGridAmp = Results.ImageGridAmp(:,TimeInd(1):TimeInd(end));
        Results.ImageGridAmp = Results.ImageGridAmp - repmat(mean(Results.ImageGridAmp,2),1,length(TimeInd(1):TimeInd(end)));
        
        
        F = F(TimeInd(1):TimeInd(end));
        F = (F - mean(F));
        Results.ImageGridAmp = (double(Results.ImageGridAmp));
        
        %Compute correlation coefficients for all sources
        
        [CorrCoeff(k) imax] = max(abs((F/norm(F))...
            *Results.ImageGridAmp'*inorcol(Results.ImageGridAmp')));
        
        SignCorr(k) = sign(((F/norm(F))...
            *Results.ImageGridAmp(imax,:)'*inorcol(Results.ImageGridAmp(imax,:)')));
        
        Source(k,:) = Results.ImageGridAmp(imax,:);
        Data(k,:) = F;
        
        %end
        
        % Plot
        %for k = 1:length(DepthScouts.CorticalScouts.Vertex)
        if abs(CorrCoeff(k)) > .4
            hfig = figure;
            TimeVec = 1000*linspace(CorrTimeWin(1),CorrTimeWin(2), length(F));
            ploot = 100*Data(k,:)/(norm(Data,'fro')/sqrt(length(DepthScouts.CorticalScouts.Vertex))); % Normalize amplitudes by average signal energy over all electrodes
            plot(TimeVec,ploot,'r','linewidth',2), hold on
            if SignCorr(k) < 0
                plot(TimeVec,-Source(k,:)*max(abs(ploot))/max(abs(Source(k,:))),'b','linewidth',2);
            else
                plot(TimeVec, Source(k,:)*max(abs(ploot))/max(abs(Source(k,:))),'b','linewidth',2);
            end
            set(gcf,'color','k')
            set(gca,'xcolor','w','ycolor','w')
            ylabel('% average total signal energy')
            title(sprintf('Channel: %s - Correlation: %3.2f %% - Area Explored: %3.2f cm2',DepthScouts.CorticalScouts.HandleLabels{k},100*abs(CorrCoeff(k)),area(k)),'color','w')
            set(hfig,'Name',CorrRefChanName)
        end
        %end
        ans = input('Plot cortical maps (y/n) ? ','s');
        switch(ans)
        case 'y'
            
            imax = patchh(imax);
            patchh = ...
                [DepthScouts.CorticalScouts.Vertex(9),patch_swell(DepthScouts.CorticalScouts.Vertex(9),VertConn{nsurf})];
            patchh = [patchh,patch_swell(patchh,VertConn{nsurf})];
            for kkkk =1 : 35
                patchh = [patchh,patch_swell(patchh,VertConn{nsurf})];
            end
            Results = get(hScoutsManagerFig.ResultFiles,'Userdata');
            tmp = Results.ImageGridAmp(patchh,TimeInd(1):TimeInd(end));
            Results.ImageGridAmp = NaN * Results.ImageGridAmp(:,TimeInd(1):TimeInd(end));
            Results.ImageGridAmp(patchh,:) = tmp; clear tmp
            
            
            % Surface smoothing
            fv.vertices = Vertices';
            fv.faces = Faces;
            load C:\Users\Silvin_E\BrainStorm2001\Subjects\epilepsie\02May02\02May02_tess_vertconn
            VertConn = VertConn{6};
            [smoothedFV,A]=smooth_cortex(fv,VertConn,.3,100);
            
            %tmpVertices = NaN * Vertices;
            %tmpVertices(:,patchh) = Vertices(:,patchh);
            tmpVertices = NaN * smoothedFV.vertices';
            tmpVertices(:,patchh) = smoothedFV.vertices(patchh,:)';
            [curvature_sigmoid,curvature] = curvature_cortex(fv,VertConn,.1,0);
            curvature(curvature>=0) = .8;
            curvature(curvature<0) = -.8;
            
            data.curDensityMapArray{1} = [zeros(17,2) (0.5:0.03:1)'
                zeros(29,1) (0.03:0.034:1)' ones(29,1)
                (0.3:0.3:1)' ones(3,2)
                1 1 1
                ones(26,1) (1:-0.039:0)' (1:-0.039:0)'
                (1:-0.021:0.5)' zeros(24,2)];
            
            data.curDensityMapArray{2} = [zeros(17,2) (0.5:0.03:1)' %17
                zeros(29,1) (0.03:0.034:1)' ones(29,1) %29
                (0.3:0.3:1)' ones(3,2) %3
                1 1 1
                ones(50,2) (1:-0.0202:0)']; %50
            
            tmp = hot(140);%flipdim(hot(100),1);
            data.curDensityMapArray{3} = tmp(1:100,:);%flipdim(hot(100),1);
            
            %             [TH,PHI,R] = cart2sph(tmpVertices(1,patchh),tmpVertices(2,patchh),tmpVertices(3,patchh));
            %             PHI2 = zeros(size(PHI));
            %             R2 = R./cos(PHI).^.2;
            %             [tmpVertices(2,patchh),tmpVertices(1,patchh)] = pol2cart(TH,R2);
            %             tmpVertices(3,patchh) = 0;
            %             
            
            [hf,hs,hl] = view_surface('test',Faces,tmpVertices,Results.ImageGridAmp(:,1));
            %             hold on
            subplot(221)
            [hf,hs,hl] = view_surface('test',Faces,tmpVertices,Results.ImageGridAmp(:,1));
            delete(findobj(gcf,'type','light'))
            view(0,0), material dull, camlight, 
            
            hold on
            hax = subplot(222);
            hs(2) = copyobj(hs,hax);
            axis image
            axis off
            view(-180,0), material dull, camlight
            
            
            set(hs,'backfacelighting','unlit')
            hold on
            iScout = DepthScouts.CorticalScouts.Vertex(k); 
            ind = k;
            hsaxes = get(hs(1),'Parent'); 
            axes(hsaxes)
            helec = plot3(tmpVertices(1,iScout),tmpVertices(2,iScout),tmpVertices(3,iScout),'ow','MarkerFaceColor','none','MarkerSize',10);
            hbest = plot3(tmpVertices(1,(imax)),tmpVertices(2,(imax)),tmpVertices(3,(imax)),'ow','MarkerFaceColor','none','MarkerSize',5);
            hsaxes = get(hs(2),'Parent'); 
            axes(hsaxes)
            helec(2) = plot3(tmpVertices(1,iScout),tmpVertices(2,iScout),tmpVertices(3,iScout),'ow','MarkerFaceColor','none','MarkerSize',10);
            hbest(2) = plot3(tmpVertices(1,(imax)),tmpVertices(2,(imax)),tmpVertices(3,(imax)),'ow','MarkerFaceColor','none','MarkerSize',5);
            set(hf,'color','k')
            %colorbar
            
            subplot(212)
            Time = 1000*Results.ImageGridTime(TimeInd(1):TimeInd(end));
            hp = plot(Time,ploot,'w','linewidth',2);
            %set(hp,'visible','off')
            
            hold on
            
            %Source(ind,:) = norcol(F);
            
            if SignCorr(ind) < 0
                plot(TimeVec,-Source(ind,:)*max(abs(ploot))/max(abs(Source(ind,:))),'w','linewidth',1);
            else
                plot(TimeVec, Source(ind,:)*max(abs(ploot))/max(abs(Source(ind,:))),'w','linewidth',1);
            end
            
            %ploot = Source(ind,:)*max(abs(ploot))/max(abs(Source(ind,:)));
            hcurr = plot(Time(1),ploot(1),'or','linewidth',2);
            
            set(gca,'xcolor','w','ycolor','w','color','k')
            xlabel('Time (ms)')
            hcurraxes = get(hcurr,'Parent');
            
            colormap flag; colormenu
            flag = get(gcf,'colormap'); %colormap([0 0 1; repmat([1 0 0],size(flag,1)-1,1)])
            flag = flag(:,[3 2 1]);
            colormap(flag)
            
            pwd
            %hsaxes = get(hs(1),'Parent'); 
            %set(hsaxes,'clim',[.98 1])
            %set(hsaxes,'clim',[0 max(abs(Results.ImageGridAmp(:)))])
            %hsaxes = get(hs(2),'Parent'); 
            %set(hsaxes,'clim',[.98 1])
            %set(hsaxes,'clim',[0 max(abs(Results.ImageGridAmp(:)))])
            
            % ------------------------------------------------------------------------------
            cd_threshold = .2;
            cd_transparency = 0.2;
            anatomy = 0;
            viewmode = 1; % Real data or 0 = abs values
            colormap_selected = 3;
            movie = 0;
            activations = 0;
            if movie 
                Avifile = uiputfile('*.avi');%TObas4_corr.avi';%[DepthScouts.CorticalScouts.HandleLabels{ind},'_8.avi'];
                mov = avifile(Avifile,'FPS',2,'QUALITY',100,'compression','indeo5');   
            end
            corrsize = 30;%50
            % ------------------------------------------------------------------------------
            
            if ~activations
                timee = corrsize/2+1:2:size(Results.ImageGridAmp,2)-corrsize/2;
            else
                timee = 1:2:size(Results.ImageGridAmp,2);%
            end
            
            corr = zeros(size(Results.ImageGridAmp(patchh,:),1),length(timee));
            tt = 0;
            hw = waitbar(0,'Computing correlations');
            cc = zeros(length(timee),1);
            for t = timee
                tt = tt+1;
                if ~rem(tt,10)
                    waitbar(tt/length(timee))
                end
                
                if ~activations
                    activ = Results.ImageGridAmp(patchh,t-corrsize/2:t+corrsize/2);
                    elec = F(t-corrsize/2:t+corrsize/2);
                    %                     corr(:,tt) = abs(((elec/norm(elec))...
                    %                         *activ'*inorcol(activ')))';
                    %                     
                    % Compute corrcoef over time and accross cortex (patchh)
                    activ = (activ)/max(abs(activ(:)));
                    energy = norlig(activ);
                    [senergy iener] = sort(energy);
                    activ = activ(iener(end-100+1:end),:)';
                    tmp = corrcoef(activ);
                    cc(tt) = mean(abs(tmp(:)));
                    
                    
                else
                    if t == timee(1)
                        tmp = Results.ImageGridAmp(patchh,timee);
                        corr = abs(tmp)/max(tmp(:));%abs(activ(:,1+corrsize/2));
                    end    
                end
                
            end
            close(hw)
            corr(corr<cd_threshold) = 0; %Show current density that is at least 'display_threshold'% the maximum absolute value
            limitValue = max(abs(corr(:)));
            
            
            %Apply alpha blending
            %create anatomy and current density colormaps
            if(anatomy==1) anatomy=0.99; end
            scale = linspace(.4,.6,100);%anatomy:(1-anatomy)/99:1;
            data.anatMap = repmat(scale',1,3);
            if(viewmode == 1) %if real values
                data.curDensityMap = data.curDensityMapArray{colormap_selected};
            else %if abs values
                data.curDensityMap = data.curDensityMapArray{colormap_selected+2};
            end
            colormap(data.curDensityMap);
            
            %Get RGB info for anatomy
            min_curv = min(curvature);
            max_curv = max(curvature);
            data.index_anatomy = (floor((curvature-min_curv)*99/(max_curv-min_curv)+1));
            data.anatRGB = data.anatMap(data.index_anatomy,:);
            
            for t = 1:length(timee)
                
                %                 tmp = NaN * Results.ImageGridAmp(:,t);
                %                 tmp(patchh) = corr(:,t);
                %                 
                if(viewmode == 1)%if real values
                    index_density = floor((corr(:,t)+limitValue)*99/(2*limitValue)+1);
                    data.curDensityMap = data.curDensityMapArray{colormap_selected};
                else %if abs values
                    corr = abs(corr);
                    index_density = floor((corr(:,t)-limitValue)*99/limitValue+100);
                    data.curDensityMap = data.curDensityMapArray{colormap_selected+2};
                end
                data.curDensityRGB = data.curDensityMap(index_density,:);
                data.toBlend = find(corr(:,t) ~= 0);
                mixedRGB = data.anatRGB;
                toBlend = data.toBlend;
                mixedRGB(patchh(toBlend),:) = cd_transparency*data.anatRGB(toBlend,:)+(1-cd_transparency)*data.curDensityRGB(toBlend,:);
                
                set(hs,'FaceVertexCData',mixedRGB);
                
                if movie
                    Frame = getframe(gcf); %get(hs,'Parent')
                    mov = addframe(mov,Frame);    
                end
                set(hcurr,'xdata',Time(timee(t)),'ydata',ploot(timee(t)))
                drawnow
                keyboard
            end
            
            if movie 
                moc = close(mov);
            end
            
            keyboard
            
        case 'n'
            if ishandle(hfig)
                close(hfig)
            end
            
        end
        
    end
    
    bst_message_window({...
            ' ',...
            sprintf('Average Correlation : %3.2f %%',100*mean(abs(CorrCoeff))),...
            sprintf('Average area explored : %3.2f cm2%%',mean(area)),...
            ' '});
    
    %----------------------------------------------------------------------------
case 'CorrelationMapFromSensor' % Compute a correlation map between source time series and a the signal from a sensor
    % Useful for instance when looking for areas correlating the most with some depth electrode signal.
    CorrCoefThres = .98;
    data.Results.Display.Handles.hPatch  = findobj(0,'tag','tessellation_window'); 
    %CorticalScouts = get(data.Results.Display.Handles.hPatch,'Userdata');
    ScoutsManagerFig = findobj(0,'Tag','tesselation_select');
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    % First ask which sensor will work as the reference 
    % and on which time window
    Prompt = {...
            'NAME of sensor to act as reference in correlation computation: ',...
            'Time window for correlation computation (in ms):',...
        };
    
    Answer = inputdlg(Prompt,'',1,{'tobas4',''});
    
    if isempty(Answer{1})
        return
    end
    
    CorrRefChan = Answer{1}; % Name of reference channel for correlation computation
    CorrTimeWin = str2num(Answer{2})/1000; % Time window for correlation computation
    clear Answer
    
    % Find sensor in channels
    Visu = get(DATAPLOT,'Userdata');
    CorrRefChan = find(strcmpi(CorrRefChan,{Visu.Channel.Name}));
    if isempty(CorrRefChan), return, end % Channel name not found
    CorrRefChanName = Visu.Channel(CorrRefChan).Name;
    F = Visu.Data.F(CorrRefChan,:); 
    clear Visu 
    % Find time indices for correlation computation
    % Trim out source and sensor time series
    Results = get(hScoutsManagerFig.ResultFiles,'Userdata');
    
    F = F(Results.Time);
    
    TimeInd = findclosest(CorrTimeWin,Results.ImageGridTime');
    Results.ImageGridAmp = (Results.ImageGridAmp(:,TimeInd(1):TimeInd(end)));
    % Apply amplitude threshold
    M = max(abs(Results.ImageGridAmp),[],2); % Take maximum over time
    MAX = max(M);
    cThres = str2num(get(hScoutsManagerFig.TruncateFactor,'string'))/100;
    % Indices of the sources above threshold
    GoodSrc = find(M >= cThres*MAX);
    Results.ImageGridAmp  = Results.ImageGridAmp(GoodSrc,:);
    
    Results.ImageGridAmp = Results.ImageGridAmp - repmat(mean(Results.ImageGridAmp,2),1,length(TimeInd(1):TimeInd(end)));
    F = F(TimeInd(1):TimeInd(end));
    F = F - mean(F);
    Results.ImageGridAmp = double(Results.ImageGridAmp);
    
    %Compute correlation coefficients for all sources
    %     CorrCoeff = abs(F/norm(F))...
    %         *abs(Results.ImageGridAmp'*inorcol(Results.ImageGridAmp'));
    
    CorrCoeff = abs((F/norm(F))...
        *Results.ImageGridAmp'*inorcol(Results.ImageGridAmp'));
    
    % Keep only sources which correlate best
    % Indices of the sources above threshold
    ind = find(CorrCoeff >= CorrCoefThres);
    CorticalScouts.Vertex = GoodSrc(ind);
    if isempty(CorticalScouts.Vertex)
        [maxx ind] = max(CorrCoeff);
        CorticalScouts.Vertex = GoodSrc(ind);    
    end
    
    Source = Results.ImageGridAmp(ind,:);
    
    %     if length( CorticalScouts.Vertex) == 1
    %         figure, plot(F/norm(F),'r'), hold on, plot(Results.ImageGridAmp(ind,:)/norm(Results.ImageGridAmp(ind,:)),'g');
    %     end
    
    sprintf('Found %d sources correlating at least up to %3.2f with sensor %s', ...
        length(CorticalScouts.Vertex),...
        CorrCoefThres, ...
        CorrRefChanName)
    
    
    % Plot
    CorrCoeff = CorrCoeff(ind);
    for k = 1:length(CorticalScouts.Vertex)
        if abs(CorrCoeff(k)) > .4
            hfig = figure;
            TimeVec = 1000*linspace(CorrTimeWin(1),CorrTimeWin(2), length(F));
            ploot = 100*F/(norm(F,'fro')/sqrt(length(CorticalScouts.Vertex))); % Normalize amplitudes by average signal energy over all electrodes
            plot(TimeVec,ploot,'r','linewidth',2), hold on
            plot(TimeVec, Source(k,:)*max(abs(ploot))/max(abs(Source(k,:))),'b','linewidth',2);
            set(gcf,'color','k')
            set(gca,'xcolor','w','ycolor','w')
            ylabel('% average total signal energy')
            title(sprintf('Correlation: %3.2f %%',100*abs(CorrCoeff(k))),'color','w')
        end
    end
    
    CorticalScouts.CorrCoefThres = CorrCoefThres;
    CorticalScouts.CorrRefChan = CorrRefChan;
    CorticalScouts.CorrTimeWin = CorrTimeWin;
    
    %CorrCoefThres
    
    % Now generate cortical scout file
    
    % What's the name of the current cortical envelope ?
    ActiveTess = get(hScoutsManagerFig.removed,'String'); % Find the active Cortical surface
    iCortex = get(hScoutsManagerFig.removed,'Value'); 
    if iscell(ActiveTess)
        ActiveTess= ActiveTess{iCortex};
    end
    
    % XYZ Coordinates of the scouts
    Visu = get(DATAPLOT,'Userdata');
    cd(Users.SUBJECTS)
    load(Visu.Tesselation,'Comment','Vertices')
    % What's the current surface ?
    imesh = find(strcmp(Comment,ActiveTess));
    Vertices = Vertices{imesh};
    CorticalScouts.Handle = Vertices(:,CorticalScouts.Vertex)';
    clear Vertices
    for k =1 : length(CorticalScouts.Vertex)
        CorticalScouts.HandleLabels{k} = int2str(k);
    end
    
    if 0%length( CorticalScouts.HandleLabels) == 1
        % Load specific cortical scout file (depth electrode)
        DepthScouts = load('I:\Data\Studies\Epilepsie\02May02\DepthElectrodes_Scouts.mat');
        % Find scout with same name
        ScoutRef = find(strcmpi(CorrRefChanName,cellstr(char(DepthScouts.CorticalScouts.HandleLabels{:}))));
        % Add to current Cortical Scout File
        CorticalScouts.Handle = [CorticalScouts.Handle; DepthScouts.CorticalScouts.Handle(ScoutRef,:)];
        CorticalScouts.HandleLabels{end+1} = DepthScouts.CorticalScouts.HandleLabels{ScoutRef};
        CorticalScouts.Vertex(end+1) = DepthScouts.CorticalScouts.Vertex(ScoutRef);
        Dist = 1000*norm(diff(CorticalScouts.Handle)); % Distance in mm
        sprintf('Distance between depth electrode and MaxCorr spot : %3.2f',Dist)
    end
    
    % What's the name of the current data set ?
    ResultFiles = get(hScoutsManagerFig.ResultFiles,'String');
    ResultFile = ResultFiles{get(hScoutsManagerFig.ResultFiles,'Value')};
    
    cd(Users.STUDIES)
    dirr = fileparts(Users.CurrentData.StudyFile);
    cd(dirr);
    
    % Create Cortical Scout File
    [filename, pathname] = uiputfile(['CorrCoefScouts_',get(hScoutsManagerFig.TruncateFactor,'string'),'.mat'], 'Save Cortical Scout in...');
    save(fullfile(pathname,filename), 'ResultFile','ActiveTess','CorticalScouts','-mat');
    
    % Visualization
    dataplot_cb('LoadCorticalSpot',fullfile(pathname,filename))
    
    
    
    %----------------------------------------------------------------------------
    
case 'CorticalScoutFromThres' % Create a cortical scout file with areas above a given %-threshold
    data.Results.Display.Handles.hPatch  = findobj(0,'tag','tessellation_window'); 
    CorticalScouts = get(data.Results.Display.Handles.hPatch,'Userdata');
    ScoutsManagerFig = findobj(0,'Tag','tesselation_select');
    hScoutsManagerFig = guihandles(ScoutsManagerFig);
    
    % Let's do it simple first
    Results = get(hScoutsManagerFig.ResultFiles,'Userdata');
    M = max(abs(Results.ImageGridAmp),[],2); % Take maximum over time
    MAX = max(M);
    cThres = str2num(get(hScoutsManagerFig.TruncateFactor,'string'))/100;
    
    % Indices of the sources above threshold
    CorticalScouts.Vertex = find(M >= cThres*MAX);
    % What's the name of the current cortical envelope ?
    ActiveTess = get(hScoutsManagerFig.removed,'String'); % Find the active Cortical surface
    iCortex = get(hScoutsManagerFig.removed,'Value'); 
    if iscell(ActiveTess)
        ActiveTess= ActiveTess{iCortex};
    end
    
    % XYZ Coordinates of the scouts
    Visu = get(DATAPLOT,'Userdata');
    cd(Users.SUBJECTS)
    load(Visu.Tesselation,'Comment','Vertices')
    % What's the current surface ?
    imesh = find(strcmp(Comment,ActiveTess));
    Vertices = Vertices{imesh};
    CorticalScouts.Handle = Vertices(:,CorticalScouts.Vertex)';
    clear Vertices
    for k =1 : length(CorticalScouts.Vertex)
        CorticalScouts.HandleLabels{k} = int2str(k);
    end
    
    
    % What's the name of the current data set ?
    ResultFiles = get(hScoutsManagerFig.ResultFiles,'String');
    ResultFile = ResultFiles{get(hScoutsManagerFig.ResultFiles,'Value')};
    
    cd(Users.STUDIES)
    dirr = fileparts(Users.CurrentData.StudyFile);
    cd(dirr);
    
    % Create Cortical Scout File
    [filename, pathname] = uiputfile(['CorticalScoutFromThresh_',get(hScoutsManagerFig.TruncateFactor,'string'),'.mat'], 'Save Cortical Scout in...');
    save(fullfile(pathname,filename), 'ResultFile','ActiveTess','CorticalScouts','-mat');
    
    % Visualization
    dataplot_cb('LoadCorticalSpot',fullfile(pathname,filename))
    
end
