tesa_plot.m

% tesa_plot() - plots TMS-evoked activity averaged over trials. Timing of
%               the TMS pulse is indicated with a red dashed line. Output
%               from peak analyses can also be plotted as well as
%               confidence intervals across trials.
%
% Usage:
%   >>  tesa_plot( EEG );
%   >>  tesa_plot( EEG, 'key1',value1... );
%
% Inputs:
%   EEG             - EEGLAB EEG structure
% 
% Optional input pairs:
%   'xlim', [min,max] - integers describing the x axis limits in ms.
%                   default = [-100,500]
%   'ylim', [min,max] - integers describing the y axis limits in ms. If
%                   left blank, the plot will automatically scale.
%                   default = []
%                   Examples: [-10,10]
%   'elec','str'  - string describing a single electrode  
%                   for plotting. If left blank, a butterfly plot of all electrodes 
%                   will be plotted.
%                   default = [];
%                   Example: 'Cz'
%   'CI','str'      - 'on' | 'off'. Plot confidence interval calculated
%                   across trials. This option is not available for GMFA
%                   and butterfly plots.
%                   default = 'off';
%   
% Options for plotting output from tesa_tepextract and tesa_peakanalysis
%   'tepType','str'  - 'data' | 'ROI' | 'GMFA'. 'Data' input extracts data
%                   for plotting from EEG.data. 'ROI' input selects a ROI
%                   generated by tesa_tepextract. 'GMFA' input selects a
%                   GMFA generated by tesa_tepextract. Note that if
%                   multiple ROIs or GMFAs are present, 'tepName' must also
%                   be included to determine which one to plot.
%                   Default = 'data';
%   'tepName', 'str' - String is either the name of a ROI or GMFA generated
%                   by tesa_tepextract. The default names generated by
%                   tesa_tepextract are 'R1','R2'... etc. however the user
%                   can also define names. This is required if multiple ROIs or GMFAs are present. 
%                   Default = []; 
%   'plotPeak','str' - 'off' | 'on'. If 'on', peaks identified with
%                   tesa_peakanalysis will also be plotted. A box defining
%                   the lower and upper time limits in which the peak was
%                   searched will also be plotted. Green indiciates that a
%                   peak was found, red a peak was not detected. Note that
%                   tesa_peakanalysis must be run prior to using this
%                   function and 'tepType' must also be called.
% 
% Examples:
%   tesa_plot(EEG); % plot a butterfly plot.
%   tesa_plot(EEG, 'xlim', [-200,600], 'ylim', [-10,10], 'elec', 'P1', 'CI', 'on'); % plot electrode P1, rescale the axes and include confidence intervals.
%   tesa_plot(EEG, 'tepType','ROI','tepName','parietal','plotPeak','on','CI','on'); % plot the output from a ROI analysis called parietal including detected peaks and confidence intervals
%   tesa_plot(EEG, 'tepType','GMFA','plotPeak','on'); % plot output from a GMFA analysis including detected peaks
% 
% See also:
%   tesa_tepextract, tesa_peakanalysis, tesa_peakoutput 

% Copyright (C) 2016  Nigel Rogasch, Monash University,
% nigel.rogasch@monash.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function tesa_plot( EEG, varargin )

if nargin < 1
	error('Not enough input arguments.');
end

%define defaults
options = struct('xlim',[-100,500],'ylim',[],'elec',[],'CI','off','tepType','data','tepName',[],'plotPeak','off');

% read the acceptable names
optionNames = fieldnames(options);

% count arguments
nArgs = length(varargin);
if round(nArgs/2)~=nArgs/2
   error('EXAMPLE needs key/value pairs')
end

for pair = reshape(varargin,2,[]) % pair is {propName;propValue}
   inpName = pair{1}; % make case insensitive

   if any(strcmpi(inpName,optionNames))%looks for known options and replaces these in options
      options.(inpName) = pair{2};
   else
      error('%s is not a recognized parameter name',inpName)
   end
end

%Check xlim input
if size(options.xlim,2)~=2
    error('Input for ''xlim'' must be in the following format [min,max] e.g. [-100,500].')
end

%Check xlim is within range of data
if options.xlim(1,1) < EEG.times(1,1) || options.xlim(1,1) > EEG.times(1,end) || options.xlim(1,2) < EEG.times(1,1) || options.xlim(1,2) > EEG.times(1,end)
    error('Input for ''xlim'' is outside of the range of the data (%d to %d).', EEG.times(1,1), EEG.times(1,end));
end

%Check ylim input
if ~isempty(options.ylim)
    if size(options.ylim,2)~=2
        error('Input for ''ylim'' must be in the following format [min,max] e.g. [-10,10].')
    end
end

%Check CI input
if ~(strcmp(options.CI,'off') || strcmp(options.CI,'on'))
    error('Input for ''CI'' must be either ''on'' or ''off''.');
end

