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,...