%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example batch script for analysing NeuroMag VectorView Raw Data (FIF) in SPM5
% Rik Henson, 23/4/07
%
% MEG data can be found in /megdata/cbu (retrieved by MaxFilter below)
% MRI data (smri.img/hdr) needs to be downloaded from WIKI and put in a directory "sMRI"
% (together with a file smri_fids.mat coding 3 fiducials in MRI space).
% Also in the WIKI WinZip archive are the event-code files "hands.txt" and "faces.txt"
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First call NeuroMag Max Filter program (e.g, piped from Matlab):
%   1. Apply Signal Space Separation (eg to remove effects of Active Shielding)
%   2. Apply Movement correction 
%   3. Detect Bad channels (these written to log, but not yet automatically 
%      incorporated into SPM)
%   4. Downsample by a factor of 4 and convert to short
% Output is "raw.fif", and movement parameters in "raw_movement.pos"

!/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

% Could insert code to read and display raw_movement.pos (and warn about bad HPI fits)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read into Matlab and convert to SPM format 

S.Fdata = 'raw.fif'
S.rawflag = 1;                            % Raw rather than averaged data
S.Fchannels = fullfile(spm('dir'),...
 'EEGtemplates', 'FIF306_setup.mat');     % Template for channel display
S.trig_chan = 'STI101';                   % Trigger channel
S.twin = [0 Inf];                         % Timewindow to read in (i.e, all!)
S.veogchan = 2;			          % Which of two EOGs (if present) is VEOG

S.sclchan = 1;				  % Scale (gradio)meters 
% (default scaling is by distance between gradiometer coil centres, ie 1/16.8, 
% but can override with S.sclfact if like, see spm_eeg_rdata_FIF.m)
%S.sclfact = ones(1,306);
%grdscl = 1/20; 	 	          % from prestim noise estimates below
%grdscl = 1/60;			          % from one part of Neuromag manuals
%grdscl = 1/100;		          % from another part of Neuromag manuals
%S.sclfact([2:3:end 3:3:end]) = grdscl;

D = spm_eeg_rdata_FIF(S);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create two versions of data (one for Demo 1; one for Demo 2)

D.fname = 'faces.mat';
D.fnamedat = 'faces.dat';
eval('!mv raw.dat faces.dat')
save(D.fname,'D');
eval('!rm raw.mat')

D.fname = 'hands.mat';
D.fnamedat = 'hands.dat';
eval('!cp faces.dat hands.dat')
save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEMO 1 (hands) - skip if already done via GUI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Epoch

clear S; S.D = D.fname;
S.events.start = -300;
S.events.stop = 300;
S.events.types = [-2304 -1280 -768 -512];   % Because trigger codes had offset error!
S.events.Inewlist = 1;
S.events.Ec = load('hands.txt')             % New codes (1 = left, 2 = right)
D = spm_eeg_epochs(S);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Artifact (limited to rejection via thresholding of EOG; more complex possible)

clear S; S.D = D.fname
S.thresholds.External_list = 0;
S.artefact.weighted = 0;
S.thresholds.Check_Threshold = 1;
S.thresholds.threshold = [ones(1,306)*Inf 200 200];	% 200 uV rejection for EOG
D = spm_eeg_artefact(S);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Average

clear S; S.D = D.fname;
D = spm_eeg_average(S)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lowpass filter to 40Hz (usually better to do before epoching)

clear S; S.D = D.fname;
S.filter.type = 'butterworth';
S.filter.band = 'low';
S.filter.PHz = 40;
D = spm_eeg_filter(S);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results (optional)

clear S; S.D = spm_eeg_ldata(D.fname);
S.chans = [1:D.Nchannels];		    % show all channels
%S.chans = 'chans_mags_eog.mat';	    % channel file specifies mags and eog only
spm_eeg_display_ui(S);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Imaging inverse (L2-type norm), using Multiple Sparse Priors (MSP)
%
% (3 localisations done below: all channels, mags only; grads only)

badchan{1}=[];		comm{1} = 'All channels';
badchan{2}=[103:306];	comm{2} = 'Mags only';
badchan{3}=[1:102];	comm{3} = 'Grads only';

