Back to the main page.
Bug 1794 - ft_prepare_leadfield output discrepancy between cm and mm units for predefined grid
| Status | CLOSED WORKSFORME |
| Reported | 2012-10-26 21:44:00 +0200 |
| Modified | 2019-08-10 12:41:36 +0200 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Linux |
| Importance: | P1 major |
| Assigned to: | Johanna |
| URL: | |
| Tags: | |
| Depends on: | 2355 |
| Blocks: | 2265 |
| See also: | http://bugzilla.fcdonders.nl/show_bug.cgi?id=686http://bugzilla.fcdonders.nl/show_bug.cgi?id=2817 |
Giorgos Michalareas - 2012-10-26 21:44:06 +0200
Created attachment 354 mat file containing the cfg to be input in ft_prepare_leadfield.m The main issue is that if a predefined standard grid is used, then if the grid positions are defined in MM, then the output leadfields have huge values. At a first glance it seems that volume and sensor structures are converted to CM units while the precomputed grid units are not and although they are MM they are treated as CM (not sure about this last part). If the grid units are converted to CM along with the volume and the sensor and then input to ft_prepare_leadfield, then the leadfields look much more realistic. So I guess the first issue is that ft_prepare_leadfield should produce the same leadfields independent of the units in the input grid, volume and sensor information. The second issue is that I compared the leadfields produced by the latest version of FT (20121025) with an older version of fieldtrip, namely 20120426. The only difference in the input was that when using the latest FT the grid,volume and sensor units were CM , while when using the old FT version the units were MM. %=================================== load bugcfg cfg; % find attached file bugcfg containing a cfg structure with a precomputed grid, ready for input in ft_prepare_leadfield restoredefaultpath; % Here Add fieltrip-20120426 directory ft_defaults; gridLFold= ft_prepare_leadfield(cfg); restoredefaultpath; % Here Add fieltrip-20121025 directory gridLFnew= ft_prepare_leadfield(cfg); %=================================== By comparing the leadfields, I observed that the leadfields were differing by 2 orders of magnitude. gridLFold.leadfield{gridLFold.inside(1)}(1,:)= [4.325760634513e-14 4.46333595600794e-13 6.86774945338438e-14]; gridLFnew.leadfield{gridLFnew.inside(1)}(1,:)= 4.32576063451312e-12 4.46333595600792e-11 6.86774945338437e-12 This is very strange as the leadfields (and consequently projected signal magnitude and power) should not be affected by the units used. So the second issue I guess is if the latest version computes leadfields correctly and if yes then does this mean that the leadfields from older versions were not correct? And if this is the case which fieldtrip versions had this problem? I have already produced results based on leadfields so I d like to know which of my results I need to recompute. Thank you Best
Giorgos Michalareas - 2012-10-26 21:54:46 +0200
(In reply to comment #0) I just post here again the code I posted in the initial post, but with the conversion from mm to cm before the use of the latest FT version (forgot to put it in the initial post - sorry). %=================================== load bugcfg cfg; % find attached file bugcfg containing a cfg structure with a precomputed grid, ready for input in ft_prepare_leadfield restoredefaultpath; % Here Add fieltrip-20120426 directory ft_defaults; gridLFold= ft_prepare_leadfield(cfg); restoredefaultpath; % Here Add fieltrip-20121025 directory ft_defaults; cfg.grid=ft_convert_units(cfg.grid,'cm'); cfg.vol=ft_convert_units(cfg.vol,'cm'); cfg.grad=ft_convert_units(cfg.grad,'cm'); gridLFnew= ft_prepare_leadfield(cfg); %===================================
Johanna - 2013-01-23 13:17:12 +0100
I can take a look.
Johanna - 2013-01-23 17:35:23 +0100
I have replicated that the results from a cfg with fields in mm gives a different answer than a cfg with fields in cm. however, the results are the same between a verison 20120630 and current (20130123) as those were two easy versions I had on hand. I will look further later.
Johanna - 2013-01-25 15:27:34 +0100
Sticking with the current svn version for the moment (see Comment 3 why), it occurs because ft_compute_leadfield only knows about sens.unit and vol.unit (and does indeed ensure they are the same) but not the units of pos which are entered in as a vector, not as a structure with units. Any attempt to discover their units with grid.pos=pos; ft_convert_units(grid); fails since it doesn't work well on a 1x3 size vector of .pos. Furthermore, ft_compute_leadfield seems to assume units for pos. Not clear to me which units, but I suspect 'cm'. It gets worse if I try giving mix/match combinations of vol.unit, sens.unit and grid.unit. Then the grid that gets computed in line 174 of ft_prepare_leadfield by the call to ft_prepare_sourcemodel has units that don't match its .pos, if the units of vol(sens).unit are different. I commit test_bug1794 which shows these scenarios. I can think to suggest 3 options: 1) Force ft_prepare_leadfield to convert grid to units that ft_compute_leadfield assumes prior to calling ft_compute_leadfield. If so, do vol.unit and sens.unit also need to be this same? I think so. 2) Force an obligatory 4th input argument to ft_compute_leadfield with known units of pos. (and then convert within the function, using this info). 3) Force the 1st input argument to ft_compute_leadfield to be a structure with grid.pos (1x3) and with grid.unit. I feel like the latter two are safer for ensuring that inputs to ft_compute_leadfield are correct irrespective of how this function is called, but are worse in terms of backwards compatibility. At a minimum, make explicit in the 'help' of ft_compute_leadfield what the units should be for pos, sens, vol.
Robert Oostenveld - 2013-01-28 12:14:50 +0100
There is the (old) option cfg.sourceunits that we used in the past to distinguish between source and mri units. mac001> grep -l sourceunits *.m ft_dipolefitting.m ft_megrealign.m ft_prepare_headmodel.m ft_prepare_leadfield.m ft_prepare_localspheres.m ft_prepare_mesh.m ft_prepare_sourcemodel.m ft_sourceanalysis.m ft_sourceinterpolate.m There also used to be a cfg.mriunits, but that is now forbidden/deprecated because those are always present or can be inferred from the mri. mac001> grep mriunits *.m ft_prepare_mesh.m:cfg = ft_checkconfig(cfg, 'forbidden', {'numcompartments', 'outputfile', 'sourceunits', 'mriunits'}); ft_prepare_singleshell.m:cfg = ft_checkconfig(cfg, 'deprecated', 'mriunits'); ft_prepare_sourcemodel.m:cfg = ft_checkconfig(cfg, 'deprecated', 'mriunits'); ft_sourceinterpolate.m:cfg = ft_checkconfig(cfg, 'deprecated', {'sourceunits', 'mriunits'}); I suggest to rethink how cfg.sourceunits is supposed to behave and whether all required high-level ft functions use it accordingly.
Johanna - 2013-02-07 13:25:16 +0100
I cc Vladimir, as I saw his comment in another bug on lead field units.
Johanna - 2013-04-23 16:35:00 +0200
Regarding the difference between grid.leadfield on the order of units (meaning, a factor of 10 or 100): If I fix prepare_headmodel so that cfg.sourceunit = 'cm' is default (rather than only set to 'cm' if vol and sens units are different), then it works that the grid.leadfield will be (*nearly) the same irrespective of units of grid, sens and vol. So that seems an initial good solution to me. I will commit this change today. However, I have traced a few problems to vol.forwpar, which is created inside ft_prepare_vol_sens (called by prepare_headmodel, which is called by ft_prepare_leadfield). Firstly, vol.forwpar will be in units of vol as it is passed into ft_prepare_vol_sens, but then this later cannot be changed with ft_convert_units (i.e. if a later higher level function receiving this vol wants to change units). Perhaps this is a bug of ft_convert_units? (not able to convert vol.forwpar ?) Secondly, the (*nearly) refers to some remaining numerical imprecision (on the order of 10*eps) bewteen grid.leadfield of inputs of different units, which can be traced to imprecision within vol.forwpar. Suggestions what to do with vol.forwpar, regarding both its instability at the level of eps, and regarding it preventing the units of vol to be changed later on?
Robert Oostenveld - 2013-04-23 17:45:36 +0200
(In reply to comment #7) thanks for looking into this. ft_convert_units was only designed to deal with geometrical objects, not with parameters related to the forward model. forwpar is for single shell, and should only be present after ft_prepare_vol_sens. Other cases of forward model parameters are for bem and fem models. I think you have identified an issue with at which moment which conversion can be made. Geometrical conversions are expected to happen prior preparing the model parameters, but now happen after. We should rethink as to which conversions can happen when. Some are harmless, like converting T to fT, but some are harmful, such as T/cm to T/mm. ft_compute_leadfield now attempts to convert with % 'unit' = string, can be 'arbitrary' or 'si' (default = 'arbitrary')
Robert Oostenveld - 2013-11-01 14:33:41 +0100
I have just now made changes to the code to manage units more consistently, in the future allowing for explicit units throughout all analysis. I have seen this bug that seems to be related, and will go through it once to see whether it reproduces and whether it is still an issue. If I cannot reproduce it, I will close it.
Robert Oostenveld - 2013-11-01 14:48:37 +0100
the initial assumption by Giorgos about the leadfields not changing if you change the units is not fully correct. If the distance between the source and sensor is 10x smaller, the field will be stronger. Of course if the leadfield is computed with knowledge about the units, and if the leadfield units (T/Am) is considered, then the leadfield should not change, just like 10cm being equal to 100mm. But up to now FT is not so smart that it knows the units throughout the code, hence 10 and 100 are different. Right now the existing script fails with Error using test_bug1794 (line 68) Assertion failed. on assert(isequalwithequalnans(rmfield(gridLFmmcur,'cfg'),rmfield(gridLFcmcur,'cfg'))) which is due to >> gridLFmmcur gridLFmmcur = pos: [3094x3 double] unit: 'mm' dim: [17 14 13] xgrid: [-5.5000 -4.5000 -3.5000 -2.5000 -1.5000 -0.5000 0.5000 1.5000 2.5000 3.5000 4.5000 5.5000 6.5000 7.5000 8.5000 9.5000 10.5000] ygrid: [-6.8000 -5.8000 -4.8000 -3.8000 -2.8000 -1.8000 -0.8000 0.2000 1.2000 2.2000 3.2000 4.2000 5.2000 6.2000] zgrid: [-0.6000 0.4000 1.4000 2.4000 3.4000 4.4000 5.4000 6.4000 7.4000 8.4000 9.4000 10.4000 11.4000] inside: [1x1488 double] outside: [1x1606 double] cfg: [1x1 struct] leadfield: {1x3094 cell} >> gridLFcmcur gridLFcmcur = pos: [3094x3 double] unit: 'cm' dim: [17 14 13] xgrid: [-5.5000 -4.5000 -3.5000 -2.5000 -1.5000 -0.5000 0.5000 1.5000 2.5000 3.5000 4.5000 5.5000 6.5000 7.5000 8.5000 9.5000 10.5000] ygrid: [-6.8000 -5.8000 -4.8000 -3.8000 -2.8000 -1.8000 -0.8000 0.2000 1.2000 2.2000 3.2000 4.2000 5.2000 6.2000] zgrid: [-0.6000 0.4000 1.4000 2.4000 3.4000 4.4000 5.4000 6.4000 7.4000 8.4000 9.4000 10.4000 11.4000] inside: [1x1488 double] outside: [1x1606 double] cfg: [1x1 struct] leadfield: {1x3094 cell} which are different because one is in cm and the other in mm. So according to me the assertion should have failed, as it does: If you numerically compare 10cm and 100mm, them 10~=100. So this mostly seems to be an inconsistency between the users' (Giorgos') expectation and what fFT does, not a bug in the code. Please see http://code.google.com/p/fieldtrip/source/detail?r=8691 and http://code.google.com/p/fieldtrip/source/detail?r=8695 which are two changes that relate to this. The test case implemented in the test_bug1794 script is interesting however, as it goes beyond ft_prepare_headmodel all the way to ft_prepare_leadfield. I will revitis this at a later time, but now have to leave.
Eelke Spaak - 2015-02-11 15:41:15 +0100
I ran into this bug's failing test script today at the bug binge. The test script doesn't really makes sense to me, perhaps it should be revisited at some point?
Johanna - 2015-02-11 16:11:26 +0100
The failure now seems to be something unrelated to the main point of the bug, but with ft_preamble. But if you have any questions on the purpose of test_bug1794, maybe I can help...?
Robert Oostenveld - 2015-02-12 10:43:23 +0100
(In reply to Johanna from comment #12) I just ran the test script, it was failing for me due to the change of the indexed inside representation (1, 2, 3, 4, ...) into the boolean representation, as I announced on the mailing list. I changed gridLFcmcur.inside(1) into find(gridLFcmcur.inside, 1) to get the index of first inside source position. Subsequently it runs further, but fails on assert(isequaln(rmfield(gridLFmmcur,'cfg'),rmfield(gridLFcmcur,'cfg'))) where the fro grids are obviously different because one is in mm and the other in cm. mac011> svn commit Sending test/test_bug1794.m Transmitting file data . Committed revision 10227.
Robert Oostenveld - 2017-01-10 15:14:37 +0100
this is an old issue and can be closed. If the input is according to SI, the output is also according to SI. This is being tested with - test_meg_leadfield_units.m - test_eeg_leadfield_units.m - test_neuromag_units.m See also http://www.fieldtriptoolbox.org/faq/units If there are concerns about automatic conversion between units, please post an example script including ALL input data that demonstrates the problem.
Robert Oostenveld - 2019-08-10 12:35:24 +0200
This closes a whole series of bugs that have been resolved (either FIXED/WONTFIX/INVALID) for quite some time. If you disagree, please file a new issue on https://github.com/fieldtrip/fieldtrip/issues.