improve documentation for dykstra2012 algorithm

Author: gpiantoni

ft_electroderealign.m

@@ -53,6 +53,9 @@
 %                        'nonlin3'         apply a 3rd order non-linear warp
 %                        'nonlin4'         apply a 4th order non-linear warp
 %                        'nonlin5'         apply a 5th order non-linear warp
+%                        'dykstra2012'     non-linear wrap only for headshape
+%                                          method useful for projecting ECoG onto
+%                                          cortex hull.
 %   cfg.channel        = Nx1 cell-array with selection of channels (default = 'all'),
 %                        see  FT_CHANNELSELECTION for details
 %   cfg.fiducial       = cell-array with the name of three fiducials used for
@@ -85,7 +88,21 @@
 %                        single triangulated boundary, or a Nx3 matrix with surface
 %                        points
 %
-% See also FT_READ_SENS, FT_VOLUMEREALIGN, FT_INTERACTIVEREALIGN
+% If you want to align ECoG electrodes to the pial surface, you first need to
+% compute the cortex hull with FT_PREPARE_MESH. dykstra2012 uses algorithm
+% described in Dykstra et al. (2012, Neuroimage) in which electrodes are
+% projected onto pial surface while minimizing the displacement of the
+% electrodes from original location and maintaining the grid shape. It relies
+% on the optimization toolbox.
+%   cfg.method         = 'headshape'
+%   cfg.warp           = 'dykstra2012'
+%   cfg.headshape      = a filename containing headshape, a structure containing a
+%                        single triangulated boundary, or a Nx3 matrix with surface
+%                        points
+%   cfg.feedback       = 'yes' or 'no' (feedback includes the output of the iteration
+%                        procedure.
+%
+% See also FT_READ_SENS, FT_VOLUMEREALIGN, FT_INTERACTIVEREALIGN, FT_PREPARE_MESH
 
 % Copyright (C) 2005-2015, Robert Oostenveld
 %
@@ -297,13 +314,13 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 if strcmp(cfg.method, 'template')
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  
+
   % determine electrode selection and overlapping subset for warping
   cfg.channel = ft_channelselection(cfg.channel, elec.label);
   for i=1:Ntemplate
     cfg.channel = ft_channelselection(cfg.channel, target(i).label);
   end
-  
+
   % make consistent subselection of electrodes
   [cfgsel, datsel] = match_str(cfg.channel, elec.label);
   elec.label = elec.label(datsel);
@@ -313,18 +330,18 @@
     target(i).label   = target(i).label(datsel);
     target(i).elecpos = target(i).elecpos(datsel,:);
   end
-  
+
   % compute the average of the target electrode positions
   average = ft_average_sens(target);
-  
+
   fprintf('warping electrodes to average template... '); % the newline comes later
   [norm.elecpos, norm.m] = ft_warp_optim(elec.elecpos, average.elecpos, cfg.warp);
   norm.label = elec.label;
-  
+
   dpre  = mean(sqrt(sum((average.elecpos - elec.elecpos).^2, 2)));
   dpost = mean(sqrt(sum((average.elecpos - norm.elecpos).^2, 2)));
   fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);
-  
+
   if strcmp(cfg.feedback, 'yes')
     % create an empty figure, continued below...
     figure
@@ -334,28 +351,28 @@
     xlabel('x')
     ylabel('y')
     zlabel('z')
-    
+
     % plot all electrodes before warping
     ft_plot_sens(elec, 'r*');
-    
+
     % plot all electrodes after warping
     ft_plot_sens(norm, 'm.', 'label', 'label');
-    
+
     % plot the template electrode locations
     ft_plot_sens(average, 'b.');
-    
+
     % plot lines connecting the input and the realigned electrode locations with the template locations
     my_line3(elec.elecpos, average.elecpos, 'color', 'r');
     my_line3(norm.elecpos, average.elecpos, 'color', 'm');
   end
-  
+
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 elseif strcmp(cfg.method, 'headshape')
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  
+
   % determine electrode selection and overlapping subset for warping
   cfg.channel = ft_channelselection(cfg.channel, elec.label);
-  
+
   % make subselection of electrodes
   [cfgsel, datsel] = match_str(cfg.channel, elec.label);
   elec.label   = elec.label(datsel);
@@ -372,14 +389,14 @@
     dpost = ft_warp_error(norm.m, elec.elecpos, headshape, cfg.warp);
     fprintf('mean distance prior to warping %f, after warping %f\n', dpre, dpost);
   end
