# Modeling and hexahedral meshing of arterial networks from centerlines

Méghane Decroocq1,2,4,5,6 1 CREATIS, Université Lyon1, CNRS UMR5220, INSERM
U1206, INSA-Lyon, 69621 Villeurbanne, France
2 LIRIS, CNRS UMR 5205, F-69621, France
3 Ecole Centrale de Lyon, France
4 ELyTMaX IRL3757, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard
Lyon 1, Tohoku University, 980-8577, Sendai, Japan
5Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku,
Sendai, Miyagi 980-8577, Japan
6 Graduate School of Biomedical Engineering, Tohoku University, 6-6 Aramaki-
aza-aoba, Aoba-ku, Sendai, Miyagi 980-8579, Japan
Corresponding Author<EMAIL_ADDRESS>Carole Frindel1,5 1
CREATIS, Université Lyon1, CNRS UMR5220, INSERM U1206, INSA-Lyon, 69621
Villeurbanne, France
2 LIRIS, CNRS UMR 5205, F-69621, France
3 Ecole Centrale de Lyon, France
4 ELyTMaX IRL3757, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard
Lyon 1, Tohoku University, 980-8577, Sendai, Japan
5Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku,
Sendai, Miyagi 980-8577, Japan
6 Graduate School of Biomedical Engineering, Tohoku University, 6-6 Aramaki-
aza-aoba, Aoba-ku, Sendai, Miyagi 980-8579, Japan
Corresponding Author<EMAIL_ADDRESS>Makoto Ohta4,5 1
CREATIS, Université Lyon1, CNRS UMR5220, INSERM U1206, INSA-Lyon, 69621
Villeurbanne, France
2 LIRIS, CNRS UMR 5205, F-69621, France
3 Ecole Centrale de Lyon, France
4 ELyTMaX IRL3757, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard
Lyon 1, Tohoku University, 980-8577, Sendai, Japan
5Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku,
Sendai, Miyagi 980-8577, Japan
6 Graduate School of Biomedical Engineering, Tohoku University, 6-6 Aramaki-
aza-aoba, Aoba-ku, Sendai, Miyagi 980-8579, Japan
Corresponding Author<EMAIL_ADDRESS>Guillaume Lavoué2,3
1 CREATIS, Université Lyon1, CNRS UMR5220, INSERM U1206, INSA-Lyon, 69621
Villeurbanne, France
2 LIRIS, CNRS UMR 5205, F-69621, France
3 Ecole Centrale de Lyon, France
4 ELyTMaX IRL3757, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard
Lyon 1, Tohoku University, 980-8577, Sendai, Japan
5Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku,
Sendai, Miyagi 980-8577, Japan
6 Graduate School of Biomedical Engineering, Tohoku University, 6-6 Aramaki-
aza-aoba, Aoba-ku, Sendai, Miyagi 980-8579, Japan
Corresponding Author<EMAIL_ADDRESS>
###### Abstract

Computational fluid dynamics (CFD) simulation provides valuable information on
blood flow from the vascular geometry. However, it requires to extract
accurate models of arteries from low resolution medical images, which remains
challenging. Centerline-based representation is widely used to model large
vascular networks with small vessels, as it enables manual editing and encodes
the topological information. In this work, we propose an automatic method to
generate an hexahedral mesh suitable for CFD directly from centerlines. The
proposed method is an improvement of the state-of-the-art in terms of
robustness, mesh quality and reproductibility.

Both the modeling and meshing tasks are addressed. A new vessel model based on
penalized splines is proposed to overcome the limitations inherent to the
centerline representation, such as noise and sparsity. Bifurcations are
reconstructed using a physiologically accurate parametric model that we
extended to planar n-furcations. Finally, a volume mesh with structured,
hexahedral and flow oriented cells is produced from the proposed vascular
network model.

The proposed method offers a better robustness and mesh quality than the
state-of-the-art methods. As it combines both modeling and meshing techniques,
it can be applied to edit the geometry and topology of vascular models
effortlessly to study the impact on hemodynamics. We demonstrate the
efficiency of our method by entirely meshing a dataset of 60 cerebral vascular
networks. 92% of the vessels and 83% of the bifurcations where mesh without
defects needing manual intervention, despite the challenging aspect of the
input data. The source code will be released publicly.

Keywords: Arterial network, Centerlines, Hexahedral mesh, Computational fluid
dynamics

##  1 Introduction

Vascular diseases, such as stroke, can cause severe disability or death [30].
The relationship between the topology and geometry of the vascular network and
the onset and the outcome of the pathology is increasingly investigated in the
literature. Computational fluid dynamics is a key tool for this type of study,
at it provides information on the hemodynamics from the vessel geometry [33,
36]. Numerical simulation requires a smooth, anatomically accurate model of
the arterial wall to give reliable results. In the finite element method, the
shape of the cells inside the volume also affects the simulations. Due to
their ability to mesh automatically very complex shapes, tetrahedral cells are
widely used. However, unstructured tetrahedral cells lead to higher
computational cost and less accurate results than flow-oriented, structured
hexahedral cells, as demonstrated by [39], [10] and [13]. Two main approaches
used to address the meshing of vascular networks are segmentation-based and
centerline-based methods.

The segmentation of magnetic resonance angiography (MRA) images is a non-
invasive way to access patient-specific vasculature. A lot of effort was put
in order to develop efficient vessel enhancing filters and to improve the
segmentation methods [22, 27]. However, it remains a challenging task when it
comes to large networks with voxel-sized vessels [25]. Besides, a good
segmentation does not guarantee the accuracy of the network topology (e.g
vessels merging due to the image resolution, disconnected vessels), and
geometry (e.g bumps caused by voxelization). The segmentation result is
generally meshed with tetrahedral elements and hexahedral remeshing is not
straightforward.

Following the tubularity assumption, vessels can be reduced to a centerline-
radius description. Segmentation-based and centerline-based models complement
each others, centerline extraction being used as a pre-processing or post-
processing of segmentation. Many methods to extract vessel skeletons from
binary or raw images were proposed in the literature [34, 17]. As opposed to
image segmentation, centerline-based representation advantageously
incorporates the network topology and enables manual extraction and editing.
This simplified representation is more suitable to the construction of big
databases of large vascular networks [40] or the creation of ideal models. As
it encodes the vessel orientation, it has a high potential for the creation of
meshes with flow-oriented, hexahedral cells. Nevertheless, the representation
of vessels by centerlines lowers the geometrical information content ;
depending on the extraction method, only a limited number of data points are
used and noise can be introduced in the dataset. It causes inaccuracy in the
shape of the vessels and in the position and geometry of bifurcations. These
limitations make it difficult to reconstruct a smooth and physiologically
accurate surface model that matches the requirements of numerical simulation.

In this work, we propose a method overcoming the limitations of centerlines to
generate a high quality mesh ready for CFD. The input data is presented in
section 3. Section 2 reviews the state-of-the-art methods in hexahedral
meshing of vessels and surface reconstruction from centerlines. In section 4,
we introduce a mathematical model to represents vessels as tubes, based on the
approximation of centerline data by penalized splines. The proposed
approximation framework allows for the combination of the spatial coordinates
and radius in a single function. It is shown to be robust to noise and low
sampling of data point, compared to other conventional approximation methods
[28, 11]. The vessel model is combined with a physiological parametric model
of bifurcation proposed by Zakaria et al. [42] to form a what we call a
parametric arterial network. In the section 5, we present a method to create a
structured hexahedral surface and volume mesh with flow oriented cells from
the arterial network model detailed in the previous section. Original
relaxation and smoothing steps are introduced to improve the quality of the
cells without deforming the model shape. The meshing is fast and fully
automatic, which opens the way to numerical simulation in large networks. As
the modeling and meshing steps are intertwined, the proposed framework enables
to easily and inexpensively edit the topology of the networks to study its
influence on blood flow. In section 7, we use our method to mesh a database of
60 large cerebral networks. Other applications such as pathology modeling and
hexahedral remeshing are also presented.

##  2 Related work

###  2.1 Hexahedral meshing

Hexahedral meshing has been investigated over the past years as it offers many
advantages. It better capture the underlying topology of the object with a
lower cell-density, and bridges the gap between representation and physical
simulation, as it simplifies the boundary layer creation and domain
discretization, and provides a basis for NURBS modeling [43]. Hexahedral
meshing is however limited by a far more complex generation process than
standard tetrahedral meshes. In the application to arterial network, the main
challenge is the treatment of bifurcations. In the literature, this task was
addressed by a two-step pipeline ; the bifurcations are first decomposed into
three branches, then the hexahedral mesh is generated using the decomposition.

