Decomposition of Velocity Field Along a Centerline Curve Using Frenet-Frames: Application to Arterial Blood Flow Simulations

This paper presents a novel method for the evaluation of three-dimensional blood-flow simulations based, on the decomposition of the velocity field into localized coordinate systems along the vessels centerline. The method is based on the computation of accurate centerlines with the Vascular Modeling Toolkit (VMTK) library, to calculate the localized Frenet-frames along the centerline and the morphological features, namely the curvature and torsion. Using the Frenet-frame unit vectors, the velocity field can be decomposed into axial, circumferential and radial components and visualized in a diagram along the centerline. This paper includes case studies with four idealized geometries resembling the carotid siphon and two patient-specific cases to demonstrate the capability of the method and the connection between morphology and flow. The proposed evaluation method presented in this paper can be easily extended to other derived quantities of the velocity fields, such as the wall shear stress field. Furthermore, it can be used in other fields of engineering with tubular cross-sections with complex torsion and curvature distribution.


Introduction
In the last two decades the number of research studies in computational hemodynamics rose substantially, due to the ever-growing computational resources and their acceptance by the medical practitioners. Since cardiovascular vessel malformations, such as narrowings (stenoses) or enlargements (aneurysms) are still on the top of the list of death causes worldwide [1], the research of these diseases is of great importance. The past two decades has seen a rapid development of computational [1][2][3][4] and imaging tools [5][6][7][8] that can help researchers avoid the tedious pre-processing work from medical imaging to the solution of the numerical problem. More recently, the research community stepped forward by conducting international research collaborations to study these computational tools and the decisions that has to be made during modeling to reach a consensus about their correct usage [9][10][11][12][13]. Computational fluid analysis has been used in this field to reach several goals; risk prediction of intracranial and aortic aneurysm rupture [14][15][16]; severity analysis of carotid and coronary stenosis [17][18][19]; modeling the formation [20][21][22][23][24] of the disease or the process of pathogenesis either with follow-up studies or numerical growth modeling [25][26][27]; assessment of treatments or planning the procedure [28][29][30][31][32]. One important aspect common among most of these investigations is the consideration of the local geometrical and morphological features, primarily as they are -through the locally emerging flow field -linked to the local mechanobiological sensing [33][34][35][36]. Vascular Modelling Toolkit (VMTK), a software package is one of the widely utilized tools for geometric analysis of vessel segments [23,[37][38][39]. Among other functionalities, the calculation of the centerline provides a number of useful quantities about the global morphological properties of the given vessel section such as the unit vectors of local Frenet-frames and the curvature and torsion metric. There are two main uses of geometrical information in blood flow analysis: 1. to correlate the flow field to morphological features [18,40,41], or 2. to utilize geometrical information in the pre-and post-processing of the simulations [19,42,43].
The latter is unequivocal, but the former is often challenging as the simulated flow field in the global coordinate system is not easily reconcilable with the local geometric and morphological properties of the vessel segment. These properties (curvature and torsion) are based upon the locally varying Frenet-frames of the centerline and their unit vectors. In this paper, we present an evaluation methodology based on the centerline calculations by VMTK for the decomposition of the velocity field in local coordinates using the Frenet-frames. The main purpose of this study is to develop a framework in which the velocity field (and the wall shear stress field) lies in the same system of reference as the geometrical properties.
The paper is structured as follows: In Section 2 the mathematical background of the Frenet-frames is introduced, including the definition of the curvature and torsion. Then the proposed decomposition technique for the velocity field is given, followed by the transformation procedure of the geometric coordinates of slices that will be used later in our evaluation approach. Section 2 is closed by a summary on the Computational Fluid Dynamics (CFD) simulations (Subsection 2.4). In Section 3 the method is demonstrated and discussed on several artificial cases and two patient-specific geometries. Finally, the paper ends with the Conclusion.

Centerlines and Frenet-frames
As mentioned in the Introduction, the procedure starts with a centerline calculation in VMTK. Without further explanations, a centerline by definition lies on the Voronoi-diagram of the triangulated surface geometry. In VMTK, it is achieved by solving the Eikonal equation numerically with the Fast-Marching method [44], which is essentially the solution of a minimization problem for the shortest path between a source and a target point through the Voronoi diagram. Next, a local coordinate system is created at every point of the centerline by the Frenet-Serret formulas [5]. These vectors form an orthonormal basis at each point on the curve in the three-dimensional space. The curve c(s) is parameterized by the arc length. The three unit vectors are referred to as tangential t, normal n and binormal b ones, with the following definitions: These calculations are automatically carried out by the VMTK package, where the derivatives are calculated using finite differences. Since this method uses second and later third derivatives of the centerline curves, a smoothing (by cubic splines) process is required to create a continuous smooth centerline from the noisy curve defined by the Voronoi-diagram.

Morphological properties
From the smoothed centerline and its Frenet-frames, two morphological parameters can be quantified along the centerline: the curvature (κ) and torsion (τ) (see Eqs. (4) and (5)). The curvature is the inverse of the radius of a fitted circle to the centerline in the osculating plane, defined by the tangential and normal vectors of the Frenet-frame at a given point. Torsion (also called second curvature) is the rate of twisting out of the osculating plane, or the speed of rotation of the binormal vector. These quantities can give a detailed insight into the morphological characteristics of the vessel and are influential in the formation of the flow field.

Velocity field decomposition
The first step in our method is to define plane slices that are perpendicular to the centerline. In these planes the field variables can be defined in an arbitrary point by interpolation from the mesh points.

Axial and secondary components
While several conclusions can be drawn about the vascular system using a global velocity field, a local decomposition is desirable to assess the relation of the axial and secondary flow features. Before going into detail, a terminological clarification is in order here. In the global coordinate system, the component showing in the direction of the tube is called tangential component. In the local system, to avoid confusion, this component is called axial component. Perpendicular to this, we call the component in this plane of the cross section the secondary component ( Fig. 1 a), c)). The secondary component is further divided into the circumferential and radial components, respectively, representing local cylindrical coordinates defined below ( Fig. 1 b), d)). In our application of the Frenet frames, the tangential unit vector was used to decompose the global velocity field into the axial component, then the secondary component, as follows. According to Eqs. (6) and (7)) an arbitrary velocity vector can be decomposed into the axial flow v ax and the secondary flow components v sec defined using the tangential part of the Frenet unit vectors t(p) in the corresponding p centerline point. Then by Eq. (7) the secondary flow vector v sec is calculated by subtracting the axial flow vector v ax from the velocity vector v.

