Back to the main page.
Bug 3037 - Inconsistent results with ft_dipolefitting
| Status | CLOSED FIXED |
| Reported | 2016-01-08 22:06:00 +0100 |
| Modified | 2016-06-14 16:14:52 +0200 |
| Product: | FieldTrip |
| Component: | core |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Windows |
| Importance: | P5 blocker |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
François Tadel - 2016-01-08 22:06:54 +0100
Created attachment 764 Example data and results I get wrong results from ft_dipolefitting after an update to a new version of FieldTrip: 20151001 gives good results, 20160107 gives wrong results. Attached: example of fitting a dipole 100ms after a simple auditory stim. - input.mat: contains the "cfg" and "data" variables in input - test_dipolefitting: illustrates how the code is executed - 20151001.jpg and 20160107.jpg: illustrates the results Is it the way to interpret the input structures that changed or is there a new bug somewhere?
Robert Oostenveld - 2016-01-09 12:11:25 +0100
thanks for reporting and the initial triage. First thing that comes to mind looking at the screenshots is that it is a unit mismatch. That might be a bug introduced not in the dipole fitting code, but elsewhere in the FT data handling sections. I'll have a look at the data as well.
Robert Oostenveld - 2016-01-09 12:38:45 +0100
(In reply to Robert Oostenveld from comment #1) >> cfg.grid ans = pos: [1524x3 double] inside: [1524x1 double] does not specify unit='m' (ideally it should), but it is consistent with grad and headmodel units. After grid search but before nonlinear optimization I get found minimum after scanning for topography 1 on grid point [-0.0724206 -0.021562 0.0484371] which is also an OK coordinate (in meter). After nonlinear I get >> dip1.dip ans = pos: [3.7980 0.2000 0.6490] mom: [3x1 double] pot: [274x1 double] rv: 0.9532 unit: 'm' which is consistent with your figure of the dipole being 3 meters away. But what strikes me is that the residual variance is nearly 100%. Also during nonlinear fitting it decreases only from 97 to 95 percent. So the fit is simply very bad. I notice you use a multiple spheres model. We don't use those in Nijmegen any more. So those parts of the code don't get used that often. If I do it with a single sphere like this cfg.headmodel = []; cfg.headmodel.r = 0.12 cfg.headmodel.o = [0 0 0.04]; cfg.headmodel.unit = 'm' dip2 = ft_dipolefitting(cfg, ftData) I get a RV of 55%. Still very bad.
Robert Oostenveld - 2016-01-09 13:07:37 +0100
Let's take a step back. with ft_multiplotER([], ftData) I see a clean ERF, which might be the result of auditory stimulation? The ERF is strong on the right side, much weaker (but not absent) on the left side. A single dipole does not seem like the most appropriate model, a dipole. Ah, messing around with the data and cfg I suddenly notice that you have cfg.channels: {1x274 cell} whereas it should be cfg.channel without the S. .... Nope, does not make a difference. Although it defaults to cfg.channel='all', the head model only has the 274, so it is still what you would want. Hmm, no luck identifying the cause of the difference yet. Summary so far: the data is perhaps not ideal for a single dipole, the localspheres model is not frequently used by us, the 2*4 coils per channel different from our default (although it is actually better than our 2 coil per channel default). The combination between the localspheres and the 2*4 coil representation of channels has never been explicitly considered when writing or updating the code, so that is a risk. But that has always been the case in FT code, and you write that it worked before, so I think that I should better be looking for something that is due to changes to FT rather than something that has always been the case. Ah, trying out a symmetric two-dipole model I just found a bug in a piece of code that did change over the last few months and that relates to localspheres. But I now have to join the family, shopping for the weekend... to be followed up.
François Tadel - 2016-01-11 18:23:11 +0100
Created attachment 765 Test data
François Tadel - 2016-01-11 18:24:40 +0100
Thanks for looking into this! I did the changes you suggested: - grid.unit = 'm' - .channel instead of .channels It's true that this bilateral response is not the best to test the single dipole fitting. I attached here a better example, a left median nerve stim, 30ms post-stim. I identified the dates of the change in behavior: - fieldtrip-lite-20151005: Works - fieldtrip-lite-20151006: Doesn't work
Robert Oostenveld - 2016-01-11 23:17:03 +0100
(In reply to François Tadel from comment #5) the double dipole should not be the solution to the bug (but would improve the interpretation of the data), but trying it (in the weekend) did give me an explicit error that I suspect to be pointing to the source of the problem. I am now tied up on a workshop in Marseille, but could you try the following (informal) patch to ft_compute_leadfield? mac011> svn diff forward/ Index: forward/ft_compute_leadfield.m =================================================================== --- forward/ft_compute_leadfield.m (revision 11049) +++ forward/ft_compute_leadfield.m (working copy) @@ -155,7 +155,7 @@ if isfield(headmodel, 'o') % shift dipole and magnetometers to origin of sphere - dippos = dippos - repmat(headmodel.o, Ndipoles, 1); + dippos = dippos - repmat(headmodel.o, Ndipoles, 1); coilpos = coilpos - repmat(headmodel.o, size(coilpos, 1), 1); end @@ -191,13 +191,13 @@ end lf = zeros(ncoils, 3*Ndipoles); - for chan=1:ncoils + for coil=1:ncoils for dip=1:Ndipoles % shift dipole and magnetometer coil to origin of sphere - dippos = dippos(dip, :) - headmodel.o(chan, :); - chnpos = sens.coilpos(chan, :) - headmodel.o(chan, :); - tmp = meg_leadfield1(dippos, chnpos, sens.coilori(chan, :)); - lf(chan, (3*dip-2):(3*dip)) = tmp; + tmppos = dippos(dip, :) - headmodel.o(coil, :); + coilpos = sens.coilpos(coil, :) - headmodel.o(coil, :); + tmp = meg_leadfield1(tmppos, coilpos, sens.coilori(coil, :)); + lf(coil, (3*dip-2):(3*dip)) = tmp; end end @@ -217,7 +217,7 @@ % tmp2 = 0.01*dippos'; %convert to cm % lf = megfield([tmp2 tmp2 tmp2], [[1 0 0]'*tmp1 [0 1 0]'*tmp1 [0 0 1]'*tmp1]); for dip=1:Ndipoles - R = 0.01*dippos(i, :)'; % convert from cm to m + R = 0.01*dippos(dip, :)'; % convert from cm to m Qx = [1 0 0]; Qy = [0 1 0]; Qz = [0 0 1];
François Tadel - 2016-01-12 01:55:18 +0100
Awesome! It looks like it's working! Thanks for the very quick update. Is this modification going to be part of fieldtrip-20160112.zip?
Robert Oostenveld - 2016-01-15 11:15:54 +0100
(In reply to François Tadel from comment #7) this is the fix. mac011> svn commit ft_compute_leadfield.m Sending ft_compute_leadfield.m Transmitting file data . Committed revision 11061. It will be in the FTP version this evening. I will also send an email to the list.
François Tadel - 2016-01-15 16:53:56 +0100
Thanks! I'll make sure it all works next week.