sig_raster

 sig_raster() - Plots 2D (electrode x time/frequency) raster diagram of 
                variables that show significant effects according to t-tests. 
                Significant positive effects are represented 
                with white squares (red if color is requested).  
                Significant negative effects are represented with
                black squares (blue if color is requested).
                
 Usage:
  >> [img, h_ax]=sig_raster(GND_GRP_specGND_or_fname,test_id,varargin);

 Required Inputs:
  GND_GRP_specGND_or_fname  - A GND/GRP/specGND structure variable or the filename of 
                      such a variable that has been saved to disk.   
                      To create a GND variable from Kutaslab ERP files (e.g.,
                      *.nrm files) use avgs2GND.m.  To do the same from 
                      EEGLAB *.set files use sets2GND.m.  To create a
                      a GRP structure use GNDs2GRP.m. See Mass Univariate ERP  
                      Toolbox documentation for detailed information about the 
                      format of GND and GRP variables. 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:
  x_ticks          - [vector] The times/frequencies you want marked with ticks
                     on the x-axis.  Note, because the EEG has been sampled 
                     at discrete time points, not all times between the 
                     minimum and maximum analyzed times/frequencies can be
                     used.  When a tick is requested that is not available, 
                     the closest possible value is used. If not
                     specified, x_ticks are automatically chosen based
                     on the time/frequency range covered by the diagram. 
                     This option has no effect if the t-test was
                     executed based on mean voltages/power within specified 
                     time windows/frequency bands.
  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}
  use_color        - ['n','rb','rgb'] If 'n', raster will be black, white,
                     and grey.  If 'rb', raster will be red (significantly
                     positive), blue (significantly negative), and white
                     (nonsignificant). If 'rgb', raster will use a 
                     temperature scale to indicate graded values of 
                     significance. {default: 'n'} 
  units            - ['t' or 'uV'] If 'use_color' is set to 'rgb', 't'
                     means that raster values will be in units of t-scores. 
                     'uV' means that raster values will be in microvolts 
                     (minus the mean of the null hypothesis). {default: 't'}
  plot_vert_lines  - [integer] If non-zero, vertical lines separating all
                     time points included in the t-test will be 
                     shown. If 0, the only vertical lines drawn will be 
                     those that separate times/frequencies included in a test 
                     from those not included. This option has no effect if 
                     the t-test was executed based on mean voltage/power 
                     within specified time windows. {default: 1}
  lr_sym           - [integer] If non-zero, left hemisphere electrodes are
                     plot from most anterior to most posterior and right 
                     hemisphere electrodes are plot from most posterior to 
                     most anterior.  This may make it easier to visualize 
                     the degree of effect left-right symmetry.  If 0, all
                     electrodes are plot from most anterior to most
                     posterior.
  scale_limits     - [cmin cmax] The minimum and maximum value of the
                     colorscale when 'use_color' option is set to 'rgb'. 
                     If 'use_color' is NOT set to 'rgb', this option has 
                     no effect. If 'scale_limits' is not specified the 
                     limits are +/- the maximum absolute value being shown.
                     This option is useful for making multiple sig_raster
                     plots on the same scale.
  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)

 Outputs:
  img              - Two-dimensional matrix (channel x time/frequency) 
                     representing the visualized raster.  This could be
                     useful for making your own personalized raster plot
                     (e.g., >> imagesc(img);).  Note that the matrix
                     includes rows for any blank rows separating groups of
                     channels (e.g., left electrodes from midline
                     electrodes).
  h_ax             - h_ax(1) is the handle of the raster axis. h_ax(2) is
                     the handle of the colorscale (if used).
         

 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.

 -Some image viewers (e.g., Apple's Preview) may have trouble viewing 
 postscript raster images since they downsample the image to make the file 
 smaller, which blurs the lines.  Ghostview  and Adobe Illustrator should
 be able to view the files just fine.  If viewing the raster as a
 postscript file is producing problems, you can save the figure in jpg 
 format.

 -Assignment of electrode to left, midline, or right grouping of
 electrodes is based on the GND.chanlocs(x).theta coordinate.  Anterior to
 posterior organization of electrodes is based on GND.chanlocs(x).radius.

 Author:
 David Groppe
 Kutaslab, 3/2010

0001 % sig_raster() - Plots 2D (electrode x time/frequency) raster diagram of
0002 %                variables that show significant effects according to t-tests.
0003 %                Significant positive effects are represented
0004 %                with white squares (red if color is requested).
0005 %                Significant negative effects are represented with
0006 %                black squares (blue if color is requested).
0007 %
0008 % Usage:
0009 %  >> [img, h_ax]=sig_raster(GND_GRP_specGND_or_fname,test_id,varargin);
0010 %
0011 % Required Inputs:
0012 %  GND_GRP_specGND_or_fname  - A GND/GRP/specGND structure variable or the filename of
0013 %                      such a variable that has been saved to disk.
0014 %                      To create a GND variable from Kutaslab ERP files (e.g.,
0015 %                      *.nrm files) use avgs2GND.m.  To do the same from
0016 %                      EEGLAB *.set files use sets2GND.m.  To create a
0017 %                      a GRP structure use GNDs2GRP.m. See Mass Univariate ERP
0018 %                      Toolbox documentation for detailed information about the
0019 %                      format of GND and GRP variables. If you specifiy a filename
0020 %                      be sure to include the file's path, unless the file is
0021 %                      in the current working directory.
0022 %  test_id           - [integer] The index # of the t-test results
0023 %                      stored in the GND/GRP/specGND variable that you wish to visualize.
0024 %                      To see what test results are available, look at
0025 %                      the "t_tests" field of your variable (e.g., GND.t_tests)
0026 %                      or use the functions headinfo.m or headinfo_spec.m.
0027 %
0028 % Optional Inputs:
0029 %  x_ticks          - [vector] The times/frequencies you want marked with ticks
0030 %                     on the x-axis.  Note, because the EEG has been sampled
0031 %                     at discrete time points, not all times between the
0032 %                     minimum and maximum analyzed times/frequencies can be
0033 %                     used.  When a tick is requested that is not available,
0034 %                     the closest possible value is used. If not
0035 %                     specified, x_ticks are automatically chosen based
0036 %                     on the time/frequency range covered by the diagram.
0037 %                     This option has no effect if the t-test was
0038 %                     executed based on mean voltages/power within specified
0039 %                     time windows/frequency bands.
0040 %  fig_id           - [integer] The index # of the MATLAB figure in which
0041 %                     the diagram will be produced.  Useful for overwriting
0042 %                     old figures. {default: lowest unused index}
0043 %  use_color        - ['n','rb','rgb'] If 'n', raster will be black, white,
0044 %                     and grey.  If 'rb', raster will be red (significantly
0045 %                     positive), blue (significantly negative), and white
0046 %                     (nonsignificant). If 'rgb', raster will use a
0047 %                     temperature scale to indicate graded values of
0048 %                     significance. {default: 'n'}
0049 %  units            - ['t' or 'uV'] If 'use_color' is set to 'rgb', 't'
0050 %                     means that raster values will be in units of t-scores.
0051 %                     'uV' means that raster values will be in microvolts
0052 %                     (minus the mean of the null hypothesis). {default: 't'}
0053 %  plot_vert_lines  - [integer] If non-zero, vertical lines separating all
0054 %                     time points included in the t-test will be
0055 %                     shown. If 0, the only vertical lines drawn will be
0056 %                     those that separate times/frequencies included in a test
0057 %                     from those not included. This option has no effect if
0058 %                     the t-test was executed based on mean voltage/power
0059 %                     within specified time windows. {default: 1}
0060 %  lr_sym           - [integer] If non-zero, left hemisphere electrodes are
0061 %                     plot from most anterior to most posterior and right
0062 %                     hemisphere electrodes are plot from most posterior to
0063 %                     most anterior.  This may make it easier to visualize
0064 %                     the degree of effect left-right symmetry.  If 0, all
0065 %                     electrodes are plot from most anterior to most
0066 %                     posterior.
0067 %  scale_limits     - [cmin cmax] The minimum and maximum value of the
0068 %                     colorscale when 'use_color' option is set to 'rgb'.
0069 %                     If 'use_color' is NOT set to 'rgb', this option has
0070 %                     no effect. If 'scale_limits' is not specified the
0071 %                     limits are +/- the maximum absolute value being shown.
0072 %                     This option is useful for making multiple sig_raster
0073 %                     plots on the same scale.
0074 %  verblevel        - An integer specifiying the amount of information you want
0075 %                     this function to provide about what it is doing during runtime.
0076 %                     Options are:
0077 %                      0 - quiet, only show errors, warnings, and EEGLAB reports
0078 %                      1 - stuff anyone should probably know
0079 %                      2 - stuff you should know the first time you start working
0080 %                          with a data set {default value}
0081 %                      3 - stuff that might help you debug (show all
0082 %                          reports)
0083 %
0084 % Outputs:
0085 %  img              - Two-dimensional matrix (channel x time/frequency)
0086 %                     representing the visualized raster.  This could be
0087 %                     useful for making your own personalized raster plot
0088 %                     (e.g., >> imagesc(img);).  Note that the matrix
0089 %                     includes rows for any blank rows separating groups of
0090 %                     channels (e.g., left electrodes from midline
0091 %                     electrodes).
0092 %  h_ax             - h_ax(1) is the handle of the raster axis. h_ax(2) is
0093 %                     the handle of the colorscale (if used).
0094 %
0095 %
0096 % Global Variable:
0097 %   VERBLEVEL - Mass Univariate ERP Toolbox level of verbosity (i.e., tells functions
0098 %               how much to report about what they're doing during
0099 %               runtime) set by the optional function argument 'verblevel'
0100 %
0101 % Notes:
0102 % -The printed/exported figure will have the same dimensions as the figure
0103 % on the display.  Thus you can undo figure clutter by re-sizing it.
0104 %
0105 % -Some image viewers (e.g., Apple's Preview) may have trouble viewing
0106 % postscript raster images since they downsample the image to make the file
0107 % smaller, which blurs the lines.  Ghostview  and Adobe Illustrator should
0108 % be able to view the files just fine.  If viewing the raster as a
0109 % postscript file is producing problems, you can save the figure in jpg
0110 % format.
0111 %
0112 % -Assignment of electrode to left, midline, or right grouping of
0113 % electrodes is based on the GND.chanlocs(x).theta coordinate.  Anterior to
0114 % posterior organization of electrodes is based on GND.chanlocs(x).radius.
0115 %
0116 % Author:
0117 % David Groppe
0118 % Kutaslab, 3/2010
0119 
0120 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
0121 % 3/27/2010-Comments revised to accommodate GRP variable
0122 %
0123 % 4/8/2010-chanloc theta values standardized to a range of +/- 180 (A2 had
0124 % a theta of 240 for some reason)
0125 %
0126 % 10/5/2010-channels lacking a theta coordinate are assigned a theta
0127 % coordinate of 0 and radius of 2
0128 %
0129 % 12/14/2010-Compatible with specGND variables now too
0130 %
0131 % 3/29/2011-Classifies any electrode with a radius of 0 as a midline
0132 % electrode
0133 %
0134 % 3/5/2012-'rgb' color option, 'units' option, and 'scale_limits' added.
0135 % h_ax output added.
0136 
0137 %%%%%%%%%%%%%%%% FUTURE WORK %%%%%%%%%%%%%%%%%
0138 % -When you click on figure, box with time and electrode appear behind edges
0139 % of figure.  Fix this somehow?
0140 % -Make it possible to further exclude channels from plot?
0141 % -Make it possible to deal with overlapping time windows when mean time
0142 % windows are used?
0143 
0144 
0145 function [img, h_ax]=sig_raster(GND_GRP_specGND_or_fname,test_id,varargin)
0146 
0147 p=inputParser;
0148 p.addRequired('GND_GRP_specGND_or_fname',@(x) isstruct(x) || ischar(x)); 
0149 p.addRequired('test_id',@(x) isnumeric(x) && (length(x)==1));
0150 p.addParamValue('x_ticks',[],@isnumeric);
0151 p.addParamValue('fig_id',[],@(x) isnumeric(x) && (length(x)==1));
0152 p.addParamValue('use_color','n',@(x) strcmpi(x,'n') || strcmpi(x,'rb') || strcmpi(x,'rgb'));
0153 p.addParamValue('units','t',@(x) strcmpi(x,'t') || strcmpi(x,'uV'));
0154 p.addParamValue('plot_vert_lines',1,@(x) isnumeric(x) || ischar(x));
0155 p.addParamValue('lr_sym',0,@(x) isnumeric(x) || ischar(x));
0156 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1));
0157 p.addParamValue('scale_limits',[],@(x) isnumeric(x) && (length(x)==2));
0158 p.addParamValue('unused_color',[1 1 1]*.7,@(x) isnumeric(x) && (length(x)==3)); % this option isn't all that useful actually
0159 
0160 p.parse(GND_GRP_specGND_or_fname,test_id,varargin{:});
0161 
0162 % Manage VERBLEVEL
0163 if isempty(p.Results.verblevel),
0164     VERBLEVEL=2; %not global, just local
0165 else
0166     VERBLEVEL=p.Results.verblevel;
0167 end
0168 
0169 %use_color=str2bool(p.Results.use_color);
0170 use_color=p.Results.use_color;
0171 if ~verLessThan('matlab','8.4') && strcmpi(use_color,'rgb')
0172     warning('rbg option for use_color is not yet supported for your version of Matlab. Reverting to grey scale. Try an older version of Matlab.');
0173     use_color='n';
0174 end
0175 plot_vert_lines=str2bool(p.Results.plot_vert_lines);
0176 lr_sym=str2bool(p.Results.lr_sym);
0177 
0178 %Load GND, GRP, or specGND struct
0179 if ischar(GND_GRP_specGND_or_fname),
0180     fprintf('Loading GND, GRP, or specGND struct variable from file %s.\n',GND_GRP_specGND_or_fname);
0181     load(GND_GRP_specGND_or_fname,'-MAT');
0182     if ~exist('GND','var') && ~exist('GRP','var') && ~exist('specGND','var')
0183         error('File %s does not contain a GND, GRP, or specGND variable.',GND_GRP_specGND_or_fname);
0184     end
0185     if exist('GRP','var'),
0186         GND=GRP; %for simplicity GRP variables are re-named GND
0187         clear GRP;
0188     elseif exist('specGND','var'),
0189         GND=specGND; %for simplicity specGND variables are re-named GND
0190         clear specGND;
0191     end
0192     VerbReport(sprintf('Experiment: %s',GND.exp_desc),2,VERBLEVEL);
0193 else
0194     GND=GND_GRP_specGND_or_fname; %for simplicity GRP and specGND variables are re-named GND
0195     clear GND_GRP_specGND_or_fname;
0196 end
0197 
0198 n_test=length(GND.t_tests);
0199 if test_id>n_test,
0200    error('There are only %d t-tests stored with these data, but you requested Test #%d',n_test,test_id); 
0201 end
0202 
0203 %FDR or permutation test correction for multiple comparisons
0204 if isnan(GND.t_tests(test_id).n_perm)
0205     fdr_crct=1;
0206     if VERBLEVEL,
0207        fprintf('Correcting for multiple comparisons via FDR procedure: %s\n', ...
0208            GND.t_tests(test_id).mult_comp_method); 
0209     end
0210 else
0211     fdr_crct=0;
0212     if VERBLEVEL,
0213        fprintf('Correcting for multiple comparisons via permutation test.\n');  
0214     end
0215 end
0216 
0217 if isfield(GND.t_tests(1),'freq_band')
0218     %create temporary fields to make frequency domain plotting
0219     % easily compatible with time domain plot
0220     GND.t_tests(test_id).time_wind=GND.t_tests(test_id).freq_band;
0221     GND.t_tests(test_id).mean_wind=GND.t_tests(test_id).mean_band;
0222     GND.t_tests(test_id).used_tpt_ids=GND.t_tests(test_id).used_freq_ids;
0223     GND.grands_t=GND.grands_pow_dB_t;
0224     freq_step=GND.freqs(2)-GND.freqs(1);
0225     freq_ord=orderofmag(freq_step);
0226     ord_pow=log10(freq_ord); %useful for formatting x-tick labels
0227     if ord_pow>=0,
0228         n_dig_past_dot=0;
0229     else
0230         n_dig_past_dot=-ord_pow;
0231     end
0232     GND.time_pts=rnd_orderofmag(GND.freqs);
0233     freq_domain=1;
0234 else
0235     freq_domain=0;
0236 end
0237 
0238 %Clean up code by copying some GND.t_tests variables to their own
0239 %variable
0240 use_chans=GND.t_tests(test_id).used_chan_ids;
0241 n_wind=size(GND.t_tests(test_id).time_wind,1);
0242 use_bin=GND.t_tests(test_id).bin;
0243 if strcmpi(GND.t_tests(test_id).mean_wind,'yes'),
0244     use_tpts=zeros(1,length(GND.time_pts)); %preallocate mem
0245     pval=repmat(use_tpts,length(GND.t_tests(test_id).used_chan_ids),1);
0246     grands_t=pval;
0247     for w=1:n_wind,
0248         ids=GND.t_tests(test_id).used_tpt_ids{w};
0249         use_tpts(ids)= ...
0250             use_tpts(ids)+1;
0251         if fdr_crct,
0252             pval(:,ids)=repmat(1-GND.t_tests(test_id).fdr_rej(:,w),1,length(ids));
0253         else
0254             pval(:,ids)=repmat(GND.t_tests(test_id).adj_pval(:,w),1,length(ids));
0255         end
0256         if strcmpi(p.Results.units,'t')
0257             grands_t(:,ids)=repmat(GND.t_tests(test_id).data_t(:,w),1,length(ids));
0258         else
0259             temp_mn=mean(GND.grands(GND.t_tests(test_id).used_chan_ids, ...
0260                 GND.t_tests(3).used_tpt_ids{w},use_bin),2);
0261             temp_mn=temp_mn-GND.t_tests(test_id).null_mean; % Deals with possible non-zero null hypothesis mean
0262             grands_t(:,ids)=repmat(temp_mn,1,length(ids));
0263         end
0264     end
0265     if max(use_tpts)>1,
0266         if freq_domain,
0267             error('Tested frequency bands overlap. sig_raster.m can''t handle this currently.');
0268         else
0269             error('Tested time windows overlap. sig_raster.m can''t handle this currently.');
0270         end
0271     end
0272     ids=find(use_tpts);
0273     use_tpts=ids;
0274     pval=pval(:,ids);
0275     grands_t=grands_t(:,ids);
0276 else
0277     use_tpts=GND.t_tests(test_id).used_tpt_ids;
0278     if fdr_crct,
0279         pval=1-GND.t_tests(test_id).fdr_rej;
0280     else
0281         pval=GND.t_tests(test_id).adj_pval;
0282     end
0283     if GND.t_tests(test_id).null_mean,
0284         %recompute t-scores because mean of null hypothesis is not zero
0285         if strcmpi(p.Results.units,'t')
0286             grands_t=squeeze( (GND.grands(GND.t_tests(test_id).used_chan_ids, ...
0287                 use_tpts,use_bin)-GND.t_tests(test_id).null_mean)./ ...
0288                 GND.grands_stder(GND.t_tests(test_id).used_chan_ids, ...
0289                 use_tpts,use_bin));
0290         else
0291             grands_t=squeeze( (GND.grands(GND.t_tests(test_id).used_chan_ids, ...
0292                 use_tpts,use_bin)-GND.t_tests(test_id).null_mean)); %Note, grands_t is storing voltage data despite its name
0293         end
0294     else
0295         if strcmpi(p.Results.units,'t')
0296             grands_t=squeeze(GND.grands_t(GND.t_tests(test_id).used_chan_ids, ...
0297                 use_tpts,use_bin));
0298         else
0299             grands_t=squeeze(GND.grands(GND.t_tests(test_id).used_chan_ids, ...
0300                 use_tpts,use_bin)); %Note, grands_t is storing voltage data despite its name
0301         end
0302     end
0303 end
0304 
0305 % Figure out electrode ordering
0306 n_use_chans=length(use_chans);
0307 
0308 % Standardize theta to -180<theta<=180 degrees
0309 for c=1:n_use_chans,
0310     theta=mod(GND.chanlocs(use_chans(c)).theta,360);
0311     if isempty(theta),
0312         watchit(sprintf('Channel #%d (Label: %s) does not have a theta coordinate. To produce raster, I''m assigning it temporary coordinate of: theta=0.',c,GND.chanlocs(use_chans(c)).labels));
0313         theta=0;
0314     end
0315     if theta>180,
0316         theta=theta-360;
0317     elseif theta<=-180,
0318         theta=360+theta;
0319     end
0320     GND.chanlocs(use_chans(c)).theta=theta;
0321     
0322     if isempty(GND.chanlocs(use_chans(c)).radius),
0323         GND.chanlocs(use_chans(c)).radius=2;
0324         watchit(sprintf('Channel #%d (Label: %s) does not have a radius coordinate. To produce raster, I''m assigning it temporary coordinate of: radius=2.',c,GND.chanlocs(use_chans(c)).labels));
0325     end
0326 end
0327 
0328 % Get midline electrodes from front to back
0329 midline=[];
0330 midline_dist=[]; %distance from MiCe
0331 for c=1:n_use_chans,
0332     theta=GND.chanlocs(use_chans(c)).theta;
0333     radius=GND.chanlocs(use_chans(c)).radius;
0334 
0335     if theta==0 || radius==0,
0336         midline=[midline use_chans(c)];
0337         midline_dist=[midline_dist radius];
0338     elseif theta==180,
0339         midline=[midline use_chans(c)];
0340         midline_dist=[midline_dist -radius];
0341     end
0342 end
0343 [midline_dist id]=sort(midline_dist,2,'descend');
0344 midline=midline(id);
0345 n_midline=length(midline);
0346 
0347 % Get left electrodes from front to back
0348 left=[];
0349 left_dist=[];
0350 for c=1:n_use_chans,
0351     %standardize theta to a range between +/-180
0352     theta=GND.chanlocs(use_chans(c)).theta;
0353     
0354     if (theta<0) && (theta>-180)
0355         if GND.chanlocs(use_chans(c)).radius,
0356             %make sure electrode isn't midline or right
0357             left=[left use_chans(c)];
0358             ang=90-abs(theta);
0359             dist=GND.chanlocs(use_chans(c)).radius*sin(ang*2*pi/360);
0360             left_dist=[left_dist dist];
0361         end
0362     end
0363 end
0364 [left_dist id]=sort(left_dist,2,'descend');
0365 left=left(id);
0366 n_left=length(left);
0367 
0368 % Get right electrodes from front to back
0369 right=[];
0370 right_dist=[];
0371 for c=1:n_use_chans,
0372     %standardize theta to a range between +/-180
0373     theta=GND.chanlocs(use_chans(c)).theta;
0374     
0375     if (theta>0) && (theta<180)
0376         if GND.chanlocs(use_chans(c)).radius,
0377             %make sure electrode isn't midline or right
0378             right=[right use_chans(c)];
0379             ang=90-abs(theta);
0380             dist=GND.chanlocs(use_chans(c)).radius*sin(ang*2*pi/360);
0381             right_dist=[right_dist dist];
0382         end
0383     end
0384 end
0385 % Make right electrodes back to front if requested
0386 if lr_sym,
0387    right_dist=-right_dist; 
0388 end
0389 [right_dist id]=sort(right_dist,2,'descend');
0390 right=right(id);
0391 n_right=length(right);
0392 
0393 % Compute time/frequency tick mark locations and labels
0394 show_tpts=use_tpts(1):use_tpts(end); %show_tpts differs from use_tpts if multiple non-contiguous time windows were used
0395 n_show_tpts=length(show_tpts);
0396 if strcmpi(GND.t_tests(test_id).mean_wind,'yes'),
0397     % Test was based on mean voltage in time window(s)/frequency band(s)
0398     xtick=NaN;
0399     n_xtick=0;
0400     for a=1:n_wind,
0401         n_xtick=n_xtick+1;
0402         mn_tpt=mean(GND.t_tests(test_id).used_tpt_ids{a});
0403         xtick(n_xtick)=mn_tpt-show_tpts(1)+1;
0404         if freq_domain,
0405             lab_form=['%.' int2str(n_dig_past_dot) 'f-%.' int2str(n_dig_past_dot)  'f'];
0406             xtick_lab{n_xtick}=sprintf(lab_form,GND.t_tests(test_id).time_wind(a,1), ...
0407                 GND.t_tests(test_id).time_wind(a,2));
0408         else
0409             xtick_lab{n_xtick}=sprintf('%d-%d',GND.t_tests(test_id).time_wind(a,1), ...
0410                 GND.t_tests(test_id).time_wind(a,2));
0411         end
0412         if a~=n_wind,
0413             skipped_wind=[(GND.t_tests(test_id).used_tpt_ids{a}(end)+1): ...
0414                 (GND.t_tests(test_id).used_tpt_ids{a+1}(1)-1)];
0415             if ~isempty(skipped_wind),
0416                 n_xtick=n_xtick+1;
0417                 mn_tpt=mean(skipped_wind);
0418                 xtick(n_xtick)=mn_tpt-show_tpts(1)+1;
0419                 if freq_domain,
0420                     if length(skipped_wind)>1,
0421                         lab_form=['%.' int2str(n_dig_past_dot) 'f-%.' int2str(n_dig_past_dot)  'f'];
0422                         xtick_lab{n_xtick}=sprintf(lab_form,GND.time_pts(skipped_wind(1)), ...
0423                             GND.time_pts(skipped_wind(end)));
0424                     else
0425                         lab_form=['%.' int2str(n_dig_past_dot) 'f'];
0426                         xtick_lab{n_xtick}=sprintf(lab_form,GND.time_pts(skipped_wind(1)));
0427                     end
0428                 else
0429                     if length(skipped_wind)>1,
0430                         xtick_lab{n_xtick}=sprintf('%d-%d',GND.time_pts(skipped_wind(1)), ...
0431                             GND.time_pts(skipped_wind(end)));
0432                     else
0433                         xtick_lab{n_xtick}=sprintf('%d',GND.time_pts(skipped_wind(1)), ...
0434                             GND.time_pts(skipped_wind(end)));
0435                     end
0436                 end
0437             end
0438         end
0439     end
0440 else
0441     % Test was based on each time point in time window(s)/frequency band(s)
0442     if isempty(p.Results.x_ticks),
0443         tm_range=GND.time_pts(use_tpts(end))-GND.time_pts(use_tpts(1));
0444         omag=orderofmag(tm_range);
0445         tm_step=round(omag*GND.srate/1000);
0446         xtick=1:tm_step:n_show_tpts;
0447         xtick_lab=cell(1,length(xtick));
0448         for a=1:length(xtick),
0449             xtick_lab{a}=num2str(GND.time_pts(show_tpts(xtick(a))));
0450         end
0451     else
0452         tk_ct=0;
0453         xtick=zeros(1,length(p.Results.x_ticks));
0454         for t=p.Results.x_ticks,
0455             tk_ct=tk_ct+1;
0456             xtick(tk_ct)=find_tpt(GND.time_pts(show_tpts),t);
0457         end
0458         xtick=unique(xtick); %get rid of any redundant ticks
0459         xtick_lab=cell(1,length(xtick));
0460         for a=1:length(xtick),
0461             xtick_lab{a}=num2str(GND.time_pts(show_tpts(xtick(a))));
0462         end
0463     end
0464 end
0465 
0466 if VERBLEVEL,
0467     if fdr_crct,
0468         fprintf('q level of critical t-scores: %f.\n',GND.t_tests(test_id).desired_alphaORq);
0469     else
0470         fprintf('Estimated alpha level of critical t-scores: %f.\n',GND.t_tests(test_id).estimated_alpha);
0471     end
0472 end
0473 
0474 %Initialize variables
0475 n_blank=(n_right>0)+(n_left>0)+(n_midline>0)-1;
0476 if ~strcmpi(use_color,'n'),
0477     img=-.25*ones(n_use_chans+n_blank,n_show_tpts);%makes skipped lines grey (non-sig rectangles are white)
0478 else
0479     img=zeros(n_use_chans+n_blank,n_show_tpts); %makes skipped lines grey
0480 end
0481 if strcmpi(use_color,'rgb')
0482     mask=img;
0483 end
0484 chan_lab=cell(1,n_use_chans);
0485 dat.chan_lab=cell(1,n_use_chans+n_blank); %channel labels including empty cells for skipped lines
0486 %dat_chan_lab is needed to make sig_raster plot interactive
0487 
0488 %Get IDs of time points/frequencies used in perm test (some shown time points/frequencies may not
0489 %have been included
0490 used_tpt_ids=zeros(1,n_show_tpts);
0491 for a=1:n_show_tpts,
0492    if ismember(show_tpts(a),use_tpts),
0493       used_tpt_ids(a)=1; 
0494    end
0495 end
0496 used_tpt_ids=find(used_tpt_ids);
0497 
0498 %Construct raster for left electrodes
0499 ct=0;
0500 skipped=[];
0501 for a=left,
0502     ct=ct+1;
0503     elec_id=find(use_chans==a);
0504     if fdr_crct,
0505         if strcmpi(use_color,'rgb')
0506             mask(ct,used_tpt_ids)=(pval(elec_id,:)<=GND.t_tests(test_id).desired_alphaORq);
0507             img(ct,used_tpt_ids)=grands_t(elec_id,:);
0508         else
0509             img(ct,used_tpt_ids)=(pval(elec_id,:)<=GND.t_tests(test_id).desired_alphaORq).*sign(grands_t(elec_id,:));
0510         end
0511     else
0512         if strcmpi(use_color,'rgb')
0513             mask(ct,used_tpt_ids)=(pval(elec_id,:)<GND.t_tests(test_id).estimated_alpha);
0514             img(ct,used_tpt_ids)=grands_t(elec_id,:);
0515         else
0516             img(ct,used_tpt_ids)=(pval(elec_id,:)<GND.t_tests(test_id).estimated_alpha).*sign(grands_t(elec_id,:));
0517         end
0518     end
0519     chan_lab{ct}=GND.chanlocs(use_chans(elec_id)).labels;
0520     dat.chan_lab{ct}=GND.chanlocs(use_chans(elec_id)).labels;
0521 end
0522 
0523 if ~isempty(midline),
0524     if ct>0,
0525         %skip line to separate left and midline electrodes
0526         ct=ct+1;
0527         skipped=ct;
0528     end
0529     for a=midline,
0530         ct=ct+1;
0531         elec_id=find(use_chans==a);
0532         if fdr_crct,
0533             if strcmpi(use_color,'rgb')
0534                 mask(ct,used_tpt_ids)=(pval(elec_id,:)<=GND.t_tests(test_id).desired_alphaORq);
0535                 img(ct,used_tpt_ids)=grands_t(elec_id,:);
0536             else
0537                 img(ct,used_tpt_ids)=(pval(elec_id,:)<=GND.t_tests(test_id).desired_alphaORq).*sign(grands_t(elec_id,:));
0538             end
0539         else
0540             if strcmpi(use_color,'rgb')
0541                 mask(ct,used_tpt_ids)=(pval(elec_id,:)<GND.t_tests(test_id).estimated_alpha);
0542                 img(ct,used_tpt_ids)=grands_t(elec_id,:);
0543             else
0544                 img(ct,used_tpt_ids)=(pval(elec_id,:)<GND.t_tests(test_id).estimated_alpha).*sign(grands_t(elec_id,:));
0545             end
0546         end
0547         chan_lab{ct-length(skipped)}=GND.chanlocs(use_chans(elec_id)).labels;
0548         dat.chan_lab{ct}=GND.chanlocs(use_chans(elec_id)).labels;
0549     end
0550 end
0551 
0552 if ~isempty(right),
0553     if ct>0,
0554         %skip line to separate right from left or midline electrodes
0555         ct=ct+1;
0556         skipped=[skipped ct];
0557     end
0558     for a=right,
0559         ct=ct+1;
0560         elec_id=find(use_chans==a);
0561         if fdr_crct,
0562             if strcmpi(use_color,'rgb')
0563                 mask(ct,used_tpt_ids)=(pval(elec_id,:)<=GND.t_tests(test_id).desired_alphaORq);
0564                 img(ct,used_tpt_ids)=grands_t(elec_id,:);
0565             else
0566                 img(ct,used_tpt_ids)=(pval(elec_id,:)<=GND.t_tests(test_id).desired_alphaORq).*sign(grands_t(elec_id,:));
0567             end
0568         else
0569             if strcmpi(use_color,'rgb')
0570                 mask(ct,used_tpt_ids)=(pval(elec_id,:)<GND.t_tests(test_id).estimated_alpha);
0571                 img(ct,used_tpt_ids)=grands_t(elec_id,:);
0572             else
0573                 img(ct,used_tpt_ids)=(pval(elec_id,:)<GND.t_tests(test_id).estimated_alpha).*sign(grands_t(elec_id,:));
0574             end
0575         end
0576         chan_lab{ct-length(skipped)}=GND.chanlocs(use_chans(elec_id)).labels;
0577         dat.chan_lab{ct}=GND.chanlocs(use_chans(elec_id)).labels;
0578     end
0579 end
0580 
0581 
0582 if isempty(p.Results.fig_id)
0583     fig_h=figure;
0584 else
0585     fig_h=figure(p.Results.fig_id); clf;
0586 end
0587 if isfield(GND,'bin_info')
0588     set(fig_h,'name',['Bin ' int2str(use_bin) ' [' GND.bin_info(use_bin).bindesc ']'],'paperpositionmode','auto');
0589 else
0590     set(fig_h,'name',['Bin ' int2str(use_bin) ' [' GND.bindesc{use_bin} ']'],'paperpositionmode','auto');
0591 end
0592 %setting paperpositionmode to 'auto' means that if the figure is manually
0593 %resized, the printed version of the figure will reflect the whatever the
0594 %shown size was (at the time of printing)
0595 
0596 if strcmpi(use_color,'rb'),
0597     %skipped rows/columns=grey, nonsig=white, +sig=red, -sig=blue
0598     h_img=imagesc(img,[-1 1]);colormap([0 0 1; .5 .5 .5; 1 1 1; 1 0 0]);
0599 elseif strcmpi(use_color,'rgb')
0600     cmap=colormap('jet');
0601     %cmap(32,:)=[1 1 1]*.7; %set 0 vals to grey
0602     zero_color=cmap(33,:);
0603     %%cmap(33,:)=[1 1 1]*.7; %set 0 vals to grey ORIG WORKS
0604     cmap(33,:)=p.Results.unused_color; %set 0 vals to grey
0605     colormap(cmap);
0606     abs_mx=max(max(abs(img)));
0607     if isempty(p.Results.scale_limits)
0608         h_img=imagesc(img.*mask,[-1 1]*abs_mx);
0609     else
0610         h_img=imagesc(img.*mask,p.Results.scale_limits);
0611     end
0612 else
0613     h_img=imagesc(img,[-1 1]);colormap([0 0 0; .6 .6 .6; 1 1 1]);
0614 end
0615 h_ax(1)=gca;
0616 if freq_domain,
0617     bdfcn=['Cp = get(gca,''CurrentPoint''); ' ...
0618         'Xp=Cp(2,1);', ...
0619         'Yp=Cp(2,2);', ...
0620         'dat=get(gcbo,''userdata'');', ...
0621         'ht=text(Xp,Yp,[num2str(dat.times(round(Xp))) '' Hz, '' dat.chan_lab{round(Yp)}]);' ...
0622         'set(ht,''backgroundcolor'',''w'',''horizontalalignment'',''center'',''verticalalignment'',''middle'',''buttondownfcn'',''delete(gcbo);'');'];
0623 else
0624     bdfcn=['Cp = get(gca,''CurrentPoint''); ' ...
0625         'Xp=Cp(2,1);', ...
0626         'Yp=Cp(2,2);', ...
0627         'dat=get(gcbo,''userdata'');', ...
0628         'ht=text(Xp,Yp,[int2str(dat.times(round(Xp))) '' ms, '' dat.chan_lab{round(Yp)}]);' ...
0629         'set(ht,''backgroundcolor'',''w'',''horizontalalignment'',''center'',''verticalalignment'',''middle'',''buttondownfcn'',''delete(gcbo);'');'];
0630 end
0631 dat.times=GND.time_pts(show_tpts); %dat.chan_lab has been set above
0632 set(h_img,'buttondownfcn',bdfcn,'userdata',dat);
0633 set(gca,'xtick',xtick,'xticklabel',xtick_lab,'tickdir','out');
0634 
0635 set(gca,'ytick',setdiff(1:(n_use_chans+n_blank),skipped),'yticklabel',chan_lab,'box','off');
0636 
0637 if freq_domain,
0638     h=xlabel('Hz');
0639 else
0640     h=xlabel('Time (ms)');
0641 end
0642 set(h,'fontsize',12,'fontweight','bold');
0643 
0644 h=ylabel('Electrode');
0645 set(h,'fontsize',12,'fontweight','bold');
0646 
0647 
0648 if isfield(GND,'bin_info')
0649     h=title(['Bin ' int2str(use_bin) ': ' GND.bin_info(use_bin).bindesc]);
0650 else
0651     h=title(['Bin ' int2str(use_bin) ': ' GND.bindesc{use_bin}]);
0652 end
0653 set(h,'fontsize',14,'fontweight','bold');
0654 
0655 hold on;
0656 v=axis;
0657 %rightmost vertical line to cover up weird white lip that peeks out from
0658 %under imagesc
0659 h=plot([1 1]*v(2),v(3:4));
0660 set(h,'color',[1 1 1]*.6);
0661 
0662 %horizontal lines
0663 show_times=GND.time_pts(show_tpts);
0664 for a=1:n_wind,
0665     start_tpt=find_tpt(GND.t_tests(test_id).time_wind(a,1),show_times);
0666     stop_tpt=find_tpt(GND.t_tests(test_id).time_wind(a,2),show_times);
0667     for c=0:(n_use_chans+n_blank-1),
0668         h=plot([start_tpt-.5 stop_tpt+.5],[1 1]*c+.5);
0669         if ismember(c+1,skipped) || ismember(c,skipped),
0670             %make lines near boundaries different??  currently not used
0671             set(h,'color',[0 0 0],'buttondownfcn',bdfcn,'userdata',dat);
0672         else
0673             set(h,'color',[0 0 0],'buttondownfcn',bdfcn,'userdata',dat);
0674         end
0675     end
0676 end
0677 
0678 %vertical lines
0679 if plot_vert_lines && ~strcmpi(GND.t_tests(test_id).mean_wind,'yes'),
0680     for t=1:(n_show_tpts),
0681         if ismember(show_tpts(t),use_tpts),
0682             if ismember(t+1,xtick) || ismember(t,xtick),
0683                 %draw vertical lines at time tick marks differently??
0684                 %currently option isn't used
0685                 c=[0 0 0];
0686             else
0687                 c=[0 0 0];
0688             end
0689             if isempty(skipped),
0690                 h=plot([1 1]*t+.5,v(3:4));
0691                 set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0692                 if t==1 || ~ismember(show_tpts(t-1),use_tpts),
0693                     h=plot([1 1]*t-.5,v(3:4));
0694                     set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0695                 end
0696             else
0697                 h=plot([1 1]*t+.5,[v(3) skipped(1)-.5]);
0698                 set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0699                 if length(skipped)==2,
0700                     h=plot([1 1]*t+.5,[skipped(1)+.5 skipped(2)-.5]);
0701                     set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0702                 end
0703                 h=plot([1 1]*t+.5,[skipped(end)+.5 v(4)]);
0704                 set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0705                 
0706                 if t==1 || ~ismember(show_tpts(t-1),use_tpts),
0707                     h=plot([1 1]*t-.5,[v(3) skipped(1)-.5]);
0708                     set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0709                     if length(skipped)==2,
0710                         h=plot([1 1]*t-.5,[skipped(1)+.5 skipped(2)-.5]);
0711                         set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0712                     end
0713                     h=plot([1 1]*t-.5,[skipped(end)+.5 v(4)]);
0714                     set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0715                 end
0716             end
0717         end
0718     end
0719 else
0720     c=[0 0 0];
0721     for w=1:n_wind,
0722         tt(1)=find_tpt(show_times,GND.t_tests(test_id).time_wind(w,1));
0723         tt(2)=find_tpt(show_times,GND.t_tests(test_id).time_wind(w,2));
0724         ct=-1;
0725         for t=tt,
0726             ct=ct+1;
0727             if isempty(skipped),
0728                 if ct,
0729                     t_plt=[1 1]*t+.5;
0730                 else
0731                      t_plt=[1 1]*t-.5;
0732                 end
0733                 h=plot(t_plt,v(3:4));
0734                 set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0735             else
0736                 if ct,
0737                     t_plt=[1 1]*t+.5;
0738                 else
0739                     t_plt=[1 1]*t-.5;
0740                 end
0741                 
0742                 h=plot(t_plt,[v(3) skipped(1)-.5]);
0743                 set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0744                 if length(skipped)==2,
0745                     h=plot(t_plt,[skipped(1)+.5 skipped(2)-.5]);
0746                     set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0747                 end
0748                 h=plot(t_plt,[skipped(end)+.5 v(4)]);
0749                 set(h,'color',c,'buttondownfcn',bdfcn,'userdata',dat);
0750             end
0751         end
0752     end
0753 end
0754 
0755 
0756 %% Add Colorbar (if using RGB)
0757 if strcmpi(use_color,'rgb')
0758     cbar_ax=colorbar;
0759     axes(cbar_ax);
0760     hold on;
0761     v=axis;
0762     dlt=2*abs_mx/length(colormap);
0763     hf=fill([v(1) v(1) v(2) v(2)],[0 dlt dlt 0],zero_color); % this fills in the gap in the colorbar that is made to keep grey for non-significant cells
0764     set(hf,'linestyle','none');
0765     v=axis;
0766     axis(v);
0767     rngX=v(2)-v(1);
0768     rngY=v(4)-v(3);
0769     if strcmpi(p.Results.units,'t')
0770         ht=text(v(2)+rngX*.8,rngY*.005,'t');
0771         %ht=ylabel('t');
0772         set(ht,'fontsize',14,'rotation',0,'fontweight','bold', ...
0773             'horizontalalignment','left','verticalalignment','middle');
0774     else
0775         ht=text(v(2)+rngX*.8,rngY*.005,'\muV');
0776         %ht=ylabel('t');
0777         set(ht,'fontsize',14,'rotation',0,'fontweight','bold', ...
0778             'horizontalalignment','left','verticalalignment','middle');
0779     end
0780     h_ax(2)=cbar_ax;
0781 end
0782 
0783 %
0784 %% %%%%%%%%%%%%%%%%%%%%% function orderofmag() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0785 %
0786 function ord=orderofmag(val)
0787 %function ord=orderofmag(val)
0788 %
0789 % Returns the order of magnitude of the value of 'val' in multiples of 10
0790 % (e.g., 10^-1, 10^0, 10^1, 10^2, etc ...)
0791 % used for computing erpimage trial axis tick labels as an alternative for
0792 % plotting sorting variable
0793 
0794 val=abs(val);
0795 if val>=1
0796     ord=1;
0797     val=floor(val/10);
0798     while val>=1,
0799         ord=ord*10;
0800         val=floor(val/10);
0801     end
0802     return;
0803 else
0804     ord=1/10;
0805     val=val*10;
0806     while val<1,
0807         ord=ord/10;
0808         val=val*10;
0809     end
0810     return;
0811 end
0812 
0813 
0814 %
0815 %% %%%%%%%%%%%%%%%%%%%%% function str2bool() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0816 %
0817 function bool=str2bool(str)
0818 %function bool=str2bool(str)
0819 
0820 if ischar(str),
0821     if strcmpi(str,'yes') || strcmpi(str,'y')
0822         bool=1;
0823     else
0824         bool=0;
0825     end
0826 else
0827    bool=str; 
0828 end