Back to the main page.

Bug 1967 - ft_prepare_vol_sens needs to be updated for simbio

Reported 2013-02-01 17:03:00 +0100
Modified 2021-10-29 12:38:34 +0200
Product: FieldTrip
Component: forward
Version: unspecified
Hardware: PC
Operating System: Windows
Importance: P3 normal
Assigned to: Robert Oostenveld
Depends on:
See also:

Lilla Magyari - 2013-02-01 17:03:32 +0100

What is happening in ft_prepare_vol_sens in case 'simbio'? The code also contains reference to 'bnd' and 'tri' but the simbio headmodel is not triangulated.

Robert Oostenveld - 2013-04-06 16:30:58 +0200

I have made a test script and resolved the issue. The skin surface is now extracted (either from the tet or hex) and for each electrode the distance to all vertices is determined. The electrode is then replaced by the vertex that is nearest. This is consistent with how it is done in the lower level simbio code, exept that there no outer surface assumption is made. I have just committed the code through github. Still to be done is to give some information about the projection, e.g. "the electrode were shifted to the skin surface average displacement %d mm". This should also be added to the BEM methods. Furthermore it might be desireable to allow the surfaceprojection to be switched on/off, i.e. make it optional (and default=yes)

Johannes Vorwerk - 2013-05-22 15:58:12 +0200

