1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   2 % Example batch script for analysing NeuroMag VectorView Raw Data (FIF) in SPM5
   3 % Rik Henson, 23/4/07
   4 %
   5 % MEG data can be found in /megdata/cbu (retrieved by MaxFilter below)
   6 % MRI data (smri.img/hdr) needs to be downloaded from WIKI and put in a directory "sMRI"
   7 % (together with a file smri_fids.mat coding 3 fiducials in MRI space).
   8 % Also in the WIKI WinZip archive are the event-code files "hands.txt" and "faces.txt"
   9 %
  10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  11 %
  12 % First call NeuroMag Max Filter program (e.g, piped from Matlab):
  13 %   1. Apply Signal Space Separation (eg to remove effects of Active Shielding)
  14 %   2. Apply Movement correction 
  15 %   3. Detect Bad channels (these written to log, but not yet automatically 
  16 %      incorporated into SPM)
  17 %   4. Downsample by a factor of 4 and convert to short
  18 % Output is "raw.fif", and movement parameters in "raw_movement.pos"
  19 
  20 !/neuro/bin/util/maxfilter -gui -f /megdata/cbu/henson/henson_rik/070412/faces_raw.fif -o raw.fif -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -format short -movecomp -hpistep 200 -hpisubt amp -hp raw_movement.pos -ds 4
  21 
  22 % Could insert code to read and display raw_movement.pos (and warn about bad HPI fits)
  23 
  24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  25 % Read into Matlab and convert to SPM format 
  26 
  27 S.Fdata = 'raw.fif'
  28 S.rawflag = 1;                            % Raw rather than averaged data
  29 S.Fchannels = fullfile(spm('dir'),...
  30  'EEGtemplates', 'FIF306_setup.mat');     % Template for channel display
  31 S.trig_chan = 'STI101';                   % Trigger channel
  32 S.twin = [0 Inf];                         % Timewindow to read in (i.e, all!)
  33 S.veogchan = 2;			          % Which of two EOGs (if present) is VEOG
  34 
  35 S.sclchan = 1;				  % Scale (gradio)meters 
  36 % (default scaling is by distance between gradiometer coil centres, ie 1/16.8, 
  37 % but can override with S.sclfact if like, see spm_eeg_rdata_FIF.m)
  38 %S.sclfact = ones(1,306);
  39 %grdscl = 1/20; 	 	          % from prestim noise estimates below
  40 %grdscl = 1/60;			          % from one part of Neuromag manuals
  41 %grdscl = 1/100;		          % from another part of Neuromag manuals
  42 %S.sclfact([2:3:end 3:3:end]) = grdscl;
  43 
  44 D = spm_eeg_rdata_FIF(S);
  45 
  46 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  47 % Create two versions of data (one for Demo 1; one for Demo 2)
  48 
  49 D.fname = 'faces.mat';
  50 D.fnamedat = 'faces.dat';
  51 eval('!mv raw.dat faces.dat')
  52 save(D.fname,'D');
  53 eval('!rm raw.mat')
  54 
  55 D.fname = 'hands.mat';
  56 D.fnamedat = 'hands.dat';
  57 eval('!cp faces.dat hands.dat')
  58 save(D.fname,'D');
  59 
  60 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  61 % DEMO 1 (hands) - skip if already done via GUI
  62 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  63 
  64 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  65 % Epoch
  66 
  67 clear S; S.D = D.fname;
  68 S.events.start = -300;
  69 S.events.stop = 300;
  70 S.events.types = [-2304 -1280 -768 -512];   % Because trigger codes had offset error!
  71 S.events.Inewlist = 1;
  72 S.events.Ec = load('hands.txt')             % New codes (1 = left, 2 = right)
  73 D = spm_eeg_epochs(S);
  74 
  75 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  76 % Artifact (limited to rejection via thresholding of EOG; more complex possible)
  77 
  78 clear S; S.D = D.fname
  79 S.thresholds.External_list = 0;
  80 S.artefact.weighted = 0;
  81 S.thresholds.Check_Threshold = 1;
  82 S.thresholds.threshold = [ones(1,306)*Inf 200 200];	% 200 uV rejection for EOG
  83 D = spm_eeg_artefact(S);
  84 
  85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  86 % Average
  87 
  88 clear S; S.D = D.fname;
  89 D = spm_eeg_average(S)
  90 
  91 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  92 % Lowpass filter to 40Hz (usually better to do before epoching)
  93 
  94 clear S; S.D = D.fname;
  95 S.filter.type = 'butterworth';
  96 S.filter.band = 'low';
  97 S.filter.PHz = 40;
  98 D = spm_eeg_filter(S);
  99 
 100 
 101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 102 % Display results (optional)
 103 
 104 clear S; S.D = spm_eeg_ldata(D.fname);
 105 S.chans = [1:D.Nchannels];		    % show all channels
 106 %S.chans = 'chans_mags_eog.mat';	    % channel file specifies mags and eog only
 107 spm_eeg_display_ui(S);
 108 
 109 
 110 
 111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 112 % Imaging inverse (L2-type norm), using Multiple Sparse Priors (MSP)
 113 %
 114 % (3 localisations done below: all channels, mags only; grads only)
 115 
 116 badchan{1}=[];		comm{1} = 'All channels';
 117 badchan{2}=[103:306];	comm{2} = 'Mags only';
 118 badchan{3}=[1:102];	comm{3} = 'Grads only';
 119 
 120 %lwin = [-300 300];	                    % time window to localise
 121 
 122 D = spm_eeg_ldata('mae_hands.mat');	    % If want to localise evoked
 123 %D = spm_eeg_ldata('ae_hands.mat');	    % If want to localise evoked+induced
 124 
 125 
 126 %If want to clear previous localisations
 127 %D = rmfield(D,'inv'); D = rmfield(D,'con');  save(D.fname,'D');
 128 
 129 for val=1:length(comm)
 130 
 131  D.channels.Bad = badchan{val};
 132  save(D.fname,'D');
 133 
 134 % This loop could be speeded up if could just increment D.val using spm_eeg_inv_copyfields.m
 135 % but this is not possible if remove channels (since forward model has been be re-estimated)
 136 
 137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise
 138 
 139  D.val = val;
 140  D.inv{val}.method = 'imaging';
 141  load('defaults_eeg_inv.mat')
 142  D.inv{val} = invdefts.imag;
 143  D.inv{val}.date = mat2str(fix(clock));
 144  D.inv{val}.comment = comm{val};
 145  save(D.fname,'D');
 146 
 147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI
 148 
 149 %%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI
 150 
 151   D.inv{val}.mesh.Msize = 4;	% ~7200 sources
 152   D.inv{val}.mesh.sMRI = 'sMRI/smri.img';
 153   [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);
 154 
 155 % No need to redo normalisation if done already...
 156   if ~exist(fullfile(pth,['wm' nam ext]))
 157    D = spm_eeg_inv_spatnorm(D);	
 158   else
 159    def_name      = [nam '_sn_1.mat'];
 160    isndef_name   = [nam '_inv_sn_1.mat'];
 161    D.inv{val}.mesh.def    = fullfile(pth,def_name);
 162    D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
 163    D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
 164    D.inv{val}.mesh.wmMRI = fullfile(pth,['wm' nam ext]);
 165   end
 166   save(D.fname,'D');
 167 
 168 %%%%%%%%%%%%%%%%%%% Meshing 2: Warp canonical mesh
 169 
 170   D.inv{val}.mesh.msk_iskull = spm_select('List','sMRI','^smri_iskull.img$');
 171 
 172 if isempty(D.inv{val}.mesh.msk_iskull);
 173    D = spm_eeg_inv_getmasks(D);	
 174 else	% assume other masks exist too!
 175    D.inv{val}.mesh.msk_iskull = fullfile('sMRI',spm_select('List','sMRI','^smri_iskull.img$'));
 176    D.inv{val}.mesh.msk_cortex = fullfile('sMRI',spm_select('List','sMRI','^smri_cortex.img$'));
 177    D.inv{val}.mesh.msk_scalp = fullfile('sMRI',spm_select('List','sMRI','^smri_scalp.img$'));
 178    D.inv{val}.mesh.msk_flags  = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
 179 end
 180 
 181   D.inv{val}.mesh.canonical = 1;
 182   D.inv{val}.mesh.template = 0;
 183 %  D = spm_eeg_inv_meshing(D);
 184 
 185   Tmesh = load('wmeshTemplate_7204d.mat');
 186   D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
 187   D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
 188   Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 189   D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
 190   D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
 191   Tmesh      = load('wmeshTemplate_scalp.mat');
 192   Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 193   D.inv{val}.mesh.tess_scalp.vert  = Tmesh.vert;
 194   D.inv{val}.mesh.tess_scalp.face  = uint16(Tmesh.face);
 195   Tmesh      = load('wmeshTemplate_skull.mat');
 196   Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 197   D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
 198   D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
 199 
 200   spm_eeg_inv_checkmeshes(D);
 201   save(D.fname,'D');
 202 
 203 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg
 204 
 205   if isempty(D.inv{val}.datareg.fid_coreg)
 206    D.inv{val}.datareg.sensors = D.channels.Loc';
 207    D.inv{val}.datareg.megorient = D.channels.Orient';
 208    D.inv{val}.datareg.fid_eeg = D.channels.fid_eeg;
 209    D.inv{val}.datareg.headshape = D.channels.headshape;
 210    cont = load('sMRI/smri_fids.mat');  name = fieldnames(cont);  D.inv{val}.datareg.fid_mri = getfield(cont,name{1});
 211    D.inv{val}.datareg.scalpvert = D.inv{val}.mesh.tess_scalp.vert;
 212    D = spm_eeg_inv_datareg(D);
 213   end
 214 
 215   spm_eeg_inv_checkdatareg(D);
 216   save(D.fname,'D');
 217 
 218 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model
 219 
 220  D.inv{val}.method = 'Imaging';
 221  D.inv{val}.forward.sphere2fit = 2;  	% inner skull
 222  [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
 223  save(D.fname,'D');
 224 
 225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert
 226 
 227  D.inv{val}.inverse.trials = [1 2];
 228  D.inv{val}.inverse.type = 'GS';	% Actually MSP (!)
 229 % D.inv{val}.inverse.smooth = 0.6;
 230 % D.inv{val}.inverse.Np = 256;
 231 % D.inv{val}.inverse.sdv = 4;
 232 % D.inv{val}.inverse.woi = lwin;
 233  D.con = 1;
 234  D = spm_eeg_invert(D);
 235  save(D.fname,'D');
 236 end
 237 
 238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fancy Display
 239 
 240 val = 1; D.val = 1;
 241 spm_eeg_inv_visu3D_api(D);
 242 
 243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time/Freq Contrast
 244 
 245 cwin = [0 50];	                  	% time window for contrast
 246 D.inv{val}.contrast.woi = cwin;
 247 D.inv{val}.contrast.fboi = 0;	   	% All frequencies	
 248 D = spm_eeg_inv_results(D); 
 249 D.con = 1;
 250 save(D.fname,'D');
 251 
 252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Write Images
 253 
 254 D.inv{val}.contrast.smooth = 12;
 255 
 256 for con = [1 2]
 257   D.con = con;
 258   spm_eeg_invert_display(D,cwin(1)+round((cwin(2)-cwin(1))/2));
 259 %  pause
 260   spm_eeg_inv_Mesh2Voxels(D);
 261 end
 262 
 263 
 264 
 265 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 266 % DEMO 2 (faces)
 267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 268 
 269 
 270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 271 % Epoch
 272 
 273 clear S; S.D = 'faces.mat';
 274 S.events.start = -200;
 275 S.events.stop = 400;
 276 S.events.resynch = 34;			% correct 2 refresh delay (60Hz) of projector
 277 S.events.types = [-255:-247];           % Because trigger codes had offset error!
 278 S.events.Inewlist = 1;
 279 S.events.Ec = load('faces.txt')         % New codes (1 = face, 2 = scrambled)
 280 
 281 D = spm_eeg_epochs(S);
 282 
 283 
 284 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 285 % Average
 286 
 287 clear S; S.D = D.fname;
 288 D = spm_eeg_average(S)
 289 
 290 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 291 % Create new contrasts
 292 
 293 clear S; S.D = D.fname;
 294 S.c = [1 -1; 1/2 1/2; 1 0; 0 1]; 	% Differential ERF, common ERF, individual ERFs
 295 S.WeightAve = 0;
 296 D = spm_eeg_weight_epochs(S);
 297 
 298 
 299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 300 % Display results (optional)
 301 
 302 clear S; S.D = spm_eeg_ldata(D.fname);
 303 S.chans = [1:D.Nchannels];		    % show all channels
 304 %S.chans = 'chans_mags_eog.mat';	    % channel file specifies mags and eog only
 305 spm_eeg_display_ui(S);
 306 
 307 
 308 
 309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 310 % Imaging inverse (L2-type norm), using Multiple Sparse Priors (MSP)
 311 %
 312 % (7 localisations:  see comments (comm) below)
 313 
 314 badchan{1}=[];	
 315 badchan{2}=[];	
 316 badchan{3}=[];	
 317 badchan{4}=[103:306];
 318 badchan{5}=[1:102];	
 319 badchan{6}=[];	
 320 badchan{7}=[];	
 321 
 322 loc_cons{1} = [1];	
 323 loc_cons{2} = [2];	
 324 loc_cons{3} = [3 4];
 325 loc_cons{4} = [3 4];
 326 loc_cons{5} = [3 4];
 327 loc_cons{6} = [1];
 328 loc_cons{7} = [3 4];
 329 
 330 clear lwin
 331 lwin{1} = [-200 400];
 332 lwin{2} = [-200 400];
 333 lwin{3} = [-200 400];
 334 lwin{4} = [-200 400];
 335 lwin{5} = [-200 400];
 336 lwin{6} = [150 200];	
 337 lwin{7} = [150 200];	
 338 
 339 comm{1} = 'F-S, all chan, all epoch';
 340 comm{2} = 'F+S, all chan, all epoch';
 341 comm{3} = 'F,S, all chan, all epoch';
 342 comm{4} = 'F,S, mags, all epoch';
 343 comm{5} = 'F,S, grads, all epoch';
 344 comm{6} = 'F-S, all chan, m170';
 345 comm{7} = 'F,S, all chan, m170';
 346 
 347 
 348 D = spm_eeg_ldata('mme_faces.mat');	    % If want to localise evoked
 349 %D = spm_eeg_ldata('e_faces.mat');	    % If want to localise evoked+induced
 350 
 351 %If want to clear previous localisations
 352 %D = rmfield(D,'inv'); D = rmfield(D,'con');  save(D.fname,'D');
 353 
 354 
 355 for val=1:length(comm)
 356 
 357  D.channels.Bad = badchan{val};
 358  save(D.fname,'D');
 359 
 360 % This loop could be speeded up if could just increment D.val using spm_eeg_inv_copyfields.m
 361 % but this is not possible if remove channels (since forward model has been be re-estimated)
 362 
 363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise
 364 
 365  D.val = val;
 366  D.inv{val}.method = 'imaging';
 367  load('defaults_eeg_inv.mat')
 368  D.inv{val} = invdefts.imag;
 369  D.inv{val}.date = mat2str(fix(clock));
 370  D.inv{val}.comment = comm{val};
 371  save(D.fname,'D');
 372 
 373 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI
 374 
 375 %%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI
 376 
 377   D.inv{val}.mesh.Msize = 4;	% ~7200 sources
 378   D.inv{val}.mesh.sMRI = 'sMRI/smri.img';
 379   [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);
 380 
 381 % No need to redo normalisation if done already...
 382   if ~exist(fullfile(pth,['wm' nam ext]))
 383    D = spm_eeg_inv_spatnorm(D);	
 384   else
 385    def_name      = [nam '_sn_1.mat'];
 386    isndef_name   = [nam '_inv_sn_1.mat'];
 387    D.inv{val}.mesh.def    = fullfile(pth,def_name);
 388    D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
 389    D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
 390    D.inv{val}.mesh.wmMRI = fullfile(pth,['wm' nam ext]);
 391   end
 392   save(D.fname,'D');
 393 
 394 %%%%%%%%%%%%%%%%%%% Meshing 2: Warp canonical mesh
 395 
 396   D.inv{val}.mesh.msk_iskull = spm_select('List','sMRI','^smri_iskull.img$');
 397 
 398 if isempty(D.inv{val}.mesh.msk_iskull);
 399    D = spm_eeg_inv_getmasks(D);	
 400 else	% assume other masks exist too!
 401    D.inv{val}.mesh.msk_iskull = fullfile('sMRI',spm_select('List','sMRI','^smri_iskull.img$'));
 402    D.inv{val}.mesh.msk_cortex = fullfile('sMRI',spm_select('List','sMRI','^smri_cortex.img$'));
 403    D.inv{val}.mesh.msk_scalp = fullfile('sMRI',spm_select('List','sMRI','^smri_scalp.img$'));
 404    D.inv{val}.mesh.msk_flags  = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
 405 end
 406 
 407   D.inv{val}.mesh.canonical = 1;
 408   D.inv{val}.mesh.template = 0;
 409 %  D = spm_eeg_inv_meshing(D);
 410 
 411   Tmesh = load('wmeshTemplate_7204d.mat');
 412   D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
 413   D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
 414   Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 415   D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
 416   D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
 417   Tmesh      = load('wmeshTemplate_scalp.mat');
 418   Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 419   D.inv{val}.mesh.tess_scalp.vert  = Tmesh.vert;
 420   D.inv{val}.mesh.tess_scalp.face  = uint16(Tmesh.face);
 421   Tmesh      = load('wmeshTemplate_skull.mat');
 422   Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
 423   D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
 424   D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
 425 
 426   spm_eeg_inv_checkmeshes(D);
 427   save(D.fname,'D');
 428 
 429 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg
 430 
 431   if isempty(D.inv{val}.datareg.fid_coreg)
 432    D.inv{val}.datareg.sensors = D.channels.Loc';
 433    D.inv{val}.datareg.megorient = D.channels.Orient';
 434    D.inv{val}.datareg.fid_eeg = D.channels.fid_eeg;
 435    D.inv{val}.datareg.headshape = D.channels.headshape;
 436    cont = load('sMRI/smri_fids.mat');  name = fieldnames(cont);  D.inv{val}.datareg.fid_mri = getfield(cont,name{1});
 437    D.inv{val}.datareg.scalpvert = D.inv{val}.mesh.tess_scalp.vert;
 438    D = spm_eeg_inv_datareg(D);
 439   end
 440 
 441   spm_eeg_inv_checkdatareg(D);
 442   save(D.fname,'D');
 443 
 444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model
 445 
 446  D.inv{val}.method = 'Imaging';
 447  D.inv{val}.forward.sphere2fit = 2;  	% inner skull
 448  [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
 449  save(D.fname,'D');
 450 
 451 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert
 452 
 453  D.inv{val}.inverse.trials = loc_cons{val};
 454  D.inv{val}.inverse.type = 'GS';	% Actually MSP (!)
 455 % D.inv{val}.inverse.smooth = 0.6;
 456 % D.inv{val}.inverse.Np = 256;
 457 % D.inv{val}.inverse.sdv = 4;
 458  D.inv{val}.inverse.woi = lwin{val};
 459  D.con = 1;
 460  D = spm_eeg_invert(D);
 461  save(D.fname,'D');
 462 
 463 end
 464 
 465 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fancy Display
 466 
 467 val = 3; D.val = 3;
 468 spm_eeg_inv_visu3D_api(D);
 469 
 470 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time/Freq Contrast
 471 
 472 cwin = [150 200];	                % time window for contrast
 473 D.inv{val}.contrast.woi = cwin;
 474 D.inv{val}.contrast.fboi = 0;	   	% All frequencies	
 475 D = spm_eeg_inv_results(D); 
 476 D.con = 1;
 477 save(D.fname,'D');
 478 
 479 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Write Images
 480 
 481 D.inv{val}.contrast.smooth = 12;
 482 
 483 for con = [1 2]
 484   D.con = con;
 485   spm_eeg_invert_display(D,cwin(1)+round((cwin(2)-cwin(1))/2));
 486 %  pause
 487   spm_eeg_inv_Mesh2Voxels(D);
 488 end