ENH - added example script to construct MNI-aligned template source model for FEM

thanks to @NataschaRoos and @vitpia

Author: robertoostenveld

example.md

@@ -55,6 +55,7 @@ See also the [tutorials](/tutorial) and [frequently asked questions](/faq).
 - [Compute EEG leadfields using a concentric spheres headmodel](/example/source/concentricspheres)
 - [Compute EEG leadfields using a BEM headmodel](/example/source/bem)
 - [Compute EEG leadfields using a FEM headmodel](/example/source/fem)
+- [Create a template source model aligned to MNI space](/example/source/sourcemodel_mnitemplate)
 - [Compute forward simulated data and apply a beamformer scan](/example/compute_forward_simulated_data_and_apply_a_beamformer_scan)
 - [Compute forward simulated data and apply a dipole fit](/example/compute_forward_simulated_data_and_apply_a_dipole_fit)
 - [Compute forward simulated data using ft_dipolesimulation](/example/compute_forward_simulated_data)

example/source/sourcemodel_mnitemplate.m

@@ -0,0 +1,118 @@
+---
+title: Create a template source model aligned to MNI space
+category: example
+tags: [meg, mri, headmodel, source]
+---
+
+# Create template source models aligned to MNI space
+
+On the [template sourcemodel](/template/sourcemodel) page we describe that we have a number of 3D grid
+source models defined in MNI space. These models can be used as templates for individual subject data,
+allowing for more accurate and consistent source reconstruction across subjects. The procedure for this
+is described in the [source model tutorial](/tutorial/sourcemodel).
+
+These source models are only available at a few resolutions, and have dipoles throughout the whole brain
+compartment, including CSF, white matter and gray matter. That is convenient for beamformer source
+reconstruction of MEG and EEG data with a BEM or single shell head model, but not so for minimum
+norm source reconstruction or for source estimates using a FEM model.
+
+The following script describes how to create your own template source model in MNI space, at 3mm
+resolution, and only with dipoles in gray matter.
+
+{% include markup/yellow %}
+Note that "gray" is the American English spelling, while "grey" is the British English spelling.
+Both are correct, but **[ft_volumesegment](/reference/ft_volumesegment)** uses the American 
+spelling, whereas SPM uses the British spelling.
+{% include markup/end %}
+
+    % read the 1 mm resolution canonical MNI template MRI
+    [ftver, ftpath] = ft_version
+    mri = ft_read_mri(fullfile(ftpath, 'template/anatomy/single_subj_T1_1mm.nii'));
+    mri.coordsys = 'mni';
+    
+    % make a tissue probability map segmentation of CSF, white and gray matter
+    cfg = [];
+    cfg.output = 'tpm';
+    cfg.spmmethod = 'old';
+    seg = ft_volumesegment(cfg,mri);
+    
+    % there are also voxels that are neither CSF, white or gray matter
+    seg.otherwise = ones(seg.dim)-seg.gray-seg.white-seg.csf;
+    
+    % make a matrix with 4 columns, each row is a voxel
+    % the probabilistic sum of all tissue types adds up to 1, or 100%
+    prob = [seg.gray(:) seg.white(:) seg.csf(:) seg.otherwise(:)];
+    
+    % identify which voxels have the highest probability of being gray matter
+    % these are considered to be inside the (binary) grey matter compartment
+    inside = zeros(seg.dim);
+    for voxel = 1:prod(seg.dim)
+    [m,tissue] = max(prob(voxel,:)); % tissue = column, m = number
+    inside(voxel) = tissue == 1;
+    end
+    
+    % these are not needed any more
+    seg = rmfield(seg,'gray');
+    seg = rmfield(seg,'white');
+    seg = rmfield(seg,'csf');
+    seg = rmfield(seg,'otherwise');
+    
+    % add the inside-gray-matter mask
+    seg.inside = inside;
+    
+    %% prepare a source model for gray matter
+    
+    cfg = [];
+    cfg.method = 'basedongrid';
+    cfg.xgrid = -90:3:90;
+    cfg.ygrid = -120:3:90;
+    cfg.zgrid = -90:3:90;
+    cfg.unit = 'mm';
+    
+    % you want the xgrid/ygrid/zgrid numbers to be such, that they cover the entire brain
+    % and that they are nicely symmetric around the origin. In this case they are chosen
+    % such that there is also a dipole exactly at [0, 0, 0].
+    
+    sourcemodel_template = ft_prepare_sourcemodel(cfg);
+    
+    % the resulting source model consists of a 3D grid of dipoles that spans the brain,
+    % but does not specify which ones are inside the grey matter (or the brain) or outside
+    
+    figure
+    ft_plot_mesh(sourcemodel_template.pos)
+    
+    % you should rotate the figure, you will see that it is a square block with many dipoles
+    
+    % we now determine for each dipole what tissue type it is in
+    sourcemodel_template.inside = zeros(prod(sourcemodel_template.dim),1);
+    
+    % the following uses the inverse homogenous transformation to go from head to voxel coordinates
+    
+    for dipole = 1:prod(sourcemodel_template.dim)       % loop over all dipoles
+        thispos = sourcemodel_template.pos(dipole,:);   % the position of this dipole in head coordinates
+        thispos = [thispos 1]';
+        thisvox = round(inv(seg.transform)*thispos);    % the indices of the nearest voxel in the segmented MRI
+        if thisvox(1) < 1 || thisvox(1) > seg.dim(1)
+            % it falls outside the segmented volume
+            sourcemodel_template.inside(dipole) = 0;
+        elseif thisvox(2) < 1 || thisvox(2) > seg.dim(2)
+            % it falls outside the segmented volume
+            sourcemodel_template.inside(dipole) = 0;
+        elseif thisvox(3) < 1 || thisvox(3) > seg.dim(3)
+            % it falls outside the segmented volume
+            sourcemodel_template.inside(dipole) = 0;
+        else
+            % look up in the segmented volume whether the nearest voxel is gray matter
+            sourcemodel_template.inside(dipole) = seg.inside(thisvox(1), thisvox(2), thisvox(3));
+        end
+    end
+    
+    % convert it into a logical array with true/false values
+    sourcemodel_template.inside = logical(sourcemodel_template.inside);
+    
+    % make a plot of the dipoles that are inside the gray matter
+    figure
+    ft_plot_mesh(sourcemodel_template.pos(sourcemodel_template.inside,:))
+    
+    % you should rotate the figure, you will see that it is a brain-shaped cloud of dipoles, only in gray matter
+    

