run_tutorial

 MarsBaR batch script to show off MarsBaR batching
 This script replicates most of the work in the MarsBaR tutorial.
 See http://marsbar.sourceforge.net

 The script assumes that the current directory is the 'batch' directory
 of the example data set files.

 $Id: run_tutorial.m,v 1.3 2004/08/15 01:19:43 matthewbrett Exp $

0001 % MarsBaR batch script to show off MarsBaR batching
0002 % This script replicates most of the work in the MarsBaR tutorial.
0003 % See http://marsbar.sourceforge.net
0004 %
0005 % The script assumes that the current directory is the 'batch' directory
0006 % of the example data set files.
0007 %
0008 % $Id: run_tutorial.m,v 1.3 2004/08/15 01:19:43 matthewbrett Exp $
0009 
0010 % Start marsbar to make sure spm_get works
0011 marsbar('on')
0012 
0013 % You might want to define the path to the example data here, as in
0014 % root_dir = '/my/path/somewhere';
0015 root_dir = spm_get(-1, '', 'Root directory of example data');
0016 sesses = {'sess1','sess2','sess3'};
0017 
0018 % Directory to store (and load) ROIs
0019 roi_dir = fullfile(root_dir, 'rois');
0020 
0021 % MarsBaR version check
0022 if isempty(which('marsbar'))
0023   error('Need MarsBaR on the path');
0024 end
0025 v = str2num(marsbar('ver'));
0026 if v < 0.35
0027   error('Batch script only works for MarsBaR >= 0.35');
0028 end
0029 marsbar('on');  % needed to set paths etc
0030 
0031 % SPM version check. We need this to guess which model directory to use and
0032 % to get the SPM configured design name.
0033 spm_ver = spm('ver');
0034 sdirname = [spm_ver '_ana'];
0035 if strcmp(spm_ver, 'SPM99') 
0036   conf_design_name = 'SPMcfg.mat';
0037 else
0038   spm_defaults;
0039   conf_design_name = 'SPM.mat';
0040 end
0041 
0042 % Set up the SPM defaults, just in case
0043 spm('defaults', 'fmri');
0044 
0045 % Subdirectory for reconfigured design
0046 mars_sdir = 'Mars_ana';
0047 
0048 % Get SPM model for session 2
0049 sess2_model_dir = fullfile(root_dir, 'sess2', sdirname);
0050 sess2_model = mardo(fullfile(sess2_model_dir, 'SPM.mat'));
0051 if ~is_spm_estimated(sess2_model)
0052   error(['Session 2 model has not been estimated. ' ...
0053     'You may need to run the run_preprocess script']);
0054 end
0055 
0056 % Get activation cluster by loading T image
0057 con_name = 'stim_hrf';
0058 t_con = get_contrast_by_name(sess2_model, con_name);
0059 if isempty(t_con)
0060   error(['Cannot find the contrast ' con_name ...
0061     ' in the design; has it been estimated?']);
0062 end
0063 % SPM2 stores contrasts as vols, SPM99 as filenames
0064 if isstruct(t_con.Vspm)
0065   t_con_fname = t_con.Vspm.fname;
0066 else
0067   t_con_fname = t_con.Vspm;
0068 end
0069 t_img_name = fullfile(sess2_model_dir, t_con_fname);
0070 if ~exist(t_img_name)
0071   error(['Cannot find t image ' t_img_name ...
0072      '; has it been estimated?']);
0073 end
0074 
0075 % Get t threshold of uncorrected p < 0.05
0076 p_thresh = 0.05;
0077 erdf = error_df(sess2_model);
0078 t_thresh = spm_invTcdf(1-p_thresh, erdf);
0079 
0080 % get all voxels from t image above threshold
0081 V = spm_vol(t_img_name);
0082 img = spm_read_vols(V);
0083 tmp = find(img(:) > t_thresh);
0084 img = img(tmp);
0085 XYZ = mars_utils('e2xyz', tmp, V.dim(1:3));
0086 
0087 % make into clusters, find max cluster
0088 cluster_nos = spm_clusters(XYZ);
0089 [mx max_index] = max(img);
0090 max_cluster = cluster_nos(max_index);
0091 cluster_XYZ = XYZ(:, cluster_nos == max_cluster);
0092 
0093 % Make ROI from max cluster
0094 act_roi = maroi_pointlist(struct('XYZ', cluster_XYZ, ...
0095                  'mat', V.mat), 'vox');
0096 
0097 % Make box ROI to do trimming
0098 box_limits = [-20 20; -66 -106; -20 7]';
0099 box_centre = mean(box_limits);
0100 box_widths = abs(diff(box_limits));
0101 box_roi = maroi_box(struct('centre', box_centre, ...
0102                'widths', box_widths));
0103 
0104 % Combine for trimmed ROI
0105 trim_stim = box_roi & act_roi;
0106 
0107 % Give it a name
0108 trim_stim = label(trim_stim, 'batch_trim_stim');
0109 
0110 % save ROI to MarsBaR ROI file, in current directory, just to show how
0111 saveroi(trim_stim, fullfile(roi_dir, 'batch_trim_stim_roi.mat'));
0112 
0113 % Save as image
0114 save_as_image(trim_stim, fullfile(root_dir, 'batch_trim_stim.img'));
0115 
0116 % We will do estimation for the trimmed functional ROI, and for two
0117 % anatomical ROIs
0118 bg_L_name = fullfile(roi_dir, 'MNI_Putamen_L_roi.mat');
0119 bg_R_name = fullfile(roi_dir, 'MNI_Putamen_R_roi.mat');
0120 roi_array{1} = trim_stim;
0121 roi_array{2} = maroi(bg_L_name);
0122 roi_array{3} = maroi(bg_R_name);
0123 
0124 % MarsBaR estimation for sessions 1 and 3 follows
0125 pwd_orig = pwd;
0126 sesses = {'sess1', 'sess3'};
0127 
0128 % event specification for getting fitted event time-courses
0129 % A bit silly here, as we only have one session per model and one event type
0130 event_session_no = 1;
0131 event_type_no = 1;
0132 event_spec = [event_session_no; event_type_no];
0133 event_duration = 0; % default SPM event duration
0134 
0135 clear model_file
0136 for roi_no = 1:length(roi_array)
0137   roi = roi_array{roi_no};
0138   for ss = 1:length(sesses)
0139     
0140     % Run SPM model configuration, just to show we don't need to do SPM
0141     % estimation before using MarsBaR
0142     % We only need to do this for the first ROI, because we reuse the
0143     % design for each ROI
0144     if roi_no == 1
0145       model_file{ss} = configure_er_model(root_dir, sesses{ss}, mars_sdir);
0146     end
0147     
0148     D = mardo(model_file{ss});
0149     
0150     % Extract data
0151     Y = get_marsy(roi, D, 'mean');
0152     
0153     % MarsBaR estimation
0154     E = estimate(D, Y);
0155     
0156     % Add contrast, return model, and contrast index
0157     [E Ic] = add_contrasts(E, 'stim_hrf', 'T', [1 0 0]);
0158     
0159     % Get, store statistics
0160     stat_struct(ss) = compute_contrasts(E, Ic);
0161     
0162     % And fitted time courses
0163     [this_tc dt] = event_fitted(E, event_spec, event_duration);
0164 
0165     % Make fitted time course into ~% signal change
0166     tc(:, ss) = this_tc / block_means(E) * 100;
0167     
0168   end
0169 
0170   % Show calculated t statistics and contrast values
0171   % NB this next line only works when we have only one stat/contrast
0172   % value per analysis
0173   vals = [ [1 3]; [stat_struct(:).con]; [stat_struct(:).stat]; ];
0174   fprintf('Statistics for %s\n', label(roi));
0175   fprintf('Session %d; contrast value %5.4f; t stat %5.4f\n', vals);
0176   
0177   % Show fitted event time courses
0178   figure
0179   secs = [0:length(tc) - 1] * dt; 
0180   plot(secs, tc)
0181   title(['Time courses for ' label(roi)], 'Interpreter', 'none');
0182   xlabel('Seconds')
0183   ylabel('% signal change');
0184   legend(sesses)
0185   
0186 end