A variety of method were proposed to obtain a robust branch decomposition. De
Santis et al. introduced semi-automatic methods, ranging from manual selection
of the most relevant slices of the input surface mesh [10], user defined
bifurcation coordinate system [8], to the generation and adjustment of a
block-structure representation of the network [9]. Automatic methods are based
on Voronoi diagram [4], resolution of the Laplace’s equation [37], random-walk
algorithm [41] or branching templates [43] or parametric models [13]. The
hexahedral meshing can then be created from the decomposition through various
techniques ; Copper scheme in the work of [4], template grid sweeping for
[37], [43] and [13], Bezier spline modeling followed by an iso-parametric
transformation of a template mesh [8], projection and refinement of block-
structures [9], Laplacian-based harmonic functions combined with Catmull-Clark
subdivision [41].

There are limitations to the application of the described methods to our
purpose. First, they often rely on manual intervention, which is limiting when
applied to large datasets of complex vascular networks. The automatic methods
proposed involve complex algorithms, and only [8] provided their code through
the user-friendly interface PyFormex, which enables to generate hexahedral
meshes semi-automatically from a single bifurcation vascular geometry.
Moreover, only [13] use centerlines as input, and they are low noise
centerlines extracted from a surface mesh using VMTK. The other methods
require a tetrahedral surface mesh as input for the branch decomposition and
the meshing steps. They can not be applied directly from realistic centerlines
(both sparse and noisy), as those tackled by this work. The application of
such hexahedral meshing methods to centerlines would require a continuous and
accurate reconstruction of the vascular surface to be produced first.

###  2.2 Vascular modeling and reconstruction from centerlines

In this section, we review the methods used to recreate an accurate vascular
surface from centerline information. The main issues to overcome in this task
arise from the defects commonly observed in the vascular centerline extracted
from medical images ; local discontinuities causing a lack of information -
especially at the bifurcation parts -, noise due to the voxelization. In this
context, the smoothness of the vessel surface and the accurate reconstruction
of the bifurcation part are important locks. The reconstruction methods can be
divided in explicit methods, where a tetrahedral mesh of the surface is
produced, and implicit methods where the surface is represented by implicit
functions. Implicit methods employs radial basis functions [18], implicit
extrusion surfaces [19] or local implicit modeling [23] to reconstruct
vascular networks from medical images. If they stand out by their ability to
reconstruct complex branching topology, they are less suited for hexahedral
meshing than explicit methods.

In explicit methods, the vessel surface is obtained by sweeping along the
centerline. The quality of the reconstruction depends on the way centerline
points are approximated, usually with Bezier or spline functions. Details of
the approximation methods, in particular the accuracy of the tangent
estimation and the robustness to noise were not given in previous studies.
Various bifurcation models were proposed. In the work of [24] and [12], the
three branches of the bifurcation are modeled separately and joined at the
bifurcation center. The junction is then blended to restore the continuity, by
a subdivision scheme for [24] and Bezier segments for [12]. This geometric
model is particularly suited for hexahedral meshing, as the junction can be
used as branch decomposition. However, it depends on the accuracy of the
position of the bifurcation center and the tangent of the branches, which is
hard to estimate correctly from centerlines.

Besides, [16] and [14] modeled bifurcations using three tubes connecting the
inlet and outlet sections. A half of each section of the tube is meshed by
sweeping and the hole left in the middle is filled afterwards. Their method
guarantees the smoothness of the model, but results in unnatural looking
bifurcations. Finally, [42] proposed a physiologic model were the bifurcations
are represented by two merged tubes. It was validated with regards to both the
accuracy of the anatomy and of the CFD simulations. It does not rely on the
geometrical center of the bifurcation, but on a set of physiological
parameters (apex, apical sections, inlet and oulet sections). The hexahedral
meshing is less straightforward from this model than the geometric model
presented in the previous paragraph but it shows more anatomical realism. In
[42], the authors extracted the model parameters from a surface mesh, they did
not suggest a way to extract them from centerlines.

In this work, our purpose is to meet the challenges arising from this state-
of-the-art with a framework integrating modeling and hexahedral meshing. Our
method does not require to generate a tetrahedral mesh at any stage of its
application. It produces an hexahedral mesh from centerlines only, with a
controlled trade-off between smoothness and proximity to the input centerline
data points, and a realistic bifurcation shape with customizable physiological
parameters as proposed by [42].

##  3 Input data

The input vessel centerlines we consider are composed of set of data points
with three spatial coordinates (x,y,z), radius value (r), and the connectivity
between points. Data points might have several successors (e.g bifurcations).
A point with $n$ successors is a $n$-furcation. The centerlines are stored
using the swc format or VMTK format of [21]. In this work, we use centerlines
from two publicly available datasets. The Aneurisk database [3] provides 3D
models of the main arteries of the circle of Willis for patients with
aneurism. High resolution centerlines were extracted from the surface meshes
using VMTK software. The BraVa database [40] gathers the centerlines of the
whole cerebral network for 60 patients. To create this dataset, the data
points were manually placed by medical doctors on medical images using imageJ
plugin Neurite Tracer [26] and the radius was automatically computed. As a
result, the data points have a lower spatial resolution and are prone to
errors and noise.

##  4 Modeling

###  4.1 Vessels

In this part, we focus on the modeling of vessels from centerline data; the
case of bifurcations is addressed in the next section. Different models of
centerlines were proposed in the literature, based on the approximation of
data points by Bezier segments [13], regression splines [24], free knot
regression splines or local polynomial smoothing [32]. Only [32] gives the
detail of the implementation of the approximation method and provides a
thorough study of the accuracy of their model in regard with the spatial
coordinates and the derivatives. Actually, the accuracy of both the first and
second derivatives is crucial because the vessel curvature impacts the
haemodynamics [31]. Moreover, meshing techniques are often based on the
normals of the centerline [24, 13]. It is important to note that the proposed
approximation methods [32, 24, 13] focus on the spatial coordinates of the
centerlines, excluding the radius. In this work, we propose a parametric model
of vessels based on approximation by penalized splines. Our approximation
method enables to combine the spatial coordinates and radius in a single
function with physiologically accurate values and derivatives and is robust to
noise and low sampling of the input data.

####  4.1.1 Penalized splines

We want to approximate a set of $m$ points $\\{D_{0},D_{1}...,D_{m-1}\\}$ with
$4$ coordinates $(x,y,z,r)$, using a spline function $s$ defined as

$s(u)=\sum_{i=0}^{n-1}N_{i,p}(u)P_{i},$ (1)

for $u\in[0,1]$, where $N_{i,p}$ is the $i$th basis spline function of order
$p$ and $\\{P_{0},P_{1}...P{n-1}\\}$ the $n$ control points of the spline. The
shape of the basis splines functions and therefore the part of the spline
controlled by a given control point is given by a set of knots.

The main challenge in the approximation of noisy data is to find the optimal
balance between the proximity of the curve to data points and the smoothness
of the curve (i.e. the accuracy of the derivatives). There are two main
approaches to control the smoothness of a spline function. The first is to
change the number of control points : a low number of control points will
result in a smoother curve. In this case, the position of the knots can be
optimized like in [32]. In the other approach, a relatively large number of
control point and a uniform knot vector are used and the smoothness is
constrained by a penalty on the second derivatives [7, 11]. For reasons
further detailed in the next paragraphs, the second approach was judged more
suitable for our task. The vessels are modeled with penalized splines, as
introduced by [11]. For penalized splines, the optimization of the control
points is based on a cost function with two terms. The first term takes into
account the closeness to the data point and the second term the smoothness of
the approximation spline. The parameter $\lambda$ controls the balance between
both closeness and smoothness. The cost function is defined as

$f(P_{0},...,P_{n-1})=\sum_{k=0}^{m}{|D_{k}-s(t_{k})|^{2}}+\lambda\sum_{j=2}^{n}{(P_{j}-2P_{j-1}+P_{j-2})^{2}},$
(2)

where $t$ is a time parametrization vector which associates each data point to
a position on the spline.

####  4.1.2 Approximation strategy

Centerline data provides both the spatial coordinate (x, y, z) and radius r,
two variables of different scale that might show different noise levels. For
this reason, they must be approximated separately. The choice of penalized
splines allows us to dissociate $\lambda$ values for the position and the
radius in a two-step approximation algorithm. With this approach, the spatial
and radius coordinates can be modeled by a single spline.

For the approximation, we use a uniform knot vector and a parametrization
obtained by the chord-length method. The number of control points is set so
that the non-penalized approximation curve (i.e produced by solving equation 2
with $\lambda=0$) has a root mean square distance from the original data lower
than a given value, which is set in this work to $10^{-1}$ for spatial
coordinates and $10^{-3}$ for the radius.