(In reply to comment #1) Should a call to ft_prepare_vols_sens be included somewhere in the test_headmodel_simbio script since the electrode positions are created outside the mesh? I couldn't find a call to this.

Johannes Vorwerk - 2013-05-23 12:18:49 +0200

(In reply to comment #2) Ok, just missed it. Forget about my comment...

Johannes Vorwerk - 2013-08-06 10:50:01 +0200

I think it would be reasonable to include a check whether an eeg transfer matrix is already computed in ft_prepare_vol_sens, my idea would be something like @@ -502,7 +502,12 @@ sens = rmfield(sens, 'chanpos'); end - vol.transfer = sb_transfer(vol,sens); + if ~(isfield(vol, 'transfer') && isfield(vol, 'elecpos') && (sens.elecpos == vol.elecpos)) + vol.transfer = sb_transfer(vol,sens); + vol.elecpos = sens.elecpos; + else + fprintf('using precomputed transfer matrix\n'); + end case 'interpolate' This would also ensure that the transfer matrix fits to the electrode positions given to ft_prepare_vol_sens, even though i am not sure if sens.elecpos == vol.elecpos works or if one should check whether the distance between the two sets is below a certain limit?

Robert Oostenveld - 2013-08-06 11:26:11 +0200

(In reply to comment #4) so you are suggesting that ft_prepare_vol_sens should behave differently if you call it the second time compared to the first time. Right now I think that the behaviour of the function is not defined when you call it a second time, i.e. in the way that the function is documented and used in fieldtrip, spm and eeglab it is never called twice. Note that it supports ~10 headmodels, and the behaviour of the function should remain consistent for all. In the simbio case: where does the transfer matrix come from? Ah, I see line 505: vol.transfer = sb_transfer(vol,sens); so that means that the transfer depends on the sens. What should happen if you call it like [vol1, sens1a] = ft_prepare_vol_sens(vol, sens1) [vol2, sens2a] = ft_prepare_vol_sens(vol1, sens2) where sens1 is completely different from sens2?

Johannes Vorwerk - 2013-08-06 11:37:57 +0200

(In reply to comment #5) my idea is not to recompute things that are already included in the headmodel. if the sensors do not change (assured by sens.elecpos == vol.elecpos, possibly a different way to check this would be more optimal), the transfer matrix stays the same so that it would be an unnecessary and long (setting up the transfer matrix is the most time consuming part!) recomputation. note that this change is in a part that is only used for a simbio calculation and that the output of the routine remains unchanged.

Robert Oostenveld - 2013-08-06 11:52:16 +0200

(In reply to comment #6) > note that this change is in a part that is only used for a simbio calculation > and that the output of the routine remains unchanged. Allowing the function to be called twice is a change that also affects other headmodels, as it would have to be reflected in higher level code (e.g. SPM). So the considerations not only apply to simbio, but also openmeeg, dipoli and other models. If higher level code change and the user is able to make a choice, then the behaviour for the second call has to be defined for all user-choices, hence also for all other head models. Right now that behaviour is not defined. But it does make sense to consider how to avoid the lengthy computation from happening twice. An alternative place to implement in the overall structure this is one level up, i.e. not in the forward functions. I do suspect that to be more problematic, though.

Johannes Vorwerk - 2013-08-06 12:29:11 +0200

(In reply to comment #7) i don't really get why it should not be possible to call the routine a second time for other dipole models? it will just overwrite the previous output. the point is that it will usually not happen, since the computational effort is very low for these approaches, while in the fem case it might be reasonable to call ft_prepare_vol_sens already prior to ft_prepare_leadfield to be able to easily save the result.

Robert Oostenveld - 2013-08-06 12:37:37 +0200

(In reply to comment #8) "dipole models" -> head models. The function has not been written with calling it twice in mind. Hence the behaviour is undefined. I expect for example the bookkeeping for multisphere MEG models with a mixture of magnetometers and gradiometers in the sens structure to get confused if you call it twice. The expected behaviour there is that it gets as input one channel per channel, and outputs one channel per coil (two coils per channel). Also for the EEG BEM methods there is some channel specific handling (computing of an surface interpolation matrix) that I expect to be problematic if you were to call it twice. The function behaviour being "undefined" does not mean that it will always fail if you call it twice. It just means that there are cases where it might fail. And if I change the high level code allowing the function to be called twice, then those cases would surface as bugs. If you want to work on it, I suggest that you go through the list ft_headmodel_asa.m ft_headmodel_bemcp.m ft_headmodel_concentricspheres.m ft_headmodel_dipoli.m ft_headmodel_fns.m ft_headmodel_halfspace.m ft_headmodel_infinite.m ft_headmodel_interpolate.m ft_headmodel_localspheres.m ft_headmodel_openmeeg.m ft_headmodel_simbio.m ft_headmodel_singleshell.m ft_headmodel_singlesphere.m ft_headmodel_slab.m and for each one check whether the function behaves on the second call as you would expect it. Or as alternative you might look at the higher level code and try to find a solution there.

Johannes Vorwerk - 2013-08-06 12:45:22 +0200

(In reply to comment #9) my plan was not to change anything in the higher level code. as i wrote my plan was to change one part in ft_prepare_vol_sens. afterwards, one could tell (fem) users that they have the possibility to call ft_prepare_vol_sens prior to calling ft_prepare_leadfield to be able to store the result of the transfer matrix computation. of course, this advice would have to be marked as experimental/not recommended and fem only.

Robert Oostenveld - 2013-08-06 13:07:31 +0200

(In reply to comment #10) at this moment ft_prepare_vol_sens will not and cannot be called twice from software such as fieldtrip, eeglab and spm, so it would not make sense to only change the section of code for simbio. Either the function always works as expected, and the code that makes use of it is updated to reflect this, or the code stays as it is (i.e. not supporting two calls).

Johannes Vorwerk - 2013-08-06 13:18:40 +0200

i don't see why it cannot be called twice? right now it can be called twice by any user willing to check the construction of the sensors. it is a seperate (not private!) function which the user can call anytime he/she wants to and it would offer a very simple opportunity to make the fem-pipeline more convenient without affecting any of the other pipelines.

Robert Oostenveld - 2013-08-06 14:01:05 +0200

(In reply to comment #12) > i don't see why it cannot be called twice Because the functions have a stable API and that API does not specify calling the function twice. Please read the design documents at If you want to change or extend the API (i.,e. the way that other code calls the functions from the forward toolbox, such as high-level fieldtrip, eeglab or spm makes use of it), then the consequences need to be thought through. Right now one consequence would be that bugs might be introduced if high level code calls the function twice with a localspheres or openmeeg vol as input. An API is the consequence of an agreement on how functions should behave, not a consequence of what technically is possible. Technically it is indeed possible to change how functions behave, but it being possible does not make it automatically desirable to make ad-hoc changes. In which end-user scenario would presently ft_prepare_vol_sens be called twice? Would changing the ft_prepare_vol_sens behaviour be an improvement for that specific scenario? If at this moment it is never called twice, then changing the function behaviour would not make a difference.

Johannes Vorwerk - 2013-08-06 14:42:53 +0200

(In reply to comment #13) scenario: user calls ft_prepare_vol_sens to check meg sensor arrangement / channel deselection / eeg electrode projection manually. afterwards user applies ft_prepare_leadfield.

Johannes Vorwerk - 2013-08-06 14:46:02 +0200

(In reply to comment #9) there will be no second call to ft_prepare_vol_sens in these scenarios, thus the function will behave as is. the only scenario i am thinking of is a user that wants to use the same transfer matrix several times (which is very probable) and wants to store it therefore. this has no effect at the moment, since a new one will be computed regardless of one being already present or not.

Robert Oostenveld - 2013-08-06 15:21:49 +0200

(In reply to comment #14) > scenario: user calls ft_prepare_vol_sens to ... the user does never call ft_prepare_vol_sens, as that is low-level code (like all other code in fieldtrip/forward). The user calls ft_prepare_mesh, ft_prepare_headmodel, ft_prepare_leadfield, ft_sourceanalysis. These functions have a non-graphical user interface (UI), whereas forward has an application programming interface (API). Did you try out the tutorial? it might help to get to the end-user perspective.

Johannes Vorwerk - 2013-08-06 15:33:29 +0200

(In reply to comment #16) that may be true for the absolute beginner that follows the tutorial, when i use fieldtrip, i nearly exclusively used low-level functions. also i am surely aware of the "regular" pipeline. nevertheless, i spoke to some people that would be interested to use this pipeline and are clearly above beginner level. none of those would be interested to use it if they are not able to store the transfer matrix. since this is nearly impossible to include in the current high level functions, i would like to introduce a workaround for these users and possibly shortly mention it in the tutorial. as i said previously, this will not change any of the high level functions, the regular workflow will just stay the same. BUT: we will offer the more experience user a convenient way to save a lot of time, surely at his own risk. an alternative would be, to create an additional function that has to be called prior to ft_prepare_leadfield for fem-users, where the transfer matrix is calculated.

Robert Oostenveld - 2013-08-06 16:09:17 +0200

(In reply to comment #17) People that want to follow the pipeline still should use the documentation. Or they can hack their own version of the code together, as it is all open source. If you want to continue working towards integrating simbio with fieldtrip and spm and not just cobble some code together, then please consider that there structure in the code that extends beyond your own needs. I am not opposing the idea of improving this code, I just mentioned that you have to keep in mind that other pieces of the code should not break as a consequence thereof. Rather than providing workarounds, why don't we just finish this project? If finishing it requires that ft_prepare_vol_sens can be called twice, then that is fine with me. As long as it does not introduce bugs elsewhere. However the suggested change to the code seems to introduce bugs, therefore it qualifies as a hack rather than a work-around.

Johannes Vorwerk - 2013-08-06 21:09:35 +0200

(In reply to comment #18) ok, probably my thinking this afternoon was a bit too complicated. my goal was mainly to give the user an easy option not to compute the transfer matrix over and over again, which seemed to be just to let him call ft_prepare_vol_sens by hand previously to ft_prepare_leadfield in the fem case and save the result. to prevent confusion for users that might accidentally do this in other (not tested) cases, how about introducing a function that mainly mimics the behaviour of ft_prepare_vol_sens, but works only for the fem part. the user can then optionally precompile a transfer matrix and save the result. additionally we add the above code, to check whether a transfer matrix was previously calculated for the fitting set of sensors. an alternative might be to add an option in ft_prepare_leadfield to save the result of the transfer matrix calculation, but this seems to be more complicated to me.

Robert Oostenveld - 2013-08-07 09:22:53 +0200

(In reply to comment #19) I think that we should look at what is presently preventing people from saving it. if I do [vol2, sens2] = ft_prepare_vol_sens(vol1, sens1) then I can simply save vol2 to disk and reuse it the next day (not putting it through ft_prepare_vol_sens again) and pass it onto ft_compute_leadfield. But if I do the high-level approach 1 ft_prepare_mesh 2 ft_prepare_headmodel 3 ft_compute_leadfield 4 ft_sourceanalysis I can save the output of each step. However, after 2 the output does not have the transfer matrix, and in step 3 it is computed on the fly, but the output also does not contain the transfer matrix. Consequently If I want to redo step 3 multiple times (on the same day or one week later), it presently has to recompute the transfer. Furthermore, in step 4 there is also an (unneeded) recomputation of the transfer (due to historical reasons). I believe that to be the actual annoyance that we want to address here. Being able to call ft_prepare_vol_sens multiple times would not help in the call-sequence above. It would help if the user were to mix and match low- and high-level code and insert a separate call to ft_prepare_vol_sens in between 2 and 3. But I don't like the mix-and-match approach. So low-level-only there is not a problem, but high-level-only there is, and it has to be addressed in multiple locations (among others 3 and 4, but also ft_dipolefitting and others).

Robert Oostenveld - 2013-08-07 09:23:33 +0200

(In reply to comment #20) now that I think of it, it is probably due to the way that the high-level functions call prepare_headmodel, which in turn calls ft_prepare_vol_sens. I.e. it affects roboos@mentat001> grep -l prepare_headmodel *.m ft_dipolefitting.m ft_dipolesimulation.m ft_megplanar.m ft_megrealign.m ft_prepare_headmodel.m ft_prepare_leadfield.m ft_sourceanalysis.m ft_volumerealign.m

Robert Oostenveld - 2013-08-07 09:24:31 +0200

(In reply to comment #21) oh, not ft_prepare_headmodel. And do note that ft_prepare_headmodel.m and private/prepare_headmodel,m are different!

Johannes Vorwerk - 2013-08-07 12:14:37 +0200

(In reply to comment #20) Right, that is exactly what I wanted to improve. Unfortunately at the moment I only see two options to solve this: a) add the output of the vol struct to the high level routines so that it can be stored afterwards -> a lot of code changes, but (hopefully) nice result b) add a seperate function for calculating the transfer matrix in between 2 and 3 for those users that want to store it -> little code changes in both cases a check whether a transfer matrix is already present in the vol struct inside ft_prepare_vol_sens would be necessary

Robert Oostenveld - 2021-10-29 12:17:31 +0200

This issue is from 2013. I don't think there is still work to be done, or if so, someone should pick it up again. I will close this for now.

Robert Oostenveld - 2021-10-29 12:38:34 +0200

Let me close these bugs, now that they have been resolved.