function varargout = bst_imaging_display(varargin)
%BST_IMAGING_DISPLAY - : display utility for cortical distributed source imaging 
% function varargout = bst_imaging_display(varargin)
%function varargout = bst_imaging_display(ImagingFigure,ImageGridLoc,ImageGridAmp,varargin)
% if nargin == 2
%       ImageGridAmp = varargin{1};
% 		OPTIONS = varargin{2};
%    
% ----------------------------------------------------------------------------------------------

%<autobegin> ---------------------- 12-Oct-2004 01:10:59 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\catci.m
%   toolbox\curvature_cortex.m
%   toolbox\dataplot_cb.m
%   toolbox\get_user_directory.m
%   toolbox\good_channel.m
%   toolbox\grayish.m
%   toolbox\interp_mail.m
%   toolbox\norcol.m
%   toolbox\norlig.m
%   toolbox\plot_dipole.m
%   toolbox\vertices_connectivity.m
%   toolbox\view_surface.m
%
% Subfunctions in this file, in order of occurrence in file:
%   [] = DataOnSensors % Sensors are color-coded according to the data at each location
%   [] = DataOnScalp % Interpolated Scalp Topography of Surface Data
%   [] = CorticalMap(TessWin,OPTIONS) % See cortical current maps interpolated on the proper cortical surface
%   [mixedRGB] = mix_curvature_current(varargin)
%
% Figure Files opened by this function:
%   'mapping.fig'
%   'tessellation_window.fig'
%
%   Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
%   <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure mapping.fig
%   figure::mapping "" uses ButtonDownFcn for dataplot_cb line_select
%   uicontrol:checkbox:ABSOLUTE "Absolute" uses Callback for dataplot_cb mutincomp_cmapping
%   uicontrol:checkbox:Colorbar "Colorbar" uses Callback for dataplot_cb Colorbar_mapping
%   uicontrol:checkbox:fit "Sensor array" uses Callback for dataplot_cb mapping_head_shape
%   uicontrol:checkbox:gradient "Gradient" uses Callback for dataplot_cb mapping_type
%   uicontrol:checkbox:magnetic "Mag. Field Estim." uses Callback for dataplot_cb mapping_type
%   uicontrol:checkbox:movie "Movie" uses Callback for dataplot_cb mapping_display_type
%   uicontrol:checkbox:raw "Raw Mapping" uses Callback for dataplot_cb mapping_type
%   uicontrol:checkbox:RELATIVE "Relative" uses Callback for dataplot_cb mutincomp_cmapping
%   uicontrol:checkbox:scalp "Scalp surface" uses Callback for dataplot_cb mapping_head_shape
%   uicontrol:checkbox:ShowSensorLabels "Sensor Labels" uses Callback for dataplot_cb mapping_type
%   uicontrol:checkbox:single "Single" uses Callback for dataplot_cb mapping_display_type
%   uicontrol:checkbox:slides "Slides" uses Callback for dataplot_cb mapping_display_type
%   uicontrol:checkbox:sphere "Sphere" uses Callback for dataplot_cb mapping_head_shape
%   uicontrol:edit:current_time "0" uses Callback for dataplot_cb set_current_time
%   uicontrol:pushbutton:go "Apply" uses Callback for dataplot_cb gomapping
%   uicontrol:pushbutton:Pushbutton1 "Quit" uses Callback for delete(gcbf)
%   uicontrol:slider:slider_time "" uses Callback for dataplot_cb mapping_slider
%
% Callbacks by figure tessellation_window.fig
%   figure::tessellation_window "" uses CreateFcn for movegui northwest
%
% At Check-in: $Author: Mosher $  $Revision: 28 $  $Date: 10/11/04 11:32p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 23-Sep-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> ------------------------ 12-Oct-2004 01:10:59 -----------------------

% EK 4-Oct-2004   Minor change in displaying the cortex: If the
%                 ImageGridAmp is a vector (like in Neural Index Case) then 
%                 display this for any time point selected in time series window. 
% ----------------------------------------------------------------------------


if nargin == 4
    OPTIONS = varargin{1};
end

Users = get_user_directory;

% Define default values for display options
Def_CutPlane_OPTIONS = struct(...
  'Top', 0, ...
  'Bottom', 0, ...
  'Front', 0, ...
  'Back', 0,...
  'Right', 0,...
  'Left', 0,...
  'X', [],...
  'Y', [],...
  'Z', [] ...
  );

Def_OPTIONS = struct(...
    'CurvatureColor1',[.4 .4 .4],...
    'CurvatureColor2',[.6 .6 .6],...
    'Dipoles',[],...
    'OrthoViews',0,...
    'Time', 1, ...
    'ScalpMaps',0,...
    'ColorSensors',0,...
    'CorticalMap',1,...
    'ColorMapping','normalized',...
    'CDThreshold',.5,...
    'currDensityTransparency',.2, ...
    'CurvatureThreshold',.16,...
    'FaceVertexAlpha',.2,...
    'hPatch', NaN,...
    'MapCurvature',1,... 
    'timeText',NaN,...
    'timeValue',[],...
    'Curvature', [], ...
	'InitialCurvature',[],...
	'FaceAlpha',1,...
    'VertConn', [],...
    'CutPlane', Def_CutPlane_OPTIONS...
    );

% Check arguments passed to and required from function call
if nargin < 3 % return default OPTIONS values
    
    if nargout > 1
        varargout{1} = Def_OPTIONS;
        varargout{2} = Def_OPTIONS;
    else
        varargout{1} = Def_OPTIONS;
    end
    if nargin == 2
        ImageGridAmp = varargin{1};
        if isempty(ImageGridAmp)
            ImageGridAmp = 0; % no current mapping - just a look at anatomy
        end
		OPTIONS = varargin{2};
    else
        return
    end
    
else
    
    ImagingFigure = varargin{1};
    ImageGrid = varargin{2};
    ImageGridAmp = varargin{3};
	OPTIONS = varargin{4};
    
end

