% Batch script (linear pipeline) for MEG analysis in SPM5   
% Rik Henson MRC CBU Sep 08
% Latest Revision July 09
%
% Note that would benefit from modularisation (eg for stop/start
% running), but wanted all steps in one file to ease exposition
%
% First start SPM with (at the CBU) "spm 5 eeg" from Linux prompt

clear

wd = '/imaging/rh01/Methods/CBUMEGDemo/Group';   % !!!REPLACE WITH YOUR DIRECTORY!!!
cd(wd)

addpath /imaging/local/spm/spm5/cbu_updates/    
% Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/spm5_cbu_updates/devel/
addpath /imaging/local/meg_misc/    
% Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/meg_misc/devel/

%%%%%%%%%% Different types of sensor
sennam{1}   = 'mags'; sennam{2} = 'grds'; sennam{3} = 'eeg';
Nsen     = length(sennam);
eegflag  = [0 0 1];         % Binary flag for which sensor is EEG

% Num channels per sensor type (if differs across subjects, can re-specify below)
Nchan(1) = 102;    Nchan(2) = 204;   Nchan(3) = 70;

% Any user-thresholds for rejection for each of above (Inf=none)
thr(1)   = Inf;    thr(2)   = Inf;   thr(3)   = 120;
EOGthr   = 120;   % EOG threshold (EOG is duplicated in each file)

%%%%%%%%%% Experiment-specific parameters
expwin      = [-100 400];                    % Epoch window (in ms)
expcodes    = [1:8];                         % Codes in FIF file to extract
offset      = [ones(1,length(expcodes))*34]; % t=0 offset (eg visual delay)
newcodes    = [1 1 1 1 2 2 2 2];             % Any new (recoding) of above
expcons     = [1 0; 0 1; 1 -1; 1/2 1/2];     % Contrasts (on NEW codes)
Ncons       = size(expcons,1); 
		 
%%%%%%%%%% Raw FIF files on network 
% (cell of cells for each session of each sub (only one session here))
% Eg if instead two sessions per subject:
%    Fifdata = { { '/megdata/cbu/lfp/meg08_0014/080117/block1_raw.fif';
%                  '/megdata/cbu/lfp/meg08_0014/080117/block2_raw.fif'} };
fifdata = {
           {'/megdata/cbu/lfp/meg08_0014/080117/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0015/080118/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0016/080118/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0038/080131/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0064/080212/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0111/080303/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0112/080303/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0113/080303/block4_raw.fif'};
           {'/megdata/cbu/lfp/meg08_0127/080310/block4_raw.fif'};
	     }


%%%%%%%%%% Raw DICOM files on network (one per subject)
mridata = {
           '/mridata/cbu/CBU080769_CBU080769/20080916_095501/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU071255_CBU071255/20071218_150851/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU071151_CBU071151/20071116_140044/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU071229_CBU071229/20071211_112249/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU080034_RIS13667/20080111_102845/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU080292_CBU080292/20080418_161843/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU080358_CBU080358/20080506_151233/Series_003_CBU_MPRAGE';
           '/mridata/cbu/CBU080413_CBU080413/20080522_124313/Series_002_CBU_MPRAGE';
           '/mridata/cbu/CBU080802_CBU080802/20080930_090538/Series_002_CBU_MPRAGE';
	  };


%%%%%%%%%% Bad channel names noticed by Operator or post hoc (eeg/meg,sub,ses)
badchans{1,1,1} = [1143];
badchans{1,2,1} = [1143 0913 1013 0313 0933 2222 0642 1232 1913];
badchans{1,3,1} = [1143];
badchans{1,4,1} = [1143 2223 0432];
badchans{1,5,1} = [1143];
badchans{1,6,1} = [1143];
badchans{1,7,1} = [1143 2223];
badchans{1,8,1} = [1143 1412];
badchans{1,9,1} = [1143 2223];

%%%%%%%%%% Bad EEG channels noticed by Operator or post hoc (in terms of NAME)
badchans{2,1,1} = [48 50 73];                 % Poor contact
badchans{2,2,1} = [7 13 16 24 46 48 57];      % Drifting
badchans{2,3,1} = [48];                       % Drifting
badchans{2,4,1} = [48];                       % Drifting
badchans{2,5,1} = [15 40 49];                 % Poor contact
badchans{2,6,1} = [48];                       % Drifting
badchans{2,7,1} = [48 72];                    % No deflection (odd N170 topo)?
badchans{2,8,1} = [65 48];                    % odd in final topo of N170?
badchans{2,9,1} = [20 31 45];                 % 45 poor contact; 20+31 highfreq noise

% (Just happens that any EEG Channels 61 onwards are numbered 65 onwards in our lab!)
for sub=1:size(badchans,2)
  for ses=1:size(badchans,3)
   f = find(badchans{2,sub,ses}>60);
   badchans{2,sub,ses}(f) = badchans{2,sub,ses}(f)-4;
  end
end

%%%%%%%%%%% subject-specific parameters
%dosubs   = [1:9];      % Subjects to run, eg for subset: dosubs = [5:8]; 
dosubs   = [1:8];       % Exclude sub9 because MRI segmentation requires manual re-
                        % positioning first (can try yourself; see MRI section below!)
                        % (in fact for many, eg sub5 too, normalisation is much
                        % better after manual repositioning)
Nsub     = length(dosubs);

% Below assumes that two EOG channels (VEOG and HEOG), which will be
% appended at end of every data-file. Note that this script does not yet
% handle concurrent ECG (see Jason Taylor for how to do this)
for s = 1:Nsub
 subnum = dosubs(s);
 for m = 1:Nsen;
   Nsubchan(s,m) = Nchan(m);  % In case different number of channels (eg EEG) for some subjects
  chan_thr{s,m} = [ones(1,Nsubchan(s,m))*thr(m) EOGthr EOGthr]; 
 end
end

VitECapsules = ones(1,max(dosubs));  % Whether MRIs have Vitamin E capsule (all here)

%%%%%%%%%%% control flags (1=do step; 0=omit step)
maxfiltflag = 1; maxmvcompflag = 1; maxtransflag = 1; 
usemaxmvcompflag = 0; usemaxtransflag = 1; 
chancheckflag = 0; respmflag = 1; reprocessflag = 1; 
getmriflag = 1;
write2Dflag = [1 0 1]; write3Dflag = 1;
invertflag = [1 1 1];  % Invert all modalities
fuseinvertflag = 1;    % Whether want to simultaneously invert (fuse) modalities (sennam)
                       % (Henson et al, 2009b)
grpinvertflag = 1;     % Whether invert using group priors (after individual inversions)
                       % (Litvak & Friston, 2008, Neuroimage)
fmriinvertflag = 1;    % Whether want to use fMRI priors on inversion
                       % (Henson et al, submitted)
		       
%%%%%%%%%% General parameters
% Maxfilter
downsample_maxf   = '-ds 4';           % maxfilter downsampling factor (eg 4 to save space)
                                       % though beware of effect of downsampling on triggers
				       % (http://imaging.mrc-cbu.cam.ac.uk/meg/CbuSpmParameters)
%downsample_maxf  = '';                % (if none)
format_maxf       = '-format float';   % maxfilter data format (native = float32)
%format_maxf      = '-format short';   % (if to save space)
trans_offset_maxf = [0 -13 +6];        % New centre for trans default
                                       % (see our maxfilter WIKI)

% SPM
downsample_spm    = 100;                % spm new sample rate (Hz; should be at
                                        % least twice filt_cut_spm below)
filt_cut_spm      = 40;                 % Lowpass filter cutoff for SPM (Hz)
dformat_spm       = 'int16';            % data format for spm (eg to save space)
sensor_smooth_spm = [5 5 20];           % Smoothing in SensorSpace-Time [mm mm ms];
source_smooth_spm = [10 10 10];         % Smoothing in Source space [mm mm mm]

%%%%%%%%%%% Initialisation
subsphere = [];
allsssbad = cell(Nsub,1);
nbadchan = cell(Nsen,Nsub);
nrejects = cell(Nsen,Nsub);
nevents  = cell(Nsen,Nsub);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Preprocess MEG/EEG data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for sub = 1:Nsub
 
 subnum    = dosubs(sub);
 Nses(sub) = size(fifdata{subnum},1);

 cd(wd)
 swd = sprintf('sub%d',subnum);
 try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end
  
 for ses = 1:Nses(sub)

   allsssbad{sub} = cell(1,Nses(sub));
   
   fnm_stem = sprintf('sub%d_ses%d',subnum,ses);

   if maxfiltflag

%%%%% Use autobad to detect bad channels from first 20s, and select from output
     rawfile  = fifdata{subnum}{ses};
     sssfile  = sprintf('%s_sss.fif',fnm_stem)
     logfile  = sprintf('%s_autobad.log',fnm_stem)
     headptsfile  = sprintf('%s_headpoints.txt',fnm_stem)

     if exist(sssfile), delete(sssfile), end
     if exist(logfile), delete(logfile), end
     [Centre,Radius]  = meg_fit_sphere(rawfile,pwd,headptsfile);
     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))
     subsphere(sub,:) = [Centre Radius];

     disp('Checking first 20s for bad channels...')
     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))
     
     badfile = sprintf('%s_badchan.txt',fnm_stem)
     if exist(badfile), delete(badfile), end
     eval(sprintf('!cat %s | sed -n ''/Static/p'' | cut -f 5- -d '' '' > %s',logfile,badfile));
     sssbadchans{sub}{ses} = unique([textread(badfile,'%d')' badchans{1,subnum,ses}]);
     for bc=1:length(sssbadchans{sub}{ses}); allsssbad{sub}{ses} = [allsssbad{sub}{ses} sprintf(' %04d',sssbadchans{sub}{ses}(bc))]; end

%%%%% Now mark bad channels and run ST 
     logfile = sprintf('%s_ssst.log',fnm_stem)
     if exist(sssfile), delete(sssfile), end    % Don't need previous file from autobad
     sssfile  = sprintf('%s_ssst.fif',fnm_stem)
     if exist(logfile), delete(logfile), end
     posfile = sprintf('%s_mvpos.txt',fnm_stem)

     disp(['Now maxfiltering with bad channels (and tSSS): ' allsssbad{sub}{ses}])
     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))

     S.logfile = logfile;
     [max_d(sub),max_r(sub)] = meg_plot_maxfilter_move(S);