faq/development/datatype.md

@@ -47,7 +47,7 @@ An example of a source structure obtained after performing a frequency domain so
          inside: [1x3415 double]           indices to the positions at which the source activity is actually estimated
         outside: [1x3317 double]           indices to the positions at which the source activity has not been estimated
 
-            dim: [xdim ydim zdim]          if the positions are described as a 3D regular grid, this contains the
+            dim: [xdim ydim zdim]          if the positions are described as a regularly spaced 3D grid, this contains the
                                              dimensionality of the 3D volume
             vol: [1x1 struct]              volume conductor model
       cumtapcnt: [10x1 double]             information about the number of tapers per original trial
@@ -69,7 +69,7 @@ This is the new definition of a data structure that represents data correspondin
             pos: [6732x3 double]       positions at which the source activity could have been estimated
          inside: [6732x1 logical]      logical vector of positions at which the source activity is actually estimated
 
-            dim: [xdim ydim zdim]      if the positions can be described as a 3D regular grid, this contains the
+            dim: [xdim ydim zdim]      if the positions can be described as a regularly spaced 3D grid, this contains the
                                          dimensionality of the 3D volume
             vol: [1x1 struct]          volume conductor model
       cumtapcnt: [120x1 double]        information about the number of tapers per original trial

faq/source/sticking_out.md

@@ -26,7 +26,7 @@ For the singleshell, singlesphere or localspheres models for MEG it is not a pro
 
 ## Solution
 