% Check field names of passed OPTIONS and fill missing ones with default values
DefFieldNames = fieldnames(Def_OPTIONS);
for k = 1:length(DefFieldNames)
    if ~isfield(OPTIONS,DefFieldNames{k})
        
        OPTIONS = setfield(OPTIONS,DefFieldNames{k},getfield(Def_OPTIONS,DefFieldNames{k}));
        
    end
end

clear Def_OPTIONS

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

% tmp  = findobj(0,'tag','subplot_tess_window');
% if ~isempty(tmp)
%     ImagingFigure = tmp; % Consider the slides or movie window
% end

if ~ishandle(OPTIONS.hPatch) % no surface available: load vertex/faces information
    % Load surface information ------------------
    try
        load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1})) % CBB - suppose there is one surface per tessellation file
    catch
        try
            load(ImageGrid.GridLoc{1})
        catch
            errordlg(sprintf('Subject file %s could not be found. Please update subject''s entry in database or manually update the subject information in selected results file.',fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1})),'Unvalid subject database entry');
            return
        end
    end
    
    FigName = Comment{ImageGrid.iGrid(1)}; clear Comment % Name of figure where to plot the surface
    %OPTIONS.SurfaceName = FigName;
    fv.faces = Faces{ImageGrid.iGrid(1)}; clear Faces
    fv.vertices = Vertices{ImageGrid.iGrid(1)}'; clear Vertices
    nverts_max = 10000; % Max number of vertices before warning 
    if size(fv.vertices,1) > 10000 
        ButtonName=questdlg(sprintf('Surface has %d vertices; would you like to downsample this number ?',size(fv.vertices,1)), ...
            'Large Surface Detected', ...
            'Yes','No','Yes');
        switch(ButtonName)
        case 'Yes'
            hmsg = msgbox('Downsampling original surface. . .','Working!','warn'); drawnow
            fv = reducepatch(fv,10000/size(fv.vertices,1),'fast','verbose');
            delete(hmsg);
        otherwise
            % do nothing
        end
    end
    
    % Load surface information----DONE----------
    figure(ImagingFigure)
    % Figure Parameters ------------------------
    set(ImagingFigure,'CurrentAxes',findobj(ImagingFigure,'Tag','MainAxis'),...
        'Renderer','OpenGL',...
        'Color','w', 'Colormap',gray,'Menubar','figure','Name',FigName);
else
    fv.vertices = get(OPTIONS.hPatch,'Vertices');
end



%         % Surface orientation
%         if ~isfield(Visu.Data,'System')
%             Visu.Data.System = 'ctf';
%         end

% Figure Parameters -----------DONE---------


if 0 % deprecated code (Sb 08-Oct-2003)
    
    %FragmentMenu = findobj(gcbf,'Tag','FragmentMenu');
    %nclust = get(FragmentMenu,'Value')-1; % Get label of clustering to visualize 
    nclust = 0;
    
    hpatch = []; % handle(s) to patch objects in figure
    if OPTIONS.OrthoViews == 0
        hpatch = findobj(findobj(ImagingFigure,'Tag','MainAxes'),'Type','patch','Visible','on'); % CBB - store these handles in figure appdata ? (SB)
    else
        if ~isempty(findobj(findobj(ImagingFigure,'Tag','axsub1'),'Type','patch'))
            hpatch(1) = findobj(findobj(ImagingFigure,'Tag','axsub1'),'type','patch');
            hpatch(2) = findobj(findobj(ImagingFigure,'Tag','axsub2'),'type','patch');
            hpatch(3) = findobj(findobj(ImagingFigure,'Tag','axsub3'),'type','patch');
            hpatch(4) = findobj(findobj(ImagingFigure,'Tag','axsub4'),'type','patch');
        else
            if isempty(hpatch)
                hpatch = [];
            end
        end
    end
end % Deprecated code -------------------------


% Work on cut planes --------------------------

% work with facelapha property only?


% Work on cut planes -------------DONE---------