We first solve the linear system arising from equation 2 for the spatial
coordinates (x,y,z) of the centerline data points. The system can be written
as

$P_{(x,y,z)}=(N^{T}N+\lambda_{s}\Delta)^{-1}N^{T}D_{(x,y,z)},$ (3)

where $N$ is the matrix of representation of the basis spline functions and
$\Delta$ is the matrix representation of the difference operator which appears
in the second term of the cost function 2. The optimal value for $\lambda_{s}$
is obtained by minimizing the Akaike criterion $AIC_{2}$. A comparison study
with other selection criteria for $\lambda$, such as the Bayesian information
criterion or cross validation, detailed in section 6, led to the choice of
$AIC$. Then, the linear system is solved for the data $(t,r)$ were $t$ is the
time parametrization of each data points and $r$ their radius value :

$P_{(t,r)}=(N^{T}N+\lambda_{r}\Delta)^{-1}N^{T}D_{(t,r)}.$ (4)

The value of $\lambda_{r}$ is also selected by minimizing AIC on the
time/radius data. The spatial coordinates and radius of the optimized control
points are then concatenated to form the 4-coordinates control points of the
final spline. Figure 1 illustrates this two-part approximation scheme. The
proposed approximation method is compared with other conventional
approximation methods regarding the robustness to noise and low sampling of
the data points in section 6.1.

Figure 1: Approximation of a noisy centerline with the proposed method. (a)
and (b) show respectively the approximation of the spatial coordinates and the
radius. Figure (c) shows the input centerline data in 3-dimensions as red dots
with radius represented by black circles. In (d), the surface defined by the
approximating spline-based model is represented.

###  4.2 Bifurcations

In this part, we focus on modeling the bifurcations from centerline data.

####  4.2.1 Zakaria’s model

[42] proposed a parametric model for non-planar bifurcations. Their model was
validated in regard with both the anatomy and numerical simulation of blood
flow and showed a good agreement with real cerebral bifurcations. It requires
only a few physiological parameters and is well suited for reconstruction of
bifurcations from sparse data. In this model, bifurcations are created by
merging two tubes which represents the daughter vessels. The tubes are defined
by a shared inlet cross section $C_{0}$, separate apical cross sections
$AC_{1}$, $AC_{2}$ and outlet sections $C_{1}$ and $C_{2}$. The apical cross
sections $AC_{1,2}$ are located at the apex point $AP$ of the bifurcation,
where both tubes merges. The outlet sections $C_{1,2}$ are cut one diameter
away from the apex. In total, five cross sections and their normals are
required to build the model. A cross section $C$ is considered circular and is
represented by the three spatial coordinates of its center $P_{c}$, a radius
$r_{c}$ and a normal $\vec{n_{c}}$. The centerline of each tube is defined by
a centerline $spl_{1,2}$. The first segment of the centerline connects the
inlet section $C_{0}$ to the apical section, the second connects the apical
section to the outlet sections. The tangent of the centerline segments matches
the normal of the joined cross sections. The radius along the segments evolves
linearly between $r_{C_{0}}$, $r_{AC_{1,2}}$ and $r_{AC_{1,2}}$,
$r_{C_{1,2}}$. The bifurcation model is illustrated in Figure 2. The
unphysiological sharp angle produced between tubes at the apex is rounded by a
segment of constant radius of curvature $R$.

Figure 2: Five cross sections bifurcation model introduced by [42].

####  4.2.2 Parameter estimation

We introduce an algorithm to estimate the parameters of the bifurcations
directly from the input centerline data. The inlet data points (in light blue
in Figure 3 (a)) are merged with each of the oulet data points (resp. in deep
blue and green in Figure 3 (a)) to form two input centerlines going through
the bifurcation, as shown in Figure 3 (b). The two vessels based on these
centerlines are modeled independently by splines using the approximation
strategy presented in section 4.1. The apex $AP$ of the bifurcation is set as
the point where the surface of the two vessel models first intersect (red dot
on the Figure 3 (c)). $AP$ is then projected on the model splines $spl_{1}$
and $spl_{2}$. The tangent and position of the obtained projection points then
define the normal and the center of the apical cross sections $AC_{1}$ and
$AC{2}$. The outlet sections $C_{1}$ (resp. $C_{2}$) are computed in the same
way from the evaluation of the spline where the length from the apex
projection point is twice the radius of the apical section.

Figure 3: Pipeline of the bifurcation parameter estimation. (a) shows the
inlet and outlet data points, (b) the independent vessel models and (c) the
parameter extraction and resulting bifurcation.

####  4.2.3 Tangent continuity

The full vascular network model is created by assembling the vessels and
bifurcations. In order to preserve the continuity of the different parts of
the network, the end tangents of the extremities of the vessel models must
match the normal of the inlet and outlet cross sections of the bifurcation
models. For this, we introduce an additional constraint on the end point and
tangent in the resolution of the approximation equations 3 and 4 used to model
the vessels. A least square spline approximation with arbitrary end
derivatives was proposed by [29]. We extend this approach to penalized splines
and we propose a weaker constraint which fixes the end tangent while the
derivative is free. We consider a spline $s$ as defined by Equation 1, and
because we work with clamped curves, $s(0)=P_{0}$ and $s(1)=P_{n-1}$.
Moreover, $s^{\prime}(0)$ (respectively $s^{\prime}(1)$) is in the same
direction as vector $P_{1}$ \- $P_{0}$ (respectively $P_{n-2}$ \- $P_{n-1}$).
If we note $S_{0}$ and $S_{n-1}$ the fixed end points and $T_{0}$ and
$T_{n-1}$ the fixed end tangents, the following new conditions are applied to
the system 3:

$\begin{cases}P_{0}=S_{0}\\\ P_{n-1}=S_{n-1}\\\ P_{1}=P_{0}+\alpha T_{0}\\\
P_{n-1}=P_{n-2}+\beta T_{n-1},\\\ \end{cases}$ (5)

where $\alpha$ and $\beta$, the end tangent magnitude, are additional
parameters to optimize. Those constraints guarantees the $G_{1}$ continuity of
the final network model.

##  5 Structured hexahedral meshing

###  5.1 Bifurcations

####  5.1.1 Decomposition

For further meshing, a decomposition scheme is needed to split the bifurcation
in three geometrical branches; one inlet branch and two oulet branches. The
Figure 4 (b) gives an example of branch splitting using three separation
planes. Antiga et al. [5] proposed a bifurcation decomposition scheme based on
Voronoi diagram of the surface mesh. The decomposition they propose is robust
to variations in input geometry and has been successfully used for hexahedral
meshing [5]. However, it was not originally designed to obtain high quality
meshes but to offer a robust mapping of bifurcations. Moreover, it requires a
surface mesh to be computed and its transposition to centerline data is not
straightforward. In this work, we introduced a decomposition scheme which
relies on the spline and bifurcation models described above.

Three separation planes are defined by a set of five points; the apex point
$AP$, which is a already a parameter of the model, two center points $CT_{0}$
and $CT_{1}$ and two separation points $SP_{1}$ and $SP_{2}$. As illustrated
in Figure 4 (a), we first define the geometric center of the bifurcation $X$,
as the barycenter of $AP$, $p_{m_{1}}$ and $p_{m_{2}}$, where $p_{m_{1,2}}$
are the projection of the key points $m_{1,2}$ located at the intersection of
one centerline with the surface of the other vessel. The separation points
$SP_{1,2}$ are obtained by projecting $X$ on the surface in the opposite
direction from $AP$.

Figure 4: Geometric decomposition of the bifurcation model. In (b), the end
cross sections are represented in red and the separations planes in back.

Finally, the position of center points $CT_{0}$ and $CT_{1}$ is obtained by
projecting $X$ on the surface of the vessels. The direction of projection is
the normal to the plane defined by the three points $AP$, $SP_{1}$ and
$SP_{2}$. The separation points $AP$, $SP_{1}$ and $SP_{2}$ are finally
connected to the center points $CT_{0}$, $CT_{1}$ by arcs, which delineate a
geometrical frontier between the branches of the bifurcation (see Figure 4
(b)), providing the desired branch decomposition. This decomposition enables
to handle large radius difference between the daughter vessels. The barycenter
$X$ is naturally closer to the vessel with the smallest radius, which relaxes
the angles between the separation planes and improves the quality of the
resulting mesh.

####  5.1.2 Surface meshing