% Do check the movement parameters - particularly goodness of fit -
% because movement compensation below will be invalidated if poor fit
% (also, the 'inter' MaxFilter option is used below in case the HPI 
% temporarily fails (which it does for sub5), but if your subject's 
% HPI failed  permanently from the start, you might not know this otherwise,
% because 'inter' would handle it fine!).
     
     infile = sssfile;
     
%%%%% Now mark bad channels and run MVCOMP
     if maxmvcompflag

      sss2file = sprintf('%s_mvcomp.fif',fnm_stem)
      logfile = sprintf('%s_mvcomp.log',fnm_stem)
      if exist(sss2file), delete(sss2file), end
      if exist(logfile), delete(logfile), end

      disp(['Now compensating movement... (every 1000ms)'])
      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))
%      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))

      infile = sss2file;
     end
     
%%%%% Now Trans Default (if requested)
     if maxtransflag
      trans_origin = Centre(1:3) + trans_offset_maxf;
      logfile = sprintf('%s_trans.log',fnm_stem)
      sss3file = sprintf('%s_trans.fif',fnm_stem)
      if exist(sss3file), delete(sss3file), end
      if exist(logfile), delete(logfile), end

      disp(['Now trans-ing to DEFAULT ' num2str(round(trans_origin))])
      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))
     end
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Olaf's channel checking (have not used for a while; should put in SVN!)
    if chancheckflag 
     addpath /imaging/olaf/MEG/artefact_scan_tool/
     global fiffs tasks values criteria sensors datapara figs;
     ini_check4artefacts; % Initialisation of variables
     epochlength = 5;    
     amplitude_threshold_mag = 5000;  
     amplitude_threshold_gra = 10000;         
     out_dir = fullfile(wd,'ChanCheck');
     addfiff(infile,fnm_stem);
     addtask('maxminamp');   % maximum minus minimum amplitude within epoch
     addtask('maxmingrad');  % maximum signal gradient (amplitude difference between successive data points in time) with epoch
     addtask('threshold');   % number of epochs with max-min amplitudes above specified threshold
     addtask('crosscorr');   % intercorrelation matrix for magnetometers and gradiometers separately
     addtask('trendamp');    % linear trend in max-min differences across all epochs
     addtask('trendgrad');   % linear trend in maximum signal gradients across all epochs
     check4artefacts;    % That's where it all happens... lean back and enjoy!
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Read into SPM (writes separate files for MEG and EEG)

    clear S
    if usemaxtransflag
      S.Fdata = sprintf('%s_trans.fif',fnm_stem);
      if maxmvcompflag
       S.HPIfile = sprintf('%s_mvcomp.fif',fnm_stem);   % If trans'ed (hpi lost)
      else
       S.HPIfile = sprintf('%s_sss.fif',fnm_stem);      % If trans'ed (hpi lost)
      end
    elseif usemaxmvcompflag
      S.Fdata = sprintf('%s_mvcomp.fif',fnm_stem);
      S.HPIfile = [];
    else
      S.Fdata = sprintf('%s_sss.fif',fnm_stem);      
      S.HPIfile = [];
    end
%    meg_plot_fifpoints(S);                              % If want to check points
    S.Fchannels = 'FIF306_setup.mat';
    S.veogchan  = [62 63]; S.veog_unipolar = 0;
    S.heogchan  = [61 64]; S.heog_unipolar = 0;
    S.trig_chan = 'STI101';
    S.twin      = [0 Inf];
    S.eeg_ref   = 'average';
    S.trig_type = 'positive_from_zero';
%    S.Dtype     = 'float32';                             % make int16 if files too big
    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
%    S.Fchannels_eeg =  fullfile(wd,swd,sprintf('eeg_montage_sub%02d.mat',subnum)); % use subject-specific (care labels may not be standard)
    
    if respmflag
      D = spm_eeg_rdata_FIF(S);
      basefile = D.fname(1:(end-4));
    else
      basefile = S.Fdata(1:(end-4));
    end

%%%%% Downsample, then split Mags and Grads (file too big otherwise)

    if reprocessflag
      clear S; S.D = sprintf('%s-eeg.mat',basefile); disp(S.D)
      