% Compute or retrieve curvature information --
if OPTIONS.MapCurvature % Color coding of surface curvature 
    
    if isempty(OPTIONS.Curvature) % No curvature mapping is available - compute it
        
        % Is vertex connectivity available ?
        %load(fullfile(Users.SUBJECTS,Users.CurrentData.SubjectFile),'VertConn'); 
        try
            load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'VertConn');
        catch
            ImageGrid.GridLoc{1} = strrep(ImageGrid.GridLoc{1},Users.SUBJECTS,'');
        end
        load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'VertConn');
       
        if ~exist('VertConn','var') % Vertex connectivity not computed yet
            load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Comment');
            VertConn = cell(length(Comment),1);
        end
          
        %VertConnFile = VertConn; % Filename containing connectivity information
        [pathname, filename] = fileparts(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}));    
        
        cd(pathname)
        %if ~exist(VertConnFile,'file')
        if isempty(VertConn{ImageGrid.iGrid(1)})
    
            % How many surfaces in tessellation file ?
            % should be only one but older BrainStorm files may contain several
            load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Comment');
            nsurfaces = length(Comment); % Number of surfaces in tessellation file
            VertConn = cell(nsurfaces,1); % Create cell array where to store vertex connectivity values
  
            VertConn{ImageGrid.iGrid(1)} = vertices_connectivity(fv); % compute vertex connectivity 
            [pathname, filename] = fileparts(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}));    
            cd(pathname)
            save(ImageGrid.GridLoc{1},'VertConn','-append');  
            
        else % load existing vertex connectivity map
        
            %load(fullfile(pathname,VertConnFile))
            
        end
        
        % Now compute curvature information
        OPTIONS.VertConn = VertConn{ImageGrid.iGrid(1)}; clear VertConn
        
        %[vertexcolorh,map]=c_activ(FV,TessWinUserData.VertConn); clear FV ; OPTIONS.Curvature = vertexcolorh;
        % Parameters for curvature computation
        sigmoid_const = .1;  
        show_sigmoid = 0;
       
		if isempty(OPTIONS.InitialCurvature)
            load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Curvature');
            if ~exist('Curvature','var') % Vertex connectivity not computed yet
                load(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Comment');
                Curvature = cell(length(Comment),1);
            end
            
            if isempty(Curvature{ImageGrid.iGrid(1)}) | size(Curvature{ImageGrid.iGrid(1)},2) > 1 % Older curvature computation
                
                bst_message_window('wrap',...
                    {' ','Computing surface curvature. . .'})
                [curvature_sigmoid,Curvature{ImageGrid.iGrid(1)}] = curvature_cortex(fv,OPTIONS.VertConn,sigmoid_const,show_sigmoid);
                bst_message_window('wrap',...
                    {'Computing surface curvature ->DONE',' '})
                clear curvature_sigmoid

                save(fullfile(Users.SUBJECTS,ImageGrid.GridLoc{1}),'Curvature','-append');
                
            end
            
		end
        OPTIONS.InitialCurvature = Curvature{ImageGrid.iGrid(1)}; clear Curvature
		
	end
	
	% apply threshold
	OPTIONS.Curvature(OPTIONS.InitialCurvature >= OPTIONS.CurvatureThreshold) = .8;
	OPTIONS.Curvature(OPTIONS.InitialCurvature < OPTIONS.CurvatureThreshold) = -.8;
	   
end

% Compute or retrieve curvature information --


% Colormap for current density---
tmp = hot(140);%flipdim(hot(100),1);
currDensityColorMap = tmp(1:100,:); % Current density colormap
%anatColorMap = repmat(linspace(.4,.6,100)',1,3); % Anatomy colormap
anatColorMap = [OPTIONS.CurvatureColor1;...
        zeros(98,3);...
        OPTIONS.CurvatureColor2]; % Anatomy colormap - articifical zero padding 
% Colormap ----------DONE--------

if ImageGridAmp == 0 % No activation is specified - look at anatomy only
    ImageGridAmp = ones(size(fv.vertices,1),length(OPTIONS.Time));
end

% Create patch object -------------------------
if ~ishandle(OPTIONS.hPatch) % No patch handle is passed (i.e. need to create one)
    figure(ImagingFigure)
    PlotHandles.SurfaceName = FigName;
    [PlotHandles.hFigure,PlotHandles.hPatch,PlotHandles.hLights] = ...
        view_surface(FigName,fv.faces,fv.vertices,ImageGridAmp(:,OPTIONS.Time));
    % Dull rendering
    delete(findobj(PlotHandles.hFigure,'type','light'))
    %view(-10,0), 
    material dull
    %PlotHandles.hLights = camlight;
    set(PlotHandles.hPatch,'backfacelighting','unlit')
    set(gca,'AmbientLightColor',[1 1 1])
    hold on
	rotate3d(PlotHandles.hFigure,'on')
    %set(hf,'color','k')
else
    PlotHandles.hPatch = OPTIONS.hPatch;
end

    
% Create patch object ----------DONE-----------


% Translate anatomy curvature in RGB values
if ~OPTIONS.MapCurvature % no color coding of surface curvature 
    if ishandle(PlotHandles.hPatch)
        fv.vertices = get(PlotHandles.hPatch,'Vertices');
    end
    
    index_anatomy = 100*ones(size(fv.vertices,1),1);
       
else
    min_curv = min(OPTIONS.Curvature);
    max_curv = max(OPTIONS.Curvature);
    index_anatomy = floor((OPTIONS.Curvature-min_curv)*99/(max_curv-min_curv)+1);
    
end

anatRGB = anatColorMap(index_anatomy,:);

% Translate current density values in RGB 

% Amplitude scaling
ImageGridAmp = double(ImageGridAmp);
switch lower(OPTIONS.ColorMapping) % is color mapping time-relative or time-normalized ?
case 'normalized'
    if size(ImageGridAmp,2) == 1 %E.g Representation of  neural activity index
       if OPTIONS.Time >1
           bst_message_window('wrap',{'Currently visualized result does not change with time',' '});
       end

        currDensity = abs(ImageGridAmp(:,1))/max(abs(ImageGridAmp(:,1)));
    else
    currDensity = abs(ImageGridAmp(:,OPTIONS.Time))/max(abs(ImageGridAmp(:,OPTIONS.Time)));
    end
    limitValue = max(currDensity); % Store maximum value 
case 'relative'
        if size(ImageGridAmp,2) == 1 %E.g Representation of  neural activity index
                   if OPTIONS.Time >1
           bst_message_window('wrap',{'Currently visualized result does not change with time',' '});
       end
    currDensity = abs(ImageGridAmp(:,1))/max(abs(ImageGridAmp(:)));

        else
    currDensity = abs(ImageGridAmp(:,OPTIONS.Time))/max(abs(ImageGridAmp(:)));
        end
    limitValue = 1; % Force it to 1
end
% Apply amplitude threshold
%tmpcurrDensity = currDensity;
currDensity(currDensity < OPTIONS.CDThreshold) = 0; 
if limitValue == min(currDensity) % no current value above threshold
    limitValue = -1 ;
end
%clear tmpcurrDensity
if 1%(viewmode == 1)%if real values
    index_density = floor((currDensity + limitValue)*99/(2*limitValue)+1);
else %if abs values
    currDensity = abs(currDensity);
    index_density = floor((currDensity-limitValue)*99/limitValue+100);
end
currDensityRGB = currDensityColorMap(index_density,:);


% Now mix curaveture and current RGBs
mixedRGB = anatRGB;
toBlend = find(currDensity ~= 0); % Find vertex indices holding non-zero activation (after thresholding)
mixedRGB(toBlend,:) = OPTIONS.currDensityTransparency*anatRGB(toBlend,:)+...
    (1-OPTIONS.currDensityTransparency)*currDensityRGB(toBlend,:);


% Assign color to surface elements

set(PlotHandles.hPatch,'FaceVertexCData',mixedRGB);
drawnow

% Add text label to indicate current time
if ~isempty(OPTIONS.timeValue)
    if ~ishandle(OPTIONS.timeText)
        OPTIONS.timeText = text(0.05,0.05,0,sprintf('%3.1f ms',1000*OPTIONS.timeValue),'units','normalized');
        set(OPTIONS.timeText,'fontweight','normal','color',[.3 .3 .3],'Fontname','helvetica','Fontsize',8,'FontUnits','Point');
    else
        set(OPTIONS.timeText,'String',sprintf('%3.1f ms',1000*OPTIONS.timeValue))
    end
end

% Now overlay dipole models if any -------------------------------------
%rotate3d on

if isempty(OPTIONS.Dipoles)
    varargout{1} = PlotHandles;
    varargout{2} = OPTIONS;
    %rotate3d(get(PlotHandles.hPatch,'Parent'),'on')
    return
end

%set(PlotHandles.hPatch,'edgecolor','none','FaceAlpha',0.3); hold on
% Make transparency decrease with depth
%fv.vertices = fv.vertices - repmat( mean(fv.vertices),size(fv.vertices,1),1);
%  [TH,PHI,R] = CART2SPH(fv.vertices(:,1),fv.vertices(:,2),fv.vertices(:,3)); clear TH PHI
%  R = R/max(R);
% set(PlotHandles.hPatch,'edgecolor','none','FaceVertexAlpha',1-R,'facealpha','flat'); hold on

axes(get(PlotHandles.hPatch,'Parent'))
set(PlotHandles.hPatch,'FaceVertexAlpha',OPTIONS.FaceVertexAlpha,'facealpha','flat'); hold on

scale = norcol(OPTIONS.Dipoles.TimeSeries); % Scale factor for amplitude display
scale = scale/max(scale); 
scale(scale <.5) = .5; % Avoid too small dipoles
scale = scale/200; % Dipole size needs to fit within cortical volume

for src = 1:length(OPTIONS.Dipoles.SourceLoc)
    axes(get(PlotHandles.hPatch,'Parent'))
    if isempty(OPTIONS.Dipoles.Color{src})
        OPTIONS.Dipoles.Color{src} = rand(1,3);
    end
    plot_dipole(OPTIONS.Dipoles.SourceLoc{src}',...
        OPTIONS.Dipoles.Orientation{src}',...
        scale(src),OPTIONS.Dipoles.Color{src});
end


varargout{1} = PlotHandles;
varargout{2} = OPTIONS;

rotate3d(PlotHandles.hFigure,'on')

return

% Deprecated code below

% Call subfunctions depending on type fo display
if OPTIONS.ColorSensors
    DataOnSensors
end
if OPTIONS.ScalpMaps
    DataOnScalp
end
if OPTIONS.CorticalMap
    %CorticalMap(ImagingFigure,hpatch)
else
    if get(OPTIONS.MapCurvature,'Value')
        set(hpatch,'FaceVertexCData',vertexcolorh)
    end
end

set(hpatch,'Visible','on')

vertexcolorh =  get(hpatch,'FaceVertexCdata');

%         if nclust > 0%k == nsurf & nclust > 0 % Some clustering has been selected - display clusters on the surface
%             set(h,'FacevertexCdata',Clusters{nsurf}.Cluster{nclust},'Userdata',Clusters{nsurf}.Cluster{nclust}')
%             colormap(rand(Clusters{nsurf}.NClusters(nclust),3))
%         end

hold on
        
set(hpatch,'edgecolor','none','facealpha',OPTIONS.FaceAlpha)

% Update envelope properties ------------------------------------------------------

ctet = get(hpatch,'vertices');
siz = size(vertexcolorh,2);
if size(vertexcolorh,1) == 1
    vertexcolorh = ones(size(ctet,1),1)*vertexcolorh;
end


if 0
    
    if Vals(k,3) == 1
        vertexcolorh(find(ctet(:,2) < 0),:) = NaN * ones(length(find(ctet(:,2) < 0)),siz);	   	
    end
    if Vals(k,4) == 1
        vertexcolorh(find(ctet(:,2) > 0),:) = NaN * ones(length(find(ctet(:,2) > 0)),siz);	   	
    end
    
    
    if Vals(k,3) == 1 | Vals(k,4) == 1 % Separate both hemispheres
        
        % centering of the brain
        ctet2(:,1)=ctet(:,1)- mean(ctet(:,1));%sum(sum(ctet(:,1)))/size(ctet,1).*ones(size(ctet,1),1);
        ctet2(:,2)=ctet(:,2)- mean(ctet(:,2));%sum(sum(vu(:,2)))/size(vu,1).*ones(size(vu,1),1);
        ctet2(:,3)=ctet(:,3)- mean(ctet(:,3));;%sum(sum(vu(:,3)))/size(vu,1).*ones(size(vu,1),1);
        %             ctet2(:,1)=ctet(:,1)- sum(sum(ctet(:,1)))/size(ctet,1).*ones(size(ctet,1),1);
        %             ctet2(:,2)=ctet(:,2)- sum(sum(ctet(:,2)))/size(ctet,1).*ones(size(ctet,1),1);
        %             ctet2(:,3)=ctet(:,3)- sum(sum(ctet(:,3)))/size(ctet,1).*ones(size(ctet,1),1);
        %             
        %             % find short and long axes
        %             larg1=max(max(ctet2(:,1)))-min(min(ctet2(:,1)));
        %             larg2=max(max(ctet2(:,2)))-min(min(ctet2(:,2)));
        %             if larg1<larg2
        %                 short_axe=1;
        %             else
        %                 short_axe=2;
        %             end
        %             
        %             cd(Users.SUBJECTS)
        %             load(Users.CurrentData.SubjectFile,'VertConn'); % Is the vertex connectivity available ?
        %             load(VertConn)
        %             xVertConn = VertConn;
        %             if isempty(xVertConn{nsurf}) % Vertex connectivity not available: compute and save.
        %                 FV.faces = get(h,'faces');
        %                 FV.vertices = get(h,'vertices');
        %                 tmp = vertices_connectivity(FV); clear FV
        %                 xVertConn{nsurf} = tmp; clear tmp
        %                 load(Users.CurrentData.SubjectFile,'VertConn'); 
        %                 VertConnFile= VertConn; VertConn = xVertConn; clear xVertConn
        %                 save(VertConnFile,'VertConn','-append')
        %                 VertConn = VertConn{nsurf};
        %             else
        %                 VertConn = xVertConn{nsurf};
        %             end
        %             
        
        %             [I_pos,I_neg] = separation(ctet2,VertConn); clear ctet2
        
        %SVD the vertex coordinates to find principal axes
        %Dowsize the vertex location matrix to run the SVD and find the axis of principal direction for the cortical volume
        nverts = size(ctet2,1);
        FV.faces = get(h,'faces');
        FV.vertices = get(h,'vertices');
        nfv = reducepatch(FV,1000); clear FV
        tmp = nfv.vertices; clear nfv
        tmp = tmp - repmat(mean(tmp),size(tmp,1),1);
        ndown = 500; % Objective for the number of vertices considered for determination of the main axes.
        rate = ceil(nverts/ndown);
        tmp = ctet2(1:rate:end,:);
        [U,s,V] = svd(tmp'/sqrt(size(tmp,1)-1),0); % Compute the SVD of the covariance coordinate matrix.
        s = sqrt(diag(s)); % Get the standard deviation in each direction
        
        % Draw orthogonal planes to illustrate the 3 orthogonal axes
        % Brain centroid
        brain_c = mean(ctet);
        % Get the axes limits;
        axh = get(h,'Parent'); % Get current surface parent
        Xlims = get(axh,'Xlim'); Ylims = get(axh,'Ylim'); Zlims = get(axh,'Zlim');
        
        % Parametric equation for all lines passing by brain_c in 3D: X = brain_c + U(:,i)*t
        fline = inline('brain_c'' + u*x','x','brain_c','u');% Create a function for the parametric expression of the line
        
        % Sagittal Line
        p1 = fline(-.1,brain_c,U(:,1))';
        p2 = fline(.1,brain_c,U(:,1))';
        
        hline(1) = line([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)]);
        
        % Sagittal Line
        p1 = fline(-.1,brain_c,U(:,2))';
        p2 = fline(.1,brain_c,U(:,2))';
        
        hline(2) = line([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)]);
        
        % Sagittal Line
        p1 = fline(-.1,brain_c,U(:,3))';
        p2 = fline(.1,brain_c,U(:,3))';
        
        hline(3) = line([p1(1) p2(1)],[p1(2) p2(2)],[p1(3) p2(3)]);
        
        %             % Select the clustering seed in one of the hemisphere: pick up the extreme point left or right vertex
        %             if short_axe == 1
        %                 [mm, imax] = min(abs(ctet2(:,1)));
        %             else
        %                 [mm, imax] = min(abs(ctet2(:,2)));
        %             end
        %             
        %             [CLASS, NumClass] = crtx_cluster(VertConn,2,imax,1);
        %             
        
        % Get corresponding t values 
        
        
        if Vals(k,3) == 1
            I_pos = find(CLASS{1} == 1);
            vertexcolorh(I_pos,:) = NaN * ones(length(I_pos),siz);	   	
        elseif Vals(k,4) == 1
            I_neg = find(CLASS{1} == 2);
            vertexcolorh(I_neg,:) = NaN * ones(length(I_neg),siz);	   	
        end
        
        if Vals(k,5) == 1
            vertexcolorh(find(ctet(:,1) > 0),:) = NaN * ones(length(find(ctet(:,1) > 0)),siz);	   	
        end
        if Vals(k,6) == 1
            vertexcolorh(find(ctet(:,1) < 0),:) = NaN * ones(length(find(ctet(:,1) < 0)),siz);	   	
        end
        if Vals(k,7) == 1
            vertexcolorh(find(ctet(:,3) > 0),:) = NaN * ones(length(find(ctet(:,3) > 0)),siz);	   	
        end
        if Vals(k,8) == 1
            vertexcolorh(find(ctet(:,3) < 0),:) = NaN * ones(length(find(ctet(:,3) < 0)),siz);	   	
        end
    end
end

clear ctet

if iscell(vertexcolorh)
    set(h,'FaceVertexCData',vertexcolorh{1},'Visible','on')
    if length(hpatch)>1 % OrthoViews
        for kk = 1:length(hpatch)
            h = hpatch(kk);
            set(h,'FaceVertexCData',vertexcolorh{1},'Visible','on') 
        end
        
    end
    
else
    if OPTIONS.OrthoViews == 1
        set(hpatch,'FaceVertexCData',vertexcolorh,'Visible','on')
    else
        set(hpatch,'FaceVertexCData',vertexcolorh,'Visible','on')
    end
end
vertexcolorh = []; % Move to next surface       


figure(ImagingFigure)
hold off
rotate3d(PlotHandles.hFigure,'on')
if ~isempty(findobj(gcf,'Tag','MainAxes')) 
    if length(hpatch)>1 & length(imesh) == 1% suplots with only one surface to display (need to generalize this)
        if ~strcmp(get(hpatch(4),'Visible'),'on')
            set(gca,'View',get(findobj(gcf,'Tag','MainAxes'),'view'))
        end
    end
    
end

%dataplot_cb mesh_lighting_props

%------------------------------------ SubFunctions

function [] = DataOnSensors % Sensors are color-coded according to the data at each location

%dataplot_cb mapping_create
Visu = get(DATAPLOT,'Userdata');
previous  = openfig('tessellation_window.fig','reuse');


% Handles to the spherical representation of the sensors 
sph = findobj(previous,'Tag','SENSORS');
if isempty(sph) | isempty(get(previous,'userdata'))
    
    if ~isempty(sph)
        if strcmp(get(sph(1),'type'),'surface')
            sph = dataplot_cb('see_sensors');
        else
            sph = dataplot_cb('see_sensors_markers');
        end
    else
        TessWin = findobj(0,'Tag','tesselation_select');
        hTessWin = guihandles(TessWin);
        if get(hTessWin.Sensors3D,'Value')
            sph = dataplot_cb('see_sensors');
        else
            sph = dataplot_cb('see_sensors_markers');    
        end
    end
    set(previous,'UserData',sph);
else
    sph = get(previous,'Userdata'); % Need this and line above to keep sensor ordering correct - sph = findobj creates an arbitrary ordering otherwise
end

% Bring the time slider - call the mapping window if available
MAPPING = openfig('mapping.fig','reuse');
hMAPPING = guihandles(MAPPING);

% Good/Bad Channels
if ~isempty(find(Visu.Data.ChannelFlag == -1))
    goodchannels = good_channel(Visu.Channel,Visu.Data.ChannelFlag,DataType{current}); % Discard bad channels
    badchannels = setdiff(good_channel(Visu.Channel,ones(length(Visu.Channel),1),DataType{current}),goodchannels);
    %      goodchannels = setdiff(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));
    %badchannels = intersect(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));