-  
+
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 elseif strcmp(cfg.method, 'fiducial')
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  
+
   % the fiducials have to be present in the electrodes and in the template set
   label = intersect(lower(elec.label), lower(target.label));
-  
+
   if ~isfield(cfg, 'fiducial') || isempty(cfg.fiducial)
     % try to determine the names of the fiducials automatically
     option1 = {'nasion' 'left' 'right'};
@@ -405,17 +422,17 @@
     end
   end
   fprintf('matching fiducials {''%s'', ''%s'', ''%s''}\n', cfg.fiducial{1}, cfg.fiducial{2}, cfg.fiducial{3});
-  
+
   % determine electrode selection
   cfg.channel = ft_channelselection(cfg.channel, elec.label);
   [cfgsel, datsel] = match_str(cfg.channel, elec.label);
   elec.label     = elec.label(datsel);
   elec.elecpos   = elec.elecpos(datsel,:);
-  
+
   if length(cfg.fiducial)~=3
     error('you must specify exaclty three fiducials');
   end
-  
+
   % do case-insensitive search for fiducial locations
   nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
   lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
@@ -426,11 +443,11 @@
   elec_nas = elec.elecpos(nas_indx,:);
   elec_lpa = elec.elecpos(lpa_indx,:);
   elec_rpa = elec.elecpos(rpa_indx,:);
-  
+
   % FIXME change the flow in the remainder
   % if one or more template electrode sets are specified, then align to the average of those
   % if no template is specified, then align so that the fiducials are along the axis
-  
+
   % find the matching fiducials in the template and average them
   tmpl_nas = nan(Ntemplate,3);
   tmpl_lpa = nan(Ntemplate,3);
@@ -449,21 +466,21 @@
   tmpl_nas = mean(tmpl_nas,1);
   tmpl_lpa = mean(tmpl_lpa,1);
   tmpl_rpa = mean(tmpl_rpa,1);
-  
+
   % realign both to a common coordinate system
   elec2common  = ft_headcoordinates(elec_nas, elec_lpa, elec_rpa);
   templ2common = ft_headcoordinates(tmpl_nas, tmpl_lpa, tmpl_rpa);
-  
+
   % compute the combined transform
   norm         = [];
   norm.m       = templ2common \ elec2common;
-  
+
   % apply the transformation to the fiducials as sanity check
   norm.elecpos(1,:) = ft_warp_apply(norm.m, elec_nas, 'homogeneous');
   norm.elecpos(2,:) = ft_warp_apply(norm.m, elec_lpa, 'homogeneous');
   norm.elecpos(3,:) = ft_warp_apply(norm.m, elec_rpa, 'homogeneous');
   norm.label        = cfg.fiducial;
-  
+
   nas_indx = match_str(lower(elec.label), lower(cfg.fiducial{1}));
   lpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{2}));
   rpa_indx = match_str(lower(elec.label), lower(cfg.fiducial{3}));
@@ -473,7 +490,7 @@
   rpa_indx = match_str(lower(norm.label), lower(cfg.fiducial{3}));
   dpost = mean(sqrt(sum((norm.elecpos([nas_indx lpa_indx rpa_indx],:) - [tmpl_nas; tmpl_lpa; tmpl_rpa]).^2, 2)));
   fprintf('mean distance between fiducials prior to realignment %f, after realignment %f\n', dpre, dpost);
-  
+
   if strcmp(cfg.feedback, 'yes')
     % create an empty figure, continued below...
     figure
@@ -483,23 +500,23 @@
     xlabel('x')
     ylabel('y')
     zlabel('z')
-    
+
     % plot the first three electrodes before transformation
     my_plot3(elec.elecpos(1,:), 'r*');
     my_plot3(elec.elecpos(2,:), 'r*');
     my_plot3(elec.elecpos(3,:), 'r*');
     my_text3(elec.elecpos(1,:), elec.label{1}, 'color', 'r');
     my_text3(elec.elecpos(2,:), elec.label{2}, 'color', 'r');
     my_text3(elec.elecpos(3,:), elec.label{3}, 'color', 'r');
-    
+
     % plot the template fiducials
     my_plot3(tmpl_nas, 'b*');
     my_plot3(tmpl_lpa, 'b*');
     my_plot3(tmpl_rpa, 'b*');
     my_text3(tmpl_nas, ' nas', 'color', 'b');
     my_text3(tmpl_lpa, ' lpa', 'color', 'b');
     my_text3(tmpl_rpa, ' rpa', 'color', 'b');