%%%%% Take opportunity to relabel EEG channels according to approx 10-20 system
%%%%% (see http://imaging.mrc-cbu.cam.ac.uk/meg/EEGChannelNames)
      D = spm_eeg_ldata(S.D); 
      D.channels.ctf = fullfile(spm('dir'),'EEGtemplates', sprintf('cbu_meg_%deeg_montage_named.mat',Nsubchan(sub,find(eegflag))));
      tmp = load(D.channels.ctf);
      D.channels.name = tmp.Cnames;           % NB: This assumes channel order 
            		                      % has not been changed from 1-70!
      spm_eeg_inv_save(D);

      
      S.Radc_new = downsample_spm;
      D = spm_eeg_downsample(S);

      clear S; S.D = basefile;
      S.Radc_new = downsample_spm;
      D = spm_eeg_downsample(S);
      basefile = D.fname(1:(end-4));

      clear S; S.D = D.fname;
      D = spm_eeg_splitFIF(S);
    end
    
%%%%% If want to view rawdata (with Danny's tool):
%     meg_viewdata(D.fname);
     
%%%%% Preprocess each sensor type

     for sen = 1:Nsen
 
       if reprocessflag
	clear S
        S.D = sprintf('%s-%s.mat',basefile,sennam{sen}); disp(S.D)
	S.filter.type = 'butterworth';
	S.filter.band = 'low';	S.filter.PHz = filt_cut_spm;
%	S.filter.band = 'stop'; S.filter.PHz = [49 51];
        S.Dtype       = dformat_spm;
	D = spm_eeg_filter(S);

%%%%%% !!!Insert ICA here if want to remove EOG/ECG artifacts!!!
%%%%%% (See Jason Taylor; assume filtered, downsampled data appropriate)
	
	clear S; S.D = D.fname;
 	S.events.start    = expwin(1);
 	S.events.stop     = expwin(2);
 	S.events.types    = expcodes;
 	S.events.Inewlist = 0;
 	S.events.resynch  = offset;
	D = spm_eeg_epochs(S);
	
	% re-classify event-codes
	for e = 1:length(expcodes);
	    fi = find(D.events.code==expcodes(e));
	    if ~isempty(fi)
		D.events.code(fi) = newcodes(e);
	    end
	end
	D.events.types  = unique(D.events.code);
	D.events.Ntypes = length(D.events.types);
	save(D.fname,'D')
	clear S; S.D = D.fname;

      else
        S.D = sprintf('e_fd%s-%s.mat',basefile,sennam{sen}); disp(S.D)
	D = spm_eeg_ldata(S.D);
      end
      
%%%%% Threshold channels (eg EOG)

      nc = D.Nchannels;
      S.thresholds.External_list   = 0;
      S.thresholds.Check_Threshold = 1;
      S.thresholds.threshold       = chan_thr{sub,sen};
      S.BadThr                     = 0.2;  
      S.artefact.weighted          = 0;
      if eegflag(sen)==0
%           S.channels.Bad = unique([D.channels.Bad strcmpi(D.channels.name....sssbadchans{sub}{ses}]);     	% If dont trust MaxFilter channel reconstruction
      else
            S.channels.Bad = unique([D.channels.Bad badchans{2,subnum,ses}]);  
      end
      D = spm_eeg_artefact(S);
      
      nbadchan{sen,sub} = [nbadchan{sen,sub} length(D.channels.Bad)];	
      nrejects{sen,sub} = [nrejects{sen,sub} length(find(D.events.reject))];

%%%%% Re-reference EEG if removed channels
      if eegflag(sen) & ~isempty(D.channels.Bad)
	   D = spm_eeg_ldata(D.fname);
           meg_reavg_ref(D);
	   save(D.fname,'D')
      end
      sesnam{sen,ses} = D.fname;
      
%%%%% Uncomment if want to average and contrast each session separately 
%%%%% (as well as after merging below)
%
%      clear S; S.D = D.fname;
%      D = spm_eeg_average(S);
%
%      clear S; S.D = D.fname;
%      S.c         = expcons;
%      S.WeightAve = 1;
%      D = spm_eeg_weight_epochs(S);
%      nevents{sen,sub} = D.events.repl;
      
    end  % end of sen loop
  end  % end of ses loop
    
%%%%% Merge sessions, average, contrast 

  for sen = 1:Nsen
    if Nses(sub)>1
      if ~usemaxtransflag
        disp(['Concatenating sessions... (NB: merged files will not have correct channel Loc/Orient info unless trans-ed to first or to default!!)']);
      end
      clear S; S.D = strvcat(sesnam{sen,:});
      S.recode = cell(1,Nses(sub));
      D = spm_eeg_merge(S);
      clear S; S.D = D.fname;
    else
      clear S; S.D = sesnam{sen,1};
      D = spm_eeg_ldata(S.D);
    end
    preavgnam{sen,sub} = fullfile(D.path,S.D);  
    Nsubchan(sub,sen) = length(D.channels.eeg) - length(D.channels.Bad);
  end
  
 
%%%% If want to take RMS of grds before averaging (for Sensor-Time image files)
%%%% Noise will not cancel as well, but results less biased by number of trials
%%%% (for more info, see end of http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm)
%
% for sen = 1:Nsen 
%    if strcmp(sennam{sen},'grds')
%      clear S; S.D = preavgnam{sen,sub};
%      S.bc = 1;
%      S.do_write = 1;
%      D = meg_grds2grms(S); % RMS of grads; just for sensor-level analyses 
%      sennam{Nsen+1} = 'grms';
%      preavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
%      nevents{Nsen+1,sub}   = nevents{sen,sub};
%      nrejects{Nsen+1,sub}  = nrejects{sen,sub};
%      nbadchan{Nsen+1,sub}  = nbadchan{sen,sub};
%      Nsubchan(sub,Nsen+1)  = Nsubchan(sub,sen);
%      write2Dflag(Nsen+1)   = 1;
%      invertflag(Nsen+1)    = 0;
%    end
%  end
  
  for sen = 1:length(sennam)

    clear S; S.D = preavgnam{sen,sub};
    D = spm_eeg_average(S);

    clear S; S.D = D.fname;
    S.c = expcons;
    S.WeightAve = 0;
    D = spm_eeg_weight_epochs(S);
    
    nevents{sen,sub}  = D.events.repl;
    avgnam{sen,sub} = fullfile(D.path,D.fname);
  end
end

cd(wd)

for sen = 1:length(sennam)
 disp(sprintf('%s...',sennam{sen}))
 for sub = 1:Nsub
   subnum = dosubs(sub);
   disp(sprintf('\tSub %d (%d)...',subnum,sub))
   disp(sprintf('\t\tNevents (per condition): %s',mat2str(nevents{sen,sub})))
   disp(sprintf('\t\tRejects (per session): %s',mat2str(nrejects{sen,sub})))
   disp(sprintf('\t\tBadchans (per session): %s',mat2str(nbadchan{sen,sub})))
 end
end

%%%% If want to take RMS after averaging (cf. commented section above)
%%%% In theory more sensitive, but biased by trial numbers and less Gaussian
%%%% (See end of http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm page)

for sen = 1:Nsen
  if strcmp(sennam{sen},'grds')
     for sub = 1:Nsub
       cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
       clear S; S.D = avgnam{sen,sub};
       S.bc = 1;
       S.do_write = 1;
       S.log = 1;   % If your RMS are skewed (see above WIKI page)
       D = meg_grds2grms(S);
       avgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
       nevents{Nsen+1,sub}  = nevents{sen,sub};
     end
     Nsen=Nsen+1; sennam{Nsen} = 'grms';
     write2Dflag(Nsen) = 1;
  end
end

%%%% Write out Sensor-Time image files for each subject

for sen = find(write2Dflag)
  for sub = 1:Nsub
         cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
	 disp(avgnam{sen,sub})
         [pth,nam,ext] = fileparts(avgnam{sen,sub});
         clear S; S.Fname = [nam ext]; 
	 S.interpolate_bad = 0;
         S.n = 32;
         S.pixsize = 3;
         S.trialtypes = [1:Ncons];
         P = spm_eeg_convertmat2ana3D(S);
	 if any(sensor_smooth_spm)   
	  if length(sensor_smooth_spm)==1; sensor_smooth_spm = ones(1,3)*sensor_smooth_spm; end
	  for con = 1:size(P,1);
	   [pth,nam,ext] = fileparts(P(con,:));
	   Pout          = fullfile(pth,['s' nam ext]);
           spm_smooth(spm_vol(P(con,:)),Pout,sensor_smooth_spm);
	   Pin = strvcat(P(con,:),Pout);
	   spm_imcalc_ui(Pin,Pout,'((i1+eps).*i2)./(i1+eps)',{[],[],'float32',0});   %Reinsert NaNs
 	   SensorTimeImgs{sen}{sub}(con,:) = fullfile(wd,swd,Pout);
	  end
	 else
	  for con = 1:size(P,1);
 	   SensorTimeImgs{sen}{sub}(con,:) = fullfile(wd,swd,P(con,:));
	  end
	 end
  end
end

%%%% Create Grand Mean across subjects
cd(wd)
for sen = 1:length(sennam)
 clear S; S.P = strvcat(avgnam{sen,:});
 S.Pout = sprintf('G%d-%s.mat',Nsub,sennam{sen})
 D = spm_eeg_grandmean(S);
end

eval(sprintf('save preproc%d.mat subsphere nrejects nsubchan nbadchan nevents allsssbad dosubs preavgnam avgnam',length(dosubs)))

disp('!!!NOW USE SPM-Display-M/EEG button to examine results, eg G8-*.mat files !!!')

%return   % uncommment if want script to finish here


% Once displaying one of the modalities, you can hold shift and select both 
% conditions 1+2, and click on a posterior right channel to see the M170
% difference between faces (1) and scrambled faces (2). You can also select
% condition 3, which is the differential ERF/ERP between faces and scrambled,
% then press "topography", enter -100ms for initial time, then press and hold
% the time slider to see topography of difference at every time point (displayed
% in Matlab window) - random noise until ~100ms, when strong differences emerge


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 2D sensor x 1D Time SPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stdir = fullfile(wd,'SensorTimeSPMs');
try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end

clear grpsubs grpcons imgfiles
grpsubs{1} = [1:Nsub];
% If want to split subjects, eg into two groups, eg:
% grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
grpcons{1} = [1:2];  % Include each condition for mags
grpcons{3} = [1:2];  % Include each condition for eeg
grpcons{4} = [1:2];  % Diff of RMS (be careful if diff number of trials for diff conditions)
%grpcons{4} = [3];    % Or RMS of diff (not equivalent to above)

Ngrp = length(grpsubs);
for sen = find(write2Dflag)
 outdir = fullfile(stdir,sennam{sen})
 for grp = 1:Ngrp
  for sub = 1:length(grpsubs{grp});
   subnum = grpsubs{grp}(sub);
   imgfiles{grp}{sub} = SensorTimeImgs{sen}{subnum}(grpcons{sen},:);
  end
 end
 meg_batch_anova(imgfiles,outdir);
end

disp('!!!NOW press Results button in SPM!!!')

%return

% You need to know a bit about the SPM Results interface, but press the Results 
% button (make sure SPM is in "EEG" mode) (or better still, use the "ns" toolbox
% under "Toolbox" button - see WIKI page below) and select the SPM.mat file for 
% one of the modalities, select the default effects of interest F-contrast
% (because we don't care about polarity, at least for Mags and EEG), choose
% 0.001 uncorrected. Then press "Volume" to get table of stats (ignore brain image). 
% You should see, at least in Mags and EEG, frontal and posterior clusters 
% (simply of different polarity) maximal around 170ms. Few individual
% "voxels" survive whole-image correction because only 8 subjects (so low
% df's, ie low power, and RFT conservative) - though the space-time extent
% of clusters around 170ms do survive correction at the cluster level
% for Mags and EEG (using "ns" toolbox); For GRMS, you can evaluate a
% T-contrast [1 -1] for Faces>Scrambled, threshold at p<001 uncorrected 
% (using ns toolbox) and find some corrected clusters around 170ms,
% though the results are not as convincing as Mags or EEG (probably
% related to RMS issues: see end of WIKI page below)
%
% Correction for height would also be significant if you used SVC because
% of an a priori timewindow of interest around 170ms (again see WIKI page below).
%
% For more info on interpreting these SPMs, see
%                   http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm
%
% Also, for an example use with EEG data, see Henson et al (2008), Neuroimage.
%
% Note that you can do SVC for a priori timewindows of interest, but the main 
% use of such SPMs is to define timewindows of interest (in the absence of 
% a priori ones) for the subsequent source localisation below...
%
% (Note that these SPMs can also handle different numbers and locations of channels
% across subjects, provided their locations are digitised in same space, because the 
% data are interpolated to same grid. For MEG data, particularly gradiometers, the
% results benefit from using 'trans -default' above, to realign across different
% headpositions for different subjects - see Smith et al (2009), HBM Abstract.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Get MRIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for sub=1:Nsub 
   subnum = dosubs(sub);
   swd = sprintf('sub%d',subnum);
   mwd{subnum} = fullfile(wd,swd,'MRI');
   try cd(mwd{subnum}); catch eval(sprintf('!mkdir %s',mwd{subnum})); cd(mwd{subnum}); end

   if getmriflag
     hdr = spm_select('FPList',mridata{subnum},'^*.dcm');
     hdr = spm_dicom_headers(hdr);
     
     patsex{subnum} = hdr{1}.PatientsSex;
     patage{subnum} = hdr{1}.PatientsAge;   % NB At time of scanning, not time of MEG!
     
     spm_dicom_convert(hdr,'all','flat','nii');

     mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-0[0-9].nii$');
     % (filter setting just to ensure picks up raw structural - assuming CBU onvention)
     mrifile{sub} = strtrim(mrifile{sub}(1,:));   % Take first if multiple previous

     disp('!!!GOOD IDEA TO MANUALLY REPOSITION STRUCTURAL VIA DISPLAY BUTTON!!!')
%     pause
     
     if VitECapsules(subnum)
      meg_headmask(mrifile{sub},'display');         % Remove any VitE capsules
      disp('!!!Check difference image!!!');
%      pause     
      eval(sprintf('!mv %s_masked.nii %s',mrifile{sub}(1:end-4),mrifile{sub}))
     end
     patage
     patsex
   
   else
     mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-0[0-9].nii$');
     mrifile{sub} = strtrim(mrifile{sub}(1,:));   % Take first if multiple previous     
     wmrifile{sub} = spm_select('FPList',mwd{subnum},'^wms.*-0[0-9].nii$');
     wmrifile{sub} = strtrim(wmrifile{sub}(1,:));   % Take first if multiple previous     
   end
end
spm_check_registration(strvcat(mrifile));             % Only feasible for <~20 images
try, spm_check_registration(strvcat(wmrifile)); end   % (sub5 bit too warped at back)

disp('!!!NOW DISPLAY MRI IMAGE in SPM AND DETERMINE FIDUCIALS!!!')

%%%%%%%%%% Coordinates of fiducials from above MRIs (nasion, left, right)
% Below are from using MRI as is (without manually reorienting) - very approximate
% For information about manual repositioning, see 
%                         http://imaging.mrc-cbu.cam.ac.uk/meg/RepositioningMRIs
smri_fids{1}  = [ 3 111  22; -78 18 -31; 77 11 -33];
smri_fids{2}  = [ 6 107 -28; -69 26 -41; 68 14 -45];
smri_fids{3}  = [-2 124 -23; -78 22 -52; 70 21 -58];
smri_fids{4}  = [ 3 109 -8;  -72 18 -41; 72 18 -47];
smri_fids{5}  = [-15 103 -33; -79 3 -38; 62 13 -46];
smri_fids{6}  = [-4 109 -20; -71 15 -32; 63 11 -35];
smri_fids{7}  = [-1 93  -5;  -72 14 -44; 70 10 -44];
smri_fids{8}  = [ 1 99  -22; -73  8 -27; 69 10 -42];
smri_fids{9}  = [-3 125 -24; -76 22 -55; 69 22 -58];
% Below are after my manual reorienting (just for your reference)
%smri_fids{1}  = [ 0.7 80.0 -31.2;   -76.9 -14.3 -45.9;   79.9 -19.2 -41.0];
%smri_fids{2}  = [ 1.3 78.7 -30.9;   -71.8  -2.5 -58.0;   71.5   1.8 -58.0];  
%smri_fids{3}  = [ 1.7 88.8 -24.6;   -74.0  -1.4 -66.1;   71.7  -2.8 -66.1];
%smri_fids{4}  = [-1.6 78.9 -22.4;   -74.8 -13.3 -57.7;   71.6  -7.2 -54.6];
%smri_fids{5} = {TBA};
%smri_fids{6}  = [-2.4 79.1 -27.1;   -72.1  -9.8 -44.6;   64.3  -9.8 -51.9];
%smri_fids{7}  = [-2.3 73.2  -7.0;   -71.8  -2.0 -46.9;   70.2  -4.9 -49.8]; 
%smri_fids{8}  = [ 0.0 78.1 -30.3;   -71.3  -1.5 -48.5;   72.8   6.0 -51.5]; 
%smri_fids{9}  = [ 0.0 88.9 -31.1;   -76.8 -15.5 -56.6;   78.0 -12.7 -55.0];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Localise on Mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addpath /imaging/local/Brainstorm/Toolbox

% !!! There are several options for source localisation, eg choice of
% forward model (eg Spheres or BEMs), and choice of inversion method (eg
% MSP, COH, MNM). Below illustrates just two possible inversions (each 
% inversion is indexed by "val"), the second of which has been commented 
% out because BEMs take a VERY LONG time (and don't work well for EEG yet). 
% Of course many more options can be tried (and compared using the model-evidence,
% recorded in mod_evd below, provided that you don't also change the
% data, eg epoch inverted). 
%
% For more information about forward models, see 
%    http://imaging.mrc-cbu.cam.ac.uk/meg/SpmForwardModels
% For more information about inversion parameters, see 
%    Friston et al (2008) Neuroimage 
%    (PDF available here http://imaging.mrc-cbu.cam.ac.uk/meg/SpmAnalysis)

%%%%%%%%%%%%%%%%%%%%
val = 1;  % Spherical forward model using canonical cortical mesh but individual head meshes
comment{val} = 'Concentric Spheres';
cortex{val}  = 'Can';   % Canonical (inverse-normalised template) cortical mesh
CtxSize(val) = 8196;    % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
iskull{val}  = 'Ind';   % Individual (SPM-created) inner skull mesh
oskull{val}  = '';      % (Outer skull irrelevant for spherical models)
scalp{val}   = 'Ind';   % Individual (SPM-created) scalp mesh
formod{val}  = 'Sph';   % Concentric spheres
surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all sensor-types

invcons{val} = [1 2];   % Conditions to invert
invtype{val} = 'GS';    % Type of inversion (This is MSP, with a Greedy Search (GS))
invtwin{val} = expwin;  % Invert whole epoch

%%%%% Options for subsequent time-freq contrast
contwin{val}{1} = [140 190];  % Time window for contrast
contwin{val}{2} = [];         % Frequency window for contrast ([]=all freqs) 
conengy{val}    = 'evoked';   % or 'induced', but spm_eeg_invert_fuse doesnt handle yet

%%%%%%%%%%%%%%%%%%%%
val = 2;  % As val=1 except IID (ie standard minimum norm) used to invert
comment{val} = 'Concentric Spheres';
cortex{val}  = 'Can';   % Canonical (inverse-normalised template) cortical mesh
CtxSize(val) = 8196;    % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
iskull{val}  = 'Ind';   % Individual (SPM-created) inner skull mesh
oskull{val}  = '';      % (Outer skull irrelevant for spherical models)
scalp{val}   = 'Ind';   % Individual (SPM-created) scalp mesh
formod{val}  = 'Sph';   % Concentric spheres
surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all data-types

invcons{val} = [1 2];   % Conditions to invert
invtype{val} = 'IID';   % Type of inversion
invtwin{val} = expwin;  % Invert whole epoch

%%%%% Options for subsequent time-freq contrast
contwin{val}{1} = [140 190];  % Time window for contrast
contwin{val}{2} = [];         % Frequency window for contrast ([]=all freqs) 
conengy{val}    = 'evoked';   % or 'induced', but spm_eeg_invert_fuse doesnt handle yet

%%%%%%%%%%%%%%%%%%%%
%val = 2;  % BEM using canonical meshes (Warning: BEMs take a long time to compute!)
%
%comment{val} = 'Boundary Element Model';
%cortex{val}  = 'Can';    % Canonical (inverse-normalised template) cortical mesh
%CtxSize(val) = 8196;     % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
%iskull{val}  = 'Can';    % Canonical (inverse-normalised template) inner skull mesh
%oskull{val}  = 'Can';    % Canonical (inverse-normalised template) outer skull mesh
%scalp{val}   = 'Can';    % Canonical (inverse-normalised template) scalp mesh
%formod{val}  = 'BEM';    % Boundary Element Model
%surfit(val,:)= [2 2 4];  % Surfaces for BEM: up to 2nd for MEG (iskull) and 4th for EEG (+oskull+scalp)
%invcons{val} = [1 2];    % Conditions to invert
%invtype{val} = 'GS';     % Type of inversion (This is MSP, with a Greedy Search (GS))
%invtwin{val} = expwin;   % Invert whole epoch

%%%%% Options for subsequent time-freq contrast
%contwin{val}{1} = [140 190];  % Time window for contrast
%contwin{val}{2} = [];         % Frequency window for contrast ([]=all freqs) 
%conengy{val}    = 'evoked';   % or 'induced', but spm_eeg_invert_fuse doesnt handle yet

% Other options you can consider are "overlapping spheres" for the
% forward model (see spm_eeg_inv_BSTfwdsol.m), minimumum norm (IID) for
% the inversion, etc. See Henson et al (2009), Neuroimage, for more options.

Nval = length(comment);
mod_evd=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If want to pass prior estimates of sensor-level error (noise), eg
% estimates from empty-room data in case of MEG(see Henson et al, 2009b)...
tmp = load('/imaging/rh01/Methods/MEG_Group/MeanCn.mat');
Cn{1} = tmp.Cn{1};
Cn{2} = tmp.Cn{2}; 
%Cn{1} = tmp.Cn{1}-eye(size(tmp.Cn{1}))*mean(diag(tmp.Cn{1}));  % Orthogonalise
%Cn{2} = tmp.Cn{2}-eye(size(tmp.Cn{2}))*mean(diag(tmp.Cn{2}));  % Orthogonalise
% !!!Care: Cn only relevant to SSSt and 40Hz lowpass (see Henson et al, 2009b)!!!

for sub = 1:Nsub

 cd(wd)
 subnum = dosubs(sub);
 swd = sprintf('sub%d',subnum);
 cd(swd)
 
% If want to pass prior estimates of sensor-level error (noise)...
% !!!Care - sensor order in sennam must match order in Cn, ie Cn{1}=mags, Cn{2}=grads!!!
 Ce{1} = {speye(Nsubchan(sub,1)) Cn{1}};   % Mags (Cn assumes no bad channels!!!!)
 Ce{2} = {speye(Nsubchan(sub,2)) Cn{2}};   % Grds (Cn assumes no bad channels!!!!)
 Ce{3} = {speye(Nsubchan(sub,3))};         % EEG  (no noise estimates yet)
% Otherwise, do not pass anything, corresponding to just IID ('white') noise only:
% Ce = cell(1,Nsen);
   
 for val = 1:Nval
   
  for sen = find(invertflag) % Nsen

   clear S;  
   S.D = char(preavgnam{sen,sub});
% This is using all trials to estimate sources, ie total (induced+evoked) power.
% This gives more sample points in total, and the pure evoked part can always be 
% extracted when estimating time-freq contrasts afterwards (see below). If you want
% to localise purely evoked power, then use instead:
%   S.D = char(avgnam{sen,sub});
% Also may prefer "mvcomp" rather than "trans" data, if don't trust "trans default"

   D = spm_eeg_ldata(S.D); disp(S.D)

%If want to clear previous localisations (safer?)
   if val==1
     if isfield(D,'inv');  D = rmfield(D,'inv'); save(D.fname,'D'); end
     if isfield(D,'con');  D = rmfield(D,'con'); save(D.fname,'D'); end
   end
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise

   disp(comment{val})
   D.val = val; D.fname
   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 = comment{val};

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

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

   D.inv{val}.mesh.sMRI = mrifile{sub};   % In case many from previous

   [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);
   mri_stem = nam;
   wmrifile{sub} = fullfile(pth,['wm' nam ext]);
   if ~exist(wmrifile{sub})
    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

%%%%%%%%%%%%%%%%%%% Meshing 2: Creating binary images for scalp (and inner skull)

   D.inv{val}.mesh.msk_scalp = spm_select('List',mwd{subnum},'^*scalp\.nii$');
 
   if isempty(D.inv{val}.mesh.msk_scalp)
    D = spm_eeg_inv_getmasks(D);	
   else	% assume other masks exist too!
    D.inv{val}.mesh.msk_scalp   = fullfile(mwd{subnum},D.inv{val}.mesh.msk_scalp);
    D.inv{val}.mesh.msk_iskull  = spm_select('FPList',mwd{subnum},'^*iskull\.nii$');
    D.inv{val}.mesh.msk_cortex  = spm_select('FPList',mwd{subnum},'^*cortex\.nii$');
    D.inv{val}.mesh.msk_flags   = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
   end
   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));

%%%%%%%%%%%%%%%%%%% Meshing 2: Creating meshes

   D.inv{val}.mesh.template  = 0;
   D.inv{val}.mesh.canonical = 1;
   D.inv{val}.mesh.Msize = CtxSize(val);
  
  % Prompt for canonical or user-specified cortical meshs (SPM does not currently
  % produce individual ones, because automatic meshing of cortex is difficult!)
 
   if strcmp(cortex{val},'Can')
    Tmesh = load(fullfile(spm('dir'),'EEGtemplates',sprintf('mni_mesh_cortex_%d.mat',CtxSize(val))));
    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);
   elseif strcmp(cortex{val},'Load')
    Tmesh = spm_select(1,'mat','Select cortex mat file containing vert and face',[],pwd,'.*mat');
    D.inv{val}.mesh.tess_ctx.vert    = Tmesh.vert;
    D.inv{val}.mesh.tess_ctx.face    = uint16(Tmesh.face);
    Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.invdef);
    D.inv{val}.mesh.tess_mni.vert    = Tmesh.vert;
    D.inv{val}.mesh.tess_mni.face    = uint16(Tmesh.face);
   else
    error('No cortical mesh specified????')
   end

   iskull_mesh_name = fullfile('MRI',[mri_stem '_iskull_mesh.mat']);
   scalp_mesh_name = fullfile('MRI',[mri_stem '_scalp_mesh.mat']);
    
  % Prompt for canonical, individual (from SPM) or user-specified inner skull
  % (Use canonical if worried about canonical cortex piercing innner skull)

   if strcmp(iskull{val},'Can')
    Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_iskull_2562.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);
   elseif strcmp(iskull{val},'Ind')
    if ~exist(iskull_mesh_name)
     D = spm_eeg_inv_getmeshes(D,1);
     vert = D.inv{val}.mesh.tess_iskull.vert;
     face = D.inv{val}.mesh.tess_iskull.face;
     save(iskull_mesh_name,'vert','face');
    else
     Tmesh = load(iskull_mesh_name);
     D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
     D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
    end
   elseif strcmp(iskull{val},'Load')
     Tmesh = spm_select(1,'mat','Select iskull mat file containing vert and face',[],pwd,'.*mat');
     D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
     D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
   else
     warning('No inner skull mesh specified')
   end

  % Prompt for canonical or user-specified outer skull (SPM does not currently
  % produce individual ones, because segmenting outer skull from scalp difficult)
  % (Note that outer skull only really necessary for EEG BEMs)

   if strcmp(oskull{val},'Can')
    Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_oskull_2562.mat'));
    Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
    D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
    D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
   elseif strcmp(oskull{val},'Load')
     Tmesh = spm_select(1,'mat','Selectoskull  mat file containing vert and face',[],pwd,'.*mat');
     D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
     D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
   else
     warning('No outer skull mesh specified')
   end

  % Prompt for canonical, individual (from SPM) or user-specified scalp
  % (Use canonical if worried about canonical cortex or inner skull piercing scalp)

   if strcmp(scalp{val},'Can')
    Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_scalp_2562.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);
   elseif strcmp(scalp{val},'Ind')
    if ~exist(scalp_mesh_name)
     D = spm_eeg_inv_getmeshes(D,2);
     vert = D.inv{val}.mesh.tess_scalp.vert;
     face = D.inv{val}.mesh.tess_scalp.face;
     save(scalp_mesh_name,'vert','face');
    else
     Tmesh = load(scalp_mesh_name);
     D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
     D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
    end
   elseif strcmp(scalp{val},'Load')
     Tmesh = spm_select(1,'mat','Select scalp mat file containing vert and face',[],pwd,'.*mat');
     D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
     D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
   else
     warning('No scalp mesh specified')
   end

  % Always individual scalp for datareg
  
   if ~exist(scalp_mesh_name)
     S = spm_eeg_inv_getmeshes(D,2);   % lastarg==1 returns iskull (not scalp)
     vert = S.inv{val}.mesh.tess_scalp.vert;
     face = S.inv{val}.mesh.tess_scalp.face;
     save(scalp_mesh_name,'vert','face');
     D.inv{val}.mesh.tess_scalp_datareg.vert  = vert;
     D.inv{val}.mesh.tess_scalp_datareg.face  = uint16(face);  
   else
     Tmesh = load(scalp_mesh_name);
     D.inv{val}.mesh.tess_scalp_datareg.vert  = Tmesh.vert;
     D.inv{val}.mesh.tess_scalp_datareg.face  = uint16(Tmesh.face);  
   end
  
   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.fid_eeg    = D.channels.fid_eeg;
    D.inv{val}.datareg.headshape  = D.channels.headshape;
    if ~eegflag(sen)
     D.inv{val}.datareg.megorient = D.channels.Orient';
    end
    D.inv{val}.datareg.scalpvert  = D.inv{val}.mesh.tess_scalp_datareg.vert;
    D.inv{val}.datareg.fid_mri    = smri_fids{subnum};
    D = spm_eeg_inv_datareg(D);
    fid_err(sub) = sqrt(mean(mean((D.inv{val}.datareg.fid_coreg - D.inv{val}.datareg.fid_mri).^2)))
   end
  
   spm_eeg_inv_checkdatareg(D);
   save(D.fname,'D');
   clear S; S.D=D.fname; 
%   meg_pretty_datareg(S);        % Prettier rendering of head and sensors

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

   rotate3d off
   D.inv{val}.method = 'Imaging';
   if ~eegflag(sen)
    if strcmp(formod{val},'Sph'), meth = 'meg_sphere';
    elseif strcmp(formod{val},'BEM'), meth = 'meg_bem'; 
    elseif strcmp(formod{val},'OS'), meth = 'meg_os'; 
    else, error('Unknown forward model!'); end
   else
    if strcmp(formod{val},'Sph'), meth = 'eeg_3sphereBerg';
    elseif strcmp(formod{val},'BEM'), meth = 'eeg_bem'; 
    elseif strcmp(formod{val},'OS'), meth = 'eeg_os'; 
    else, error('Unknown forward model!'); end
   end
   D.inv{val}.forward.method =  meth;
   D.inv{val}.forward.surface2fit = surfit(val,sen);
%   D.inv{val}.forward.NVertMax = 1000;   % If want to save some time!
   gain_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatrix_' upper(sennam{sen}) ...
 		    '_' formod{val} '_' num2str(surfit(val,sen)) '.mat'])
   gxyz_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatxyz_' upper(sennam{sen}) ...
		    '_' formod{val} '_' num2str(surfit(val,sen)) '.mat'])
  
