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.You are not allowed to attach a file to this page.