Back to the main page.
Bug 3119 - fitting 2 moving dipoles with constr.sequential = 1 probably wrong
| Status | ASSIGNED |
| Reported | 2016-04-29 17:22:00 +0200 |
| Modified | 2016-06-02 16:00:34 +0200 |
| Product: | FieldTrip |
| Component: | core |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Linux |
| Importance: | P5 normal |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Diego Lozano Soldevilla - 2016-04-29 17:22:46 +0200
I simulate EEG data to fit two symmetric moving dipoles for every time point of a time interval of interest. I used the constr.sequential = 1 to explore for what might be useful (no documentation and hard to get from the code) and I notice an inconsistency in /inverse/dipole_fit.m. Withing lines 318-330 the sequential thing happens but the 'framesel' variable will be always 0.5, yielding a moment zero because data selection gives and empty matrix: line 329: pinv(lf(:,dipsel))*dat(:,framesel); The thing is that the fitting is done sample-by-sample (ft_dipolefitting.m line 506) so previous line of code should be: pinv(lf(:,dipsel))*dat; However, what I don't understand is that later in ft_dipolefitting.m line 567, the dipole moment is recomputed ignoring the sequential constrain dip(t).mom = pinv(lf)*Vdata(:,t); Looking forward to know how the sequential constrain should behave (fitting one dipole in one frame and another in the next?) and why this second computation of the dipole moment. SUGESTION: would it make sense to add a dipole option on where one can fit two dipoles in a given sample without symmetry constraints. It could be possible by selecting the two lowest error values in the grid.error (around ft_dipolefitting.m line 421) Below you can reproduce the problem using the attached data. % load template mri ftdir = fileparts(which('ft_defaults')); mri = ft_read_mri(fullfile(ftdir,'template','headmodel','standard_mri.mat')); % used later to plot dipoles cfg = []; cfg.resolution = 1; cfg.xrange = [-100 100]; cfg.yrange = [-140 140]; cfg.zrange = [-80 120]; mris = ft_volumereslice(cfg, mri); mris = ft_convert_units(mris, 'cm'); % load template segmented mri and rearrange fields seg = ft_read_mri(fullfile(ftdir,'template','headmodel','standard_seg.mat')); brain = zeros(size(seg.seg)); bri = find(seg.seg==3); brain(bri)=1; clear bri; skull = zeros(size(seg.seg)); ski = find(seg.seg==2); skull(ski)=1; clear ski; scalp = zeros(size(seg.seg)); sci = find(seg.seg==1); scalp(sci)=1; clear sci; seg.brain = brain; clear brain; seg.scalp = scalp; clear scalp; seg.skull = skull; clear skull; seg = rmfield(seg,'seg'); vol = ft_read_vol(fullfile(ftdir,'template','headmodel','standard_bem.mat')); figure; ft_plot_mesh(vol.bnd(1),'edgecolor', 'none', 'facecolor', 'r'); % scalp ft_plot_mesh(vol.bnd(2),'edgecolor', 'none', 'facecolor', 'g'); % skull ft_plot_mesh(vol.bnd(3),'edgecolor', 'none', 'facecolor', 'b'); % brain alpha 0.3; view(132, 14) % get custom electrodes label={'Fp1';'AF7';'AF3';'F1';'TP9';'F5h';'F7';'FT7';'FC5h';'PO9';'FC1';'C1';'C3';'C5';'T7';'TP7';'CP5';'CP3';'CP1';'P1';'I1';'P5h';'P7';'P9';'PO7';'PO3';'O1';'Iz';'Oz';'POz';'Pz';'CPz';'Fpz';'Fp2';'AF8';'AF4';'AFz';'Fz';'F2';'TP10';'F6h';'F8';'FT8';'FC6h';'PO10';'FC2';'FCz';'Cz';'C2';'C4';'C6';'T8';'TP8';'CP6';'CP4';'CP2';'P2';'I2';'P6h';'P8';'P10';'PO8';'PO4';'O2'}; elec = ft_read_sens(fullfile(ftdir,'template','electrode','standard_1005.elc')); [sel1,sel2] = match_str(label,elec.label); % isequal(elec.label(sel2),label(sel1)) chanpos = elec.chanpos(sel2,:); elecpos = elec.elecpos(sel2,:); elec.chanpos = chanpos; elec.elecpos = elecpos; elec.label = label; figure; hold on ft_plot_sens(elec, 'style', 'ob','label','label'); ft_plot_vol(vol, 'facealpha', 0.5, 'edgecolor', 'none'); % "lighting phong" does not work with opacity material dull; camlight; % cfg = []; cfg.elec = ft_convert_units(elec,'cm'); cfg.headmodel = ft_convert_units(vol,'cm'); cfg.reducerank = 3; cfg.channel = 'all'; cfg.grid.resolution = 1; % use a 3-D grid with a 1 cm resolution cfg.grid.unit = 'cm'; grid = ft_prepare_leadfield(cfg); %% load simulated data load datat1.mat cfg = []; cfg.channel = 'all'; cfg.elec = ft_convert_units(elec,'cm'); cfg.senstype = 'eeg'; cfg.latency = [0.2 0.3]; cfg.reducerank = 3; cfg.numdipoles = 1; cfg.grid = ft_convert_units(grid,'cm'); cfg.gridsearch = 'yes'; cfg.nonlinear = 'yes'; cfg.vol = ft_convert_units(vol,'cm'); cfg.model = 'moving'; cfg.numdipoles = 2; cfg.symmetry = 'x'; cfg.dipfit.constr.sequential = 1; dipole = ft_dipolefitting(cfg, datat1); %% figure; ft_plot_slice(mris.anatomy, 'transform', mris.transform, 'location', [-3.1 -4.6 -1.4], 'orientation', [0 1 0], 'resolution', 0.1); ft_plot_slice(mris.anatomy, 'transform', mris.transform, 'location', [-3.1 -4.6 -1.4], 'orientation', [1 0 0], 'resolution', 0.1); ft_plot_slice(mris.anatomy, 'transform', mris.transform, 'location', [-3.1 -4.6 -1.4], 'orientation', [0 0 1], 'resolution', 0.1); axis tight; axis off; for frame=1:17; ft_plot_dipole(dipole.dip(frame).pos(1,:), dipole.dip(frame).mom(1:3,:),'color','r'); ft_plot_dipole(dipole.dip(frame).pos(2,:), dipole.dip(frame).mom(4:6,:),'color','b'); end
Diego Lozano Soldevilla - 2016-04-29 17:25:12 +0200
find the data below https://www.dropbox.com/s/9ftdl2kz14xmb8h/datat1.mat?dl=0
Robert Oostenveld - 2016-06-01 17:17:14 +0200
(In reply to Diego Lozano Soldevilla from comment #1) Hi Diego, Since I am working on something related, I came across your bug report. The sequential option was implemented by me and Arjen (now CC) to fit multiple dipoles at the same time, where each dipole has a different topography. This only makes sense (as far as I know) in combination with the rigidbody option. The specific use case is this: cfg = []; cfg.channel = 'MEGMAG'; cfg.headmodel = headmodel; cfg.dip.unit = 'm'; cfg.dip.pos = polhemus_hpipos; cfg.model = 'regional'; cfg.dipfit.constr.rigidbody = true; cfg.dipfit.constr.sequential = true; cfg.frequency = hpicoilfreq; cfg.gridsearch = 'no'; hpifit_jointly = ft_dipolefitting(cfg, freq); where (in my data) I have 10 HPI head coils, each at its own frequency, for which the fourier topographies are contained in the freq structure. I want to fit 10 magnetic dipoles, whose relative position is fixed and whose starting point is given by the polhemus recording, but the topography is not a linear mixture of the 10 dipole topographies. Rather, each dipole has its own topography (which is spectrally separated). So you should think of the data as Nchan*Ncoil, and the dipole model is Ncoil*3, and I want to move all dipoles jointly (hence rigid body) to minimize the overall error between model and data. Your use case is rather complicated, but I don't think the "sequential" option should be part of it. Actually, sequential and moving don't match. I will add a check and error to the ft_dipolefitting function.
Robert Oostenveld - 2016-06-01 17:20:26 +0200
(In reply to Robert Oostenveld from comment #2) mac011> git commit -a [bug3119-constr-sequential c2dde19] ENH - give an error with incorrect configuration for dipole fitting, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3119 3 files changed, 130 insertions(+), 1 deletion(-) create mode 100644 test/test_bug3119.m
Robert Oostenveld - 2016-06-01 17:43:25 +0200
specifying cfg.symmetry=x together with a cfg.grid that is NOT a symmetric dipole grid is also weird. I followed it through the code, and accidentally it does work correctly. Actually, it is fine I realize... the symmetry is implemented on the 2nd call to ft_prepare_sourcemodel (from ft_dipolefitting). I have improved this a bit. mac011> git commit -a [bug3119-constr-sequential 97a4578] ENH - added some sanity checks, extended test script. See http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3119 3 files changed, 67 insertions(+), 9 deletions(-)
Robert Oostenveld - 2016-06-01 17:50:16 +0200
Created attachment 793 screen shot of data and model this is how I think your simulated data should be fitted cfg = []; cfg.channel = 'all'; cfg.elec = ft_convert_units(elec,'cm'); cfg.grid = ft_convert_units(grid,'cm'); cfg.headmodel = ft_convert_units(vol,'cm'); cfg.senstype = 'eeg'; cfg.latency = [0.2 0.3]; cfg.reducerank = 3; cfg.gridsearch = 'yes'; cfg.nonlinear = 'yes'; cfg.model = 'regional'; cfg.numdipoles = 2; cfg.symmetry = 'x'; dipole = ft_dipolefitting(cfg, datat1); The fit is not perfect, but see original = rmfield(dipole, 'dip'); original.avg = dipole.Vdata; fitted = rmfield(dipole, 'dip'); fitted.avg = dipole.Vmodel; cfg = []; cfg.layout = 'elec1005.lay'; ft_multiplotER(cfg, original, fitted);
Robert Oostenveld - 2016-06-01 17:51:24 +0200
please note that I do not rule out that there are errors due to inconsistent model handling in ft_dipolefitting in these sections: 1) the grid search 2) the nonlinear part 3) the post-hoc recomputing of the fit
Robert Oostenveld - 2016-06-01 17:52:12 +0200
(In reply to Robert Oostenveld from comment #5) mac011> git commit -a [bug3119-constr-sequential dc689a0] ENH - added the "best" configuration to the test script, and added visualisation of the result 1 file changed, 31 insertions(+)
Robert Oostenveld - 2016-06-01 17:54:31 +0200
I merged the changes with https://github.com/fieldtrip/fieldtrip/pull/168
Diego Lozano Soldevilla - 2016-06-02 00:12:43 +0200
(In reply to Robert Oostenveld from comment #8) Hi Robert, Thanks for your thorough answer. Now I see the rigid body option for the headlocalizer application. Regarding my proposal of fitting two dipoles without symmetry constraints I still think makes sense. I think BESA has a way to fit a single dipole that explains as much variance as possible and then it fits a second one to increase the variance explained (here the reason I reported about the cfg.dipfit.constr.sequential = true;). This is the reason I thought about selecting the two lowest error values in the grid.error. Do you have any suggestions about how to fit two dipoles without symmetry constrains? I thought about it for a while but no worries if it's too much nuisance
Robert Oostenveld - 2016-06-02 09:23:21 +0200
(In reply to Diego Lozano Soldevilla from comment #9) it is possible to fit two dipoles without symmetry constraint, but not to do a grid search with two dipoles without symmetry. You'll have to position them at suitable starting locations. One strategy is to call ft_dipolefitting twice, in the first gridsearch and nonlinear with symmetry, and then in the second call to release the symmetry and further nonlinear optimize (without grid search, but using the result from the first). The low-level dipole_fit function does currently not support the fitting of one dipole while keeping another dipole at a fixed position. It would be possible to extend the dipole_fit->dipolemodel2param and param2dipolemodel function to support. Let me give that a try... mac011> git push origin bug3119-constr-sequential X11 forwarding request failed on channel 0 Counting objects: 4, done. Delta compression using up to 4 threads. Compressing objects: 100% (4/4), done. Writing objects: 100% (4/4), 859 bytes | 0 bytes/s, done. Total 4 (delta 3), reused 0 (delta 0) To git@github.com:robertoostenveld/fieldtrip.git dc689a0..dcc246c bug3119-constr-sequential -> bug3119-constr-sequential This branch (in robertoostenveld/fieldtrip) contains dipole_fit function with an extended constraint. Look in the code for "constr.nofit". You should be able to make your own "bug3119-constr-sequential" branch locally and to pull the changes on my branch into yours. Then you can try out. After fitting the first dipole, you should proceed with a 2nd call to ft_dipolefitting with cfg.dip.pos(1,:) = position from first fit cfg.dip.pos(2,:) = your best guess for the 2nd dipole cfg.dipfit.constr.nofit = [1 1 1 0 0 0] I still don't know how to deal with the interaction between the (simple) gridsearch, the (elaborate) nonlinear model and the strength/moment optimization at the end.