Source Estimation in SPM8

This is a basic script for source estimation of individual subject data in SPM8. It uses as input the averaged data from SPM8 pre-processing, as well as individual MRI data.

% Root directory for EMEG data
bwd = '/MyDataDir/';
data_name = 'Mmacespm8_block1_raw_ssstf_raw.mat'; % File name for results after averaging (input for source estimation)

%% Define your subject data and MRIs
cnt = 0;

cnt = cnt + 1;
subjects{cnt} = {'meg10_1000', '101224'};   % letter strings for subject-specific paths
MRI_dir{cnt} = '/MRIdirectory/thisone.nii'; % directory with this participant's MRI (*.nii)

cnt = cnt + 1;
subjects{cnt} = {'meg10_1001', '101231'};   % letter strings for subject-specific paths
MRI_dir{cnt} = '/MRIdirectory/thisone.nii'; % directory with this participant's MRI (*.nii) 
% Note: you may have to copy your MRI first (e.g. from /imaging/local/structurals/cbu)

nr_subs = length(subjects);
fprintf(1, 'Going to process %d data sets\n', nr_subs);

%% Define Analysis Parameters 

% Names of conditions after averaging
inv_trls = {'cond1';'cond2';'control'};  

% Time (ms) and Freq (Hz) contrast window (0=all freqs)
con_win = {[0 400]; 0};

% Forward model options, for_typ{EEG/MEG}
for_typ{1} = 'EEG BEM'; % EEG
for_typ{2} = 'Single Shell'; % MEG (mags+grads)

%% Processing options
segment_only = 0;    % Flag whether only MRI segmentation (up to, not including, coregistration) shall be done (1: only segment; 0: do it all)
redoreg_flag = 0;    % Flag whether to redo MRI and registration model (eg if want to change fids)
redoformod_flag = 1; % Flag whether to redo forward modelling (eg if only one for_typ above)
display_flag = 1;    % display results on the fly
write3Dflag = 0;     % write SPM volumes of results

%% Define Inversion Options

% Which sensor types to use
inv_mods = {'MEG';'MEGPLANAR';'EEG'};   % Fusion of mulitple modalities

% Which inversion method to use
inv_typ = 'IID';    % Minimum Norm Least-Squares
%inv_typ = 'GS';   % (MSP)

% Cortical surface option
mesh_size = 2;   % [1-3] for coarse-medium-fine

