% sig_topo() - Plots scalp topographies of EEG efffects that have been
% tested for significance in a family of tests. Only works for
% EEG effects derived from averaging time points within
% time windows (for tests of every single time point within
% time windows use gui_erp.m or gui_pow.m). Topographies can
% be visualized in microvolts (for ERPs), decibels (for spectral
% power), or t-scores.
%
% Usage:
% >> sig_topo(GND_GRP_specGND_or_fname,test_id,varargin)
%
% Required Inputs:
% GND_GRP_specGND_or_fname - A GND, GRP, or specGND structure variable or the
% filename of such a variable that has
% been saved to disk. If you specifiy a filename be
% sure to include the file's path, unless the file is
% in the current working directory.
% test_id - [integer] The index # of the t-test results
% stored in the GND/GRP/specGND variable that you wish to visualize.
% To see what test results are available, look at
% the "t_tests" field of your variable (e.g., GND.t_tests)
% or use the functions headinfo.m or headinfo_spec.m.
%
% Optional Inputs:
% fig_id - [integer] The index # of the MATLAB figure in which
% the diagram will be produced. Useful for overwriting
% old figures. {default: lowest unused index}
% units - ['uV','dB', or 't'] If 'uV' topographies will be
% visualized in units of microvolts. If 'dB' spectral
% power will be shown in decibels (10*log10((uV^2)/Hz)).
% If 't', they will be shown in t-scores. {default: 't'}
% title_on - [0 or 1] If 1, the bin # and bin descriptor will be
% printed at the top of the figure. {default: 1}
% one_scale - [0 or 1] If 1, all topographies will be plot using
% the same color scale (useful for comparing effect
% sizes across topographies). If 0, each topography
% will be plot using a color scale tailored to it
% (provides maximal resolution of topography
% distribution). {default: 1 if plotting in uV or dB,
% 0 if plotting t-scores}
% scale_limits - [min_value max_value] Minimum and maximum value of
% color scale used to illustrate topographies (e.g.,
% [-10 10]). If specified, one_scale parameter is set
% to 1. If not specified, color scales are set to +/-
% the maximum absolute value of the topography.
% {default: not used}
% verblevel - An integer specifiying the amount of information you want
% this function to provide about what it is doing during runtime.
% Options are:
% 0 - quiet, only show errors, warnings, and EEGLAB reports
% 1 - stuff anyone should probably know
% 2 - stuff you should know the first time you start working
% with a data set {default value}
% 3 - stuff that might help you debug (show all
% reports)
%
% Global Variable:
% VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells
% functions how much to report about what they're doing during
% runtime) set by the optional function argument 'verblevel'
%
% Notes:
% -The printed/exported figure will have the same dimensions as the figure
% on the display. Thus you can undo figure clutter by re-sizing it.
%
% -Note you need data at least three channels to interpolate a scalp topography.
%
% Author:
% David Groppe
% Kutaslab, 3/2010
%%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
% 5/6/2010-Made compatible with FDR control
%
% 12/14/2010-Compatible with specGND variables now too
%
% 9/10/2012-scale_limits option added
%%%%%%%%%%%%%%%% FUTURE WORK %%%%%%%%%%%%%%%%%
% -
function sig_topo(GND_GRP_specGND_or_fname,test_id,varargin)
p=inputParser;
p.addRequired('GND_GRP_specGND_or_fname',@(x) ischar(x) || isstruct(x));
p.addRequired('test_id',@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('units','t',@ischar);
p.addParamValue('fig_id',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('title_on',1,@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('one_scale',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1));
p.addParamValue('scale_limits',[],@(x) isnumeric(x) && (length(x)==2));
p.parse(GND_GRP_specGND_or_fname,test_id,varargin{:});
% Manage VERBLEVEL
if isempty(p.Results.verblevel),
VERBLEVEL=2; %not global, just local
else
VERBLEVEL=p.Results.verblevel;
end
if ~(strcmpi(p.Results.units,'t') || strcmpi(p.Results.units,'uV') || strcmpi(p.Results.units,'dB'))
error('Value of optional input ''units'' must be ''t'', ''uV'', or ''dB''.');
end
%Load GND, GRP, or specGND struct
if ischar(GND_GRP_specGND_or_fname),
fprintf('Loading GND, GRP, or specGND struct variable from file %s.\n',GND_GRP_specGND_or_fname);
load(GND_GRP_specGND_or_fname,'-MAT');
if ~exist('GND','var') && ~exist('GRP','var') && ~exist('specGND','var')
error('File %s does not contain a GND, GRP, or specGND variable.',GND_GRP_specGND_or_fname);
end
if exist('GRP','var'),
GND=GRP; %for simplicity GRP variables are re-named GND
clear GRP;
elseif exist('specGND','var'),
GND=specGND; %for simplicity specGND variables are re-named GND
clear specGND;
end
VerbReport(sprintf('Experiment: %s',GND.exp_desc),2,VERBLEVEL);
else
GND=GND_GRP_specGND_or_fname; %for simplicity GRP and specGND variables are re-named GND
clear GND_GRP_specGND_or_fname;
end
n_test=length(GND.t_tests);
if test_id>n_test,
error('There are only %d t-tests stored with these data, but you requested Test #%d',n_test,test_id);
end
if isfield(GND.t_tests(1),'freq_band')
% We're dealing with spectral power analysis. Create temporary fields
% to make frequency domain plotting compatible easily compatible with time domain plot
GND.t_tests(test_id).time_wind=GND.t_tests(test_id).freq_band;
GND.t_tests(test_id).mean_wind=GND.t_tests(test_id).mean_band;
GND.t_tests(test_id).used_tpt_ids=GND.t_tests(test_id).used_freq_ids;
GND.grands_t=GND.grands_pow_dB_t;
GND.grands=GND.grands_pow_dB;
freq_step=GND.freqs(2)-GND.freqs(1);
freq_ord=orderofmag(freq_step);
ord_pow=log10(freq_ord); %useful for formatting topo titles
if ord_pow>=0,
n_dig_past_dot=0;
else
n_dig_past_dot=-ord_pow;
end
GND.time_pts=rnd_orderofmag(GND.freqs);
freq_domain=1;
if strcmpi(p.Results.units,'uV'),
watchit('You can''t plot spectral power topographies in units of microvolts. They will be plot in decibels (i.e., 10*log10((uV^2)/Hz))) instead.');
end
else
freq_domain=0;
if strcmpi(p.Results.units,'dB'),
watchit('You can''t plot ERP topographies in units of decibels. They will be plot in microvolts instead.');
end
end
if ~strcmpi(GND.t_tests(test_id).mean_wind,'yes'),
error('sig_topo.m only works for significance tests applied to data averaged across multiple time points.');
end
n_topos=length(GND.t_tests(test_id).used_tpt_ids);
wdth=floor(sqrt(n_topos));
ht=ceil(n_topos/wdth);
bin=GND.t_tests(test_id).bin;
chans=GND.t_tests(test_id).used_chan_ids;
if length(chans)<3,
watchit('sig_topo.m cannot plot scalp topographies with less than 3 electrodes.');
return
end
if isempty(p.Results.fig_id)
fig_h=figure;
else
fig_h=figure(p.Results.fig_id); clf;
end
if isfield(GND,'bin_info'),
set(fig_h,'name',['Bin ' int2str(bin) ' [' GND.bin_info(bin).bindesc ']'],'paperpositionmode','auto');
else
set(fig_h,'name',['Bin ' int2str(bin) ' [' GND.bindesc{bin} ']'],'paperpositionmode','auto');
end
%setting paperpositionmode to 'auto' means that if the figure is manually
%resized, the printed version of the figure will reflect whatever the
%shown size was (at the time of printing)
if strcmpi(p.Results.units,'uV') || strcmpi(p.Results.units,'dB') ,
if freq_domain,
units='dB';
else
units='\muV';
end
%Compute mean voltage or spectral power in each time window
mns=zeros(length(chans),n_topos);
for w=1:n_topos,
mns(:,w)=mean(GND.grands(chans,GND.t_tests(test_id).used_tpt_ids{w},bin),2);
end
if ~isempty(p.Results.scale_limits)
one_scale=1;
maplimits=p.Results.scale_limits;
VerbReport(sprintf('Using color scale limits of: %f to %f',maplimits(1),maplimits(2)),2,VERBLEVEL);
VerbReport(sprintf('All topographies will have the same color scale.'),2,VERBLEVEL);
elseif isempty(p.Results.one_scale) || (p.Results.one_scale<=0),
one_scale=0;
maplimits='absmax';
else
one_scale=1;
mx=max(max(abs(mns)));
maplimits=[-1 1]*mx;
end
else
units='t';
if ~isempty(p.Results.scale_limits)
one_scale=1;
maplimits=p.Results.scale_limits;
VerbReport(sprintf('Using color scale limits of: %f to %f',maplimits(1),maplimits(2)),2,VERBLEVEL);
VerbReport(sprintf('All topographies will have the same color scale.'),2,VERBLEVEL);
elseif isempty(p.Results.one_scale) || (p.Results.one_scale>0),
one_scale=1;
mx=max(max(abs(GND.t_tests(test_id).data_t)));
maplimits=[-1 1]*mx;
else
one_scale=0;
maplimits='absmax';
end
end
if VERBLEVEL,
if isnan(GND.t_tests(test_id).n_perm),
fprintf('q level of critical t-scores: %f., FDR method: %s\n', ...
GND.t_tests(test_id).desired_alphaORq,GND.t_tests(test_id).mult_comp_method);
else
fprintf('Estimated alpha level of critical t-scores: %f.\n', ...
GND.t_tests(test_id).estimated_alpha);
end
end
for a=1:n_topos,
subplot(ht,wdth,a);
time_wind=GND.t_tests(test_id).time_wind(a,:);
if strcmpi(p.Results.units,'uV') || strcmpi(p.Results.units,'dB') ,
mn=mns(:,a);
else
mn=GND.t_tests(test_id).data_t(:,a);
end
if isnan(GND.t_tests(test_id).n_perm),
sig_chans=find(GND.t_tests(test_id).fdr_rej(:,a));
else
sig_chans=find(GND.t_tests(test_id).adj_pval(:,a)<GND.t_tests(test_id).estimated_alpha);
end
topoplotMK(mn,GND.chanlocs(chans),'emarker2',{sig_chans,'o',[1 1 1],4}, ...
'maplimits',maplimits);
if freq_domain,
lab_form=['%.' int2str(n_dig_past_dot) 'f-%.' int2str(n_dig_past_dot) 'f Hz'];
else
lab_form='%d-%d ms';
end
hh=title(sprintf(lab_form,time_wind(1),time_wind(2)));
if ~one_scale,
hcb=colorbar;
hy=ylabel(hcb,units);
set(hy,'rotation',0,'verticalalignment','middle','position',[4.9+n_topos*.2 0.03 1.0001]);
end
end
if one_scale,
hcb=cbar_nudge('vert',0,maplimits);
hy=ylabel(hcb,units);
set(hy,'rotation',0,'verticalalignment','middle');
end
% Plot title of entire figure
if p.Results.title_on,
if isfield(GND,'bin_info'),
h=textsc(['Bin ' int2str(bin) ': ' GND.bin_info(bin).bindesc],'title');
else
h=textsc(['Bin ' int2str(bin) ': ' GND.bindesc{bin}],'title');
end
set(h,'fontsize',14,'fontweight','bold');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% cbar_nudge (plots colorbar slightly to left of cbar.m
% search for "DG change" to find the line I edited
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [handle]=cbar_nudge(arg,colors,minmax, grad)
if nargin < 2
colors = 0;
end
posscale = 'off';
if nargin < 1
arg = 'vert';
ax = [];
else
if isempty(arg)
arg = 0;
end
if arg(1) == 0
ax = [];
arg = 'vert';
elseif strcmpi(arg, 'pos')
ax = [];
arg = 'vert';
posscale = 'on';
else
if isstr(arg)
ax = [];
else
ax = arg;
arg = [];
end
end
end
if nargin>2
if size(minmax,1) ~= 1 | size(minmax,2) ~= 2
help cbar
fprintf('cbar() : minmax arg must be [min,max]\n');
return
end
end
if nargin < 4
grad = 5;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Choose colorbar position
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (length(colors) == 1) & (colors == 0)
t = caxis;
else
t = [0 1];
end
if ~isempty(arg)
if strcmp(arg,'vert')
cax = gca;
pos = get(cax,'Position');
stripe = 0.04;
edge = 0.01;
space = -.02;
set(cax,'Position',[pos(1) pos(2) pos(3) pos(4)])
rect = [pos(1)+pos(3)+space pos(2) stripe*pos(3) pos(4)];
ax = axes('Position', rect);
elseif strcmp(arg,'horiz')
cax = gca;
pos = get(cax,'Position');
stripe = 0.075;
space = .1;
set(cax,'Position',...
[pos(1) pos(2)+(stripe+space)*pos(4) pos(3) (1-stripe-space)*pos(4)])
rect = [pos(1) pos(2) pos(3) stripe*pos(4)];
ax = axes('Position', rect);
end
else
pos = get(ax,'Position');
if pos(3) > pos(4)
arg = 'horiz';
else
arg = 'vert';
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw colorbar using image()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
map = colormap;
n = size(map,1);
if length(colors) == 1
if strcmp(arg,'vert')
if strcmpi(posscale, 'on')
image([0 1],[0 t(2)],[ceil(n/2):n-colors]');
else
image([0 1],t,[1:n-colors]');
end;
set(ax,'xticklabelmode','manual')
set(ax,'xticklabel',[],'YAxisLocation','right')
else
image(t,[0 1],[1:n-colors]);
set(ax,'yticklabelmode','manual')
set(ax,'yticklabel',[],'YAxisLocation','right')
end
set(ax,'Ydir','normal','YAxisLocation','right')
else % length > 1
if max(colors) > n
error('Color vector excedes size of colormap')
end
if strcmp(arg,'vert')
image([0 1],t,[colors]');
set(ax,'xticklabelmode','manual')
set(ax,'xticklabel',[])
else
image([0 1],t,[colors]);
set(ax,'yticklabelmode','manual')
set(ax,'yticklabel',[],'YAxisLocation','right')
end
set(ax,'Ydir','normal','YAxisLocation','right')
end
%%%%%%%%%%%%%%%%%%%%%%%%%
% Adjust cbar ticklabels
%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin > 2
if strcmp(arg,'vert')
Cax = get(ax,'Ylim');
else
Cax = get(ax,'Xlim');
end;
CBTicks = [Cax(1):(Cax(2)-Cax(1))/(grad-1):Cax(2)]; % caxis tick positions
CBLabels = [minmax(1):(minmax(2)-minmax(1))/(grad-1):minmax(2)]; % tick labels
dec = floor(log10(max(abs(minmax)))); % decade of largest abs value
CBLabels = ([minmax]* [ linspace(1,0, grad);linspace(0, 1, grad)]);
if dec<1
CBLabels = round(CBLabels*10^(1-dec))*10^(dec-1);
elseif dec == 1
CBLabels = round(CBLabels*10^(2-dec))*10^(dec-2);
else
CBLabels = round(CBLabels);
end
if strcmp(arg,'vert')
set(ax,'Ytick',CBTicks);
set(ax,'Yticklabel',CBLabels);
else
set(ax,'Xtick',CBTicks);
set(ax,'Xticklabel',CBLabels);
end
end
handle = ax;
%%%%%%%%%%%%%%%%%%
% Adjust cbar tag
%%%%%%%%%%%%%%%%%%
set(ax,'tag','cbar')
