Merge pull request #1919 from schoffelen/femstuff

ENH - working on a reported issue w.r.t. ill-defined reference for si…

external/simbio/sb_calc_vecx.m

@@ -1,4 +1,4 @@
-function vecx = sb_calc_vecx(stiff,vecb,ref);
+function vecx = sb_calc_vecx(stiff,vecb,ref)
 
 % SB_CALC_VECX
 %
@@ -8,6 +8,6 @@
 vecdi(ref) = 1;
 vecva = zeros(size(stiff,1),1);
 [stiff, vecb] = sb_set_bndcon(stiff,vecb,vecdi,vecva);
-clear vecdi, vecva;
+clear vecdi vecva;
 vecx = sb_solve(stiff,vecb);
 end

external/simbio/sb_set_bndcon.m

@@ -6,7 +6,7 @@
 
 dia = diag(stiff);
 stiff = stiff - diag(dia);
-[indexi indexj s] = find(stiff);
+[indexi, indexj, s] = find(stiff);
 clear stiff;
 dind = find(dirinode > 0);
 indi = find(ismember(indexi,dind));

external/simbio/sb_solve.m

@@ -1,14 +1,14 @@
-function x = sb_solve(sysmat,vecb);
+function x = sb_solve(sysmat,vecb)
 
 % SB_SOLVE
 %
 % $Id$
 
 %scalen
-disp('Scaling stiffnes matrix...')
+disp('Scaling stiffness matrix...')
 dkond = 1./(sqrt(diag(sysmat)));
 vecb = vecb.*dkond;
-[indexi indexj s] = find(sysmat);
+[indexi, indexj, s] = find(sysmat);
 sys_size = size(sysmat,1);
 clear sysmat;
 s = (s.*dkond(indexi)).*dkond(indexj);

external/simbio/sb_transfer.m

@@ -33,6 +33,7 @@
 %TODO: add possibility for parallel computation
 %matlabpool local 2;
 for i=2:length(vol.elecnodes)
+    % NOTE: the loop starts at 2, which seems intentional, because the reference is the first electrode
     str = ['Electrode ',num2str(i),' of ',num2str(size(vol.elecnodes,1))];
     disp(str)
     vecb = zeros(size(vol.stiff,1),1);

ft_electroderealign.m

@@ -189,21 +189,9 @@
 cfg = ft_checkconfig(cfg, 'renamedval', {'warp', 'homogenous', 'rigidbody'});
 cfg = ft_checkconfig(cfg, 'renamedval', {'warp', 'homogeneous', 'rigidbody'});
 
-% set the defaults
-cfg.warp          = ft_getopt(cfg, 'warp', 'rigidbody');
-cfg.channel       = ft_getopt(cfg, 'channel',  'all');
-cfg.keepchannel   = ft_getopt(cfg, 'keepchannel', 'no');
-cfg.feedback      = ft_getopt(cfg, 'feedback', 'no');
-cfg.casesensitive = ft_getopt(cfg, 'casesensitive', 'no');
-cfg.headshape     = ft_getopt(cfg, 'headshape', []);     % for triangulated head surface, without labels
+% set these defaults already here, the rest will follow later
 cfg.target        = ft_getopt(cfg, 'target',  []);       % for electrodes or fiducials, always with labels
 cfg.coordsys      = ft_getopt(cfg, 'coordsys');          % this allows for automatic template fiducial placement
-
-if isempty(cfg.target)
-  % remove the field, otherwise ft_checkconfig will complain
-  cfg = rmfield(cfg, 'target');
-end
-
 if ~isempty(cfg.coordsys) && isempty(cfg.target)
   % set the template fiducial locations according to the coordinate system
   switch lower(cfg.coordsys)
@@ -219,20 +207,32 @@
     otherwise
       ft_error('the %s coordinate system is not automatically supported, please specify fiducial details in cfg.target')
   end
+elseif isempty(cfg.target)
+  % remove the field, otherwise ft_checkconfig will complain depending on
+  % the settings for checkconfig
+  cfg = rmfield(cfg, 'target');
 end
 
 % ensure that the right cfg options have been set corresponding to the method
 switch cfg.method