%   delete(gain_name);   % if want to be safe (uncomment if want to save time)
   if ~exist(gain_name) %| ~exist(gxyz_name)
    [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
    eval(sprintf('!mv %s %s',D.inv{val}.forward.gainmat,gain_name))
    eval(sprintf('!mv %s %s',D.inv{val}.forward.gainxyz,gxyz_name))
   end

   D.inv{val}.forward.gainmat = gain_name;
   D.inv{val}.forward.gainxyz = gxyz_name;
   save(D.fname,'D');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert each modality separately

   D.inv{val}.inverse = [];
   D.inv{val}.inverse.trials = invcons{val};
   D.inv{val}.inverse.type   = invtype{val};
   D.inv{val}.inverse.woi    = invtwin{val};
   D.inv{val}.inverse.lpf    = 0;   % Confusing, but hpf refers to cutoff
   D.inv{val}.inverse.hpf    = 40;  % of lowpass (and lpf to cutoff of highpass)
% Some extra parameters that you could play with...(but defaults seem ok!)
%   D.inv{val}.inverse.smooth = 0.8;
%   D.inv{val}.inverse.Np = 256;
%   D.inv{val}.inverse.Han = 0;
%   D.inv{val}.inverse.Nm = [];
%   D.inv{val}.inverse.NmTOL = -Inf;    
%   D.inv{val}.inverse.NrTOL = -1;
   D.inv{val}.inverse.pQe = Ce{sen};
   D.inv{val}.inverse.pQp = [];

   D = spm_eeg_invert(D);
%   D = spm_eeg_invert_fuse(D);   % If want to compare exactly with fused below
   save(D.fname,'D');
%   mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end);  % if want to include model complexity
   mod_evd(sub,val,sen) = D.inv{val}.inverse.F(end);
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast

   D.inv{val}.contrast.woi  = contwin{val}{1};
   D.inv{val}.contrast.fboi = contwin{val}{2};
   D.inv{val}.contrast.type = conengy{val};
   D.val = val;
   D = spm_eeg_inv_results(D); 
   spm_eeg_inv_results_display(D); 

   if write3Dflag
    D.inv{val}.contrast.smooth = source_smooth_spm;
    D.inv{val}.contrast.rms = 1;