else
    goodchannels = setdiff(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));
    badchannels = [];
end

% Sampling rate
srate = abs(Visu.Data.Time(1)-Visu.Data.Time(2) );   % Assuming time is in sec.
ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000-Visu.Data.Time(1))/srate)+1;
axx = get(sph(1),'Parent');

data = Visu.Data.F(goodchannels,:); clear Visu
%set(axx,'Clim',[-max(abs(data(:))) max(abs(data(:)))]);
caxis(axx,[-max(abs(data(:))) max(abs(data(:)))]);
data = data(:,ctime);

C = [flipud(fliplr(hot(128)));hot(128)];
colormap(axx,C);

% Manipulation to get the color properly for each sensor sphere

set(sph,'visible','off')
TessWin = findobj(0,'Tag','tesselation_select');
hTessWin = guihandles(TessWin);

if get(hTessWin.Sensors3D,'Value')
    zdata = get(sph(1),'Zdata');                
    for k=1:length(sph)
        cdata = data(k)*ones(size(get(sph(k),'zdata')));
        set(sph(k),'CData',cdata,'facecolor','flat','CDataMapping','scaled');
    end
else
    for k=1:length(sph)
        set(sph(k),'Markerfacecolor',C( max([round(size(C,1)*(data(k)-min(caxis(axx)))/(max(caxis(axx))-min(caxis(axx)))),1]),:));
    end
    