In this step, we create the mesh grid connecting the end cross sections to the
separation planes with a set of successive sections, as illustrated in Figure
6. Each section of the mesh have a number $N$ of nodes where $N$ can be any
multiple of $4$. Figure 5 illustrates the process for constructing them. We
first compute the $N$ nodes of the end sections $C_{0}$, $C_{1}$ and $C_{2}$.
A normalized reference vector $\vec{ref_{C}}$ which minimizes the rotation
with the separation points $SP_{1,2}$ is defined for each end cross section
$C$. The nodes of the end sections are placed on the outline of the cross
section, with evenly spaced angles starting by $\vec{ref_{C}}$ and rotated
counterclockwise. The nodes of the separation planes are positioned with
equally sampled angles along the arcs which join the separation point $AP$,
$SP_{1}$ and $SP_{2}$ to both center points $CT_{0}$ and $CT_{1}$.

Figure 5: Computation of the nodes (black dots) of the end cross sections and
the separation planes for $N=8$ and splines $tspl_{1,2,3}$.

Once the nodes of the end section and the separation half-sections are
computed, they are connected to form a surface mesh. The nodes of $C0$ are
connected to the nodes of the half-sections defined by $SP_{1}$ and $SP_{2}$,
and the nodes of $C_{1,2}$ are connected respectively to the nodes of the
half-sections $SP_{1,2}$ and $AP$.

We first define an initialization of the 3D trajectory that connects two
nodes, as shown in the left column of Figure 6. This initialization is an
approximation which is used to control the topology and geometry of the final
mesh grid, but it does not necessary lies on the exact surface of the
bifurcation at this point. The initial trajectories are evenly sampled with
$n$ nodes, where $n$ determines the number of cross sections to compute along
a given branch. This number is proportional to the radius of the end section
of the branch, by a coefficient $d$ which can be adjusted to obtain the
intended density of faces in the mesh.

The nodes are then projected radially to the surface of the two vessels, as
illustrated in the right column of Figure 6. The direction of the projection
is important to maintain the quality of the faces of the initial grid after
projection. Ideally, the nodes of the initial trajectory must be displaced
only radially from the center of the branch vessel. However, the shape splines
$spl_{1}$ and $spl_{2}$ does not constitute a good approximation of the
centerline of the branches joining the geometric center of the bifurcation.
For this reason, we create another set of splines $tspl_{1,2,3}$ connecting
the center of each end section to the center $X$ of the bifurcation,
represented in blue in Figure 5. The nodes are projected to the surface of the
bifurcation model according to the normal of this new set of splines.

Figure 6: Initial surface mesh and mesh after projection for the two types of
initialization considered. The red squares emphasize the impact of the two
types of initialization on the final mesh.

The properties of the resulting mesh depends on the initial trajectory
approximation. Figure 6 illustrates the meshes obtained after projection
considering two types of initialization. The first row shows the simple case
where the nodes of the end sections are linearly connected to the nodes of the
separation geometry. In the second row, connection trajectories are computed
so that the normal of the end sections is preserved in the output surface
mesh. Both approximations allow to preserve the topology of the grid and the
quality of the faces after projection. The initial trajectories with normal
preservation are closer to the actual surface of the bifurcation, there is
less displacement of the nodes during projection. Moreover, the preservation
of the normal of the end sections enables to include the bifurcation mesh in
larger arterial network; the connecting curves can be smoothly extended to
downstream vessels. In the rest of this work, we use the normal preserving
initialization.

####  5.1.3 Relaxation

The projection step of the meshing method results in an uneven sampling of the
nodes along the trajectories and can lead to faces with heterogeneous size or
important skewness. Moreover, a rupture of continuity is observed when the
curves crosses the separation between two branches. Those unwanted features
are corrected by relaxation of the nodes of the surface mesh. Mesh smoothing
methods are an easy way to reduce the skewness of faces, however it triggers
important deformations of the general shape of the model. To avoid
deformations, [38] proposed to combine smoothing with a back projection on the
surface. Following this approach, an iteration of Laplacian smooth (relaxation
factor of $0.8$) is first applied to the bifurcation mesh, then the nodes are
projected back to the original surface. To prevent cross sections from
intersecting, the projection is made in the direction of the line connecting
the center of the cross section to the node to project. This process can be
repeated until the relaxation is satisfying. Figure 7 displays a bifurcation
mesh after 1 and 5 relaxation iterations. The faces are colored according to
their geometric quality, measure by the scaled Jacobian. We observe that while
the shape of the model is preserved, the quality of the faces near the
separation plan is improved, and the trajectories smoothly crosses the
separation plans. Based on the average quality of the faces, we estimate that
5 relaxation iterations are sufficient.

Figure 7: Original bifurcation mesh and mesh after 1 and 5 relaxation
iterations. The scaled Jacobian is used to measure the quality of the cells
between $-1$ (poor quality) and $1$ (high quality).

####  5.1.4 Apex smoothing

The last step of the bifurcation meshing is the smoothing of the apical
region. The model presents an unwanted sharp angle were the two vessels merge.
The curvature in the apex regions impacts the pressure and velocity fields
obtained by numerical simulation, as shown by [15]. Conventional mesh
smoothing methods (e.g Laplacian, Taubin smoothing) are fast and can produce
smooth meshes with high quality faces. However as they are global methods,
they struggle to generate important local deformations. [42] proposed to
smooth the apex region by projecting the nodes on a sphere of given radius,
rolling on the surface. This method is accurate, but it is computationally
expensive and might not preserve the quality of the cells in the case of
hexahedral meshes. Taking advantage of the topology of the proposed surface
mesh, we reduced this complex 3-dimensional problem to a 2-dimensional
problem.

Figure 8: Illustration of the apex smoothing pipeline.

Figure 8 illustrates the proposed smoothing method. The 3D polylines
connecting two nodes of the end sections of the bifurcation are extracted
(e.g, the curve in red on Figure 8). They are then projected on the 2D plane
defined by the normal of the mesh at the separation point and the normal of
the separation plane (resp. green and blue arrows on Figure 8). A circle whose
radius correspond to the desired apex radius of curvature is rolled along the
2D curves. The position of the circle where it is in contact with a further
part of the curve is mathematically computed. The points located under the
circle are moved to its outline, while preserving their original sampling.
Finally, the new coordinates of the points are projected back on the original
3D referential to form the output surface mesh.

Figure 9: Apex smoothing with different radius of curvature values. The
colormap encodes the local distance to the original mesh, on the left.

The described smoothing method enables to control the direction of projection
and the sampling of the projected nodes. Therefore, the quality of the faces
is preserved. As shown in Figure 9, the smoothing is very local and does not
affect the shape of the vessels outside of the apical region.

####  5.1.5 Planar n-furcations

If the cerebral arterial network is composed of a majority of bifurcations,
multifurcations may also be present (e.g trifurcations are frequently found on
the basilar artery). To address this requirement, we generalized the model of
[42] to planar n-furcations. The generalized n-furcation model is built with
$n-1$ splines, $2n+1$ cross sections and $n-1$ apex points, as illustrated for
the case $n=3$ in Figure 10 (a). We adapted the decomposition scheme presented
in section 5.1 to compute $n+1$ separation plans, as in Figure 10 (b). Figure
10 (c) shows an example of planar trifurcation mesh obtained with this
generalization.

Figure 10: (a) and (b) respectively illustrates the parametric model and the
branch decomposition scheme for a trifurcation. (c) shows an example of
trifurcation mesh.

###  5.2 Vessels

For the vessels, we adapted the meshing method proposed by [13] to the
parametric model proposed in section 4.1 in order to obtain surface and volume
meshes.

####  5.2.1 Surface meshing

To mesh the surface of a vessel, the spline model is evaluated at a set of
time values equally sampled in the $[0,1]$ interval. Those values are used to
set the center position and radius of the cross sections along the vessel (i.e
the longitudinal resolution of the mesh). The density of cross sections
(number of sections per $mm$) is proportional to the mean radius of the
vessel, with a proportional coefficient $d$ which is set by the user. From
each center positions, $N$ nodes are radially projected on the model surface
to form a circular cross section. The projection vector is swept along the
centerline and avoids twisting between the sections. The successive sections
are connected to each other to form the mesh faces. In the case of vessels
connecting one bifurcation to another, an extra rotation is smoothly applied
to the cross sections along the vessels so that the last vessel section is
aligned with the first section of the next bifurcation.

####  5.2.2 Volume meshing

The volume of the vessel is meshed following the method of [13]. From each
cross section of the surface mesh, a structured O-grid pattern is created. It
has 3 different areas, including the boundary layers, intermediary layers and
a central block. The relative size $\alpha$, $\beta$, $\gamma$ of the areas,
the number $N_{\alpha}$ of boundary layers and the number $N_{\beta}$ of
intermediary layers can be adjusted. The separation planes of the n-furcations
are handles by combining $n+1$ halves grids. The successive O-grid patterns
are connected to form the hexahedral cells of the volume mesh, as showed in
Figure 11.

