diff --git a/external/simbio/sb_calc_vecx.m b/external/simbio/sb_calc_vecx.m index cf633b61a1..c998a5116f 100644 --- a/external/simbio/sb_calc_vecx.m +++ b/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 diff --git a/external/simbio/sb_set_bndcon.m b/external/simbio/sb_set_bndcon.m index 092736f883..dd98eb0658 100644 --- a/external/simbio/sb_set_bndcon.m +++ b/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)); diff --git a/external/simbio/sb_solve.m b/external/simbio/sb_solve.m index 3a52d7119f..33c2c560a7 100644 --- a/external/simbio/sb_solve.m +++ b/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); diff --git a/external/simbio/sb_transfer.m b/external/simbio/sb_transfer.m index bc8a1d21fb..e92f727b0e 100644 --- a/external/simbio/sb_transfer.m +++ b/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); diff --git a/ft_electroderealign.m b/ft_electroderealign.m index ea2adef6eb..f21e124210 100644 --- a/ft_electroderealign.m +++ b/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; diff --git a/test/test_example_fem.m b/test/test_example_fem.m new file mode 100644 index 0000000000..7d86d8339f --- /dev/null +++ b/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(:)