%lwin = [-300 300];	                    % time window to localise

D = spm_eeg_ldata('mae_hands.mat');	    % If want to localise evoked
%D = spm_eeg_ldata('ae_hands.mat');	    % If want to localise evoked+induced


%If want to clear previous localisations
%D = rmfield(D,'inv'); D = rmfield(D,'con');  save(D.fname,'D');

for val=1:length(comm)

 D.channels.Bad = badchan{val};
 save(D.fname,'D');

% This loop could be speeded up if could just increment D.val using spm_eeg_inv_copyfields.m
% but this is not possible if remove channels (since forward model has been be re-estimated)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise

 D.val = val;
 D.inv{val}.method = 'imaging';
 load('defaults_eeg_inv.mat')
 D.inv{val} = invdefts.imag;
 D.inv{val}.date = mat2str(fix(clock));
 D.inv{val}.comment = comm{val};
 save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI

%%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI

  D.inv{val}.mesh.Msize = 4;	% ~7200 sources
  D.inv{val}.mesh.sMRI = 'sMRI/smri.img';
  [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);

% No need to redo normalisation if done already...
  if ~exist(fullfile(pth,['wm' nam ext]))
   D = spm_eeg_inv_spatnorm(D);	
  else
   def_name      = [nam '_sn_1.mat'];
   isndef_name   = [nam '_inv_sn_1.mat'];
   D.inv{val}.mesh.def    = fullfile(pth,def_name);
   D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
   D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
   D.inv{val}.mesh.wmMRI = fullfile(pth,['wm' nam ext]);
  end
  save(D.fname,'D');

%%%%%%%%%%%%%%%%%%% Meshing 2: Warp canonical mesh

  D.inv{val}.mesh.msk_iskull = spm_select('List','sMRI','^smri_iskull.img$');

if isempty(D.inv{val}.mesh.msk_iskull);
   D = spm_eeg_inv_getmasks(D);	
else	% assume other masks exist too!
   D.inv{val}.mesh.msk_iskull = fullfile('sMRI',spm_select('List','sMRI','^smri_iskull.img$'));
   D.inv{val}.mesh.msk_cortex = fullfile('sMRI',spm_select('List','sMRI','^smri_cortex.img$'));
   D.inv{val}.mesh.msk_scalp = fullfile('sMRI',spm_select('List','sMRI','^smri_scalp.img$'));
   D.inv{val}.mesh.msk_flags  = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
end

  D.inv{val}.mesh.canonical = 1;
  D.inv{val}.mesh.template = 0;
