Back to the main page.
Bug 1818 - ft_read_headshape has to be made consistent with new segmentation/parcellation
| Status | CLOSED FIXED |
| Reported | 2012-11-07 15:53:00 +0100 |
| Modified | 2013-06-30 15:16:24 +0200 |
| Product: | FieldTrip |
| Component: | forward |
| Version: | unspecified |
| Hardware: | PC |
| Operating System: | Mac OS |
| Importance: | P3 normal |
| Assigned to: | Jakob Ludewig |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | 18221840 |
| See also: |
Robert Oostenveld - 2012-11-07 15:53:48 +0100
at this moment there is case 'vista' ft_hastoolbox('simbio', 1); [nodes,elements,labels] = read_vista_mesh(filename); shape.pnt = nodes; if size(elements,2)==8 shape.hex = elements; elseif size(elements,2)==4 shape.tet = elements; else error('unknown elements format') end shape.index = labels; that should be changed such that it is compatible with FT_DATATYPE_PARCELLATION.
Robert Oostenveld - 2012-11-07 18:33:24 +0100
It might not be possible to make it fully compatible. As altrnative it can be considered to do the reformatting of the mesh structure into a parcellation in ft_read_atlas.
Robert Oostenveld - 2012-11-08 12:09:46 +0100
to do: make a test script. This test script can only be made if we have a variety of test data. So to do first: get test data.
Johannes Vorwerk - 2012-11-09 10:41:34 +0100
After a closer look at the parcellation datatype it seems to me as a reasonable option to store the meshes according to it. The extension of the datatype to add the option of assigning anatomical labels to elements (triangles, tetrahedrons, hexahedrons etc.) instead of assigning them to points makes very much sense in my opinion (it even has the advantage that you gain a unique coverage of your domain). The question, which I cannot overlook, is, whether there would occur problems in other routines relying on this datatype when adding this possibility.
Robert Oostenveld - 2012-11-09 11:18:59 +0100
(In reply to comment #3) Other functions calling ft_read_headshape expect the output to contain a pnt field, not a pos field. That is a small and silly difference, but still. So a mesh has a pnt, and a parcellation has a pos. Although the ft_datatype_xxx definitions are meant to standardize "objects" and ensure them to be consistent in different parts of the code, there are still "objects" that are not properly defined. That is largely due to there not being enough usage of the objects and consensus on how they should look. So here we have this problem, i.e. ft_datatype_mesh or ft_datatype_headshape are simply not defined (as function), since the mesh or headshape object are still ambiguous. But what we can easily do is to extend ft_datatype_parcellation such that pnt is renamed in pos. Other FT functions that expect a parcellation will always be passing their input data through ft_datatype_parcellation anyway. This on-the-fly conversion (in this case trivial) is something that we also do in other cases, e.g converting between source and volume, raw and timelock, etc. I suggest Lilla stays in charge of this bug. Could you send her some different example files that she can work on? The test case would then look like mesh = ft_read_headshape('yourfile.xxx'); parcellation1 = ft_datatype_parcellation(mesh) % should rename ont into pos; parcellation2 = ft_datatype(mesh, 'datatype', 'parcellation;) % should also work and then check for the different files and that the MATLAB structures are as expected.
Robert Oostenveld - 2012-11-12 18:08:47 +0100
Lilla and I were confused about file formats and software packages that write them. We spent some time with google... ------- Vista employs a data file format that allows any collection of images, edge vectors, and other structured objects to be stored in a single file. This allows, for example, both an edge set and the image from which it was derived to be kept in one file. http://www.cs.ubc.ca/nest/lci/vista/vista.html http://www.cs.ubc.ca/nest/lci/vista/file.html -------- Vgrid is a fast and robust octree-based 3D mesh generator for unstructured gridswhich works directly on segmented (labeled) 3D images http://www.rheinahrcampus.de/~medsim/vgrid/ http://www.rheinahrcampus.de/~medsim/vgrid/manual.html Vgrid saves the data in one of these formats: - hf selects Vista graph format (default). - vm and ipe are variations of the vm format - ascii selects a simple home-grown ASCII mesh format. Volume meshes only. - gmv selects the GMV ASCII format (http://www-xdiv.lanl.gov/XCM/gmv/GMVHome.html). For both surface and volume meshes. - gpp selects the GPPView/GeoFEM format (http://geofem.tokyo.rist.or.jp/download_en/gppview.html). Volume meshes only. - complex2d selects the complex2d format of GrAL (http://www.math.tu-cottbus.de/~berti/gral). For surface meshes only. - dx selects the format of OpenDX (http://www.opendx.org). For surface meshes only. -------- TetGen is a program to generate tetrahedral meshes of any 3D polyhedral domains. TetGen generates exact constrained Delaunay tetrahedralizations, boundary conforming Delaunay meshes, and Voronoi partitions. http://tetgen.berlios.de/ TetGen supports a bunch of file formats: .node files .poly files .smesh files .ele files .face files Other supported file formats (recognized by file extensions): .off files .mesh files .ply files .stl files http://tetgen.berlios.de/fformats.html In general it seems that there is always a *.node file for the vertices with a corresponding file that describes how the vertices are combined. ------------ It would be good to have more example files from tetgen and from vista. Vgrid does not seem to have its own file format, but uses formats defined elsewhere. Also for those example files might be useful.
Johannes Vorwerk - 2012-11-13 10:45:09 +0100
I'm not sure where this should lead, even though there is a variety of possible output formats for both programs, I would think that it should be sufficient to be able to read their native/standard output format? This would be the vista format for vgrid and *.ele with a *.node for tetgen. Of course, it is always advantageous to be able to handle as many file formats as possible, but I am not sure whether there is any need for the handling of further file formats at this moment. If necessary, read functions for further file formats should be quite easy to implement.
Robert Oostenveld - 2012-11-13 11:20:23 +0100
(In reply to comment #6) > I'm not sure where this should lead, ... It should lead to a complete set of test files and a clear description of what file can act as input and/or output to which software in case an excursion out of the ideal fieldtrip pipeline is needed. For example the "vista format" seems to be a container that can contain multiple geometrical representations. How does a hypothetical user get an anatomical MRI that was segmented with AFNI and written in the Analyze format into the SIMBIO pipeline? This might involve segmentation = ft_read_mri(analyze_file) ft_write_mri(vista_file, segmentation) call vgrid on the command line, have it write the result to the hf(?) vista format mesh = ft_read_headshape(hf_vista_file) headmodel = ft_headmodel_simbio(mesh) etc. So here I am seeing an interaction between analyze (could also have been nifti, or minc, or dicom) and two vista formats. Reading from "the" vista format would not be sufficient to allow this excursion. That is why I want to have an inventory of formats that ft_read_headshape could support. Some formats seem to be more appropriate for ft_read_mri. For some formats we need the respective write functions. These functions does not have to support all formats, but we have to be able to communicate to users how to implement a particular pipeline. This is not only about programming, but as much about documentation. Remember that someone has to be using the code, and that that someone does not have your expertise. Comment 5 could at a certain point be changed into a FAQ: "how do I make an excursion to software XX for the SIMBIO pipeline?". The answer to that needs to address both the "why" and the "how".
Johannes Vorwerk - 2012-11-13 11:51:07 +0100
(In reply to comment #7) "It should lead to a complete set of test files and a clear description of what file can act as input and/or output to which software in case an excursion out of the ideal fieldtrip pipeline is needed." The question is, whether it should be a description of what can act as input and/or output or a description what should act as input and/or output. There is a bunch of different formats, but most of theme are needless in our case, so that a documentation of these would only confuse the user. Thus, from my point of view a complete set of test files for the parts of ft_read_headshape needed for simbio consists of a *.v-mesh and a *.ele with *.node. The other file formats should not be used, so that it should be clearly stated in the documentation to use the standard file formats when taking an excursion to get a simbio-headmodel. Do we agree here? Additionally, there should be test files for the input of vgrid and tetgen, which would be a MRI/segmentation in the vista format for vgrid and a set of surfaces in the *.poly format for tetgen, but this should not be handled by ft_read_headshape of course.
Robert Oostenveld - 2012-11-13 12:23:10 +0100
(In reply to comment #8) whether a format is needed depends on the users needs, not on the developer's E.g. George probably has his rat brain files in *.stl format, not in a vista format. So he should be enabled to convert them to vista format. The SIMBIO developers can indeed work with the limited perspective that they don't have to care about the users and only worry about their own specific formats, the FieldTrip developers should care about bridging the gaps. To repeat myself: I don't want to make read/write code for all formats. But in order to determine usefulness and which formats need support, we have to know which formats exist. So far there was no list, now there is. I suggest we move on with getting things done. @Lilla, how far are you with making test scripts and assigning bugs to the others? @Johannes if you have more files that you consider to be useful for testing any read/write code, please pass them on to Lilla.
Johannes Vorwerk - 2012-11-13 15:16:10 +0100
(In reply to comment #9) Ok, then maybe I didn't get your point correctly before. Including the read/write functions for the vista MRI format to ft_read/write_mri might make sense so that we would need a test file for that too, I'll search for one...
Lilla Magyari - 2012-11-13 17:34:00 +0100
hi Jakob, I assigned this bug for you. You can find a test script in fieldtrip/test/test_bug1818.m This script is testing if the ft_read_headshape function can read in vista and tetgen meshes in an appropriate FieldTrip format. When I run the test script, I get an error because ft_read_headshape can not read in the tetgen mesh (but it works well with vista mesh). When test script can run without an error, the meshes are also in the appropriate format, and this bug is solved. Please, take a look on bug1819 as well, because also labels for the elements should be added for the tetgen mesh. You will probably need to change the path to the example files in the test script. Thanks and let us know if you have questions. Lilla
Johannes Vorwerk - 2012-11-14 10:17:57 +0100
(In reply to comment #11) Hey Lilla, there where some smaller mistakes in the test script, if you correct them it should work fine (at least it did on my pc...). 1. ft_read_headshape doesn't recognize the tetgen-format automatically so that you have to add ft_read_headshape(...,'format','tetgen') to the call. 2. When reading in tetgen meshes, you shouldn't give the full filename but only the name without ending, so without the .1.ele/.1.node part. Furthermore, you don't read .ele and .node seperately but at once. 3. To create a parcellation afterwards, you already need the information which element belongs to which compartment. To achieve this, you have to include what I posted for bug 1819. A corrected script for the tetgen part looks like this: % tetgen mesh = ft_read_headshape('/home/common/matlab/fieldtrip/data/test/bug1818/tet_4layer_127_127_127','format','tetgen'); parcellation1 = ft_datatype_parcellation(mesh); parcellation2 = ft_datatype_parcellation(mesh,'parcellationstyle','probabilistic'); if~(ft_datatype(parcellation1,'parcellation')) error('the conversion to a parcellation failed'); end if~(ft_datatype(parcellation2,'parcellation')) error('the conversion to a parcellation failed'); end Best, Johannes
Robert Oostenveld - 2012-11-14 11:09:28 +0100
(In reply to comment #12) > 1. ft_read_headshape doesn't recognize the tetgen-format automatically so that > you have to add ft_read_headshape(...,'format','tetgen') to the call. Then that is something to be fixed in fileio/ft_filetype if possible. Please have a look there.
Robert Oostenveld - 2012-11-14 11:14:20 +0100
(In reply to comment #12) > 2. When reading in tetgen meshes, you shouldn't give the full filename but only > the name without ending, so without the .1.ele/.1.node part. Furthermore, you > don't read .ele and .node seperately but at once. This has to be made robust to user specification of either one of the files. I suggest something along the line of [p, f, x] = fileparts(inputfilename) elefile = fullfile(p, [f '.ele']); nodefile = fullfile(p, [f '.node']); for consistency with other fieldtrip behavior, the following should be made to work: [f, p] = uigetfile('*', 'please select one of the tetgen files'); and then allow the user to do mesh = ft_read_mesh(fullfile(p, f))
Johannes Vorwerk - 2012-11-14 12:06:33 +0100
Ok, if I think about it it's perhaps a bit more complicated. The problem is, that a .ele file necessarily has to come with a .node file, but not the other way round. A .node might also be the only file, just containing a node list, or come with a .face file containing surface definitions. Would it be reasonable to let the user enter the name of the .ele file if he wants to read the mesh?
Robert Oostenveld - 2012-11-14 12:16:12 +0100
(In reply to comment #15) This is not unique to the problem we are solving here, there are many cases like this. E.g. (known for Lilla) the brainvision and ctf dataset formats have a variety of files that belong together. ft_filetype can do something like if node and ele exist, then it is tetgen_ele elseif node and surf exist, then it is tetgen_surf elseif only node exist, then it is tetgen_node In the tetgen_ele section in ft_read_headshape you know that you need a ele and a node file and hence can do [p, f, x] = fileparts(inputfilename) elefile = fullfile(p, [f '.ele']); nodefile = fullfile(p, [f '.node']); In the tetgen_surf section in ft_read_headshape (*) you know that you need a surf and a node file and hence can do [p, f, x] = fileparts(inputfilename) surffile = fullfile(p, [f '.surf']); nodefile = fullfile(p, [f '.node']); and in the tetgen_node section just read the nodes, there is no filename confusion. Note that we do not need the tetgen_surf section for now, I am just mentioning it for illustrative purposes.
Johannes Vorwerk - 2012-11-14 14:08:58 +0100
(In reply to comment #16) Ok, I have now added elseif (filetype_check_extension(filename, '.ele') || filetype_check_extension(filename, '.node')) && exist([p f '.node']) && exist([p f '.ele']) type = 'tetgen_tet'; manufacturer = 'TetGen - A Quality Tetrahedral Mesh Generator and a 3D Delaunay Triangulator'; content = 'tetgen tetrahedral mesh'; elseif (filetype_check_extension(filename, '.face') || filetype_check_extension(filename, '.node')) && exist([p f '.node']) && exist([p f '.face']) type = 'tetgen_fac'; manufacturer = 'TetGen - A Quality Tetrahedral Mesh Generator and a 3D Delaunay Triangulator'; content = 'tetgen surface mesh'; elseif filetype_check_extension(filename, '.node') type = 'tetgen_nod'; manufacturer = 'TetGen - A Quality Tetrahedral Mesh Generator and a 3D Delaunay Triangulator'; content = 'tetgen nodes '; end to ft_filetype and case 'tetgen_tet' % reads in the tetgen format and rearranges according to FT conventions % tetgen files also return a 'faces' field (not used here) [p, f, x] = fileparts(filename); shape = rmfield(shape,'fid'); IMPORT = importdata([p f '.node'],' ',1); shape.pnt = IMPORT.data(:,2:4); IMPORT = importdata([p f '.ele'],' ',1); shape.tet = IMPORT.data(:,2:5); if size(IMPORT.data,2)==6 shape.tissue = IMPORT.data(:,6); end % add label (see above at vista) % shape.tissue % shape.tissuelabel case 'tetgen_fac' error('TetGen .face not yet implemented!') case 'tetgen_nod' error('TetGen .node not yet implemented!') to ft_read_headshape, which works fine on my machine.
Johannes Vorwerk - 2012-11-14 14:16:43 +0100
One more point for Lilla: If you assign the labels by % representation of data is compatible with ft_datatype_parcellation shape.tissue = zeros(size(labels)); numlabels = size(unique(labels),1); shape.tissuelabel = {}; for i = 1:numlabels ulabel = unique(labels); shape.tissue(labels == ulabel(i)) = i; shape.tissuelabel{i} = num2str(ulabel(i)); end it might happen, that you change the labels of the different tissues when reading a mesh. Is there a need for the parcellation datatype, that the tissues start at 1 and are increasing in steps of one? Otherwise it would probably be better to keep the original labels to avoid confusion, which would mean shape.tissue(labels == ulabel(i)) = ulabel(i); or easier shape.tissue = labels;
Robert Oostenveld - 2012-11-14 16:03:21 +0100
(In reply to comment #18) > Is there a need for the parcellation datatype, that the tissues > start at 1 and are increasing in steps of one? mesh.tissue contains the indices into mesh.tissuelabel, so the following should work: for i=1:length(mesh.tissue) fprintf('the tissue is %s\n' mesh.tissuelabel{mesh.tissue(i)}); end It would be possible to represent it like mesh.tissue = [10 11 12] mesh.tissuelabel = {[], [], [], [], [], [], [], [], [], 'label10', 'label11', 'label12'} but that would not be as nice. The index numbers should not be interpreted to start with, they serve to look it up. The consequence of the empty tissies later in the code is that we would also have to represent conductivities like conductivity = [ inf, inf, inf, inf, inf, inf, inf, inf, inf, cond10, cond11, cond12 ] Note that due to the requirement of mesh.tissue to be used to index other fields that negative numbers are not possible. Fractional numbers (floats) are also not possible, nor is zero.
Jakob Ludewig - 2012-11-20 16:16:14 +0100
(In reply to comment #11) Ok, I'll look into it. Greetings Jakob
Lilla Magyari - 2012-11-21 17:17:25 +0100
(In reply to comment #19) hi Robert, can I close this bug? We created this bug primarily for reading in fem tetrahedral and hexahedral meshes and converting them to FT-style matlab structures. The test script (test_bug1818) for this bug runs without error. I can file a new bug for testing other meshes as well. Lilla