end

set(sph,'visible','on')

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

function [] = DataOnScalp % Interpolated Scalp Topography of Surface Data

%dataplot_cb mapping_create

TessSelect = findobj(0,'Tag','tesselation_select');
hSelect = guihandles(TessSelect);
ActiveTess = get(hSelect.removed,'String'); % Find the active scalp surface
% Check if the CORTEX keyword is present
for k = 1:length(ActiveTess)
    flag(k) = ~isempty(findstr(lower(ActiveTess{k}),'head'));%'scalp'));
end
iScalp = find(flag);
iScalp = iScalp(1);
if isempty(iScalp)
    h = msgbox('Please select a scalp envelope from the tessellated surfaces available.');
    return
end
TessWin = openfig('tessellation_window.fig','reuse');

Visu = get(DATAPLOT,'Userdata');
if ~isempty(find(Visu.Data.ChannelFlag == -1))
    goodchannels = good_channel(Visu.Channel,Visu.Data.ChannelFlag,DataType{current}); % Discard bad channels
    badchannels = setdiff(good_channel(Visu.Channel,ones(length(Visu.Channel),1),DataType{current}),goodchannels);
    %      goodchannels = setdiff(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));% - min(Visu.ChannelID{current})+1; % Discard bad channels
    %badchannels = intersect(Visu.ChannelID{current},find(Visu.Data.ChannelFlag == -1));% - min(Visu.ChannelID{current})+1;