% Whether to use headshape for coregistration (in datareg)
use_headshape = 1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Source localisation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ss = 1:nr_subs,
    
    % Create individual output directory if necessary
    cd(bwd)
    swd = fullfile(bwd,subjects{ss}{1},subjects{ss}{2});
    try
        eval(sprintf('!mkdir %s',swd))
    end
    cd(swd)

    % current subject's data directory
    subj_dir = fullfile(bwd, subjects{ss}{1}, subjects{ss}{2});
    fprintf(1, '\n Current subject directory: %s\n\n', subj_dir);

    % Load result of data averaging
    data_file = fullfile(subj_dir, data_name);
    D = spm_eeg_load( data_file );

    % Initialise... (if want to be safe!)
    D.inv = {struct('mesh', [], 'gainmat', [])};

    current_formod = 1;

    indF = [];
    val = 0;

    %% MRI processing and coregistration
    if redoreg_flag == 1 | ~isfield(D.inv{1},'mesh') | ~isfield(D.inv{1},'datareg')

        val = val + 1;
        D.val = val;
        D.inv{val}.date    = strvcat(date,datestr(now,15));
        D.inv{val}.comment = {sprintf('%s_%s',inv_typ,char(inv_mods)')};    % remember inversion type and sensor configuration    
        D.inv{val}.mesh    = [];
        D.inv{val}.datareg = [];
        D.inv{val}.forward = [];

        %% Check whether MRI present (care if more than one such *img!)        
        D.inv{val}.mesh.sMRI = spm_select('FPList',MRI_dir{ss},'^s.*\.nii$');
        if isempty(D.inv{val}.mesh.sMRI)
            error('MRI not found')
        end

        %% Check whether inverse-normalised surfaces present
        D.inv{val}.mesh.def  = spm_select('FPList',MRI_dir{ss},'^y_sM.*\.nii$');
        if isempty(D.inv{val}.mesh.def)
            warning('No inverse normalisation parameters found - segmenting may take a while!')
        end
        %% Normalise sMRI (if not done already), and create inverse-normalised surfaces
        D.inv{val}.mesh = spm_eeg_inv_mesh(D.inv{val}.mesh.sMRI, mesh_size);
        if display_flag, spm_eeg_inv_checkmeshes(D); end
        % D.save;
            
        % If coregistration not yet done...
        if segment_only,
            continue;   % Stop processing for this subject here, continue with next one
        end;
              
        % If coregistration done manually and saved...
        newmrifid           = [];
        newmrifid.fid.pnt   = D.inv{val}.mesh.fid.fid.pnt(1:3,:);  % I think these are the fids after "Save"...
        newmrifid.fid.label = {'Nasion';'LPA';'RPA'};
        newmrifid.pnt       = D.inv{val}.mesh.fid.pnt;        % Scalp mesh points from MRI above
    
        meegfid = D.fiducials;
        %% Remove nose points (assume those for which y>0 and z<0)
        meegfid.pnt(find(meegfid.pnt(:,2)>0 & meegfid.pnt(:,3)<0),:) = [];
    
        fprintf(1, 'Coregistering\n');
        D = spm_eeg_inv_datareg_ui(D, val, meegfid, newmrifid,use_headshape);
    
        fprintf(1, 'Done Coregistering\n');

        for ind = 1:length(D.inv{val}.datareg)
            fprintf(1, '%d ', ind);
            d = D.inv{val}.datareg(ind).fid_eeg.fid.pnt - D.inv{val}.datareg(ind).fid_mri.fid.pnt;
            err(ss,ind,val) = mean(sqrt(sum(d.^2,2)));
        end 
        fprintf(1, '\n');

        redoreg_flag = 0;

    else                
        D.inv{val}.mesh    = D.inv{1}.mesh;
        D.inv{val}.datareg = D.inv{1}.datareg;                
    end % if redoreg_flag...

    %% Computing forward model/leadfield
    if redoformod_flag == 1 | ~isfield(D.inv{1},'forward') | ~exist(D.inv{1}.gainmat)

        fprintf(1, 'Creating forward model\n');
        %% Create forward model (BEM) (could conditionalise this bit on modality inverted...)

        D.inv{val}.forward = struct([]);
        for ind = 1:length(D.inv{val}.datareg)
            D.inv{val}.forward(ind).voltype = for_typ{ind};
        end

        fprintf(1, 'Computing leadfield\n');
        D = spm_eeg_inv_forward(D);

        if display_flag
             for ind = 1:length(D.inv{val}.datareg)
                 spm_eeg_inv_checkforward(D, val, ind); pause(3);
             end
        end
        current_formod = val;
        redoformod_flag = 0;
        D.save;
    else
        D.inv{val}.forward = D.inv{current_formod}.forward;
        D.inv{val}.gainmat = D.inv{current_formod}.gainmat;
    end     % redoformod_flag = ...


    %% Invert & Contrast

    D.inv{val}.inverse = [];   % Clear to be safe!
    D.inv{val}.inverse.trials = inv_trls;
    D.inv{val}.inverse.type   = inv_typ;
    %D.inv{val}.inverse.lpf    = 0;
    %D.inv{val}.inverse.hpf    = 40;
    %D.inv{val}.inverse.woi    = [-100 300];
    %D.inv{val}.inverse.Han    = 1;

    D.inv{val}.inverse.modality = inv_mods;

    D = spm_eeg_invert(D);

    D.inv{val}.contrast.woi  = con_win{1};
    D.inv{val}.contrast.fboi = con_win{2};
    D.inv{val}.contrast.type = 'evoked';

    D = spm_eeg_inv_results(D);

    D.con=1; spm_eeg_inv_results_display(D);    % Display first condition

    % Write result to SPM volumes
    if write3Dflag
        D.inv{val}.contrast.smooth = 12;
        % D.inv{val}.contrast.rms = 1;
        % D.inv{val}.contrast.scalefactor = 1;
        D = spm_eeg_inv_Mesh2Voxels(D);
        SourceImgs{val}{ss} = strvcat(D.inv{val}.contrast.fname);
    end

    indF(ss,val) = D.inv{val}.inverse.F;

    D.save;

end;    % for ss = 1:nr_subs,...