VTK support

[jtgreen]

Another big ask; would you mind loading an example where a file, such as that at

https://github.com/UK-Digital-Heart-Project/Statistical-Shape-Model/blob/master/Both_ED_mean.vtk

is loaded in such a way that is compatible with the previous biventricular instantiation?

[finsberg]

You could try to convert the mesh using meshio?

[jtgreen]

I imported the vtk but it doesn’t expose the relevant methods you use in the biv example; I’m a coder and a physician but in no way a veteran cardiac modeler, so the complex geometries are still new to me. If it’s a trivial task for you it might be of general interest to see how other types like vtk (which for some reason seems to be the preferred publication method in open source repositories for many cardiac meshes) are imported. I know I’d appreciate it!

[jtgreen]

That is I imported the vtk with meshio, but it doesn’t expose the same methods on the geometry object that the raw mesh does.

[jtgreen]

I’ll try a few other simple conversions with meshio and post back

[finsberg]

I was not able to read the meshes in the link you provided. Seems like meshio does not support reading vtk polydata. You can check out some of the meshes that I haved listed here: https://computationalphysiology.github.io/cardiac_geometries/third_party.html

(for example https://kcl.figshare.com/articles/dataset/A_Virtual_Cohort_of_Twenty-four_Left-ventricular_Models_of_Ischemic_Cardiomyopathy_Patients/16473903)

[jtgreen]

OK, I will. Not that it's this project's problem, but what I'm trying to do is populate a biventricular model (healthy subject, no atria) with your algorithm. In the third party link above there is a zenodo link I have also been working with, namely https://zenodo.org/record/4419784#.ZAUUK-zMI1J
which contains the vtu/vtp. I'm having a hard time importing the files to dolfin. I've come up with the following:

###########################

import meshio
import vtk

# Convert vtp to vtu for meshio
reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName("heart_sur.vtp")
reader.Update()
vtp = reader.GetOutput()
append_filter = vtk.vtkAppendFilter()
append_filter.AddInputData(vtp)
append_filter.Update()
vtu_surface = append_filter.GetOutput()
# Set filename and input
writer = vtk.vtkXMLUnstructuredGridWriter() # is it a vtu or vtp?
writer.SetFileName("heart_sur_vtp.vtu")
writer.SetInputData(vtu_surface)
writer.Update()

# Write
writer.Write()

# Read vtp and vtu as vtu's
msh_vtu = meshio.read("heart_vol.vtu")
msh_vtp = meshio.read("heart_sur_vtp.vtu")

# Convert volume mesh, equivalent of the commandline:
# > meshio-convert inputfile.vtu outputfile.xdmf
meshio.write("heart_vol.xdmf", msh_vtu)

##################################

but mapping the points and element types (RV/LV etc) from the VTP 'class' element back to the VTU is proving difficult. The number of element tags matches the number of points in the VTP, but not the number of cells. You could use the points to disambiguate, but there seems to be truncations in the VTP not present in the VTU, so you'd have to match rounded values. This may be an ignorance thing, but I can't find a good way to map the element tags to the triangular boundary faces back to the overall points.

Again, I know this isn't an issue with your project, but would love to use your application here to populate that model.

Thanks in advance,
John

[finsberg]

I think you are doing the right thing. For this particular mesh it looks like the markers are defined in the nodes on the boundary mesh, so you could perhaps do something like

import dolfin

meshio.write("heart_vol.xdmf", msh_vtu)

mesh = dolfin.Mesh()

with dolfin.XDMFFile("heart_vol.xdmf") as infile: 
    infile.read(mesh)

bmesh = dolfin.BoundaryMesh(mesh, "exterior")

Now msh_vtp.point_data["class"] contains the markers for the points in bmesh (but I think this is a bit strange... ).

[jtgreen]

Do you think your code will work with markers defined at the nodes?

I think ideally what I'd like to commit here are some examples using vtk/vtu/vtp and not only that I'd like to show an example of taking the output and populating it with a purkinje network like from here: https://github.com/fsahli/fractal-tree

The contribution of course being upstream input variation and downstream augmentation, like with adding a purkinje system.

I'm working on both, but I think they'd be good examples to add to the demos section.

[jtgreen]

As an aside, do you know how to get paraview to show the glyphs demonstrating the angles? I see the angles represented as the color variation, but the vector of the glyph arrows is all down, instead of representing the fiber angle.

[finsberg]

Do you think your code will work with markers defined at the nodes?

No this would not work in the current version. I think what you need to do is to loop through all the point and then set a value for the facets in the facet function that is appropriate. Of course this can be a bit problematic when a node is shared across multiple surfaces. I can try to come up with a strategy for this.

I think ideally what I'd like to commit here are some examples using vtk/vtu/vtp and not only that I'd like to show an example of taking the output and populating it with a purkinje network like from here: https://github.com/fsahli/fractal-tree

The contribution of course being upstream input variation and downstream augmentation, like with adding a purkinje system.

I'm working on both, but I think they'd be good examples to add to the demos section.

As an aside, do you know how to get paraview to show the glyphs demonstrating the angles? I see the angles represented as the color variation, but the vector of the glyph arrows is all down, instead of representing the fiber angle.

Not entirely sure what you mean, but if you are using the option 3D glyphs then this might resulting in all arrows pointing down. In stead you should select the Glyph filter.

[finsberg]

I too am struggling to use Paraview to create a view like the one shown in the attached picture. I did try the Glyph filter, but not sure how to compute the fibre_angle in Paraview. Any corresponding script for Paraview to generate this result is highly appreciated.

I too was initially using 3D Glyph till I read this post. Now, I do see angles at +60 and -60, but dont quite see the picture as attached..

Thanks!

PS: attaching my output below.. Cant seem to make it prettier than this. Also looks like arrows are not marked on all cells...

Take a look at this video

Screen.Recording.2024-08-07.at.22.22.11.mp4