Radial and circumferential components
Furthermore, the secondary flow components can be further separated into a radial and circumferential component by a local coordinate system transformation by switching to a polar representation around the centerline point ( r p ), lying in the given cross section. To calculate these properties, first a radial unit vector r rad is defined by Eq. (8) between an arbitrary point (r) on the plane and ( r p ). Then, using Eqs. (9) and (10), the secondary velocity vector can be decomposed into a radial v rad and circumferential v cir components, respectively.

Transformation of the geometric coordinates
Having discussed, how to construct the velocity components in Frenet-frames along the centerline, the final part of our evaluation procedure is to organize the data into a comprehensible form. In this final step the slices are transformed to the origin of the global coordinate system, offset, and then rotated to match the Frenet-frame unit vectors with the global axis directions (n -x, b -y, t -z as shown in Fig. 2.). To express it in an illustrative way: the vessel is "straightened". After transforming the slices to the global coordinate system, the slices were partitioned into an inner and outer section, based on the transformed Frenet-frames. In the global coordinate system, an arbitrary point on a slice lies in the outer region if its x coordinate value is less than zero, therefore, those points where x > 0, are in the inner partition. The inner and outer partition represent the inside and outside part of a bend. It follows from the above definition that if there is an inflection point on the centerline (as it is often the case in practice) then the coordinate system makes a sudden 180° turn. These partitions were averaged for the calculated flow components along the arterial section. These averaged vector components (axial, radial, circumferential) provide a basis for understanding the influence of the geometry and morphology on the flow structure.
Here, we want to remark that this partitioning can serve as example but different and more complex partitionings can be implemented easily.

Computational fluid dynamics
Case studies were used to demonstrate the capability of the previously described evaluation approach. First, four idealized geometries were created for the CFD simulations that distantly resemble an internal carotid artery (ICA) section with the so-called carotid siphon. Each artificial vessel consists of two bends with arc radii of 4.2 and 3.6 mm, respectively. The diameter is 4 mm along the whole geometry, while the inlet and outlet sections are 30 mm long which ensures proper velocity profile development [45]. Fig. 3 a) shows that the only difference between the cases is the out-of-plane rotation at the beginning of idealized carotid siphon, separating the two bends. Apart from the section separating the two bends, curvature and torsion along the centerline are the same for all geometries. At that intermediate section a minor change in curvature but a sudden change of torsion occurs. Steady-state simulations were prepared within the ANSYS software package. The numerical solutions were carried out with the finite volume solver of ANSYS CFX. The unstructured numerical meshes consisted of about 2.5 million cells with seven prismatic layers adjacent to the wall. Second order discretization was used with a convergence criterion of 10 −6 but the simulations were only stopped when a monitored physical quantity  converged to constant value. A Newtonian approach was used for the blood flow with a density and viscosity of 1055 kg/m 3 and 3.4 mPas, respectively. A parabolic velocity profile with a mean velocity of 0.121 m/s was imposed on the inlet corresponding to a Reynolds number of 150. 0 Pa static pressure was prescribed on the outlets and the walls were assumed to be rigid.
Furthermore, two patient-specific geometries (Case01 and Case02) were investigated to demonstrate the real capability of the proposed evaluation technique, depicted in Fig. 3 b) and c). These geometries are the part of an ongoing investigation on side-wall aneurysm initiation, in which the aneurysmal bulge is virtually removed from the vessel by an objective algorithm. The simulations associated with these geometries were transient but for the demonstration of the evaluation approach only a certain time instant is selected in this paper. The numerical meshes contain about three million cells. The inlet boundary condition is a time-varying parabolic velocity profile.
The waveform of the heart cycle is the same as in [30] in a non-dimensional form and is scaled according to the inlet diameter as proposed by [46].
Flow rate splitting boundary conditions are set on the outlets according to a modified Murray's law [47], closed by a 0 Pa static pressure prescribed only on the smallest outlet. Adaptive time-stepping was used to ensure a proper Courant number criterion and a second order discretization scheme for the time integration. After the solution the result files are exported into a format that is suitable for ParaView and our code.

