Back to the main page.
Bug 1815 - ft_prepare_mesh should make hexahedral mesh
| Status | CLOSED FIXED |
| Reported | 2012-11-07 15:46:00 +0100 |
| Modified | 2013-09-24 14:44:04 +0200 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P4 enhancement |
| Assigned to: | Lilla Magyari |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | 1822 |
| See also: |
Robert Oostenveld - 2012-11-07 15:46:07 +0100
this is needed, since ft_prepare_headmodel with method=simbio needs the hex-mesh.
Robert Oostenveld - 2012-11-08 12:05:32 +0100
I suggest to make a test script that looks like this seg = zeros(101, 91, 81); % slightly different numbers in x, y, z transform = eye(4); % update this so that the center of the volume is in [0 0 0] create a sphere in this seg with ones. put them together into something compatible with ft_datatype_segmentation call ft_prepare_mesh to make a triangulated surface mesh call ft_prepare_mesh to make a hexahedral mesh -> this requires code changes call ft_prepare_mesh to make a tetrahedral mesh -> this requires code changes use ft_plot_mesh for plotting all of them (ft_plot_mesh cannot do the hexahedral mesh yet, so it also needs to be extended) seg = zeros(101, 91, 81); % slightly different numbers in x, y, z create three sphere in the seg with value 1, 2 and 3. repeat the mesh generation and plotting.
Jakob Ludewig - 2012-11-20 16:17:54 +0100
(In reply to comment #0) I talked to Johannes and I will focus on this bug for now. Cheers
Lilla Magyari - 2012-11-21 16:48:59 +0100
(In reply to comment #2) hi Jakub and Johannes, I am trying to use the prepare_volume_mesh_segmentation function that Johannes sent me a few days ago, but I couldn't run it. I did the following: I copied this function and also vgrid.m to fieldtrip/external/simbio, I created an example segmentation structure (fieldtrip style): disp(singlesphere) dim: [91 104 111] transform: [4x4 double] coordsys: 'ctf' unit: 'mm' seg: [91 104 111] % this is binary mask of a sphere and I called: cfg=[]; cfg.tissue = 'seg'; mesh = prepare_volume_mesh_segmentation(cfg,singlesphere) but I get an error: vgrid is writing the wireframe mesh file, this may take some time ... ./tp26314fb9_afbb_43e8_9e1e_19dbd6020202.sh: line 2: /home/mrphys/lilmag/matlab/fieldtrip/external/simbio/vgrid: No such file or directory elapsed time: 0.015383 Trying to read FE mesh from file tpd1aede87_4d9e_4059_a32f_43a868ad17a1.v. Could not open file tpd1aede87_4d9e_4059_a32f_43a868ad17a1.v for reading. ??? Error using ==> read_vista_mesh Error while reading mesh! Error in ==> prepare_volume_mesh_segmentation at 121 [mesh.pnt,mesh.hex,labels] = read_vista_mesh(meshfile); Could you help me why it does not work out for me? What did you use an input for the function? Or is this a problem with the path? Lilla
Johannes Vorwerk - 2012-11-21 19:56:24 +0100
(In reply to comment #3) Hi Lilla, you have to put the vgrid.m exactly into the folder where the vgrid binary is located, since it is just a dummy-file so that matlab can locate the binary. Usually this would be /external/simbio/vgrid1.3.1/program/. You have to obey that the binary is compiled under linux at the moment, I might have a try to compile it on a mac the next days. Best, Johannes
Lilla Magyari - 2012-11-22 13:13:48 +0100
(In reply to comment #4) hi Johannes, Thanks. I am always using the latest svn version, so the binary file is located for me where it is located in the latest fieldtrip version. I put the vgrid.m to /external/simbio/vgrid1.3.1/program/ (where I could find also vgrid.C, is this the binary file? sorry, I do not really know which one that should be.) and I got the following error message: >> mesh = prepare_volume_mesh_segmentation([],singlesphere) vgrid is writing the wireframe mesh file, this may take some time ... /home/mrphys/lilmag/matlab/fieldtrip/external/simbio/vgrid1.3.1/program/vgrid: Option -min has incorrect value -max. Usage: /home/mrphys/lilmag/matlab/fieldtrip/external/simbio/vgrid1.3.1/program/vgrid
Johannes Vorwerk - 2012-11-22 14:00:31 +0100
(In reply to comment #5) Ok, that was an error I already fixed but forgot to mail you, see e-mail. Sorry, my bad.
Lilla Magyari - 2012-11-22 14:53:52 +0100
(In reply to comment #6) hey, no problem. that's why I am testing it. and thanks! it worked! Lilla
Lilla Magyari - 2012-11-30 16:46:48 +0100
(In reply to comment #7) hi Johannes and Jakob, I am working on the prepare_volume_mesh_segmentation function that Johannes sent me a week ago. I try to fit this among the other fieldtrip mesh functions. We (me and Robert) planned to have two separate functions: 1. prepare_mesh_hexahedral for preparing hexahedral meshes from a segmentation 2. and prepare_mesh_tetrahedral for tetrahedral meshes. But I am realizing now during running prepare_volume_mesh that it depends on the resolution what kind of mesh (tetrahedral or hexahedral) I get. How is the resolution specification implemented? When resolution =1, is the output always a hexahedral mesh independent of the size of the volume? and when resolution = 2 (or larger) is the mesh always tetrahedral? Or does the output depend also on the volume size? Is it possible to specify in vgrid to output a hexahedral mesh in any case independently of the resolution? In general: what does resolution mean? thanks for the help! Lilla
Johannes Vorwerk - 2012-12-03 17:40:41 +0100
(In reply to comment #8) Resolution in this case means how vgrid proceeds with the given vista-image, resolution 1 means a direct transformation, where each voxel is transformed to a hexaeder. Resolution 2 will then mean 8 voxels to a hexaeder, non-integer steps would force a resampling. But, since I am not totally sure, how this resampling is executed, I think we should mainly not use it and better do a resampling inside matlab with fieldtrip-functions and then create a mesh with vgrid afterwards, which would probably then need to be transformed in matlab. Forcing vgrid to always create cubic elements should already be done by the switch "-elem cube" (line 118), did you experience cases, where this did not happen? I never experienced that before, I would have to have a look at that. Regarding splitting creation of tetrahedral and hexahdreal meshes: This function is meant for creating hexahedral meshes only. vgrid is not a program one should use to create tetrahedral meshes, since it simply splits up hexahedras into tetrahedras to achieve this, which mainly means a much higher computational effort, without having the advantages of tetrahedral meshes (better shape of compartment boundaries etc.). For tetrahedral meshes indeed a second function needs to be written, which should then implement a tetgen call.
Lilla Magyari - 2012-12-04 18:15:15 +0100
(In reply to comment #9) (In reply to comment #9) hi Johannes, thanks a lot for your answer. "Resolution 2 will then mean 8 voxels to a hexaeder, non-integer steps would force a resampling." when I run prepare_volume_mesh_segmentation with no resolution specified, or with with cfg.resolution = 1 the output structure looks like this: disp(mesh) pnt: [283008x3 double] hex: [267731x8 double] tissue: [267731x1 double] tissuelabel: {'101' '102' '103'} tissuename: {'tissue_1' 'tissue_2' 'tissue_3'} but when cfg.resolution = 2, the output is like this: mesh4=prepare_volume_mesh_segmentation(cfg,seg); disp(mesh4) pnt: [37224x3 double] hex: [166855x4 double] tissue: [166855x1 double] tissuelabel: {'1'} tissuename: {'tissue_3'} As you can see, in the hex field, there are actually 4 points defining a polyhedron. Therefore, I was asking if it makes tetrahedral meshes. And actually, I also get a warning on the screen: % vgrid is writing the wireframe mesh file, this may take some time ... % WARNING: max_maxdim must divide nx,ny,nz! % WARNING: enlarging image from size 91, 104, 111 to 92,106,112 % WARNING: gridGenerator: Changing grid generation to 'tetra5' % elapsed time: 3.231 % Trying to read FE mesh from file tp4f96190f_6895_4d28_b487_f938a19f7b37.v. % Successfully read mesh with 37224 nodes and 166855 elements. which says that it is changing the generation to 'tetra...'. I am wondering that if you do not suggest to use the resolution option from the vgrid, then maybe we could just have resolution = 1, without offering the option to change it. My other question is: is it important what are the tissuelabels, right? Are these numbers only for differentiating the different tissue-types (in case there are more) or do they represent which hexahedron belongs to skin, brain... etc.? Is it ok for example if the input contains only the segmentation of the headsurface (here: tissue_3) and then it gets tissuelabel 101 (and not 103... for example)? cfg=[]; cfg.tissue = {'tissue_3'}; mesh3=prepare_volume_mesh_segmentation(cfg,seg); disp(mesh3) % pnt: [37224x3 double] % hex: [33371x8 double] % tissue: [33371x1 double] % tissuelabel: {'101'} % tissuename: {'tissue_3'} I tested this function with test_volume_mesh_segmentation which you can find in the fieldtrip/test directory (it contains the example segmentation structure). Best, Lilla
Robert Oostenveld - 2012-12-04 21:17:41 +0100
not to go into any details, but just as a hint: you can do tmp.pnt = bnd.pnt; tmp.hex = bnd.hex(1,:); % get the first one and idem for tet and tri. Subsequently you can do ft_plot_mesh(tmp) to get a quick view of the first (or any other) element. You can also plot all elements with ft_plot_mesh(bnd) where facecolor=none and edgecolor=b would be suitable for a tetrahedral or hexahedral wireframe. But it is very slow in MATLAB if there are many elements. It works much better if you do ft_write_headshape(filename, bnd, 'format', 'ply') and use the free Paraview or Meshlab software for the display (see google). That allows you to look at the whole mesh. But looking at an individual element (see the first lines) is also useful.
Johannes Vorwerk - 2012-12-05 12:04:06 +0100
(In reply to comment #10) Ok, I found the error, will fix it soon. For the second part: The tissuelabels do not directly fix which values in the tissues vector correspond to which compartment, it is more important to identify the different compartments clearly so that we are able to assign the correct conductivities later on.
Johannes Vorwerk - 2012-12-05 12:18:56 +0100
(In reply to comment #12) If you replace lines 107+108 in prepare_volume_mesh_segmentation.m by sb_write_materials(materialsfile,[1:numel(cfg.tissue)],tissue,cfg.resolution); it should work now.
Robert Oostenveld - 2012-12-13 13:04:47 +0100
mac001> svn commit prepare_mesh_hexahedral.m Sending prepare_mesh_hexahedral.m Transmitting file data . Committed revision 7176. changed to reflect the new handling of teh vista executable, do not change into tempdir: that causes problems when it fails (and does not cd back).
Lilla Magyari - 2012-12-18 17:49:14 +0100
Created attachment 395 hexagonal mesh with 1 hexagon
Lilla Magyari - 2012-12-18 17:50:01 +0100
Created attachment 396 hexahedral mesh with 27 hexahedrons
Lilla Magyari - 2012-12-18 17:54:42 +0100
(In reply to comment #13) hi Johannes, I attached the screenshots of two hexahedral meshes (sorry for the file descriptions, I always mix up hexagons with hexahedrals). I created an example segmentation (script for it in test/test_prepare_mesh_hexahedral) with really small radius (1 and 2 mms) and so I got meshes which contained only a few elements. Then I wrote them with ft_write_headshape into ply format and I opened them in paraview. I am wondering if the hexahedrons should look like this... Could you take a look on the images, please? thanks. Best, Lilla
Lilla Magyari - 2012-12-21 14:38:49 +0100
(In reply to comment #17) hi all, I plotted the meshes also with ft_plot_mesh, and they looked fine. So, it seems that problem arises when ft_write_headshape writes the mesh into 'ply' format. I will make an other bug for this, and I think we can conclude now that the hexahedral meshes are fine. Therefore, I will close this bug. P.s. I also inserted a ft_hastoolbox('simbio',1) in prepare_mesh_hexahedral, otherwise it doesn't run. Best, Lilla
Lilla Magyari - 2013-09-13 16:24:15 +0200
I created a new test script: test_ft_prepare_mesh. The content of the test script of this bug (test/test_bug1815) is moved to there. I will delete test_bug1815.