-  case 'template'        % realign the sensors to match a template set
+  case 'template'      % realign the sensors to match a template set
     cfg = ft_checkconfig(cfg, 'required', 'target', 'forbidden', 'headshape');
   case 'headshape'     % realign the sensors to fit the head surface
     cfg = ft_checkconfig(cfg, 'required', 'headshape', 'forbidden', 'target');
-  case 'fiducial'        % realign using the NAS, LPA and RPA fiducials
+  case 'fiducial'      % realign using the NAS, LPA and RPA fiducials
     cfg = ft_checkconfig(cfg, 'required', 'target', 'forbidden', 'headshape');
-  case 'moveinward'      % moves eletrodes inward
+  case 'moveinward'    % moves eletrodes inward
     cfg = ft_checkconfig(cfg, 'required', 'moveinward');
 end % switch cfg.method
 
+% set the rest of the defaults, this is to avoid confusion in ft_checkconfig w.r.t. to the method specific checks
+cfg.warp          = ft_getopt(cfg, 'warp', 'rigidbody');
+cfg.channel       = ft_getopt(cfg, 'channel',  'all');
+cfg.keepchannel   = ft_getopt(cfg, 'keepchannel', 'no');
+cfg.feedback      = ft_getopt(cfg, 'feedback', 'no');
+cfg.casesensitive = ft_getopt(cfg, 'casesensitive', 'no');
+cfg.headshape     = ft_getopt(cfg, 'headshape', []);     % for triangulated head surface, without labels
+
 if strcmp(cfg.method, 'fiducial') && isfield(cfg, 'warp') && ~isequal(cfg.warp, 'rigidbody')
   ft_warning('The method ''fiducial'' implies a rigid body tramsformation. See also http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1722');
   cfg.warp = 'rigidbody';
@@ -298,9 +298,14 @@
   end
   clear tmp
   for i=1:Ntemplate
-    % ensure up-to-date sensor description
-    % ensure that the units are consistent with the electrodes
-    tmp(i) = ft_convert_units(ft_datatype_sens(target(i)), elec.unit);
+    try
+      % ensure up-to-date sensor description
+      % ensure that the units are consistent with the electrodes
+      tmp(i) = ft_convert_units(ft_datatype_sens(target(i)), elec.unit);
+    catch
+      ft_warning('cannot check the consistency of the units in the fiducials and the electrodes, you need to ensure this yourself');
+      tmp(i) = target(i);
+    end
   end
   target = tmp;
 end
@@ -561,19 +566,19 @@
 
   tmpcfg = [];
   tmpcfg.individual.elec = elec;
-  if isfield(cfg, 'headshape') && ~isempty(cfg.headshape)
-    tmpcfg.template.headshape = cfg.headshape;
+  if useheadshape
+    tmpcfg.template.headshape = headshape;
   end
-  if isfield(cfg, 'target') && ~isempty(cfg.target)
-    if iscell(cfg.target)
-      if numel(cfg.target)>1
+  if usetarget
+    if iscell(target)
+      if numel(target)>1
         ft_notice('computing the average electrode positions');
-        tmpcfg.template.elec = ft_average_sens(cfg.target);
+        tmpcfg.template.elec = ft_average_sens(target);
       else
-        tmpcfg.template.elec = cfg.target{1};
+        tmpcfg.template.elec = target{1};
       end
     elseif isstruct(cfg.target)
-      tmpcfg.template.elec = cfg.target;
+      tmpcfg.template.elec = target;
     end
     tmpcfg.template.elecstyle = {'facecolor', 'blue'};
     ft_info('plotting the target electrodes in blue');
@@ -665,7 +670,7 @@
       elec_realigned.(cfg.warp) = norm.m;
     end
 
-  case  {'fiducial' 'interactive'}
+  case  {'fiducial', 'interactive'}
     % the transformation is a 4x4 homogenous matrix
     % apply the transformation to the original complete set of sensors
     elec_realigned = ft_transform_geometry(norm.m, elec_original);
@@ -735,10 +740,6 @@
   elec_realigned.elecpos = [elec_realigned.elecpos; elec_original.elecpos(idx,:)];
 end
 
