[SOLVED] How to import volume mesh from Matlab iso2mesh toolkit?

[VincentG]

With the Matlab iso2mesh toolkit, I am able to segment and process DICOM data to a volumetric mesh and save it as .stl, .dxf, .vrml or .medit.
It seems that there is only a STL loader available in SOFA. This works for visuals, but not for mechanics, apparently because in a STL file the vertex data is not shared among voxels.
The .dxf, .vrml and .medit are all ASCII textfiles with proper shared vertices, but can’t seem to find a suitable loader in SOFA.
How to proceed? Thanks for the help!

[lionel]

I think medit files have .mesh or .meshb extension.
To my knowledge there is no .dxf, vrml or .mesh/.meshb in SOFA.
I would propose to save it in meshb format and use a converter to vtk files for example.
Maybe medit or parawiew can do the job.
If not, I know that this library has a convert executable from mesh/meshb to vtk files.

[VincentG]

Thanks for the quick reply! I found a VTK export module (vtkwrite) for Matlab, which can write the iso2mesh data arrays to a VTK file. This seems to work fine.

[VincentG]

Looks like I was a bit too quick. While I can produce VTK output files, I consistently get the error:
[WARNING] [TetrahedralCorotationalFEMForceField] ERROR(TetrahedralCorotationalFEMForceField): object must have a Tetrahedral Set Topology.

Comparing my phantom_volume.vkt file contents with the provided HeartTetra_200k.vtk (which works), I see a lot of differences. Apparently, more information is needed in the output. Any specific documentation on this point? I hope there is still a way to import Matlab-meshed data into SOFA.

######### phantom_volume.vtk #########

# vtk DataFile Version 2.0
VTK from Matlab
ASCII
DATASET POLYDATA
POINTS 591 float
-81.13 92.76 -32.01 -92.76 85.00 -20.98 -81.17 92.33 -32.19
…

POLYGONS 2282 11410
4 422 544 528 79
…

######### HeartTetra_200k.vtk #########

# vtk DataFile Version 3.0
vtk output
ASCII
DATASET UNSTRUCTURED_GRID
POINTS 37923 float
71.0287 51.4427 26.4789 31.9127 65.3531 31.8623 72.5554 39.0259 56.9355
…

CELLS 204841 1024205
4 8170 8171 8172 8173
…

CELL_TYPES 204841
10
…

CELL_DATA 204841
FIELD FieldData 1
Zones 1 204841 int
6 6 6 21 6 21 21 25 7
…

POINT_DATA 37923
FIELD FieldData 1
DTI_Fibres 3 37923 double
0.204005 -0.779598 0.592123
…

[lionel]

Your phantom_volume.vtk is a surface mesh. You can see it using paraview.
You should find a vtk exporter for matlab that exports unstructured grids (e.g. volumetric meshes)

[VincentG]

I see that the word “POLYGONS” in phantom_volume.vtk makes the file reader think that it is a surface mesh, but the data itself is actually a volumetric mesh. Now I understand why Paraview consistently shows only two faces of each tetrahedra! I’m now trying to get the field names right to get it work…

[VincentG]

Importing meshes now works. I wrote my own Matlat scripts to export in VTK format with correct field names. Because there were still some errors, I also wrote a script for gmesh export, but this didn’t matter because it turned out that some tetras were too small. Now I can focus on making good quality meshes from DICOM data.

In case anyone else is interested in Matlab to VTK/gmesh export, here is the code:

function saveTetraVtk(filename,nodes,elements)
fileID = fopen(filename,'w');
fprintf(fileID,'# vtk DataFile Version 3.0\n');
fprintf(fileID,'VTK by Vincent\n');
fprintf(fileID,'ASCII\n');
fprintf(fileID,'DATASET UNSTRUCTURED_GRID\n');
fprintf(fileID,'POINTS %d float\n',size(nodes,1));
fprintf(fileID,'%3.3f %3.3f %3.3f\n',nodes');
fprintf(fileID,'\n');
fprintf(fileID,'CELLS %d %d\n',size(elements,1),5*size(elements,1));
fprintf(fileID,'4 %d %d %d %d\n',elements(:,1:4)'-1);
fprintf(fileID,'\n');
fprintf(fileID,'CELL_TYPES %d\n',size(elements,1));
fprintf(fileID,'%d\n',10*ones(size(elements,1),1));
fclose(fileID);
end

function saveTetraMesh(filename,nodes,elements)
% first column is 1..#nodes
nodes(:,2:4)=nodes(:,1:3);
nodes(:,1)=1:size(nodes,1);
elements(:,2:5)=elements(:,1:4);
elements(:,1)=1:size(elements,1);
fileID = fopen(filename,'w');
fprintf(fileID,'$NOD\n');
fprintf(fileID,'%d\n',size(nodes,1));
fprintf(fileID,'%d %3.3f %3.3f %3.3f\n',nodes');
fprintf(fileID,'$ENDNOD\n');
fprintf(fileID,'$ELM\n');
fprintf(fileID,'%d\n',size(elements,1));
fprintf(fileID,'%d 4 1 1 4 %d %d %d %d\n',elements(:,1:5)');
fprintf(fileID,'$ENDELM\n');
fclose(fileID);
end

[Hugo]

Wahoo, nice work Vincent!
Thank you so much for sharing!

[Author]

Posts