Applicable version(s):
Latest stable release & Development version
This tutorial demonstrates how to analyze the 3D geometry of a vascular segment.
Before proceeding, make sure you went through the centerlines tutorial and the branch splitting tutorial.
The theory behind this tutorial can be found in
- Antiga L, Steinman DA. Robust and objective decomposition and mapping of bifurcating vessels. IEEE Transactions on Medical Imaging, 23(6), June 2004.
- Piccinelli M, Veneziani A, Steinman DA, Remuzzi A and Antiga L. A framework for geometric analysis of vascular structures: application to cerebral aneurysms. IEEE Transactions on Medical Imaging 2009 Aug;28(8):1141-55.
You can see it in action in
- Thomas JB, Antiga L, Che S, Milner JS, Hangan Steinman DA, Spence JD, Rutt BK and Steinman DA. Variation in the carotid bifurcation geometry of young vs older adults: Implications forgeometric risk of atherosclerosis. Stroke, 36(11), Nov 2005.
- Lee SW, Antiga L, Spence JD and Steinman DA. Geometry of the carotid bifurcation predicts its exposure to disturbed flow. Stroke, 39(8): 2341-2347, Aug 2008.
Overview
Let’s assume we have a surface representing a branching vascular segment, with multiple branches each with N vessels departing from the parent vessel. As in the branch splitting tutorial, let’s assume the network has a tree-like topology.
Using vmtk, you should have been able to compute centerlines, to identify bifurcations and branches on the centerlines and to clip the surface according to the branches. Now we’ll face the problem of quantifying geometric features of the vascular segment, those associated to bifurcations, such as bifurcation planes and bifurcation angles, and those associated to branches, such as curvature and torsion.
However, before characterizing geometry, we introduce a few geometric concepts which the subsequent analysis relies on.
Centerline attributes
Centerlines are defined by the position of their points. In order to define the position of a point along the centerline or to quantify the angular position of a point on the surface around the centerline, it is necessary to equip centerlines with two attributes, Abscissas and Normals.
Abscissas are easy to define: they measure the distances along a line. The point with abscissa 0 is not necessarily the first point of the centerline, but can be chosen to represent a landmark along a centerline, such that, if we have a population of models, all the abscissas will be referred to the same anatomical location. In absence of such location, abscissas can be generated relative to the starting point and can be offset to a different location at a later time.
Measuring the angular position of a generic point around a line in 3D is not so trivial. It depends on the construction of a reference system along the line that defines the reference orientations in each point. Unlike a surface, on which defining a normal orientation is quite intuitive, defining a normal orientation on a line is not straightforward, and several approaches exist. The most popular is the Frenet reference system, in which the normal to a line points towards the center of the osculating plane, and the binormal is the normal to the osculating plane. Although this system has several nice properties, it has also a few disadvantages: the normal can experience abrupt changes in its orientation whenever the osculating plane changes, and can be undefined where the line is straight. In contrast, our intuitive requirements for a reference normal along a line is that the normal has to follow the line geometry avoiding too much twisting around the line. This is offered by parallel transport reference systems: given a starting reference system (the analogous to the 0 abscissa point), the next reference system is generated by moving along the centerline of a small amount and rotating the system on the osculating plane of an amount equal to the change in the orientation of the line tangent. This operation is repeated until all points are equipped with a reference system. Since reference systems are rotated within the osculating plane no artificial torsion is introduced proceeding along the centerline. As with abscissas, the choice of the initial reference system influences the orientation of the resulting normals. In case of a population of vessels, if a common orientation is chosen at a landmark, orientations can be compared across the population. If a landmark is not available, the starting orientation can be given arbitrarily and the reference systems can eventually be rigidly rotated to conform to a landmark orientation at a later time.
The script that computes abscissas and parallel transport normals is vmtkcenterlineattributes. It computes abscissas and normals relative to the first point of each centerline (they can be offset to landmarks at a later time). As such, it should be called before vmtkbranchextractor, otherwise each tract into which centerlines are split will have an abscissa starting from 0 at the first point. The usage is simply
vmtkcenterlineattributes -ifile foo_cl.vtp -ofile foo_clat.vtp
Before learning how to offset attributes at a landmark, let’s first define what landmarks are in vmtk and how they are defined.
Bifurcation reference systems
Since we are dealing with bifurcating vessels, it is quite natural that the position and orientation of bifurcations become important features of a vascular segment. If one is looking for a point where to place the 0 abscissa and the normal, the position of a bifurcation and its spatial orientations are good candidates for such a job, also because it is quite common to find the same bifurcations across a population of subjects, while the geometry of branches can change remarkably (indeed, several parts of the vasculature present variability event in the configuration of bifurcations, but this is a problem that can be dealt with segment-specific ad-hoc strategies).
Although at a first thought defining the position and orientation of a bifurcation is an easy task, formulating an objective criterion to robustly perform such operation given the surface of a bifurcating vessel is not trivial. Fortunately, we already learnt how to compute centerlines, decompose them into branches and split the surface accordingly. We now take a step further and define a bifurcation reference system.
In order to define a reference system we need the position of the origin and the orientation of three normals. Well, two normals, if we consider that the third can be obtained by a cross product. In our bifurcation reference system, the first normal can be taken as the normal defining the bifurcation plane, i.e. the plane onto which the daughter vessels depart from the parent vessel, while the second normal (which we’ll call upnormal) describes the direction pointing from the parent vessel to the daughters.
We here consider a single bifurcation: the extension to a more general case is rather straightforward. From the centerline splitting and grouping phase we’re left with a first group that ends with the second reference points (one for each tract), a second group representing the bifurcation whose tracts start with the second reference point and end with the first, and two last groups representing the daughter branches that start with the first reference points (where each centerline exits the other’s tube). We define the bifurcation reference system the following way:
- the origin of the bifurcation is defined as the barycenter of the four reference points weighted by the surface of the maximum inscribed sphere defined on the reference points. The reason of the weighting is that small branches have less impact on the position of the bifurcation origin
- the normal to the bifurcation plane is defined as the normal to the polygon defined by the four reference points computed in the bifurcation origin
- the bifurcation upnormal is the similarly weighted average vector between the vectors pointing from the second to the first reference point on each centerline
Let’s see how to generate reference systems with vmtk scripts. Given the centerlines, we first split them into branches and then compute the reference systems:
vmtkbranchextractor -ifile foo_cl.vtp -radiusarray@ MaximumInscribedSphereRadius --pipe vmtkbifurcationreferencesystems -ofile foo_rs.vtp
A note: see how in the pipe above I made use of argument pushing, described at the end of the basic PypeS tutorial. The radius array parameter is taken in input but normally not sent in output along the pipe by vmtkbranchextractor, and so vmtkbifurcationreferencesystems would not be able to use it. By appending @ to -radiusarray, I forced its value to be pusehd along the pipe and thus to be used by vmtkbifurcationreferencesystems. Thanks to Barry Hogan for pointing this aspect out.
Back to the tutorial, the file foo_rs.vtp will contain a set of points, one for each bifurcation, located at the bifurcation origin, with two point data arrays defined, Normal and UpNormal.
If you need the positions and orientations in a human readable format, use the pointdata format (extension dat)
vmtkbranchextractor -ifile foo_cl.vtp -radiusarray@ MaximumInscribedSphereRadius --pipe vmtkbifurcationreferencesystems -ofile foo_rs.dat
Offsetting centerline attributes to bifurcations
As anticipated, centerline attributes (i.e. abscissas and parallel transport normals), which were computed relative to the first point of the centerline (not a very meaningful landmark), can be offset to bifurcations. This way, a population of vessels will have the 0 abscissa corresponding to the origin of the same bifurcation, which sounds like a reasonable thing to do. Similarly, parallel transport normals at the bifurcation will be oriented like the bifurcation plane normal.
To do this, you have to first choose a bifurcation to use as reference. Bifurcations are identified by their group id: to inspect group ids, use vmtkcenterlineviewer
vmtkbranchextractor -ifile foo_cl.vtp -radiusarray@ MaximumInscribedSphereRadius --pipe vmtkcenterlineviewer -cellarray GroupIds
Once you identify the reference bifurcation, use vmtkcenterlineoffsetattributes, which will take the attributes computed by vmtkcenterlineattributes and offset them to the reference bifurcation
vmtkcenterlineattributes -ifile foo_cl.vtp --pipe vmtkbranchextractor -radiusarray@ MaximumInscribedSphereRadius --pipe vmtkbifurcationreferencesystems --pipe vmtkcenterlineoffsetattributes -referencegroupid 1 -ofile foo_clatof.vtp
In the above example, -referencegroupid 1 means that the group id of the bifurcation we chose as a landmark is 1 (0 will likely be the first tract). In case we only have one bifurcation in the segment we can omit the option and the only bifurcation in the dataset will be taken as the reference.
One interesting option is -replaceattributes: if 1, the arrays containing the original abscissas and normals will be replaced by the offset ones, otherwise those will be kept and new arrays (with new names) will be created. This has an implication for using the script within a pipe. If -replaceattributes is 1 (the default) the script will be easier to pipe, otherwise the name of the new arrays will have to be explicitly specified to the following scripts in the pipe.
Note that, while vmtkcenterlineattributes expects to be run before vmtkbranchextractor, vmtkcenterlineoffsetattributes runs on centerlines that are already split and grouped in branches, and therefore its output will be as such.
Computing bifurcation geometry
We will now see how to compute the geometry of a bifurcation, in terms of the description of how the parent and daughter vessels are geometrically related. First we will describe such relationships are quantified in the case of a generic bifurcation with N branches. Then we will derive from the general approach more intuitive quantities for bifurcations with one parent and two daughter branches.
Each branch, parent or daughter, arrives or leaves the bifurcation region with a certain orientation. Such orientation can in general be decomposed in in-plane and out-of-plane components, where the plane is the bifurcation plane. The geometry of the bifurcation will therefore be described by the in-plane and out-of-plane angles that each branch forms with the bifurcation reference system. Taking the upnormal as an in-plane reference orientation, each branch at the bifurcation is described by its in-plane angle with the upnormal direction and by its out-of-plane angle with the bifurcation normal.
From such information, it is easy to compute branch angles, which we define as the angles formed by parent vessels with each of the daughter vessels.
Operatively, let’s start by showing the pipe:
vmtksurfacereader -ifile foo.vtp --pipe vmtkcenterlines -seedselector openprofiles --pipe vmtkbranchextractor --pipe vmtkbifurcationreferencesystems --pipe vmtkbifurcationvectors -ofile foo_bv.vtp
The pipe generates a vtp file containing the data on the bifurcation vectors with the following layout: one point each bifurcation reference point, equipped with the following point data:
- BifurcationVectors
- InPlaneBifurcationVectors
- OutOfPlaneBifurcationVectors
- InPlaneBifurcationVectorAngles (the angle between the InPlaneBifurcationVectors and the bifurcation UpNormal, in radians, from -pi to pi, zero for a UpNormal oriented vector, positive in the clockwise direction with respect to the bifurcation Normal)
- OutOfPlaneBifurcationVectorAngles (the angle between the BifurcationVectors and the bifurcation plane, in radians, positive if the OutOfPLaneBifurcationVector is directed as the bifurcation Normal)
- BifurcationVectorsOrientation (flag accounting for the role played by the branch in the bifurcation, 0 for upstream, 1 for downstream the bifurcation)
- GroupIds (the groupId of the branch described by the vector)
- BifurcationGroupIds (the groupId of the bifurcation)
You can visualize all the above in Paraview by loading the vtp file and running it through a Glyph filter (which will show you arrows representing the vectors), or by selecting a Spreadsheet view and viewing the row numbers in tabular format.
As usual, you can get the above information in a text file using the pointdata format
vmtkcenterlines -ifile foo.vtp -seedselector openprofiles --pipe vmtkbranchextractor --pipe vmtkbifurcationreferencesystems --pipe vmtkbifurcationvectors -ofile foo_bv.dat
The above data should be sufficient to extract whatever geometric characteristic of the bifurcation. For instance, how do you compute the bifurcation angle (usually referred to as the angle between the daughter branches of a Y-shaped bifurcation?): just take the difference between the in-plane angles of the two daughter branches. And what about the ICA angle in the Stroke 2005 paper cited at the top of the page? Take the supplementary angle of the difference between in-plane angle of the CCA and the ICA. A quick look at the figures of that paper will clarify these suggestions.
Computing centerline geometry
The geometry of a line in space is generally described by curvature (the inverse of the radius of the local osculating circle) and torsion (the amount by which the osculating plane rotates along the line). A quick search on Wikipedia or the Piccinelli et al paper referenced above will tell you more about the formal expressions. Curvature and torsion are tightly linked to the definition of the Frenet line frame, constituted by a tangent, a normal (pointing towards the center of the osculating circle) and the binormal (the normal to the osculating plane). All these quantities are computed in vmtk by
vmtkcenterlinegeometry -ifile foo_cl.vtp -ofile foo_clgm.vtp
Note that the input centerlines don’t have to be split in branches. Since the computation of the above quantities depend on first, second and third derivatives of the line coordinates, and since such derivatives are approximated using a simple finite difference scheme along the line, it is very likely that such derivatives will be affected by noise that is not appreciable when looking at the line itself. For this reason, it is possible to instruct the script to first run the centerline through a Laplacian smoothing filter first and then compute the derivatives and the related quantities. This is done by
vmtkcenterlinegeometry -ifile foo_cl.vtp -smoothing 1 -ofile foo_clgm.vtp
You can control the number of iterations and the relaxation factor of the Laplacian filter using
vmtkcenterlinegeometry -ifile foo_cl.vtp -smoothing 1 -iterations 100 -factor 0.1 -ofile foo_clgm.vtp
You can also choose to output the smoothed version of the centerlines (i.e. not just using it to compute derivatives):
vmtkcenterlinegeometry -ifile foo_cl.vtp -smoothing 1 -iterations 100 -factor 0.1 -outputsmoothed 1 -ofile foo_clgm.vtp
In addition to curvature, torsion and Frenet frames, the script also outputs centerline tortuosity (the ratio between the centerline length and the distance of the line endpoints).
In case you are dealing with split centerlines, it could be useful to compute these quantities for each branch, not for all centerline tracts composing the branch individually, as it would happen using vmtkcenterlinegeometry. For this task, you can use vmtkbranchgeometry, which will compute an abscissa-based group average to provide the same branch-wise geometric quantities. For instance:
vmtkcenterlines -ifile foo.vtp -seedselector openprofiles --pipe vmtkbranchextractor --pipe vmtkbranchgeometry -ofile foo_clcg.vtp
The optional parameters are the same as for vmtkcenterlinegeometry, except that you cannot output the smoothed version of the centerlines.
