Back to the main page.
Bug 1317 - implement forward solution based on leadfield interpolation
| Status | CLOSED FIXED |
| Reported | 2012-02-08 13:13:00 +0100 |
| Modified | 2021-10-29 12:38:38 +0200 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P3 normal |
| Assigned to: | Vladimir Litvak |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
Robert Oostenveld - 2012-02-08 13:13:09 +0100
this is a desired feature for SPM12, and also useful for FT in general. TODO - It has been proposed in at least one, but if I recall correctly two published papers. I should look them up. - I have done this in the past (pre-FT time) and might have some relevant code lying around. I have to dig this up. - design a test script that demonstrates how it should work DONE - Christope sent me the files, they are 4D nifti files, one each for x, y, z. The position is coded inside the file. I will add them to FT test directory.
Robert Oostenveld - 2012-11-27 16:55:39 +0100
I have implementer the first farmework, including a test script. Please start with test_headmodel_interpolate and see how it is supposed to work. Concrete works still needs to be done in ft_prepare_vol_sens.m and in leadfield_interpolate.m. mbp> svn commit ft_voltype.m ft_compute_leadfield.m ft_prepare_vol_sens.m private/leadfield_interpolate.m ft_headmodel_interpolate.m ../test/test_headmodel_interpolate.m Sending forward/ft_compute_leadfield.m Adding forward/ft_headmodel_interpolate.m Sending forward/ft_prepare_vol_sens.m Sending forward/ft_voltype.m Adding forward/private/leadfield_interpolate.m Adding test/test_headmodel_interpolate.m Transmitting file data ...... Committed revision 7000.
Robert Oostenveld - 2012-11-30 16:50:30 +0100
I have made a first full implementation. It turned out to be quite difficult to get the montage handling right, i.e. the mapping between electrode contacts and channels (e.g. bipolar channels). This is specifically relevant for this method, because we want to use the more complex models for iEEG and ECoG, where the EEG default common-average-reference will not suffice. Still to do: - replace the interpolation (now nearest neighbor) - prevent the double handling of sens.tra in the interpolate section and the general section of ft_compute_leadfield - ensure that the test script (also updated) runs and add some (even more) challenging test cases mac001> svn commit forward/ test/ Sending forward/ft_apply_montage.m Sending forward/ft_headmodel_interpolate.m Sending forward/ft_prepare_vol_sens.m Sending forward/private/ft_datatype_headmodel.m Sending forward/private/leadfield_interpolate.m Sending test/test_headmodel_interpolate.m Transmitting file data ...... Committed revision 7069.
Robert Oostenveld - 2012-12-04 14:19:26 +0100
Moved the channel interpolation from ft_prepare_vol_sens to test_headmodel_interpolate, moved the nifti memmapping from ft_compute_leadfield to ft_prepare_vol_sens (i.e. both one level up). Updated the test script. mac001> svn commit test forward Sending forward/ft_compute_leadfield.m Sending forward/ft_headmodel_interpolate.m Sending forward/ft_prepare_vol_sens.m Sending forward/private/leadfield_interpolate.m Sending test/test_headmodel_interpolate.m Transmitting file data ..... Committed revision 7078. What still fails in the test script is the on-the-fly average referencing (or not). This needs a conceptual discussion. What also still needs to be done is the leadfield interpolation. At the moment it is (as a proof of concept) a nearest neighbor interpolation (i.e. simply return the leadfield of the nearest neighbor), but that should become a proper interpolation.
Robert Oostenveld - 2012-12-04 14:23:04 +0100
@Christophe and Johannes, Can you both have a look at the code, reading yourself into the code in the order 1) fieldtrip/forward/ft_headmodel_interpolate 2) fieldtrip/forward/ft_prepare_vol_sens 3) fieldtrip/forward/ft_compute_leadfield -> which calls private/leadfield_forward In leadfield_forward the source space interpolation should be done, please come with a suggestion. Should it be a tri-linear one, or one with a Gaussian smoothing?
Johannes Vorwerk - 2012-12-05 16:26:16 +0100
(In reply to comment #4) I once tested trilinear and higher order Bezier interpolation, which showed good results and is quite easy to implement. The most complicated part (especially for higher order interpolations) seemed to be the management of the grid underlying the interpolation.
Robert Oostenveld - 2012-12-05 17:25:36 +0100
if you have >> [x, y, z] = ndgrid(1:30, 1:40, 1:50); >> whos x y z Name Size Bytes Class Attributes x 30x40x50 480000 double y 30x40x50 480000 double z 30x40x50 480000 double and val = rand(30, 40, 50) or another scalar function that is defined on the 3-D grid, how do you assign for a given position that is inside the space spanned by the grid pos(1) = ceil(rand(1)*30); pos(2) = ceil(rand(1)*40); pos(3) = ceil(rand(1)*50); the non-zero elements of a weight matrix >> weight = zeros(30, 40, 50) such that the interpolated value at the desired position is >> vali = weight(:)' * val(:) i.e. the sum of the product of the elements.
Robert Oostenveld - 2013-05-24 14:11:52 +0200
@vladimir, please see forward/private/leadfield_interpolate.m forward/ft_headmodel_interpolate.m
Vladimir Litvak - 2013-06-18 11:52:44 +0200
Created attachment 487 1-r^2 for correlation between interpolated and bemcp lf
Vladimir Litvak - 2013-06-18 12:23:18 +0200
I did some work on interpolated lead fields. It seems to work but there are several outstanding issues: 1) @Robert, please check my changes to filetype_check_header and prepare_vol_sens and make sure they are OK with you. I remember that a long time ago you opposed allowing for vol to be a mat-file name in prepare_vol_sens but it seems to be the only way for me pass in the path where the lead fields are located. Since you put absolute paths in the vol struct the code won't work if you copy the nifti files to a different location unless you update those paths. That would be unfortunate as sharing pre-computed leadfields and downloading them from the web might be one of the main applications. If you have a different idea about how to achieve that, I'm open to it. 2) The code is currently quite slow. It takes tens of minutes to compute lead-fields for the SPM mesh 8000 vertices on my PC. This is better than running a FEM but still quite inconvenient. Most of the time is spent loading the nifti files into memory. The way things are implemented now the code must load all the channel files for every point. This is not very efficient. SPM interpolation allows passing all the positions at once but the code in ft_compute_leadfield does not support that at the moment. It would be much better if the SPM code could loop over channels rather than over vertices. 3) When comparing the mesh lead-fields computed by interpolation with those of the original bemcp I surprisingly found very weak correlations. I first thought that this was an implementation bug but investigating further I found that the correlation is only low close to the surface but very good in depth (see attached image). This is something we knew about in principle but I didn't think the problem was that severe. I now want to compute lead fields with OPENMEEG to see what happens there but that will take a while to run so I'm checking in the code for you to look at without waiting for it.
Robert Oostenveld - 2013-06-18 14:23:48 +0200
(In reply to comment #9) thanks for your work. I also got your text message. At the moment I am in Aspet for the EEGLAB workshop with rather crappy wifi, but I'll give it a try to work on it. Re 1: good point. I did not put any thought into absolute or relative filenames yet, but it should be relative. Re 2: looping over positions is supported in ft_compute_leadfield, see at line 117. However, I am not sure whether it is actually supported for all the voltypes. I see that in quite some cases it should be dealt with at a lower level... I checked, private/leadfield_interpolate indeed ignores it. Needs to be fixed. Re 2: I consider looping over channels not a good idea. It would interfere with the on-the-fly referencing/montage which ft_compute_leadfield does towards the end. Re 3: the results look nice. It getting more inaccurate at the boundaries is not surprising, the leadfields change more rapidly there. We have to find out a decent resolution that is not too slow when precomputing , and not too inaccurate when interpolating. I will check your code and get back to you.
Vladimir Litvak - 2013-06-18 15:30:42 +0200
(In reply to comment #10) OK. Looping over channels is not a must. It's just to have the interface update and not freeze. My code in leadfield_interpolate already supports multiple points at once so if you fix your part it should work.
Robert Oostenveld - 2013-06-18 15:53:10 +0200
enhancement - allow reading of any file, not only of a mat file. mbp> svn commit forward/ft_prepare_vol_sens.m Sending forward/ft_prepare_vol_sens.m Transmitting file data . Committed revision 8273.
Robert Oostenveld - 2013-06-18 15:55:04 +0200
enhancement - replaced 0/1 with false/true, added some documentation. Functionality did not change. mbp> svn commit fileio/private/filetype_check_header.m Sending fileio/private/filetype_check_header.m Transmitting file data . Committed revision 8274.
Vladimir Litvak - 2013-07-01 13:16:49 +0200
Created attachment 493 1-r^2 for correlation between interpolated and openmeeg lf
Vladimir Litvak - 2013-07-01 13:18:23 +0200
(In reply to comment #14) I tried the same thing for OPENMEEG and it's not much better (see the second attachment). I'm quite puzzled by this. The grid resolution was 5mm. Do you think it might be insufficient? We should probably get to the bottom of this before we roll it out. Vladimir
Vladimir Litvak - 2013-07-30 16:18:27 +0200
(In reply to comment #15) Also see this unfinished issue. There are two things here: 1) Changing ft_compute_leadfield to consistently handle multiple points at once. 2) Understanding why there are differences between interpolated and full computation close to the surface.
Vladimir Litvak - 2013-10-29 16:43:49 +0100
(In reply to comment #16) I tried again to compare original to interpolated leadfields with BEMCP in a 2x2 cm cube with 1mm spacing. Still for voxels around the boundary there are low correlations. So it might be the case that it's just a BEM problem and if we use a FEM we'd be able to sample the lead fields much coarser. In any case BEM gives rubbish values there so there is no point to sample them better. You can try with dipoli and see what you get. Vladimir
Robert Oostenveld - 2013-10-29 17:25:38 +0100
BEM is known to become inaccurate when the dipoles get "close" to the surface, where close is defined as "in the order of 1/2 to 1 edge-length". In the comparison of BEM versus interpolated-BEM that should ideally not matter, but if the BEM leadfields start behaving erratically, that might be a reason for the interpolation to get into trouble.
Vladimir Litvak - 2013-10-30 14:04:27 +0100
(In reply to comment #18) I repeated the comparison for 3-sphere model and it looks like things only break down at the boundary so inside the brain everything looks OK and there is no bug. My evaluation method is somewhat coarse so if someone does it more carefully everything should look even better. I also looked at the HBM paper from 2001 which seems to be the only previous paper about lead field interpolation. They also encountered this problems with BEMs and used 'truncated grids' to get around it. They also describe some comparisons of grid resolutions (they got 8mm as a good compromise) and interpolation methods. However, they did not interpolate on the scalp so this is a new development that should probably be written up. Also they didn't look at FEMs.
Vladimir Litvak - 2013-10-30 14:06:39 +0100
(In reply to comment #19) Given that the problem with points close to the surface is quite severe should we perhaps also have some code to automatically correct grids and meshes to make sure they don't suffer from BEM errors?
Vladimir Litvak - 2013-11-06 15:55:14 +0100
(In reply to comment #20) I added code to ft_headmodel_interpolate to perform extrapolation by smoothing and normalizing by smoothed mask (this is a trick used in SPM that Guillaume showed me). At the moment this uses SPM smoothing functions and creates and deletes temporary nifti files, which is slightly inelegant (although the code depends on SPM anyway). If you want to implement the same thing cleaner, please do. I now convinced myself that IF everything is done carefully on both interpolated and original BEM sides then there is good match between original BEM and interpolation. However, that IF is not trivial and requires using the right inwardshift and also removing/correcting mesh vertices that stick out. I empirically determined the critical inwardshift for BEMCP and SPM meshes and it turns out to be 6mm which is the same as the average triangle side length (but not half of that as in some references).
Robert Oostenveld - 2013-11-06 16:03:36 +0100
(In reply to comment #21) I think the "1/2 triangle length" depends on the BEM details. BEMCP is not the most sophisticated implementation. But it is good to know that the rule-of-thumb is to be used with caution.
Vladimir Litvak - 2013-11-06 16:08:43 +0100
(In reply to comment #22) When do you think you will finish implementing everything we planned last week?
Vladimir Litvak - 2013-11-08 13:01:20 +0100
(In reply to comment #23) Just to note some things I thought of when writing the paper: 1) It might be good to save the extra information for the volume not as a .mat file but as something that can be easily read by other software because we want to sell this as a way to establish a repository of lead-fields for any software not just Matlab-based. One easy possibility is to save it as an .m file that when run generates the vol. SPM has some code for generating code from structs, not sure about FT. 2) The way the code is written now is that the spline coefficients are computed from the beginning in ft_headmodel_inteprolate. As those coefficients can be computed just from the lead-field volumes for distributing lead-fields via the web it'd be better not to compute them at that stage but if ft_prepare_vol_sens. That'd mean that something would have to be done even if there is match between sensors. 3) There should be an option to say where to store the re-interpolated vol. Some users might have problems with the temp directory or they want to keep it for subsequent use.
Robert Oostenveld - 2013-11-08 17:11:50 +0100
(In reply to comment #24) Robert will do 1, Vladimir will do 2. 3 is not needed
Robert Oostenveld - 2013-11-26 23:00:44 +0100
to do: consider the units, see http://bugzilla.fcdonders.nl/show_bug.cgi?id=2377#c27
Robert Oostenveld - 2013-12-09 14:24:12 +0100
I revisited the test_headmodel_interpolate script. I found an issue with file naming, inside ft_headmodel_interpolate it was assumed that files with name xxx were located in directory xxx. However, that was not consistently implemented. Furthermore, temporary files were created in a location that (potentially) did not have write permission. I changed it so that teh filename passed to ft_headmodel_interpolate is a file name. It can include the directory, and in that directory (or in pwd) the temp files are written. So if you want xxx/xxx_elec1.nii, then you have to specify xxx/xxx as filename. Another issue is that the leadfields are not exactly the same any more. I suspect that to be due to the smoothing that is now implemented. I wonder whether that should be always applied, or whether it should be optimal. The electrode close to the source now gets a mush smaller potential (as there the potential changes most rapidly as function of dipole position). For now I have added an option "smooth" with default=true. When false, the script runs through as expected (i.e. leadfields are exactly the same). mac001> svn commit Sending forward/ft_compute_leadfield.m Sending forward/ft_headmodel_interpolate.m Sending forward/ft_prepare_vol_sens.m Sending forward/private/leadfield_interpolate.m Sending test/test_headmodel_interpolate.m Sending utilities/ft_datatype.m Transmitting file data ...... Committed revision 8997. The present code does not deal correctly with the referencing of the EEG, i.e. channel selections messes up the referencing.
Robert Oostenveld - 2021-10-29 12:30:28 +0200
I came across this thread while reviewing the "interpolate" head model in relation to https://www.fieldtriptoolbox.org/development/project/femfuns/. I thin that the interpolate head model can be reused for FEMFUNS. There was a test script "failed_headmodel_interpolate" which, as its name suggested, failed to run. The errors were relatively minor, partially to some bookkeeping, and partially to non-standard reference schemes in elec.tra. The commits from yesterday https://github.com/fieldtrip/fieldtrip/commit/bd8e8d58a55b393c8043647ac3cb951c0b1ec9dc and https://github.com/fieldtrip/fieldtrip/commit/0c2c4df798ca674e6eecf4d1f13353879b8e336b address the bookkeeping updates. It turned out that fixing the reference handling was also not too hard. Rather than initializing the elec.tra as eye(N), it should be initialized as eye(N)-1/N, since in thaht case it corresponds to an average reference forward solution. In case elec.tra is not present, ft_compute_leadfield will also compute the average reference over the lead field. See https://github.com/fieldtrip/fieldtrip/pull/1909 for the pull request.
Robert Oostenveld - 2021-10-29 12:35:57 +0200
as far as I can tell, this implementation is complete and correct now that https://github.com/fieldtrip/fieldtrip/pull/1909 has been merged. Hence I will close this thread.