else
    goodchannels = Visu.ChannelID{current};
    badchannels = [];
end

Scalp = findobj(TessWin,'Type','patch','Tag',ActiveTess{iScalp}); % Need to find a better way to identify scalp when multiple surfaces are present in tessellation_win
if isempty(Scalp), tesselation_select_done, return, end

if isempty(get(Scalp,'Userdata'))
    figure(TessWin)
    C = [flipud(fliplr(grayish(hot(128),.3)));grayish(hot(128),.3)];
    colormap(C)
    
    scalp.vertices = get(Scalp,'vertices');
    % Channels Locations
    j = 0;
    for chan = goodchannels 
        j = j+1;
        %Projection of sensors on surface
        celec = ones(size(scalp.vertices,1),1)*Visu.Channel(chan).Loc(:,1)';
        dist = norlig(scalp.vertices-celec);
        [minn(j) imin(j)] = min(dist);
    end
    scalp.chanloc = [scalp.vertices(imin,1),scalp.vertices(imin,2),scalp.vertices(imin,3)];
    
    % Sampling rate
    scalp.srate = abs((Visu.Data.Time(2)-Visu.Data.Time(1))); 
    data = Visu.Data.F(goodchannels,:);
    axx= get(Scalp,'Parent');
    set(axx,'Clim',[-max(abs(data(:))) max(abs(data(:)))]);    
    set(Scalp,'Userdata',scalp)
else
    scalp = get(Scalp,'Userdata');
end


% Get the data from the mapping window
MAPPING = openfig('mapping.fig','reuse');
hMAPPING = guihandles(MAPPING);
ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000-Visu.Data.Time(1))/scalp.srate)+1;
vertxcolor = interp_mail(scalp.vertices,scalp.chanloc,Visu.Data.F(goodchannels,ctime));
set(Scalp,'FaceVertexCData',vertxcolor,'facecolor','interp');

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

function [] = CorticalMap(TessWin,OPTIONS) % See cortical current maps interpolated on the proper cortical surface

% TessSelect = findobj(0,'Tag','tesselation_select');
% hSelect = guihandles(TessSelect);
% ActiveTess = get(hSelect.removed,'String'); % Find the active scalp surface
% iCortex = get(hSelect.removed,'Value'); % Find the active scalp surface    

% if isempty(iCortex)
%     h = msgbox('Please select a cortex surface from the tessellations available.');
%     return
% end

%TessWin = findobj(0,'Tag','tessellation_window');
SlidesWin = findobj(0,'Tag','subplot_tess_window');
if ~isempty(SlidesWin)
    TessWin = SlidesWin; % Subplot or movie window
end

% if isempty(TessWin)
%     TessWin = openfig('tessellation_window.fig');
%     dataplot_cb('loaddata',Users.CurrentData.StudyFile,bst_static_taskbar('GET','DATA')); % Refresh by reloading data 
% end