Figure 11: Illustration of the O-grid pattern and volume meshing method.

##  6 Results

In this section, we evaluate both the modeling and the meshing methods
proposed. The robustness and accuracy of the proposed vessel modeling method
is assessed in an comparative study performed on a synthetic dataset of
distorded centerlines. The performance of the meshing method is measured in
terms of cell quality and computational time.

###  6.1 Vessel model evaluation

In this part, we evaluate the robustness of the approximation method presented
in section 5.2 to noise and low sampling of the data points.

####  6.1.1 Validation dataset

For this evaluation, we built a dataset of ground truth vessel models. Four
surface meshes of cerebral arteries from the Aneurisk database were selected.
For each mesh, a single vessel starting from the inlet of the network and
ending at an outlet was selected so that is does not include pathologies but
go through bifurcations, where we generally observe high curvature and big
radius change. The selected vessel centerlines were extracted with a good
resolution and low noise using the VMTK software. The obtained high quality
centerline data points were then approximated by a 4-coordinate spline $s$
that constitute the ground truth. The control points were manually added and
the accuracy of the fitting of the spatial coordinates, radius and first
derivatives were checked visually until the approximation was judged
satisfying. Ground truth vessels and their creation process are illustrated in
supplementary material, section 1.1.

To evaluate the robustness of our approach, those ground truth data were
distorded to mimic defects commonly observed in realistic centerline data; low
sampling and noise. Spatial noise and radius noise were applied separately, as
they might differ in level. To generate spatial noise, the data points were
displaced from their original position. The magnitude of displacement is
randomly picked from a zero-centered Gaussian distribution with standard
deviation $\sigma_{spatial}$. In order not to affect the radius values, the
direction of the displacement is normal to the ground-truth spline $s$. Random
radius noise is generated randomly from a zero-centered Gaussian distribution
of standard deviation $\sigma_{radius}$ and added to the ground truth radius.
In both case, the applied standard deviation value is proportional to the
point radius, as indicated in Table 1, in order to keep a similar levels of
noise between big and small vessels. Finally, low sampling is obtained by
removing data points along the centerline to reach target point densities.

Table 1: Parameters used for the distortion of the ground-truth centerlines density $(mm^{-1})$ | $2$ | $4$ | $10$ | $16$ | $20$
---|---|---|---|---|---
$\sigma_{radius}(mm)$ | $0.01r$ | $0.05r$ | $0.1r$ | $0.3r$ | $0.5r$
$\sigma_{spatial}(mm)$ | $0.01r$ | $0.05r$ | $0.1r$ | $0.3r$ | $0.5r$

For each density value in Table 1, ten combination of noise parameters are
used, spatial and radius noise being added separately. Radius noise is applied
to the ground-truth data with parameters $\sigma_{radius}$ as given in Table 1
while the spatial noise is set to $0$. Then spatial noise is applied with
parameters $\sigma_{spatial}$ as given in Table 1 while the radius noise is
nul. Each noise combination is repeated three times to account for stochastic
effect; we get $30$ data per density value, thus $150$ in total. This is done
for the four vessels of the ground-truth dataset, bringing the number of data
in the distorded dataset to $600$.

####  6.1.2 Approximation methods

To demonstrate the robustness and the accuracy of the proposed approximation
strategy presented in section 4.1, we compared it to approximation methods
commonly used in the literature. If they claim to use splines or Bezier
functions, few authors provided details on the approximation method employed.
For this reason, we implemented four commonly used approximations methods with
incremental complexity - including the original penalized splines proposed by
[11] \- in order to emphasize the contributions of the proposed method.

  * •

Global Non-Penalized (GNP) : In this basic approach, the control points are
optimized without smoothness penalty in the cost function (Equation 2 with
$\lambda=0$). The number of control points is set to match the RMSE threshold
given in section 4.1, paragraph 2. We call it global because spatial and
radius dimensions are not addressed separately.

  * •

Global Non-Penalized with Akaike criterion (GNP-AIC) : Optimizing the number
of control points to obtain the desired spline smoothness is a common
approximation method in the literature. In this approach, the optimal number
of control points minimizes the Akaike information criterion ([1]) $AIC_{1}$ :

$AIC_{1}=m\times log(SSE)+8(n+p).$ (6)

where $m$ is the number of data points, $p$ the degree of the spline, $n$ the
number of control points and SSE is the sum squared error from the data
points, including their four coordinates.

  * •

Global Penalized with Akaike criterion (GP-AIC) : This approach corresponds to
the original approximation by penalized splines described in [11]. It uses the
same global approach as in GNP, but with a smoothing penalty defined with a
parameter $\lambda\neq 0$ as in Equation 2.

  * •

Spatial coordinates and Radius Penalized with Akaike criterion (SRP-AIC) : The
approximation strategy that we propose in this work penalizes spatial and
radius dimensions separately. The comparison of our strategy with GP-AIC
allows to evaluate the contribution of treating the spatial and radius
coordinates individually.

In methods GP-AIC and SRP-AIC, the criterion used to optimize the $\lambda$
values is another formulation of the Akaike information criterion ($AIC_{2}$),
adapted to penalized splines, as proposed by [11]:

$AIC_{2}=m\times log(SSE/m)+2tr,$ (7)

where $tr$ is the trace of the matrix $H=N(N^{t}N+\lambda\Delta)^{-1}N^{t}$.

Unlike $AIC_{1}$, it is not employed to choose an optimal number of control
points but to select the optimal value for the smoohting parametre $\lambda$.
This criterion was compared to other criteria of the literature ; the
corrected Akaike information criterion of [20], the Schwarz’s Bayesian
criterion of [35], the cross validation criterion and generalized cross
validation criterion of [7]. According to this study, given in supplementary
material (section 1.3), the Akaike criterion yielded the best results on our
data.

####  6.1.3 Quality metrics

A total of six quality metrics were selected to evaluate the approximation
strategies presented in the previous paragraph. To build an accurate measure
of distance between the ground-truth spline $s$ and the approximation spline
$\hat{s}$, we project one curve on the other. As illustrated in Figure 12, two
matched sets of time parameters are built. The spline $s$ is equally sampled
with a time vector $t$, then projected on $\hat{s}$ according to the minimum
distance to form the matched time vector $T$.

Figure 12: Matching time parameters by minimum distance projection from $s$
onto $\hat{s}$

Once the projection is performed, the matched values can be compared. We use
the root mean squared error (RMSE) as a measure of the closeness of the
approximation spline to the ground-truth spline. The spatial coordinates and
the radius values are treated separately in the evaluation. We note
$\mathrm{RMSE_{radius}}$ (respectively $\mathrm{RMSE_{spatial}}$) the root
mean squared error of the radius (respectively the spatial coordinates). In
order to have a robust comparison between the curves, the projection is
computed in both ways (from $s$ to $\hat{s}$ and from $\hat{s}$ to $s$) and
the final RMSE value is the average of the RMSE yielded by both projections.

The accuracy of the first derivatives of the model is evaluated by the metrics
$\mathrm{RMSEder_{spatial}}$ and $\mathrm{RMSEder_{radius}}$. As curvature is
commonly considered in hemodynamic studies, the model performance with regard
to the centerline curvature is also measured, by the metric
$\mathrm{RMSEcurv}$. Finally, the length of the vessel affects the delay of
blood arrival between inlet and oulet of the vascular tree in numerical
simulations. Therefore, the difference $\mathrm{L_{diff}}$ of length between
the ground-truth and the approximated centerline was considered.

####  6.1.4 Results

As the spatial and radius distortions are not comparable in nature and
magnitude, the evaluation results are presented in two different tables. Table
2 (respectively Table 3) shows the mean values of the six quality criteria for
the four methods after radius noise (respectively spatial noise) addition. As
expected, the non-penalized model (GNP) is sensible to the added noise and
performs poorly for all radius-related metrics. In Figure 13, the radius
estimation error is clearly visible on the vessel produced by this method. In
the same way, the spatial-related metrics are impacted when spatial noise is
added (Table 3). In addition, a tendency to overfit the data is observed in
Table 2, causing surprisingly high spatial error. The overfitting and noise
problems are partially solved by optimizing the number of control points with
method GNP-AIC. However, this approach still yields a poor approximation of
the derivatives : as the number of control points is lower, the space between
data points might not be correctly interpolated, which particularly impacts
the curvature values.

Figure 13: Mesh resulting from the approximation of distorded data
($density=1mm^{-1}$, $\sigma_{radius}=0.1$) by three of the methods compared
in section 6.1.