%    D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
    D = spm_eeg_inv_Mesh2Voxels(D);
    SourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
   end  

%  If want to save on filesize
%   if val>2
%    canvert = D.inv{val-1}.mesh.tess_mni.vert;
%    canface = D.inv{val-1}.mesh.tess_mni.face;
%    D.inv{val-1}.mesh = [];
%    D.inv{val-1}.mesh.tess_mni.vert = canvert;
%    D.inv{val-1}.mesh.tess_mni.face = canface;
%    D.inv{val-1}.datareg = [];
%    clear canvert canface
%   end 
 
   disp(sprintf('Sub=%d, %s, val=%d: %s',sub,sennam{sen},val,D.inv{val}.comment))
   save(D.fname,'D');
   
   invertedflag(sen) = 1;
   
  end % end sen

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fusion

  if fuseinvertflag
   DD = {};
   for sen = find(invertflag)
     DD{sen} = spm_eeg_ldata(preavgnam{sen,sub});
     DD{sen}.val = val;
     DD{sen}.inv{val}.inverse.pQe = Ce{sen};
     DD{sen}.inv{val}.inverse.pQp = [];
   end
   DD{1}.inv{val}.inverse.Nm = [];   % Reset any parameters from last inversion
   DD{1}.inv{val}.inverse.Nr = [];   % Reset any parameters from last inversion
