Create MNI-aligned grids in individual head-space

When combining the source-level data of multiple subjects this data is typically first interpolated (e.g. using ft_sourceinterpolate) and then spatially normalized to a template brain (e.g. using ft_volumenormalise). However it is also possible to define the source-reconstruction grid for each individual subject in such a way that all these grids are already aligned in MNI-space. The combination or statistic of source-level data across subjects can then directly be computed within the source-structure without the need to interpolate and normalize each volume.

You first make a template grid based on a template MRI in MNI coordinates. Subsequently, each individual MRI is warped to the template MRI, and the inverse of that warp is applied to the dipole grid. Hereby the individual subjects' grids are not regularly spaced anymore (meaning the distance between 2 grid points can vary) See figure 1. But because of this warping, a specific grid point is located in the same area of the brain in all subjects and the MNI template.

fig 1a: Example: the MNI template brain and the brains of 3 subjects

fig 1b: Example: the MNI grid and the grid of 3 subjects, note that the each grid point points to the same location in all brains, and that the subjects' grids are not regularly spaced

In the figures above, the spatial deformation of the individual subjects' brains relative to the template brain is exemplified. However, other relevant differences are that the CTF coordinate system (in which the grid positions are expressed) is defined relative to the three coils that are placed on the nose and the ears, whereas the MNI/SPM coordinate system is defined relative to the ventricles. There is a small rotation around translation involved to match the CTF and MNI/SPM origin. Another difference is that in CTF the positive x-axis is pointing towards the nose, and in MNI/SPM the positive x-axis is approximately pointing towards the right ear. Furthermore, the CTF positions are expressed in cm and the MNI positions in mm.

The following code shows how to construct such 'normalized' grids. For each individual, the script also constructs a single-shell volume-conductor model that can be used for the computation of the leadfields (i.e. the forward solution) in ft_sourceanalysis or in ft_prepare_leadfield. The following code is using some standard FieldTrip functions, combined with some low-level fieldtrip functions and some standard Matlab. You should add the fieldtrip/private directory to your Matlab path in order for this to work (explained here). And don't forget to read the How to proceed from here part, totally at the bottom of this page.

% define the location of the single subject anatomical MRI files
mrifile = {
  '/home/common/meg_mri/nass_c/nass_c.mri'
  '/home/common/meg_mri/jokisch_d/jokisch_d.mri'
  '/home/common/meg_mri/henrich_e/henrich_e.mri'
  '/home/common/meg_mri/snoeijs_h/snoeijs_h.mri'
  '/home/common/meg_mri/grothe_i/grothe_i.mri'
  '/home/common/meg_mri/udden_j2/udden_j2.mri'
  '/home/common/meg_mri/vink_k/vink_k.mri'
  '/home/common/meg_mri/weber_k/weber_k.mri'
  '/home/common/meg_mri/bisschops_l/bisschops_l.mri'
  '/home/common/meg_mri/brinkmann_m/brinkmann_m.mri'
  '/home/common/meg_mri/moebius_m/moebius_m.mri'
  '/home/common/meg_mri/ploechl_m/ploechl_m.mri'
  '/home/common/meg_mri/schmidt_m/schmidt_m.mri'
  '/home/common/meg_mri/hoogenboom_n/hoogenboom_n.mri'
  '/home/common/meg_mri/hoemke_p/hoemke_p.mri'
  '/home/common/meg_mri/schelbergen_r/schelbergen_r.mri'
  '/home/common/meg_mri/kommattan_s/kommattan_s.mri'
  '/home/common/meg_mri/nispel_s/nispel_s.mri'
  '/home/common/meg_mri/verhoeven_t/verhoeven_t.mri'
  '/home/common/meg_mri/blankers_v_2/blankers_v.mri'
};
 
