Back to the main page.
Bug 1954 - ft_prepare_head_model has singular matrices with 'bemcp'
| Status | CLOSED FIXED |
| Reported | 2013-01-24 18:46:00 +0100 |
| Modified | 2015-07-15 13:31:05 +0200 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | All |
| Operating System: | All |
| Importance: | P3 normal |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: | http://bugzilla.fcdonders.nl/show_bug.cgi?id=2817 |
Tilmann Sander-Thommes - 2013-01-24 18:46:19 +0100
Running ft_defaults load('~/segmentedmri.mat'); cfg = []; cfg.method = 'bemcp'; vol = ft_prepare_headmodel(cfg, segmentedmri); yields ... weight = -0.155225 , defl = 0.000000 Nverta = 1600 Ntrib = 4796 Warning: Matrix is singular to working precision. > In ft_headmodel_bemcp at 132 In ft_prepare_headmodel at 206 In test_vol_mod_seg at 78 ... and the resulting vol contains NaNs. It works with "cfg.method = 'dipoli';" I am not an expert on BEM, but there might really be something wrong or at least a clarification needed. Many thanks, Till
Johanna - 2013-01-25 12:05:41 +0100
Hi Till, Thanks for spotting this. I'm able to confirm your finding, and also that these NaN values occur in the 'old' way of doing bemcp headmodels as well. The dipoli method does not have a similar .mat field in which to assess if it went ok or not. In all cases, ft_prepare_mesh is called and this results in a bnd(1) which is not optimal (spikes coming out of the head). Here is my test function for anyone else to try: %% load standard_mri cfg = []; cfg.output = {'brain','skull','scalp'}; segmentedmri = ft_volumesegment(cfg, mri); cfg = []; cfg.method = 'bemcp'; vol1 = ft_prepare_headmodel(cfg, segmentedmri); % vol1.mat are all NaN cfg = []; cfg.method = 'dipoli'; vol2 = ft_prepare_headmodel(cfg, segmentedmri); % looks ok, but no .mat to assess % old way cfg=[]; cfg.numvertices=3000; bnd=ft_prepare_mesh(cfg,segmentedmri); % ft_plot_mesh(bnd(1)) shows that the segmentation did not go well. segmentedmri.bnd=bnd; cfg=[]; cfg.method='bemcp'; vol1o=ft_prepare_bemmodel(cfg,segmentedmri); % vol1o.mat also all NaN
Johanna - 2013-01-25 12:57:41 +0100
I trace it to line 131 of ft_headmodel_bemcp where C21 becomes a matrix of all NaN. I feel like I've seen this before, but can't remember when/where/why. My guess is an imperfect segmentation leading to an imperfect mesh. I cc Robert.
Tilmann Sander-Thommes - 2013-01-25 13:08:15 +0100
In external/bemcp is a test script called bemcp_example.m . This works fine until after ft_prepare_bemmodel, but the meshes are just spheres geometrically well separated ("onion"). The problem is that the code from the website http://fieldtrip.fcdonders.nl/tutorial/headmodel_eeg really works (it uses dipoli). But you said, one does not have a success/fail flag for dipoli ? So if the mashes have sikes coming out of the head dipoli should fail or maybe it has some safeguards. This is a bit worrying altogether. Actually a chinese collaborator of mine found this problem. And we both are NOT experts on bem. Many thanks, Till
Robert Oostenveld - 2013-01-25 13:40:12 +0100
(In reply to comment #1) Hi Johanna, Could you post the specific segmentedmri an the meshes on home/common with the test script? I wonder whether there is something wrong with the segmentation, like some voxels outside the head. Lilla (CC) and me are working on improving the handling of the segmentations, but I hope we did not temporarily introduce a bug. If the segmentation has some rogue voxels, one of the inner meshes may stick out the outer meshes. The BEM (any one of them) requires the meshes to be nested, which is then not the case any more. At this moment there is only a very limited sanity check on the meshes, as those checks take quite some time. It is on the todo list to implement these sanity checks, at least to allow the user to do them him/herself. Robert
Johanna - 2013-01-25 14:05:51 +0100
I have commited test/test_bug1954. I did not upload the segmented mri or meshes, to save space, as they are computed from loading /home/common/matlab/fieldtrip/template/headmodel/standard_mri.mat But indeed figure;imagesc(squeeze(segmentedmri.scalp(:,110,:))) shows voxels outside (on top of) the head. And ft_plot_mesh(bnd(1)) shows voxels on top of head, as well as that the bottom of the head got cut off (so indeed, probably the inner brain bnd extends below this outer scalp bnd.).
Lilla Magyari - 2013-01-25 15:09:50 +0100
(In reply to comment #5) interestingly, in line 513 of ft_volumesegment scalpmask still shows the little bar at the head. However, volumethreshold(scalpmask) in line 509 should remove the "little parts" outside of the head. Is it possible that something is not working right at volumethreshold and volumefillholes just before that? Jan-Mathijs, could you, please, look at it? Lilla
Lilla Magyari - 2013-01-25 15:12:19 +0100
and now it goes to Jan-Mathijs, too. And I have just assigned this bug to myself, but please, feel free to assigned it to yourself if you know what the problem is.
Jan-Mathijs Schoffelen - 2013-01-25 16:38:45 +0100
Created attachment 415 aliasing in standard mri
Jan-Mathijs Schoffelen - 2013-01-25 16:46:07 +0100
the aliased stuff is attached to the scalp mask, due to the relatively low (default) threshold. Increasing the scalpthreshold to 0.3 seems to get rid of the problem. Should we adjust the default?
Jan-Mathijs Schoffelen - 2013-01-25 16:46:38 +0100
...but there still seem to be nans in the bemcp bem
Jan-Mathijs Schoffelen - 2013-01-25 16:50:07 +0100
the nans are now probably caused by holes in the scalp surface. this situation calls for my new volumeedit function ;-)
Robert Oostenveld - 2013-01-28 12:02:09 +0100
(In reply to comment #11) The bar at the top is indeed a nasty one, but a good example of what can go wrong. I recall having edited away in 2001 manually with some matlab code. Too bad I did not have the JM editor back then. It also calls for some checks when the mesh is made from the segmentation. E.g. with bwlabeln -> how many clusters are there? (should be one per tissue type).
Jan-Mathijs Schoffelen - 2013-01-28 13:06:53 +0100
bash-4.1$ svn commit -m "enhancement - manually edited anatomy to remove the aliasing above the head that caused scalp segmentation issues" standard_mri.mat Sending standard_mri.mat Transmitting file data . Committed revision 7414. I used private/volumeedit.m to erase the aliased part. Scalp morphology now looks better, but still the bemcp complains about matrix singularity. All three compartments consist of 1 blob after sgmentation.
Robert Oostenveld - 2013-01-28 13:09:27 +0100
(In reply to comment #13) might the order of the compartments be incorrect?
Johanna - 2013-01-28 13:29:46 +0100
I think the order is ok, as this is the printout I get from calling: cfg = []; cfg.method = 'bemcp'; vol1 = ft_prepare_headmodel(cfg, segmentedmri); triangulating the outer boundary of compartment 1 (scalp) with 3000 vertices triangulating the outer boundary of compartment 2 (brain) with 3000 vertices triangulating the outer boundary of compartment 3 (skull) with 3000 vertices the call to "ft_prepare_mesh" took 9 seconds and required the additional allocation of an estimated 0 MB reordering the boundaries to: 2 3 1 I also checked that each was inside each other, just by quick visual inspection and they look ok (though can't promise I checked every .pnt!)
Robert Oostenveld - 2013-01-28 13:48:41 +0100
(In reply to comment #13) mac001> svn commit -m "restructuring - added a README file to each of the template directories with a pointer to the wiki (where appropriate). I also moved the headmodel readme.txt content over to the wiki, reformatted it and extended it." Adding template/atlas/README Adding template/electrode/README Adding template/headmodel/README Deleting template/headmodel/README.txt Adding template/layout/README Adding template/sourcemodel/README Transmitting file data ..... Committed revision 7416. @JM can you please document the change you made at http://fieldtrip.fcdonders.nl/template/headmodel#standard_mri I suggest something like === Revision information === Up to r7413 (Jan 2013) the file content corresponded to the original anatomical MRI. From r7414 onwards the anatomical MRI was updated with ... (describe changes).
Tilmann Sander-Thommes - 2013-01-28 14:43:06 +0100
Just a reminder: The test script bemcp_example.m in external/bemcp, shows that bemcp basically works with the chosen spherical geometry. Therefore if the segmentation is OK for a real geometry, but bemcp fails, then something very subtle is going on.
Lilla Magyari - 2013-07-23 22:03:57 +0200
(In reply to comment #17) hi all, first of all I re-run bemcp on the standard_mri in the template directory and I got the matrix singularity problem and the NaNs in the matrix. The earlier artifact of the segmentation (caused by the bar in the mri) is gone, but the surfaces are really close to each other (see attached jpg.) I was wondering if this caused the problem. Therefore, I run bemcp also on the tutorial segmentation (in home/common/matlab/fieldtrip/data/ftp/tutorial/headmodel_eeg/segmentedmri) where I am almost sure that the segmentation is fine. However, I got the same error! I got the warning message and the matrix is full of NaNs. Is it possible that somehow bemcp is just not working properly? I assign this bug back to joint development user. Lilla
Lilla Magyari - 2013-07-23 22:07:55 +0200
Created attachment 497 brain skull segmentation of standard mri (see lowest point)
Johanna - 2015-01-29 16:36:53 +0100
Perhaps Robert's recent change: https://code.google.com/p/fieldtrip/source/detail?spec=svn10143&r=9863 may have partly or fully fixed this bug. The test_bug1954 now produces: vol1.mat all as numbers, but vol1o.mat are all still NaN. I'm not sure I understand the difference of why it works for vol1 and not vol1o.
Robert Oostenveld - 2015-02-12 10:30:46 +0100
I made some updates to the test script and code. I don't know what was initially causing the problem, but it is all working fine now. I indeed suspect that the fix that Johanna pointed to is the solution. mac011> svn commit Sending forward/ft_headmodel_bemcp.m Sending private/prepare_mesh_segmentation.m Sending private/volumethreshold.m Sending test/test_bug1954.m Transmitting file data .... Committed revision 10225.
Tilmann Sander-Thommes - 2015-02-12 11:28:20 +0100
The bug was originally communicated to me from a chinese colleague. I will inform him about the solution. Many thanks, Till
Johanna - 2015-02-16 18:19:34 +0100
Hi Robert, Sorry about this, but I just check test_bug1954 and it crashed for me, as well as reported on the dashboard. I think some change in private/volumethreshold.m wasn't right. http://fieldtrip.fcdonders.nl/development/dashboard/r10235/test_bug1954 This is in the call to ft_volumesegment (thus perhaps unrelated to the bemcp of this bug).
Robert Oostenveld - 2015-02-16 22:28:30 +0100
(In reply to Johanna from comment #23) You are right, there was a stupid error in volumethreshold. I fixed it. Furthermore, the dipole part was failing because of a correctly detected problem with the meshes. At the bottom of the head the brain and skull were going too far down, causing the meshes to touch the skin surface. I updated the segmentation code, a bit like the strategy at http://fieldtrip.fcdonders.nl/tutorial/natmeg/dipolefitting#construct_the_eeg_volume_conduction_model roboos@mentat001> svn commit Sending forward/ft_headmodel_dipoli.m Sending private/volumethreshold.m Sending test/test_bug1954.m Transmitting file data ... Committed revision 10238.
Johanna - 2015-02-17 10:57:05 +0100
seems to work now, and great addition to check the overlapping boundaries. thanks!
Robert Oostenveld - 2015-02-22 14:02:39 +0100
*** Bug 2849 has been marked as a duplicate of this bug. ***