Summary of the algorithm
Before moving on to Section 3, the algorithm is summarized, as follows. Based on the results of a CFD simulation and a centerline geometry computation in VMTK, plane slices perpendicular to the centerline are defined at every point of the centerline. The field variables in these planes are calculated by interpolation from the mesh points. The velocity components are calculated with Eqs. (6)- (10). At this point the flow field, although decomposed into components along a centerline, is defined in terms of the global coordinate system. Afterwards, slice coordinates are transformed, so that the Frenet reference system belonging to each slice will align with the global coordinate systems as shown in Fig. 2 c).

Results and discussion
As pointed out in Subsection 2.6, case studies with a series of idealized geometries and two patient-specific ones were completed to demonstrate the capability of the method. In order to introduce some order into the vast amount of data, also to determine the relative importance of the secondary flow in a certain vessel segment, the velocity ratio was defined as the ratio of the secondary and axial velocity magnitude:

Idealized geometry
The results of the four idealized geometries can be seen in Fig. 4. Each column represents one variant of the idealized geometries starting from the in-plane geometry to the 90° out of plane one. The averaged velocity fields plotted along the centerline length with yellow and blue curves symbolize the inner and outer segments, respectively. The rows represent the axial, secondary, circumferential, and radial components of the velocity field in this order, while the fifth one depicts the aforementioned velocity ratio. The last row displays the morphological properties, the curvature and torsion distributions with a shaded gray area and green line, respectively. In order to ease the understanding, the same format is used for the latter figures depicting the results of the realistic cases. The last row of Fig. 4 demonstrates that curvatures along the centerline are almost identical in all cases with a slight increase between the two bends with an increasing out-of-plane angle. Torsion, on the other hand, has an abrupt change at this location. Since the centerline of the first planar geometry is not C 3 continuous at the intermediate point between the two bends, the torsion here has a singularity and is practically infinite. Similarly, the torsion has a singular but not infinite value in the other geometries as well, decreasing with an increasing out-of-plane angle.

Overview of the flow features
In order to understand the effect of morphological change on the velocity field in this idealized carotid siphon, an overview is granted on the emergent flow phenomenon due to the subsequent bends (first diagram in Fig. 4). Due to inertial forces, the velocity profile starts to distort before the bend at the end of the inlet section and leans toward the inner segment where the inner axial velocity increases in the first shaded region, before reaching the outer wall. Here, the outer axial velocity gradually reaches and eventually overcomes the inner one (second shaded region). In the intermediate section between the second and third shaded region, the coordinate system switches direction because of the inflection point, which also means a torsion discontinuity and a change of the inner and outer sides at this point. Subsequently, the outer axial velocity takes over between the third and fourth region again as the same phenomenon occurs as in the first bend. In the fourth region the velocity profile leans again toward the outer part before eventually converging back to a symmetric profile. The effect of the aforementioned singular point can be observed on all the velocity components.

Effect on the velocity components
After the aforementioned point the axial velocities in the inner segment increase with the torsion, while the outer velocity decreases. The second row of Fig. 4 shows that up until 45°, with decreasing torsion the secondary velocity at the second bend increases in the inner segment and much stronger decreases in the outer segment.
In the outer segment, at 90°, the secondary velocity still increases right after the turning point, but the maximum value slightly decreases in the second bend, compared to the 45° case.
By further decomposition, the results show that the radial velocity component at 90° decreases more than the circumferential one. It indicates that a weaker transition from one bend to the next (i.e. lower torsion) distort the velocity field less as it can be seen in Fig. 5. All the studied parameters show a certain trend between 0 and 45°. This trend, however, seems to change between 45° and 90° and to find the point of change, as well as its cause need a better resolution in this range of angles. The explanation of the above observations needs a deeper analysis with the consideration of the various forces, arising in curved coordinate systems. This is beyond the scope of the present paper. Only a slight, almost negligible change can be observed in the velocity ratio in the inner segment. However, the decreasing torsion had a significant and intensifying effect in the second bend of the outer segment, almost tripling the maximum value from around 0.2 to about 0.6.

