Back to the main page.
Bug 1756 - ft_voltype not work with bem, or 'bnd' field not detected correctly
| Status | CLOSED FIXED |
| Reported | 2012-10-01 12:46:00 +0200 |
| Modified | 2014-02-24 10:56:34 +0100 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Windows |
| Importance: | P3 normal |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | 937 |
| See also: |
Johanna - 2012-10-01 12:46:35 +0200
since r5287, bug was introduced with variable 'type' being probed but without existing, on line 105. elseif isfield(vol, 'bnd') && strcmp(type, {'dipoli', 'asa', 'bemcp', 'openmeeg'}) type = any(strcmp(type, {'dipoli', 'asa', 'bemcp', 'openmeeg'})); how should in fact a vol of type BEM(subtype one of 4 bems) be established? Not all vol with .bnd are BEM.
Robert Oostenveld - 2012-10-02 11:52:42 +0200
I made a test script manzana> svn commit test_bug1756.m Adding test_bug1756.m Transmitting file data . Committed revision 6629. But I don't understand the problem. If you want to know that it is a subtype, you ask for the specific subtype. Can you explain once more what is bugging you?
Johanna - 2012-10-02 13:49:24 +0200
There is no problem if vol.type exists. And it does in the bemcp and dipoli created in the test_bug1756.m But try >> bemcp1=bemcp; >> bemcp1=rmfield(bemcp1,'type') and then >> ft_voltype(bemcp1) will now crash on the 'elseif' line mentioned in initial bug Description, since that line is currently illogical. Do all 'bem' types have both a .bnd and .mat? Can that be used as a check if a BEM? (I initially stumbled on this bug with a vol that had .bnd field but was not a BEM).
Johanna - 2012-10-02 13:51:35 +0200
Related question: is it ok/intended that in the output of the test function: >> size(bemcp.bnd(1).pnt) ans = 900 3 >> size(bemcp.bnd(2).pnt) ans = 600 3 >> size(bemcp.bnd(3).pnt) ans = 300 3 and >> size(dipoli.bnd(1).pnt) ans = 300 3 >> size(dipoli.bnd(2).pnt) ans = 600 3 >> size(dipoli.bnd(3).pnt) ans = 900 3 In other words, the order of the number of points in the bnd is swapped between dipoli and bemcp.
Robert Oostenveld - 2012-10-02 14:52:33 +0200
I just committed the change that I also mentioned to you in person: if ~isempty(unit) % use the user-specified units for the output - vol.unit = geometry.unit; + vol.unit = unit; elseif isfield(geometry, 'unit') % copy the geometrical units into he volume conductor - vol.unit = geometry.unit; + % assume that in case of multiple meshes that they have the same units + vol.unit = geometry(1).unit; end Sending forward/ft_headmodel_concentricspheres.m Transmitting file data . Committed revision 6631.
Robert Oostenveld - 2012-10-02 14:53:30 +0200
(In reply to comment #3) that is allowed if you ask me: dipoli has its own preferred order, whereas ft_prepare_mesh resurns them according to the users request.
Robert Oostenveld - 2012-10-02 18:08:42 +0200
manzana> svn commit test/test_bug1756.m forward/ft_voltype.m Sending forward/ft_voltype.m Sending test/test_bug1756.m Transmitting file data .. Committed revision 6635.
Johanna - 2012-10-03 09:18:55 +0200
Hi Robert, I'm reopening as I don't think your change [&& isfield(vol,'type')] fixed the problem. If isfield(vol,'type'), then that is already taken care of on line 85, so this line on 105 will never be reached. Also, 'type' still doesn't exist on line 105; perhaps you meant to change it to && strcmp(vol.type,....)? But still my question is: what should be done with a 1) vol that has no .type and no .mat but does have .bnd? 2) vol that has no .type but does have both .bnd and .mat? (sorry I didn't make this clear in the first description of the bug).
Robert Oostenveld - 2012-10-03 10:05:40 +0200
(In reply to comment #7) I indeed misread/misunderstood the code. So line 108 should deal with - it does not contain a type - based on the structure content, what does it look like? We could do elseif isfield(vol, 'bnd') && isfield(vol, 'mat') type = 'bem'; % it could be dipoli, asa, bemcp or openmeeg and then change around line 120 switch desired case 'bem' type = any(strcmp(type, {'dipoli', 'asa', 'bemcp', 'openmeeg'})); into type = any(strcmp(type, {'bem', 'dipoli', 'asa', 'bemcp', 'openmeeg'})); However, not that >> t_bemcp t_bemcp = bnd: [1x3 struct] cond: [0.3300 0.0041 0.3300] skin_surface: 3 source: 1 mat: [300x1800 double] unit: 'mm' cfg: [1x1 struct] >> t_dipoli t_dipoli = bnd: [1x3 struct] cond: [0.3300 0.0041 0.3300] mat: [1800x1800 double] unit: 'mm' cfg: [1x1 struct] the structures are different. The respective vol.mat should only be used with the corresponding code in th_compute_leadfield->eeg_leadfieldb. So the bemcp "mat" should not be used with the "openmeeg" code. That means that the bem type should not get lost. If it gets lost, there is a problem that cannot be fixed. The only that can be said after the vol.type got lost is that it is a bem model, but forward computations are impossibile.
Robert Oostenveld - 2012-10-03 10:07:05 +0200
the relevant section in the test code then becomes assert(ft_voltype(t_bemcp, 'bem')); assert(ft_voltype(t_dipoli, 'bem')); assert(~ft_voltype(t_singleshell, 'bem'));
Robert Oostenveld - 2012-10-03 10:09:19 +0200
(In reply to comment #9) manzana> svn commit test/test_bug1756.m forward/ft_voltype.m forward/ft_compute_leadfield.m Sending forward/ft_compute_leadfield.m Sending forward/ft_voltype.m Sending test/test_bug1756.m Transmitting file data ... Committed revision 6642. -------- Are you satisfied with this solution? If not, drop me a call.
Johanna - 2012-10-03 16:42:29 +0200
I will commit the change we discussed regarding if isfield(vol,'forwpar'). Now however, ft_voltype still cannot solve the artifical situation suggested at the beginning of test_bug937 (which itself had bugs in, but after fixing and simplification, still reduces to): [pnt, tri] = icosahedron162; ssvol.bnd(1).pnt = 10*pnt; ssvol.bnd(1).tri = tri; ssvol.bnd(2).pnt = 50*pnt; ssvol.bnd(2).tri = tri; ssvol.bnd(3).pnt = 60*pnt; ssvol.bnd(3).tri = tri; then ft_voltype(ssvol) will give 'unknown'. But of course tcfg=[]; tcfg.headshape=ssvol.bnd; ssvolcs=ft_prepare_concentricspheres(tcfg) results in ft_voltype(ssvolsc) with correct output 'concentricspheres'. Question: 1) Should a 'multishell' description exist from ft_voltype to name the 'ssvol'? This example is a weird hybrid between BEM (3 shells contained within each other) and concentricspheres (except defined by .pnt rather than .r/.o).
Robert Oostenveld - 2012-10-03 16:49:57 +0200
(In reply to comment #11) no, such a multishell should not exist. The ssvol is a structure that contains a struct-array with three meshes. The bnd field is a geometry, but not yet a volume conduction model. For it to become a volume conduction model, more information/data needs to be added, such as tissue conductivities and potentially other numerical parameters.
Johanna - 2012-10-04 09:53:45 +0200
Now committing to ft_voltype and updated the test script here to reflect. svn commit 6671. elseif isfield(vol, 'bnd') && isfield(vol, 'forwpar') type = 'singleshell'; elseif isfield(vol, 'bnd') && numel(vol.bnd)==1 type = 'singleshell'; I close this bug.
Johanna - 2012-10-05 09:04:55 +0200
reopened test_bug1756 fails now.
Robert Oostenveld - 2012-10-05 10:04:52 +0200
(In reply to comment #14) ah, this is because I removed the default conductivity values yesterday. What do you think about having default values? For 3-shell bem and 3-sphere they are doable perhaps, because the skin and brain are equal and hence reordering the meshes (inside/outside ot the other way around) does not affect the outcome.
Johanna - 2012-10-05 18:15:19 +0200
(In reply to comment #15) Yes, I think default conductivity values for 3shell (BEM or concentric spheres, and where the order can safely be reversed) is not a problem/confound, and is in fact welcome for those researchers who forget exactly what those default values from the literature should be. (e.g. I remember the ratio of 1/80 but not the actual values). I guess then for 4shell that cfg.conductivity should be required so as to ensure the order is correct. But perhaps a useful error message could be given, suggesting values to be used if the order is 'scalp, skull, csf, brain' for example.
Robert Oostenveld - 2012-10-07 11:34:33 +0200
(In reply to comment #16) We should check whether the order has the expected consequences in the code. I see a potential problem with the user giving bnd(1) to bnd(4) as input and cond=1x4, but not knowing about the explicit order of the surfaces. Explicit feedback would be useful, e.g. created sphere 1 (radius = xxx) with conductivity c1 (innermost, i.e. brain) created sphere 2 (radius = xxx) with conductivity c2 created sphere 3 (radius = xxx) with conductivity c3 created sphere 4 (radius = xxx) with conductivity c4 (outermost, i.e. scalp) this calls for a test script with all possible permutations of the surfaces ....
Robert Oostenveld - 2012-10-07 11:37:40 +0200
(In reply to comment #17) regarding default values there is now mbp> grep 'conductivity.*\[' ft_prepare_* ft_prepare_bemmodel.m:% cfg.conductivity = [Cskin_surface Couter_skull_surface Cinner_skull_surface] ft_prepare_concentricspheres.m:% cfg.conductivity = conductivity values for the model (default = [0.3300 1 0.0042 0.3300]) ft_prepare_concentricspheres.m:% cfg.conductivity = [0.3300 1 0.0042 0.3300] ft_prepare_concentricspheres.m:if ~isfield(cfg, 'conductivity'), cfg.conductivity = []; end % this should be specified by the user ft_prepare_headmodel.m: cfg.conductivity = ft_getopt(cfg, 'conductivity', []); ft_prepare_headmodel.m: cfg.conductivity = ft_getopt(cfg, 'conductivity', []); and mbp> grep 'conductivity.*\[' forward/ft_headmodel_* forward/ft_headmodel_bemcp.m: conductivity = [1 1/80 1] * 0.33; forward/ft_headmodel_bemcp.m: conductivity = [1 1/80 1 1] * 0.33; forward/ft_headmodel_dipoli.m: conductivity = [1 1/80 1] * 0.33; forward/ft_headmodel_dipoli.m: conductivity = [1 1/80 1 1] * 0.33; forward/ft_headmodel_fns.m:% conductivity = matrix C [9XN tissue types] forward/ft_headmodel_openmeeg.m: conductivity = [1 1/80 1] * 0.33; forward/ft_headmodel_openmeeg.m: conductivity = [1 1/80 1 1] * 0.33; which is clearly not consistent with the desired behavior.
Johanna - 2012-10-19 14:04:12 +0200
Hi Robert, I'm working now on the updates to ft_prepare_headmodel (bug 1745), and in the case of concentric spheres, if the user gives 3 boundaries but no conductivity values, then as we discussed previously, it should be ok to assume default values of conductivity = [1 1/80 1] * 0.33; Where should this default be placed? Within the top of ft_headmodel_concentricspheres.m, or in the higher function ft_prepare_headmodel that calls it?
Robert Oostenveld - 2012-10-20 08:39:16 +0200
(In reply to comment #19) I suggest to put it in ft_headmodel_concentricspheres, and in the respective BEM versions. In those functions the rationale can be explained (with one sentence of comments). And in the high-level function it would have be be conditional on a variety of characteristics, whereas in the low-level function it is much easier to check that numel(conductivity) is consistent with the number of compartments. E.g. for concentric spheres set teh default at line 29 and then around line 90 if numel(conductivity)~=length(vol.r) give an error end
Johanna - 2012-10-26 11:15:53 +0200
Hi Robert, I'm working on this change now (defaults for 3 boundaries) but have noticed another thing: an inconsistency in naming of the conductivity subfield name in the vol: concentricspheres uses vol.c but the BEM use vol.cond. (Old versions of concentricspheres used vol.cond but this was changed at some point). The low-level code (ft_compute_leadfield or eeg_leadfieldb.m) know the correct name for a given method, thus it runs correctly, but 1) should we make it consistent over all methods? (I vote yes.) 2) If yes, which should it be? .c or .cond? (I actually like .cond better as it's more descriptive, but .c is more consistent with the short .r and .o subfields).
Robert Oostenveld - 2013-12-05 09:43:57 +0100
(In reply to comment #21) I believe that all voltype detection issues are resolved. That leaves this thread with a one request, which is to use vol.cond instead of vol.cI agree that it is better for consistency. I have changed vol.c into vol.cond. I have also fixed a minor bug in ft_datatype_headmodel (it had strcmp, whereas it should be isfield) and added '2013' as the latest version. If all goes well, the user specified cfg.vol.c is now automatically upgraded to cfg.vol.cond. I tried to find all code that uses it (including test code) and updated it. mac001> svn commit Sending forward/ft_compute_leadfield.m Sending forward/ft_headmodel_concentricspheres.m Sending forward/ft_headmodel_singlesphere.m Sending forward/private/eeg_leadfield1.m Sending forward/private/eeg_leadfield4.m Sending forward/private/eeg_leadfield4_prepare.m Sending ft_prepare_concentricspheres.m Sending test/test_bug1042.m Sending test/test_bug1988.m Sending test/test_bug2148.m Sending utilities/ft_datatype_headmodel.m Transmitting file data ........... Committed revision 8963.