-If you are working with regular 3D grids as source model, you don't have to do anything. FieldTrip will check for every source whether it is inside or outside the innermost compartment of the volume conductor, and all sources outside will be flagged and not used in the source reconstruction. FieldTrip will actually do the same check for source models based on a triangulated cortical sheets. You can use the `cfg.inwardshift` option to ensure that sources that are on the inside - but too close to the surface for BEM - are also flagged as outside.
+If you are working with regularly spaced 3D grids as source model, you don't have to do anything. FieldTrip will check for every source whether it is inside or outside the innermost compartment of the volume conductor, and all sources outside will be flagged and not used in the source reconstruction. FieldTrip will actually do the same check for source models based on a triangulated cortical sheets. You can use the `cfg.inwardshift` option to ensure that sources that are on the inside - but too close to the surface for BEM - are also flagged as outside.
 
 If you do not want to flag any of your dipoles as outside and exclude them from the source estimation, you have to modify your source model. First of all, you should ensure that your cortical sheet is properly aligned with your volume conduction model of the head. The `cfg.moveinward` option in **[ft_prepare_sourcemodel](/reference/ft_prepare_sourcemodel)** takes all dipoles that are outside, and projects then on an inward-shifted surface of the inside-skull. I.e. it deforms the cortical sheet by pushing some vertices inward. This should work both for spherical models (i.e. concentric spheres) and BEM models. In the case of spherical models, it will result in a cortical mesh that is squeezed in to a sphere shape.
 

template/sourcemodel.md

@@ -20,7 +20,7 @@ When you do source reconstruction with dipole fit methods (as implemented in **[
 When doing source reconstruction with beamformers, people typically scan the brain volume where dipoles are defined on a regular 3D grid, with a regular spacing between the dipole locations. These grids are usually optimized to the individual anatomy of the participant. To facilitate group analysis, however, a clever strategy is to use a template grid (based on a template anatomical volume) that will be (linearly or non-linearly) warped to the individual participant's anatomy. Although this may lead to irregular spacing between the dipole locations (both across and potentially within participants), the dipole locations are directly comparable across participants, because they coincide in standard space.
 The FieldTrip template directory provides a set of sourcemodels defined on regular 3D-grids that are constructed from the MNI-template anatomy from SPM. These template sourcemodels can subsequently be used to be inverse normalized to the individual participant's anatomy. See this [example](/example/sourcemodel_aligned2mni) for more information.
 
-Template source models with the varying dipole spacing (4, 5, 6, 7.5, 8 and 10 mm) on a regular 3-D grid are available in fieldtrip/template/sourcemodel director
+Template source models with the varying dipole spacing (4, 5, 6, 7.5, 8 and 10 mm) on a regular 3-D grid are released along with the FieldTrip toolbox and available in the `fieldtrip/template/sourcemodel` directory:
 
 - standard_sourcemodel3d10mm.mat
 - standard_sourcemodel3d4mm.mat
@@ -29,7 +29,7 @@ Template source models with the varying dipole spacing (4, 5, 6, 7.5, 8 and 10 m
 - standard_sourcemodel3d7point5mm.mat
 - standard_sourcemodel3d8mm.mat
 
-To load and visualize the 3D regular grids, you can do for example
+To load and visualize the regularly spaced 3D grids, you can do
 
     load standard_sourcemodel3d5mm.mat
 
@@ -44,7 +44,7 @@ or you can use
     figure
     ft_plot_mesh(sourcemodel)
 
-You will notice that the 3D regular grids are not that interesting to look at; since there are no connecting elements (triangles), **[ft_plot_mesh](/reference/plotting/ft_plot_mesh)** will only show the vertices.
+You will notice that the regularly spaced 3D grids are not that interesting to look at; since there are no connecting elements (triangles), **[ft_plot_mesh](/reference/plotting/ft_plot_mesh)** will only show the vertices.
 
 ### Distributed source models with MNE
 

workshop/neuroimaging2-2425/sr4.md

@@ -354,7 +354,7 @@ row in the 'avg’ field belongs to which channel).
 Next we will make a source reconstruction using the 'mne’ method of FieldTrip’s
 ft_sourceanalysis function. Before we can do this, we need to define our source
 model, i.e., the set of locations that we assume to be active. For now we assume
-that the active dipoles are distributed on a 3D regular grid, with a spacing
+that the active dipoles are distributed on a regular 3D grid, with a spacing
 of 1 cm between the dipoles:
 
     sourcemodel = ni2_sourcemodel('type', 'grid', 'resolution', 1);