% Check if loaded Results match this surface (ie number of sources = number of vertices of selected cortical surface)
Results = get(hSelect.ResultFiles,'Userdata');
nSources = size(Results.ImageGridAmp,1);

if nargin == 1 & ~OPTIONS.OrthoViews
    if iscell(ActiveTess)
        Cortex = findobj(get(TessWin,'CurrentAxes'),'Type','patch','Tag',ActiveTess{iCortex}); %CHEAT - need to be improved ?
    else
        Cortex = findobj(get(TessWin,'CurrentAxes'),'Type','patch','Tag',ActiveTess); %CHEAT - need to be improved ?    
    end
elseif nargin == 1 & get(hSelect.OrthoViews,'Value') == 1 % Orthogonal views are requested
    Cortex = findobj(get(TessWin,'CurrentAxes'),'Type','patch','Tag',ActiveTess,'Visible','on'); %CHEAT - need to be improved ?    
    Cortex = Cortex(1);
else
    Cortex = varargin{1}(1);
end

%set(Cortex,'Visible','off')
if ~isfield(get(Cortex,'Userdata'),'srate')
    Visu = get(DATAPLOT,'Userdata');
    figure(TessWin)
    
    cortex = get(Cortex,'Userdata');
    
    cortex.vertices = get(Cortex,'vertices');
    % Sampling rate
    cortex.srate = abs(Visu.Data.Time(2)-Visu.Data.Time(1)); 
    data = Results.ImageGridAmp;
    axx= get(Cortex,'Parent');
    set(axx,'Clim',[-max(abs(data(:))) max(abs(data(:)))]);    
    set(Cortex,'Userdata',cortex)
    
    % Result time window has priority over the original time window for the data - replace
    Visu.Data.Time = Results.ImageGridTime; % Need to do better than that - CHEAT
    try
        Visu.Data.F = Visu.Data.F(:,Results.Time(1):Results.Time(end)); % Should work only the first time the result time series are called - otherwise do nothing
    catch
        Visu.Data.F = Visu.Data.F; % Do nothing
    end
    
    set(DATAPLOT,'Userdata',Visu);
    
    M = max(abs(Results.ImageGridAmp(:)));
    set(hSelect.Colorbar,'Userdata',M);
    
else
    cortex = get(Cortex,'Userdata');
end

%dataplot_cb mapping_create

% Get the data from the mapping window
MAPPING = findobj(0,'Tag','mapping');
hMAPPING = guihandles(MAPPING);

if iscell(cortex)
    ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000- Results.ImageGridTime(1))/cortex{1}.srate)+1;
else
    ctime = round((get(hMAPPING.hDataplot.TimeSlider,'Value')/1000- Results.ImageGridTime(1))/cortex.srate)+1;
end
if ctime < 0, ctime =1; end

if get(hSelect.AbsoluteCurrent,'Value')
    set(Cortex,'FaceVertexCData',abs(Results.ImageGridAmp(:,ctime)));
else
    set(Cortex,'FaceVertexCData',Results.ImageGridAmp(:,ctime));
end

figure(TessWin)
hTessWin = guihandles(TessSelect);
if get(hTessWin.Normalize,'Value')
    TessWin = findobj(0,'Tag','tessellation_window');
    axx = findobj(TessWin,'tag','MainAxes');
    set(axx,'ClimMode','auto')  
    M = max(abs(Results.ImageGridAmp(:,ctime)));
    set(hSelect.Normalize,'Userdata',M)
    
    SlidesWin = findobj(0,'Tag','subplot_tess_window');
    if ~isempty(SlidesWin)
        TessWin = SlidesWin; % Subplot or movie window
        axx = findobj(TessWin,'type','axes');
        set(axx,'ClimMode','auto')    
    end
    
else
    if get(hTessWin.FreezeColormap,'Value')
        data = get(Cortex,'FaceVertexCData');
        FreezeTime = get(hTessWin.FreezeColormap,'Userdata');
        M = max(abs(Results.ImageGridAmp(:,FreezeTime)));
    else
        M = max(abs(Results.ImageGridAmp(:)));
    end
    
    if get(hSelect.AbsoluteCurrent,'Value')
        set(findobj(TessWin,'type','axes'),'ClimMode','manual','Clim',[0 M])
    else
        set(findobj(TessWin,'type','axes'),'ClimMode','manual','Clim',[-M M]) 
    end
end

if isempty(get(hSelect.ColorMAP,'Userdata'))
    %if ~isfield(cortex,'CDepth') & ~get(hSelect.MapCurvature,'Value')% No curvature mapping
    if ~get(hSelect.MapCurvature,'Value')% No curvature mapping
        if ~get(hSelect.AbsoluteCurrent,'Value')
            C = [flipud(fliplr(grayish(hot(128),.3)));grayish(hot(128),.3)];
        else
            %C = [(grayish(hot(128),.3))];
            load bst_cactivcmap
            C = map;
        end
    else
        [FVC,C]=catci(get(Cortex,'FaceVertexCData'),cortex.CDepth,M);
        set(Cortex,'FaceVertexCData',FVC);
    end
    
else % Truncated Colormap
    if ~isfield(cortex,'CDepth') % No curvature mapping
        C = get(hSelect.ColorMAP,'Userdata'); % Truncated Colormap
    else
        if get(hSelect.ZScoreThresholdApply,'Value') & get(hSelect.ZScore,'Value') 
            cThres = 100*get(hSelect.ZScoreThreshold,'Userdata');
        else
            cThres = str2num(get(hTessWin.TruncateFactor,'String'));
        end
        
        [FVC,C]=catci(get(Cortex,'FaceVertexCData'),cortex.CDepth,M,cThres);
        if get(hSelect.OrthoViews,'Value') == 0
            set(Cortex,'FaceVertexCData',FVC);
        else
            Cortex = findobj(TessWin,'Type','patch','Tag',ActiveTess,'Visible','on'); %CHEAT - need to be improved ?    
            set(Cortex,'FaceVertexCData',FVC);
        end
        
    end
end

colormap(C)