%   DD{1}.inv{val}.inverse.NrTOL = -1;

   DD{1}.inv{val}.inverse.OutExt = '_fused';
   D = spm_eeg_invert_fuse(DD);  
%   mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
   mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F(end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast

   D.inv{val}.contrast.woi  = contwin{val}{1};
   D.inv{val}.contrast.fboi = contwin{val}{2};
   D.inv{val}.contrast.type = conengy{val};
   D.val = val;
   D = spm_eeg_inv_results(D); 
   spm_eeg_inv_results_display(D); 

   if write3Dflag
    D.inv{val}.contrast.smooth = source_smooth_spm;
    D.inv{val}.contrast.rms = 1;
%    D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
    D = spm_eeg_inv_Mesh2Voxels(D);
    SourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
   end

   sennam{Nsen+1} = 'fused';
   preavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
   invertedflag(Nsen+1) = 1;
  end
  
 end % end val
 
end % end subs

% Review model-evidences (note cannot compare fused to individual modality
% inversions, since data differs)
for sub = 1:Nsub
   subnum = dosubs(sub);
   disp(sprintf('\tSub %d (%d)...',subnum,sub))
   disp(sprintf('\t\tFiducial error (RMS mm): %3.2f',fid_err(sub)))
   for val=1:Nval
     disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(mod_evd(sub,val,:))))))
   end
end

eval(sprintf('save preproc%d.mat subsphere nrejects Nsubchan nbadchan nevents allsssbad dosubs preavgnam sennam invertflag SourceImgs Nval mod_evd fid_err Ce',length(dosubs)))

% Fid errors may be quite large (~8mm) for some subjects, simply because
% I did not bother to mark them carefully on the MRI. Note that the
% actual coregistration is determined primarily by the digitised headshape;
% the fiducials are only really used to check the quality of registration.

%spm_check_registration(strvcat(wmrifile))   % double-check normaliation worked

disp('!!!Review single-subject localisations if you wish!!!')

% At this point, you can review a single-subject's localisation by
% pressing the "3D source reconstruction" button, and then loading one
% of the ae* files (for one modality, or the "*_fused.mat" for the
% fusion of the modalities). You can then press the "mip" button to see
% the overall reconstruction, or the "display" button under "window" to
% see the power in the time-freq contrast (for the selected
% condition). You can also display a movie and change the start-stop
% times, etc. Don't press the other buttons for inversion, fusion, group
% inversion etc because they will be covered below...

%%%%%%%%%%%%%%%%% Or a quick look at mean over subjects...

clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
S.fig = 2;
for val=1:Nval
 S.val = val; 
 for sen=find(invertedflag)
  S.fig = S.fig+1;
  S.P = strvcat(preavgnam{sen,:});
  meg_inv_display_con(S);
 end
end
lastfig = S.fig;

% You should see bilateral fusiform, maximal anteriorly in MEG but more
% posteriorly in EEG (more occipital), plus a mixture after fusion
% Note that this is just the mean across subjects, and stats (below) differ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 3D  SPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stdir = fullfile(wd,'SourceSPMs');
try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end

clear grpsubs imgfiles grpcons
grpsubs{1} = [1:Nsub];
% If want to split subjects, eg into two groups, eg:
% grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];

grpcons = [1:2];
Ngrp = length(grpsubs);
for val=1:Nval
 valdir = fullfile(stdir,sprintf('Inv%d',val));
 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
 for sen = find(invertedflag)
  outdir = fullfile(valdir,sennam{sen})
  for grp = 1:Ngrp
   for sub = 1:length(grpsubs{grp});
    subnum = grpsubs{grp}(sub);
    imgfiles{grp}{sub} = SourceImgs{sen}{val}{subnum}(grpcons,:);
   end
  end
  meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
 end
end
cd(stdir)
disp('!!!Review the group 3D stats if you wish!!!')

%return