Patient-specific cases
Although the morphology of the patient-specific cases is significantly more complex along the carotid siphon, some similarities can be discussed in Subsection 3.2 between the patient-specific and idealized cases, since a patient-specific carotid siphon also resembles an S-bend with an inflection point.
The results of the patient-specific cases can be seen in Figs. 6 and 7. In both cases, a similar inversion between the components of the inner and outer segment can be observed, where the torsion has a spike, as depicted at the largest occurrence with the first dashed line in the figures.
This aforementioned location also marks the beginning of the carotid siphon segment with the largest curvature (up to the second dashed line in the figures) throughout the ICA vessel. Furthermore, significant variations in the flow components evolve after this point, in both cases. In both cases, the inner segment dominates for the axial and circumferential components until the apex point in the curvature. Subsequently, at the distal side of the bend the outer segment starts to dominate the axial component as the flow impinges on to the outer wall of the vessel. Up to this point along the centerline, clear similarities could be identified between the two patient-specific and the idealized cases, as the internal carotid artery shares a number of morphological landmarks in every person. However, distal to carotid siphon section, leading to one of the junction points of the Willis circle (middle and anterior cerebral artery) a more in-depth characterization is needed with further cases which we postpone to the next paper. Here we want to remark, that quantitative comparisons will not be given here, as the purpose of this paper is to present the evaluation methodology.

Potential applications
In this current study, we proposed a methodology for the decomposition of the velocity field into local coordinate systems defined at the points of the centerline. Several reports have shown the influence of the geometry on the emerging flow field [18,22,24,42,48,49]. However, systematic relationships between flow and geometry have not been established. In this paper, we demonstrated that by using the Frenet-frames, the flow field can be analyzed in a fused reference system following the vessel morphology. The results of the idealized and patient-specific cases indicate that the quantitative effect of morphological changes -such as the curvature and torsion along the centerline -can be exposed with the proposed decomposition technique. Furthermore, another convenience of depicting the flow field with these components is the direct linkage with mechanobiological sensing. The directions of the axial, circumferential and radial components are in line with main types of directional load sensing capacity of a vessel wall [50].
In this paper, only a simple decomposition into an inner and outer segment was presented to show the capabilities of the evaluation methodology. In general, a more detailed analysis can be easily implemented by using a different partition strategy. The purpose of this paper is not the formulation of any hypothesis about the generation of aneurysm, but to introduce a toolkit.
Here only the effects of time dependence were considered but time-dependent analysis can be easily carried out in the form of a three-dimensional waterfall or a spectrogram-like diagram for the purpose. Moreover, the vector decomposition could be used for the evaluation of the wall   shear stress field and other quantities as well. The only difference being that on the wall, the vector quantities will only decompose into a circumferential and an axial component. Finally, the proposed method is not limited to the ICA vessel segment, it can also be used for other tubular geometries with complex morphology.

Remark on the morphological properties
Although the methodology presented in this paper has satisfied our objectives, it has certain points that need to be considered cautiously. As mentioned before, to calculate the torsion and the curvature, the curve is required to be C 3 continuous. The centerline based solely on the Voronoi diagram is noisy by definition and needs smoothing. With sufficient resolution and smoothing, the Frenetframes will follow the geometry continuously. However, the smoothing procedure demands caution as it can introduce a minor difference in the centerline location in order to provide better differentiability properties. Another source of uncertainty may arise at the inflection point of an S-bend like section with a sudden change in Frenet frame orientation. The switch is highly dependent on the centerline smoothing and resolution, therefore, special care has to be taken there as well.
Since in the VMTK separate whole centerlines are calculated from the inlet to each opening, separate evaluations can be performed for each branch vessel, if needed. Currently our method uses a single centerline picked out after the centerline calculation. Caution has to be exercised during an evaluation at the bifurcation regions, as different Frenet-frames can be identified with centerlines, leading to different vessel branches.
A final remark concerns the fact that the flow features are not only influenced by the local morphology but also the evolution of the flow upstream of the region of interest. In the artificial geometries the upstream region was therefore identical and this way we could make the study independent of the upstream effects. The patient-specific cases were too small in number to make a systematic study on the upstream effects.

Conclusions
This project was undertaken to design an evaluation framework to integrate the hemodynamic field into a frame of reference given by the centerline geometry. The study establishes a quantitative framework for finding the connection between geometric and hemodynamic variables. The partitioning of the segments by the Frenet frames gives a more complete overview of the emerging flow features at given sections of the vessel. The technique offers a wealth of further opportunities for analysis. For example, it might help land-marking certain points of the artery that are prone to disease formation and potential progression.