%Check tepType input
if ~(strcmp(options.tepType,'data') || strcmp(options.tepType,'ROI') || strcmp(options.tepType,'GMFA'))
    error('Input for ''tepType'' must be either ''data'', ''ROI'' or ''GMFA''.');
end    

%If tepType is ROI, check whether tepName is needed
if strcmp(options.tepType,'ROI')
    if ~isfield(EEG,'ROI')
        error('There are no ROI analyses present in the data. Please run tesa_tepextract.')
    elseif isempty(options.tepName) && size(fieldnames(EEG.ROI),1) > 1
        error('There are multiple ROIs present in the data. Please enter a specific ROI using ''tepName'', ''str'' where str is the name of the specific ROI.')
    end
end

%If tepType is GMFA, check whether tepName is needed
if strcmp(options.tepType,'GMFA')
    if ~isfield(EEG,'GMFA')
        error('There are no GMFA analyses present in the data. Please run tesa_tepextract.')
    elseif isempty(options.tepName) && size(fieldnames(EEG.GMFA),1) > 1
        error('There are multiple GMFAs present in the data. Please enter a specific GMFA using ''tepName'', ''str'' where str is the name of the specific GMFA.')
    end
end

%Check tepName input exists
if ~isempty(options.tepName)
    if ~(strcmp(options.tepType,'ROI') || strcmp(options.tepType,'GMFA'))
        error('Please indicate which type of TEP you would like to perform analysis on using ''tepType'', ''str'', where str is either ''ROI'' or ''GMFA''.');
    elseif strcmp(options.tepType,'ROI')
        if ~isfield(EEG,'ROI')
            error('There are no ROI analyses present in the data. Please run tesa_tepextract.')
        elseif ~isfield(EEG.ROI,options.tepName)
            error('''tepName'' ''%s'' does not exist for tepType ROI. Please revise.',options.tepName);
        end
    elseif strcmp(options.tepType,'GMFA')
        if ~isfield(EEG,'GMFA')
            error('There is no GMFA analysis present in the data. Please run tesa_tepextract.')
        elseif ~isfield(EEG.GMFA,options.tepName)
            error('''tepName'' ''%s'' does not exist for tepType GMFA. Please revise.',options.tepName);
        end
    end
end

%Check plotPeak input
if ~(strcmp(options.plotPeak,'off') || strcmp(options.plotPeak,'on'))
    error('Input for ''plotPeak'' must be either ''on'' or ''off''.');
end

t = figure;

%Plot figure
if strcmpi(options.tepType,'data')
    if isempty(options.elec)
        plot(EEG.times,mean(EEG.data,3),'b'); hold on;
    elseif ~isempty(options.elec);
        for z = 1:EEG.nbchan;
            chan{1,z} = EEG.chanlocs(1,z).labels;
        end;
        elecNum = find(strcmpi(options.elec,chan));
        plot(EEG.times,mean(EEG.data(elecNum,:,:),3),'b'); hold on;
    end
elseif strcmpi(options.tepType,'ROI')
    if isempty(options.tepName)
        tempN = fieldnames(EEG.ROI);
        tempName = tempN{1,1};
    else
        tempName = options.tepName;
    end
    plot(EEG.ROI.(tempName).time,EEG.ROI.(tempName).tseries,'b'); hold on;
elseif strcmpi(options.tepType,'GMFA')
    if isempty(options.tepName)
        tempN = fieldnames(EEG.GMFA);
        tempName = tempN{1,1};
    else
        tempName = options.tepName;
    end
    plot(EEG.GMFA.(tempName).time,EEG.GMFA.(tempName).tseries,'b'); hold on;
end

%Figure settings
if isempty(options.ylim)
    set(gca,'xlim',options.xlim,'box','off','tickdir','out')
elseif ~isempty(options.ylim)
    set(gca,'xlim',options.xlim,'ylim',options.ylim,'box','off','tickdir','out')
end

if strcmpi(options.tepType,'GMFA')
    ylabel('GMFA (\muV)');
else
    ylabel('Amplitude (\muV)');
end
xlabel('Time (ms)');

if strcmp(options.tepType,'data') && isempty(options.elec)
    title('All electrodes');
elseif strcmp(options.tepType,'data') && ~isempty(options.elec)
    title(options.elec);
elseif strcmp(options.tepType,'ROI')
    if isempty(options.tepName)
        tempName = fieldnames(EEG.ROI);
        titleIn = ['Region of interest (', tempName{1,1}, ')'];
    else
        titleIn = ['Region of interest (', options.tepName, ')'];
    end
    title(titleIn)
elseif strcmp(options.tepType,'GMFA')
    if isempty(options.tepName)
        tempName = fieldnames(EEG.GMFA);
        titleIn = ['Global mean field amplitude (', tempName{1,1}, ')'];
    else
        titleIn = ['Global mean field amplitude (', options.tepName, ')'];
    end
    title(titleIn)
end

%Confidence intervals (if requested)
if strcmp(options.CI,'on')
    if strcmpi(options.tepType,'data')
        if isempty(options.elec)
            error('Confidence intervals are only available for single channel plots. Please use ''elec'',''str'' to enter channel name.')
        end
        M = mean(EEG.data(elecNum,:,:),3); 
        CI = 1.96*(std(EEG.data(elecNum,:,:),0,3)./(sqrt(size(EEG.data,3))));
        f = fill([EEG.times,fliplr(EEG.times)],[M-CI,fliplr(M+CI)],'b');
        set(f,'FaceAlpha',0.3);set(f,'EdgeColor', 'none');
    elseif strcmpi(options.tepType,'ROI')
        if isempty(options.tepName)
            tempN = fieldnames(EEG.ROI);
            tempName = tempN{1,1};
        else
            tempName = options.tepName;
        end
        M = EEG.ROI.(tempName).tseries; 
        CI = EEG.ROI.(tempName).CI;
        f = fill([EEG.times,fliplr(EEG.times)],[M-CI,fliplr(M+CI)],'b');
        set(f,'FaceAlpha',0.3);set(f,'EdgeColor', 'none');
    elseif strcmpi(options.tepType,'GMFA')
        error('Confidence intervals are not available for GMFA');
    end
end

%Plot peaks found
if strcmpi(options.plotPeak,'on')
    if strcmpi(options.tepType,'data')
        error('Peaks are unavailable for butterfly plots and single electrodes. Please plot ROI or GMFA data.');
    else
        if isempty(options.tepName)
            tempN = fieldnames(EEG.(options.tepType));
            tempName = tempN{1,1};
        else
            tempName = options.tepName;
        end
        peaksTemp = fieldnames(EEG.(options.tepType).(tempName));
        peakLog = strncmp('P',peaksTemp,1) | strncmp('N',peaksTemp,1);
        if sum(peakLog)==0
            error('Peaks are unavailable. Please run pop_tesa_peakanalysis first.');
        end
        peak = peaksTemp(peakLog)';
    end

    Ysize = get(gca,'ylim');
    Yheight = (Ysize(1,2)-Ysize(1,1))./10;
    for z=1:size(peak,2);
        if strcmp(options.tepType,'ROI')
            if strcmp(EEG.(options.tepType).(tempName).(peak{1,z}).found,'yes')
                tempColour = 'g';
                tempLat = EEG.(options.tepType).(tempName).(peak{1,z}).lat;
            elseif strcmp(EEG.(options.tepType).(tempName).(peak{1,z}).found,'no')
                tempColour = 'r';
                tempLat = EEG.(options.tepType).(tempName).(peak{1,z}).peak;
            end
        elseif strcmp(options.tepType,'GMFA')
            if strcmp(EEG.(options.tepType).(tempName).(peak{1,z}).found,'yes')
                tempColour = 'g';
                tempLat = EEG.(options.tepType).(tempName).(peak{1,z}).lat;
            elseif strcmp(EEG.(options.tepType).(tempName).(peak{1,z}).found,'no')
                tempColour = 'r';
                tempLat = EEG.(options.tepType).(tempName).(peak{1,z}).peak;
            end
        end        

        %Plot cross at peak
        plot(tempLat,EEG.(options.tepType).(tempName).(peak{1,z}).amp,'+','MarkerEdgeColor',tempColour,'MarkerSize',20,'LineWidth',2);%hold on;
        
        %Plot box
        rectangle('position',[EEG.(options.tepType).(tempName).(peak{1,z}).minWin, EEG.(options.tepType).(tempName).(peak{1,z}).amp-(Yheight./2),EEG.(options.tepType).(tempName).(peak{1,z}).maxWin-EEG.(options.tepType).(tempName).(peak{1,z}).minWin,Yheight], 'LineStyle','--','EdgeColor',tempColour);
        
        %Plot text
        if strncmp(peak{1,z},'P',1)
            text(double(tempLat),double(EEG.(options.tepType).(tempName).(peak{1,z}).amp + Yheight),peak{1,z});
        elseif strncmp(peak{1,z},'N',1)
            text(double(tempLat),double(EEG.(options.tepType).(tempName).(peak{1,z}).amp - Yheight),peak{1,z});
        end
%         centreLine = EEG.(options.tepType).(tempName).(peak{1,z}).amp-(Yheight/2):0.01:EEG.(options.tepType).(tempName).(peak{1,z}).amp+(Yheight/2);
%         centreLat = ones(1,size(centreLine,2))*EEG.(options.tepType).(tempName).(peak{1,z}).peak;
%         plot(centreLat,centreLine,'LineStyle','--','Color',tempColour);
    end
end

%Plot timing of TMS pulse
plot([0 0], get(gca,'ylim'),'r--');

fprintf('TEP plot generated. \n');

end