% Again you need to know a bit about the SPM Results interface, but press
% Results and select the SPM.mat file for the "mags" directory in the
% "Inv1" (MSP(GS)) directory, and enter a new T-contrast [1 -1] to find 
% voxels in which greater source amplitude for faces than scrambled faces
% Choose 0.001 uncorrected and then press "Volume" to get table of stats.
% You should see a swathe of ventral temporal activity, particularly on right, 
% with two clusters that survive correction for multiple comparisons.
% For GRDS, you get smaller, more anterior fusiform clusters; for EEG
% you get much more posterior occipital clusters; for fused results, you
% get a more focal lateral midfusiform. Note that the differences across
% these SPMs are likely more apparent than real, given thresholding and
% the relatively small number of subjects (see Henson et al, 2007,
% Neuroimage and Taylor & Henson, in prep, for more information). MSP
% should also benefit more from group inversion (below).
%
% For the IID inversions (in the "Inv2" directory), the results are more
% consistent across modalities, but also more superficial (eg lateral);
% a well known bias of standard minimum norm.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 3D  PPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This uses SPMs above to calculate Posterior Probability Maps (PPMs),
% which do not have any mutliple comparison issues (eg no need for RFT)
% However, it will take a long time...

%for val=1:Nval
% valdir = fullfile(stdir,sprintf('Inv%d',val));
% cd(valdir); 
% for sen = find(invertedflag)
%  outdir = fullfile(valdir,sennam{sen})
%  cd(outdir)
%  load('SPM.mat');
%  spm_spm_Bayes(SPM);
% end
%end
%disp('!!!Review the group 3D PPMs if you wish!!!')

% Again you need to know a bit about the PPM Results interface, but press
% Results and select an SPM.mat file for one of the modalities and enter
% a new T-contrast [1 -1], but this time, select "Bayesian" rather than
% "classical" when prompted. Then change the default effect size threshold 
% to 0, and probability threshold to 0.99. This will show voxels with a
% >99% probability of the increase for faces relative to scrambled being
% greater than zero.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Group inversion (Litvak & Friston, 2008)
% This pools MSP priors across subjects, in theory to make localisations
% more consistent (it has no effect on IID inversions).

grp_mod_evd=[];
if grpinvertflag
 for sen = find(invertedflag)
  clear DD  
  for sub=1:Nsub
   DD{sub} = spm_eeg_ldata(preavgnam{sen,sub});
  end
  for val=1:Nval   % NB: Group inversion makes no difference for IID (val=2 here)
   for sub=1:Nsub
    Ce{3} = {speye(Nsubchan(sub,3))};         % EEG  (no noise estimates yet)
    DD{sub}.val = val;
    DD{sub}.inv{val}.comment    = [DD{sub}.inv{val}.comment 'Group'];
    DD{sub}.inv{val}.inverse.Nm = [];   % Reset any parameters from last inversion
    DD{sub}.inv{val}.inverse.Nr = [];   % Reset any parameters from last inversion
%    DD{sub}.inv{val}.inverse.NrTOL = -1;
    DD{sub}.inv{val}.inverse.pQe = Ce{sen};
    DD{sub}.inv{val}.inverse.pQp = [];
   end
   DD = spm_eeg_invert(DD);

   % Save outputs
   for sub=1:Nsub
    D       = DD{sub};
    D.fname = [D.fname(1:(end-4)) '_grp.mat'];
    grppreavgnam{sen,sub} = fullfile(D.path, D.fname);
    save(fullfile(D.path, D.fname), 'D');
    grp_mod_evd(sub,val,sen) = D.inv{val}.inverse.F(end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast

    D.inv{val}.contrast.woi  = contwin{val}{1};
    D.inv{val}.contrast.fboi = contwin{val}{2};
    D.inv{val}.contrast.type = conengy{val};
    D.val = val;
    D = spm_eeg_inv_results(D); 
    spm_eeg_inv_results_display(D); 

    if write3Dflag
     D.inv{val}.contrast.smooth = source_smooth_spm;
     D.inv{val}.contrast.rms = 1;
%     D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
     D = spm_eeg_inv_Mesh2Voxels(D);
     GrpSourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
    end 
    save(fullfile(D.path, D.fname), 'D');
   end % end sub
  end % end val
 end % end sen
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fuse any Group inversion
 if fuseinvertflag
  for sub=1:Nsub
   Ce{3} = {speye(Nsubchan(sub,3))};         % EEG  (no noise estimates yet)
   clear DD;
   for sen = find(invertflag)
     DD{sen} = spm_eeg_ldata(grppreavgnam{sen,sub});
   end

   for val=1:Nval   % NB: Group inversion makes no difference for IID (val=2 here)
    for sen = find(invertflag)
     DD{sen}.val = val;
     DD{sen}.inv{val}.inverse.Nm = [];   % Reset any parameters from last inversion
     DD{sen}.inv{val}.inverse.Nr = [];   % Reset any parameters from last inversion
%     DD{sen}.inv{val}.inverse.NrTOL = -1;
     DD{sen}.inv{val}.inverse.pQe = Ce{sen};
     DD{sen}.inv{val}.inverse.OutExt = '_fused';
     DD{sen}.inv{val}.inverse.pQp = {DD{sen}.inv{val}.inverse.QP};   % Use group source priors...
    end
    DD{1}.inv{val}.inverse.grpopt = 0;  % ...so no need to (re)optimise source priors
    
    D = spm_eeg_invert_fuse(DD);  
    grp_mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F(end);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast

    D.inv{val}.contrast.woi  = contwin{val}{1};
    D.inv{val}.contrast.fboi = contwin{val}{2};
    D.inv{val}.contrast.type = conengy{val};
    D.val = val;
    D = spm_eeg_inv_results(D); 
    spm_eeg_inv_results_display(D); 

    if write3Dflag
     D.inv{val}.contrast.smooth = source_smooth_spm;
     D.inv{val}.contrast.rms = 1;
%     D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
     D = spm_eeg_inv_Mesh2Voxels(D);
     GrpSourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
    end
    save(fullfile(D.path, D.fname), 'D');
   end  % end val
   
   grppreavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
  end   % end sub
  sennam{Nsen+1} = 'fused';
 end
end

% Review model-evidences (note cannot compare fused to individual modality
% inversions, since data differs)
for sub = 1:Nsub
   subnum = dosubs(sub);
   disp(sprintf('\tSub %d (%d)...',subnum,sub))
   for val=1:Nval
     disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(grp_mod_evd(sub,val,:))))))
   end
end

eval(sprintf('save preproc%d.mat subsphere nrejects Nsubchan nbadchan nevents allsssbad dosubs preavgnam sennam invertflag SourceImgs Nval mod_evd fid_err Ce grp_mod_evd GrpSourceImgs grppreavgnam',length(dosubs)))

%%%%%%%%%%%%%%%%% Quick look at mean over subjects...

clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
S.fig = lastfig;
for val=1:Nval   % NB: Group inversion makes no difference for IID (val=2 here)
 S.val = val;
 for sen = find(invertedflag)
  S.fig = S.fig+1;
  S.P = strvcat(grppreavgnam{sen,:});
  meg_inv_display_con(S);
  drawnow
 end
end
lastfig = S.fig;

% The group MSP results look similar to before, but with higher maximal values.
% The IID results are (as expected) unchanged by group inversion, apart
% from when the group inversions are fused (when there is a small change,
% as expected owing to the third inversion step).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 3D  SPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stdir = fullfile(wd,'GrpSourceSPMs');
try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end

clear grpsubs imgfiles grpcons
grpsubs{1} = [1:Nsub];
% If want to split subjects, eg into two groups, eg:
% grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
grpcons = [1:2];
Ngrp = length(grpsubs);
for val=1:Nval   % NB: Group inversion makes no difference for IID (val=2 here)
 valdir = fullfile(stdir,sprintf('Inv%d',val));
 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
 for sen = find(invertedflag)
  outdir = fullfile(valdir,sennam{sen})
  for grp = 1:Ngrp
   for sub = 1:length(grpsubs{grp});
    subnum = grpsubs{grp}(sub);
    imgfiles{grp}{sub} = GrpSourceImgs{sen}{val}{subnum}(grpcons,:);
   end
  end
  meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
 end
end
cd(stdir)
disp('!!!Review the group 3D stats if you wish!!!')

%return

% Again you need to know a bit about the SPM Results interface, but press
% Results and select an SPM.mat file for one of the modalities under MSP
% (GS) and enter a new T-contrast [1 -1] to find voxels in which greater source
% amplitude for faces than scrambled faces at 0.001 uncorrected. Results
% are similar, though a bit more reliable, after group inversion than
% previously (though additional orbitofrontal in MAGS and GRDS unclear!).
% The IID results will not have changed (except for the fused results).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 3D  PPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This uses SPMs above to calculate Posterior Probability Maps (PPMs),
% which do not have any mutliple comparison issues (eg no need for RFT)
% However, it will take a long time...

%for val=1:Nval   % NB: Group inversion makes no difference for IID (val=2 here)
% valdir = fullfile(stdir,sprintf('Inv%d',val));
% cd(valdir); 
% for sen = find(invertedflag)
%  outdir = fullfile(valdir,sennam{sen})
%  cd(outdir)
%  load('SPM.mat');
%  spm_spm_Bayes(SPM);
% end
%end
%disp('!!!Review the group 3D PPMs if you wish!!!')

