ui_plot

 method plots data in various formats
 FORMAT r = ui_plot(o, plot_spec, plot_params)

 Input
 o            - marsy object

 plot_spec    - can be a
                string, in which case it has the same meaning as
                    the structure version below, with only field 'types'
                    set to this value, OR 
                structure, with none ore more of fields
                (in all cases, field missing is the same as the field
                being empty)
                r_nos      - region number(s) 
                             (if empty, -> all regions)
                b_nos      - session/subject numbers 
                             (if empty, -> all blocks)  
                types      - string, or cell arrays of strings
                             specifying plot type(s).  Plot types can be
                             one or more of: 
                  'raw'    - plots raw time series 
                  'acf     - plots autocorrelation function
                  'fft'    - plots fourier analysis of time series
                  'all'    - all of the above 
                  (defaults to 'raw')
                graphics_obj - graphics object(s) for plot.  
                            If empty, becomes handle SPM graphics
                            window, otherwise, can be a handle to
                            figure, or enough handles to axes for all
                            requested plots.  This is useful for putting
                            the plot on a convenient figure or axis 
   
 plot_params  - (optional) Parameters to pass to plots
                Can be empty (giving defaults) or structure, or
                cell array of structures, one per requested plot
                Relevant fields of structure differ for different plot
                types;
                Plot 'fft': fields 'bin_length' distance between time
                            points in seconds
                     'acf': fields 'lags' no of lags to plot [10]

 Returns
 r            - results that were plotted; all the data for 'raw', all
                the acf values for 'acf', all the fft coefficients for
                'fft'.  Cell array, one cell per plot.

 Examples
 ui_plot(Y, 'all')  plots all plot types for all regions and sessions 
 ui_plot(Y, 'acf')  just plots ACF, for all regions
 ui_plot(Y, struct('types', 'fft', 'r_nos', 1))
                    will plot fft for region 1 only
 f = figure; ui_plot(Y, struct('graphics_obj', f))
                    will plot all plot types to the new figure
 ui_plot(Y, 'acf', struct('lags', 8))
                    plots ACF for all regions, with lag of 8 rather than
                    the default of 10

 $Id$ 

