sig_topo

 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

0001 % sig_topo() - Plots scalp topographies of EEG efffects that have been
0002 %              tested for significance in a family of tests.  Only works for
0003 %              EEG effects derived from averaging time points within
0004 %              time windows (for tests of every single time point within
0005 %              time windows use gui_erp.m or gui_pow.m).  Topographies can
0006 %              be visualized in microvolts (for ERPs), decibels (for spectral
0007 %              power), or t-scores.
0008 %
0009 % Usage:
0010 %  >> sig_topo(GND_GRP_specGND_or_fname,test_id,varargin)
0011 %
0012 % Required Inputs:
0013 %  GND_GRP_specGND_or_fname - A GND, GRP, or specGND structure variable or the
0014 %                             filename of such a variable that has
0015 %                             been saved to disk. If you specifiy a filename be
0016 %                             sure to include the file's path, unless the file is
0017 %                             in the current working directory.
0018 %  test_id          - [integer] The index # of the t-test results
0019 %                     stored in the GND/GRP/specGND variable that you wish to visualize.
0020 %                     To see what test results are available, look at
0021 %                     the "t_tests" field of your variable (e.g., GND.t_tests)
0022 %                     or use the functions headinfo.m or headinfo_spec.m.
0023 %
0024 % Optional Inputs:
0025 %  fig_id           - [integer] The index # of the MATLAB figure in which
0026 %                     the diagram will be produced.  Useful for overwriting
0027 %                     old figures. {default: lowest unused index}
0028 %  units            - ['uV','dB', or 't'] If 'uV' topographies will be
0029 %                     visualized in units of microvolts.  If 'dB' spectral
0030 %                     power will be shown in decibels (10*log10((uV^2)/Hz)).
0031 %                     If 't', they will be shown in t-scores. {default: 't'}
0032 %  title_on         - [0 or 1] If 1, the bin # and bin descriptor will be
0033 %                     printed at the top of the figure. {default: 1}
0034 %  one_scale        - [0 or 1] If 1, all topographies will be plot using
0035 %                     the same color scale (useful for comparing effect
0036 %                     sizes across topographies).  If 0, each topography
0037 %                     will be plot using a color scale tailored to it
0038 %                     (provides maximal resolution of topography
0039 %                     distribution).  {default: 1 if plotting in uV or dB,
0040 %                     0 if plotting t-scores}
0041 %  scale_limits     - [min_value max_value] Minimum and maximum value of
0042 %                     color scale used to illustrate topographies (e.g.,
0043 %                     [-10 10]). If specified, one_scale parameter is set
0044 %                     to 1. If not specified, color scales are set to +/-
0045 %                     the maximum absolute value of the topography.
0046 %                     {default: not used}
0047 %  verblevel        - An integer specifiying the amount of information you want
0048 %                     this function to provide about what it is doing during runtime.
0049 %                     Options are:
0050 %                      0 - quiet, only show errors, warnings, and EEGLAB reports
0051 %                      1 - stuff anyone should probably know
0052 %                      2 - stuff you should know the first time you start working
0053 %                          with a data set {default value}
0054 %                      3 - stuff that might help you debug (show all
0055 %                          reports)
0056 %
0057 % Global Variable:
0058 %   VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells
0059 %               functions how much to report about what they're doing during
0060 %               runtime) set by the optional function argument 'verblevel'
0061 %
0062 % Notes:
0063 % -The printed/exported figure will have the same dimensions as the figure
0064 % on the display.  Thus you can undo figure clutter by re-sizing it.
0065 %
0066 % -Note you need data at least three channels to interpolate a scalp topography.
0067 %
0068 % Author:
0069 % David Groppe
0070 % Kutaslab, 3/2010
0071 
0072 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
0073 % 5/6/2010-Made compatible with FDR control
0074 %
0075 % 12/14/2010-Compatible with specGND variables now too
0076 %
0077 % 9/10/2012-scale_limits option added
0078 
0079 %%%%%%%%%%%%%%%% FUTURE WORK %%%%%%%%%%%%%%%%%
0080 % -
0081 
0082 function sig_topo(GND_GRP_specGND_or_fname,test_id,varargin)
0083 
0084 p=inputParser;
0085 p.addRequired('GND_GRP_specGND_or_fname',@(x) ischar(x) || isstruct(x));
0086 p.addRequired('test_id',@(x) isnumeric(x) && (length(x)==1));
0087 p.addParamValue('units','t',@ischar);
0088 p.addParamValue('fig_id',[],@(x) isnumeric(x) && (length(x)==1));
0089 p.addParamValue('title_on',1,@(x) isnumeric(x) && (length(x)==1));
0090 p.addParamValue('one_scale',[],@(x) isnumeric(x) && (length(x)==1));
0091 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1));
0092 p.addParamValue('scale_limits',[],@(x) isnumeric(x) && (length(x)==2));
0093 
0094 p.parse(GND_GRP_specGND_or_fname,test_id,varargin{:});
0095 
0096 % Manage VERBLEVEL
0097 if isempty(p.Results.verblevel),
0098     VERBLEVEL=2; %not global, just local
0099 else
0100     VERBLEVEL=p.Results.verblevel;
0101 end
0102 
0103 if ~(strcmpi(p.Results.units,'t') || strcmpi(p.Results.units,'uV') || strcmpi(p.Results.units,'dB'))
0104     error('Value of optional input ''units'' must be ''t'', ''uV'', or ''dB''.');
0105 end
0106     
0107 %Load GND, GRP, or specGND struct
0108 if ischar(GND_GRP_specGND_or_fname),
0109     fprintf('Loading GND, GRP, or specGND struct variable from file %s.\n',GND_GRP_specGND_or_fname);
0110     load(GND_GRP_specGND_or_fname,'-MAT');
0111     if ~exist('GND','var') && ~exist('GRP','var') && ~exist('specGND','var')
0112         error('File %s does not contain a GND, GRP, or specGND variable.',GND_GRP_specGND_or_fname);
0113     end
0114     if exist('GRP','var'),
0115         GND=GRP; %for simplicity GRP variables are re-named GND
0116         clear GRP;
0117     elseif exist('specGND','var'),
0118         GND=specGND; %for simplicity specGND variables are re-named GND
0119         clear specGND;
0120     end
0121     VerbReport(sprintf('Experiment: %s',GND.exp_desc),2,VERBLEVEL);
0122 else
0123     GND=GND_GRP_specGND_or_fname; %for simplicity GRP and specGND variables are re-named GND
0124     clear GND_GRP_specGND_or_fname;
0125 end
0126 
0127 n_test=length(GND.t_tests);
0128 if test_id>n_test,
0129    error('There are only %d t-tests stored with these data, but you requested Test #%d',n_test,test_id); 
0130 end
0131 
0132 if isfield(GND.t_tests(1),'freq_band')
0133     % We're dealing with spectral power analysis.  Create temporary fields
0134     % to make frequency domain plotting compatible easily compatible with time domain plot
0135     GND.t_tests(test_id).time_wind=GND.t_tests(test_id).freq_band;
0136     GND.t_tests(test_id).mean_wind=GND.t_tests(test_id).mean_band;
0137     GND.t_tests(test_id).used_tpt_ids=GND.t_tests(test_id).used_freq_ids;
0138     GND.grands_t=GND.grands_pow_dB_t;
0139     GND.grands=GND.grands_pow_dB;
0140     freq_step=GND.freqs(2)-GND.freqs(1);
0141     freq_ord=orderofmag(freq_step);
0142     ord_pow=log10(freq_ord); %useful for formatting topo titles
0143     if ord_pow>=0,
0144         n_dig_past_dot=0;
0145     else
0146         n_dig_past_dot=-ord_pow;
0147     end
0148     GND.time_pts=rnd_orderofmag(GND.freqs);
0149     freq_domain=1;
0150     if strcmpi(p.Results.units,'uV'),
0151        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.');
0152     end
0153 else
0154     freq_domain=0;
0155     if strcmpi(p.Results.units,'dB'),
0156         watchit('You can''t plot ERP topographies in units of decibels.  They will be plot in microvolts instead.');
0157     end
0158 end
0159 
0160 if ~strcmpi(GND.t_tests(test_id).mean_wind,'yes'),
0161    error('sig_topo.m only works for significance tests applied to data averaged across multiple time points.'); 
0162 end
0163 
0164 n_topos=length(GND.t_tests(test_id).used_tpt_ids);
0165 wdth=floor(sqrt(n_topos));
0166 ht=ceil(n_topos/wdth);
0167 
0168 bin=GND.t_tests(test_id).bin;
0169 chans=GND.t_tests(test_id).used_chan_ids;
0170 if length(chans)<3,
0171    watchit('sig_topo.m cannot plot scalp topographies with less than 3 electrodes.');
0172    return
0173 end
0174 if isempty(p.Results.fig_id)
0175     fig_h=figure;
0176 else
0177     fig_h=figure(p.Results.fig_id); clf;
0178 end
0179 if isfield(GND,'bin_info'),
0180     set(fig_h,'name',['Bin ' int2str(bin) ' [' GND.bin_info(bin).bindesc ']'],'paperpositionmode','auto');
0181 else
0182     set(fig_h,'name',['Bin ' int2str(bin) ' [' GND.bindesc{bin} ']'],'paperpositionmode','auto');
0183 end
0184 %setting paperpositionmode to 'auto' means that if the figure is manually
0185 %resized, the printed version of the figure will reflect whatever the
0186 %shown size was (at the time of printing)
0187 
0188 if strcmpi(p.Results.units,'uV') || strcmpi(p.Results.units,'dB') ,
0189     if freq_domain,
0190         units='dB';
0191     else
0192         units='\muV';
0193     end
0194     %Compute mean voltage or spectral power in each time window
0195     mns=zeros(length(chans),n_topos);
0196     for w=1:n_topos,
0197         mns(:,w)=mean(GND.grands(chans,GND.t_tests(test_id).used_tpt_ids{w},bin),2);
0198     end
0199     if ~isempty(p.Results.scale_limits)
0200         one_scale=1;
0201         maplimits=p.Results.scale_limits;
0202         VerbReport(sprintf('Using color scale limits of: %f to %f',maplimits(1),maplimits(2)),2,VERBLEVEL);
0203         VerbReport(sprintf('All topographies will have the same color scale.'),2,VERBLEVEL);
0204     elseif isempty(p.Results.one_scale) || (p.Results.one_scale<=0),
0205         one_scale=0;
0206         maplimits='absmax';
0207     else
0208         one_scale=1;
0209         mx=max(max(abs(mns)));
0210         maplimits=[-1 1]*mx;
0211     end
0212 else
0213     units='t';
0214     if ~isempty(p.Results.scale_limits)
0215         one_scale=1;
0216         maplimits=p.Results.scale_limits;
0217         VerbReport(sprintf('Using color scale limits of: %f to %f',maplimits(1),maplimits(2)),2,VERBLEVEL);
0218         VerbReport(sprintf('All topographies will have the same color scale.'),2,VERBLEVEL);
0219     elseif isempty(p.Results.one_scale) || (p.Results.one_scale>0),
0220         one_scale=1;
0221         mx=max(max(abs(GND.t_tests(test_id).data_t)));
0222         maplimits=[-1 1]*mx;
0223     else
0224         one_scale=0;
0225         maplimits='absmax';
0226     end
0227 end
0228 
0229 if VERBLEVEL,
0230     if isnan(GND.t_tests(test_id).n_perm),
0231         fprintf('q level of critical t-scores: %f., FDR method: %s\n', ...
0232             GND.t_tests(test_id).desired_alphaORq,GND.t_tests(test_id).mult_comp_method);
0233     else
0234         fprintf('Estimated alpha level of critical t-scores: %f.\n', ...
0235             GND.t_tests(test_id).estimated_alpha);
0236     end
0237 end
0238 
0239 for a=1:n_topos,
0240     subplot(ht,wdth,a);
0241     time_wind=GND.t_tests(test_id).time_wind(a,:);
0242     if strcmpi(p.Results.units,'uV') || strcmpi(p.Results.units,'dB') ,
0243         mn=mns(:,a);
0244     else
0245         mn=GND.t_tests(test_id).data_t(:,a);
0246     end
0247     
0248     if isnan(GND.t_tests(test_id).n_perm),
0249         sig_chans=find(GND.t_tests(test_id).fdr_rej(:,a));
0250     else
0251         sig_chans=find(GND.t_tests(test_id).adj_pval(:,a)<GND.t_tests(test_id).estimated_alpha);
0252     end
0253     topoplotMK(mn,GND.chanlocs(chans),'emarker2',{sig_chans,'o',[1 1 1],4}, ...
0254         'maplimits',maplimits);
0255     if freq_domain,
0256         lab_form=['%.' int2str(n_dig_past_dot) 'f-%.' int2str(n_dig_past_dot)  'f Hz'];
0257     else
0258         lab_form='%d-%d ms';
0259     end
0260     
0261     hh=title(sprintf(lab_form,time_wind(1),time_wind(2)));
0262     if ~one_scale,
0263         hcb=colorbar;
0264         hy=ylabel(hcb,units);
0265         set(hy,'rotation',0,'verticalalignment','middle','position',[4.9+n_topos*.2 0.03 1.0001]);
0266     end
0267 end
0268 
0269 if one_scale,
0270     hcb=cbar_nudge('vert',0,maplimits);
0271     hy=ylabel(hcb,units);
0272     set(hy,'rotation',0,'verticalalignment','middle');
0273 end
0274 
0275 % Plot title of entire figure
0276 if p.Results.title_on,      
0277     if isfield(GND,'bin_info'),
0278         h=textsc(['Bin ' int2str(bin) ': ' GND.bin_info(bin).bindesc],'title');
0279     else
0280         h=textsc(['Bin ' int2str(bin) ': ' GND.bindesc{bin}],'title');
0281     end
0282     set(h,'fontsize',14,'fontweight','bold');
0283     if ~verLessThan('matlab','8')
0284         set(h,'position',[.5 .975 0]);
0285     end
0286 end
0287 
0288 
0289 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0290 %% cbar_nudge (plots colorbar slightly to left of cbar.m
0291 % search for "DG change" to find the line I edited
0292 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0293 function [handle]=cbar_nudge(arg,colors,minmax, grad)
0294 
0295 if nargin < 2
0296   colors = 0;
0297 end
0298 posscale = 'off';
0299 if nargin < 1
0300   arg = 'vert';
0301   ax = [];
0302 else
0303   if isempty(arg)
0304     arg = 0;
0305   end
0306   if arg(1) == 0
0307     ax = [];
0308     arg = 'vert';
0309   elseif strcmpi(arg, 'pos')
0310     ax = [];
0311     arg = 'vert';
0312     posscale = 'on';
0313   else      
0314     if isstr(arg)
0315       ax = [];
0316     else
0317       ax = arg;
0318       arg = [];
0319     end
0320   end
0321 end
0322 
0323 if nargin>2
0324   if size(minmax,1) ~= 1 | size(minmax,2) ~= 2
0325     help cbar
0326     fprintf('cbar() : minmax arg must be [min,max]\n');
0327     return
0328   end
0329 end
0330 if nargin < 4
0331     grad = 5;
0332 end;
0333 
0334 
0335 %%%%%%%%%%%%%%%%%%%%%%%%%%%
0336 % Choose colorbar position
0337 %%%%%%%%%%%%%%%%%%%%%%%%%%%
0338 
0339 if (length(colors) == 1) & (colors == 0)
0340   t = caxis;
0341 else
0342   t = [0 1];
0343 end
0344 if ~isempty(arg)
0345   if strcmp(arg,'vert')  
0346     cax = gca;
0347     pos = get(cax,'Position');
0348     stripe = 0.04; 
0349     edge = 0.01;
0350     space = -.02;
0351 
0352     set(cax,'Position',[pos(1) pos(2) pos(3) pos(4)])
0353     rect = [pos(1)+pos(3)+space pos(2) stripe*pos(3) pos(4)];
0354     ax = axes('Position', rect);
0355   elseif strcmp(arg,'horiz')
0356     cax = gca;
0357     pos = get(cax,'Position');
0358     stripe = 0.075; 
0359     space = .1;  
0360     set(cax,'Position',...
0361         [pos(1) pos(2)+(stripe+space)*pos(4) pos(3) (1-stripe-space)*pos(4)])
0362     rect = [pos(1) pos(2) pos(3) stripe*pos(4)];
0363     ax = axes('Position', rect);
0364   end
0365 else
0366   pos = get(ax,'Position');
0367   if pos(3) > pos(4)
0368     arg = 'horiz';
0369   else
0370     arg = 'vert';
0371   end
0372 end
0373   
0374 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0375 % Draw colorbar using image()
0376 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0377 
0378 map = colormap;
0379 n = size(map,1);
0380 
0381 if length(colors) == 1
0382   if strcmp(arg,'vert')
0383       if strcmpi(posscale, 'on')
0384           image([0 1],[0 t(2)],[ceil(n/2):n-colors]');
0385       else
0386           image([0 1],t,[1:n-colors]');
0387       end;
0388       set(ax,'xticklabelmode','manual')
0389       set(ax,'xticklabel',[],'YAxisLocation','right')
0390       
0391   else
0392     image(t,[0 1],[1:n-colors]);
0393     set(ax,'yticklabelmode','manual')
0394     set(ax,'yticklabel',[],'YAxisLocation','right')
0395   end
0396   set(ax,'Ydir','normal','YAxisLocation','right')
0397 
0398 else % length > 1
0399 
0400   if max(colors) > n
0401     error('Color vector excedes size of colormap')
0402   end
0403   if strcmp(arg,'vert')
0404     image([0 1],t,[colors]');
0405     set(ax,'xticklabelmode','manual')
0406     set(ax,'xticklabel',[])
0407   else
0408     image([0 1],t,[colors]);
0409     set(ax,'yticklabelmode','manual')
0410     set(ax,'yticklabel',[],'YAxisLocation','right')
0411   end  
0412   set(ax,'Ydir','normal','YAxisLocation','right')
0413 end
0414 
0415 %%%%%%%%%%%%%%%%%%%%%%%%%
0416 % Adjust cbar ticklabels
0417 %%%%%%%%%%%%%%%%%%%%%%%%%
0418 
0419 if nargin > 2 
0420   if strcmp(arg,'vert')
0421       Cax = get(ax,'Ylim');
0422   else
0423       Cax = get(ax,'Xlim');
0424   end;
0425   CBTicks = [Cax(1):(Cax(2)-Cax(1))/(grad-1):Cax(2)]; % caxis tick positions
0426   CBLabels = [minmax(1):(minmax(2)-minmax(1))/(grad-1):minmax(2)]; % tick labels
0427   
0428   dec = floor(log10(max(abs(minmax)))); % decade of largest abs value
0429   CBLabels = ([minmax]* [ linspace(1,0, grad);linspace(0, 1, grad)]);
0430   if dec<1
0431     CBLabels = round(CBLabels*10^(1-dec))*10^(dec-1);
0432   elseif dec == 1
0433     CBLabels = round(CBLabels*10^(2-dec))*10^(dec-2);
0434   else
0435     CBLabels = round(CBLabels);
0436   end
0437 
0438   if strcmp(arg,'vert')
0439     set(ax,'Ytick',CBTicks);
0440     set(ax,'Yticklabel',CBLabels);
0441   else
0442     set(ax,'Xtick',CBTicks);
0443     set(ax,'Xticklabel',CBLabels);
0444   end
0445 end
0446 handle = ax;
0447 
0448 %%%%%%%%%%%%%%%%%%
0449 % Adjust cbar tag
0450 %%%%%%%%%%%%%%%%%%
0451 
0452 set(ax,'tag','cbar')