% Again you need to know a bit about the PPM Results interface, but press
% Results and select an SPM.mat file for one of the modalities and enter
% a new T-contrast [1 -1], but this time, select "Bayesian" rather than
% "classical" when prompted. Then change the default effect size threshold 
% to 0, and probability threshold to 0.99. This will show voxels with a
% >99% probability of the increase for faces relative to scrambled being
% greater than zero.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% adding fMRI priors (Henson et al, submitted)

fmri_grp_mod_evd=[];
if fmriinvertflag

% Create fMRI priors from an SPM image....
% Here we assume a group prior in MNI space, but priors could be
% defined separately for each subject based on their own fMRI data

 clear S
 S.D = grppreavgnam{1,1};    % Take first subject/modality; doesnt really matter
 S.fmri = '/imaging/rh01/Methods/fMRIPriors/fMRI_ANOVA_flip0/UFvsS_05_cor_15.img';
 S.gm = '/imaging/local/spm/spm5/canonical/c1single_subj_T1.nii';
 S.space = 1;  % SPM and GM images are in MNI space
 S.hthr = 0.5;  S.ethr = 0;
 S.ncomp = Inf; S.smooth = 0.2;
 S.pflag = 1;

 D0 = spm_eeg_inv_fmripriors(S);
 fMRI = D0.inv{D0.val}.fmri.priors;   % Priors (cell of vectors)
 
% Note that here, we simply add the fMRI priors to the previous
% group-optimised source priors. One could also invert each subject
% afresh, re-optimising the MSP priors at the same time as the fMRI
% priors (either individually, or as a group), by using the preavgnam 
% files (not grppreavgnam ones), and not setting grpopt=0.

 for sen = find(invertflag)
  clear DD
  for sub=1:Nsub
%   DD{sub} = spm_eeg_ldata(preavgnam{sen,sub});
   DD{sub} = spm_eeg_ldata(grppreavgnam{sen,sub});
  end
  for val=1:Nval
   for sub=1:Nsub 
    Ce{3} = {speye(Nsubchan(sub,3))};         % EEG  (no noise estimates yet)
    DD{sub}.inv{val}.comment    = [DD{sub}.inv{val}.comment 'Group+fMRI priors'];
    DD{sub}.inv{val}.inverse.Nm = [];   % Reset any parameters from last inversion
    DD{sub}.inv{val}.inverse.Nr = [];   % Reset any parameters from last inversion
%    DD{sub}.inv{val}.inverse.NrTOL = -1;
    DD{sub}.val = val;
    DD{sub}.inv{val}.inverse.pQe = Ce{sen};
    DD{sub}.inv{val}.inverse.pQp = [];
   end
   DD{1}.inv{val}.inverse.pQp = fMRI;
   DD{1}.inv{val}.inverse.pQp{end+1} = DD{1}.inv{val}.inverse.QP;
%   DD{1}.inv{val}.inverse.grpopt = 1;  % Do (re)optimise source priors
   DD{1}.inv{val}.inverse.grpopt = 0;  % Do not (re)optimise source priors
   
   DD = spm_eeg_invert(DD);

   % Save outputs
   for sub=1:Nsub
    D       = DD{sub};
%    D.fname = [D.fname(1:(end-4)) '_grp_fmri.mat'];
    D.fname = [D.fname(1:(end-4)) '_fmri.mat'];
    fmrigrppreavgnam{sen,sub} = fullfile(D.path, D.fname);
    save(fullfile(D.path, D.fname), 'D');
    fmri_grp_mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast

    D.inv{val}.contrast.woi  = contwin{val}{1};
    D.inv{val}.contrast.fboi = contwin{val}{2};
    D.inv{val}.contrast.type = conengy{val};
    D.val = val;
    D = spm_eeg_inv_results(D); 
    spm_eeg_inv_results_display(D); 

    if write3Dflag
     D.inv{val}.contrast.smooth = source_smooth_spm;
     D.inv{val}.contrast.rms = 1;
%     D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
     D = spm_eeg_inv_Mesh2Voxels(D);
     FmriGrpSourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
    end 
    save(fullfile(D.path, D.fname), 'D');
   end % end val
  end % end sub
 end % end sen
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fuse using Group and fMRI priors
 if fuseinvertflag
  for sub=1:Nsub
   Ce{3} = {speye(Nsubchan(sub,3))};         % EEG  (no noise estimates yet)
   clear DD;
   for sen = find(invertflag)
    DD{sen} = spm_eeg_ldata(fmrigrppreavgnam{sen,sub});
   end
   
   for val=1:Nval
    for sen = find(invertflag)
     DD{sen}.val = val;
     DD{sen}.inv{val}.inverse.Nm = [];   % Reset any parameters from last inversion
     DD{sen}.inv{val}.inverse.Nr = [];   % Reset any parameters from last inversion
%     DD{sen}.inv{val}.inverse.NrTOL = -1;
     DD{sen}.inv{val}.inverse.OutExt = '_fused';
     DD{sen}.inv{val}.inverse.pQe = Ce{sen};
     DD{sen}.inv{val}.inverse.pQp = {DD{sen}.inv{val}.inverse.QP};   % Use group source priors...
    end
    DD{1}.inv{val}.inverse.grpopt = 0;  % ...so no need to (re)optimise source priors
    
    D = spm_eeg_invert_fuse(DD);  
    fmri_grp_mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast

    D.inv{val}.contrast.woi  = contwin{val}{1};
    D.inv{val}.contrast.fboi = contwin{val}{2};
    D.inv{val}.contrast.type = conengy{val};
    D.val = val;
    D = spm_eeg_inv_results(D); 
    spm_eeg_inv_results_display(D); 

    if write3Dflag
     D.inv{val}.contrast.smooth = source_smooth_spm;
     D.inv{val}.contrast.rms = 1;
%     D.inv{val}.contrast.scalefactor = 1;     % if want to compare across dif size windows
     D = spm_eeg_inv_Mesh2Voxels(D);
     FmriGrpSourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
    end

    save(fullfile(D.path, D.fname), 'D');
   end  % end val

   sennam{Nsen+1} = 'fused';
   fmrigrppreavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
  end
 end
end


% Review model-evidences (note cannot compare fused to individual modality
% inversions, since data differs)
for sub = 1:Nsub
   subnum = dosubs(sub);
   disp(sprintf('\tSub %d (%d)...',subnum,sub))
   for val=1:Nval
     disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(fmri_grp_mod_evd(sub,val,:))))))
   end
end

%%%%%%%%%%%%%%%%% Quick look at mean over subjects...

clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
S.fig = lastfig;
for val=1:Nval
 S.val = val;
 for sen = find(invertedflag)
  S.fig = S.fig+1;
  S.P = strvcat(fmrigrppreavgnam{sen,:});
  meg_inv_display_con(S);
  drawnow
 end
end
lastfig = S.fig;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 3D  SPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stdir = fullfile(wd,'fMRIGrpSourceSPMs');
try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end

clear grpsubs imgfiles grpcons
grpsubs{1} = [1:Nsub];
% If want to split subjects, eg into two groups, eg:
% grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
grpcons = [1:2];
Ngrp = length(grpsubs);
for val=1:Nval 
 valdir = fullfile(stdir,sprintf('Inv%d',val));
 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
 for sen = find(invertedflag)
  outdir = fullfile(valdir,sennam{sen})
  for grp = 1:Ngrp
   for sub = 1:length(grpsubs{grp});
   subnum = grpsubs{grp}(sub);
   imgfiles{grp}{sub} = FmriGrpSourceImgs{sen}{val}{subnum}(grpcons,:);
   end
  end
  meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
 end
end
cd(stdir)
disp('!!!Review the group+fMRI 3D stats if you wish!!!')

% Again you need to know a bit about the SPM Results interface, but press
% Results and select an SPM.mat file for one of the modalities under MSP
% (GS) and enter a new T-contrast [1 -1] to find voxels in which greater source
% amplitude for faces than scrambled faces at 0.001 uncorrected. MSP results
% are not changed much by the addition of fMRI priors (one possibility
% being that the MSP are already optimal; Henson et al, submitted). The
% fMRI priors have a more noticeable effect on the IID results, moving
% them closer to the fMRI priors.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 3D  PPMs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This uses SPMs above to calculate Posterior Probability Maps (PPMs),
% which do not have any mutliple comparison issues (eg no need for RFT)
% However, it will take a long time...

%for val=1:Nval 
% valdir = fullfile(stdir,sprintf('Inv%d',val));
% cd(valdir); 
% for sen = find(invertedflag)
%  outdir = fullfile(valdir,sennam{sen})
%  cd(outdir)
%  load('SPM.mat');
%  spm_spm_Bayes(SPM);
% end
%end
%disp('!!!Review the group+fMRI 3D PPMs if you wish!!!')

% Again you need to know a bit about the PPM Results interface, but press
% Results and select an SPM.mat file for one of the modalities and enter
% a new T-contrast [1 -1], but this time, select "Bayesian" rather than
% "classical" when prompted. Then change the default effect size threshold 
% to 0, and probability threshold to 0.99. This will show voxels with a
% >99% probability of the increase for faces relative to scrambled being
% greater than zero.

% Phew! You got to here! Well done!