0001 function r = ui_plot(o, plot_spec, plot_params)
0002 % method plots data in various formats
0003 % FORMAT r = ui_plot(o, plot_spec, plot_params)
0004 %
0005 % Input
0006 % o            - marsy object
0007 %
0008 % plot_spec    - can be a
0009 %                string, in which case it has the same meaning as
0010 %                    the structure version below, with only field 'types'
0011 %                    set to this value, OR
0012 %                structure, with none ore more of fields
0013 %                (in all cases, field missing is the same as the field
0014 %                being empty)
0015 %                r_nos      - region number(s)
0016 %                             (if empty, -> all regions)
0017 %                b_nos      - session/subject numbers
0018 %                             (if empty, -> all blocks)
0019 %                types      - string, or cell arrays of strings
0020 %                             specifying plot type(s).  Plot types can be
0021 %                             one or more of:
0022 %                  'raw'    - plots raw time series
0023 %                  'acf     - plots autocorrelation function
0024 %                  'fft'    - plots fourier analysis of time series
0025 %                  'all'    - all of the above
0026 %                  (defaults to 'raw')
0027 %                graphics_obj - graphics object(s) for plot.
0028 %                            If empty, becomes handle SPM graphics
0029 %                            window, otherwise, can be a handle to
0030 %                            figure, or enough handles to axes for all
0031 %                            requested plots.  This is useful for putting
0032 %                            the plot on a convenient figure or axis
0033 %
0034 % plot_params  - (optional) Parameters to pass to plots
0035 %                Can be empty (giving defaults) or structure, or
0036 %                cell array of structures, one per requested plot
0037 %                Relevant fields of structure differ for different plot
0038 %                types;
0039 %                Plot 'fft': fields 'bin_length' distance between time
0040 %                            points in seconds
0041 %                     'acf': fields 'lags' no of lags to plot [10]
0042 %
0043 % Returns
0044 % r            - results that were plotted; all the data for 'raw', all
0045 %                the acf values for 'acf', all the fft coefficients for
0046 %                'fft'.  Cell array, one cell per plot.
0047 %
0048 % Examples
0049 % ui_plot(Y, 'all')  plots all plot types for all regions and sessions
0050 % ui_plot(Y, 'acf')  just plots ACF, for all regions
0051 % ui_plot(Y, struct('types', 'fft', 'r_nos', 1))
0052 %                    will plot fft for region 1 only
0053 % f = figure; ui_plot(Y, struct('graphics_obj', f))
0054 %                    will plot all plot types to the new figure
0055 % ui_plot(Y, 'acf', struct('lags', 8))
0056 %                    plots ACF for all regions, with lag of 8 rather than
0057 %                    the default of 10
0058 %
0059 % $Id$
0060   
0061 % Get, check data from object
0062 [n_rows n_cols] = summary_size(o);
0063 if ~prod([n_rows n_cols]), warning('No data to plot'), return, end
0064 rows = block_rows(o);
0065 Y    = summary_data(o);
0066 N    = region_name(o);
0067 S    = sumfunc(o);
0068 info = summary_info(o);
0069 if strcmp(S, 'unknown'), S = ''; end
0070 if ~isempty(S), S = [' - ' S]; end
0071 
0072 def_spec = struct('types', {{'raw'}}, ...
0073           'r_nos', 1:n_cols, ...
0074           'b_nos', 1:length(rows),...
0075           'graphics_obj', []);
0076 
0077 % Process plot type arguments
0078 if nargin < 2
0079   plot_spec = [];
0080 end
0081 if ischar(plot_spec), plot_spec = struct('types', plot_spec); end
0082 plot_spec = mars_struct('fillafromb', plot_spec, def_spec);
0083 
0084 plot_types = plot_spec.types;
0085 if ~iscell(plot_types)
0086   plot_types = {plot_types};
0087 end
0088 if strcmp('all', plot_types{1})
0089   plot_types = {'raw','acf','fft'};
0090 end
0091 n_p_t   = length(plot_types);
0092 n_plots = n_p_t * length(plot_spec.r_nos); 
0093 
0094 % Check passed graphics object(s)
0095 gr_ob = plot_spec.graphics_obj;
0096 if isempty(gr_ob)
0097   gr_ob = spm_figure('GetWin','Graphics');
0098   spm_results_ui('Clear', gr_ob, 0);
0099 end
0100 if ~all(ishandle(gr_ob))
0101   error('One or more graphics objects are not valid');
0102 end
0103 o_t = get(gr_ob, 'type');
0104 if iscell(o_t)
0105   if any(diff(strvcat(o_t{:}),[],1))
0106     error('All graphics objects must be of the same type');
0107   end
0108   o_t = o_t{1};
0109 end
0110 switch o_t
0111  case 'figure'
0112   figure(gr_ob);
0113  case 'axes'
0114   if n_plots > length(gr_ob)
0115     error('Not enough axes for planned number of plots');
0116   end
0117  otherwise
0118   error(['Completely confused by graphics object type: ' o_t]);
0119 end
0120 
0121 rows     = rows(plot_spec.b_nos);
0122 n_blocks = length(rows);
0123 
0124 % Process plot_param arguments
0125 if nargin < 3
0126   plot_params = cell(1,n_p_t);
0127 end
0128 if ~iscell(plot_params)
0129   plot_params = {plot_params};
0130 end
0131 if length(plot_params) == 1
0132   plot_params = repmat(plot_params, 1, n_p_t);
0133 elseif length(plot_params) < n_p_t
0134   error(sprintf('You need %d plot_param entries', n_p_t));
0135 end
0136  
0137 % Default string, bin_length for fft plots
0138 if mars_struct('isthere', info, 'TR')
0139   bin_length = info.TR;
0140   bin_str = 'Hz';    
0141 else
0142   bin_length = 1;
0143   bin_str = 'cycles per time point';
0144 end
0145 
0146 p_ctr = 1;
0147 for c = plot_spec.r_nos
0148   for p = 1:n_p_t
0149     switch o_t
0150      case 'figure'
0151       subplot(n_plots, 1, p_ctr); 
0152      case 'axes'
0153       axes(gr_ob(p_ctr));
0154     end
0155 
0156     y = Y(:,c);
0157     
0158     switch lower(plot_types{p})
0159      case 'raw'
0160       a_r = [];
0161       for s = 1:n_blocks
0162     a_r = [a_r; rows{s}(:)];
0163       end
0164       plot(y(a_r));
0165       axis tight
0166       hold on
0167       % add session indicators
0168       yrg = [min(y) max(y)];
0169       for s = 1:n_blocks
0170     rw = rows{s};
0171     if s > 1 % session divider
0172       plot([rw(1) rw(1)]-0.5, yrg, 'k:');
0173     end
0174       end
0175       ylabel(['Signal intensity' S])
0176       xlabel('Time point');
0177       r{p_ctr} = y;
0178       
0179      case 'acf'
0180       if isfield(plot_params{p}, 'lags')
0181     lags = plot_params{p}.lags;
0182       else lags = 10;
0183       end
0184       mn_len = Inf;
0185       for s = 1:n_blocks
0186     if length(rows{s}) < mn_len, mn_len = length(rows{s}); end
0187       end
0188       lags = max([1 min([mn_len-4 lags])]);
0189 
0190       C = [];
0191       for s = 1:n_blocks
0192     by = y(rows{s});
0193     nb_rows = length(by);
0194     ty = toeplitz(by, [by(1) zeros(1, lags)]);
0195     ty([1:lags n_rows+1:end], :) = [];
0196     Cs      = corrcoef(ty);
0197     C       = [C Cs(1,:)];
0198     n       = nb_rows - 2;                   % df for correlation
0199     t_th    = spm_invTcdf(1-0.025, n);       % t for two tailed p=0.05
0200     r_th(s) = sqrt(1/(n/t_th.^2+1));         % r equivalent
0201       end
0202       stem(C);
0203       axis tight
0204       hold on
0205       % add confidence intervals, labels, session indicators
0206       Crg = [min(C) max(C)];
0207       for s = 1:n_blocks
0208     s_vals = [0 lags] + (s-1)*(lags+1) + 1;
0209     plot(s_vals,  [r_th(s) r_th(s)], 'r:');
0210     plot(s_vals, -[r_th(s) r_th(s)], 'r:');
0211     if s > 1 % session divider line
0212       plot([s_vals(1) s_vals(1)]-0.5, Crg, 'k:');
0213     end
0214       end
0215       xtl = get(gca,'xticklabel');
0216       xtl = num2str(mod(str2num(xtl)-1, lags+1));
0217       set(gca, 'xticklabel', xtl);
0218 
0219       ylabel('Correlation coefficient')
0220       xlabel('Lag');
0221       r{p_ctr} = C;
0222      case 'fft'
0223       if isfield(plot_params{p}, 'bin_length')
0224     b_len = plot_params{p}.bin_length;
0225     b_str = 'Hz';
0226       else
0227     b_len = bin_length;
0228     b_str = bin_str;
0229       end
0230       P = []; H = []; St = [];
0231       for s = 1:n_blocks
0232     by    = y(rows{s});
0233     gX    = abs(fft(by)).^2;
0234     gX    = gX*diag(1./sum(gX));
0235     q     = size(gX,1);
0236     Hz    = [0:(q - 1)]/(q * b_len);
0237     q     = 2:fix(q/2);
0238     P     = [P gX(q)'];
0239     St(s) = length(H)+1;
0240     H     = [H; Hz(q)'];
0241       end
0242       plot(P)
0243       hold on
0244       Prg = [min(P) max(P)];
0245       for s = 2:n_blocks
0246     plot([St(s) St(s)]-0.5, Prg, 'k-');
0247       end
0248       axis tight
0249       % Rename tick labels
0250       xt  = get(gca, 'xtick');
0251       xt  =  xt(xt==fix(xt));
0252       set(gca, 'xtick', xt);
0253       for t = 1:length(xt)
0254     xtl_fft{t} = sprintf('%5.3f', H(xt(t)));
0255       end
0256       set(gca, 'xticklabel', xtl_fft);
0257       xlabel(sprintf('Frequency (%s)', b_str))
0258       ylabel('Relative spectral density')
0259       axis tight
0260       r{p_ctr} = P;
0261      otherwise
0262       error(['What is this plot type: ' plot_types{p} '?']);
0263     end
0264     title(N{c}, 'interpreter', 'none');
0265   
0266     p_ctr = p_ctr + 1;
0267   end
0268 end
0269 
0270 return
0271