attachment:cbu_meeg_spm5_pipeline_5-3-09.m of SpmDemo - Meg Wiki
location: attachment:cbu_meeg_spm5_pipeline_5-3-09.m of SpmDemo

Attachment 'cbu_meeg_spm5_pipeline_5-3-09.m'

Download

   1 % Batch script (linear pipeline) for MEG analysis in SPM5   
   2 % Rik Henson MRC CBU Sep 08
   3 %
   4 % Note that would benefit from modularisation (eg for stop/start
   5 % running), but wanted all steps in one file to ease exposition
   6 %
   7 % First start SPM with (at the CBU) "spm 5 eeg" from Linux prompt
   8 
   9 clear
  10 
  11 wd = '/imaging/rh01/Methods/CBUMEGDemo/Group';   % !!!REPLACE WITH YOUR DIRECTORY!!!
  12 cd(wd)
  13 
  14 addpath /imaging/local/spm/spm5/cbu_updates/    
  15 % Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/spm5_cbu_updates/devel/
  16 addpath /imaging/local/meg_misc/    
  17 % Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/meg_misc/devel/
  18 
  19 %%%%%%%%%% Different types of sensor
  20 typ{1}   = 'mags'; typ{2} = 'grds'; typ{3} = 'eeg';
  21 Ntyp     = length(typ);
  22 eegflag  = [0 0 1];         % Binary flag for which typ is EEG
  23 
  24 % Num channels per type (if differs across subjects, can re-specify below)
  25 Nchan(1) = 102;    Nchan(2) = 204;   Nchan(3) = 70;
  26 
  27 % Any user-thresholds for rejection for each of above (Inf=none)
  28 thr(1)   = Inf;    thr(2)   = Inf;   thr(3)   = 120;
  29 EOGthr   = 120;   % EOG threshold (EOG is duplicated in each file)
  30 
  31 %%%%%%%%%% Experiment-specific parameters
  32 expwin      = [-100 400];                    % Epoch window (in ms)
  33 expcodes    = [1:8];                         % Codes in FIF file to extract
  34 offset      = [ones(1,length(expcodes))*34]; % t=0 offset (eg visual delay)
  35 newcodes    = [1 1 1 1 2 2 2 2];             % Any new (recoding) of above
  36 expcons     = [1 0; 0 1; 1 -1; 1/2 1/2];     % Contrasts (on NEW codes)
  37 Ncons       = size(expcons,1); 
  38 		 
  39 %%%%%%%%%% Raw FIF files on network 
  40 % (cell of cells for each session of each sub (only one session here))
  41 % Eg if instead two sessions per subject:
  42 %    Fifdata = { { '/megdata/cbu/lfp/meg08_0014/080117/block1_raw.fif';
  43 %                  '/megdata/cbu/lfp/meg08_0014/080117/block2_raw.fif'} };
  44 fifdata = {
  45            {'/megdata/cbu/lfp/meg08_0014/080117/block4_raw.fif'};
  46            {'/megdata/cbu/lfp/meg08_0015/080118/block4_raw.fif'};
  47            {'/megdata/cbu/lfp/meg08_0016/080118/block4_raw.fif'};
  48            {'/megdata/cbu/lfp/meg08_0038/080131/block4_raw.fif'};
  49            {'/megdata/cbu/lfp/meg08_0064/080212/block4_raw.fif'};
  50            {'/megdata/cbu/lfp/meg08_0111/080303/block4_raw.fif'};
  51            {'/megdata/cbu/lfp/meg08_0112/080303/block4_raw.fif'};
  52            {'/megdata/cbu/lfp/meg08_0113/080303/block4_raw.fif'};
  53            {'/megdata/cbu/lfp/meg08_0127/080310/block4_raw.fif'};
  54 	     }
  55 
  56 
  57 %%%%%%%%%% Raw DICOM files on network (one per subject)
  58 mridata = {
  59            '/mridata/cbu/CBU080769_CBU080769/20080916_095501/Series_002_CBU_MPRAGE';
  60            '/mridata/cbu/CBU071255_CBU071255/20071218_150851/Series_002_CBU_MPRAGE';
  61            '/mridata/cbu/CBU071151_CBU071151/20071116_140044/Series_002_CBU_MPRAGE';
  62            '/mridata/cbu/CBU071229_CBU071229/20071211_112249/Series_002_CBU_MPRAGE';
  63            '/mridata/cbu/CBU080034_RIS13667/20080111_102845/Series_002_CBU_MPRAGE';
  64            '/mridata/cbu/CBU080292_CBU080292/20080418_161843/Series_002_CBU_MPRAGE';
  65            '/mridata/cbu/CBU080358_CBU080358/20080506_151233/Series_003_CBU_MPRAGE';
  66            '/mridata/cbu/CBU080413_CBU080413/20080522_124313/Series_002_CBU_MPRAGE';
  67            '/mridata/cbu/CBU080802_CBU080802/20080930_090538/Series_002_CBU_MPRAGE';
  68 	  };
  69 
  70 
  71 %%%%%%%%%% Bad channel names noticed by Operator or post hoc (eeg/meg,sub,ses)
  72 badchans{1,1,1} = [1143];
  73 badchans{1,2,1} = [1143 0913 1013 0313 0933 2222 0642 1232 1913];
  74 badchans{1,3,1} = [1143];
  75 badchans{1,4,1} = [1143 2223 0432];
  76 badchans{1,5,1} = [1143];
  77 badchans{1,6,1} = [1143];
  78 badchans{1,7,1} = [1143 2223];
  79 badchans{1,8,1} = [1143 1412];
  80 badchans{1,9,1} = [1143 2223];
  81 
  82 %%%%%%%%%% Bad EEG channels noticed by Operator or post hoc (in terms of NAME)
  83 badchans{2,1,1} = [48 50 73];                 % Poor contact
  84 badchans{2,2,1} = [7 13 16 24 46 48 57];      % Drifting
  85 badchans{2,3,1} = [48];                       % Drifting
  86 badchans{2,4,1} = [48];                       % Drifting
  87 badchans{2,5,1} = [15 40 49];                 % Poor contact
  88 badchans{2,6,1} = [48];                       % Drifting
  89 badchans{2,7,1} = [48 72];                    % No deflection (odd N170 topo)?
  90 badchans{2,8,1} = [65 48];                    % odd in final topo of N170?
  91 badchans{2,9,1} = [20 31 45];                 % 45 poor contact; 20+31 highfreq noise
  92 
  93 % (Just happens that any EEG Channels 61 onwards are numbered 65 onwards in our lab!)
  94 for sub=1:size(badchans,2)
  95   for ses=1:size(badchans,3)
  96    f = find(badchans{2,sub,ses}>60);
  97    badchans{2,sub,ses}(f) = badchans{2,sub,ses}(f)-4;
  98   end
  99 end
 100 
 101 %%%%%%%%%%% subject-specific parameters
 102 %dosubs   = [1:9];      % Subjects to run, eg for subset: dosubs = [5:8]; 
 103 dosubs   = [1:8];       % Exclude sub9 because MRI segmentation requires manual re-
 104                         % positioning first (can try yourself; see MRI section below!)
 105                         % (in fact for many, eg sub5 too, normalisation is much
 106                         % better after manual repositioning)
 107 Nsub     = length(dosubs);
 108 
 109 % Below assumes that two EOG channels (VEOG and HEOG), which will be
 110 % appended at end of every data-file. Note that this script does not yet
 111 % handle concurrent ECG (see Jason Taylor for how to do this)
 112 for s = 1:Nsub
 113  subnum = dosubs(s);
 114  for m = 1:Ntyp;
 115    Nsubchan(s,m) = Nchan(m);  % In case different number of channels (eg EEG) for some subjects
 116   chan_thr{s,m} = [ones(1,Nsubchan(s,m))*thr(m) EOGthr EOGthr]; 
 117  end
 118 end
 119 
 120 VitECapsules = ones(1,max(dosubs));  % Whether MRIs have Vitamin E capsule (all here)
 121 
 122 %%%%%%%%%%% control flags (1=do step; 0=omit step)
 123 maxfiltflag = 1; maxmvcompflag = 1; maxtransflag = 1; 
 124 usemaxmvcompflag = 0; usemaxtransflag = 1; 
 125 chancheckflag = 0; respmflag = 1; reprocessflag = 1; 
 126 getmriflag = 1;
 127 write2Dflag = [1 0 1]; write3Dflag = 1;
 128 grpinvertflag = 1;     % Whether invert using group priors (after individual inversions)
 129                        % (Litvak & Friston, 2008, Neuroimage)
 130 fusionflag = 1;        % Whether want to simultaneously invert (fuse) modalities (typs)
 131                        % (Henson et al, submitted)
 132 		       
 133 %%%%%%%%%% General parameters
 134 % Maxfilter
 135 downsample_maxf   = '-ds 4';           % maxfilter downsampling factor (eg 4 to save space)
 136                                        % though beware of effect of downsampling on triggers
 137 				       % (http://imaging.mrc-cbu.cam.ac.uk/meg/CbuSpmParameters)
 138 %downsample_maxf  = '';                % (if none)
 139 format_maxf       = '-format float';   % maxfilter data format (native = float32)
 140 %format_maxf      = '-format short';   % (if to save space)
 141 trans_offset_maxf = [0 -13 +6];        % New centre for trans default
 142                                        % (see our maxfilter WIKI)
 143 
 144 % SPM
 145 downsample_spm    = 100;                % spm new sample rate (Hz; should be at
 146                                         % least twice filt_cut_spm below)
 147 filt_cut_spm      = 40;                 % Lowpass filter cutoff for SPM (Hz)
 148 dformat_spm       = 'int16';            % data format for spm (eg to save space)
 149 sensor_smooth_spm = [5 5 10];           % Smoothing in SensorSpace-Time [mm mm ms];
 150 source_smooth_spm = [10 10 10];         % Smoothing in Source space [mm mm mm]
 151 
 152 %%%%%%%%%%% Initialisation
 153 subsphere = [];
 154 allsssbad = cell(Nsub,1);
 155 nbadchan = cell(Ntyp,Nsub);
 156 nrejects = cell(Ntyp,Nsub);
 157 nevents  = cell(Ntyp,Nsub);
 158 
 159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 160 %%%%%%%%%%% Preprocess MEG/EEG data
 161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 162 
 163 for sub = 1:Nsub
 164  
 165  subnum    = dosubs(sub);
 166  Nses(sub) = size(fifdata{subnum},1);
 167 
 168  cd(wd)
 169  swd = sprintf('sub%d',subnum);
 170  try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end
 171   
 172  for ses = 1:Nses(sub)
 173 
 174    allsssbad{sub} = cell(1,Nses);
 175    
 176    fnm_stem = sprintf('sub%d_ses%d',subnum,ses);
 177 
 178    if maxfiltflag
 179 
 180 %%%%% Use autobad to detect bad channels from first 20s, and select from output
 181      rawfile  = fifdata{subnum}{ses};
 182      sssfile  = sprintf('%s_sss.fif',fnm_stem)
 183      logfile  = sprintf('%s_autobad.log',fnm_stem)
 184      headptsfile  = sprintf('%s_headpoints.txt',fnm_stem)
 185 
 186      if exist(sssfile), delete(sssfile), end
 187      if exist(logfile), delete(logfile), end
 188      [Centre,Radius]  = meg_fit_sphere(rawfile,pwd,headptsfile);
 189      disp(sprintf('Subject %d (%d): Session %d: Sphere centre [%d %d %d] and radius %d mm',subnum,sub,ses,Centre(1),Centre(2),Centre(3),Radius))
 190      subsphere(sub,:) = [Centre Radius];
 191 
 192      disp('Checking first 20s for bad channels...')
 193      eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d -autobad 20 -skip 21 9999999 -v | tee %s',rawfile,sssfile,Centre(1),Centre(2),Centre(3),logfile))
 194      
 195      badfile = sprintf('%s_badchan.txt',fnm_stem)
 196      if exist(badfile), delete(badfile), end
 197      eval(sprintf('!cat %s | sed -n ''/Static/p'' | cut -f 5- -d '' '' > %s',logfile,badfile));
 198      sssbadchans{sub}{ses} = unique([textread(badfile,'%d')' badchans{1,subnum,ses}]);
 199      for bc=1:length(sssbadchans{sub}{ses}); allsssbad{sub}{ses} = [allsssbad{sub}{ses} sprintf(' %04d',sssbadchans{sub}{ses}(bc))]; end
 200 
 201 %%%%% Now mark bad channels and run ST 
 202      logfile = sprintf('%s_ssst.log',fnm_stem)
 203      if exist(sssfile), delete(sssfile), end    % Don't need previous file from autobad
 204      sssfile  = sprintf('%s_ssst.fif',fnm_stem)
 205      if exist(logfile), delete(logfile), end
 206      posfile = sprintf('%s_mvpos.txt',fnm_stem)
 207 
 208      disp(['Now maxfiltering with bad channels (and tSSS): ' allsssbad{sub}{ses}])
 209      eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -headpos -hpistep 1000 -hpisubt amp -hp %s %s -v | tee %s',rawfile,sssfile,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},posfile,downsample_maxf,logfile))
 210 
 211      S.logfile = logfile;
 212      [max_d(sub),max_r(sub)] = meg_plot_maxfilter_move(S);
 213 % Do check the movement parameters - particularly goodness of fit -
 214 % because movement compensation below will be invalidated if poor fit
 215 % (also, the 'inter' MaxFilter option is used below in case the HPI 
 216 % temporarily fails (which it does for sub5), but if your subject's 
 217 % HPI failed  permanently from the start, you might not know this otherwise,
 218 % because 'inter' would handle it fine!).
 219      
 220      infile = sssfile;
 221      
 222 %%%%% Now mark bad channels and run MVCOMP
 223      if maxmvcompflag
 224 
 225       sss2file = sprintf('%s_mvcomp.fif',fnm_stem)
 226       logfile = sprintf('%s_mvcomp.log',fnm_stem)
 227       if exist(sss2file), delete(sss2file), end
 228       if exist(logfile), delete(logfile), end
 229 
 230       disp(['Now compensating movement... (every 1000ms)'])
 231       eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -movecomp inter -hpistep 200 -hpisubt amp %s -v | tee %s',rawfile,sss2file,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},downsample_maxf,logfile))
 232 %      eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -movecomp -hpistep 1000 -hpisubt amp %s -v | tee %s',rawfile,sss2file,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},downsample_maxf,logfile))
 233 
 234       infile = sss2file;
 235      end
 236      
 237 %%%%% Now Trans Default (if requested)
 238      if maxtransflag
 239       trans_origin = Centre(1:3) + trans_offset_maxf;
 240       logfile = sprintf('%s_trans.log',fnm_stem)
 241       sss3file = sprintf('%s_trans.fif',fnm_stem)
 242       if exist(sss3file), delete(sss3file), end
 243       if exist(logfile), delete(logfile), end
 244 
 245       disp(['Now trans-ing to DEFAULT ' num2str(round(trans_origin))])
 246       eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -autobad off -trans default -frame head -origin %d %d %d %s -v -force | tee %s',infile,sss3file,trans_origin(1),trans_origin(2),trans_origin(3),format_maxf,logfile))
 247      end
 248     end
 249 
 250 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 251 %%%%% Olaf's channel checking (have not used for a while; should put in SVN!)
 252     if chancheckflag 
 253      addpath /imaging/olaf/MEG/artefact_scan_tool/
 254      global fiffs tasks values criteria sensors datapara figs;
 255      ini_check4artefacts; % Initialisation of variables
 256      epochlength = 5;    
 257      amplitude_threshold_mag = 5000;  
 258      amplitude_threshold_gra = 10000;         
 259      out_dir = fullfile(wd,'ChanCheck');
 260      addfiff(infile,fnm_stem);
 261      addtask('maxminamp');   % maximum minus minimum amplitude within epoch
 262      addtask('maxmingrad');  % maximum signal gradient (amplitude difference between successive data points in time) with epoch
 263      addtask('threshold');   % number of epochs with max-min amplitudes above specified threshold
 264      addtask('crosscorr');   % intercorrelation matrix for magnetometers and gradiometers separately
 265      addtask('trendamp');    % linear trend in max-min differences across all epochs
 266      addtask('trendgrad');   % linear trend in maximum signal gradients across all epochs
 267      check4artefacts;    % That's where it all happens... lean back and enjoy!
 268     end
 269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 270 
 271 %%%%% Read into SPM (writes separate files for MEG and EEG)
 272 
 273     clear S
 274     if usemaxtransflag
 275       S.Fdata = sprintf('%s_trans.fif',fnm_stem);
 276       if maxmvcompflag
 277        S.HPIfile = sprintf('%s_mvcomp.fif',fnm_stem);   % If trans'ed (hpi lost)
 278       else
 279        S.HPIfile = sprintf('%s_sss.fif',fnm_stem);      % If trans'ed (hpi lost)
 280       end
 281     elseif usemaxmvcompflag
 282       S.Fdata = sprintf('%s_mvcomp.fif',fnm_stem);
 283       S.HPIfile = [];
 284     else
 285       S.Fdata = sprintf('%s_sss.fif',fnm_stem);      
 286       S.HPIfile = [];
 287     end
 288 %    meg_plot_fifpoints(S);                              % If want to check points
 289     S.Fchannels = 'FIF306_setup.mat';
 290     S.veogchan  = [62 63]; S.veog_unipolar = 0;
 291     S.heogchan  = [61 64]; S.heog_unipolar = 0;
 292     S.trig_chan = 'STI101';
 293     S.twin      = [0 Inf];
 294     S.eeg_ref   = 'average';
 295     S.trig_type = 'positive_from_zero';
 296 %    S.Dtype     = 'float32';                             % make int16 if files too big
 297     S.Fchannels_eeg = fullfile(spm('dir'),'EEGtemplates', sprintf('cbu_meg_%deeg_montage_old.mat',Nsubchan(sub,find(eegflag)))); % use default montage (which simply explains how to display in 2D), prior to new CBU caps in mid-2008
 298 %    S.Fchannels_eeg =  fullfile(wd,swd,sprintf('eeg_montage_sub%02d.mat',subnum)); % use subject-specific (care labels may not be standard)
 299     
 300     if respmflag
 301       D = spm_eeg_rdata_FIF(S);
 302       basefile = D.fname(1:(end-4));
 303     else
 304       basefile = S.Fdata(1:(end-4));
 305     end
 306 
 307 %%%%% Downsample, then split Mags and Grads (file too big otherwise)
 308 
 309     if reprocessflag
 310       clear S; S.D = sprintf('%s-eeg.mat',basefile); disp(S.D)
 311       S.Radc_new = downsample_spm;
 312       D = spm_eeg_downsample(S);
 313 
 314       clear S; S.D = basefile;
 315       S.Radc_new = downsample_spm;
 316       D = spm_eeg_downsample(S);
 317       basefile = D.fname(1:(end-4));
 318 
 319       clear S; S.D = D.fname;
 320       D = spm_eeg_splitFIF(S);
 321     end
 322     
 323 %%%%% If want to view rawdata (with Danny's tool):
 324 %     meg_viewdata(D.fname);
 325      
 326 %%%%% Preprocess each data-type
 327 
 328      for typn = 1:Ntyp
 329  
 330        if reprocessflag
 331 	clear S
 332         S.D = sprintf('%s-%s.mat',basefile,typ{typn}); disp(S.D)
 333 	S.filter.type = 'butterworth';
 334 	S.filter.band = 'low';	S.filter.PHz = filt_cut_spm;
 335 %	S.filter.band = 'stop'; S.filter.PHz = [49 51];
 336         S.Dtype       = dformat_spm;
 337 	D = spm_eeg_filter(S);
 338 
 339 %%%%%% !!!Insert ICA here if want to remove EOG/ECG artifacts!!!
 340 %%%%%% (See Jason Taylor; assume filtered, downsampled data appropriate)
 341 	
 342 	clear S; S.D = D.fname;
 343  	S.events.start    = expwin(1);
 344  	S.events.stop     = expwin(2);
 345  	S.events.types    = expcodes;
 346  	S.events.Inewlist = 0;
 347  	S.events.resynch  = offset;
 348 	D = spm_eeg_epochs(S);
 349 	
 350 	% re-classify event-codes
 351 	for e = 1:length(expcodes);
 352 	    fi = find(D.events.code==expcodes(e));
 353 	    if ~isempty(fi)
 354 		D.events.code(fi) = newcodes(e);
 355 	    end
 356 	end
 357 	D.events.types  = unique(D.events.code);
 358 	D.events.Ntypes = length(D.events.types);
 359 	save(D.fname,'D')
 360 	clear S; S.D = D.fname;
 361 
 362       else
 363         S.D = sprintf('e_fd%s-%s.mat',basefile,typ{typn}); disp(S.D)
 364 	D = spm_eeg_ldata(S.D);
 365       end
 366       
 367 %%%%% Threshold channels (eg EOG)
 368 
 369       nc = D.Nchannels;
 370       S.thresholds.External_list   = 0;
 371       S.thresholds.Check_Threshold = 1;
 372       S.thresholds.threshold       = chan_thr{sub,typn};
 373       S.BadThr                     = 0.2;  
 374       S.artefact.weighted          = 0;
 375       if eegflag(typn)==0
 376 %           S.channels.Bad = unique([D.channels.Bad strcmpi(D.channels.name....sssbadchans{sub}{ses}]);     	% If dont trust MaxFilter channel reconstruction
 377       else
 378             S.channels.Bad = unique([D.channels.Bad badchans{2,subnum,ses}]);  
 379       end
 380       D = spm_eeg_artefact(S);
 381       
 382       nbadchan{typn,sub} = [nbadchan{typn,sub} length(D.channels.Bad)];	
 383       nrejects{typn,sub} = [nrejects{typn,sub} length(find(D.events.reject))];
 384 
 385 %%%%% Re-reference EEG if removed channels
 386       if eegflag(typn) & ~isempty(D.channels.Bad)
 387 	   D = spm_eeg_ldata(D.fname);
 388            meg_reavg_ref(D);
 389 	   save(D.fname,'D')
 390       end
 391       sesnam{typn,ses} = D.fname;
 392       
 393 %%%%% Uncomment if want to average and contrast each session separately 
 394 %%%%% (as well as after merging below)
 395 %
 396 %      clear S; S.D = D.fname;
 397 %      D = spm_eeg_average(S);
 398 %
 399 %      clear S; S.D = D.fname;
 400 %      S.c         = expcons;
 401 %      S.WeightAve = 1;
 402 %      D = spm_eeg_weight_epochs(S);
 403 %      nevents{typn,sub} = D.events.repl;
 404       
 405     end  % end of typ loop
 406   end  % end of session loop
 407     
 408 %%%%% Merge sessions, average, contrast 
 409 
 410   for typn = 1:Ntyp
 411 
 412     if Nses>1
 413       if ~maxtransflag
 414         disp(['Concatenating sessions... (NB: merged files will not have correct channel Loc/Orient info unless trans-ed to first or to default!!)']);
 415       end
 416       clear S; S.D = strvcat(sesnam{typn,:});
 417       S.recode = cell(1,Nses);
 418       D = spm_eeg_merge(S);
 419       clear S; S.D = D.fname;
 420     else
 421       clear S; S.D = sesnam{typn,1};
 422     end
 423 
 424     D = spm_eeg_average(S);
 425 
 426     clear S; S.D = D.fname;
 427     S.c = expcons;
 428     S.WeightAve = 0;
 429     D = spm_eeg_weight_epochs(S);
 430     
 431     nevents{typn,sub}  = D.events.repl;
 432     finalnam{typn,sub} = fullfile(D.path,D.fname);
 433   end
 434 end
 435 
 436 cd(wd)
 437 
 438 for typn = 1:Ntyp
 439  disp(sprintf('%s...',typ{typn}))
 440  for sub = 1:Nsub
 441    subnum = dosubs(sub);
 442    disp(sprintf('\tSub %d (%d)...',subnum,sub))
 443    disp(sprintf('\t\tNevents (per condition): %s',mat2str(nevents{typn,sub})))
 444    disp(sprintf('\t\tRejects (per session): %s',mat2str(nrejects{typn,sub})))
 445    disp(sprintf('\t\tBadchans (per session): %s',mat2str(nbadchan{typn,sub})))
 446  end
 447 end
 448 
 449 %%%% Create RMS of Grads (mainly for Sensor-Time image files)
 450 
 451 for typn = 1:Ntyp
 452   if strcmp(typ{typn},'grds')
 453      for sub = 1:Nsub
 454        cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
 455        clear S; S.D = finalnam{typn,sub};
 456        S.do_write = 1;
 457        D = meg_grds2grms(S); % RMS of grads; just for sensor-level analyses 
 458        finalnam{Ntyp+1,sub} = fullfile(D.path,D.fname);
 459        nevents{Ntyp+1,sub}  = nevents{typn,sub};
 460      end
 461      Ntyp=Ntyp+1; typ{Ntyp} = 'grms';
 462      write2Dflag(Ntyp) = 1;
 463   end
 464 end
 465 
 466 %%%% Write out Sensor-Time image files for each subject
 467 
 468 for typn = find(write2Dflag)
 469   for sub = 1:Nsub
 470          cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
 471 	 disp(finalnam{typn,sub})
 472          [pth,nam,ext] = fileparts(finalnam{typn,sub});
 473          clear S; S.Fname = [nam ext]; 
 474 	 S.interpolate_bad = 0;
 475          S.n = 32;
 476          S.pixsize = 3;
 477          S.trialtypes = [1:Ncons];
 478          P = spm_eeg_convertmat2ana3D(S);
 479 	 if any(sensor_smooth_spm)   
 480 	  if length(sensor_smooth_spm)==1; sensor_smooth_spm = ones(1,3)*sensor_smooth_spm; end
 481 	  for con = 1:size(P,1);
 482 	   [pth,nam,ext] = fileparts(P(con,:));
 483 	   Pout          = fullfile(pth,['s' nam ext]);
 484            spm_smooth(spm_vol(P(con,:)),Pout,sensor_smooth_spm);
 485 	   Pin = strvcat(P(con,:),Pout);
 486 	   spm_imcalc_ui(Pin,Pout,'((i1+eps).*i2)./(i1+eps)',{[],[],'float32',0});   %Reinsert NaNs
 487  	   SensorTimeImgs{typn}{sub}(con,:) = fullfile(wd,swd,Pout);
 488 	  end
 489 	 else
 490 	  for con = 1:size(P,1);
 491  	   SensorTimeImgs{typn}{sub}(con,:) = fullfile(wd,swd,P(con,:));
 492 	  end
 493 	 end
 494   end
 495 end
 496 
 497 %%%% Create Grand Mean across subjects (inc RMS (grms))
 498 cd(wd)
 499 for typn = 1:Ntyp
 500  clear S; S.P = strvcat(finalnam{typn,:});
 501  S.Pout = sprintf('G%d-%s.mat',Nsub,typ{typn})
 502  D = spm_eeg_grandmean(S);
 503 end
 504 
 505 
 506 %return
 507 
 508 
 509 eval(sprintf('save preproc%d.mat subsphere nrejects nbadchan nevents allsssbad dosubs',length(dosubs)))
 510 
 511 disp('!!!NOW USE SPM-Display-M/EEG button to examine results, eg G8-*.mat files !!!')
 512 
 513 % Once displaying one of the modalities, you can hold shift and select both 
 514 % conditions 1+2, and click on a posterior right channel to see the M170
 515 % difference between faces (1) and scrambled faces (2). You can also select
 516 % condition 3, which is the differential ERF/ERP between faces and scrambled,
 517 % then press "topography", enter -100ms for initial time, then press and hold
 518 % the time slider to see topography of difference at every time point (displayed
 519 % in Matlab window) - random noise until ~100ms, when strong differences emerge
 520 
 521 
 522 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 523 %%%%%%%%%%%%%%%%% 2D sensor x 1D Time SPMs 
 524 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 525 
 526 stdir = fullfile(wd,'SensorTimeSPMs');
 527 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
 528 
 529 clear grpsubs grpcons imgfiles
 530 grpsubs{1} = [1:Nsub];
 531 % If want to split subjects, eg into two groups, eg:
 532 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
 533 grpcons{1} = [1:2];  % Include each condition for mags
 534 grpcons{3} = [1:2];  % Include each condition for eeg
 535 grpcons{4} = [1:2];  % Diff of RMS (be careful if diff number of trials for diff conditions)
 536 %grpcons{4} = [3];    % Or RMS of diff (not equivalent to above)
 537 
 538 Ngrp = length(grpsubs);
 539 for typn = find(write2Dflag)
 540  outdir = fullfile(stdir,typ{typn})
 541  for grp = 1:Ngrp
 542   for sub = 1:length(grpsubs{grp});
 543    subnum = grpsubs{grp}(sub);
 544    imgfiles{grp}{sub} = SensorTimeImgs{typn}{subnum}(grpcons{typn},:);
 545   end
 546  end
 547  meg_batch_anova(imgfiles,outdir);
 548 end
 549 
 550 disp('!!!NOW press Results button in SPM!!!')
 551 
 552 %return
 553 
 554 % You need to know a bit about the SPM Results interface, but press the Results 
 555 % button (make sure SPM is in "EEG" mode) (or better still, use the "ns" toolbox
 556 % under "Toolbox" button - see WIKI page below) and select the SPM.mat file for 
 557 % one of the modalities, select the default effects of interest F-contrast
 558 % (because we don't care about polarity, at least for Mags and EEG), choose
 559 % 0.001 uncorrected. Then press "Volume" to get table of stats (ignore brain image). 
 560 % You should see, at least in Mags and EEG, frontal and posterior clusters 
 561 % (simply of different polarity) maximal around 170ms. Not much survives whole-
 562 % image correction because only 8 subjects, so few df's ie low power and RFT 
 563 % conservative) - though effects around 170ms do for corrected extent-level 
 564 % for Mags and EEG (using "ns" toolbox); GRMS results are poor though. 
 565 % Correction for height would also be significant if you used SVC because
 566 % of an a priori timewindow of interest around 170ms (again see WIKI page below).
 567 %
 568 % For more info on interpreting these SPMs, see
 569 %                   http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm
 570 %
 571 % Also, for an example use with EEG data, see Henson et al (2008), Neuroimage.
 572 %
 573 % Note that you can do SVC for a priori timewindows of interest (see end of above WIKI)
 574 % But one use of this SPMs is to define timewindows of interest (in the absence of 
 575 % a priori ones) for the subsequent source localisation below...
 576 %
 577 % (Note that these SPMs can also handle different numbers and locations of channels
 578 % across subjects, provided their locations are digitised in same space, because the 
 579 % data are interpolated to same grid. For MEG data, particularly gradiometers, the
 580 % results benefit from using 'trans -default' above, to realign across different
 581 % headpositions for different subjects - see Smith et al (in press), HBM Abstract.
 582 
 583 
 584 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 585 %%%%%%%%%%%%%%%%% Get MRIS
 586 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 587 
 588 for sub=1:Nsub 
 589    subnum = dosubs(sub);
 590    swd = sprintf('sub%d',subnum);
 591    mwd{subnum} = fullfile(wd,swd,'MRI');
 592    try cd(mwd{subnum}); catch eval(sprintf('!mkdir %s',mwd{subnum})); cd(mwd{subnum}); end
 593 
 594    if getmriflag
 595      hdr = spm_select('FPList',mridata{subnum},'^*.dcm');
 596      hdr = spm_dicom_headers(hdr);
 597      
 598      patsex{subnum} = hdr{1}.PatientsSex;
 599      patage{subnum} = hdr{1}.PatientsAge;   % NB At time of scanning, not time of MEG!
 600      
 601      spm_dicom_convert(hdr,'all','flat','nii');
 602 
 603      mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-01.nii$');
 604      mrifile{sub} = strtrim(mrifile{sub}(1,:));   % Take first if multiple previous
 605 
 606      disp('!!!GOOD IDEA TO MANUALLY REPOSITION STRUCTURAL VIA DISPLAY BUTTON!!!')
 607 %     pause
 608      
 609      if VitECapsules(subnum)
 610       meg_headmask(mrifile{sub},'display');         % Remove any VitE capsules
 611       disp('!!!Check difference image!!!');
 612 %      pause     
 613       eval(sprintf('!mv %s_masked.nii %s',mrifile{sub}(1:end-4),mrifile{sub}))
 614      end
 615      patage
 616      patsex
 617    
 618    else
 619      mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-01.nii$');
 620      mrifile{sub} = strtrim(mrifile{sub}(1,:));   % Take first if multiple previous     
 621      wmrifile{sub} = spm_select('FPList',mwd{subnum},'^wms.*-01.nii$');
 622      wmrifile{sub} = strtrim(wmrifile{sub}(1,:));   % Take first if multiple previous     
 623    end
 624 end
 625 spm_check_registration(strvcat(mrifile));    % Only feasible for <~12 images
 626 spm_check_registration(strvcat(wmrifile));   % (sub5 bit too warped at back)
 627 
 628 disp('!!!NOW DISPLAY MRI IMAGE in SPM AND DETERMINE FIDUCIALS!!!')
 629 
 630 %%%%%%%%%% Coordinates of fiducials from above MRIs (nasion, left, right)
 631 % Below are from using MRI as is (without manually reorienting) - very approximate
 632 smri_fids{1}  = [ 3 111  22; -78 18 -31; 77 11 -33];
 633 smri_fids{2}  = [ 6 107 -28; -69 26 -41; 68 14 -45];
 634 smri_fids{3}  = [-2 124 -23; -78 22 -52; 70 21 -58];
 635 smri_fids{4}  = [ 3 109 -8;  -72 18 -41; 72 18 -47];
 636 smri_fids{5}  = [-15 103 -33; -79 3 -38; 62 13 -46];
 637 smri_fids{6}  = [-4 109 -20; -71 15 -32; 63 11 -35];
 638 smri_fids{7}  = [-1 93  -5;  -72 14 -44; 70 10 -44];
 639 smri_fids{8}  = [ 1 99  -22; -73  8 -27; 69 10 -42];
 640 smri_fids{9}  = [-3 125 -24; -76 22 -55; 69 22 -58];
 641 % Below are after my manual reorienting (just for your reference)
 642 %smri_fids{1}  = [ 0.7 80.0 -31.2;   -76.9 -14.3 -45.9;   79.9 -19.2 -41.0];
 643 %smri_fids{2}  = [ 1.3 78.7 -30.9;   -71.8  -2.5 -58.0;   71.5   1.8 -58.0];  
 644 %smri_fids{3}  = [ 1.7 88.8 -24.6;   -74.0  -1.4 -66.1;   71.7  -2.8 -66.1];
 645 %smri_fids{4}  = [-1.6 78.9 -22.4;   -74.8 -13.3 -57.7;   71.6  -7.2 -54.6];
 646 %smri_fids{5} = {TBA};
 647 %smri_fids{6}  = [-2.4 79.1 -27.1;   -72.1  -9.8 -44.6;   64.3  -9.8 -51.9];
 648 %smri_fids{7}  = [-2.3 73.2  -7.0;   -71.8  -2.0 -46.9;   70.2  -4.9 -49.8]; 
 649 %smri_fids{8}  = [ 0.0 78.1 -30.3;   -71.3  -1.5 -48.5;   72.8   6.0 -51.5]; 
 650 %smri_fids{9}  = [ 0.0 88.9 -31.1;   -76.8 -15.5 -56.6;   78.0 -12.7 -55.0];
 651 
 652 
 653 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 654 %%%%%%%%%%%%%%%%% Localise on Mesh
 655 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 656 
 657 addpath /imaging/local/Brainstorm/Toolbox
 658 
 659 clear typ;
 660 typ{1}   = 'mags'; typ{2} = 'grds'; typ{3} = 'eeg';
 661 Ntyp     = length(typ);
 662 		      
 663 % !!! There are several options for source localisation, eg choice of
 664 % forward model (eg Spheres or BEMs), and choice of inversion method (eg
 665 % MSP, COH, MNM). Below illustrates just two possible inversions (each 
 666 % inversion is indexed by "val"), the second of which has been commented 
 667 % out because BEMs take a VERY LONG time (and don't work well for EEG yet). 
 668 % Of course many more options can be tried (and compared using the model-evidence,
 669 % recorded in mod_evds below, provided that you don't also change the
 670 % data, eg epoch inverted). 
 671 %
 672 % For more information about forward models, see 
 673 %    http://imaging.mrc-cbu.cam.ac.uk/meg/SpmForwardModels
 674 % For more information about inversion parameters, see 
 675 %    Friston et al (2008) Neuroimage 
 676 %    (PDF available here http://imaging.mrc-cbu.cam.ac.uk/meg/SpmAnalysis)
 677 
 678 %%%%%%%%%%%%%%%%%%%%
 679 val = 1;  % Spherical forward model using canonical cortical mesh but individual head meshes
 680 comment{val} = 'Concentric Spheres';
 681 cortex{val}  = 'Can';   % Canonical (inverse-normalised template) cortical mesh
 682 CtxSize(val) = 8196;    % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
 683 iskull{val}  = 'Ind';   % Individual (SPM-created) inner skull mesh
 684 oskull{val}  = '';      % (Outer skull irrelevant for spherical models)
 685 scalp{val}   = 'Ind';   % Individual (SPM-created) scalp mesh
 686 formod{val}  = 'Sph';   % Concentric spheres
 687 surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all data-types
 688 
 689 invcons{val} = [1 2];   % Conditions to invert
 690 invtype{val} = 'MSP';   % Type of inversion
 691 invtwin{val} = expwin;  % Invert whole epoch
 692 
 693 %%%%% Options for subsequent time-freq contrast
 694 contwin{val}{1} = [140 190];  % Time window for contrast
 695 contwin{val}{2} = [];         % Frequency window for contrast ([]=all freqs) 
 696 conengy{val}    = 'evoked';   % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
 697 
 698 %%%%%%%%%%%%%%%%%%%%
 699 %val = 2;  % BEM using canonical meshes (Warning: BEMs take a long time to compute!)
 700 %
 701 %comment{val} = 'Boundary Element Model';
 702 %cortex{val}  = 'Can';    % Canonical (inverse-normalised template) cortical mesh
 703 %CtxSize(val) = 8196;     % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
 704 %iskull{val}  = 'Can';    % Canonical (inverse-normalised template) inner skull mesh
 705 %oskull{val}  = 'Can';    % Canonical (inverse-normalised template) outer skull mesh
 706 %scalp{val}   = 'Can';    % Canonical (inverse-normalised template) scalp mesh
 707 %formod{val}  = 'BEM';    % Boundary Element Model
 708 %surfit(val,:)= [2 2 4];  % Surfaces for BEM: up to 2nd for MEG (iskull) and 4th for EEG (+oskull+scalp)
 709 %invcons{val} = [1 2];     % Conditions to invert
 710 %invtype{val} = 'MSP';     % Type of inversion
 711 %invtwin{val} = expwin;    % Invert whole epoch
 712 
 713 %%%%% Options for subsequent time-freq contrast
 714 %contwin{val}{1} = [140 190];  % Time window for contrast
 715 %contwin{val}{2} = [];         % Frequency window for contrast ([]=all freqs) 
 716 %conengy{val}    = 'evoked';   % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
 717 
 718 % Other options you can consider are "overlapping spheres" for the
 719 % forward model (see spm_eeg_inv_BSTfwdsol.m), minimumum norm (IID) for
 720 % the inversion, etc. See Henson et al (2009), Neurimage, for more options.
 721 
 722 Nval = length(comment);
 723 
 724 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 725 
 726 for sub = 1:Nsub
 727 
 728  cd(wd)
 729  subnum = dosubs(sub);
 730  swd = sprintf('sub%d',subnum);
 731  cd(swd)
 732  
 733  for val = 1:Nval
 734    
 735   for typn = 1:Ntyp
 736 
 737    clear S;  
 738    S.D = char(finalnam{typn,sub});
 739 % If use finalnam, then will be localising evoked energy only; if want to
 740 % localise total energy (evoked and induced), then select (c)ae_* file instead
 741 %     e.g S.D = sprintf('ae_fdsub%d_ses1_trans-%s.mat',subnum,typ{typn}) 
 742 % Also may prefer "mvcomp" rather than "trans" data, if don't trust "trans default"
 743 
 744    D = spm_eeg_ldata(S.D); disp(S.D)
 745 
 746 %If want to clear previous localisations (safer?)
 747    if isfield(D,'inv');  D = rmfield(D,'inv'); save(D.fname,'D'); end
 748    if isfield(D,'con');  D = rmfield(D,'con'); save(D.fname,'D'); end
 749 
 750 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise
 751 
 752    disp(comment{val})
 753    D.val = val; D.fname
 754    D.inv{val}.method = 'imaging';
 755    load('defaults_eeg_inv.mat')
 756    D.inv{val} = invdefts.imag;
 757    D.inv{val}.date = mat2str(fix(clock));
 758    D.inv{val}.comment = comment{val};
 759 
 760 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI
 761 
 762 %%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI
 763 
 764    D.inv{val}.mesh.sMRI = mrifile{sub};   % In case many from previous
 765 
 766    [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);
 767    mri_stem = nam;
 768    wmrifile{sub} = fullfile(pth,['wm' nam ext]);
 769    if ~exist(wmrifile{sub})
 770     D = spm_eeg_inv_spatnorm(D);	
 771    else
 772     def_name      = [nam '_sn_1.mat'];
 773     isndef_name   = [nam '_inv_sn_1.mat'];
 774     D.inv{val}.mesh.def    = fullfile(pth,def_name);
 775     D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
 776     D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
 777     D.inv{val}.mesh.wmMRI  = fullfile(pth,['wm' nam ext]);
 778    end
 779 
 780 %%%%%%%%%%%%%%%%%%% Meshing 2: Creating binary images for scalp (and inner skull)
 781 
 782    D.inv{val}.mesh.msk_scalp = spm_select('List',mwd{subnum},'^*scalp\.nii$');
 783  
 784    if isempty(D.inv{val}.mesh.msk_scalp)
 785     D = spm_eeg_inv_getmasks(D);	
 786    else	% assume other masks exist too!
 787     D.inv{val}.mesh.msk_scalp   = fullfile(mwd{subnum},D.inv{val}.mesh.msk_scalp);
 788     D.inv{val}.mesh.msk_iskull  = spm_select('FPList',mwd{subnum},'^*iskull\.nii$');
 789     D.inv{val}.mesh.msk_cortex  = spm_select('FPList',mwd{subnum},'^*cortex\.nii$');
 790     D.inv{val}.mesh.msk_flags   = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
 791    end
 792    spm_check_registration(char(D.inv{val}.mesh.msk_scalp,D.inv{val}.mesh.msk_iskull,D.inv{val}.mesh.msk_cortex,D.inv{val}.mesh.sMRI));
 793 
 794 %%%%%%%%%%%%%%%%%%% Meshing 2: Creating meshes
 795 
 796    D.inv{val}.mesh.template  = 0;
 797    D.inv{val}.mesh.canonical = 1;
 798    D.inv{val}.mesh.Msize = CtxSize(val);
 799   
 800   % Prompt for canonical or user-specified cortical meshs (SPM does not currently
 801   % produce individual ones, because automatic meshing of cortex is difficult!)
 802  
 803    if strcmp(cortex{val},'Can')
 804     Tmesh = load(fullfile(spm('dir'),'EEGtemplates',sprintf('mni_mesh_cortex_%d.mat',CtxSize(val))));
 805     D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
 806     D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
 807     Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 808     D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
 809     D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
 810    elseif strcmp(cortex{val},'Load')
 811     Tmesh = spm_select(1,'mat','Select cortex mat file containing vert and face',[],pwd,'.*mat');
 812     D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
 813     D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
 814     Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.invdef);
 815     D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
 816     D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
 817    else
 818     error('No cortical mesh specified????')
 819    end
 820 
 821    iskull_mesh_name = fullfile('MRI',[mri_stem '_iskull_mesh.mat']);
 822    scalp_mesh_name = fullfile('MRI',[mri_stem '_scalp_mesh.mat']);
 823     
 824   % Prompt for canonical, individual (from SPM) or user-specified inner skull
 825   % (Use canonical if worried about canonical cortex piercing innner skull)
 826 
 827    if strcmp(iskull{val},'Can')
 828     Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_iskull_2562.mat'));
 829     Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 830     D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
 831     D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
 832    elseif strcmp(iskull{val},'Ind')
 833     if ~exist(iskull_mesh_name)
 834      D = spm_eeg_inv_getmeshes(D,1);
 835      vert = D.inv{val}.mesh.tess_iskull.vert;
 836      face = D.inv{val}.mesh.tess_iskull.face;
 837      save(iskull_mesh_name,'vert','face');
 838     else
 839      Tmesh = load(iskull_mesh_name);
 840      D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
 841      D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
 842     end
 843    elseif strcmp(iskull{val},'Load')
 844      Tmesh = spm_select(1,'mat','Select iskull mat file containing vert and face',[],pwd,'.*mat');
 845      D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
 846      D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
 847    else
 848      warning('No inner skull mesh specified')
 849    end
 850 
 851   % Prompt for canonical or user-specified outer skull (SPM does not currently
 852   % produce individual ones, because segmenting outer skull from scalp difficult)
 853   % (Note that outer skull only really necessary for EEG BEMs)
 854 
 855    if strcmp(oskull{val},'Can')
 856     Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_oskull_2562.mat'));
 857     Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 858     D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
 859     D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
 860    elseif strcmp(oskull{val},'Load')
 861      Tmesh = spm_select(1,'mat','Selectoskull  mat file containing vert and face',[],pwd,'.*mat');
 862      D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
 863      D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
 864    else
 865      warning('No outer skull mesh specified')
 866    end
 867 
 868   % Prompt for canonical, individual (from SPM) or user-specified scalp
 869   % (Use canonical if worried about canonical cortex or inner skull piercing scalp)
 870 
 871    if strcmp(scalp{val},'Can')
 872     Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_scalp_2562.mat'));
 873     Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 874     D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
 875     D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
 876    elseif strcmp(scalp{val},'Ind')
 877     if ~exist(scalp_mesh_name)
 878      D = spm_eeg_inv_getmeshes(D,2);
 879      vert = D.inv{val}.mesh.tess_scalp.vert;
 880      face = D.inv{val}.mesh.tess_scalp.face;
 881      save(scalp_mesh_name,'vert','face');
 882     else
 883      Tmesh = load(scalp_mesh_name);
 884      D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
 885      D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
 886     end
 887    elseif strcmp(scalp{val},'Load')
 888      Tmesh = spm_select(1,'mat','Select scalp mat file containing vert and face',[],pwd,'.*mat');
 889      D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
 890      D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
 891    else
 892      warning('No scalp mesh specified')
 893    end
 894 
 895   % Always individual scalp for datareg
 896   
 897    if ~exist(scalp_mesh_name)
 898      S = spm_eeg_inv_getmeshes(D,2);   % lastarg==1 returns iskull (not scalp)
 899      vert = S.inv{val}.mesh.tess_scalp.vert;
 900      face = S.inv{val}.mesh.tess_scalp.face;
 901      save(scalp_mesh_name,'vert','face');
 902      D.inv{val}.mesh.tess_scalp_datareg.vert  = vert;
 903      D.inv{val}.mesh.tess_scalp_datareg.face  = uint16(face);  
 904    else
 905      Tmesh = load(scalp_mesh_name);
 906      D.inv{val}.mesh.tess_scalp_datareg.vert  = Tmesh.vert;
 907      D.inv{val}.mesh.tess_scalp_datareg.face  = uint16(Tmesh.face);  
 908    end
 909   
 910    spm_eeg_inv_checkmeshes(D);
 911    save(D.fname,'D');
 912 
 913 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg
 914 
 915    if isempty(D.inv{val}.datareg.fid_coreg) 
 916     
 917     D.inv{val}.datareg.sensors    = D.channels.Loc';
 918     D.inv{val}.datareg.fid_eeg    = D.channels.fid_eeg;
 919     D.inv{val}.datareg.headshape  = D.channels.headshape;
 920     if ~eegflag(typn)
 921      D.inv{val}.datareg.megorient = D.channels.Orient';
 922     end
 923     D.inv{val}.datareg.scalpvert  = D.inv{val}.mesh.tess_scalp_datareg.vert;
 924     D.inv{val}.datareg.fid_mri    = smri_fids{subnum};
 925     D = spm_eeg_inv_datareg(D);
 926     fid_err(sub) = sqrt(mean(mean((D.inv{val}.datareg.fid_coreg - D.inv{val}.datareg.fid_mri).^2)))
 927    end
 928   
 929    spm_eeg_inv_checkdatareg(D);
 930    save(D.fname,'D');
 931    clear S; S.D=D.fname; 
 932 %   meg_pretty_datareg(S);        % Prettier rendering of head and sensors
 933 
 934 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model
 935 
 936    rotate3d off
 937    D.inv{val}.method = 'Imaging';
 938    if ~eegflag(typn)
 939     if strcmp(formod{val},'Sph'), meth = 'meg_sphere';
 940     elseif strcmp(formod{val},'BEM'), meth = 'meg_bem'; 
 941     elseif strcmp(formod{val},'OS'), meth = 'meg_os'; 
 942     else, error('Unknown forward model!'); end
 943    else
 944     if strcmp(formod{val},'Sph'), meth = 'eeg_3sphereBerg';
 945     elseif strcmp(formod{val},'BEM'), meth = 'eeg_bem'; 
 946     elseif strcmp(formod{val},'OS'), meth = 'eeg_os'; 
 947     else, error('Unknown forward model!'); end
 948    end
 949    D.inv{val}.forward.method =  meth;
 950    D.inv{val}.forward.surface2fit = surfit(val,typn);
 951 %   D.inv{val}.forward.NVertMax = 1000;   % If want to save some time!
 952    gain_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatrix_' upper(typ{typn}) ...
 953  		    '_' formod{val} '_' num2str(surfit(val,typn)) '.mat'])
 954    gxyz_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatxyz_' upper(typ{typn}) ...
 955 		    '_' formod{val} '_' num2str(surfit(val,typn)) '.mat'])
 956   
 957 %  delete(gain_name);   % if want to be safe 
 958    if ~exist(gain_name) %| ~exist(gxyz_name)
 959     [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
 960     eval(sprintf('!mv %s %s',D.inv{val}.forward.gainmat,gain_name))
 961     eval(sprintf('!mv %s %s',D.inv{val}.forward.gainxyz,gxyz_name))
 962    end
 963 
 964    D.inv{val}.forward.gainmat = gain_name;
 965    D.inv{val}.forward.gainxyz = gxyz_name;
 966    save(D.fname,'D');
 967 
 968 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert each modality separately
 969 
 970    D.inv{val}.inverse = [];
 971    D.inv{val}.inverse.trials = invcons{val};
 972    D.inv{val}.inverse.type   = invtype{val};
 973    D.inv{val}.inverse.woi    = invtwin{val};
 974    D.inv{val}.inverse.lpf    = 0;   % Confusing, but hpf refers to cutoff
 975    D.inv{val}.inverse.hpf    = 40;  % of lowpass (and lpf to cutoff of highpass)
 976 % Some extra parameters that you could play with...(but defaults seem ok!)
 977 %   D.inv{val}.inverse.smooth = 0.8;
 978 %   D.inv{val}.inverse.Np = 256;
 979 %   D.inv{val}.inverse.Han = 0;
 980 
 981    D = spm_eeg_invert(D);
 982 %   D = spm_eeg_invert_fuse(D);   % If want to compare exactly with fused below
 983    save(D.fname,'D');
 984    mod_evds(sub,typn) = D.inv{val}.inverse.F1;
 985    
 986 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
 987 
 988    D.inv{val}.contrast.woi  = contwin{val}{1};
 989    D.inv{val}.contrast.fboi = contwin{val}{2};
 990    D.inv{val}.contrast.type = conengy{val};
 991    D.val = val;
 992    D = spm_eeg_inv_results(D); 
 993    spm_eeg_inv_results_display(D); 
 994 
 995    if write3Dflag
 996     D.inv{val}.contrast.smooth = source_smooth_spm;
 997     D.inv{val}.contrast.rms = 1;
 998 %    D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
 999     D = spm_eeg_inv_Mesh2Voxels(D);
1000     SourceImgs{typn}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1001    end  
1002 
1003 %  If want to save on filesize
1004 %   if val>2
1005 %    canvert = D.inv{val-1}.mesh.tess_mni.vert;
1006 %    canface = D.inv{val-1}.mesh.tess_mni.face;
1007 %    D.inv{val-1}.mesh = [];
1008 %    D.inv{val-1}.mesh.tess_mni.vert = canvert;
1009 %    D.inv{val-1}.mesh.tess_mni.face = canface;
1010 %    D.inv{val-1}.datareg = [];
1011 %    clear canvert canface
1012 %   end 
1013  
1014    disp(sprintf('Sub=%d, %s, val=%d: %s',sub,typ{typn},val,D.inv{val}.comment))
1015    save(D.fname,'D');
1016    
1017   end % end typ
1018 
1019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1020 % Fusion
1021 
1022   if fusionflag
1023    DD = {};
1024    for typn=1:Ntyp
1025      DD{typn} = spm_eeg_ldata(finalnam{typn,sub});
1026    end
1027    DD{1}.inv{DD{1}.val}.inverse.Nm = [];   % Reset any parameters from last inversion
1028    DD{1}.inv{DD{1}.val}.inverse.Nr = [];   % Reset any parameters from last inversion
1029    try rmfield(DD{1}.inv{DD{1}.val}.inverse,'QP'), end
1030 
1031    D = spm_eeg_invert_fuse(DD);  % note necessary so not overwritten by next val
1032 
1033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1034 
1035    D.inv{val}.contrast.woi  = contwin{val}{1};
1036    D.inv{val}.contrast.fboi = contwin{val}{2};
1037    D.inv{val}.contrast.type = conengy{val};
1038    D.val = val;
1039    D = spm_eeg_inv_results(D); 
1040    spm_eeg_inv_results_display(D); 
1041 
1042    if write3Dflag
1043     D.inv{val}.contrast.smooth = source_smooth_spm;
1044     D.inv{val}.contrast.rms = 1;
1045 %    D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
1046     D = spm_eeg_inv_Mesh2Voxels(D);
1047     SourceImgs{Ntyp+1}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1048    end
1049 
1050    typ{Ntyp+1} = 'fused';
1051    finalnam{Ntyp+1,sub} = fullfile(D.path,D.fname);
1052   end
1053   
1054  end % end val
1055  
1056 end % end subs
1057 
1058 for sub = 1:Nsub
1059    subnum = dosubs(sub);
1060    disp(sprintf('\tSub %d (%d)...',subnum,sub))
1061    disp(sprintf('\t\tFiducial error (RMS mm): %3.2f',fid_err(sub)))
1062    for val=1:Nval
1063      disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(mod_evds(sub,:))))
1064    end
1065 end
1066 
1067 % Fid errors may be quite large (~8mm) for some subjects, simply because
1068 % I did not bother to mark them carefully on the MRI. Note that the
1069 % actual coregistration is determined primarily by the digitised headshape;
1070 % the fiducials are only really used to check the quality of registration.
1071 
1072 %spm_check_registration(strvcat(wmrifile))   % double-check normaliation worked
1073 
1074 disp('!!!Review single-subject localisations if you wish!!!')
1075 
1076 % At this point, you can review a single-subject's localisation by
1077 % pressing the "3D source reconstruction" button, and then loading one
1078 % of the mmca* files (for one modality, or the "*_fused.mat" for the
1079 % fusion of the modalities). You can then press the "mip" button to see
1080 % the overall reconstruction, or the "display" button under "window" to
1081 % see the power in the time-freq contrast (for the selected
1082 % condition). You can also display a movie and change the start-stop
1083 % times, etc. Don't press the other buttons for inversion, fusion, group
1084 % inversion etc because they will be covered below...
1085 
1086 %%%%%%%%%%%%%%%%% Or a quick look at mean over subjects...
1087 
1088 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1089 S.fig = 2;
1090 for val=1:Nval
1091  S.val = val; 
1092  for typn=1:(Ntyp+fusionflag)
1093   S.fig = S.fig+1;
1094   S.P = strvcat(finalnam{typn,:});
1095   meg_inv_display_con(S);
1096  end
1097 end
1098 
1099 % You should see bilateral fusiform and anterior inferior temporal
1100 % cortex, plus some more posterior occipitotemporal stuff on the right,
1101 % in most of the modalities, particularly the fused result
1102   
1103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1104 %%%%%%%%%%%%%%%%% 3D  SPMs 
1105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1106 
1107 stdir = fullfile(wd,'SourceSPMs');
1108 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1109 
1110 clear grpsubs imgfiles grpcons
1111 grpsubs{1} = [1:Nsub];
1112 % If want to split subjects, eg into two groups, eg:
1113 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1114 
1115 grpcons = [1:2];
1116 Ngrp = length(grpsubs);
1117 for typn = 1:(Ntyp+fusionflag)
1118  outdir = fullfile(stdir,typ{typn})
1119  for grp = 1:Ngrp
1120   for sub = 1:length(grpsubs{grp});
1121    subnum = grpsubs{grp}(sub);
1122    imgfiles{grp}{sub} = SourceImgs{typn}{subnum}(grpcons,:);
1123   end
1124  end
1125  meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1126 end
1127 
1128 disp('!!!Review the group 3D stats if you wish!!!')
1129 
1130 %return
1131 
1132 % Again you need to know a bit about the SPM Results interface, but press
1133 % Results and select an SPM.mat file for one of the modalities and enter
1134 % a new T-contrast [1 -1] to find voxels in which greater source
1135 % amplitude for faces than scrambled faces, choose 0.001 uncorrected (not
1136 % much survives whole-image correction because only 8 subjects, so few
1137 % df's ie low power and RFT conservative) and then press "Volume" to get
1138 % table of stats. It will be very rough (speckled) because of the low
1139 % dfs, but you should see a cloud around right  fusiform for example for 
1140 % the fused stats. These statistical maps look better with more subjects, 
1141 % or when using SNPM or PPM's (see Henson et al, 2007, Neuroimage and
1142 % Taylor & Henson, in prep, for more information).
1143 
1144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1145 % Group inversion
1146 
1147 if grpinvertflag
1148  for typn=1:Ntyp 
1149   clear DD
1150   
1151   for val=1:Nval
1152    for sub=1:Nsub
1153     DD{sub} = spm_eeg_ldata(finalnam{typn,sub});
1154     DD{sub}.inv{val}.comment    = [DD{sub}.inv{val}.comment 'Group'];
1155     DD{sub}.inv{val}.inverse.Nm = [];   % Reset any parameters from last inversion
1156     DD{sub}.inv{val}.inverse.Nr = [];   % Reset any parameters from last inversion
1157     DD{sub}.val = val;
1158    end
1159    DD = spm_eeg_invert(DD);
1160   end 
1161 
1162   % Save outputs
1163   for sub=1:Nsub
1164     D       = DD{sub};
1165     D.fname = [D.fname(1:(end-4)) '_grp.mat'];
1166     grpfinalnam{typn,sub} = fullfile(D.path, D.fname);
1167     save(fullfile(D.path, D.fname), 'D');
1168     
1169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1170 
1171     D.inv{val}.contrast.woi  = contwin{val}{1};
1172     D.inv{val}.contrast.fboi = contwin{val}{2};
1173     D.inv{val}.contrast.type = conengy{val};
1174     D.val = val;
1175     D = spm_eeg_inv_results(D); 
1176     spm_eeg_inv_results_display(D); 
1177 
1178     if write3Dflag
1179      D.inv{val}.contrast.smooth = source_smooth_spm;
1180      D.inv{val}.contrast.rms = 1;
1181 %     D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
1182      D = spm_eeg_inv_Mesh2Voxels(D);
1183      GrpSourceImgs{typn}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1184     end 
1185     save(fullfile(D.path, D.fname), 'D');
1186   end % end sub
1187  end % end typ
1188  
1189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1190 % Fuse any Group inversion
1191  if fusionflag
1192   for sub=1:Nsub
1193    clear DD;
1194    for typn=1:Ntyp
1195     DD{typn} = spm_eeg_ldata(grpfinalnam{typn,sub});
1196    end
1197    DD{1}.inv{DD{1}.val}.inverse.Nm = [];   % Reset any parameters from last inversion
1198    DD{1}.inv{DD{1}.val}.inverse.Nr = [];   % Reset any parameters from last inversion
1199    % (Important not to reset the source priors from the grp inversion above)
1200    D = spm_eeg_invert_fuse(DD);  % note necessary so not overwritten by next val
1201   
1202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1203 
1204    D.inv{val}.contrast.woi  = contwin{val}{1};
1205    D.inv{val}.contrast.fboi = contwin{val}{2};
1206    D.inv{val}.contrast.type = conengy{val};
1207    D.val = val;
1208    D = spm_eeg_inv_results(D); 
1209    spm_eeg_inv_results_display(D); 
1210 
1211    if write3Dflag
1212     D.inv{val}.contrast.smooth = source_smooth_spm;
1213     D.inv{val}.contrast.rms = 1;
1214 %    D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
1215     D = spm_eeg_inv_Mesh2Voxels(D);
1216     GrpSourceImgs{Ntyp+1}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1217    end
1218    save(fullfile(D.path, D.fname), 'D');
1219 
1220    typ{Ntyp+1} = 'fused';
1221    grpfinalnam{Ntyp+1,sub} = fullfile(D.path,D.fname);
1222   end
1223  end
1224 end
1225 
1226 %%%%%%%%%%%%%%%%% Quick look at mean over subjects...
1227 
1228 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1229 S.fig = 6;
1230 for val=1:Nval
1231  S.val = val;
1232  for typn=1:(Ntyp+fusionflag)
1233   S.fig = S.fig+1;
1234   S.P = strvcat(grpfinalnam{typn,:});
1235   meg_inv_display_con(S);
1236   drawnow
1237  end
1238 end
1239 
1240 % The group fused results look a bit more focal than before, even though
1241 % gradiometer result look a bit less likely. Note that this is just the
1242 % mean over subjects, so a better test is statistics over subjects below...
1243 
1244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1245 %%%%%%%%%%%%%%%%% 3D  SPMs 
1246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1247 
1248 stdir = fullfile(wd,'GrpSourceSPMs');
1249 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1250 
1251 clear grpsubs imgfiles grpcons
1252 grpsubs{1} = [1:Nsub];
1253 % If want to split subjects, eg into two groups, eg:
1254 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1255 grpcons = [1:2];
1256 Ngrp = length(grpsubs);
1257 for typn = 1:(Ntyp+fusionflag)
1258  outdir = fullfile(stdir,typ{typn})
1259  for grp = 1:Ngrp
1260   for sub = 1:length(grpsubs{grp});
1261    subnum = grpsubs{grp}(sub);
1262    imgfiles{grp}{sub} = GrpSourceImgs{typn}{subnum}(grpcons,:);
1263   end
1264  end
1265  meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1266 end
1267 
1268 % Again you need to know a bit about the SPM Results interface, but press
1269 % Results and select an SPM.mat file for one of the modalities and enter
1270 % a new T-contrast [1 -1] to find voxels in which greater source
1271 % amplitude for faces than scrambled faces, choose 0.001 uncorrected (not
1272 % much survives whole-image correction because only 8 subjects, so few
1273 % df's ie low power and RFT conservative) and then press "Volume" to get
1274 % table of stats. It will be very rough (speckled) because of the low
1275 % dfs. The results for mags, for example, seem better (eg more
1276 % suprathreshold voxels) for this group inversion than previously (above),
1277 % though it is less clear that the fused results are better following
1278 % group inversion of each modality alone (this is an area of current
1279 % investigation). You can decide not to use group inversion (or fusion)
1280 % of course. Again, all these statistical maps look better with more subjects, 
1281 % or when using SNPM or PPM's (see Henson et al, 2007, Neuroimage and
1282 % Taylor & Henson, in prep, for more information).

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2009-07-24 16:05:00, 69.8 KB) [[attachment:cbu_meeg_spm5_pipeline.m]]
  • [get | view] (2009-04-16 12:18:36, 64.2 KB) [[attachment:cbu_meeg_spm5_pipeline_13_4_09.m]]
  • [get | view] (2009-05-01 15:22:55, 64.4 KB) [[attachment:cbu_meeg_spm5_pipeline_24-7-09.m]]
  • [get | view] (2009-04-11 11:14:19, 51.7 KB) [[attachment:cbu_meeg_spm5_pipeline_5-3-09.m]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.