% the datasets are only used here for plotting the gradiometer positions, 
% and are not required for the construction of the grids
dataset = {
  'CindyNass_double_20070213_20070213_01.ds'
  'DanielJokisch_double_flash_20070118_01.ds'
  'EmanuelHenrich_double_flash_20070122_01.ds'
  'HarmenSnoeije_double_20070118_20070118_01.ds'
  'IrisGrothe_double_20070119_20070119_02.ds'
  'JuanVidal_double_flash_20061102_02.ds'
  'KarinaVink_double_flash_20070108_01.ds'
  'Kerstin_Weber_double_flash_20061205_01.ds'
  'LonnekeBisschops_double_20070108_20070108_03.ds'
  'MareikeSchmidt_double_20070209_20070209_01.ds'
  'MariaBrinkmann_double_20070209_20070209_03.ds'
  'MartinMoebius_double_flash_20070122_01.ds'
  'MichaelPloechl_double_flash_20061208_01.ds'
  'NienkeHoogenboom_double_flash_20061204__01.ds'
  'PaulHoemke_double_flash_20070110_01.ds'
  'RikSchelberg_double_20070122_01.ds'
  'SinaNispel_double_20070209_20070209_02.ds'
  'SumaKommattam_double_20070213_20070213_02.ds'
  'TjerryVerhoeven_double_flash_20070110_01.ds'
  'VivianBlankers_double_flash_20070119_01.ds'
};
 
% the single subject spherical headmodels are only used here for plotting, 
% and are not required for the construction of the grids
hdmfile = {
  '/home/common/meg_mri/nass_c/nass_c.hdm'
  '/home/common/meg_mri/jokisch_d/jokisch_d.hdm'
  '/home/common/meg_mri/henrich_e/henrich_e.hdm'
  '/home/common/meg_mri/snoeijs_h/snoeijs_h.hdm'
  '/home/common/meg_mri/grothe_i/grothe_i.hdm'
  '/home/common/meg_mri/udden_j2/udden_j2.hdm'
  '/home/common/meg_mri/vink_k/vink_k.hdm'
  '/home/common/meg_mri/weber_k/weber_k.hdm'
  '/home/common/meg_mri/bisschops_l/bisschops_l.hdm'
  '/home/common/meg_mri/brinkmann_m/brinkmann_m.hdm'
  '/home/common/meg_mri/moebius_m/moebius_m.hdm'
  '/home/common/meg_mri/ploechl_m/ploechl_m.hdm'
  '/home/common/meg_mri/schmidt_m/schmidt_m.hdm'
  '/home/common/meg_mri/hoogenboom_n/hoogenboom_n.hdm'
  '/home/common/meg_mri/hoemke_p/hoemke_p.HDM'
  '/home/common/meg_mri/schelbergen_r/schelbergen_r.HDM'
  '/home/common/meg_mri/kommattan_s/kommattan_s.hdm'
  '/home/common/meg_mri/nispel_s/nispel_s.hdm'
  '/home/common/meg_mri/verhoeven_t/verhoeven_t.hdm'
  '/home/common/meg_mri/blankers_v_2/blankers_v.hdm'
};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This template is in MNI coordinates , but from a single subject, which is needed for proper segmentation
template = '/home/common/matlab/spm2/canonical/single_subj_T1.mnc'; %template is distributed with spm 
template_mri = ft_read_mri(template);
 
% segment the template brain and construct a volume conduction model (i.e. head model)
cfg = [];
cfg.coordinates = 'spm';
template_seg = ft_volumesegment(cfg, template_mri); %NB is in mm
cfg = [];
template_hdm = ft_prepare_singleshell(cfg, template_seg); %NB transforms to cm!
 
% this dummy gradiometer definition is required for prepare_dipole_grid and headmodelplot
template_grad     = [];
template_grad.pnt = [];
template_grad.ori = [];
template_grad.tra = [];
template_grad.label = {};
 
% construct the dipole grid in the template brain coordinates
% the source units are in cm
% the negative inwardshift means an outward shift of the brain surface for inside/outside detection
cfg = [];
cfg.grid.xgrid = -20:1:20;
cfg.grid.ygrid = -20:1:20;
cfg.grid.zgrid = -10:1:20;
cfg.grid.tight = 'yes'; %new way
cfg.inwardshift = -1.5;
template_grid = prepare_dipole_grid(cfg, template_hdm, template_grad);
 
