Back to the main page.
Bug 1820 - ft_headmodel_simbio needs to be updated
| Status | CLOSED FIXED |
| Reported | 2012-11-07 15:59:00 +0100 |
| Modified | 2019-08-10 12:31:23 +0200 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P3 normal |
| Assigned to: | Johannes Vorwerk |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | 15901822 |
| See also: |
Robert Oostenveld - 2012-11-07 15:59:19 +0100
The mesh generation will move to ft_prepare_mesh, the other stuff(*) is done with a mex file. the output would be something like vol.cond vol.pnt vol.tet or hex vol.index or label or so, preferably compatible with parcellations vol.stiffness After ft_prepare_vol_sens it would look like vol.cond vol.pnt vol.tet or hex vol.index or label or so, preferably compatible with parcellations vol.stiffness vol.transfer <=== this one is added *) stuff here refers to the computation of the stiffness matrix
Robert Oostenveld - 2012-11-07 15:59:38 +0100
and vol.type should be 'simbio'
Robert Oostenveld - 2012-11-07 16:16:06 +0100
the output should also have vol.unit, which is m, cm or mm
Robert Oostenveld - 2012-11-08 12:18:34 +0100
ideas for test script: 1) make (or read) a hexahedral mesh and pass it into ft_headmodel_simbio 2) make (or read) a tetrahedral mesh and pass it into ft_headmodel_simbio (not so important initially)
Robert Oostenveld - 2012-11-08 12:31:06 +0100
(In reply to comment #3) the input here could in principle be the output of the test script for bug 1815
Robert Oostenveld - 2012-11-15 08:56:01 +0100
*** Bug 986 has been marked as a duplicate of this bug. ***
Robert Oostenveld - 2012-11-15 08:58:04 +0100
I have removed the simbio source code from the repository, saving 11MB. The compilation will not be done inside fieldtrip. Only mex files have to be provided, with a pointer as to where to find teh code and compilation instructions. mac001> svn commit Deleting simbio/simbio-src Committed revision 6927.
Robert Oostenveld - 2012-11-15 09:16:36 +0100
added the code for simbio as suggested. Also use ft_hastoolbox to ensure that the path is correct. mac001> svn commit ft_prepare_vol_sens.m Sending ft_prepare_vol_sens.m Transmitting file data . Committed revision 6928.
Robert Oostenveld - 2012-11-15 09:19:16 +0100
I replaced these two by the ones I received from Johannes: mac001> svn commit Sending forward/ft_headmodel_simbio.m Sending forward/private/leadfield_simbio.m Transmitting file data .. Committed revision 6929.
Robert Oostenveld - 2012-11-15 09:40:45 +0100
I have copied the test data to /home/common/matlab/fieldtrip/data/test/bug1820, it contains bug1820/cube2mm3layervorwerk_ns_127_127_127.v bug1820/cube2mm3layervorwerk_ns_127_127_127.v.ascii bug1820/cube2mm3layervorwerk_ns_127_127_127.v.potential bug1820/elc.mat bug1820/fort.6 bug1820/matlab.mat bug1820/phantom_1.0mm_normal_crisp.mnc bug1820/phantom_1.0mm_normal_crisp.mnc.gz bug1820/phantom_1.0mm_normal_fuzzy.tar bug1820/test_forward.mat bug1820/tet_4layer_127_127_127.1.ele bug1820/tet_4layer_127_127_127.1.node bug1820/threeshell_mr.v I think that some of these test files have to be moved elsewhere, e.g. to bug 1818. -------- I adapted the test script from Johanens to the following cd /home/common/matlab/fieldtrip/data/test/bug1820 % load the mesh load test_forward % first calculate stiffness matrix, will be stored in test_stiff.stiff test_stiff = ft_headmodel_simbio(test,'conductivity',[0.33,0.0042,0.33]); % now compute transfer matrix, will be stored in test_stiff.transfer test_transfer = ft_prepare_vol_sens(test_stiff,sens); % finally compute the leadfield lf = ft_compute_leadfield(pos,sens,test_transfer); and get the error Undefined function 'calc_stiff_matrix_val' for input arguments of type 'int32'. Error in sb_calc_stiff (line 65) [diinsy,cols,sysmat] = calc_stiff_matrix_val(node,elem,cond,mele); Error in ft_headmodel_simbio (line 80) vol.stiff = sb_calc_stiff(vol); Ah, I though it was due to the int32, but it is because I am trying on a maci64 and the mex file only exists (so far) for 64 bit linux. I will move the test to a linux machine. I have added the following to give a meaningful error mac001> svn commit calc_stiff_matrix_val.m Adding calc_stiff_matrix_val.m Transmitting file data . Committed revision 6930.
Robert Oostenveld - 2012-11-15 09:44:49 +0100
I have added the test script as fieldtrip/test/test_bug1820, see https://twitter.com/fieldtriptoolbx/status/268997645752741888
Johannes Vorwerk - 2012-11-15 10:24:41 +0100
(In reply to comment #9) svn commit calc_stiff_matrix_val.mexmaci64 Hinzuf. (bin) calc_stiff_matrix_val.mexmaci64 Übertrage Daten . Revision 6932 übertragen. Now it should also work on a maci64.
Robert Oostenveld - 2012-11-15 10:46:46 +0100
(In reply to comment #11) thanks! The mac is still running, on linux I get Error using forward/private/leadfield_simbio Too many input arguments. Error in ft_compute_leadfield (line 454) lf = leadfield_simbio(pos, sens, vol); This was easy to fix roboos@mentat001> svn commit Sending forward/ft_compute_leadfield.m Transmitting file data . Committed revision 6936. Then I get >> lf = ft_compute_leadfield(pos,sens,test_transfer); Warning: an error occurred while running simbio > In forward/private/leadfield_simbio at 26 In ft_compute_leadfield at 455 Undefined function 'sb_calc_ven_loads' for input arguments of type 'double'. Error in sb_rhs_venant (line 13) loads = sb_calc_ven_loads(pos,dir,ven_nd,vol.pos); Error in leadfield_simbio (line 21) rhs = sb_rhs_venant(locpos,dir,vol); Error in ft_compute_leadfield (line 455) lf = leadfield_simbio(pos, vol);
Johannes Vorwerk - 2012-11-15 10:54:22 +0100
(In reply to comment #12) svn commit sb_calc_ven_loads.m Hinzufügen sb_calc_ven_loads.m Übertrage Daten . Revision 6937 übertragen.
Robert Oostenveld - 2012-11-16 10:14:10 +0100
The test runs, see http://fieldtrip.fcdonders.nl/development/dashboard and http://fieldtrip.fcdonders.nl/development/dashboard/r6943/test_bug1820
Robert Oostenveld - 2013-04-10 21:04:38 +0200
The test script broke recently, due to the update of the electrode projection onto the skin surface. >> test_bug1820 Warning: adding /home/mrphys/roboos/matlab/fieldtrip/external/simbio toolbox to your Matlab path Undefined function 'mesh2edge' for input arguments of type 'struct'. Error in ft_prepare_vol_sens (line 484) bnd = mesh2edge(vol); Error in test_bug1820 (line 17) test_transfer = ft_prepare_vol_sens(test_stiff,sens); 484 bnd = mesh2edge(vol); The problem is simply a missing file. roboos@mentat001> svn commit Adding forward/private/mesh2edge.m Sending plotting/private/mesh2edge.m Committed revision 7775. fixed
Lilla Magyari - 2013-09-15 13:08:48 +0200
(In reply to comment #15) hi Robert, The test_bug1820 broke again. test_bug1820 Error using sb_calc_stiff (line 80) Elements have wrong orientation or are degenerated Error in ft_headmodel_simbio (line 89) vol.stiff = sb_calc_stiff(vol); Error in test_bug1820 (line 14) test_stiff = ft_headmodel_simbio(test,'conductivity',[0.33,0.0042,0.33]); This test script is using a mesh (called test) from /home/common/matlab/fieldtrip/data/test/bug1820 in test_forward. Interestingly, when I create a mesh on my own using the script in fieldtrip/test/test_ft_prepare_mesh (see attached the mesh) this new mesh is not failing with ft_headmodel_simbio. I could not figure out what the problem is with the mesh in test_forward. Is it possible that its orientation is actually wrong and it should be updated with a new mesh? Lilla
Lilla Magyari - 2013-09-15 13:12:09 +0200
Created attachment 518 example mesh
Johannes Vorwerk - 2013-09-15 14:08:35 +0200
(In reply to comment #16) That's probably due to a change I made in sb_calc_stiff.m. I moved some checks that were previously only made in the mex-files into the matlab-part to prevent possible matlab crashes. Since it worked previously it will not be due to the orientation since that was also checked before, it is probably due to a too strong mesh adaptation which may lead to degenerated hexaeder. This is for some reason not checked in the mex-files (and thus in the SimBio-part) but should also be controled. Probably it is easiest to change the test model so that it does not contain degenerated elements any more.
Robert Oostenveld - 2013-09-16 11:18:46 +0200
(In reply to comment #18) I checked the mesh, which models a sphere. The vertices are all OK, there is no overlap and they are regularly spaced. There are 11459 hexahedra, each is described properly with a tissue type (1, 2 or 3). >> for i=1:11459 ; assert(length(unique(meshB.hex(1,:)))==8); end returns no error, so all hexahedra consist of 8 unique vertices. >> length(unique(meshB.hex(:))) ans = 13368 indicates that all vertices are being used at least once. To summarize, I see no evidence of the mesh being degenerate. Given the error "Elements have wrong orientation or are degenerated", I therefore suspect them to have an incorrect orientation. How should the orientation be defined? Say that I have a unit cube with vertex "1" at the origin and furthermore: 5----8 / /| 6----7 | | | 4 | |/ 2----3 what would be the correct sequences of vertices?
Johannes Vorwerk - 2013-09-16 11:36:16 +0200
(In reply to comment #19) Seems like I should have taken a look at the mesh before I write something... With degenerated I meant containing missshaped hexaeder, i.e. with points that lie inside (and not on) the convex hull, but this can obviously not be the case with this mesh. The orientation should be like in your figure, Robert. When I run ft_headmodel_simbio with the mesh Lilla posted on my pc, I do not get any error, which is kind of strange.
Jakob Ludewig - 2013-09-16 11:38:58 +0200
(In reply to comment #20) Right now the elements are created in the following fashion in build_mesh_hexahedral.m(b+c is just the offset for the nodes of the element): elements(:,1) = b + c; elements(:,2) = b + c + 1; elements(:,3) = b + c + (x_dim+1) + 1; elements(:,4) = b + c + (x_dim+1); elements(:,5) = b + c + (x_dim + 1)*(y_dim+1); elements(:,6) = b + c + (x_dim + 1)*(y_dim+1) + 1; elements(:,7) = b + c + (x_dim + 1)*(y_dim+1) + (x_dim+1) + 1; elements(:,8) = b + c + (x_dim + 1)*(y_dim+1) + (x_dim+1); which would correspond to something like 8----7 / /| 5----6 | | | 3 | |/ 1----2 if I'm not mistaken. I just numbered the vertices counter-clockwise for the top and bottom face(which are squares). I think we had this kind of problem before, fixed it and arrived at this solution. It is rather strange that it should reappear. Nevertheless, this should be easy to fix. To arrive at the numbering that you proposed we would just have to change the above lines to: elements(:,1) = b + c + (x_dim+1); elements(:,2) = b + c; elements(:,3) = b + c + 1; elements(:,4) = b + c + (x_dim+1) + 1; elements(:,5) = b + c + (x_dim + 1)*(y_dim+1) + (x_dim+1); elements(:,6) = b + c + (x_dim + 1)*(y_dim+1); elements(:,7) = b + c + (x_dim + 1)*(y_dim+1) + 1; elements(:,8) = b + c + (x_dim + 1)*(y_dim+1) + (x_dim+1) + 1; Actually I'm not sure what kind of orientation/numbering the simbio routines require but once we found out it is just a matter of permuting the vertice numbering(I hope...). Greetings Jakob
Johannes Vorwerk - 2013-09-16 11:42:21 +0200
(In reply to comment #21) This is surely the same orientation as in Roberts example, since only the local orientation of the nodes is important and the 90° rotation of the whole cube shouldn't matter.
Robert Oostenveld - 2013-09-16 13:10:43 +0200
when I wrote 5----8 / /| 6----7 | | | 4 | |/ 2----3 I only thought of labeling the vertices, not of mapping them on the 1x8 vector of indices. But if I understand Johannes, then 8 vertices like this should be described with a hexaheder like mesh.hex(1,:) = [1 2 3 4 5 6 7 8] That means that (with isotropic voxels) the vertex distances can be checked. E.g. dist(1-2)=1, dist(2-7)=sqrt(2) and dist(2-8)=sqrt(3). So there are a number of checks (including the ones I did in my previous comment) which could be implemented in a stand-alone function.
Johannes Vorwerk - 2013-09-16 13:30:05 +0200
(In reply to comment #23) That's right Robert, I thought you already meant is as the numbering in the .hex vector. Acutally, these checks are implemented in sb_check_ori and this should be the function that caused the output which Lilla posted. The question is, why the above posted mesh gave an error for Lilla while it works fine here.
Lilla Magyari - 2013-09-18 13:17:43 +0200
(In reply to comment #24) hi all, I am really sorry for causing all the misunderstanding. probably, my message above was not clear. The mesh I attached (called meshB) worked totally fine. I created it last week. The mesh which is called "test" and is located in /home/common/matlab/fieldtrip/data/test/bug1820 in test_forward.mat gave an error. And I tried but could not figure out what was the difference between them. I did not attach the wrong mesh ("test") because I thought that can be reached from the network (although I realize now that Johannes and Jakob probably can not reach it.) Moreover, the quesiton was if we should just change the mesh (called "test") in test_forward to a newly created one which probably won't give an error. This is just important because the test script fails at the moment because of the wrong mesh. But I will just change it then! And thanks for all the reactions.
Lilla Magyari - 2013-09-18 13:19:08 +0200
Created attachment 519 this contains the mesh called test which gives an error
Lilla Magyari - 2013-09-24 16:53:08 +0200
(In reply to comment #26) I've fixed the test script. Lilla
Eelke Spaak - 2015-02-11 13:58:47 +0100
The test script started failing again: Error using sb_calc_stiff (line 64) Dimensions of vol.cond and entries of vol.tissuelabel do not fit! Error in ft_headmodel_simbio (line 89) vol.stiff = sb_calc_stiff(vol); Error in test_bug1820 (line 86) test_stiff = ft_headmodel_simbio(mesh,'conductivity',[0.33,0.0042,0.33]); 64 error('Dimensions of vol.cond and entries of vol.tissuelabel do not fit!');
Johannes Vorwerk - 2015-02-20 15:33:33 +0100
(In reply to Eelke Spaak from comment #28) I wrote a fix for this, unfortunately I have problems with my svn access at the moment. Could someone replace line 28 in .private/prepare_mesh_hexahedral.m for i=1:numel(fn),if numel(mri.(fn{i}))==prod(mri.dim), segfield=fn{i};end;end by for i=1:numel(fn),if (numel(mri.(fn{i}))==prod(mri.dim))&(~strcmp(fn{i},'inside')), segfield=fn{i};end;end and check this in? Thanks!
Eelke Spaak - 2015-03-02 11:15:48 +0100
(In reply to Johannes Vorwerk from comment #29) Thanks, just committed. I can confirm this fixes the error I mentioned above. Note that I have not run the entire test script again yet, as this takes quite a while. We will see whether it is flagged as passing again on the next full test run.