The penalized approximations GP-AIC and SRP-AIC enabled to drastically improve
the estimation of the derivatives and curvature. Finally, the advantage of
SRP-AIC over GP-AIC is demonstrated both in the result tables 2 and 3 and in
Figure 13. The global smoothing penalty used in GP-AIC forces a trade-off
between the radius and spatial accuracy. In Figure 13, the radius of the
vessel produced by GP-AIP is very similar to the ground-truth vessel, but in
return the trajectory of the centerline appears too smoothed. On the other
hand, both the radius and trajectory of the vessel produced with SRP-AIC is
closer to the ground-truth. As a conclusion of this comparison, the proposed
approximation method shows a very good robustness to the defects of the input
data while enabling to simultaneously and accurately model the vessel
centerline and radius. More results are provided in supplementary material,
section 1.2.

Table 2: Overall evaluation of the approximation methods : mean values of the quality criteria for all the centerlines distorded by radius noise addition. The cells in gray corresponds to the lowest error for each metric. | GNP | GNP-AIC | GP-AIC | SRP-AIC
---|---|---|---|---
$\mathrm{RMSE_{spatial}}$ | 8.462 | 0.034 | 0.053 | 0.029
$\mathrm{RMSE_{radius}}$ | 17.523 | 0.095 | 0.042 | 0.043
$\mathrm{RMSEder_{spatial}}$ | 0.218 | 0.118 | 0.042 | 0.009
$\mathrm{RMSEder_{radius}}$ | 0.391 | 0.214 | 0.032 | 0.032
$\mathrm{RMSEcurv}$ | 1919.428 | 190.531 | 0.060 | 0.035
$\mathrm{L_{diff}}$ | 718.906 | 0.057 | 0.207 | 0.004
Table 3: Overall evaluation of the approximation methods : mean values of the quality criteria for all the centerlines distorded by spatial noise addition. The cells in gray corresponds to the lowest error for each metric. | GNP | GNP-AIC | GP-AIC | SRP-AIC
---|---|---|---|---
$\mathrm{RMSE_{spatial}}$ | 0.511 | 0.152 | 0.099 | 0.096
$\mathrm{RMSE_{radius}}$ | 0.008 | 0.009 | 0.018 | 0.007
$\mathrm{RMSEder_{spatial}}$ | 0.314 | 0.343 | 0.075 | 0.076
$\mathrm{RMSEder_{radius}}$ | 0.015 | 0.019 | 0.021 | 0.013
$\mathrm{RMSEcurv}$ | 1.524 | 2.362 | 0.085 | 0.091
$\mathrm{L_{diff}}$ | 50.180 | 15.071 | 0.252 | 0.207

###  6.2 Mesh quality

In CFD, the accuracy and stability of the simulation is affected by the
quality of the mesh. In order to evaluate this quality, we computed the scaled
Jacobian of the cells in the meshes generated with the proposed method. The
scaled Jacobian ranges from -1 (worst quality) and 1 (best quality). Negative
values indicate invalid cells. The volume meshes for 60 patients from the
BraVa database were generated (see section 7.3 for details), with the
following parameters ; $N=24$, $d=0.2$, $\alpha=0.2$, $\beta=0.3$ ,
$\gamma=0.5$, $N_{\alpha}=10$, $N_{\beta}=10$. The cells of the bifurcations
and vessels are evaluated separately. Failed bifurcations and vessels (see
section 7.3) were excluded from the study. The histograms of scaled Jacobian
for the 60 patients is given Figure 14, together with an example branch.

Figure 14: Distribution of the scaled Jacobian values of the mesh cells. The
histogram (b) represents the quality of bifurcations cells and the histogram
(c) the quality of vessels cells. The image (a) illustrates the location of
high and low quality cells within a mesh.

As shown in the image (a) of Figure 14, the bifurcations are the most
challenging structures to mesh. The lower quality cells are mainly localized
in the bifurcation separation planes. Nevertheless, we achieved a very good
overall quality for bifurcation cells, with $71\%$ of the cells with scaled
Jacobian value higher than $0.9$. The vessel cells have a even better quality,
with $95.7\%$ of the cells having a scaled Jacobian higher than $0.9$. In
terms of mesh quality, our method improves the state of the art. Indeed, only
$49\%$ of the cells have a scaled Jacobian above $0.9$ in average on the
distributions given for three large cerebral networks in [13]. This proportion
goes up to $62\%$ of the cells of the abdominal aortic artery geometry meshed
by the method of [41]. Finally, in [9], between $65\%$ and $82\%$ \- depending
on the case and the cell density - of the cells of the aortic arch meshed have
a scaled Jacobian value between $0.8$ and $1$. Quantitatively, our method
gives better results, especially given that the histograms for the other
methods were computed on all the cells, both bifurcations and vessels.
However, we bear in mind that the study of [9] and [41] focuses on arterial
geometries that differs from our study.

###  6.3 Computation time

The computational time of the modeling and meshing steps for five patients of
the BraVa database were computed. The results for three of them are given in
Table 4. The average time for modeling a large cerebral vascular network is
about 16 minutes. The time for the volume meshing step is given in Table 4 for
different cell densities. The average meshing time goes from 24.6 minutes for
a the coarse mesh to 49.7 minutes for a fine mesh. We want to stress that this
study was performed on large networks, with a high number of bifurcations
(around $100$) and vessels (around $200$). The meshing time increases with the
number of bifurcations and vessels, while the modeling time is affected by the
number of data points.