% % ------------------------
% case 'mapping_create' % Generates SLIDES and/or MOVIES
%    
%    Visu = get(DATAPLOT,'Userdata');
%    MEG = get(findobj(DATAPLOT,'Tag','MEG'),'value');
%    EEG = get(findobj(DATAPLOT,'Tag','EEG'),'value');
%    OTHER = get(findobj(DATAPLOT,'Tag','OTHER'),'value');
%    current = find([MEG,EEG,OTHER] == 1);
%    
%    mapping_win = findobj(0,'Tag','mapping');
%    if isempty(mapping_win)
%       open mapping.fig; 
%       mapping_win = findobj(0,'Tag','mapping');
%    else
%       return
%    end
%    
%    RAW = findobj(mapping_win ,'Tag','raw');
%    GRADIENT = findobj(mapping_win ,'Tag','gradient');
%    MAGNETIC = findobj(mapping_win,'Tag','magnetic');
%    
%    if current == 1 % If MEG
%       set([RAW,GRADIENT,MAGNETIC],'enable','on')
%       set(RAW,'Value',1)
%    else
%       set([RAW,GRADIENT,MAGNETIC],'enable','off')
%    end
%    
%    start = findobj(mapping_win,'Tag','start');
%    step = findobj(mapping_win,'Tag','step');
%    stop = findobj(mapping_win,'Tag','stop');
%    current_time = findobj(mapping_win,'Tag','current_time');
%    hDataplot.TimeSlider = findobj(mapping_win,'Tag','hDataplot.TimeSlider');
%    go = findobj(mapping_win,'Tag','go');
%    SPHERE = findobj(mapping_win,'Tag','sphere');
%    fit = findobj(mapping_win,'Tag','fit');
%    ROTATE = findobj(mapping_win,'Tag','rotate');
%    fit = findobj(mapping_win,'Tag','fit');
%    slides = findobj(mapping_win,'Tag','slides');
%    single = findobj(mapping_win,'Tag','single');
%    MOVIE = findobj(mapping_win,'Tag','movie');
%    
%    handles = guihandles(mapping_win);
%    
%    set(mapping_win,'Name',[get(mapping_win,'Name'), ' - ', modality{current}])
%    
%    set(start,'String', num2str(1000*Visu.Data.Time(1)))
%    set(stop,'String', num2str(1000*Visu.Data.Time(end)))
%    if length(Visu.Data.Time) > 1
%       set(step,'String', num2str(1000*(Visu.Data.Time(2)-Visu.Data.Time(1)),2))
%    else
%       set(step,'String', 0)
%    end
%    
%    
%    data.TIME_MIN = str2num(get(findobj(DATAPLOT,'Tag','time_min'),'String'));
%    data.TIME_MAX = str2num(get(findobj(DATAPLOT,'Tag','time_max'),'String'));
%    
%    if data.TIME_MAX > Visu.Data.Time(end) * 1000
%       data.TIME_MAX = Visu.Data.Time(end) * 1000;
%       set(findobj(DATAPLOT,'Tag','time_max'),'String',num2str(data.TIME_MAX,5))
%    end
%    
%    if data.TIME_MIN < Visu.Data.Time(1) * 1000
%       data.TIME_MIN = Visu.Data.Time(1) * 1000;
%       set(findobj(DATAPLOT,'Tag','time_min'),'String',num2str(data.TIME_MIN,5))
%    end
%    
%    if ~isempty(mapping_win)
%       hDataplot.TimeSlider = findobj(mapping_win,'Tag','hDataplot.TimeSlider');
%       if length(Visu.Data.Time) > 1
%          set(hDataplot.TimeSlider,'enable','on')
%          set(hDataplot.TimeSlider,'Min',data.TIME_MIN,'Max',data.TIME_MAX)
%       else
%          set(hDataplot.TimeSlider,'Min',data.TIME_MIN,'Max',2*data.TIME_MAX)
%          set(hDataplot.TimeSlider,'enable','off')
%       end
%       
%    end
%    set(hDataplot.TimeSlider,'Value',data.TIME_MIN,'enable','on')
%    
%    set(current_time,'String',num2str(get(hDataplot.TimeSlider,'Value'),4))
%    set(single,'Value',1)
%    mutincomp([single,slides,MOVIE])
%    set(SPHERE,'Value',1)
%    mutincomp([SPHERE,fit,handles.scalp])
%    
%    % Draw time cursor
%    figsingle = findobj(0,'Tag',['waves_single',int2str(current)]); % Figure of overlaping plots for the current modality
%    if ~isempty(figsingle) 
%       figure(figsingle)
%       ctime = get(hDataplot.TimeSlider,'Value');
%       cursor = line([ctime ctime],get(gca,'Ylim'));
%       set(cursor,'color','r','linewidth',3,'Tag','cursor','erasemode','Xor')
%    end
%    
%    fignplot = findobj(0,'Tag','nplot_win'); % superimposed n-plots
%    if ~isempty(fignplot) 
%       figure(fignplot)
%       ctime = get(hDataplot.TimeSlider,'Value');
%       cursor = line([ctime ctime],get(gca,'Ylim'));
%       set(cursor,'color','r','linewidth',3,'Tag','cursor','erasemode','Xor')
%    end
%    %-----------------------------------------------------------------------------------------------   



function [mixedRGB] = mix_curvature_current(varargin)


% Various colormaps ---------------------------------------
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);
% Various colormaps ------------------DONE----------------

[hf,hs,hl] = view_surface('test',Faces,tmpVertices,Results.ImageGridAmp(:,1));
delete(findobj(gcf,'type','light'))
%view(-10,0), 
material dull, camlight
set(hs,'backfacelighting','unlit')
hold on
%set(hf,'color','k')

colormap flag; colormenu
flag = get(gcf,'colormap'); 
flag = flag(:,[3 2 1]);
colormap(flag)

timee = 1:timestep:size(Results.ImageGridAmp,2)-corrsize/2;%
corr = zeros(size(Results.ImageGridAmp(patchh,:),1),length(timee));
tt = 0;
tmp = Results.ImageGridAmp(patchh,timee);
corr = abs(tmp)/max(tmp(:));%abs(activ(:,1+corrsize/2));

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);
    
    drawnow
end