% make a figure with the template head model and dipole grid
cfg = [];
cfg.vol = template_hdm;
cfg.grad = template_grad;
cfg.plotsensors = 'no';
cfg.inwardshift = -1.5;
figure
ft_headmodelplot(cfg);
triplot(template_grid.pos(template_grid.inside,:), [], [], 'nodes')
 
% mni space is in mm, thus *10, and save
template_grid.pos = template_grid.pos *10 % cm->mm
template_grid.xgrid = template_grid.xgrid *10 % cm->mm
template_grid.ygrid = template_grid.ygrid *10 % cm->mm
template_grid.zgrid = template_grid.zgrid *10 % cm->mm
 
save('template_grid', 'template_grid'); %You need the .xyzgrid .pos and .dim field later

fig 2a: template grid and headmodel, top view

fig2b: template grid and headmodel, right frontal view

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
nsubj = length(mrifile);
for i.html">i=1:nsubj
 
  % read the single subject anatomical MRI
  mri = ft_read_mri(mrifile{i.html">i});
 
  % segment the anatomical MRI
  cfg = [];
  cfg.downsample = 2;
  cfg.coordinates = 'ctf';
  seg = ft_volumesegment(cfg, mri);
 
  % flip dimensions because of bug in volumesegment
  seg.gray=flipdim(seg.gray,1);
  seg.gray=flipdim(seg.gray,2);
  seg.gray=flipdim(seg.gray,3);
  seg.white=flipdim(seg.white,1);
  seg.white=flipdim(seg.white,2);
  seg.white=flipdim(seg.white,3);
  seg.csf=flipdim(seg.csf,1);
  seg.csf=flipdim(seg.csf,2);
  seg.csf=flipdim(seg.csf,3);
 
  % construct volume conductor model (i.e. head model) for each subject
  % this is optional, you can also use another model, e.g. a single sphere model
  cfg = [];
  hdm = ft_prepare_singleshell(cfg, seg);
 
  % normalize this subjects MRI to the template MRI, 
  % the purpose is to obtain the (inverse) transformation matrix
  cfg = [];
  cfg.template = 'T1.mnc' % default is 'T1.mnc' (spm2) or 'T1.nii' (spm8) 
  %is a smoothed template based on multiple subjects, to deal with between subject variation in anatomy. This is needed for proper normalization. Template distributed with spm (in templates folder) 
  cfg.downsample = 2;
  cfg.coordinates = 'ctf';
  cfg.nonlinear = 'no';
  norm = ft_volumenormalise(cfg,mri); % subjects mri to MNI template
 
  % reverse-normalize the template grid to this individual subject
  grid         = [];
  grid.pos     = warp_apply(inv(norm.cfg.final), template_grid.pos, 'homogenous')/10; %in cm
  grid.inside  = template_grid.inside;
  grid.outside = template_grid.outside;
 
  % also remember the normalization transformation matrix
  transform_ctf2spm = norm.cfg.final;
 
  % 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 = dataset{i.html">i};
  figure
  ft_headmodelplot(cfg);
  triplot(grid.pos(grid.inside,:), [], [], 'nodes')

fig 3: template grid in single-subject head coordinates, single-subject headmodel and sensor positions (right frontal view)

 
  [p, f, x] = fileparts(mrifile{i.html">i});
  filename = [f '.mat'];
  save(filename, 'seg', 'hdm', 'grid', 'transform_ctf2spm');
 
end

The above script also saves the affine transformation matrix for spatial normalization of each subject's brain. This matrix can e.g. be used to project ROI positions defined in MNI-space into the individual head-space.

Keep in mind that: The .pos field in the single subjects' grid points to the location (in cm) in the subjects' MRI which is in CTF coordinates. The single subject grid's are not spaced regularly, therefore the .dim and the .xyzgrid fields are not present in the grid.

You can use these single subjects' MNI grids in source analysis. After that you should put .pos, .xgrid, .ygrid, .zgrid, .dim field from the template_grid onto the subjects's source, which is thereby in MNI coordinates. So interpolating and normalising to compare over subjects is not necessary anymore.