Attachment 'cbu_meeg_spm5_pipeline_24-7-09.m'
Download 1 P% 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 sennam{1} = 'mags'; sennam{2} = 'grds'; sennam{3} = 'eeg';
21 Nsen = length(sennam);
22 eegflag = [0 0 1]; % Binary flag for which sensor is EEG
23
24 % Num channels per sensor 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:Nsen;
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 fuseinvertflag = 1; % Whether want to simultaneously invert (fuse) modalities (sennam)
129 % (Henson et al, submitted)
130 grpinvertflag = 1; % Whether invert using group priors (after individual inversions)
131 % (Litvak & Friston, 2008, Neuroimage)
132 fmriinvertflag = 1; % Whether want to use fMRI priors on inversion
133 % (Flandin et al, submitted)
134
135 %%%%%%%%%% General parameters
136 % Maxfilter
137 downsample_maxf = '-ds 4'; % maxfilter downsampling factor (eg 4 to save space)
138 % though beware of effect of downsampling on triggers
139 % (http://imaging.mrc-cbu.cam.ac.uk/meg/CbuSpmParameters)
140 %downsample_maxf = ''; % (if none)
141 format_maxf = '-format float'; % maxfilter data format (native = float32)
142 %format_maxf = '-format short'; % (if to save space)
143 trans_offset_maxf = [0 -13 +6]; % New centre for trans default
144 % (see our maxfilter WIKI)
145
146 % SPM
147 downsample_spm = 100; % spm new sample rate (Hz; should be at
148 % least twice filt_cut_spm below)
149 filt_cut_spm = 40; % Lowpass filter cutoff for SPM (Hz)
150 dformat_spm = 'int16'; % data format for spm (eg to save space)
151 sensor_smooth_spm = [5 5 10]; % Smoothing in SensorSpace-Time [mm mm ms];
152 source_smooth_spm = [10 10 10]; % Smoothing in Source space [mm mm mm]
153
154 %%%%%%%%%%% Initialisation
155 subsphere = [];
156 allsssbad = cell(Nsub,1);
157 nbadchan = cell(Nsen,Nsub);
158 nrejects = cell(Nsen,Nsub);
159 nevents = cell(Nsen,Nsub);
160
161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
162 %%%%%%%%%%% Preprocess MEG/EEG data
163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164
165 for sub = 1:Nsub
166
167 subnum = dosubs(sub);
168 Nses(sub) = size(fifdata{subnum},1);
169
170 cd(wd)
171 swd = sprintf('sub%d',subnum);
172 try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end
173
174 for ses = 1:Nses(sub)
175
176 allsssbad{sub} = cell(1,Nses(sub));
177
178 fnm_stem = sprintf('sub%d_ses%d',subnum,ses);
179
180 if maxfiltflag
181
182 %%%%% Use autobad to detect bad channels from first 20s, and select from output
183 rawfile = fifdata{subnum}{ses};
184 sssfile = sprintf('%s_sss.fif',fnm_stem)
185 logfile = sprintf('%s_autobad.log',fnm_stem)
186 headptsfile = sprintf('%s_headpoints.txt',fnm_stem)
187
188 if exist(sssfile), delete(sssfile), end
189 if exist(logfile), delete(logfile), end
190 [Centre,Radius] = meg_fit_sphere(rawfile,pwd,headptsfile);
191 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))
192 subsphere(sub,:) = [Centre Radius];
193
194 disp('Checking first 20s for bad channels...')
195 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))
196
197 badfile = sprintf('%s_badchan.txt',fnm_stem)
198 if exist(badfile), delete(badfile), end
199 eval(sprintf('!cat %s | sed -n ''/Static/p'' | cut -f 5- -d '' '' > %s',logfile,badfile));
200 sssbadchans{sub}{ses} = unique([textread(badfile,'%d')' badchans{1,subnum,ses}]);
201 for bc=1:length(sssbadchans{sub}{ses}); allsssbad{sub}{ses} = [allsssbad{sub}{ses} sprintf(' %04d',sssbadchans{sub}{ses}(bc))]; end
202
203 %%%%% Now mark bad channels and run ST
204 logfile = sprintf('%s_ssst.log',fnm_stem)
205 if exist(sssfile), delete(sssfile), end % Don't need previous file from autobad
206 sssfile = sprintf('%s_ssst.fif',fnm_stem)
207 if exist(logfile), delete(logfile), end
208 posfile = sprintf('%s_mvpos.txt',fnm_stem)
209
210 disp(['Now maxfiltering with bad channels (and tSSS): ' allsssbad{sub}{ses}])
211 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))
212
213 S.logfile = logfile;
214 [max_d(sub),max_r(sub)] = meg_plot_maxfilter_move(S);
215 % Do check the movement parameters - particularly goodness of fit -
216 % because movement compensation below will be invalidated if poor fit
217 % (also, the 'inter' MaxFilter option is used below in case the HPI
218 % temporarily fails (which it does for sub5), but if your subject's
219 % HPI failed permanently from the start, you might not know this otherwise,
220 % because 'inter' would handle it fine!).
221
222 infile = sssfile;
223
224 %%%%% Now mark bad channels and run MVCOMP
225 if maxmvcompflag
226
227 sss2file = sprintf('%s_mvcomp.fif',fnm_stem)
228 logfile = sprintf('%s_mvcomp.log',fnm_stem)
229 if exist(sss2file), delete(sss2file), end
230 if exist(logfile), delete(logfile), end
231
232 disp(['Now compensating movement... (every 1000ms)'])
233 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))
234 % 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))
235
236 infile = sss2file;
237 end
238
239 %%%%% Now Trans Default (if requested)
240 if maxtransflag
241 trans_origin = Centre(1:3) + trans_offset_maxf;
242 logfile = sprintf('%s_trans.log',fnm_stem)
243 sss3file = sprintf('%s_trans.fif',fnm_stem)
244 if exist(sss3file), delete(sss3file), end
245 if exist(logfile), delete(logfile), end
246
247 disp(['Now trans-ing to DEFAULT ' num2str(round(trans_origin))])
248 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))
249 end
250 end
251
252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
253 %%%%% Olaf's channel checking (have not used for a while; should put in SVN!)
254 if chancheckflag
255 addpath /imaging/olaf/MEG/artefact_scan_tool/
256 global fiffs tasks values criteria sensors datapara figs;
257 ini_check4artefacts; % Initialisation of variables
258 epochlength = 5;
259 amplitude_threshold_mag = 5000;
260 amplitude_threshold_gra = 10000;
261 out_dir = fullfile(wd,'ChanCheck');
262 addfiff(infile,fnm_stem);
263 addtask('maxminamp'); % maximum minus minimum amplitude within epoch
264 addtask('maxmingrad'); % maximum signal gradient (amplitude difference between successive data points in time) with epoch
265 addtask('threshold'); % number of epochs with max-min amplitudes above specified threshold
266 addtask('crosscorr'); % intercorrelation matrix for magnetometers and gradiometers separately
267 addtask('trendamp'); % linear trend in max-min differences across all epochs
268 addtask('trendgrad'); % linear trend in maximum signal gradients across all epochs
269 check4artefacts; % That's where it all happens... lean back and enjoy!
270 end
271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
272
273 %%%%% Read into SPM (writes separate files for MEG and EEG)
274
275 clear S
276 if usemaxtransflag
277 S.Fdata = sprintf('%s_trans.fif',fnm_stem);
278 if maxmvcompflag
279 S.HPIfile = sprintf('%s_mvcomp.fif',fnm_stem); % If trans'ed (hpi lost)
280 else
281 S.HPIfile = sprintf('%s_sss.fif',fnm_stem); % If trans'ed (hpi lost)
282 end
283 elseif usemaxmvcompflag
284 S.Fdata = sprintf('%s_mvcomp.fif',fnm_stem);
285 S.HPIfile = [];
286 else
287 S.Fdata = sprintf('%s_sss.fif',fnm_stem);
288 S.HPIfile = [];
289 end
290 % meg_plot_fifpoints(S); % If want to check points
291 S.Fchannels = 'FIF306_setup.mat';
292 S.veogchan = [62 63]; S.veog_unipolar = 0;
293 S.heogchan = [61 64]; S.heog_unipolar = 0;
294 S.trig_chan = 'STI101';
295 S.twin = [0 Inf];
296 S.eeg_ref = 'average';
297 S.trig_type = 'positive_from_zero';
298 % S.Dtype = 'float32'; % make int16 if files too big
299 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
300 % S.Fchannels_eeg = fullfile(wd,swd,sprintf('eeg_montage_sub%02d.mat',subnum)); % use subject-specific (care labels may not be standard)
301
302 if respmflag
303 D = spm_eeg_rdata_FIF(S);
304 basefile = D.fname(1:(end-4));
305 else
306 basefile = S.Fdata(1:(end-4));
307 end
308
309 %%%%% Downsample, then split Mags and Grads (file too big otherwise)
310
311 if reprocessflag
312 clear S; S.D = sprintf('%s-eeg.mat',basefile); disp(S.D)
313 S.Radc_new = downsample_spm;
314 D = spm_eeg_downsample(S);
315
316 clear S; S.D = basefile;
317 S.Radc_new = downsample_spm;
318 D = spm_eeg_downsample(S);
319 basefile = D.fname(1:(end-4));
320
321 clear S; S.D = D.fname;
322 D = spm_eeg_splitFIF(S);
323 end
324
325 %%%%% If want to view rawdata (with Danny's tool):
326 % meg_viewdata(D.fname);
327
328 %%%%% Preprocess each sensor type
329
330 for sen = 1:Nsen
331
332 if reprocessflag
333 clear S
334 S.D = sprintf('%s-%s.mat',basefile,sennam{sen}); disp(S.D)
335 S.filter.type = 'butterworth';
336 S.filter.band = 'low'; S.filter.PHz = filt_cut_spm;
337 % S.filter.band = 'stop'; S.filter.PHz = [49 51];
338 S.Dtype = dformat_spm;
339 D = spm_eeg_filter(S);
340
341 %%%%%% !!!Insert ICA here if want to remove EOG/ECG artifacts!!!
342 %%%%%% (See Jason Taylor; assume filtered, downsampled data appropriate)
343
344 clear S; S.D = D.fname;
345 S.events.start = expwin(1);
346 S.events.stop = expwin(2);
347 S.events.types = expcodes;
348 S.events.Inewlist = 0;
349 S.events.resynch = offset;
350 D = spm_eeg_epochs(S);
351
352 % re-classify event-codes
353 for e = 1:length(expcodes);
354 fi = find(D.events.code==expcodes(e));
355 if ~isempty(fi)
356 D.events.code(fi) = newcodes(e);
357 end
358 end
359 D.events.types = unique(D.events.code);
360 D.events.Ntypes = length(D.events.types);
361 save(D.fname,'D')
362 clear S; S.D = D.fname;
363
364 else
365 S.D = sprintf('e_fd%s-%s.mat',basefile,sennam{sen}); disp(S.D)
366 D = spm_eeg_ldata(S.D);
367 end
368
369 %%%%% Threshold channels (eg EOG)
370
371 nc = D.Nchannels;
372 S.thresholds.External_list = 0;
373 S.thresholds.Check_Threshold = 1;
374 S.thresholds.threshold = chan_thr{sub,sen};
375 S.BadThr = 0.2;
376 S.artefact.weighted = 0;
377 if eegflag(sen)==0
378 % S.channels.Bad = unique([D.channels.Bad strcmpi(D.channels.name....sssbadchans{sub}{ses}]); % If dont trust MaxFilter channel reconstruction
379 else
380 S.channels.Bad = unique([D.channels.Bad badchans{2,subnum,ses}]);
381 end
382 D = spm_eeg_artefact(S);
383
384 nbadchan{sen,sub} = [nbadchan{sen,sub} length(D.channels.Bad)];
385 nrejects{sen,sub} = [nrejects{sen,sub} length(find(D.events.reject))];
386
387 %%%%% Re-reference EEG if removed channels
388 if eegflag(sen) & ~isempty(D.channels.Bad)
389 D = spm_eeg_ldata(D.fname);
390 meg_reavg_ref(D);
391 save(D.fname,'D')
392 end
393 sesnam{sen,ses} = D.fname;
394
395 %%%%% Uncomment if want to average and contrast each session separately
396 %%%%% (as well as after merging below)
397 %
398 % clear S; S.D = D.fname;
399 % D = spm_eeg_average(S);
400 %
401 % clear S; S.D = D.fname;
402 % S.c = expcons;
403 % S.WeightAve = 1;
404 % D = spm_eeg_weight_epochs(S);
405 % nevents{sen,sub} = D.events.repl;
406
407 end % end of sen loop
408 end % end of ses loop
409
410 %%%%% Merge sessions, average, contrast
411
412 for sen = 1:Nsen
413
414 if Nses(sub)>1
415 if ~usemaxtransflag
416 disp(['Concatenating sessions... (NB: merged files will not have correct channel Loc/Orient info unless trans-ed to first or to default!!)']);
417 end
418 clear S; S.D = strvcat(sesnam{sen,:});
419 S.recode = cell(1,Nses(sub));
420 D = spm_eeg_merge(S);
421 clear S; S.D = D.fname;
422 else
423 clear S; S.D = sesnam{sen,1};
424 D = spm_eeg_ldata(S.D);
425 end
426 preavgnam{sen,sub} = fullfile(D.path,S.D);
427
428 Nsubchan(sub,sen) = length(D.channels.eeg) - length(D.channels.Bad);
429
430 D = spm_eeg_average(S);
431
432 clear S; S.D = D.fname;
433 S.c = expcons;
434 S.WeightAve = 0;
435 D = spm_eeg_weight_epochs(S);
436
437 nevents{sen,sub} = D.events.repl;
438 avgnam{sen,sub} = fullfile(D.path,D.fname);
439 end
440 end
441
442 cd(wd)
443
444 for sen = 1:Nsen
445 disp(sprintf('%s...',sennam{sen}))
446 for sub = 1:Nsub
447 subnum = dosubs(sub);
448 disp(sprintf('\tSub %d (%d)...',subnum,sub))
449 disp(sprintf('\t\tNevents (per condition): %s',mat2str(nevents{sen,sub})))
450 disp(sprintf('\t\tRejects (per session): %s',mat2str(nrejects{sen,sub})))
451 disp(sprintf('\t\tBadchans (per session): %s',mat2str(nbadchan{sen,sub})))
452 end
453 end
454
455 %%%% Create RMS of Grads (mainly for Sensor-Time image files)
456
457 for sen = 1:Nsen
458 if strcmp(sennam{sen},'grds')
459 for sub = 1:Nsub
460 cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
461 clear S; S.D = avgnam{sen,sub};
462 S.do_write = 1;
463 D = meg_grds2grms(S); % RMS of grads; just for sensor-level analyses
464 avgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
465 nevents{Nsen+1,sub} = nevents{sen,sub};
466 end
467 Nsen=Nsen+1; sennam{Nsen} = 'grms';
468 write2Dflag(Nsen) = 1;
469 end
470 end
471
472 %%%% Write out Sensor-Time image files for each subject
473
474 for sen = find(write2Dflag)
475 for sub = 1:Nsub
476 cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
477 disp(avgnam{sen,sub})
478 [pth,nam,ext] = fileparts(avgnam{sen,sub});
479 clear S; S.Fname = [nam ext];
480 S.interpolate_bad = 0;
481 S.n = 32;
482 S.pixsize = 3;
483 S.trialtypes = [1:Ncons];
484 P = spm_eeg_convertmat2ana3D(S);
485 if any(sensor_smooth_spm)
486 if length(sensor_smooth_spm)==1; sensor_smooth_spm = ones(1,3)*sensor_smooth_spm; end
487 for con = 1:size(P,1);
488 [pth,nam,ext] = fileparts(P(con,:));
489 Pout = fullfile(pth,['s' nam ext]);
490 spm_smooth(spm_vol(P(con,:)),Pout,sensor_smooth_spm);
491 Pin = strvcat(P(con,:),Pout);
492 spm_imcalc_ui(Pin,Pout,'((i1+eps).*i2)./(i1+eps)',{[],[],'float32',0}); %Reinsert NaNs
493 SensorTimeImgs{sen}{sub}(con,:) = fullfile(wd,swd,Pout);
494 end
495 else
496 for con = 1:size(P,1);
497 SensorTimeImgs{sen}{sub}(con,:) = fullfile(wd,swd,P(con,:));
498 end
499 end
500 end
501 end
502
503 %%%% Create Grand Mean across subjects (inc RMS (grms))
504 cd(wd)
505 for sen = 1:Nsen
506 clear S; S.P = strvcat(avgnam{sen,:});
507 S.Pout = sprintf('G%d-%s.mat',Nsub,sennam{sen})
508 D = spm_eeg_grandmean(S);
509 end
510
511 eval(sprintf('save preproc%d.mat subsphere nrejects nbadchan nevents allsssbad dosubs preavgnam avgnam',length(dosubs)))
512
513 disp('!!!NOW USE SPM-Display-M/EEG button to examine results, eg G8-*.mat files !!!')
514
515 %return % uncommment if want script to finish here
516
517
518 % Once displaying one of the modalities, you can hold shift and select both
519 % conditions 1+2, and click on a posterior right channel to see the M170
520 % difference between faces (1) and scrambled faces (2). You can also select
521 % condition 3, which is the differential ERF/ERP between faces and scrambled,
522 % then press "topography", enter -100ms for initial time, then press and hold
523 % the time slider to see topography of difference at every time point (displayed
524 % in Matlab window) - random noise until ~100ms, when strong differences emerge
525
526
527 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
528 %%%%%%%%%%%%%%%%% 2D sensor x 1D Time SPMs
529 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
530
531 stdir = fullfile(wd,'SensorTimeSPMs');
532 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
533
534 clear grpsubs grpcons imgfiles
535 grpsubs{1} = [1:Nsub];
536 % If want to split subjects, eg into two groups, eg:
537 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
538 grpcons{1} = [1:2]; % Include each condition for mags
539 grpcons{3} = [1:2]; % Include each condition for eeg
540 grpcons{4} = [1:2]; % Diff of RMS (be careful if diff number of trials for diff conditions)
541 %grpcons{4} = [3]; % Or RMS of diff (not equivalent to above)
542
543 Ngrp = length(grpsubs);
544 for sen = find(write2Dflag)
545 outdir = fullfile(stdir,sennam{sen})
546 for grp = 1:Ngrp
547 for sub = 1:length(grpsubs{grp});
548 subnum = grpsubs{grp}(sub);
549 imgfiles{grp}{sub} = SensorTimeImgs{sen}{subnum}(grpcons{sen},:);
550 end
551 end
552 meg_batch_anova(imgfiles,outdir);
553 end
554
555 disp('!!!NOW press Results button in SPM!!!')
556
557 %return
558
559 % You need to know a bit about the SPM Results interface, but press the Results
560 % button (make sure SPM is in "EEG" mode) (or better still, use the "ns" toolbox
561 % under "Toolbox" button - see WIKI page below) and select the SPM.mat file for
562 % one of the modalities, select the default effects of interest F-contrast
563 % (because we don't care about polarity, at least for Mags and EEG), choose
564 % 0.001 uncorrected. Then press "Volume" to get table of stats (ignore brain image).
565 % You should see, at least in Mags and EEG, frontal and posterior clusters
566 % (simply of different polarity) maximal around 170ms. Not much survives whole-
567 % image correction because only 8 subjects, so few df's ie low power and RFT
568 % conservative) - though effects around 170ms do for corrected extent-level
569 % for Mags and EEG (using "ns" toolbox); GRMS results are poor though.
570 % Correction for height would also be significant if you used SVC because
571 % of an a priori timewindow of interest around 170ms (again see WIKI page below).
572 %
573 % For more info on interpreting these SPMs, see
574 % http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm
575 %
576 % Also, for an example use with EEG data, see Henson et al (2008), Neuroimage.
577 %
578 % Note that you can do SVC for a priori timewindows of interest (see end of above WIKI)
579 % But one use of this SPMs is to define timewindows of interest (in the absence of
580 % a priori ones) for the subsequent source localisation below...
581 %
582 % (Note that these SPMs can also handle different numbers and locations of channels
583 % across subjects, provided their locations are digitised in same space, because the
584 % data are interpolated to same grid. For MEG data, particularly gradiometers, the
585 % results benefit from using 'trans -default' above, to realign across different
586 % headpositions for different subjects - see Smith et al (in press), HBM Abstract.
587
588
589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
590 %%%%%%%%%%%%%%%%% Get MRIS
591 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
592
593 for sub=1:Nsub
594 subnum = dosubs(sub);
595 swd = sprintf('sub%d',subnum);
596 mwd{subnum} = fullfile(wd,swd,'MRI');
597 try cd(mwd{subnum}); catch eval(sprintf('!mkdir %s',mwd{subnum})); cd(mwd{subnum}); end
598
599 if getmriflag
600 hdr = spm_select('FPList',mridata{subnum},'^*.dcm');
601 hdr = spm_dicom_headers(hdr);
602
603 patsex{subnum} = hdr{1}.PatientsSex;
604 patage{subnum} = hdr{1}.PatientsAge; % NB At time of scanning, not time of MEG!
605
606 spm_dicom_convert(hdr,'all','flat','nii');
607
608 mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-0[0-9].nii$');
609 % (filter setting just to ensure picks up raw structural - assuming CBU onvention)
610 mrifile{sub} = strtrim(mrifile{sub}(1,:)); % Take first if multiple previous
611
612 disp('!!!GOOD IDEA TO MANUALLY REPOSITION STRUCTURAL VIA DISPLAY BUTTON!!!')
613 % pause
614
615 if VitECapsules(subnum)
616 meg_headmask(mrifile{sub},'display'); % Remove any VitE capsules
617 disp('!!!Check difference image!!!');
618 % pause
619 eval(sprintf('!mv %s_masked.nii %s',mrifile{sub}(1:end-4),mrifile{sub}))
620 end
621 patage
622 patsex
623
624 else
625 mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-0[0-9].nii$');
626 mrifile{sub} = strtrim(mrifile{sub}(1,:)); % Take first if multiple previous
627 wmrifile{sub} = spm_select('FPList',mwd{subnum},'^wms.*-0[0-9].nii$');
628 wmrifile{sub} = strtrim(wmrifile{sub}(1,:)); % Take first if multiple previous
629 end
630 end
631 spm_check_registration(strvcat(mrifile)); % Only feasible for <~20 images
632 try, spm_check_registration(strvcat(wmrifile)); end % (sub5 bit too warped at back)
633
634 disp('!!!NOW DISPLAY MRI IMAGE in SPM AND DETERMINE FIDUCIALS!!!')
635
636 %%%%%%%%%% Coordinates of fiducials from above MRIs (nasion, left, right)
637 % Below are from using MRI as is (without manually reorienting) - very approximate
638 smri_fids{1} = [ 3 111 22; -78 18 -31; 77 11 -33];
639 smri_fids{2} = [ 6 107 -28; -69 26 -41; 68 14 -45];
640 smri_fids{3} = [-2 124 -23; -78 22 -52; 70 21 -58];
641 smri_fids{4} = [ 3 109 -8; -72 18 -41; 72 18 -47];
642 smri_fids{5} = [-15 103 -33; -79 3 -38; 62 13 -46];
643 smri_fids{6} = [-4 109 -20; -71 15 -32; 63 11 -35];
644 smri_fids{7} = [-1 93 -5; -72 14 -44; 70 10 -44];
645 smri_fids{8} = [ 1 99 -22; -73 8 -27; 69 10 -42];
646 smri_fids{9} = [-3 125 -24; -76 22 -55; 69 22 -58];
647 % Below are after my manual reorienting (just for your reference)
648 %smri_fids{1} = [ 0.7 80.0 -31.2; -76.9 -14.3 -45.9; 79.9 -19.2 -41.0];
649 %smri_fids{2} = [ 1.3 78.7 -30.9; -71.8 -2.5 -58.0; 71.5 1.8 -58.0];
650 %smri_fids{3} = [ 1.7 88.8 -24.6; -74.0 -1.4 -66.1; 71.7 -2.8 -66.1];
651 %smri_fids{4} = [-1.6 78.9 -22.4; -74.8 -13.3 -57.7; 71.6 -7.2 -54.6];
652 %smri_fids{5} = {TBA};
653 %smri_fids{6} = [-2.4 79.1 -27.1; -72.1 -9.8 -44.6; 64.3 -9.8 -51.9];
654 %smri_fids{7} = [-2.3 73.2 -7.0; -71.8 -2.0 -46.9; 70.2 -4.9 -49.8];
655 %smri_fids{8} = [ 0.0 78.1 -30.3; -71.3 -1.5 -48.5; 72.8 6.0 -51.5];
656 %smri_fids{9} = [ 0.0 88.9 -31.1; -76.8 -15.5 -56.6; 78.0 -12.7 -55.0];
657
658
659 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
660 %%%%%%%%%%%%%%%%% Localise on Mesh
661 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662
663 addpath /imaging/local/Brainstorm/Toolbox
664
665 clear sennam;
666 sennam{1} = 'mags'; sennam{2} = 'grds'; sennam{3} = 'eeg';
667 Nsen = length(sennam);
668
669 % !!! There are several options for source localisation, eg choice of
670 % forward model (eg Spheres or BEMs), and choice of inversion method (eg
671 % MSP, COH, MNM). Below illustrates just two possible inversions (each
672 % inversion is indexed by "val"), the second of which has been commented
673 % out because BEMs take a VERY LONG time (and don't work well for EEG yet).
674 % Of course many more options can be tried (and compared using the model-evidence,
675 % recorded in mod_evd below, provided that you don't also change the
676 % data, eg epoch inverted).
677 %
678 % For more information about forward models, see
679 % http://imaging.mrc-cbu.cam.ac.uk/meg/SpmForwardModels
680 % For more information about inversion parameters, see
681 % Friston et al (2008) Neuroimage
682 % (PDF available here http://imaging.mrc-cbu.cam.ac.uk/meg/SpmAnalysis)
683
684 %%%%%%%%%%%%%%%%%%%%
685 val = 1; % Spherical forward model using canonical cortical mesh but individual head meshes
686 comment{val} = 'Concentric Spheres';
687 cortex{val} = 'Can'; % Canonical (inverse-normalised template) cortical mesh
688 CtxSize(val) = 8196; % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
689 iskull{val} = 'Ind'; % Individual (SPM-created) inner skull mesh
690 oskull{val} = ''; % (Outer skull irrelevant for spherical models)
691 scalp{val} = 'Ind'; % Individual (SPM-created) scalp mesh
692 formod{val} = 'Sph'; % Concentric spheres
693 surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all sensor-types
694
695 invcons{val} = [1 2]; % Conditions to invert
696 invtype{val} = 'GS'; % Type of inversion (This is MSP, with a Greedy Search (GS))
697 invtwin{val} = expwin; % Invert whole epoch
698
699 %%%%% Options for subsequent time-freq contrast
700 contwin{val}{1} = [140 190]; % Time window for contrast
701 contwin{val}{2} = []; % Frequency window for contrast ([]=all freqs)
702 conengy{val} = 'evoked'; % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
703
704 %%%%%%%%%%%%%%%%%%%%
705 val = 2; % As val=1 except IID (ie standard minimum norm) used to invert
706 comment{val} = 'Concentric Spheres';
707 cortex{val} = 'Can'; % Canonical (inverse-normalised template) cortical mesh
708 CtxSize(val) = 8196; % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
709 iskull{val} = 'Ind'; % Individual (SPM-created) inner skull mesh
710 oskull{val} = ''; % (Outer skull irrelevant for spherical models)
711 scalp{val} = 'Ind'; % Individual (SPM-created) scalp mesh
712 formod{val} = 'Sph'; % Concentric spheres
713 surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all data-types
714
715 invcons{val} = [1 2]; % Conditions to invert
716 invtype{val} = 'IID'; % Type of inversion
717 invtwin{val} = expwin; % Invert whole epoch
718
719 %%%%% Options for subsequent time-freq contrast
720 contwin{val}{1} = [140 190]; % Time window for contrast
721 contwin{val}{2} = []; % Frequency window for contrast ([]=all freqs)
722 conengy{val} = 'evoked'; % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
723
724 %%%%%%%%%%%%%%%%%%%%
725 %val = 2; % BEM using canonical meshes (Warning: BEMs take a long time to compute!)
726 %
727 %comment{val} = 'Boundary Element Model';
728 %cortex{val} = 'Can'; % Canonical (inverse-normalised template) cortical mesh
729 %CtxSize(val) = 8196; % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
730 %iskull{val} = 'Can'; % Canonical (inverse-normalised template) inner skull mesh
731 %oskull{val} = 'Can'; % Canonical (inverse-normalised template) outer skull mesh
732 %scalp{val} = 'Can'; % Canonical (inverse-normalised template) scalp mesh
733 %formod{val} = 'BEM'; % Boundary Element Model
734 %surfit(val,:)= [2 2 4]; % Surfaces for BEM: up to 2nd for MEG (iskull) and 4th for EEG (+oskull+scalp)
735 %invcons{val} = [1 2]; % Conditions to invert
736 %invtype{val} = 'GS'; % Type of inversion (This is MSP, with a Greedy Search (GS))
737 %invtwin{val} = expwin; % Invert whole epoch
738
739 %%%%% Options for subsequent time-freq contrast
740 %contwin{val}{1} = [140 190]; % Time window for contrast
741 %contwin{val}{2} = []; % Frequency window for contrast ([]=all freqs)
742 %conengy{val} = 'evoked'; % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
743
744 % Other options you can consider are "overlapping spheres" for the
745 % forward model (see spm_eeg_inv_BSTfwdsol.m), minimumum norm (IID) for
746 % the inversion, etc. See Henson et al (2009), Neuroimage, for more options.
747
748 Nval = length(comment);
749 mod_evd=[];
750
751 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
752
753 % If want to pass prior estimates of sensor-level error (noise), eg
754 % estimates from empty-room data in case of MEG, such as below...
755 %tmp = load('/imaging/rh01/Methods/MEG_Group/MeanCn.mat')
756 %Cn{1} = tmp.Cn{1}-eye(size(tmp.Cn{1}))*mean(diag(tmp.Cn{1})); % Orthogonalise
757 %Cn{2} = tmp.Cn{2}-eye(size(tmp.Cn{2}))*mean(diag(tmp.Cn{2})); % Orthogonalise
758 % !!!Care: Cn only relevant to SSSt and 40Hz lowpass (see Henson et al, submitted)!!!
759
760 for sub = 1:Nsub
761
762 cd(wd)
763 subnum = dosubs(sub);
764 swd = sprintf('sub%d',subnum);
765 cd(swd)
766
767 % If want to pass prior estimates of sensor-level error (noise)...
768 % !!!Care - sensor order in sennam must match order in Cn, ie Cn{1}=mags, Cn{2}=grads!!!
769 % Ce{1} = {speye(Nsubchan(sub,1)) Cn{1}}; % Mags (Cn assumes no bad channels!!!!)
770 % Ce{2} = {speye(Nsubchan(sub,2)) Cn{2}}; % Grds (Cn assumes no bad channels!!!!)
771 % Ce{3} = {speye(Nsubchan(sub,3))}; % EEG (no noise estimates yet)
772 % Otherwise, do not pass anything, corresponding to just IID ('white') noise only:
773 Ce = cell(1,Nsen);
774
775 for val = 1:Nval
776
777 for sen = 1:Nsen
778
779 clear S;
780 S.D = char(preavgnam{sen,sub});
781 % This is using all trials to estimate sources, ie total (induced+evoked) power.
782 % This gives more sample points in total, and the pure evoked part can always be
783 % extracted when estimating time-freq contrasts afterwards (see below). If you want
784 % to localise purely evoked power, then use instead:
785 % S.D = char(avgnam{sen,sub});
786 % Also may prefer "mvcomp" rather than "trans" data, if don't trust "trans default"
787
788 D = spm_eeg_ldata(S.D); disp(S.D)
789
790 %If want to clear previous localisations (safer?)
791 if val==1
792 if isfield(D,'inv'); D = rmfield(D,'inv'); save(D.fname,'D'); end
793 if isfield(D,'con'); D = rmfield(D,'con'); save(D.fname,'D'); end
794 end
795
796 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise
797
798 disp(comment{val})
799 D.val = val; D.fname
800 D.inv{val}.method = 'imaging';
801 load('defaults_eeg_inv.mat')
802 D.inv{val} = invdefts.imag;
803 D.inv{val}.date = mat2str(fix(clock));
804 D.inv{val}.comment = comment{val};
805
806 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI
807
808 %%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI
809
810 D.inv{val}.mesh.sMRI = mrifile{sub}; % In case many from previous
811
812 [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);
813 mri_stem = nam;
814 wmrifile{sub} = fullfile(pth,['wm' nam ext]);
815 if ~exist(wmrifile{sub})
816 D = spm_eeg_inv_spatnorm(D);
817 else
818 def_name = [nam '_sn_1.mat'];
819 isndef_name = [nam '_inv_sn_1.mat'];
820 D.inv{val}.mesh.def = fullfile(pth,def_name);
821 D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
822 D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
823 D.inv{val}.mesh.wmMRI = fullfile(pth,['wm' nam ext]);
824 end
825
826 %%%%%%%%%%%%%%%%%%% Meshing 2: Creating binary images for scalp (and inner skull)
827
828 D.inv{val}.mesh.msk_scalp = spm_select('List',mwd{subnum},'^*scalp\.nii$');
829
830 if isempty(D.inv{val}.mesh.msk_scalp)
831 D = spm_eeg_inv_getmasks(D);
832 else % assume other masks exist too!
833 D.inv{val}.mesh.msk_scalp = fullfile(mwd{subnum},D.inv{val}.mesh.msk_scalp);
834 D.inv{val}.mesh.msk_iskull = spm_select('FPList',mwd{subnum},'^*iskull\.nii$');
835 D.inv{val}.mesh.msk_cortex = spm_select('FPList',mwd{subnum},'^*cortex\.nii$');
836 D.inv{val}.mesh.msk_flags = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
837 end
838 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));
839
840 %%%%%%%%%%%%%%%%%%% Meshing 2: Creating meshes
841
842 D.inv{val}.mesh.template = 0;
843 D.inv{val}.mesh.canonical = 1;
844 D.inv{val}.mesh.Msize = CtxSize(val);
845
846 % Prompt for canonical or user-specified cortical meshs (SPM does not currently
847 % produce individual ones, because automatic meshing of cortex is difficult!)
848
849 if strcmp(cortex{val},'Can')
850 Tmesh = load(fullfile(spm('dir'),'EEGtemplates',sprintf('mni_mesh_cortex_%d.mat',CtxSize(val))));
851 D.inv{val}.mesh.tess_mni.vert = Tmesh.vert;
852 D.inv{val}.mesh.tess_mni.face = uint16(Tmesh.face);
853 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
854 D.inv{val}.mesh.tess_ctx.vert = Tmesh.vert;
855 D.inv{val}.mesh.tess_ctx.face = uint16(Tmesh.face);
856 elseif strcmp(cortex{val},'Load')
857 Tmesh = spm_select(1,'mat','Select cortex mat file containing vert and face',[],pwd,'.*mat');
858 D.inv{val}.mesh.tess_ctx.vert = Tmesh.vert;
859 D.inv{val}.mesh.tess_ctx.face = uint16(Tmesh.face);
860 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.invdef);
861 D.inv{val}.mesh.tess_mni.vert = Tmesh.vert;
862 D.inv{val}.mesh.tess_mni.face = uint16(Tmesh.face);
863 else
864 error('No cortical mesh specified????')
865 end
866
867 iskull_mesh_name = fullfile('MRI',[mri_stem '_iskull_mesh.mat']);
868 scalp_mesh_name = fullfile('MRI',[mri_stem '_scalp_mesh.mat']);
869
870 % Prompt for canonical, individual (from SPM) or user-specified inner skull
871 % (Use canonical if worried about canonical cortex piercing innner skull)
872
873 if strcmp(iskull{val},'Can')
874 Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_iskull_2562.mat'));
875 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
876 D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
877 D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
878 elseif strcmp(iskull{val},'Ind')
879 if ~exist(iskull_mesh_name)
880 D = spm_eeg_inv_getmeshes(D,1);
881 vert = D.inv{val}.mesh.tess_iskull.vert;
882 face = D.inv{val}.mesh.tess_iskull.face;
883 save(iskull_mesh_name,'vert','face');
884 else
885 Tmesh = load(iskull_mesh_name);
886 D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
887 D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
888 end
889 elseif strcmp(iskull{val},'Load')
890 Tmesh = spm_select(1,'mat','Select iskull mat file containing vert and face',[],pwd,'.*mat');
891 D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
892 D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
893 else
894 warning('No inner skull mesh specified')
895 end
896
897 % Prompt for canonical or user-specified outer skull (SPM does not currently
898 % produce individual ones, because segmenting outer skull from scalp difficult)
899 % (Note that outer skull only really necessary for EEG BEMs)
900
901 if strcmp(oskull{val},'Can')
902 Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_oskull_2562.mat'));
903 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
904 D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
905 D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
906 elseif strcmp(oskull{val},'Load')
907 Tmesh = spm_select(1,'mat','Selectoskull mat file containing vert and face',[],pwd,'.*mat');
908 D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
909 D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
910 else
911 warning('No outer skull mesh specified')
912 end
913
914 % Prompt for canonical, individual (from SPM) or user-specified scalp
915 % (Use canonical if worried about canonical cortex or inner skull piercing scalp)
916
917 if strcmp(scalp{val},'Can')
918 Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_scalp_2562.mat'));
919 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
920 D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
921 D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
922 elseif strcmp(scalp{val},'Ind')
923 if ~exist(scalp_mesh_name)
924 D = spm_eeg_inv_getmeshes(D,2);
925 vert = D.inv{val}.mesh.tess_scalp.vert;
926 face = D.inv{val}.mesh.tess_scalp.face;
927 save(scalp_mesh_name,'vert','face');
928 else
929 Tmesh = load(scalp_mesh_name);
930 D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
931 D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
932 end
933 elseif strcmp(scalp{val},'Load')
934 Tmesh = spm_select(1,'mat','Select scalp mat file containing vert and face',[],pwd,'.*mat');
935 D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
936 D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
937 else
938 warning('No scalp mesh specified')
939 end
940
941 % Always individual scalp for datareg
942
943 if ~exist(scalp_mesh_name)
944 S = spm_eeg_inv_getmeshes(D,2); % lastarg==1 returns iskull (not scalp)
945 vert = S.inv{val}.mesh.tess_scalp.vert;
946 face = S.inv{val}.mesh.tess_scalp.face;
947 save(scalp_mesh_name,'vert','face');
948 D.inv{val}.mesh.tess_scalp_datareg.vert = vert;
949 D.inv{val}.mesh.tess_scalp_datareg.face = uint16(face);
950 else
951 Tmesh = load(scalp_mesh_name);
952 D.inv{val}.mesh.tess_scalp_datareg.vert = Tmesh.vert;
953 D.inv{val}.mesh.tess_scalp_datareg.face = uint16(Tmesh.face);
954 end
955
956 spm_eeg_inv_checkmeshes(D);
957 save(D.fname,'D');
958
959 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg
960
961 if isempty(D.inv{val}.datareg.fid_coreg)
962
963 D.inv{val}.datareg.sensors = D.channels.Loc';
964 D.inv{val}.datareg.fid_eeg = D.channels.fid_eeg;
965 D.inv{val}.datareg.headshape = D.channels.headshape;
966 if ~eegflag(sen)
967 D.inv{val}.datareg.megorient = D.channels.Orient';
968 end
969 D.inv{val}.datareg.scalpvert = D.inv{val}.mesh.tess_scalp_datareg.vert;
970 D.inv{val}.datareg.fid_mri = smri_fids{subnum};
971 D = spm_eeg_inv_datareg(D);
972 fid_err(sub) = sqrt(mean(mean((D.inv{val}.datareg.fid_coreg - D.inv{val}.datareg.fid_mri).^2)))
973 end
974
975 spm_eeg_inv_checkdatareg(D);
976 save(D.fname,'D');
977 clear S; S.D=D.fname;
978 % meg_pretty_datareg(S); % Prettier rendering of head and sensors
979
980 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model
981
982 rotate3d off
983 D.inv{val}.method = 'Imaging';
984 if ~eegflag(sen)
985 if strcmp(formod{val},'Sph'), meth = 'meg_sphere';
986 elseif strcmp(formod{val},'BEM'), meth = 'meg_bem';
987 elseif strcmp(formod{val},'OS'), meth = 'meg_os';
988 else, error('Unknown forward model!'); end
989 else
990 if strcmp(formod{val},'Sph'), meth = 'eeg_3sphereBerg';
991 elseif strcmp(formod{val},'BEM'), meth = 'eeg_bem';
992 elseif strcmp(formod{val},'OS'), meth = 'eeg_os';
993 else, error('Unknown forward model!'); end
994 end
995 D.inv{val}.forward.method = meth;
996 D.inv{val}.forward.surface2fit = surfit(val,sen);
997 % D.inv{val}.forward.NVertMax = 1000; % If want to save some time!
998 gain_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatrix_' upper(sennam{sen}) ...
999 '_' formod{val} '_' num2str(surfit(val,sen)) '.mat'])
1000 gxyz_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatxyz_' upper(sennam{sen}) ...
1001 '_' formod{val} '_' num2str(surfit(val,sen)) '.mat'])
1002
1003 % delete(gain_name); % if want to be safe (uncomment if want to save time)
1004 if ~exist(gain_name) %| ~exist(gxyz_name)
1005 [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
1006 eval(sprintf('!mv %s %s',D.inv{val}.forward.gainmat,gain_name))
1007 eval(sprintf('!mv %s %s',D.inv{val}.forward.gainxyz,gxyz_name))
1008 end
1009
1010 D.inv{val}.forward.gainmat = gain_name;
1011 D.inv{val}.forward.gainxyz = gxyz_name;
1012 save(D.fname,'D');
1013
1014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert each modality separately
1015
1016 D.inv{val}.inverse = [];
1017 D.inv{val}.inverse.trials = invcons{val};
1018 D.inv{val}.inverse.type = invtype{val};
1019 D.inv{val}.inverse.woi = invtwin{val};
1020 D.inv{val}.inverse.lpf = 0; % Confusing, but hpf refers to cutoff
1021 D.inv{val}.inverse.hpf = 40; % of lowpass (and lpf to cutoff of highpass)
1022 % Some extra parameters that you could play with...(but defaults seem ok!)
1023 % D.inv{val}.inverse.smooth = 0.8;
1024 % D.inv{val}.inverse.Np = 256;
1025 % D.inv{val}.inverse.Han = 0;
1026 D.inv{val}.inverse.Qe = Ce{sen};
1027
1028 D = spm_eeg_invert(D);
1029 % D = spm_eeg_invert_fuse(D); % If want to compare exactly with fused below
1030 save(D.fname,'D');
1031 mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end);
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{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1048 end
1049
1050 % If want to save on filesize
1051 % if val>2
1052 % canvert = D.inv{val-1}.mesh.tess_mni.vert;
1053 % canface = D.inv{val-1}.mesh.tess_mni.face;
1054 % D.inv{val-1}.mesh = [];
1055 % D.inv{val-1}.mesh.tess_mni.vert = canvert;
1056 % D.inv{val-1}.mesh.tess_mni.face = canface;
1057 % D.inv{val-1}.datareg = [];
1058 % clear canvert canface
1059 % end
1060
1061 disp(sprintf('Sub=%d, %s, val=%d: %s',sub,sennam{sen},val,D.inv{val}.comment))
1062 save(D.fname,'D');
1063
1064 end % end sen
1065
1066 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1067 % Fusion
1068
1069 if fuseinvertflag
1070 DD = {};
1071 for sen=1:Nsen
1072 DD{sen} = spm_eeg_ldata(preavgnam{sen,sub});
1073 DD{sen}.val = val;
1074 DD{sen}.inv{val}.inverse.Qe = Ce{sen};
1075 DD{sen}.inv{val}.inverse.QP = []; % Reset any parameters from last inversion
1076 end
1077 DD{1}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1078 DD{1}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1079
1080 DD{1}.inv{val}.inverse.OutExt = '_fused';
1081 D = spm_eeg_invert_fuse(DD);
1082 mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
1083
1084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1085
1086 D.inv{val}.contrast.woi = contwin{val}{1};
1087 D.inv{val}.contrast.fboi = contwin{val}{2};
1088 D.inv{val}.contrast.type = conengy{val};
1089 D.val = val;
1090 D = spm_eeg_inv_results(D);
1091 spm_eeg_inv_results_display(D);
1092
1093 if write3Dflag
1094 D.inv{val}.contrast.smooth = source_smooth_spm;
1095 D.inv{val}.contrast.rms = 1;
1096 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1097 D = spm_eeg_inv_Mesh2Voxels(D);
1098 SourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1099 end
1100
1101 sennam{Nsen+1} = 'fused';
1102 preavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
1103 end
1104
1105 end % end val
1106
1107 end % end subs
1108
1109 % Review model-evidences (note cannot compare fused to individual modality
1110 % inversions, since data differs)
1111 for sub = 1:Nsub
1112 subnum = dosubs(sub);
1113 disp(sprintf('\tSub %d (%d)...',subnum,sub))
1114 disp(sprintf('\t\tFiducial error (RMS mm): %3.2f',fid_err(sub)))
1115 for val=1:Nval
1116 disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(mod_evd(sub,val,:))))))
1117 end
1118 end
1119
1120 % Fid errors may be quite large (~8mm) for some subjects, simply because
1121 % I did not bother to mark them carefully on the MRI. Note that the
1122 % actual coregistration is determined primarily by the digitised headshape;
1123 % the fiducials are only really used to check the quality of registration.
1124
1125 %spm_check_registration(strvcat(wmrifile)) % double-check normaliation worked
1126
1127 disp('!!!Review single-subject localisations if you wish!!!')
1128
1129 % At this point, you can review a single-subject's localisation by
1130 % pressing the "3D source reconstruction" button, and then loading one
1131 % of the mmca* files (for one modality, or the "*_fused.mat" for the
1132 % fusion of the modalities). You can then press the "mip" button to see
1133 % the overall reconstruction, or the "display" button under "window" to
1134 % see the power in the time-freq contrast (for the selected
1135 % condition). You can also display a movie and change the start-stop
1136 % times, etc. Don't press the other buttons for inversion, fusion, group
1137 % inversion etc because they will be covered below...
1138
1139 %%%%%%%%%%%%%%%%% Or a quick look at mean over subjects...
1140
1141 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1142 S.fig = 2;
1143 for val=1:Nval
1144 S.val = val;
1145 for sen=1:(Nsen+fuseinvertflag)
1146 S.fig = S.fig+1;
1147 S.P = strvcat(preavgnam{sen,:});
1148 meg_inv_display_con(S);
1149 end
1150 end
1151 lastfig = S.fig;
1152
1153 % You should see bilateral fusiform and anterior inferior temporal
1154 % cortex, plus some more posterior occipitotemporal stuff on the right,
1155 % in most of the modalities, particularly the fused result
1156
1157 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1158 %%%%%%%%%%%%%%%%% 3D SPMs
1159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1160
1161 stdir = fullfile(wd,'SourceSPMs');
1162 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1163
1164 clear grpsubs imgfiles grpcons
1165 grpsubs{1} = [1:Nsub];
1166 % If want to split subjects, eg into two groups, eg:
1167 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1168
1169 grpcons = [1:2];
1170 Ngrp = length(grpsubs);
1171 for val=1:Nval
1172 valdir = fullfile(stdir,sprintf('Inv%d',val));
1173 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
1174 for sen = 1:(Nsen+fuseinvertflag)
1175 outdir = fullfile(valdir,sennam{sen})
1176 for grp = 1:Ngrp
1177 for sub = 1:length(grpsubs{grp});
1178 subnum = grpsubs{grp}(sub);
1179 imgfiles{grp}{sub} = SourceImgs{sen}{val}{subnum}(grpcons,:);
1180 end
1181 end
1182 meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1183 end
1184 end
1185 cd(stdir)
1186 disp('!!!Review the group 3D stats if you wish!!!')
1187
1188 %return
1189
1190 % Again you need to know a bit about the SPM Results interface, but press
1191 % Results and select an SPM.mat file for one of the modalities and enter
1192 % a new T-contrast [1 -1] to find voxels in which greater source
1193 % amplitude for faces than scrambled faces, choose 0.01 uncorrected (not
1194 % much survives whole-image correction because only 8 subjects, so few
1195 % df's ie low power and RFT conservative) and then press "Volume" to get
1196 % table of stats. It will be very rough (speckled) because of the low
1197 % dfs, but you should see a cloud around right fusiform for example for
1198 % the fused stats. These statistical maps look better with more subjects,
1199 % or when using SNPM or PPM's (see Henson et al, 2007, Neuroimage and
1200 % Taylor & Henson, in prep, for more information).
1201
1202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1203 %%%%%%%%%%%%%%%%% 3D PPMs
1204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1205
1206 % This uses SPMs above to calculate Posterior Probability Maps (PPMs),
1207 % which do not have any mutliple comparison issues (eg no need for RFT)
1208 % However, it will take a long time...
1209
1210 for val=1:Nval
1211 valdir = fullfile(stdir,sprintf('Inv%d',val));
1212 cd(valdir);
1213 for sen = 1:(Nsen+fuseinvertflag)
1214 outdir = fullfile(valdir,sennam{sen})
1215 cd(outdir)
1216 load('SPM.mat');
1217 spm_spm_Bayes(SPM);
1218 end
1219 end
1220 disp('!!!Review the group 3D PPMs if you wish!!!')
1221
1222 % Again you need to know a bit about the PPM Results interface, but press
1223 % Results and select an SPM.mat file for one of the modalities and enter
1224 % a new T-contrast [1 -1], but this time, select "Bayesian" rather than
1225 % "classical" when prompted. Then change the default effect size threshold
1226 % to 0, and probability threshold to 0.99. This will show voxels with a
1227 % >99% probability of the increase for faces relative to scrambled being
1228 % greater than zero.
1229
1230
1231 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1232 % Group inversion
1233
1234 grp_mod_evd=[];
1235 if grpinvertflag
1236 for sen=1:Nsen
1237 clear DD
1238 for sub=1:Nsub
1239 DD{sub} = spm_eeg_ldata(preavgnam{sen,sub});
1240 end
1241 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1242 for sub=1:Nsub
1243 DD{sub}.val = val;
1244 DD{sub}.inv{val}.comment = [DD{sub}.inv{val}.comment 'Group'];
1245 DD{sub}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1246 DD{sub}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1247 end
1248
1249 DD = spm_eeg_invert(DD);
1250
1251 % Save outputs
1252 for sub=1:Nsub
1253 D = DD{sub};
1254 D.fname = [D.fname(1:(end-4)) '_grp.mat'];
1255 grppreavgnam{sen,sub} = fullfile(D.path, D.fname);
1256 save(fullfile(D.path, D.fname), 'D');
1257 grp_mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end);
1258
1259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1260
1261 D.inv{val}.contrast.woi = contwin{val}{1};
1262 D.inv{val}.contrast.fboi = contwin{val}{2};
1263 D.inv{val}.contrast.type = conengy{val};
1264 D.val = val;
1265 D = spm_eeg_inv_results(D);
1266 spm_eeg_inv_results_display(D);
1267
1268 if write3Dflag
1269 D.inv{val}.contrast.smooth = source_smooth_spm;
1270 D.inv{val}.contrast.rms = 1;
1271 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1272 D = spm_eeg_inv_Mesh2Voxels(D);
1273 GrpSourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1274 end
1275 save(fullfile(D.path, D.fname), 'D');
1276 end % end val
1277 end % end sub
1278 end % end sen
1279
1280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1281 % Fuse any Group inversion
1282 if fuseinvertflag
1283 for sub=1:Nsub
1284 clear DD;
1285 for sen=1:Nsen
1286 DD{sen} = spm_eeg_ldata(grppreavgnam{sen,sub});
1287 end
1288
1289 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1290 for sen=1:Nsen
1291 DD{sen}.val = val;
1292 DD{sen}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1293 DD{sen}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1294 % DD{sen}.inv{val}.inverse.QP = []; % (Important not to reset the source priors QP from the grp inversion above)
1295 DD{sen}.inv{val}.inverse.OutExt = '_fused';
1296 end
1297
1298 D = spm_eeg_invert_fuse(DD);
1299 grp_mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
1300
1301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1302
1303 D.inv{val}.contrast.woi = contwin{val}{1};
1304 D.inv{val}.contrast.fboi = contwin{val}{2};
1305 D.inv{val}.contrast.type = conengy{val};
1306 D.val = val;
1307 D = spm_eeg_inv_results(D);
1308 spm_eeg_inv_results_display(D);
1309
1310 if write3Dflag
1311 D.inv{val}.contrast.smooth = source_smooth_spm;
1312 D.inv{val}.contrast.rms = 1;
1313 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1314 D = spm_eeg_inv_Mesh2Voxels(D);
1315 GrpSourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1316 end
1317 save(fullfile(D.path, D.fname), 'D');
1318 end % end val
1319
1320 sennam{Nsen+1} = 'fused';
1321 grppreavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
1322 end % end sub
1323 end
1324 end
1325
1326 % Review model-evidences (note cannot compare fused to individual modality
1327 % inversions, since data differs)
1328 for sub = 1:Nsub
1329 subnum = dosubs(sub);
1330 disp(sprintf('\tSub %d (%d)...',subnum,sub))
1331 for val=1:Nval
1332 disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(grp_mod_evd(sub,val,:))))))
1333 end
1334 end
1335
1336 %%%%%%%%%%%%%%%%% Quick look at mean over subjects...
1337
1338 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1339 S.fig = lastfig;
1340 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1341 S.val = val;
1342 for sen=1:(Nsen+fuseinvertflag)
1343 S.fig = S.fig+1;
1344 S.P = strvcat(grppreavgnam{sen,:});
1345 meg_inv_display_con(S);
1346 drawnow
1347 end
1348 end
1349 lastfig = S.fig;
1350
1351 % The group fused results look a bit more focal than before, even though
1352 % gradiometer result look a bit less likely. Note that this is just the
1353 % mean over subjects, so a better test is statistics over subjects below...
1354
1355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1356 %%%%%%%%%%%%%%%%% 3D SPMs
1357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1358
1359 stdir = fullfile(wd,'GrpSourceSPMs');
1360 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1361
1362 clear grpsubs imgfiles grpcons
1363 grpsubs{1} = [1:Nsub];
1364 % If want to split subjects, eg into two groups, eg:
1365 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1366 grpcons = [1:2];
1367 Ngrp = length(grpsubs);
1368 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1369 valdir = fullfile(stdir,sprintf('Inv%d',val));
1370 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
1371 for sen = 1:(Nsen+fuseinvertflag)
1372 outdir = fullfile(valdir,sennam{sen})
1373 for grp = 1:Ngrp
1374 for sub = 1:length(grpsubs{grp});
1375 subnum = grpsubs{grp}(sub);
1376 imgfiles{grp}{sub} = GrpSourceImgs{sen}{val}{subnum}(grpcons,:);
1377 end
1378 end
1379 meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1380 end
1381 end
1382 cd(stdir)
1383 disp('!!!Review the group 3D stats if you wish!!!')
1384
1385 %return
1386
1387 % Again you need to know a bit about the SPM Results interface, but press
1388 % Results and select an SPM.mat file for one of the modalities and enter
1389 % a new T-contrast [1 -1] to find voxels in which greater source
1390 % amplitude for faces than scrambled faces, choose 0.001 uncorrected (not
1391 % much survives whole-image correction because only 8 subjects, so few
1392 % df's ie low power and RFT conservative) and then press "Volume" to get
1393 % table of stats. It will be very rough (speckled) because of the low
1394 % dfs. The results for mags, for example, seem better (eg more
1395 % suprathreshold voxels) for this group inversion than previously (above),
1396 % though it is less clear that the fused results are better following
1397 % group inversion of each modality alone (this is an area of current
1398 % investigation). You can decide not to use group inversion (or fusion)
1399 % of course. Again, all these statistical maps look better with more subjects,
1400 % or when using SNPM or PPM's (see Henson et al, 2007, Neuroimage and
1401 % Taylor & Henson, in prep, for more information).
1402
1403
1404 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1405 %%%%%%%%%%%%%%%%% 3D PPMs
1406 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1407
1408 % This uses SPMs above to calculate Posterior Probability Maps (PPMs),
1409 % which do not have any mutliple comparison issues (eg no need for RFT)
1410 % However, it will take a long time...
1411
1412 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1413 valdir = fullfile(stdir,sprintf('Inv%d',val));
1414 cd(valdir);
1415 for sen = 1:(Nsen+fuseinvertflag)
1416 outdir = fullfile(valdir,sennam{sen})
1417 cd(outdir)
1418 load('SPM.mat');
1419 spm_spm_Bayes(SPM);
1420 end
1421 end
1422 disp('!!!Review the group 3D PPMs if you wish!!!')
1423
1424 % Again you need to know a bit about the PPM Results interface, but press
1425 % Results and select an SPM.mat file for one of the modalities and enter
1426 % a new T-contrast [1 -1], but this time, select "Bayesian" rather than
1427 % "classical" when prompted. Then change the default effect size threshold
1428 % to 0, and probability threshold to 0.99. This will show voxels with a
1429 % >99% probability of the increase for faces relative to scrambled being
1430 % greater than zero.
1431
1432 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1433 % adding fMRI priors (only really helps IID (MNM) inversions; little effect on MSP)
1434
1435 fmri_grp_mod_evd=[];
1436 if fmriinvertflag
1437
1438 % Create fMRI priors from an SPM image....
1439 % Here we assume a group prior in MNI space, but priors could be
1440 % defined separately for each subject based on their own fMRI data
1441
1442 clear S
1443 S.D = grppreavgnam{1,sub}; % Take first subject/modality; doesnt really matter
1444 S.fmri = '/imaging/rh01/Methods/fMRIPriors/fMRI_ANOVA_flip0/UFvsS_05_cor_15.img';
1445 S.gm = '/imaging/local/spm/spm5/canonical/c1single_subj_T1.nii';
1446 S.space = 1; % SPM and GM images are in MNI space
1447 S.hthr = 0.5; S.ethr = 0;
1448 S.ncomp = Inf; S.smooth = 0.2;
1449 S.pflag = 1;
1450
1451 D0 = spm_eeg_inv_fmripriors(S);
1452 fMRI = D0.inv{D0.val}.fmri.priors; % Priors (cell of vectors)
1453
1454 for sen=1:Nsen
1455 clear DD
1456 for sub=1:Nsub
1457 DD{sub} = spm_eeg_ldata(preavgnam{sen,sub});
1458 end
1459 for val=1:Nval
1460 for sub=1:Nsub
1461 DD{sub}.inv{val}.comment = [DD{sub}.inv{val}.comment 'Group+fMRI priors'];
1462 DD{sub}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1463 DD{sub}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1464 DD{sub}.val = val;
1465 DD{sub}.inv{val}.inverse.pQ = fMRI;
1466 end
1467
1468 DD = spm_eeg_invert(DD);
1469
1470 % Save outputs
1471 for sub=1:Nsub
1472 D = DD{sub};
1473 D.fname = [D.fname(1:(end-4)) '_grp_fmri.mat'];
1474 fmrigrppreavgnam{sen,sub} = fullfile(D.path, D.fname);
1475 save(fullfile(D.path, D.fname), 'D');
1476 fmri_grp_mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end);
1477
1478 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1479
1480 D.inv{val}.contrast.woi = contwin{val}{1};
1481 D.inv{val}.contrast.fboi = contwin{val}{2};
1482 D.inv{val}.contrast.type = conengy{val};
1483 D.val = val;
1484 D = spm_eeg_inv_results(D);
1485 spm_eeg_inv_results_display(D);
1486
1487 if write3Dflag
1488 D.inv{val}.contrast.smooth = source_smooth_spm;
1489 D.inv{val}.contrast.rms = 1;
1490 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1491 D = spm_eeg_inv_Mesh2Voxels(D);
1492 FmriGrpSourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1493 end
1494 save(fullfile(D.path, D.fname), 'D');
1495 end % end val
1496 end % end sub
1497 end % end sen
1498
1499 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1500 % Fuse using Group and fMRI priors (will take a long time!)
1501 if fuseinvertflag
1502 for sub=1:Nsub
1503 clear DD;
1504 for sen=1:Nsen
1505 DD{sen} = spm_eeg_ldata(fmrigrppreavgnam{sen,sub});
1506 end
1507
1508 for val=1:Nval
1509 for sen=1:Nsen
1510 DD{sen}.val = val;
1511 DD{sen}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1512 DD{sen}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1513 % DD{sen}.inv{val}.inverse.QP = []; % (Important not to reset the source priors QP from the grp inversion above)
1514 DD{sen}.inv{val}.inverse.OutExt = '_fused';
1515 end
1516
1517 D = spm_eeg_invert_fuse(DD);
1518 fmri_grp_mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
1519
1520 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1521
1522 D.inv{val}.contrast.woi = contwin{val}{1};
1523 D.inv{val}.contrast.fboi = contwin{val}{2};
1524 D.inv{val}.contrast.type = conengy{val};
1525 D.val = val;
1526 D = spm_eeg_inv_results(D);
1527 spm_eeg_inv_results_display(D);
1528
1529 if write3Dflag
1530 D.inv{val}.contrast.smooth = source_smooth_spm;
1531 D.inv{val}.contrast.rms = 1;
1532 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1533 D = spm_eeg_inv_Mesh2Voxels(D);
1534 FmriGrpSourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1535 end
1536
1537 save(fullfile(D.path, D.fname), 'D');
1538 end % end val
1539
1540 sennam{Nsen+1} = 'fused';
1541 fmrigrppreavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
1542 end
1543 end
1544 end
1545
1546
1547 % Review model-evidences (note cannot compare fused to individual modality
1548 % inversions, since data differs)
1549 for sub = 1:Nsub
1550 subnum = dosubs(sub);
1551 disp(sprintf('\tSub %d (%d)...',subnum,sub))
1552 for val=1:Nval
1553 disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(fmri_grp_mod_evd(sub,val,:))))))
1554 end
1555 end
1556
1557 %%%%%%%%%%%%%%%%% Quick look at mean over subjects...
1558
1559 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1560 S.fig = lastfig;
1561 for val=1:Nval
1562 S.val = val;
1563 for sen=1:(Nsen+fuseinvertflag)
1564 S.fig = S.fig+1;
1565 S.P = strvcat(fmrigrppreavgnam{sen,:});
1566 meg_inv_display_con(S);
1567 drawnow
1568 end
1569 end
1570 lastfig = S.fig;
1571
1572 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1573 %%%%%%%%%%%%%%%%% 3D SPMs
1574 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1575
1576 stdir = fullfile(wd,'fMRIGrpSourceSPMs');
1577 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1578
1579 clear grpsubs imgfiles grpcons
1580 grpsubs{1} = [1:Nsub];
1581 % If want to split subjects, eg into two groups, eg:
1582 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1583 grpcons = [1:2];
1584 Ngrp = length(grpsubs);
1585 for val=1:Nval
1586 valdir = fullfile(stdir,sprintf('Inv%d',val));
1587 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
1588 for sen = 1:(Nsen+fuseinvertflag)
1589 outdir = fullfile(valdir,sennam{sen})
1590 for grp = 1:Ngrp
1591 for sub = 1:length(grpsubs{grp});
1592 subnum = grpsubs{grp}(sub);
1593 imgfiles{grp}{sub} = FmriGrpSourceImgs{sen}{val}{subnum}(grpcons,:);
1594 end
1595 end
1596 meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1597 end
1598 end
1599 cd(stdir)
1600 disp('!!!Review the group+fMRI 3D stats if you wish!!!')
1601
1602
1603 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1604 %%%%%%%%%%%%%%%%% 3D PPMs
1605 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1606
1607 % This uses SPMs above to calculate Posterior Probability Maps (PPMs),
1608 % which do not have any mutliple comparison issues (eg no need for RFT)
1609 % However, it will take a long time...
1610
1611 for val=1:Nval
1612 valdir = fullfile(stdir,sprintf('Inv%d',val));
1613 cd(valdir);
1614 for sen = 1:(Nsen+fuseinvertflag)
1615 outdir = fullfile(valdir,sennam{sen})
1616 cd(outdir)
1617 load('SPM.mat');
1618 spm_spm_Bayes(SPM);
1619 end
1620 end
1621 disp('!!!Review the group+fMRI 3D PPMs if you wish!!!')
1622
1623
1624 % Phew! You got to here! Well done!
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.You are not allowed to attach a file to this page.