%  D = spm_eeg_inv_meshing(D);

  Tmesh = load('wmeshTemplate_7204d.mat');
  D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
  D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
  Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
  D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
  D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
  Tmesh      = load('wmeshTemplate_scalp.mat');
  Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
  D.inv{val}.mesh.tess_scalp.vert  = Tmesh.vert;
  D.inv{val}.mesh.tess_scalp.face  = uint16(Tmesh.face);
  Tmesh      = load('wmeshTemplate_skull.mat');
  Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
  D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
  D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);

  spm_eeg_inv_checkmeshes(D);
  save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg

  if isempty(D.inv{val}.datareg.fid_coreg)
   D.inv{val}.datareg.sensors = D.channels.Loc';
   D.inv{val}.datareg.megorient = D.channels.Orient';
   D.inv{val}.datareg.fid_eeg = D.channels.fid_eeg;
   D.inv{val}.datareg.headshape = D.channels.headshape;
   cont = load('sMRI/smri_fids.mat');  name = fieldnames(cont);  D.inv{val}.datareg.fid_mri = getfield(cont,name{1});
   D.inv{val}.datareg.scalpvert = D.inv{val}.mesh.tess_scalp.vert;
   D = spm_eeg_inv_datareg(D);
  end

  spm_eeg_inv_checkdatareg(D);
  save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model

 D.inv{val}.method = 'Imaging';
 D.inv{val}.forward.sphere2fit = 2;  	% inner skull
 [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
 save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert

 D.inv{val}.inverse.trials = [1 2];
 D.inv{val}.inverse.type = 'GS';	% Actually MSP (!)
% D.inv{val}.inverse.smooth = 0.6;
% D.inv{val}.inverse.Np = 256;
% D.inv{val}.inverse.sdv = 4;
% D.inv{val}.inverse.woi = lwin;
 D.con = 1;
 D = spm_eeg_invert(D);
 save(D.fname,'D');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fancy Display

val = 1; D.val = 1;
spm_eeg_inv_visu3D_api(D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time/Freq Contrast

cwin = [0 50];	                  	% time window for contrast
D.inv{val}.contrast.woi = cwin;
D.inv{val}.contrast.fboi = 0;	   	% All frequencies	
D = spm_eeg_inv_results(D); 
D.con = 1;
save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Write Images

D.inv{val}.contrast.smooth = 12;

for con = [1 2]
  D.con = con;
  spm_eeg_invert_display(D,cwin(1)+round((cwin(2)-cwin(1))/2));
%  pause
  spm_eeg_inv_Mesh2Voxels(D);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEMO 2 (faces)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Epoch

clear S; S.D = 'faces.mat';
S.events.start = -200;
S.events.stop = 400;
S.events.resynch = 34;			% correct 2 refresh delay (60Hz) of projector
S.events.types = [-255:-247];           % Because trigger codes had offset error!
S.events.Inewlist = 1;
S.events.Ec = load('faces.txt')         % New codes (1 = face, 2 = scrambled)

D = spm_eeg_epochs(S);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Average

clear S; S.D = D.fname;
D = spm_eeg_average(S)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create new contrasts

clear S; S.D = D.fname;
S.c = [1 -1; 1/2 1/2; 1 0; 0 1]; 	% Differential ERF, common ERF, individual ERFs
S.WeightAve = 0;
D = spm_eeg_weight_epochs(S);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results (optional)

clear S; S.D = spm_eeg_ldata(D.fname);
S.chans = [1:D.Nchannels];		    % show all channels
%S.chans = 'chans_mags_eog.mat';	    % channel file specifies mags and eog only
spm_eeg_display_ui(S);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Imaging inverse (L2-type norm), using Multiple Sparse Priors (MSP)
%
% (7 localisations:  see comments (comm) below)

badchan{1}=[];	
badchan{2}=[];	
badchan{3}=[];	
badchan{4}=[103:306];
badchan{5}=[1:102];	
badchan{6}=[];	
badchan{7}=[];	

loc_cons{1} = [1];	
loc_cons{2} = [2];	
loc_cons{3} = [3 4];
loc_cons{4} = [3 4];
loc_cons{5} = [3 4];
loc_cons{6} = [1];
loc_cons{7} = [3 4];

clear lwin
lwin{1} = [-200 400];
lwin{2} = [-200 400];
lwin{3} = [-200 400];
lwin{4} = [-200 400];
lwin{5} = [-200 400];
lwin{6} = [150 200];	
lwin{7} = [150 200];	

comm{1} = 'F-S, all chan, all epoch';
comm{2} = 'F+S, all chan, all epoch';
comm{3} = 'F,S, all chan, all epoch';
comm{4} = 'F,S, mags, all epoch';
comm{5} = 'F,S, grads, all epoch';
comm{6} = 'F-S, all chan, m170';
comm{7} = 'F,S, all chan, m170';


D = spm_eeg_ldata('mme_faces.mat');	    % If want to localise evoked
%D = spm_eeg_ldata('e_faces.mat');	    % If want to localise evoked+induced

%If want to clear previous localisations
%D = rmfield(D,'inv'); D = rmfield(D,'con');  save(D.fname,'D');


for val=1:length(comm)

 D.channels.Bad = badchan{val};
 save(D.fname,'D');

% This loop could be speeded up if could just increment D.val using spm_eeg_inv_copyfields.m
% but this is not possible if remove channels (since forward model has been be re-estimated)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise

 D.val = val;
 D.inv{val}.method = 'imaging';
 load('defaults_eeg_inv.mat')
 D.inv{val} = invdefts.imag;
 D.inv{val}.date = mat2str(fix(clock));
 D.inv{val}.comment = comm{val};
 save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI

%%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI

  D.inv{val}.mesh.Msize = 4;	% ~7200 sources
  D.inv{val}.mesh.sMRI = 'sMRI/smri.img';
  [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);

% No need to redo normalisation if done already...
  if ~exist(fullfile(pth,['wm' nam ext]))
   D = spm_eeg_inv_spatnorm(D);	
  else
   def_name      = [nam '_sn_1.mat'];
   isndef_name   = [nam '_inv_sn_1.mat'];
   D.inv{val}.mesh.def    = fullfile(pth,def_name);
   D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
   D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
   D.inv{val}.mesh.wmMRI = fullfile(pth,['wm' nam ext]);
  end
  save(D.fname,'D');

%%%%%%%%%%%%%%%%%%% Meshing 2: Warp canonical mesh

  D.inv{val}.mesh.msk_iskull = spm_select('List','sMRI','^smri_iskull.img$');

if isempty(D.inv{val}.mesh.msk_iskull);
   D = spm_eeg_inv_getmasks(D);	
else	% assume other masks exist too!
   D.inv{val}.mesh.msk_iskull = fullfile('sMRI',spm_select('List','sMRI','^smri_iskull.img$'));
   D.inv{val}.mesh.msk_cortex = fullfile('sMRI',spm_select('List','sMRI','^smri_cortex.img$'));
   D.inv{val}.mesh.msk_scalp = fullfile('sMRI',spm_select('List','sMRI','^smri_scalp.img$'));
   D.inv{val}.mesh.msk_flags  = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
end

  D.inv{val}.mesh.canonical = 1;
  D.inv{val}.mesh.template = 0;
%  D = spm_eeg_inv_meshing(D);

  Tmesh = load('wmeshTemplate_7204d.mat');
  D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
  D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
  Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
  D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
  D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
  Tmesh      = load('wmeshTemplate_scalp.mat');
  Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
  D.inv{val}.mesh.tess_scalp.vert  = Tmesh.vert;
  D.inv{val}.mesh.tess_scalp.face  = uint16(Tmesh.face);
  Tmesh      = load('wmeshTemplate_skull.mat');
  Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
  D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
  D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);

  spm_eeg_inv_checkmeshes(D);
  save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg

  if isempty(D.inv{val}.datareg.fid_coreg)
   D.inv{val}.datareg.sensors = D.channels.Loc';
   D.inv{val}.datareg.megorient = D.channels.Orient';
   D.inv{val}.datareg.fid_eeg = D.channels.fid_eeg;
   D.inv{val}.datareg.headshape = D.channels.headshape;
   cont = load('sMRI/smri_fids.mat');  name = fieldnames(cont);  D.inv{val}.datareg.fid_mri = getfield(cont,name{1});
   D.inv{val}.datareg.scalpvert = D.inv{val}.mesh.tess_scalp.vert;
   D = spm_eeg_inv_datareg(D);
  end

  spm_eeg_inv_checkdatareg(D);
  save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model

 D.inv{val}.method = 'Imaging';
 D.inv{val}.forward.sphere2fit = 2;  	% inner skull
 [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
 save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert

 D.inv{val}.inverse.trials = loc_cons{val};
 D.inv{val}.inverse.type = 'GS';	% Actually MSP (!)
% D.inv{val}.inverse.smooth = 0.6;
% D.inv{val}.inverse.Np = 256;
% D.inv{val}.inverse.sdv = 4;
 D.inv{val}.inverse.woi = lwin{val};
 D.con = 1;
 D = spm_eeg_invert(D);
 save(D.fname,'D');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fancy Display

val = 3; D.val = 3;
spm_eeg_inv_visu3D_api(D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time/Freq Contrast

cwin = [150 200];	                % time window for contrast
D.inv{val}.contrast.woi = cwin;
D.inv{val}.contrast.fboi = 0;	   	% All frequencies	
D = spm_eeg_inv_results(D); 
D.con = 1;
save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Write Images

D.inv{val}.contrast.smooth = 12;

for con = [1 2]
  D.con = con;
  spm_eeg_invert_display(D,cwin(1)+round((cwin(2)-cwin(1))/2));
%  pause
  spm_eeg_inv_Mesh2Voxels(D);
end