Back to the main page.
Bug 2387 - ft_prepare_leadfield does not correctly handle units
| Status | ASSIGNED |
| Reported | 2013-11-22 12:50:00 +0100 |
| Modified | 2015-03-17 09:10:08 +0100 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Windows |
| Importance: | P3 normal |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: | http://bugzilla.fcdonders.nl/show_bug.cgi?id=2853 |
Eelke Spaak - 2013-11-22 12:50:46 +0100
(making a bug out of it anyway to facilitate referencing this discussion from the code). My post on FieldTrip-dev mailing list: Dear FT dev team, As you may have read in my post on the general list, the test script test_tutorial_beamforming_extended no longer returns the expected results (i.e. plots are very different from http://fieldtrip.fcdonders.nl/tutorial/beamformingextended), *unless* one first does hdm = ft_convert_units(hdm, 'm'); sourcemodel = ft_convert_units(sourcemodel, 'm'); freq_cmb.grad = ft_convert_units(freq_cmb.grad, 'm'); before the call to ft_prepare_leadfield. (Converting all three to 'cm' also works fine, so I guess the problem is in whether the three agree or not.) Since the code used to work fine without this explicit conversion, I am considering this a bug. It might be related to bugzilla entries http://bugzilla.fcdonders.nl/show_bug.cgi?id=1794 and http://bugzilla.fcdonders.nl/show_bug.cgi?id=2265, although I have not read those two in detail. I am hesitant to file a new bug as it might be related to those two. For now, is it a good idea to add an error to ft_prepare_leadfield if the units for the inputs structures are not in agreement? We can later investigate why it used to work without this condition being met. I can add this error if people (especially JM) agree. Best, Eelke
Eelke Spaak - 2013-11-22 12:50:58 +0100
Robert's reply: agreed, it indeed seems a bug due to my recent work on the units. Please add the error, and I'll get to work on reproducing it and find a solution around it. Robert
Eelke Spaak - 2013-11-22 12:57:38 +0100
added the warning: bash-4.1$ svn commit ft_prepare_leadfield.m Sending ft_prepare_leadfield.m Transmitting file data . Committed revision 8835.
Eelke Spaak - 2013-11-22 13:15:28 +0100
For reference by users finding this bug: my educated guess is that this bug was introduced on 1 nov 2013 with revision 8695.
Robert Oostenveld - 2013-11-22 13:26:52 +0100
the problem is that ft_compute_leadfield (i.e. the low-level function, which expects the units to be consistent) is getting grad and vol in mm, but a position in cm. So the leadfields are computed for a very small region somewhere in the middle of the brain. This is clearly a bug in the handling of the geometrical units. The private/prepare_headmodel used to have if ~strcmp(vol.unit, sens.unit) % ensure that the geometrical units are the same warning('convering the sensor array units to "%s"', vol.unit); sens = ft_convert_units(sens, vol.unit); end to deal with the units of the sens and vol, but it ignores the units of the source model. In the past the units of the sourcemodel were specified in cfg.sourceunits and had a default of cm. That default is now removed, but the sourcemodel specifies its own units instead in grid.unit (like the sens and vol). I have changed it to the following % ensure that the geometrical units are the same if isfield(cfg, 'grid') && isfield(cfg.grid, 'unit') % convert it to the units of the source model sens = ft_convert_units(sens, cfg.grid.unit); vol = ft_convert_units(vol, cfg.grid.unit); else % convert it to the units of the head model sens = ft_convert_units(sens, vol.unit); end However, this is only a temporary fix. To ensure that the units of the forward and inverse computation are correct, Vladimir and I have been working on a set of changes. These involve that leadfields are only guaranteed to be correct if SI units are used throughout (i.e. meter, Volt, Tesla, Ampere, Siemens, and T/m, A*m, S/m etc.). Now that I think of it, the following is a better solution, it uses the same cfg.siunits option as ft_prepare_headmodel. if istrue(cfg.siunits) % ensure that the geometrical units are in SI units sens = ft_convert_units(sens, 'm'); vol = ft_convert_units(vol, 'm'); if isfield(cfg, 'grid') cfg.grid = ft_convert_units(cfg.grid, 'm'); end else % ensure that the geometrical units are the same if isfield(cfg, 'grid') && isfield(cfg.grid, 'unit') % convert it to the units of the source model sens = ft_convert_units(sens, cfg.grid.unit); vol = ft_convert_units(vol, cfg.grid.unit); else % convert it to the units of the head model sens = ft_convert_units(sens, vol.unit); end end Not that this solves everything, because ft_convert_units is only able to convert the geometrical aspects of the volume conduction model, not the conductivity and system matrix etc.