-% channel positions are identical to the electrode positions (this was checked at the start)
-elec_realigned.chanpos = elec_realigned.elecpos;
-elec_realigned.tra = eye(numel(elec_realigned.label));
-
 % copy over unit, chantype, chanunit, and tra information in case this was not already done
 if ~isfield(elec_realigned, 'unit') && isfield(elec_original, 'unit')
   elec_realigned.unit = elec_original.unit;

test/test_example_fem.m

@@ -0,0 +1,77 @@
+function test_example_fem
+
+% WALLTIME 08:00:00
+% MEM 12gb
+% DEPENDENCY ft_prepare_headmodel ft_prepare_mesh ft_datatype_segmentation
+
+[ftver, ftpath] = ft_version;
+
+%% read mri
+mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/Subject01.mri'));
+
+%% segmentation
+cfg          = [];
+cfg.output   = {'gray', 'white', 'csf', 'skull', 'scalp'};
+segmentedmri = ft_volumesegment(cfg, mri);
+
+%% mesh
+cfg        = [];
+cfg.shift  = 0.3;
+cfg.method = 'hexahedral';
+mesh       = ft_prepare_mesh(cfg,segmentedmri);
+
+%% volume conductor
+cfg              = [];
+cfg.method       = 'simbio';
+cfg.conductivity = [0.33 0.14 1.79 0.01 0.43];
+headmodel        = ft_prepare_headmodel(cfg, mesh);
+
+%% electrode alignment
+elec = ft_read_sens(fullfile(ftpath,'template/electrode/standard_1020.elc'));
+
+nas = mri.hdr.fiducial.head.nas;
+lpa = mri.hdr.fiducial.head.lpa;
+rpa = mri.hdr.fiducial.head.rpa;
+
+%
+fiducials.pos     = [nas; lpa; rpa];
+fiducials.label   = {'Nz','LPA','RPA'};
+fiducials.unit    = 'mm';
+
+cfg          = [];
+cfg.method   = 'fiducial';
+cfg.target   = fiducials;
+cfg.elec     = elec;
+cfg.fiducial = {'Nz', 'LPA', 'RPA'};
+elec_align   = ft_electroderealign(cfg);
+
+% add 12 mm to x-axis
+n=size(elec_align.chanpos,1);
+for i=1:n
+ elec_align.chanpos(i,1)=elec_align.chanpos(i,1)+12;
+ elec_align.elecpos(i,1)=elec_align.elecpos(i,1)+12;
+end
+
+%% make the sourcemodel/grid
+
+% At the moment the sourcemodel is defined prior
+% to the leadfield because ft_prepare_sourcemodel does not automatically create
+% a sourcemodel based on a hexahedral headmodel.
+
+cfg                 = [];
+cfg.mri             = mri;
+cfg.resolution      = 3;
+sourcemodel         = ft_prepare_sourcemodel(cfg);
+sourcemodel         = ft_convert_units(sourcemodel, headmodel.unit);
+sourcemodel.inside(500:end) = false; % do only a few dipoles
+
+cfg            = [];
+cfg.headmodel  = headmodel;
+cfg.elec       = elec_align;
+cfg.sourcemodel = sourcemodel;
+cfg.channel    = elec_align.label(4:6); % do only a few channels, to save time
+lf             = ft_prepare_leadfield(cfg);
+
+m = cellfun(@mean,lf.leadfield(lf.inside),'uniformoutput',false)';
+m = cat(1, m{:});
+assert(all(m(:)<eps));

test/test_tutorial_headmodel_eeg_fem.m

