run_tutorial

 MarsBaR batch script to show off MarsBaR batching
 This script replicates most of the work in the MarsBaR tutorial.
 See https://marsbar-toolbox.github.io

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