-    
+
     % plot all electrodes after transformation
     my_plot3(norm.elecpos, 'm.');
     my_plot3(norm.elecpos(1,:), 'm*');
@@ -509,7 +526,7 @@
     my_text3(norm.elecpos(2,:), norm.label{2}, 'color', 'm');
     my_text3(norm.elecpos(3,:), norm.label{3}, 'color', 'm');
   end
-  
+
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 elseif strcmp(cfg.method, 'interactive')
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -540,14 +557,14 @@
   clear global norm
   norm = tmp;
   clear tmp
-  
+
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 elseif strcmp(cfg.method, 'project')
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   [dum, prj] = project_elec(elec.elecpos, headshape.pos, headshape.tri);
   % replace the electrodes with the projected version
   elec.elecpos = prj;
-  
+
 else
   error('unknown method');
 end % if method
@@ -560,7 +577,7 @@
     if strcmp(lower(cfg.warp), 'dykstra2012')
       elec_realigned = norm;
       elec_realigned.unit = elec_original.unit;
-      
+
     else
       % the transformation is a linear or non-linear warp, i.e. a vector
       try
@@ -574,20 +591,20 @@
       end
       % remember the transformation
       elec_realigned.(cfg.warp) = norm.m;
-      
+
     end
-    
+
   case  {'fiducial' 'interactive'}
     % the transformation is a 4x4 homogenous matrix
     % apply the transformation to the original complete set of sensors
     elec_realigned = ft_transform_sens(norm.m, elec_original);
     % remember the transformation
     elec_realigned.homogeneous = norm.m;
-    
+
   case 'project'
     % nothing to be done
     elec_realigned = elec;
-    
+
   otherwise
     error('unknown method');
 end

ft_prepare_mesh.m

@@ -3,13 +3,16 @@
 % FT_PREPARE_MESH creates a triangulated surface mesh for the volume
 % conduction model. The mesh can either be selected manually from raw
 % mri data or can be generated starting from a segmented volume
-% information stored in the mri structure. The result is a bnd
-% structure which contains the information about all segmented surfaces
-% related to mri and are expressed in world coordinates.
+% information stored in the mri structure. FT_PREPARE_MESH can be used
+% to create a cortex hull, i.e. the smoothed envelope around the pial
+% surface created by freesurfer. The result is a bnd structure which
+% contains the information about all segmented surfaces related to mri
+% sand are expressed in world coordinates.
 %
 % Use as
 %   bnd = ft_prepare_mesh(cfg, mri)
 %   bnd = ft_prepare_mesh(cfg, seg)
+%   bnd = ft_prepare_mesh(cfg)  # for cortexhull
 %
 % Configuration options:
 %   cfg.method      = string, can be 'interactive', 'projectmesh', 'iso2mesh', 'isosurface',
@@ -20,6 +23,12 @@
 %   cfg.headshape   = (optional) a filename containing headshape, a Nx3 matrix with surface
 %                     points, or a structure with a single or multiple boundaries
 %
+% Method 'cortexhull' has its own specific configuration options:
+%   cfg.headshape   = a filename containing the pial surface computed by
+%                     freesurfer recon-all ('/path/to/surf/lh.pial')
+%   cfg.
+%
+%
 % To facilitate data-handling and distributed computing you can use
 %   cfg.inputfile   =  ...
 %   cfg.outputfile  =  ...
@@ -40,6 +49,11 @@
 %   cfg.numvertices = [800, 1600, 2400];
 %   bnd             = ft_prepare_mesh(cfg, segmentation);
 %
+%   cfg             = [];
+%   cfg.method      = 'cortexhull';
+%   cfg.headshape   = '/path/to/surf/lh.pial';
+%   cortex_hull     = ft_prepare_mesh(cfg);
+%
 % See also FT_VOLUMESEGMENT, FT_PREPARE_HEADMODEL, FT_PLOT_MESH
 
 % Undocumented functionality: at this moment it allows for either
@@ -175,8 +189,8 @@
     end
 
   case 'cortexhull'
-    bnd = prepare_cortexhull(cfg);      
-      
+    bnd = prepare_cortexhull(cfg);
+
   otherwise
     error('unsupported cfg.method')
 end