@@ -0,0 +1,84 @@
+function test_tutorial_headmodel_eeg_fem
+
+% WALLTIME 08:00:00
+% MEM 12gb
+% DEPENDENCY ft_prepare_headmodel ft_prepare_mesh ft_datatype_segmentation
+
+mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/Subject01.mri'));
+
+% this needs to be done before reslicing
+nas = mri.hdr.fiducial.mri.nas;
+lpa = mri.hdr.fiducial.mri.lpa;
+rpa = mri.hdr.fiducial.mri.rpa;
+
+vox2head = mri.transform;
+nas = ft_warp_apply(vox2head, nas, 'homogenous');
+lpa = ft_warp_apply(vox2head, lpa, 'homogenous');
+rpa = ft_warp_apply(vox2head, rpa, 'homogenous');
+
+cfg     = [];
+cfg.dim = mri.dim;
+mri     = ft_volumereslice(cfg,mri);
+
+cfg           = [];
+cfg.output    = {'gray','white','csf','skull','scalp'};
+segmentedmri  = ft_volumesegment(cfg, mri);
+
+seg_i = ft_datatype_segmentation(segmentedmri,'segmentationstyle','indexed');
+
+cfg        = [];
+cfg.shift  = 0.3;
+cfg.method = 'hexahedral';
+mesh = ft_prepare_mesh(cfg,seg_i);
+
+cfg        = [];
+cfg.method ='simbio';
+cfg.conductivity = [0.33 0.14 1.79 0.01 0.43];   % order follows mesh.tissyelabel
+headmodel  = ft_prepare_headmodel(cfg, mesh);
+
+ft_plot_mesh(mesh, 'surfaceonly', 'yes');
+
+elec = ft_read_sens('standard_1020.elc');
+
+figure
+hold on
+ft_plot_mesh(mesh,'surfaceonly','yes','vertexcolor','none','edgecolor','none','facecolor',[0.5 0.5 0.5],'face alpha',0.7)
+camlight
+
+ft_plot_sens(elec);
+
+
+% Then, we determine the translation and rotation that is needed to get the position of the fiducials in the electrode structure (defined with labels 'Nz', 'LPA', 'RPA') to their counterparts in the CTF head coordinate system that we acquired from the anatomical mri (nas, lpa, rpa).
+%
+% create a structure similar to a template set of electrodes
+fid.pos           = [nas; lpa; rpa];       % CTF head coordinates of fiducials
+fid.label         = {'Nz','LPA','RPA'};    % use the same labels as those in elec
+fid.unit          = 'mm';                  % use the same units as those in mri
+
+% alignment
+cfg               = [];
+cfg.method        = 'fiducial';
+cfg.elec          = elec;                  % the electrodes we want to align
+cfg.target        = fid;                   % the template we want to align to
+cfg.fiducial      = {'Nz', 'LPA', 'RPA'};  % labels of fiducials in fid and in elec
+elec_aligned      = ft_electroderealign(cfg);
+
+% We can check the alignment by plotting together the scalp surface with the electrodes.
+%
+figure;
+hold on;
+ft_plot_mesh(mesh,'surfaceonly','yes','vertexcolor','none','edgecolor','none','facecolor',[0.5 0.5 0.5],'face alpha',0.5)
+camlight
+ft_plot_sens(elec_aligned);
+
+%
+%
+% The alignment is much better, but not perfect. Some of the electrodes are below the head surface in the front, while the electrodes in the back do not fit tightly to the head. The remaining misalignment is due to the use of different conventions to [define the fiducials](/faq/how_are_the_lpa_and_rpa_points_defined). We can improve the alignment of the electrodes interactively.
+%
+%% ## Interactive alignment
+%
+cfg          = [];
+cfg.method   = 'interactive';
+cfg.elec     = elec_aligned;
+cfg.headshape = headmodel;
+elec_aligned = ft_electroderealign(cfg);

utilities/ft_datatype_sens.m

@@ -76,7 +76,7 @@
 % (2010) Added support for bipolar or otherwise more complex linear combinations
 %  of EEG electrodes using sens.tra, similar to MEG.
 %
-% (2009) Noice reduction has been added for MEG systems in the balance field.
+% (2009) Noise reduction has been added for MEG systems in the balance field.
 %
 % (2006) The optional fields sens.type and sens.unit were added.
 %
@@ -417,7 +417,18 @@
         sens.elecpos = sens.pnt; sens = rmfield(sens, 'pnt');
       end
     end
-
+
+    if isfield(sens, 'pos')
+      if ismeg
+        % sensor description is a MEG sensor-array, containing oriented coils
+        sens.coilpos = sens.pos; sens = rmfield(sens, 'pos');
+        sens.coilori = sens.ori; sens = rmfield(sens, 'ori');
+      else
+        % sensor description is something else, EEG/ECoG/sEEG, etc
+        sens.elecpos = sens.pos; sens = rmfield(sens, 'pos');
+      end
+    end
+
     if ~isfield(sens, 'chanpos')
       if ismeg
         % sensor description is a MEG sensor-array, containing oriented coils