VTK and VMTK Deciphered: A Step-by-Step Tutorial and Hands-On MICCAI-inspired Example
This is a short tutorial about some of the VTK and VMTK tools, mainly the ones I used for the preprocessing of the data in the paper I’m presenting at MICCAI 2023. While VTK and VMTK are widely known for their visualization capabilities, this tutorial ventures beyond visual representation. Even though this tutorial delves into a specific application, it will give you an understanding of the inner workings of VTK and VMTK. My goal is to equip you with the necessary skills to navigate their documentation and apply these powerful tools to your projects successfully. All code snippets are written in Python, as both frameworks have Python-interpreted wrappers and can be installed easily with Pip
VTK
VTK, short for the Visualization Toolkit, is a powerful open-source software designed to manipulate and visualize scientific data. It has state-of-the-art tools for 3D rendering, with tools for 3D interaction and 2D graphics capacity. The platform is used both in commercial applications and for research and development. What makes it really powerful is the pipeline it is built around, its open-source nature, and its active community. New tools are developed daily, and every user can create their own. Although it can be hard to navigate at first, the various data types it has are key features, as they allow the processing and visualization of many kinds of data, whether structured or unstructured.
VTK architecture is built around a “pipeline” concept, this pipeline has two main blocks: the first one is for importing and processing data, and the second block is for visualizing and exporting data.
Let´s see in detail each step of the pipeline:
- The first element in the pipeline consists of the sources, serving as input components. In my case, I utilize readers as input, specifically vtkObjReader and vtkPolyDataReader.
- Subsequently, the filters come into play. VTK encompasses a wide range of processing algorithms, all of which are implemented as filters. These filters can be chained together, using the output from a filter as the input to another one. Each filter has an output port that stores the information and data produced by the filter. The input to the filter can be recovered from the input port. An example of two filters chained would be:
filter2.SetInputConnection(filter1.GetOutputPort())- Mappers transform data into graphic primitives. They are an abstract way of specifying what we are going to visualize.
- Actors represent an object in the scene (geometry plus display properties). With them, you can specify stuff such as color, opacity, shading, and orientation.
- Renderers and windows describe the rendering onto the screen in a platform-independent way. [1]
The update mechanism is a vital component that gives coherence to the entire pipeline. When the settings of filters or sources undergo changes, all subsequent filters, mappers, and actors that depend on them are automatically updated. Conversely, if an object downstream in the pipeline requires previous information from an object upstream, it can access it through input and output ports. This streamlined communication between pipeline components allows for efficient data flow and ensures that the visualization remains synchronized with any modifications made to the input.
VTK Interface
File interface
Readers are source objects, they digest data from a file and create a data object that can go through the pipeline. Writers are mappers, similarly, they ingest data and write it into a file
Importers and exporters read and write files that have more than one object. They are generally used to save or read a whole scene.
VTK Data representation
Data object
A data object can be perceived as an assembly of data lacking any specific structure. Only when organized into a particular arrangement they take on a form, allowing us to apply visualization algorithms for further processing.
Dataset
Data objects with an organizing structure and associated data attributes form datasets. Most algorithms in VTK operate on datasets.
The structure has two parts: topology and geometry. Topology is the set of properties invariant under certain geometric transformations [2]. Geometry represents the concrete realization of topology, providing the specific positioning of elements within three-dimensional space.
The dataset model assumes that the structure consists of cells and points. The cells are the ones that specify the topology, while the points specify the geometry.
Cell Types
A dataset comprises one or more cells, analogous to biology, cells serve as the fundamental elements. A cell is defined by a type and an ordered list of points, referred to as the connectivity list. The cell's topology is implicitly defined by combining the type and the connectivity list. The x-y-z coordinates of the points specify the geometry of the cell.
Dataset types
A dataset’s nature can be classified based on its structure, either as regular or irregular. A dataset is considered regular when it exhibits a single consistent mathematical relationship among its composing points and cells. When the points follow a regular pattern, the dataset’s geometry is also classified as regular. If the topological relationship between the cells is regular, then the topology of the dataset is regular. The main advantage of regular (or structured) is that it can be implicitly represented, saving memory and computation resources
Irregular (or unstructured) data requires explicit representation because there is no inherent pattern that can be succinctly described. VTK offers five implemented datasets to handle various data types: vtkPolyData, vtkImageData, vtkStructuredGrid, vtkRectilinearGrid, and vtkUnstructuredGrid. Each dataset is designed to cater to different data structures and characteristics, providing flexibility for representing and visualizing a wide range of data types. You can find detailed information about each of them in [3]. However, for this example, I will provide a summary of polygonal data, which I’ll be using.
The polygonal dataset (vtkPolyData) consists of vertices, polyvertices, lines, polylines, polygons, and triangle strips. Polygonal data exhibits an unstructured nature in both its topology and geometry. The dataset’s cells can vary in their topological dimensions, making it a versatile and flexible representation of complex and diverse datasets.
VMTK
The Vascular Modeling Toolkit (VMTK) is a collection of libraries and tools for 3D reconstruction, geometric analysis, mesh generation, and surface data analysis for image-based modeling of blood vessels. The toolkit provides essential functionalities for tasks such as vessel segmentation, centerline extraction, surface reconstruction, and quantitative analysis of blood vessels. Its user-friendly interface and extensible architecture make it an accessible and valuable resource for medical professionals, researchers, and developers working in the domain of vascular imaging and modeling.
The toolkit empowers users to efficiently process and analyze various vascular imaging modalities, including CT angiography (CTA) and magnetic resonance angiography (MRA). Moreover, its support for 3D visualization and interactive exploration of vascular models allows for a comprehensive understanding of vascular morphology and facilitates effective communication of findings in a clinical or research setting.
VMTK is based on C++ classes (VTK-based algorithms) and has high-level functionality through Python classes (each class is a script). Additionally, there exists PypeS — Python pipeable scripts, a framework that facilitates interaction between vmtk scripts. This framework enables seamless communication and data flow between different vmtk scripts, enhancing the capabilities and versatility of the toolkit.
VMTK It inherits VTK file formats, adding a few ones
Let´s get practical
You can download the Intra 3D dataset at their GitHub repository (https://github.com/intra3d2019/IntrA). I will be using the meshes on the folder generated/vessel.
Let´s go through the step-by-step process of extracting the centerline from a blood vessel mesh:
The first thing we need is a reader. According to the file type to read, we´ll choose the reader. For example with an obj file:
filename = "test.obj"
reader = vtk.vtkOBJReader()
reader.SetFileName(filename)
reader.Update()
data = reader.GetOutput()The usage is quite straightforward, we create the reader object, provide it with the name of the file to read, use update() to run the pipeline, and obtain the polydata with GetOutput().
We can also create a function that chooses the correct reader according to the file´s format:
def read_vtk(filename):
if filename.endswith('xml') or filename.endswith('vtp'):
polydata_reader = vtkXMLPolyDataReader()
else:
polydata_reader = vtkPolyDataReader()
polydata_reader.SetFileName(filename)
polydata_reader.Update()
polydata = polydata_reader.GetOutput()
return polydataNow that you understand how it works, let´s see how to read the Intra dataset (https://github.com/intra3d2019/IntrA). The mesh files are .obj so we can read them with an obj reader and export them with a polydata writer as a .vtp file (remember the dataset type “polygonal data”)
reader = vtk.vtkOBJReader()
reader.SetFileName(filename)
reader.Update()
data = reader.GetOutput()
writer = vtk.vtkXMLPolyDataWriter()
writer.SetFileName("ArteryObjAN1–14.vtp");
writer.SetInputData(data);
writer.Write()This code works as a .obj to .vtp file converter, which is exactly what we need as input for the next step.
Centerline extraction
To extract the centerline of the mesh we can use the vmtk script network extraction, with the vmtk pypes framework. It looks like this:
from vmtk import pypes
myPype=pypes.PypeRun("vmtknetworkextraction -ifile ArteryObjAN129–10.vtp -advancementratio 1 -ofile centerlines/ArteryObjAN129–10.vtp")This way we obtain the centerline for each mesh in the dataset. As you can see, pypes is a great tool that allows achieving complex tasks with very few lines of code. All you have to do is specify the script, in this case, “vmtknetworkextraction”, input and output file with “-ifile” and “-ofile”, and finally, the advancement ratio which is “the ratio between the sphere step and the local maximum radius” according to the documentation. The output file is a vtk polydata containing the centerline coordinates. You can find detailed documentation of network extraction at [4]
Cross-section extraction
Although there is a cross-section calculation script in vmtk, it didn’t work with my data, likely because of the closed endings of most meshes. As a result, we will explore an alternative, more manual approach to accomplish the task.
We already have the centerline coordinates from the previous step, we now need the mesh’s radius at each of these centerline points. The idea is to cut the mesh with planes at each centerline point, and obtain the cross-section area at each plane cut. Subsequently, we calculate the radius using the cross-section area obtained, approximating it as a circle.
Let’s go step by step:
- Iterate through the points of each branch: We’re gonna use the tangent of the point with its immediate prior and posterior as the cut plane normal. Not separating branches for this calculation may result in problems at the bifurcations
The file containing the centerline is organized in a way that each cell (remember vtk cells) contains one centerline branch. With centerline.GetNumberOfCells() we get the number of branches and with centerline.GetCell(j) the j-th branch.
So, to iterate through the mesh we´ll need two nested for loops, one for the branches and one for the points within the branch.
for j in range(centerline.GetNumberOfCells()):
numberOfCellPoints = centerline.GetCell(j).GetNumberOfPoints()
for i in range (numberOfCellPoints):
#Now I´m on the i-th point of the j-th branch. 2. Calculate the point´s tangent:
The first thing to do at each point is calculate the tangent, as the average vector between this point and its prior and posterior (only prior/posterior at the beginning or end of the branch)
tangent = np.zeros((3))
weightSum = 0.0;
#point i at branch j
if (i>0):
point0 = centerline.GetPoint(points_Acum-1);
point1 = centerline.GetPoint(points_Acum);
distance= np.sqrt (vtk.vtkMath.Distance2BetweenPoints(point0,point1));
tangent[0] += (point1[0] - point0[0]) / distance;
tangent[1] += (point1[1] - point0[1]) / distance;
tangent[2] += (point1[2] - point0[2]) / distance;
weightSum += 1.0;
if (i<numberOfCellPoints-1):
point0 = centerline.GetPoint(points_Acum);
point1 = centerline.GetPoint(points_Acum+1);
distance = np.sqrt(vtk.vtkMath.Distance2BetweenPoints(point0,point1));
tangent[0] += (point1[0] - point0[0]) / distance;
tangent[1] += (point1[1] - point0[1]) / distance;
tangent[2] += (point1[2] - point0[2]) / distance;
weightSum += 1.0;
tangent[0] /= weightSum;
tangent[1] /= weightSum;
tangent[2] /= weightSum;3. Cut the mesh with the plane
Whit the calculated tangent vector, I create a plane with which I´m gonna cut the mesh. A great tool I found that saved me from diving deeper into vtk is the vedo library(https://vedo.embl.es/), it is a higher-level way to use vtk. I´m also using this library for the visualizations, as it is a great high-level tool, but remember all its functionalities are coded with VTK
from vedo import *
cutplane = Grid(pos = centerline.GetPoint(points_Acum),s = (5,5),
res = (150,150), normal = tangent).cutWithMesh(msh).triangulate()“Grid” is a vedo class and uses “vtk.vtkPlaneSource()”, as you can see in the vedo source code.
It’s a straightforward process. First, we create a plane using the centerline point coordinates as the origin and the point tangent as the plane normal. The size of the plane is not crucial at this stage; I’ll delve into that explanation later. Utilizing the cutWithMesh() method from the Grid class, we extract the section of the plane that lies inside the mesh.
Frequently, the plane may be too large, causing it to intersect the mesh in multiple areas. Automating the plane size presents a challenge, as each mesh varies, demanding a specific size for accuracy. Thus, relying on the incorrect cross-section area can result in inaccuracies, making it necessary to introduce an additional step for radius calculation.
I first panicked when I ran into this issue, but fortunately, there is an easy fix with a vtk filter. The idea is, to choose the region that is closest to the centerline point at which we are centering the cut plane, and the vtk filter that does so is vtkConnectivityFilter, with the extraction mode set to ClosestPointRegion. Below is the function I wrote
def get_closest_region(cutplane, point):
connectivityFilter = vtk.vtkConnectivityFilter()
connectivityFilter.SetExtractionModeToClosestPointRegion();
connectivityFilter.SetInputData(cutplane._data);
connectivityFilter.SetClosestPoint(point);
connectivityFilter.Update();
m = Mesh(connectivityFilter.GetOutput())
pr = vtk.vtkProperty()
pr.DeepCopy(cutplane.property)
m.SetProperty(pr)
m.property = pr
m.SetOrigin(cutplane.GetOrigin())
m.SetScale(cutplane.GetScale())
m.SetOrientation(cutplane.GetOrientation())
m.SetPosition(cutplane.GetPosition())
vis = cutplane._mapper.GetScalarVisibility()
m.mapper().SetScalarVisibility(vis)
return mThe final complete loop looks like this:
for j in range(centerline.GetNumberOfCells()):
numberOfCellPoints = centerline.GetCell(j).GetNumberOfPoints();
for i in range (numberOfCellPoints):
tangent = np.zeros((3))
weightSum = 0.0;
if (i>0):
point0 = centerline.GetPoint(points_Acum-1);
point1 = centerline.GetPoint(points_Acum);
distance = np.sqrt(vtk.vtkMath.Distance2BetweenPoints(point0,point1));
##vector between the two points divided by the distance
tangent[0] += (point1[0] - point0[0]) / distance;
tangent[1] += (point1[1] - point0[1]) / distance;
tangent[2] += (point1[2] - point0[2]) / distance;
weightSum += 1.0;
##tangent line with the next point (not calculated at the last one)
if (i<numberOfCellPoints-1):
point0 = centerline.GetPoint(points_Acum);
point1 = centerline.GetPoint(points_Acum+1);
distance = np.sqrt(vtk.vtkMath.Distance2BetweenPoints(point0,point1));
tangent[0] += (point1[0] - point0[0]) / distance;
tangent[1] += (point1[1] - point0[1]) / distance;
tangent[2] += (point1[2] - point0[2]) / distance;
weightSum += 1.0;
tangent[0] /= weightSum;
tangent[1] /= weightSum;
tangent[2] /= weightSum;
##cut the plane with the mesh
cutplane = nGrid(pos = centerline.GetPoint(points_Acum),s = (5,5), normal = tangent).cutWithMesh(msh).triangulate()
##create mesh with the section polygons
vert = cutplane.points()
faces = cutplane.faces()
m = Mesh([vert, faces])
m = get_closest_region(m, centerline.GetPoint(points_Acum)) #keep only one region if it cuts the mesh multiple times
vert = m.points()
faces = m.faces()
area = m.area()
ceradius = np.sqrt(area / np.pi)
radius_array.append(ceradius)
line_array.append(j)
points_Acum += 1 #keep track of points relative to the whole centerline not just the specific branchGreat! Now we have successfully obtained the radius for the entire centerline.
To recap, we´ve seen a basic overview of both VTK and VMTK with detailed practical usage for a real scientific paper. After reading this blog post, I hope you are not scared to create your own scripts using these valuable tools
References:
[2] K. J. Weiler. Topological Structures for Geometric Modeling. PhD thesis, Rensselaer Polytechnic Institute, Troy, NY, May 1986.
[3] https://examples.vtk.org/site/VTKBook/
[4] http://www.vmtk.org/vmtkscripts/vmtknetworkextraction.html