This example script relies on the example script Create mni-aligned grids in individual head_space. But for Neuromag data there are some differences. First make a mni template as is done in the above mentioned example script.
%%batch create_headmodel_neuromag
%==================================================================
%Loading mri of single subject and make a single shell head model
%==================================================================
filename = 'mri.fif';
mri = ft_read_mri(filename);
save mri mri
%The following makes the mri into ctf coordinates.. you have to select the
%fiducials again yourself
cfg = [];
cfg.interactive = 'yes';
mri = ft_volumerealign(cfg,mri);
save mri_volreal mri
%segment the brain into gray, white and csf matter to later make a single
%shell model.
cfg = [];
cfg.template = '/data/apps/spm/spm2/templates/single_subject_T1.mnc';
cfg.coordinates = 'ctf';
cfg.write = 'no';
cfg.name = 'temp';
[segmentedmri] = ft_volumesegment(cfg, mri)
save segmentedmri segmentedmri
%check how it looks; does the segmented mri fit into the mri? Probably not
%because of neuromag coordinates (x and y are swapped) and a bug in volume
%segment.
cfg = [];
test = segmentedmri;
test.avg.pow = test.gray+test.white+test.csf;
load mri_volreal
test.anatomy = mri.anatomy;
cfg.funparameter = avg.pow';
cfg.interactive = no';
ft_sourceplot(cfg,test);
%------NEUROMAG SPECIFIC
% Swap the x and y coordinates back
test.gray = permute(segmentedmri.gray,[3 2 1]);
test.white = permute(segmentedmri.white,[3 2 1]);
test.csf = permute(segmentedmri.csf,[3 2 1]);
test.dim = [256 256 192];
%------END
%and flip the dimensions back
for t = 1:3;
test.gray = flipdim(test.gray,t);
test.csf = flipdim(test.csf,t);
test.white = flipdim(test.white,t);
end
%check if it worked
test.avg.pow = test.gray+test.white+test.csf;
figure;
cfg.interactive = 'no';
ft_sourceplot(cfg,test);
test = rmfield(test,'anatomy');
test = rmfield(test,'avg');
segmentedmri = test;
save flippedsegmentedmri segmentedmri
%make the single_shell headmodel
cfg = [];
hdm = ft_prepare_singleshell(cfg,segmentedmri);
%------NEUROMAG SPECIFIC
%Flip x and y coordinates to fit to the elekcta neuromag MEG helmet!
xcoords = hdm.bnd.pnt(:,2);
ycoords = hdm.bnd.pnt(:,1);
hdm.bnd.pnt(:,1) = xcoords;
hdm.bnd.pnt(:,2) = ycoords;
%------END
%make the leadfield normalised to the mni template
%normalise the mri first
cfg = [];
cfg.template = '/data/apps/spm/spm2/templates/T1.mnc'; %is in MNI coordinates, from templates/spm2
load mri_volreal
cfg.downsample = 2;
cfg.coordinates = 'ctf';
cfg.nonlinear = 'no';
norm = ft_volumenormalise(cfg,mri); % subjects m
%make the leadfield
%Use the final transform matrix of norm to make the grid.pos (and the template_grid)
load Template_brain/template_grid.mat
grid = [];
grid.pos = warp_apply(inv(norm.cfg.final), template_grid.pos, 'homogenous')/10; %in cm
%------NEUROMAG SPECIFIC
%Flip x and y coordinates to fit the elekcta neuromag MEG helmet!
x = grid.pos(:,2);
y = grid.pos(:,1);
grid.pos(:,1) = x;
grid.pos(:,2) = y;
%------END
grid.inside = template_grid.inside;
grid.outside = template_grid.outside;
grid.dim = template_grid.dim;
clear template_grid
% also remember the normalization transformation matrix
grid.transform = norm.cfg.final;
%to see if everything has worked:
%make a figure of the single subject headmodel, sensor positions and grid positions
cfg = [];
cfg.vol = hdm;
cfg.grid = grid;
% cfg.hdmfile = hdmfile{i}; % an alternative to cfg.vol, this uses the single sphere model
cfg.gradfile = 'grad.mat';
figure;
ft_headmodelplot(cfg);
When you have then estimated the sources which happens in NM or CTF space, you have to replace the .pos field of the source or the result of ft_sourcedescriptives with the template_grid.pos to get it back into MNI space because the origins of the two spaces are different. When you then sourceinterpolate to the normalised mri, this should work!
load f cfg = []; cfg.grid = grid; cfg.frequency = 10; cfg.vol = hdm; cfg.gradfile = 'grad.mat'; cfg.projectnoise = 'yes'; cfg.keeptrials = 'no'; cfg.keepfilter = 'yes'; cfg.keepcsd = 'yes'; cfg.keepmom = 'yes'; cfg.lambda = 0.1 * mean(f.powspctrm(:,9),1); cfg.method = 'dics'; cfg.feedback = 'textbar'; source = ft_sourceanalysis(cfg,f); save source source load grid cfg = []; source.dim = grid.dim; sd = ft_sourcedescriptives(cfg,source); %-----MNI SPECIFIC %Because sourceanalysis worked in NM coordinates and the interpolation goes %in MNI coordinates, we have to replace the sd.pos by the %template_grid.pos load Template_brain/template_grid sd.pos = template_grid.pos; %-----END load norm cfg = []; cfg.parameter = 'avg.nai'; cfg.voxelcoord = 'no'; cfg.interpmethod = 'linear'; sdint = ft_sourceinterpolate(cfg, sd, norm); %plotting the results cfg = []; cfg.surfdownsample = 2; cfg.downsample = 2; cfg.funparameter = 'avg.nai'; cfg.anaparameter = 'anatomy'; cfg.funcolormap = 'jet'; cfg.method = 'ortho'; figure; ft_sourceplot(cfg,sdint)