Table 4: Computational time required to model and mesh large vascular networks from the BraVa dataset. id |  | furcation
---
(#)
| vessel
---
(#)
| data point
---
(#)
| modeling
---
time (min)
| cells
---
(#)
| meshing
---
time (min)
P1 | 96 | 194 | 2816 | 11.3 |  | 1389k
---
1853k
2316k
2779k
| 20.4
---
25.7
31.2
38.4
P2 | 101 | 203 | 3531 | 18.3 |  | 1916k
---
2555k
3193k
3832k
| 27.5
---
38.2
49.1
67.5
P3 | 107 | 216 | 3474 | 16.8 |  | 1737k
---
2316k
2895k
3474k
| 26.3
---
36.1
44.2
55.9

Besides, a large part of the meshing time corresponds to the computation of
the surface nodes ; in average 17.4 minutes for a coarse mesh and 34.8 minutes
for a fine mesh. The volume mesh is generated directly from the nodes of the
surface mesh without recomputing them. Finally, meshing can be run in
parallel, by splitting the network into parts to be meshed on different CPUs.
Using 12 CPUs, we were able to reduce the meshing computational times given in
the Table 4 by a factor 5.

##  7 Applications

Several applications of our framework are proposed in this section.

###  7.1 Deformation

The proposed model is based on the assumption that vessel cross sections are
circular, which is limiting when dealing with pathological vessels. A way to
address this limitation is to deform the cross sections to match a target
surface as a post processing. If the user input data is a surface mesh, we
propose the following alternative use of our meshing framework:

  1. 1.

Extract the centerline from the surface mesh (using VMTK software for
example),

  2. 2.

Create a tubular mesh from the centerline using the proposed method,

  3. 3.

Deform the tubular mesh to match the original surface.

Figure 15 illustrates an example of this pipeline to mesh arteries with
aneurysms. In the deformation step, the nodes are individually projected onto
the surface of the target mesh. To prevent the sections from intersecting, the
nodes are projected radially from the section center.

Figure 15: Structured hexahedral meshing of cerebral arteries with a fusiform
or saccular aneurysm by deformation. On the left, the tubular mesh, obtained
by our framework, is superimposed on the target surface. On the right, the
mesh after projection.

As shown in Figure 15, saccular aneurysms are initially modeled as bifurcating
vessels and then deformed. Because the shape of the volume mesh pattern
depends on the position of the section nodes (cf Section 4.1), the deformation
of the surface mesh is smoothly conveyed to the cells inside the mesh, as
illustrated by Figure 16.

Figure 16: Cross section pattern before and after deformation. The
corresponding slice of the target surface mesh is represented in grey.

This pipeline is not limited to pathological vessels. It can be extended to
remesh any vessel surface mesh with hexahedral cells.

###  7.2 Topology and geometry editing

The relationship between the vascular tree topology and geometry (e.g the
different configuration of the circle of Willis, vessel angle) and the
hemodynamics have been studied extensively in the literature, using ideal or
patient specific models [6, 2]. In this context, the proposed meshing
framework finds applications in creating and editing vascular models. Because
only a few data points are required for the meshing, the bifurcation angles,
the radius or trajectory of a vessel can be modified effortlessly. Figure 17
provides examples of such modifications. Topological and geometrical
information are linked within our parametric vascular tree representation:
this facilitates the identification and modification of the data points of a
branch of interest. Moreover, the bifurcations are based on a parametric model
whose parameters (cross sections, apex smoothing) can be adjusted. This model
guarantees the physiological realism of the bifurcations even when artificial
data points are used (e.g modeling vessels as straight lines). Thanks to the
modeling of bifurcation as merging vessels, the removal of a bifurcation
branch does not modify the trajectory of the other branch through the
bifurcation part, as illustrated on the right of Figure 17. Branches and
arterial territories inexpensively with a local re-computation of the model
and mesh parts.

Figure 17: Editing of a model of basilar artery using our framework. The
bifurcation angle and the radius of the original vertebral artery are
modified, one of the inlet vertebral arteries is removed. The centerline
points and radius used to produce the meshes are represented on the top left
for each case in which they were modified.

###  7.3 Large cerebral arterial network meshing

Figure 18: Top view of 8 meshes among the 60 generated from the patients of
the BraVa database. The bifurcations were the meshing algorithm has failed
(i.e at least of of the cells has a negative Jacobian) are represented in red.
The cross sections of the vessels with cells of negative Jacobian values are
also represented in red.

To demonstrate further the robustness and fully automatic nature of our
method, we applied it to 60 patients of the BraVa dataset. The meshes produced
for 8 of the patients are given in Figure 18. This dataset is considered as
very challenging for several reasons. As the centerlines were extracted
manually by medical doctors, they are noisy and have a low sampling. The
superimposition of the centerline data points on the magnetic resonance
angiography image in Figure 19 shows the high level of noise encountered in
the input data, both in the radius estimation and the spatial positions.
Besides, by computing the ratio of the number of data points on the total
length of the connecting polyline, we estimated the average point density in
the database to $0.45$ $mm^{-1}$, which is very low.

We evaluate the percentage of successfully meshed vessels and bifurcations
separately, as the meshing method is different. The vessels or bifurcations
with at least one cell with negative scaled Jacobian score are considered
failed. With this strict definition, a total of $83\%$ of the bifurcations and
$92\%$ of the vessels were successfully meshed.

The main reason for the failure of the vessel mesh is a too high curvature -
mainly in the arteries with high tortuosity such as the internal carotid
arteries - caused by a sharp angle in the input centerline. The main causes of
failure for the bifurcations were very low bifurcation angles and misplacement
of bifurcation points in the input data. This last case is illustrated on the
image (b) of Figure 19. We can see that the bifurcation point in the
centerline data was positioned too far downstream in the main vessel, causing
one of the daughter vessel to go backwards from the direction of the flow with
a sharp angle. As we use an oriented bifurcation model, it failed to correctly
represent the geometry. The image (a) of Figure 19, on the other hand,
illustrates a successful reconstruction of the trajectory of the vessel.
Although the input centerline was very imprecise both in the radius estimation
and point positions, we were able to produce a smooth model, closest to the
vessel geometry as given by the medical image. Moreover, as shown in the
insert of Figure 19 (a), even challenging topologies (e.g short connecting
segment between bifurcations) can be successfully meshed with hexahedral
elements. An image of all the meshes of the database, with failure areas
highlighted is given in supplementary materials, section 1.3.

Figure 19: Example of success and failure of our method for one patient of the
BraVa database. The whole brain mesh is represented in the middle with focus
on two parts of the network. For each focus, the original centerline data
points are represented by red dots (center) and black circles (radius). The
mesh obtained is superimposed on the data points, with an highlight on the
relevant parts. On the right image the original centerline data points are
overlayed on the original MRA image.

##  8 Conclusion

In this article, we addressed the problem of the reconstruction and meshing of
large vascular networks from noisy, sparse centerlines. The proposed method is
robust, accurate and automatic. It opens to way to CFD simulations in large
vascular networks manually or semi-automatically extracted by medical doctors,
with a minimal manual intervention. An original approximation method unifying
the spatial and radius information in a single function is proposed to model
the vessels. The use of a bifurcation model based on physiological parameters
is associated with a new hexahedral meshing and smoothing techniques to
produce bifurcations with a realistic shape and high quality cells in a
reasonable time. Our method finds application in automatic meshing of large
databases of vascular centerlines and hexahedral remeshing of non-tubular or
pathologic vessels. It is well suited for the creation of realistic ideal
vascular network models and the study of the impact of topological (branch
removal) and geometrical (branch angle) on blood flow.

We acknowledge some limitations to this work. The pipeline was mainly tested
on cerebral vascular networks so far, and non-planar n-furcations ($n>3$) that
are common in other vessels (e.g aorta) were not addressed. In addition, the
robustness of the modeling and meshing method could be further improved as it
failed on some cases. For this, we would like to integrate more physiological
constraints on the bifurcation and vessel models such as a maximum curvature
or maximum vessel angle for instance. Besides, we want to emphasize that our
objective with this work is not to improve the performance of the segmentation
or centerline extraction algorithms but to acknowledge the limitations of the
realistic data and generate meshes as close as possible to the real anatomy
from flawed centerlines and existing databases. Hence, the accuracy of the
reconstruction depends on the accuracy of the input centerlines and might not
be in total accordance with the real patient anatomy has showed in medical
images. In order to address the issue of the agreement between the medical
image and the centerline points position, we plan to develop a vascular
network editing software, with a relevant user interface. This interface will
integrate the different methods (modeling, meshing, deformation, branch
removal) described in the article and could also integrate editing
functionalities enabling to modify the input centerline data points in
accordance with the medical image in a 3D viewer. It would open vascular
modeling and hexahedral meshing to medical doctors and non-expert users.

## Acknowledgments

We would like to thank Simon Tupin, who initiated this work and Erwan Maury
for his contribution of the development apex smoothing method. Founding is
acknowledged through AURA region (SIMAVC project) and ANR 20-CE45-0011
(PreSpin project).

## References

  * [1] Hirotogu Akaike.  Information theory and an extension of the maximum likelihood principle.  Proceedings of the Second International Symposium on Information Theory, pages 267–281, 1973.
  * [2] Martin Sandve Alnæs, Jørgen Isaksen, Kent-André Mardal, Bertil Romner, Michael K Morgan, and Tor Ingebrigtsen.  Computation of hemodynamics in the circle of willis.  Stroke, 38(9):2500–2505, 2007.
  * [3] Aneurisk-Team.  AneuriskWeb project website, http://ecm2.mathcs.emory.edu/aneuriskweb.  Web Site, 2012.
  * [4] Luca Antiga, Bogdan Ene-Iordache, Lionello Caverni, Gian Paolo Cornalba, and Andrea Remuzzi.  Geometric reconstruction for computational mesh generation of arterial bifurcations from ct angiography.  Computerized Medical Imaging and Graphics, 26(4):227–235, 2002\.
  * [5] Luca Antiga and David A Steinman.  Robust and objective decomposition and mapping of bifurcating vessels.  IEEE transactions on medical imaging, 23(6):704–713, 2004.
  * [6] BMW Cornelissen, JJ Schneiders, ME Sprengers, René van den Berg, P Van Ooij, AJ Nederveen, E Van Bavel, WP Vandertop, CH Slump, HA Marquering, et al.  Aneurysmal parent artery–specific inflow conditions for complete and incomplete circle of willis configurations.  American journal of neuroradiology, 39(5):910–915, 2018.
  * [7] Peter Craven and Grace Wahba.  Smoothing noisy data with spline functions.  Numerische mathematik, 31(4):377–403, 1978.
  * [8] Gianluca De Santis, Matthieu De Beule, Patrick Segers, Pascal Verdonck, and Benedict Verhegghe.  Patient-specific computational haemodynamics: generation of structured and conformal hexahedral meshes from triangulated surfaces of vascular bifurcations.  Computer methods in biomechanics and biomedical engineering, 14(9):797–802, 2011.
  * [9] Gianluca De Santis, Matthieu De Beule, Koen Van Canneyt, Patrick Segers, Pascal Verdonck, and Benedict Verhegghe.  Full-hexahedral structured meshing for image-based computational vascular modeling.  Medical engineering & physics, 33(10):1318–1325, 2011.
  * [10] Gianluca De Santis, Peter Mortier, Matthieu De Beule, Patrick Segers, Pascal Verdonck, and Benedict Verhegghe.  Patient-specific computational fluid dynamics: structured mesh generation from coronary angiography.  Medical & biological engineering & computing, 48(4):371–380, 2010\.
  * [11] Paul HC Eilers and Brian D Marx.  Flexible smoothing with b-splines and penalties.  Statistical science, 11(2):89–121, 1996.
  * [12] Mahsa Ghaffari, Chih-Yang Hsu, and Andreas A Linninger.  Automatic reconstruction and generation of structured hexahedral mesh for non-planar bifurcations in vascular networks.  In Computer Aided Chemical Engineering, volume 37, pages 635–640. Elsevier, 2015.
  * [13] Mahsa Ghaffari, Kevin Tangen, Ali Alaraj, Xinjian Du, Fady T Charbel, and Andreas A Linninger.  Large-scale subject-specific cerebral arterial tree modeling using automated parametric mesh generation for blood flow simulation.  Computers in biology and medicine, 91:353–365, 2017.
  * [14] Jixiang Guo, Shun Li, Yim Pan Chui, Jing Qin, and Pheng Ann Heng.  Mesh quality oriented 3d geometric vascular modeling based on parallel transport frame.  Computers in biology and medicine, 43(7):879–888, 2013.
  * [15] IV Haljasmaa, AM Robertson, and GP Galdi.  On the effect of apex geometry on wall shear stress and pressure in two-dimensional models of arterial bifurcations.  Mathematical Models and Methods in Applied Sciences, 11(03):499–520, 2001.
  * [16] Xiaoxiao Han, Richard Bibb, and Russell Harris.  Design of bifurcation junctions in artificial vascular vessels additively manufactured for skin tissue engineering.  Journal of Visual Languages & Computing, 28:238–249, 2015.
  * [17] Jiafa He, Chengwei Pan, Can Yang, Ming Zhang, Yang Wang, Xiaowei Zhou, and Yizhou Yu.  Learning hybrid representations for automatic 3d vessel centerline extraction.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 24–34. Springer, 2020.
  * [18] Qingqi Hong, Qingde Li, Beizhan Wang, Kunhong Liu, Fan Lin, Juncong Lin, Xuan Cheng, Zhihong Zhang, and Ming Zeng.  Accurate geometry modeling of vasculatures using implicit fitting with 2d radial basis functions.  Computer Aided Geometric Design, 62:206–216, 2018.
  * [19] Qingqi Hong, Qingde Li, Beizhan Wang, Jie Tian, Fei Xu, Kunhong Liu, and Xuan Cheng.  High-quality vascular modeling and modification with implicit extrusion surfaces for blood flow computations.  Computer Methods and Programs in Biomedicine, 196:105598, 2020.
  * [20] Clifford M Hurvich, Jeffrey S Simonoff, and Chih-Ling Tsai.  Smoothing parameter selection in nonparametric regression using an improved akaike information criterion.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(2):271–293, 1998.
  * [21] Richard Izzo, David Steinman, Simone Manini, and Luca Antiga.  The vascular modeling toolkit: a python library for the analysis of tubular structures in medical images.  Journal of Open Source Software, 3(25):745, 2018.
  * [22] Tim Jerman, Franjo Pernuš, Boštjan Likar, and Žiga Špiclin.  Beyond frangi: an improved multiscale vesselness filter.  In Medical Imaging 2015: Image Processing, volume 9413, page 94132A. International Society for Optics and Photonics, 2015.
  * [23] Erwan Kerrien, Ahmed Yureidini, Jeremie Dequidt, Christian Duriez, René Anxionnat, and Stéphane Cotin.  Blood vessel modeling for interactive simulation of interventional neuroradiology procedures.  Medical image analysis, 35:685–698, 2017.
  * [24] Marek Kociński, Andrzej Materka, Andreas Deistung, and Jürgen R Reichenbach.  Centerline-based surface modeling of blood-vessel trees in cerebral 3d mra.  In 2016 Signal Processing: Algorithms, Architectures, Arrangements, and Applications (SPA), pages 85–90. IEEE, 2016.
  * [25] David Lesage, Elsa D Angelini, Isabelle Bloch, and Gareth Funka-Lea.  A review of 3d vessel lumen segmentation techniques: Models, features and extraction schemes.  Medical image analysis, 13(6):819–845, 2009.
  * [26] Mark H Longair, Dean A Baker, and J Douglas Armstrong.  Simple neurite tracer: open source software for reconstruction, visualization and analysis of neuronal processes.  Bioinformatics, 27(17):2453–2454, 2011.
  * [27] Odyssée Merveille, Hugues Talbot, Laurent Najman, and Nicolas Passat.  Curvilinear structure analysis by ranking the orientation responses of path operators.  IEEE transactions on pattern analysis and machine intelligence, 40(2):304–317, 2017.
  * [28] Les Piegl and Wayne Tiller.  The NURBS book.  Springer Science & Business Media, 1996.
  * [29] Les A Piegl and Wayne Tiller.  Least-squares b-spline curve approximation with arbitary end derivatives.  Engineering with Computers, 16(2):109–116, 2000.
  * [30] Maria José Melo Ramos-Lima, Ismênia de Carvalho Brasileiro, Tamires Layane de Lima, and Pedro Braga-Neto.  Quality of life after stroke: impact of clinical and sociodemographic factors.  Clinics, 73, 2018.
  * [31] Laura M Sangalli, Piercesare Secchi, Simone Vantini, and Alessandro Veneziani.  A case study in exploratory functional data analysis: geometrical features of the internal carotid artery.  Journal of the American Statistical Association, 104(485):37–48, 2009.
  * [32] Laura M Sangalli, Piercesare Secchi, Simone Vantini, and Alessandro Veneziani.  Efficient estimation of three-dimensional curves and their derivatives by free-knot regression splines, applied to the analysis of inner carotid artery centrelines.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 58(3):285–306, 2009.
  * [33] Khalid M Saqr, Sherif Rashad, Simon Tupin, Kuniyasu Niizuma, Tamer Hassan, Teiji Tominaga, and Makoto Ohta.  What does computational fluid dynamics tell us about intracranial aneurysms? a meta-analysis and critical review.  Journal of Cerebral Blood Flow & Metabolism, 40(5):1021–1039, 2020\.
  * [34] Michiel Schaap, Coert T Metz, Theo van Walsum, Alina G van der Giessen, Annick C Weustink, Nico R Mollet, Christian Bauer, Hrvoje Bogunović, Carlos Castro, Xiang Deng, et al.  Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms.  Medical image analysis, 13(5):701–714, 2009.
  * [35] Gideon Schwarz.  Estimating the dimension of a model.  The annals of statistics, pages 461–464, 1978.
  * [36] Shin-ichiro Sugiyama, Kuniyasu Niizuma, Kenichi Sato, Sherif Rashad, Misaki Kohama, Hidenori Endo, Toshiki Endo, Yasushi Matsumoto, Makoto Ohta, and Teiji Tominaga.  Blood flow into basilar tip aneurysms: a predictor for recanalization after coil embolization.  Stroke, 47(10):2541–2547, 2016.
  * [37] Chaman Singh Verma, Paul F Fischer, Seung E Lee, and Francis Loth.  An all-hex meshing strategy for bifurcation geometries in vascular flow simulation.  In Proceedings of the 14th international meshing roundtable, pages 363–375. Springer, 2005.
  * [38] Vincent Vidal, Guillaume Lavoué, and Florent Dupont.  Low budget and high fidelity relaxed 567-remeshing.  Computers & Graphics, 47:16–23, 2015.
  * [39] Samir Vinchurkar and P Worth Longest.  Evaluation of hexahedral, prismatic and hybrid mesh styles for simulating respiratory aerosol dynamics.  Computers & Fluids, 37(3):317–331, 2008.
  * [40] Susan N Wright, Peter Kochunov, Fernando Mut, Maurizio Bergamino, Kerry M Brown, John C Mazziotta, Arthur W Toga, Juan R Cebral, and Giorgio A Ascoli.  Digital reconstruction and morphometric analysis of human brain arterial vasculature from magnetic resonance angiography.  Neuroimage, 82:170–181, 2013.
  * [41] Guanglei Xiong, Suraj Musuvathy, and Tong Fang.  Automated structured all-quadrilateral and hexahedral meshing of tubular surfaces.  In Proceedings of the 21st international meshing roundtable, pages 103–120. Springer, 2013.
  * [42] Hasballah Zakaria, Anne M Robertson, and Charles W Kerber.  A parametric model for studies of flow in arterial bifurcations.  Annals of biomedical Engineering, 36:1515, 2008.
  * [43] Yongjie Zhang, Yuri Bazilevs, Samrat Goswami, Chandrajit L Bajaj, and Thomas JR Hughes.  Patient-specific vascular nurbs modeling for isogeometric analysis of blood flow.  Computer methods in applied mechanics and engineering, 196(29-30):2943–2959, 2007.