# Modeling and hexahedral meshing of cerebral arterial networks from centerlines

Méghane Decroocq, Carole Frindel, Pierre Rougé, Makoto Ohta, Guillaume Lavoué

► To cite this version:

Méghane Decroocq, Carole Frindel, Pierre Rougé, Makoto Ohta, Guillaume Lavoué. Modeling and hexahedral meshing of cerebral arterial networks from centerlines. Medical Image Analysis, 2023, 89, pp.102912. 10.1016/j.media.2023.102912. hal-04692570

HAL Id: hal-04692570

https://hal.science/hal-04692570v1

Submitted on 28 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Medical Image Analysis

## Modeling and hexahedral meshing of cerebral arterial networks from centerlines

--Manuscript Draft--

<table><tr><td>Manuscript Number:</td><td>MEDIA-D-22-00002R2</td></tr><tr><td>Article Type:</td><td>Research Paper</td></tr><tr><td>Keywords:</td><td>Cerebral Arterial network; Centerlines, Hexahedral mesh; Computational fluid dynamics</td></tr><tr><td>Corresponding Author:</td><td>Méghane Decroocq<br>FRANCE</td></tr><tr><td>First Author:</td><td>Méghane Decroocq</td></tr><tr><td>Order of Authors:</td><td>Méghane Decroocq<br>Carole Frindel<br>Pierre Rougé<br>Makoto Ohta<br>Guillaume Lavoué</td></tr><tr><td>Abstract:</td><td>Computational fluid dynamics (CFD) simulation provides valuable information on blood flow from the vascular geometry. However, it requires extracting 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 a 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 reproducibility.<br>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.<br>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 were meshed without defects needing manual intervention, despite the challenging aspect of the input data. The source code will be released publicly.</td></tr><tr><td>Additional Information:</td><td></td></tr><tr><td>Question</td><td>Response</td></tr><tr><td>To complete your submission you must select a statement which best reflects the availability of your research data/code. IMPORTANT: this statement will be published alongside your article. If you have selected "Other", the explanation text will be published verbatim in your article (online and in the PDF).<br>(If you have not shared data/code and wish to do so, you can still return to Attach</td><td>Other (please explain: e.g. 'I have shared the link to my data/code at the Attach File step').</td></tr></table>
Files. Sharing or referencing research data and code helps other researchers to evaluate your findings, and increases trust in your article. Find a list of supported data repositories in [Author Resources](#), including the free-to-use multidisciplinary open Mendeley Data Repository.)

You have selected other. Please explain: as follow-up to "To complete your submission you must select a statement which best reflects the availability of your research data/code. IMPORTANT: this statement will be published alongside your article. If you have selected "Other", the explanation text will be published verbatim in your article (online and in the PDF).

(If you have not shared data/code and wish to do so, you can still return to Attach Files. Sharing or referencing research data and code helps other researchers to evaluate your findings, and increases trust in your article. Find a list of supported data repositories in [Author Resources](#), including the free-to-use multidisciplinary open Mendeley Data Repository.)"

The code is publicly available. The data that has been used will be made available on request.
**Declaration of interests**

[x] The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[ ] The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
- Modeling the surface of arterial networks from centerlines

- Hexahedral meshing of large networks

- Application to mesh the open database BraVa
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Medical Image Analysis

journal homepage: www.elsevier.com/locate/media

<figure><img src="image_4.png" /></figure>

# Modeling and hexahedral meshing of cerebral arterial networks from centerlines

Méghane Decroocq<sup>a,b,d,e,f</sup>, Carole Frindel<sup>a,e,*</sup>, Pierre Rougé<sup>d</sup>, Makoto Ohta<sup>d,e</sup>, Guillaume Lavoué<sup>b,c</sup>

## ARTICLE INFO

Article history:

2000 MSC: 41A05, 41A10, 65D05, 65D17

Keywords: Cerebral arterial network, Centerlines, Hexahedral mesh, Computational fluid dynamics

## ABSTRACT

Computational fluid dynamics (CFD) simulation provides valuable information on blood flow from the vascular geometry. However, it requires extracting 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 a 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 reproducibility.

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 were meshed without defects needing manual intervention, despite the challenging aspect of the input data. The source code will be released publicly.

© 2022 Elsevier B. V. All rights reserved.

## 1. Introduction

Vascular diseases, such as stroke, can cause severe disability or death (Ramos-Lima et al., 2018). 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, as it provides information on the hemodynamics from the vessel geometry (Saqr et al., 2020; Sugiyama et al., 2016). 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 complex shapes, tetrahe-

<footnote><sup>a</sup>CREATIS, Université Lyon1, CNRS UMR5220, INSERM U1206, INSA-Lyon, 69621 Villeurbanne, France</footnote>

<footnote><sup>b</sup>LIRIS, CNRS UMR 5205, F-69621, France</footnote>

<footnote><sup>c</sup>Ecole Centrale de Lyon, France</footnote>

<footnote><sup>d</sup>ELyTMaX IRL3757, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard Lyon 1, Tohoku University, 980-8577, Sendai, Japan</footnote>

<footnote><sup>e</sup>Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan</footnote>

<footnote><sup>f</sup>Graduate School of Biomedical Engineering, Tohoku University, 6-6 Aramaki-aza-aoba, Aoba-ku, Sendai, Miyagi 980-8579, Japan</footnote>

<footnote>*Corresponding author: Carole Frindel
e-mail: carole.frindel@creatis.insa-lyon.fr (Carole Frindel)</footnote>
dral 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 Vinchurkar and Longest (2008), De Santis et al. (2010) and Ghaffari et al. (2017). 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 (Jerman et al., 2015; Merveille et al., 2017) and to improve the segmentation methods. In particular, the rise of deep learning-based segmentation methods resulted in significant progress in vascular segmentation (Tetteh et al., 2020; Livne et al., 2019). However, the accuracy of the segmentation does not guarantee the accuracy of the mesh it entails (e.g vessels merging due to the image resolution, disconnected vessels, bumps), nor its usability for numerical simulation. Besides, the vascular network 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 other, 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 (Zhang et al., 2021; He et al., 2020). 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 for the construction of big databases of large vascular networks (Wright et al., 2013) or the creation of ideal models. It also offers more editing flexibility than segmentation-based meshes. As it encodes the vessel topology and orientation, it has a high potential for the creation of meshes with high-quality, flow-oriented 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 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.

## 1.1. Contributions

In this work, we propose a method that overcomes the limitations of centerlines to generate a high-quality mesh ready for CFD. Our main contributions to the current state-of-the-art are:

* We introduced an original vessel model and its approximation algorithm based on penalized splines, which enable to model both the spatial coordinates and the radius in a single function and offer good robustness to noise and low-sampling.

* The vessel model is combined with a physiological parametric model of bifurcation proposed by Zakaria et al.

(2008) to form a parametric model of the entire vascular network. A method to extract the parameters of the bifurcation model directly from centerlines is proposed, and the model was generalized to planar n-furcations. If the bifurcation model itself is not new, the use of this type of physiologic bifurcation model, as opposed to geometric bifurcation models, to reconstruct a realistic vascular shape has not been investigated in previous studies.

* A parametric method to create a structured hexahedral volume mesh with flow-oriented cells from vessel and bifurcation models is proposed. It includes relaxation and smoothing steps to improve the quality of the cells without deforming the model shape. This meshing method gives more control over the distribution and density of the cells than the commonly used tetrahedral meshing.

* The model and the mesh are both stored in a single graph structure which enables to easily and inexpensively edit the topology and geometry of the vascular networks.

The proposed framework is fast, fully automatic, and produces high-quality meshes, therefore opening the way to numerical simulation in large networks. It was evaluated qualitatively and quantitatively against other explicit and implicit centerline-based meshing methods, as well as segmentation-based meshing methods. Finally, several practical applications are presented, including the meshing of a large database of 60 large cerebral networks, pathology modeling and hexahedral remeshing.

## 2. Related work

### 2.1. *Segmentation-based meshing*

Segmentation of medical images is the most common method to obtain patient-specific meshes for CFD. In recent years, deep learning-based models have led to significant advances in vascular segmentation. More specifically, convolutional neural networks (CNNs) have achieved very good performances (Jiang et al., 2018; Tetteh et al., 2020). The popular U-net architecture (Ronneberger et al., 2015) has been successfully applied to the segmentation of intracranial vessels in Quon et al. (2020) and Livne et al. (2019). Hilbert et al. (2020) proposed an extended U-net architecture using context aggregation and deep supervision for brain vessel segmentation. Besides, the attention mechanisms have been used to help the network to better learn global dependencies and increase the receptive field in Mou et al. (2021), Ni et al. (2020) and Li et al. (2021).

For medical applications such as CFD, more than the segmentation itself, the smoothness and the topological accuracy of the mesh it entails are critical. However, in the literature, there was very little focus on the conversion of the segmented volumes to mesh. Recently, Wickramasinghe et al. (2020) and Kong et al. (2021) introduced novel neural network architectures to go directly from a 3D volume to a 3D surfaces. Despite those recent advances, the meshing largely relies on algorithms such as the Marching Cubes to produce a surface mesh with tetrahedral elements, followed by a smoothing step
(Watanabe et al., 2018; Misaki et al., 2021). However, This type of segmentation-based meshes commonly suffers from topological inaccuracies (e.g. merging or disconnected vessels, bulges, missing vessels) and requires a burdensome manual post-processing (Glaßer et al., 2015), as we demonstrated in Section 6.2. Such problems are not correctly captured by the image-based metrics (e.g DICE score) used to evaluate the segmentation methods. To overcome those challenges, the centerline representation of vascular networks has recently gained interest.

Some recent segmentation approaches propose to integrate the vessel centerline information to build more topology-oriented metrics. (Keshwani et al., 2020) proposed to segment the vascular network from its skeleton by learning a connectivity metric between center-voxels. Besides, (Shit et al., 2021) introduced a novel topology-preserving loss for the training of neural networks, which relies on the centerlines of the predicted segmentations. The information provided by the centerlines allowed the neural network to improve the topology correctness of the segmentations. In this context, we believe that the use of centerlines in the meshing process can offer many advantages for CFD applications.

## 2.2. Centerline-based meshing

In this part, 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 - and 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 into 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 (Hong et al., 2018), implicit extrusion surfaces (Hong et al., 2020) or local implicit modeling (Kerrien et al., 2017) to reconstruct vascular networks from medical images, or centerlines only for Abdellah et al. (2020). If they stand out by their ability to reconstruct complex branching topology, they do not allow as much control on the final mesh as explicit methods (e.g. hexahedral meshing) which makes them less suited for CFD using the finite element method.

In explicit methods, the vessel surface is obtained by sweeping along the centerlines. 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 Kociński et al. (2016) and Ghaffari et al. (2015), 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 Kociński et al. (2016) and Bezier segments for Ghaffari et al. (2015). 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, Han et al. (2015) and Guo et al. (2013) 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 afterward. Their method guarantees the smoothness of the model but results in unnatural-looking bifurcations. Finally, Zakaria et al. (2008) proposed a physiologic model where the bifurcations are represented by two merged tubes. It was validated with regards to both the accuracy of the anatomy and 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 outlet sections), and shows more anatomical realism. In Zakaria et al. (2008), the authors extracted the model parameters from a surface mesh, they did not suggest a way to extract them from centerlines.

## 2.3. Hexahedral meshing

For applications such as CFD with the finite element method, the inside 3D surface mesh reconstructed - with segmentation-based or centerline-based methods - must be discretized. The shape of the cells in the discretized volume mesh influence the cost and the stability of the numerical simulation. More specifically, hexahedral structured meshes were showed to reduce the cost and increase the convergence and stability of the numerical simulation results, compared to the commonly used tetrahedral meshes (Vinchurkar and Longest, 2008; Ghaffari et al., 2017; De Santis et al., 2010). Hexahedral meshing has been investigated over the past years as it offers many advantages. It better captures 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 (Zhang et al., 2007). Hexahedral meshing is however limited by a far more complex generation process than standard tetrahedral meshes. In the application to the arterial networks, 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 methods were proposed to obtain a robust branch decomposition. De Santis et al. introduced semi-automatic methods, ranging from the manual selection of the most relevant slices of the input surface mesh (De Santis et al. (2010)), user-defined bifurcation coordinate system (De Santis et al., 2011a), to the generation and adjustment of a block-structure representation of the network (De Santis et al., 2011b). Automatic methods are based on Voronoi diagram (Antiga et al., 2002), resolution of the Laplace's equation (Verma et al., 2005), random-walk algorithm (Xiong et al., 2013) or branching templates (Zhang et al., 2007) or parametric models (Ghaffari et al., 2017). The hexahedral meshing can then be created from the decomposition through various techniques; Copper scheme in the work of Antiga et al. (2002), template grid sweeping for
Verma et al. (2005), Zhang et al. (2007) and Ghaffari et al. (2017), Bezier spline modeling followed by an iso-parametric transformation of a template mesh (De Santis et al., 2011a), projection and refinement of block-structures (De Santis et al., 2011b), Laplacian-based harmonic functions combined with Catmull-Clark subdivision (Xiong et al., 2013).

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 De Santis et al. (2011a) 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 Ghaffari et al. (2017) 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 extracted from medical images (both sparse and noisy).

In this work, our purpose is to meet the challenges arising from this state-of-the-art with a framework integrating a modeling and meshing step. The shortcomings of segmentation-based meshing are addressed by developing a method based on centerlines, which integrates topologic information and geometric a priori. The realism of the surface reconstructed from the centerlines is improved compared to other state-of-the-art methods by two main contributions: a centerline approximation method which enables to control the trade-off between smoothness and proximity to the input data points, and a physiological bifurcation model offering more realism than the geometrical bifurcation model commonly used. Finally, an original meshing algorithm is proposed to create obtain high-quality hexahedral meshes suitable for CFD simulations, as opposed to the largely used tetrahedral meshes.

### 3. Input data

The input vessel centerlines we consider are composed of a 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 Izzo et al. (2018). In this work, we use centerlines from two publicly available datasets. The Aneurisk database (Aneurisk-Team, 2012) provides 3D models of the main arteries of the circle of Willis for patients with an aneurism. High-resolution centerlines were extracted from the surface meshes using VMTK software. The BraVa database (Wright et al., 2013) 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 the ImageJ plugin Neurite Tracer (Longair et al., 2011) 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 (Ghaffari et al., 2017), regression splines (Kociński et al., 2016), free knot regression splines or local polynomial smoothing (Sangalli et al., 2009b). Only Sangalli et al. (2009b) gives the detail of the implementation of the approximation method and provides a thorough study of the accuracy of their model regarding the spatial coordinates and the derivatives. The accuracy of both the first and second derivatives is crucial because the vessel curvature impacts the hemodynamics (Sangalli et al., 2009a). Moreover, meshing techniques are often based on the normals of the centerline (Kociński et al., 2016; Ghaffari et al., 2017). It is important to note that the proposed approximation methods (Sangalli et al., 2009b; Kociński et al., 2016; Ghaffari et al., 2017) 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 combining 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, \quad (1)$$

for $u \in [0, 1]$, where $N_{i,p}$ is the ith basis spline function of order p and {$P_0, P_1, ..., Pn-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 Sangalli et al. (2009b). In the other approach, a relatively large number of control points and a uniform knot vector are used and the smoothness is constrained by a penalty on the second derivatives Craven and Wahba (1978); Eilers and Marx (1996). 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 Eilers and Marx (1996). 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, \dots, 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, \quad (2)$$

where $t$ is a time parametrization vector that 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 scales 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)}, \quad (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)$ where $t$ is the time parametrization of each data point and $r$ their radius value:

$$P_{(t,r)} = (N^T N + \lambda_r \Delta)^{-1} N^T D_{(t,r)}. \quad (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.

## 4.2. Bifurcations

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

<figure><img src="image_5.png" /><figcaption>Fig. 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.</figcaption></figure>

### 4.2.1. Zakaria's model

Zakaria et al. (2008) proposed a parametric model for non-planar bifurcations. Their model was validated regarding 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 the reconstruction of bifurcations from sparse data. In this model, bifurcations are created by merging two tubes that represent 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 merge. 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, and 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_{C0}, 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$.

### 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 outlet data points (resp. in deep blue and green in Figure 3 (a)) to form two input centerlines going through the bifurca-
<figure><img src="image_6.png" /><figcaption>Fig. 2. Five cross-sections bifurcation model introduced by Zakaria et al. (2008).</figcaption></figure>

tion, 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><img src="image_7.png" /><figcaption>Fig. 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.</figcaption></figure>

### 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 Piegl and Tiller (2000). We extend this approach to penalized splines and we propose a weaker constraint that 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'(0)$ (respectively $s'(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:

$$ \left\{ \begin{array}{l} 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{array} \right. \qquad (5) $$

where $\alpha$ and $\beta$, the end tangent magnitude, are additional parameters to optimize. Those constraints guarantee 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 into three geometrical branches; one inlet branch and two outlet branches. The Figure 4 (b) gives an example of branch splitting using three separation planes. Antiga et al. Antiga and Steinman (2004) 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 (Antiga and Steinman, 2004). 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 that 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 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_{m1}$ and $p_{m2}$, where $p_{m1,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><img src="image_8.png" /><figcaption>Fig. 4. Geometric decomposition of the bifurcation model. In (b), the end cross-sections are represented in red and the separations planes in black.</figcaption></figure>

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₁ and SP₂. The separation points AP, SP₁ and SP₂ are finally connected to the center points CT₀, CT₁ 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 differences 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 has 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₀, C₁ and C₂. A normalized reference vector $\vec{r}_{fc}$ which minimizes the rotation with the separation points SP₁,₂ 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{r}_{fc}$ 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₁ and SP₂ to both center points CT₀ and CT₁.

<figure><img src="image_9.png" /><figcaption>Fig. 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}$.</figcaption></figure>

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 C₀ are connected to the nodes of the half-sections defined by SP₁ and SP₂, and the nodes of C₁,₂ are connected respectively to the nodes of the half-sections SP₁,₂ 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 that is used to control the topology and geometry of the final mesh grid, but it does not necessarily lie 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$ do 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><img src="image_10.png" /><figcaption>Fig. 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.</figcaption></figure>

The properties of the resulting mesh depend 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 networks; 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. More-
over, a rupture of continuity is observed when the curves cross 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 but it triggers important deformations of the general shape of the model. To avoid deformations, Vidal et al. (2015) 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, measured 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 cross the separation plans. Based on the average quality of the faces, we estimate that 5 relaxation iterations are sufficient.

<figure><img src="image_11.png" /><figcaption>Fig. 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).</figcaption></figure>

### 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 where the two vessels merge. The curvature in the apex regions impacts the pressure and velocity fields obtained by numerical simulation, as shown by Haljasmaa et al. (2001). 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. Zakaria et al. (2008) proposed to smooth the apex region by projecting the nodes on a sphere of a 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 hexagonal 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 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

<figure><img src="image_12.png" /><figcaption>Fig. 8. Illustration of the apex smoothing pipeline.</figcaption></figure>

plane (resp. green and blue arrows on Figure 8). A circle whose radius corresponds 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 the 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><img src="image_13.png" /><figcaption>Fig. 9. Apex smoothing with different radius of curvature values. The colormap encodes the local distance to the original mesh, on the left.</figcaption></figure>

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 Zakaria et al. (2008) 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.

## 5.2. Vessels

For the vessels, we adapted the meshing method proposed by Ghaffari et al. (2017) to the parametric model proposed in
<figure><img src="image_14.png" /><figcaption>Fig. 10. (a) and (b) respectively illustrates the parametric model and the branch decomposition scheme for a trifurcation. (c) shows an example of trifurcation mesh.</figcaption></figure>

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 position, 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 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 Ghaffari et al. (2017). 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 central block. The relative size α, β, γ of the areas, the number Nα of boundary layers and the number Nβ of intermediary layers can be adjusted. The separation planes of the n-furcations are handled by combining n + 1 halves grids. The successive O-grid patterns are connected to form the hexahedral cells of the volume mesh, as shown in Figure 11.

<figure><img src="image_15.png" /><figcaption>Fig. 11. Illustration of the O-grid pattern and volume meshing method.</figcaption></figure>

## 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 are assessed in a comparative study performed on a synthetic dataset of distorted centerlines. Then, our meshing pipeline is compared quantitatively and qualitatively with two concurrent state-of-the-art methods: deep learning-based segmentation and implicit meshing. Finally, we provide additional performance indicators 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 it does not include pathologies but goes 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 constitutes 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 distorted 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 σ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 σradius and added to the ground truth radius. In both cases, the applied standard deviation value is proportional to the point radius, as indicated in Table 1, in order to keep 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><caption>Table 1. Parameters used for the distortion of the ground truth centerlines</caption><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>density (mm<sup>-1</sup>)</td><td>2</td><td>4</td><td>10</td><td>16</td><td>20</td></tr><tr><td>σ<sub>radius</sub>(mm)</td><td>0.01r</td><td>0.05r</td><td>0.1r</td><td>0.3r</td><td>0.5r</td></tr><tr><td>σ<sub>spatial</sub>(mm)</td><td>0.01r</td><td>0.05r</td><td>0.1r</td><td>0.3r</td><td>0.5r</td></tr></tbody></table>
For each density value in Table 1, ten combinations 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 null. Each noise combination is repeated three times to account for the 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 distorted dataset to 600.

### 6.1.2. Approximation methods

To demonstrate the robustness and the accuracy of the approximation strategy used to reconstruct the surface of the vessels presented in Section 4.1, we compared it to other explicit centerline-based meshing methods of the literature (Kociński et al., 2016; Ghaffari et al., 2017). Those methods also rely on splines or Bezier curves to approximate the centerline data points and reconstruct the vessel surface. As most of the authors did not provide the details of the fitting method employed, we implemented four commonly used spline-based approximations methods with incremental complexity 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 the 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 (Akaike (1973)) $AIC_1$:

$$AIC_1 = m \times \log(\text{SSE}) + 8(n + p). \quad (6)$$

where $m$ is the number of data points, $p$ is the degree of the spline, $n$ is 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 Eilers and Marx (1996). 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 us 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 Eilers and Marx (1996):

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

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

Unlike $AIC_1$, it is not employed to choose an optimal number of control points but to select the optimal value for the smoothing parameter $\lambda$. This criterion was compared to other criteria of the literature; the corrected Akaike information criterion of Hurvich et al. (1998), the Schwarz's Bayesian criterion of Schwarz (1978), the cross-validation criterion and generalized cross-validation criterion of Craven and Wahba (1978). 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><img src="image_16.png" /><figcaption>Fig. 12. Matching time parameters by minimum distance projection from $s$ onto $\hat{s}$</figcaption></figure>

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 RMSEradius (respectively RMSEspatial) the root mean squared error of the radius (respectively the spatial coordinates). 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 RMSEderspatial and RMSEderradius. As curvature is commonly considered in hemodynamic studies, the model performance concerning the centerline curvature is also measured, by the metric RMSEcurv. Finally, the length of the vessel affects the delay of blood arrival between the inlet and the outlet of the vascular tree in numerical simulations. Therefore, the difference $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 a surprisingly high spatial error. The overfitting and noise problems are partially solved by optimizing the number of control points with the 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><img src="image_17.png" /><figcaption>Fig. 13. Mesh resulting from the approximation of distorted data (density = $1mm^{-1}$, $\sigma_{radius}$ = 0.1) by three of the methods compared in section 6.1.</figcaption></figure>

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-AIC is very similar to the ground truth vessel, but in return, the trajectory of the centerline was too smoothed. On the other hand, both the radius and trajectory of the vessel produced with SRP-AIC are closer to the ground truth. As a conclusion, the proposed approximation method shows 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.

### 6.2. Comparison with state-of-the-art methods

In this section, meshes obtained with our method are visually and quantitatively compared to meshes produced by state-of-the-art deep learning-based segmentation methods (Tetteh et al., 2020; Livne et al., 2019), as well as a recent implicit centerline-based meshing method (Abdellah et al., 2020).

<table><caption>Table 2. Overall evaluation of the approximation methods : mean values of the quality criteria for all the centerlines distorted by radius noise addition. The cells in gray correspond to the lowest error for each metric.</caption><thead><tr><th></th><th>GNP</th><th>GNP-AIC</th><th>GP-AIC</th><th>SRP-AIC</th></tr></thead><tbody><tr><td>RMSE<sub>spatial</sub></td><td>8.462</td><td>0.034</td><td>0.053</td><td>0.029</td></tr><tr><td>RMSE<sub>radius</sub></td><td>17.523</td><td>0.095</td><td>0.042</td><td>0.043</td></tr><tr><td>RMSE<sub>der</sub><sub>spatial</sub></td><td>0.218</td><td>0.118</td><td>0.042</td><td>0.009</td></tr><tr><td>RMSE<sub>der</sub><sub>radius</sub></td><td>0.391</td><td>0.214</td><td>0.032</td><td>0.032</td></tr><tr><td>RMSE<sub>curv</sub></td><td>1919.428</td><td>190.531</td><td>0.060</td><td>0.035</td></tr><tr><td>L<sub>diff</sub></td><td>718.906</td><td>0.057</td><td>0.207</td><td>0.004</td></tr></tbody></table>

<table><caption>Table 3. Overall evaluation of the approximation methods : mean values of the quality criteria for all the centerlines distorted by spatial noise addition. The cells in gray correspond to the lowest error for each metric.</caption><thead><tr><th></th><th>GNP</th><th>GNP-AIC</th><th>GP-AIC</th><th>SRP-AIC</th></tr></thead><tbody><tr><td>RMSE<sub>spatial</sub></td><td>0.511</td><td>0.152</td><td>0.099</td><td>0.096</td></tr><tr><td>RMSE<sub>radius</sub></td><td>0.008</td><td>0.009</td><td>0.018</td><td>0.007</td></tr><tr><td>RMSE<sub>der</sub><sub>spatial</sub></td><td>0.314</td><td>0.343</td><td>0.075</td><td>0.076</td></tr><tr><td>RMSE<sub>der</sub><sub>radius</sub></td><td>0.015</td><td>0.019</td><td>0.021</td><td>0.013</td></tr><tr><td>RMSE<sub>curv</sub></td><td>1.524</td><td>2.362</td><td>0.085</td><td>0.091</td></tr><tr><td>L<sub>diff</sub></td><td>50.180</td><td>15.071</td><td>0.252</td><td>0.207</td></tr></tbody></table>

#### 6.2.1. Comparison pipeline

Centerlines can be extracted either from the grayscale image directly or from a segmented image. In this way, centerline-based meshing can be used either as a substitute or a complement to segmentation. In our comparison study, we investigated both approaches, as illustrated in Figure 14. In what follows, the centerlines manually extracted by experts from MRA images of the BraVa database (Wright et al., 2013) are considered as reference "expert centerlines".

For the first part of our comparison pipeline (in blue in Figure 14), the MRA images are segmented by state-of-the-art segmentation methods (Tetteh et al., 2020; Livne et al., 2019). A surface mesh is produced from the segmented images by the marching cube algorithm and smoothed using a Taubin filter. The very small components of this mesh are removed to keep only the largest connected parts. A set of centerlines, which we call "segmentation-based centerlines", are extracted from the segmentation. In Section 6.2.3, they are quantitatively compared to the expert centerlines to evaluate their accuracy, and thus the topological and geometrical correctness of the segmentation-based mesh.

For the second part of our comparison pipeline (in red in Figure 14), both the expert centerlines and the segmentation-based centerlines are used as input for the centerline-based meshing methods. To improve their quality and to match the input requirements of our meshing method, the segmentation-based centerlines underwent some automatic post-processing before meshing; the small ending segments are cut out, the cycles are removed by computing a maximum spanning tree of the network, and the edges connecting the data points are re-oriented in the flow direction.

Two centerline-based methods are used to create meshes from those post-processed centerlines and the expert centerline; our explicit meshing method and the method of Abdellah
<figure><img src="image_18.png" /><figcaption>Fig. 14. Pipeline used to compare segmentation-based meshing and centerline-based meshing.</figcaption></figure>

et al. (2020), which is based on implicit structures and for which we used the Blender plug-in provided by the authors, in meta-balls reconstruction mode. In Section 6.2.4, those centerline-based meshes are visually compared to the segmentation-based meshes produced by deep learning methods Tetteh et al. (2020); Livne et al. (2019).

### 6.2.2. Datasets and segmentation models training

Two databases of whole-brain MRA images for healthy patients are used for this evaluation: The BraVa database gives access to 62 expert centerlines and 50 MRA, and the TubeTK database Bullitt et al. (2005) is composed of 34 MRA and the 34 in-house expert segmentations associated. These databases were chosen because they include the same type of images while offering complementary information (resp. expert centerlines and expert segmentation). For the vessels segmentation, two state-of-the-art methods were implemented, both based on neural networks: Deep-VesselNet Tetteh et al. (2020) and U-net Livne et al. (2019). The DeepVesselNet architecture provided by the author was used (https://github.com/giesekow/deepvesselnet), and the U-net neural network was re-implemented. Both models were trained on the expert segmentations of the TubeTK database. 27 images were included in the training set and 7 in the test set.

The loss function used during training was a combination of dice loss and cross-entropy loss. The stochastic gradient descent algorithm was used for the optimization, with a learning rate of 0.01 for U-net and 0.001 for DeepVesselNet. The batch size was set to 5 for U-net and 10 for DeepVesselNet due to memory constraints. These hyperparameters were set empirically by testing a large selection of values for each hyperparameter. The U-net was trained for 200 epochs and DeepVesselNet for 300 epochs.

For this comparison study, we considered the segmentations produced for the 7 images of the test set of TubeTK and the 50 MRA of the BraVa database - segmented with the models trained on TubeTK -. Table 4 summarizes the data type and the number of patients of the different datasets created for this comparison study.

<table><caption>Table 4. Description of the different datasets used in our comparison study. The name of the original database (BraVa or TubeTK), the method used, the nature of the data and the number of patients are given.</caption><thead><tr><th>Database</th><th>Method</th><th>Data</th><th>nPatients</th></tr></thead><tbody><tr><td>BraVa</td><td>DeepVesselNet</td><td>Segmentation</td><td>50</td></tr><tr><td>BraVa</td><td>Expert</td><td>Centerlines</td><td>62</td></tr><tr><td>BraVa</td><td>Unet</td><td>Segmentation</td><td>50</td></tr><tr><td>TubeTK</td><td>DeepVesselNet</td><td>Segmentation</td><td>7</td></tr><tr><td>TubeTK</td><td>Expert</td><td>Segmentation</td><td>34</td></tr><tr><td>TubeTK</td><td>Unet</td><td>Segmentation</td><td>7</td></tr></tbody></table>

### 6.2.3. Quantitative evaluation

In this section, we present the results of the quantitative evaluation of the segmentation-based centerlines compared to the expert centerlines according to the 9 topologic and geometric features described hereafter. *nBulges* corresponds to the number of bulges in the geometry. It is obtained by counting the number of ending segments smaller than the vessel diameter. *nBranch* is the number of branches in the entire network. *nCC* is the number of connected components and *nBranchMaxCC* is the size of the largest connected component, evaluated by its number of branches. This metric highlights the disconnected vessels and small isolated parts in the mesh. *nCycle* is the number of cycles of the network. The only cycle in the cerebral vascular system is the circle of Willis, so the number of cycles should be either 1 - complete circle of Willis - or 0 - incomplete circle of Willis -. Finally, the branching topology of the network is analyzed via the number of bifurcations *nBif*, the number of trifurcations or more *nTrif+* and the minimum (resp. maximum) furcation degree *minDeg* (resp. *maxDeg*), i.e the number of in and out branches (bifurcations = 3). These metrics are reported in Table 5 for the different datasets considered.

We observe in Table 5 that the expert centerlines (in gray) does not have any small ending segments (*nBulges* = 0). The network forms a single connected component, they are no isolated vessels and no cycles. Besides, the branchings are mainly bifurcations, as expected in the cerebral vascular system where trifurcations are rare. The expert centerlines show no branching with a degree superior to 4 (= trifurcations).

On the other side, segmentation-based meshes present a lot of bulges (> 12) and cycles (> 26), mainly because closed vessels are merged in the resulting mesh. The number of trifurcation and higher degree branching is high (> 15), and furcations with up to 6 branches were observed. These metrics bring in light some inaccuracies in the topology of the meshes produced by segmentation, which will affect the mesh geometry and therefore the CFD simulation results. It is interesting to see that such problems - disconnected or merged vessels, bulges - are observed even in the meshes based on the ground truth segmentation made by medical doctors (see "TubeTK expert" row in Table 5). They are not only caused by the segmentation method but also by the nature of the meshing process itself, as it relies on the segmentation of low-resolution images. Moreover, no distinction can be made between veins and arteries in the segmentation process, which might cause peculiar topology in the network. To run numerical simulations in such segmentation-based meshes, important post-processing is required to isolate
<table><caption>Table 5. Topologic and geometric features of the segmentation-based meshes for the different datasets. For each dataset, the median value between all patients is given.</caption><thead><tr><th>Database</th><th>Method</th><th>nBulges</th><th>nBranch</th><th>nCC</th><th>nBranchMaxCC</th><th>nCycle</th><th>nBif</th><th>nTrif+</th><th>minDeg</th><th>maxDeg</th></tr></thead><tbody><tr><td>BraVa</td><td>DeepVesselNet</td><td>43</td><td>369</td><td>451</td><td></td><td>59</td><td>26</td><td>76</td><td>15</td><td>0</td><td>5</td></tr><tr><td>BraVa</td><td>Expert</td><td>0</td><td>205</td><td>1</td><td></td><td>205</td><td>0</td><td>102</td><td>1</td><td>1</td><td>4</td></tr><tr><td>BraVa</td><td>Unet</td><td>38</td><td>504</td><td>380</td><td></td><td>200</td><td>44</td><td>136</td><td>30</td><td>0</td><td>6</td></tr><tr><td>TubeTK</td><td>DeepVesselNet</td><td>52</td><td>552</td><td>557</td><td></td><td>212</td><td>46</td><td>140</td><td>28</td><td>0</td><td>6</td></tr><tr><td>TubeTK</td><td>Expert</td><td>12</td><td>551</td><td>26</td><td></td><td>508</td><td>98</td><td>230</td><td>48</td><td>0</td><td>6</td></tr><tr><td>TubeTK</td><td>Unet</td><td>34</td><td>626</td><td>300</td><td></td><td>446</td><td>85</td><td>215</td><td>46</td><td>0</td><td>6</td></tr></tbody></table>

the arterial system and reconnect or separate vessels. The topologic problems highlighted here are illustrated by enhanced visualizations of the meshes (Figures 15 and 16) in the next section.

### 6.2.4. Visual evaluation

In this part, the results from our method are visually compared against centerline-based meshes produced by the method of Abdellah et al. (2020) and the segmentation-based meshes produced by the method of Tetteh et al. (2020) and Livne et al. (2019). A BraVa patient was selected for this visual evaluation to have access to the expert centerlines associated. We selected a segmentation made with Unet as it provided better results than DeepVesselNet on this database. Figure 15 shows the whole brain meshes obtained by different methods, and Figure 16 shows enhanced visualization of some relevant parts. As shown in Figure 15 (a), the deep learning segmentation algorithm demonstrates a good ability to segment a large part of the vascular network, including small segments. However, the algorithm does not guarantee the connectivity of the network, and some post-processing filtering is needed to remove the small isolated parts (Fig. 15 (b)). The centerline-based methods (images (d) and (e)) were able to produce meshes with a topology similar to the segmentation-based mesh from the centerlines automatically extracted from it. We can see in this figure that the mesh produced by our method is smoother and more geometrically and physiologically accurate than the other meshes for the same network. The geometric quality of the meshes obtained by different methods will be further discussed below. As illustrated in the image (f), the manually extracted expert centerline allowed to reconstruct a larger arterial network with smaller vessels than the centerlines based on the segmentation results. Note that, overall, the radius of the expert centerlines is smaller, due to the extraction method Longair et al. (2011).

In View 1 of Figure 16, we observe that our method, which relies on the vessel tubularity assumption, efficiently cleaned the vascular network from the bulges observed in the segmentation-based mesh. The radius and trajectory smoothing allows for reconstructing the disconnected parts in a very natural way. View 2 highlights the merging vessels and cycles observed in segmentation-based meshes. Our automatic post-treatment of the centerlines allowed to remove unwanted cycles in the network. The implicit method VessMorphoVis Abdellah et al. (2020) offer good flexibility to mesh complex geometry such as pathology or complex branching patterns, however, the surface of the vessels looks bumpy and irregular. Moreover, as shown in both views of the Figure 16, this method appears to

be sensitive to noise on the centerline geometry and radius. As it is robust to noise and data sparsity, our method improves the smoothness and the realism of the vascular geometry.

Overall, those results demonstrate the ability of our algorithm to produce high-quality meshes not only from manually extracted centerlines but also to integrate fully-automated pipelines.

### 6.3. Mesh quality

In CFD, the accuracy and stability of the simulation is affected by the quality of the mesh. 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 are given in Figure 17, together with an example branch.

As shown in the image (a) of Figure 17, 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 a scaled Jacobian value higher than 0.9. The vessel cells have 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 on average on the distributions given for three large cerebral networks in Ghaffari et al. (2017). This proportion goes up to 62% of the cells of the abdominal aortic artery geometry meshed by the method of Xiong et al. (2013). Finally, in De Santis et al. (2011b), 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 De Santis et al. (2011b) and Xiong et al. (2013) focuses on arterial geometries that differ from our study.

### 6.4. Computation time

The computational time of the modeling and meshing steps for five patients of the BraVa database was computed. The re-
<figure><img src="image_19.png" /><figcaption>Fig. 15. Meshes produced with segmentation-based and centerline-based methods for a patient of the BraVa database. From left to right and top to bottom: original mesh created from the segmentation with Unet, the same mesh after filtering of the smallest components and smoothing; centerlines extracted from the mesh after post-processing; mesh produced by VessMorphoVis from these centerlines; mesh produced by our method from these centerlines; mesh produced by our method from manually extracted expert centerlines.</figcaption></figure>

sults for three of them are given in Table 6. 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 6 for different cell densities. The average meshing time goes from 24.6 minutes for a 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.

Besides, a large part of the meshing time corresponds to the computation of the surface nodes; on 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 Table 6 by a factor of 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:
<figure><img src="image_20.png" /><figcaption>Fig. 16. Enhanced visualization of the segmentation-based mesh, the centerline-based mesh obtained with our method and the centerline-based mesh obtained with VessMorphoVis. The red squares highlight the parts where the meshes show important differences.</figcaption></figure>

<table><caption>Table 6. Computational time required to model and mesh large vascular networks from the BraVa dataset.</caption><thead><tr><th>id</th><th>furcation (#)</th><th>vessel (#)</th><th>data point (#)</th><th>modeling time (min)</th><th>cells (#)</th><th>meshing time (min)</th></tr></thead><tbody><tr><td rowspan="4">P1</td><td rowspan="4">96</td><td rowspan="4">194</td><td rowspan="4">2816</td><td rowspan="4">11.3</td><td>1389k</td><td>20.4</td></tr><tr><td>1853k</td><td>25.7</td></tr><tr><td>2316k</td><td>31.2</td></tr><tr><td>2779k</td><td>38.4</td></tr><tr><td rowspan="4">P2</td><td rowspan="4">101</td><td rowspan="4">203</td><td rowspan="4">3531</td><td rowspan="4">18.3</td><td>1916k</td><td>27.5</td></tr><tr><td>2555k</td><td>38.2</td></tr><tr><td>3193k</td><td>49.1</td></tr><tr><td>3832k</td><td>67.5</td></tr><tr><td rowspan="4">P3</td><td rowspan="4">107</td><td rowspan="4">216</td><td rowspan="4">3474</td><td rowspan="4">16.8</td><td>1737k</td><td>26.3</td></tr><tr><td>2316k</td><td>36.1</td></tr><tr><td>2895k</td><td>44.2</td></tr><tr><td>3474k</td><td>55.9</td></tr></tbody></table>

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

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

3. Deform the tubular mesh to match the original surface.

Figure 18 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.

As shown in Figure 18, 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 19.

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 (Cornelissen et al., 2018; Alnæs et al., 2007). 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 the trajectory of a vessel can be modified effortlessly. Figure 20 provides examples of such modifications. Topological and geometrical information is 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 20. Branches and arterial territories inexpensively with a local re-computation of the model and mesh parts.

## 7.3. Large cerebral arterial network meshing

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 21. This dataset is considered 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 22 shows the high level of noise encountered in the input data, both in the radius estimation and the spatial positions. Besides, by computing the
<figure><img src="image_21.png" /><figcaption>**Fig. 17.** 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 vessel's cells. The image (a) illustrates the location of high and low-quality cells within a mesh.</figcaption></figure>

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 be $0.45 \text{ 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 a 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 in image (b) of Figure 22. 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 vessels to go backward from the direction of the flow with a sharp angle. As we use an

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** 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 is shown.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 19.** Cross-section pattern before and after deformation. The corresponding slice of the target surface mesh is represented in grey.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 20.** Editing of a model of the basilar artery using our framework. The bifurcation angle and the radius of the original vertebral artery are modified, and one of the 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.</figcaption></figure>

oriented bifurcation model, it failed to correctly represent the geometry. The image (a) of Figure 22, 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 22
<figure><img src="image_25.png" /><figcaption>Fig. 21. Top view of 8 meshes among the 60 generated from the patients of the BraVa database. The bifurcations where the meshing algorithm has failed (i.e at least one 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.</figcaption></figure>

(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.

## 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 to noise, accurate and automatic. It opens to way to CFD simulations in large vascular networks manually or semi-automatically extracted by medical doctors, with 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 the 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 originally developed for cerebral vascular networks, and

non-planar n-furcations ($n > 3$) that are common in other vessels (e.g aorta, lung vessels) were not addressed yet, which might limit its use. In addition, the robustness of the modeling and meshing method could be further improved as it failed in 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. Besides, we want to emphasize that our objective with this work was 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 some manual post-treatment may still be required before simulation. In this way, our framework offers more editing flexibility than other meshing methods. To take advantage of this flexibility, we developed a vascular network editing software, with a user-friendly interface. This interface integrates the modeling and meshing methods described in this article as well as other editing functionalities such as centerline editing, branch removal or angle modification. It opens 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 to the development apex
<figure><img src="image_26.png" /><figcaption>Fig. 22. 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 a 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 a highlight on the relevant parts. On the right image, the original centerline data points are overlayed on the original MRA image.</figcaption></figure>

smoothing method. Founding is acknowledged through AURA region (SIMAVC project) and ANR 20-CE45-0011 (PreSpin project).

References

Abdellah, M., Guerrero, N.R., Lapere, S., Coggan, J.S., Keller, D., Coste, B., Dagar, S., Courcol, J.D., Markram, H., Schürmann, F., 2020. Interactive visualization and analysis of morphological skeletons of brain vasculature networks with vessmorphovis. Bioinformatics 36, i534–i541.

Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. Proceedings of the Second International Symposium on Information Theory, 267–281.

Alnæs, M.S., Isaksen, J., Mardal, K.A., Romner, B., Morgan, M.K., IngebrigtSEN, T., 2007. Computation of hemodynamics in the circle of willis. Stroke 38, 2500–2505.

Aneurisk-Team, 2012. AneuriskWeb project website, http://ecm2.mathcs. memory.edu/aneuriskweb. Web Site. URL: http://ecm2.mathcs. memory.edu/aneuriskweb.

Antiga, L., Ene-Iordache, B., Caverni, L., Cornalba, G.P., Remuzzi, A., 2002. Geometric reconstruction for computational mesh generation of arterial bifurcations from ct angiography. Computerized Medical Imaging and Graph- ics 26, 227–235.

Antiga, L., Steinman, D.A., 2004. Robust and objective decomposition and mapping of bifurcating vessels. IEEE transactions on medical imaging 23, 704–713.

Bullitt, E., Zeng, D., Gerig, G., Aylward, S., Joshi, S., Smith, J.K., Lin, W., Ewend, M.G., 2005. Vessel tortuosity and brain tumor malignancy: a blinded study1. Academic radiology 12, 1232–1240.

Cornelissen, B., Schneiders, J., Sprengers, M., van den Berg, R., Van Ooij, P., Nederveen, A., Van Bavel, E., Vandertop, W., Slump, C., Marquering, H., et al., 2018. Aneurysmal parent artery–specific inflow conditions for complete and incomplete circle of willis configurations. American journal of neuroradiology 39, 910–915.

Craven, P., Wahba, G., 1978. Smoothing noisy data with spline functions. Numerische mathematik 31, 377–403.

De Santis, G., De Beule, M., Segers, P., Verdonck, P., Verhegghe, B., 2011a. 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, 797–802.

De Santis, G., De Beule, M., Van Canneyt, K., Segers, P., Verdonck, P., Verhegghe, B., 2011b. Full-hexahedral structured meshing for image-based computational vascular modeling. Medical engineering & physics 33, 1318–1325.

De Santis, G., Mortier, P., De Beule, M., Segers, P., Verdonck, P., Verhegghe, B., 2010. Patient-specific computational fluid dynamics: structured mesh generation from coronary angiography. Medical & biological engineering & computing 48, 371–380.

Eilers, P.H., Marx, B.D., 1996. Flexible smoothing with b-splines and penalties. Statistical science 11, 89–121.

Ghaffari, M., Hsu, C.Y., Linninger, A.A., 2015. Automatic reconstruction and generation of structured hexahedral mesh for non-planar bifurcations in vascular networks, in: Computer Aided Chemical Engineering. Elsevier. volume 37, pp. 635–640.

Ghaffari, M., Tangen, K., Alaraj, A., Du, X., Charbel, F.T., Linninger, A.A., 2017. 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.

Glaßer, S., Berg, P., Neugebauer, M., Preim, B., 2015. Reconstruction of 3d surface meshes for blood flow simulations of intracranial aneurysms, in: Proceedings of the Conference of the German Society for Computer and Robotic Assisted Surgery, pp. 163–168.

Guo, J., Li, S., Chui, Y.P., Qin, J., Heng, P.A., 2013. Mesh quality oriented 3d geometric vascular modeling based on parallel transport frame. Computers in biology and medicine 43, 879–888.

Haljasmaa, I., Robertson, A., Galdi, G., 2001. 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, 499–520.

Han, X., Bibb, R., Harris, R., 2015. Design of bifurcation junctions in artificial vascular vessels additively manufactured for skin tissue engineering. Journal of Visual Languages & Computing 28, 238–249.

He, J., Pan, C., Yang, C., Zhang, M., Wang, Y., Zhou, X., Yu, Y., 2020. Learning hybrid representations for automatic 3d vessel centerline extraction, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 24–34.

Hilbert, A., Madai, V.I., Akay, E.M., Aydin, O.U., Behland, J., Sobesky, J., Galinovic, I., Khalil, A.A., Taha, A.A., Wuerfel, J., et al., 2020. BraveNet: fully automated arterial brain vessel segmentation in patients with cerebrovascular disease. Frontiers in artificial intelligence, 78.

Hong, Q., Li, Q., Wang, B., Liu, K., Lin, F., Lin, J., Cheng, X., Zhang, Z., Zeng, M., 2018. Accurate geometry modeling of vasculatures using implicit fitting with 2d radial basis functions. Computer Aided Geometric Design 62, 206–216.

Hong, Q., Li, Q., Wang, B., Tian, J., Xu, F., Liu, K., Cheng, X., 2020. High-quality vascular modeling and modification with implicit extrusion surfaces for blood flow computations. Computer Methods and Programs in Biomedicine 196, 105598.
Hurvich, C.M., Simonoff, J.S., Tsai, C.L., 1998. Smoothing parameter selection in nonparametric regression using an improved akaike information criterion. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 60, 271–293.

Izzo, R., Steinman, D., Manini, S., Antiga, L., 2018. The vascular modeling toolkit: a python library for the analysis of tubular structures in medical images. Journal of Open Source Software 3, 745.

Jerman, T., Pernuš, F., Likar, B., Špiclin, Ž., 2015. Beyond frangi: an improved multiscale vesselness filter, in: Medical Imaging 2015: Image Processing, International Society for Optics and Photonics, p. 94132A.

Jiang, Z., Zhang, H., Wang, Y., Ko, S.B., 2018. Retinal blood vessel segmentation using fully convolutional network with transfer learning. Computerized Medical Imaging and Graphics 68, 1–15.

Kerrien, E., Yureidini, A., Dequidt, J., Duriez, C., Anxionnat, R., Cotin, S., 2017. Blood vessel modeling for interactive simulation of interventional neuroradiology procedures. Medical image analysis 35, 685–698.

Keshwani, D., Kitamura, Y., Ihara, S., Iizuka, S., Simo-Serra, E., 2020. Top-net: Topology preserving metric learning for vessel tree reconstruction and labelling, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, pp. 14–23.

Kociński, M., Materka, A., Deistung, A., Reichenbach, J.R., 2016. Centerline-based surface modeling of blood-vessel trees in cerebral 3d mra, in: 2016 Signal Processing: Algorithms, Architectures, Arrangements, and Applications (SPA), IEEE, pp. 85–90.

Kong, F., Wilson, N., Shadden, S., 2021. A deep-learning approach for direct whole-heart mesh reconstruction. Medical image analysis 74, 102222.

Li, Y., Yang, J., Ni, J., Elazab, A., Wu, J., 2021. Ta-net: Triple attention network for medical image segmentation. Computers in Biology and Medicine 137, 104836.

Livne, M., Rieger, J., Aydin, O.U., Taha, A.A., Akay, E.M., Kossen, T., Sobesky, J., Kelleher, J.D., Hildebrand, K., Frey, D., et al., 2019. A u-net deep learning framework for high performance vessel segmentation in patients with cerebrovascular disease. Frontiers in neuroscience 13, 97.

Longair, M.H., Baker, D.A., Armstrong, J.D., 2011. Simple neurite tracer: open source software for reconstruction, visualization and analysis of neuronal processes. Bioinformatics 27, 2453–2454.

Merveille, O., Talbot, H., Najman, L., Passat, N., 2017. Curvilinear structure analysis by ranking the orientation responses of path operators. IEEE transactions on pattern analysis and machine intelligence 40, 304–317.

Misaki, K., Futami, K., Uno, T., Nambu, I., Yoshikawa, A., Kamide, T., Nakada, M., 2021. Inflow hemodynamics of intracranial aneurysms: A comparison of computational fluid dynamics and 4d flow magnetic resonance imaging. Journal of Stroke and Cerebrovascular Diseases 30, 105685.

Mou, L., Zhao, Y., Fu, H., Liu, Y., Cheng, J., Zheng, Y., Su, P., Yang, J., Chen, L., Frangi, A.F., et al., 2021. Cs2-net: Deep learning segmentation of curvilinear structures in medical imaging. Medical image analysis 67, 101874.

Ni, J., Wu, J., Wang, H., Tong, J., Chen, Z., Wong, K.K., Abbott, D., 2020. Global channel attention networks for intracranial vessel segmentation. Computers in biology and medicine 118, 103639.

Piegl, L.A., Tiller, W., 2000. Least-squares b-spline curve approximation with arbitrary end derivatives. Engineering with Computers 16, 109–116.

Quon, J.L., Chen, L.C., Kim, L., Grant, G.A., Edwards, M.S., Cheshier, S.H., Yeom, K.W., 2020. Deep learning for automated delineation of pediatric cerebral arteries on pre-operative brain magnetic resonance imaging. Frontiers in Surgery, 89.

Ramos-Lima, M.J.M., Brasileiro, I.d.C., Lima, T.L.d., Braga-Neto, P., 2018. Quality of life after stroke: impact of clinical and sociodemographic factors. Clinics 73.

Ronneberger, O., Fischer, P., Brox, T., 2015. U-net: Convolutional networks for biomedical image segmentation, in: International Conference on Medical image computing and computer-assisted intervention, Springer, pp. 234–241.

Sangalli, L.M., Secchi, P., Vantini, S., Veneziani, A., 2009a. A case study in exploratory functional data analysis: geometrical features of the internal carotid artery. Journal of the American Statistical Association 104, 37–48.

Sangalli, L.M., Secchi, P., Vantini, S., Veneziani, A., 2009b. 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, 285–306.

Saqr, K.M., Rashad, S., Tupin, S., Niizuma, K., Hassan, T., Tominaga, T., Ohta, M., 2020. What does computational fluid dynamics tell us about intracranial aneurysms? a meta-analysis and critical review. Journal of Cerebral Blood Flow & Metabolism 40, 1021–1039.

Schwarz, G., 1978. Estimating the dimension of a model. The annals of statistics, 461–464.

Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P., Bauer, U., Menze, B.H., 2021. cldice-a novel topology-preserving loss function for tubular structure segmentation, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 16560–16569.

Sugiyama, S.i., Niizuma, K., Sato, K., Rashad, S., Kohama, M., Endo, H., Endo, T., Matsumoto, Y., Ohta, M., Tominaga, T., 2016. Blood flow into basilar tip aneurysms: a predictor for recanalization after coil embolization. Stroke 47, 2541–2547.

Tetteh, G., Efremov, V., Forkert, N.D., Schneider, M., Kirschke, J., Weber, B., Zimmer, C., Piraud, M., Menze, B.H., 2020. Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-d angiographic volumes. Frontiers in Neuroscience, 1285.

Verma, C.S., Fischer, P.F., Lee, S.E., Loth, F., 2005. An all-hex meshing strategy for bifurcation geometries in vascular flow simulation, in: Proceedings of the 14th international meshing roundtable, Springer, pp. 363–375.

Vidal, V., Lavoué, G., Dupont, F., 2015. Low budget and high fidelity relaxed 567-remeshing. Computers & Graphics 47, 16–23.

Vinchurkar, S., Longest, P.W., 2008. Evaluation of hexahedral, prismatic and hybrid mesh styles for simulating respiratory aerosol dynamics. Computers & Fluids 37, 317–331.

Watanabe, T., Isoda, H., Takehara, Y., Terada, M., Naito, T., Kosugi, T., Onishi, Y., Tanoi, C., Izumi, T., 2018. Hemodynamic vascular biomarkers for initiation of paraclinoid internal carotid artery aneurysms using patient-specific computational fluid dynamic simulation based on magnetic resonance imaging. Neuroradiology 60, 545–555.

Wickramasinghe, U., Remelli, E., Knott, G., Fua, P., 2020. Voxel2mesh: 3d mesh model generation from volumetric data, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, pp. 299–308.

Wright, S.N., Kochunov, P., Mut, F., Bergamino, M., Brown, K.M., Mazziotta, J.C., Toga, A.W., Cebral, J.R., Ascoli, G.A., 2013. Digital reconstruction and morphometric analysis of human brain arterial vasculature from magnetic resonance angiography. Neuroimage 82, 170–181.

Xiong, G., Musuvathy, S., Fang, T., 2013. Automated structured all-quadrilateral and hexahedral meshing of tubular surfaces, in: Proceedings of the 21st international meshing roundtable. Springer, pp. 103–120.

Zakaria, H., Robertson, A.M., Kerber, C.W., 2008. A parametric model for studies of flow in arterial bifurcations. Annals of biomedical Engineering 36, 1515.

Zhang, Y., Bazilevs, Y., Goswami, S., Bajaj, C.L., Hughes, T.J., 2007. Patient-specific vascular nurbs modeling for isogeometric analysis of blood flow. Computer methods in applied mechanics and engineering 196, 2943–2959.

Zhang, Z., Marin, D., Drangova, M., Boykov, Y., 2021. Confluent vessel trees with accurate bifurcations, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9573–9582.
January 2023

**Subject: Second revision of our submission for the Medical Image Analysis journal**

Dear Reviewers, Dear Associate Editor,

Thank you for sharing your time again for this second round of revision of our article. The remarks and suggestions on our work were relevant and helpful. We have addressed those comments in a revised version of the article. The modifications are highlighted in red.

We hope that those modifications address your concerns and meet your approval standards for the Medical Image Analysis journal.

Sincerely,

*The comments of the reviewers are addressed point by point below.*

**Associate Editor:**

> *"I have read the comments raised and the revised manuscripts. I would like to thank the authors for revising the manuscript and re-submitting it. After reading, as raised, there is a need for more clarification about the performance of unstructured tetrahedral cells, the advantages of hexahedral meshes over tetrahedral meshes and related experiments, the possible effects of centerline extraction methods and post processing methods, and CFD applications. It would be better if the authors could further address the concerns."*

**Response:** Thank you for those relevant suggestions on our work. We revised the article in light of the comments made by the reviewers. The article has been modified as follows :

1) A new experiment has been conducted to demonstrate the advantages of hexahedral meshes over tetrahedral meshes. These results are consistent with the conclusion of previous work cited in the paper to justify the choice of hexahedral cells over tetrahedral cells in our framework. The results of this experiment were summarized in a new Application section (Section 7.4) of the article. The details of the methods and results of this additional experiment were judged too long to be added to the main article, and were therefore reported in detail in Supplementary materials.

2) As requested by Reviewer 1, a new section (Section 7.4) was added to the article to show a computational fluid dynamics application of our framework. In this section, we first summarize the conclusion of the comparison experiments (further detailed in the Supplementary materials). Then, we propose a new application of our framework in which a patient-specific mesh of the middle carotid artery (MCA) and downstream vessel was reconstructed by our method, with the addition of a pathology (stenosis). CFD simulations were successfully run in the healthy and pathological cases, and the resulting hemodynamic values (blood velocity) were compared. With this application, we demonstrated the potential of our meshing method to easily design and conduct blood flow studies of cerebrovascular pathologies by CFD.

3) The distinction between the different types of meshes (hexahedral, tetrahedral, structured, unstructured) was clarified in the Introduction and State-of-the-art sections, and the results of previous hex/tet comparison studies were reported more precisely.

4) We addressed the concerns of Reviewer 1 concerning the accuracy of the centerline extraction in the present letter by discussing this point in light of the challenging aspect of the reconstruction of complex vascular networks with small arteries, the current limitations of segmentation-based methods and the editing flexibility offered by our framework. The last point was emphasized by adding a new Section (5.2.3) in the article and showing editing examples in our answer to Question 4 and 5 of Reviewer 1.
**Reviewer 1:**

"I thank the authors made a great effort in addressing my concerns raised in last round of review and revising the manuscript. Many of my concerns have been well addressed. But some of them have to be further clarified. "

"1. The authors stated that, by citing several references, unstructured tetrahedral cells lead to higher computational cost and less accurate results than flow-oriented, structured hexahedral cells. This statement is still doubtful. First, what are "unstructured tetrahedral cells"? How about "structured tetrahedral cells"? Why is its performance worse than flow-oriented, structured hexahedral cells?"

**Response:**

Structured meshes are defined by an organization in which the neighbor cells can be identified without needing a connectivity matrix. This is the case with regular grids for instance. This property, which can speed up the numerical simulation, can only be achieved with hexahedral cells. By writing "unstructured tetrahedral cells", we wanted to emphasize the fact that in tetrahedral meshing, the cells do not have any particular organization and are randomly distributed in the domain. We understand that the term "unstructured tetrahedral cells" was misleading as it suggests that "structured tetrahedral cells" are also used. We replaced this term by "tetrahedral cells" in the revised version of the article. A short paragraph was also added in the Hexahedral meshing part of the State-of-the-art (Section 2.3), to clarify this notion of structured mesh :

"The commonly used cell shapes include tetrahedral, prismatic and hexahedral cells. Hexahedral meshes can be further divided into two categories; the structured meshes, where the neighborhood relationships between the cells are defined in the mesh structure (e.g. regular grid), and the unstructured mesh. In the case of blood vessels, structured and unstructured hexahedral meshing allows for the creation of flow-oriented cells. Studies of the literature show that both the shape of the cells (tetrahedral, hexahedral) and the type of mesh (structured or unstructured) influence the cost and the stability of the numerical simulation."

We addressed the second part of the remark concerning the comparison of tetrahedral and hexahedral meshes
in the answer to the next question.

"**2. As mentioned in the introduction, many FE-based CFD models employ tetrahedral meshes. To my knowledge, it is probably more widely used than hexahedral meshes. In this case, the advantages of hexahedral meshes over tetrahedral meshes, particularly in terms of CFD, should be further analyzed.**"

**Response:**

Tetrahedral meshes are indeed more widely used than hexahedral meshes, mainly because of the difficulty of the creation of a hexahedral mesh. Structured hexahedral meshing reduces the computation time, and improves the convergence of the results. Less cells are required to reach mesh independence. The advantages of flow-oriented, hexahedral meshes over tetrahedral meshes based on the existing literature were further detailed in the State-of-the-art (Section 2.3) of the revised article :

"Vinchurkar and Longest (2008), De Santis et al. (2010) and Ghaffari et al. (2017) compared the performances of hexahedral and tetrahedral/prismatic meshes in different models (airways, coronary tree and cerebral arteries) and applications. Those studies demonstrated that hexahedral meshes in general, and more specifically structured hexahedral meshes, converge better for the same accuracy of the result. De Santis et al. (2010) and Ghaffari et al. (2017) reported that 6 times less cells (resp. 10 times) and 14 times (resp. 27 times) less computational time were required. Finally, Vinchurkar and Longest (2008) insisted on the importance of having hexahedral flow-oriented cells for near wall measurement (e.g particle deposition, wall shear stress)."

The reason for this difference is thought to be two-fold :

(1) The cells of hexahedral structured mesh are aligned in the flow direction which reduces the computational error due to numerical diffusion. On the other hand, unstructured meshes with randomly oriented cells are reported to introduce relevant numerical diffusion in the solution.
(2) As the vessels are tubular, the flow does not show strong variation in the longitudinal direction, but mainly in the radial direction. Therefore, the ideal mesh for numerical simulation has a low longitudinal cell density but a high radial cell density, This can be achieved easily with hexahedral meshes, but not with tetrahedral meshes.

The results reported in the literature were confirmed within our framework by designing a set of computational fluid dynamic experiments to compare tetrahedral meshes with hexahedral meshes created with our method concerning the simulation accuracy, stability, and computational cost. This experiment is further detailed in our answer to question 3. Besides, as we mention in the State-of-the-art section (Section 2.3), the advantages of hexahedral meshes are not limited to CFD : “*The advantages of hexahedral cells are not limited to CFD; this type of mesh simplifies the boundary layer creation, bridges the gap between representation and physical simulation, and provides a basis for NURBS modeling (Zhang et al., 2007) and isogeometric analysis*”

**3. Is it possible to design a set of experiments to demonstrate the proposed hexahedral meshes can achieve "more accurate results with less computational cost" than tetrahedral meshes?**

**Response:** A new experiment has been conducted to demonstrate the advantages of hexahedral meshes over tetrahedral meshes. For this, we created a dataset of tetrahedral and hexahedral meshes with increasing cell density, as shown in the Fig.1 below (Fig.1. was added to the Supplementary materials).

<figure><img src="image_27.png" /><figcaption>Fig.1. Illustration of the different meshes produced for this study. Tetrahedral and hexahedral meshes of different densities are represented. The meshes were cut using an oblique plane to reveal the inside cells.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig.2. Graph of the sectional maximum velocity as a function of the number of cells in the mesh for both tetrahedral meshes and hexahedral meshes. The maximum velocity was averaged on three cross-sections along the tube model. The analytical value expected is shown by the black dotted line.</figcaption></figure>
CFD simulations were run in conditions that mimic blood flow in cerebral arteries and tetrahedral and hexahedral meshes were compared according to different CFD-related criteria; the meshing time, the simulation residuals, the simulation time, the number of iterations until convergence, and the number of cells required to reach the mesh independence. A simple tubular model was first used to allow the comparison to the analytical velocity profile provided by the Poiseuille equation in a pipe. The experiment results revealed that "the mesh independence was reached faster using hexahedral meshes than tetrahedral meshes, for a more accurate sectional maximum velocity value. The convergence of the simulation was also improved, as 4 times fewer iterations were necessary to obtain convergence of the results with hexahedral meshes. The simulation time was reduced on average by a factor 3, which adds to the fact that fewer cells are required to reach accurate results with hexahedral meshes, reducing the computational cost even more." These results are shown in the main article by the graph that we copied below in Fig. 2. More graphs and visual results are provided in Supplementary materials, including the Fig. 3 below.

<figure><img src="image_29.png" /><figcaption>Fig. 3. Velocity field and velocity profile for a cross-section of the tube model for coarse and fine tetrahedral (a) and hexahedral meshes (b). The analytical velocity profile is plotted alongside the simulation velocity profile. The meshes labeled "coarse" and "fine" in this were selected so that the number of cells is similar between hexahedral and tetrahedral meshes, to facilitate the comparison.</figcaption></figure>
This experiment was then reproduced in a realistic model of middle carotid artery bifurcation, showing similar results. These results are consistent with the conclusion of previous work cited in the paper to justify the choice of hexahedral cells over tetrahedral cells in our framework. The details of the methods and results of this additional experiment were judged relevant, but too long to be added to the main article. Therefore, they were summarized in a new section of the Applications (Section 7.4), and the details of the mesh creation as well as more results figures are provided in the Supplementary materials (Section 3.2).

> **4. While you started that the performance of the segmentation models may affect the performance of segmentation-based methods, I wonder if the performance of centreline extraction methods affects the performance of the proposed methods. How can you guarantee the accuracy of the mesh based on the proposed centreline-based model?**

**Response:**

Thank you for this very relevant comment. About the first question “I wonder if the performance of centreline extraction methods affects the performance of the proposed methods” Indeed, we acknowledge that the performance of centerline extraction method affects the accuracy of the mesh produced by our method, as shown in Figure 22. As a matter of fact, we acknowledged this limitation in the previous round of revision, in the Conclusion section :

> "Besides, we want to emphasize that our objective with this work was 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 some manual post-treatment may still be required before simulation. In this way, our framework offers more editing flexibility than other meshing methods. To take advantage of this flexibility, we developed a vascular network editing software, with a user-friendly interface. This interface integrates the modeling and meshing methods described in this article as well as other editing functionalities such as centerline editing, branch removal or angle modification. It opens vascular modeling and hexahedral meshing to medical doctors and non-expert users."

In conclusion, we acknowledge this limitation, and we designed our method precisely to try to overcome the limits of centerline extraction method, either manual or automatic, based on a segmentation or not. Moreover, the automatic algorithm proposed was complemented by an editing software published in a separate article, which enables to correct manually in 3D the extracted centerlines on the basis of the medical image and/or expert knowledge, as shown in Figure 4. This software was presented in an conference proceeding article:

M. Decroocq, G. Lavoué, M. Ohta and C. Frindel, "A Software to Visualize, Edit, Model and Mesh Vascular Networks," 2022 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), Glasgow, Scotland, United Kingdom, 2022, pp. 2208-2214, doi: 10.1109/EMBC48229.2022.9871365.

https://ieeexplore.ieee.org/document/9871365

This article does not include any of the methodological aspects explained in the article that we wish to publish in
Medical Image Analysis Journal, and focuses only on the software functionalities.
<figure><img src="image_30.png" /><figcaption>Fig. 4. Software use case for the manual correction of a vascular network of the BraVa database on the basis of the medical image.</figcaption></figure>

Concerning the second question “How can you guarantee the accuracy of the mesh based on the proposed centreline-based model?”, there are indeed some clear limitations in the evaluation of the accuracy of the reconstructed mesh, as the vessels reconstructed are very small and often hardly visible on the medical image. Moreover, we demonstrated in the experiments of Section 6.2 that even the expert segmentations, which are generally the golden-standard in this type of studies, resulted in geometrical and topological inaccuracies in the final reconstructed mesh. However, we would like to put emphasis on two points concerning the accuracy of the proposed model :

(1) From a topological point of view, we showed in the experiments of Section 6.2 that the mesh reconstructed from centerlines are closer in to the real anatomy than the commonly used segmentation-based methods (no cycles, no high order branchings). The topological mistakes can be manually corrected effortlessly with this representation type (and the software we developed).

(2) From a geometrical point of view, one of the clear limitations of the proposed model is the tubularity assumption. However, an observation of the medical images (MRA) and the segmentation-based meshes led us to believe that the accuracy of the non-circular shape of vessel cross sections is not guaranteed either by segmentation-based methods. We could report for instance missing vessel parts or bumps in the vessel wall. In this type of case, we believe that our model, even though it relies on the tubularity assumption, can produce more realistic vessel shapes. Besides, as shown in Figure 16, the trajectory of the vessels in our model matches
the trajectory of the vessels in segmentation-based meshes, which leads us to think that the centerline-based method proposed can provide accurate meshes.

In conclusion, we can not give strong guarantees concerning the accuracy of the mesh reconstructed from centerlines. However, given the challenging aspect of the meshing of small cerebral arteries, our method enables to reconstruct larger cerebral vascular networks than what is feasible using segmentation. The time-consuming (and user-dependent) manual correction of segmentation-based meshes is also facilitated, which constitute a step towards the reconstruction of the whole brain arterial network. Besides, many applications do not require a perfectly accurate patient-specific mesh, but only a realistic surface.

> **5. As mentioned, segmentation-based approaches may result in topologic and geometric inaccuracies in the final mesh (bumps, disconnected or merging vessels, cycles). Is it possible to use some post-processing methods to address these shortcomings? For example, could we employ some interactive segmentation techniques to incorporate human interaction to enhance the segmentation accuracy.**

**Response:** The main limitation of segmentation based methods lies in the low resolution of medical images. As demonstrated in the experiments of Section 6.2, even the expert segmentation resulted in geometrical and topological inaccuracies in the final reconstructed mesh. This suggests that the problem is not only related to the performance of the segmentation method but also to the idea of segmentation-based meshes itself. It is this observation which led us to step away from the constraints related to segmentation by basing our method on centerlines. As for the post-processing of segmentation results or segmentation-based meshes, some methods were already proposed in this way in the literature (e.g. the tools provided in the 3D slicer software).

However, we would like to argue that such an interactive post-processing method still requires to work pixel-by-pixel in the case of segmented images or node-by-node in the case of meshing. It remains more burdensome and time consuming than the modification of centerlines, which are a lighter representation of vascular networks. Time can be a considerable obstacle when we want to reconstruct networks as large and complex as the cerebral vascular network.

<figure><img src="image_31.png" /><figcaption>Fig. 5. Segmentation-based mesh produced by Unet and centerline-based mesh obtain with our method.</figcaption></figure>

If we take the example proposed in Figure 5, two vessels are merged in the segmentation-based mesh. There are two options to correct it; cut the edges of the mesh between both vessels and reconstruct the wall to fill the hole created, or modify the segmentation, by removing the pixels in between the vessels, which may also alter the shape of both vessels. In the case of centerline-based mesh (on the left), this merging problem might also result in errors on the centerlines (as explained in the previous question). On the case shown in Figure 2 there was an inversion in the trajectory of the vessels in the extracted centerline. Nevertheless, correcting this problem is straightforward with our framework. By modifying the connection between two data points of the centerline and moving a few other data points, the vessels can be easily separated and the original topology recovered. This editing process was made easier by the editing software that we developed to supplement our method, as mentioned in the answer to question 4. The editing applications of our framework were emphasized in the revised manuscript by some modifications in Section 7.2 "Topology and geometry editing" and by adding a new subsection 5.2.3. "Data encoding".
"6. The authors employ the proposed techniques to several applications to demonstrate their effectiveness, which is good. As the main objective of this method is to improve the performance of CFD models. In this case, could we apply the proposed techniques to some CFD applications to demonstrate their merits compared with existing meshing methods? "

**Response:** As suggested, a new section (Section 7.4) was added to the article to show a computational fluid dynamics application of our framework. In this section, we first summarize the conclusion of the comparison experiments described in our answer to question 3. This demonstrates the merits of our hexahedral meshing method compared to other existing methods. Then, we propose a new application of our framework in which we compare a healthy and pathological case in a model of the middle carotid artery (MCA) and downstream vessels. The vascular network (with several branches and bifurcations) was reconstructed by our method, in a case where the segmentation-based mesh produced by Unet failed to provide a mesh usable for CFD. Flow extensions were automatically added to the inlet and outlets, and a 50% stenosis was automatically added to the MCA using our framework, as shown in Figure 6 below (this Figure is provided in the Supplementary materials).

<figure><img src="image_32.png" /><figcaption>Fig. 6. a) Segmentation-based mesh produced by Unet for the arterial network of interest, i.e the right middle cerebral artery (MCA) territory. b) Mesh produced for the same arterial territory as (a). This mesh was reconstructed with our method from the centerlines of the database BraVa. c) Mesh after automatically adding extensions to the inlet and outlets. d) Mesh after automatically adding a 50% stenosis in the MCA. In c) and d), the hexahedral mesh created in the stenosis region is highlighted.</figcaption></figure>

CFD simulations were successfully run in the healthy and pathological cases, and the resulting hemodynamic values (blood velocity) compared in Figure 7 below (see Section 7.4 of the main article). With this application, we demonstrated the potential of our meshing method to easily design and conduct blood flow studies of cerebrovascular pathologies by CFD, for arteries that could not be meshed successfully by segmentation.
<figure><img src="image_33.png" /><figcaption>Fig.7. CFD simulation results for the mesh without stenosis (a) and with stenosis (b). In both cases, the velocity streamlines were rendered and the velocity field is shown by a cut on the stenosis region.</figcaption></figure>

"7. There are still a few typos and grammatic errors and please make more effort in proofreading."

Response:

Thank you for this remark, there were indeed typos and grammatical errors that we corrected in the final version.

Reviewer 2:

"In the reviewed version, the authors have addressed concerns raised in reviews. A more comprehensive literature review is performed. Experiments on one more database are added. Besides, deep learning based methods are also considered in the experiments."

"Please double check these typos:
the 3rd line under figure 5, 'mesh. The nodes of C0 are connected to the nodes of the half-'.'C0'.
($C_{0}$);
the superscript 'm' in the equation (2) seems to be 'm-1'. (You only define points from 0 to m-1);
In the figure 14, 'quantitative comparison' (quantitative)."

Response: We thank the reviewer for these comments and the time spent reviewing our article. The typos mentioned have been corrected as suggested.
<figure><img src="image_34.png" /></figure>

# Modeling and hexahedral meshing of cerebral arterial networks from centerlines

Méghane Decroocq<sup>a,b,d,e,f</sup>, Carole Frindel<sup>a,e,*</sup>, Pierre Rougé<sup>d</sup>, Makoto Ohta<sup>d,e</sup>, Guillaume Lavoué<sup>b,c</sup>

## ARTICLE INFO

Article history:

2000 MSC: 41A05, 41A10, 65D05, 65D17

Keywords: Cerebral arterial network, Centerlines, Hexahedral mesh, Computational fluid dynamics

## 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 a 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 reproducibility.

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 to the common defects of vascular centerlines and an increased mesh quality compared to other state-of-the-art methods. As it combines both modeling and meshing techniques, it can be applied to edit the vascular models effortlessly to study the impact of vascular geometry and topology 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 were meshed without defects needing manual intervention, despite the challenging aspect of the input data. The source code will be released publicly.

© 2023 Elsevier B. V. All rights reserved.

## 1. Introduction

Cerebrovascular diseases, such as stroke, can cause severe disability or death (Ramos-Lima et al., 2018). 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, as it provides information on the hemodynamics from the vessel geometry (Saqr et al., 2020; Sugiyama et al., 2016). **The main limitation in the use of CFD is the creation of the computational mesh.** Indeed, numerical simulation requires a smooth, anatomically realistic mesh of the arterial wall to provide reliable results. **In pathologies like is-

<footnote><sup>a</sup>CREATIS, Université Lyon1, CNRS UMR5220, INSERM U1206, INSA-Lyon, 69621 Villeurbanne, France</footnote>

<footnote><sup>b</sup>LIRIS, CNRS UMR 5205, F-69621, France</footnote>

<footnote><sup>c</sup>Ecole Centrale de Lyon, France</footnote>

<footnote><sup>d</sup>ELyTMaX IRL3757, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard Lyon 1, Tohoku University, 980-8577, Sendai, Japan</footnote>

<footnote><sup>e</sup>Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan</footnote>

<footnote><sup>f</sup>Graduate School of Biomedical Engineering, Tohoku University, 6-6 Aramaki-aza-aoba, Aoba-ku, Sendai, Miyagi 980-8579, Japan</footnote>

<footnote>*Corresponding author: Carole Frindel
e-mail: carole.frindel@creatis.insa-lyon.fr (Carole Frindel)</footnote>
chemic stroke, the distribution of the vessels in the different vascular territories of the brain impacts the position and evolution of the lesion (Hodneland et al., 2019). It requires reconstructing large and complex cerebral arterial networks with small vessels whose radius is close to the image resolution, which remains very challenging. Besides, in the finite element method, the shape of the cells inside the volume also affects the simulations. In particular, the flow-oriented, structured hexahedral cells were shown to improve the stability of the simulation while lowering the computational cost (Vinchurkar and Longest, 2008; De Santis et al., 2010; Ghaffari et al., 2017). Those results were confirmed in this work by running a CFD experiment whose results are presented in Section 7. In spite of this, tetrahedral cells remain widely used due to their ability to automatically mesh any complex shape. The approaches investigated in the literature to address the meshing of vascular networks can be divided into two categories: the segmentation-based and the 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 to develop efficient vessel-enhancing filters (Jerman et al., 2015; Merveille et al., 2017) and to improve the segmentation methods. In particular, the rise of deep learning-based segmentation methods resulted in significant progress in vascular segmentation (Tetteh et al., 2020; Livne et al., 2019). However, the accuracy of the segmentation does not guarantee the accuracy of the mesh it entails (e.g vessels merging due to the image resolution, disconnected vessels, bumps), nor its usability for numerical simulation. Besides, the vascular network 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 other, 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 (Zhang et al., 2021; He et al., 2020). 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 for the construction of big databases of large vascular networks (Wright et al., 2013) or the creation of ideal models. It also offers more editing flexibility than segmentation-based meshes. As it encodes the vessel topology and orientation, it has a high potential for the creation of meshes with high-quality, flow-oriented 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 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 article, we propose a method to overcome those limitations and create a high-quality, structured hexahedral mesh for CFD from centerlines only, opening the way to CFD in large

cerebral arterial networks.

## 2. Related work

### 2.1. *Segmentation-based meshing*

Segmentation of medical images is the most common method to obtain patient-specific meshes for CFD. In recent years, deep learning-based models have led to significant advances in vascular segmentation. More specifically, convolutional neural networks (CNNs) have achieved very good performances (Jiang et al., 2018; Tetteh et al., 2020). The popular U-net architecture (Ronneberger et al., 2015) has been successfully applied to the segmentation of intracranial vessels in Quon et al. (2020) and Livne et al. (2019). Hilbert et al. (2020) proposed an extended U-net architecture using context aggregation and deep supervision for brain vessel segmentation. Besides, the attention mechanisms have been used to help the network to better learn global dependencies and increase the receptive field in Mou et al. (2021), Ni et al. (2020) and Li et al. (2021).

For medical applications such as CFD, more than the segmentation itself, the smoothness and the topological accuracy of the mesh it entails are critical. However, in the literature, there was very little focus on the conversion of the segmented volumes to mesh. Recently, Wickramasinghe et al. (2020) and Kong et al. (2021) introduced new neural network architectures to reconstruct 3D meshes directly from 3D image volumes. Despite those recent advances, the meshing largely relies on algorithms such as the marching cubes, followed by a smoothing step, to produce a surface mesh with tetrahedral elements (Watanabe et al., 2018; Misaki et al., 2021). However, this type of segmentation-based meshes commonly suffers from topological inaccuracies (e.g. merging or disconnected vessels, bulges, missing vessels) and requires a burdensome manual post-processing (Glaßer et al., 2015), as we demonstrated in Section 6.2. Such problems are not correctly captured by the image-based metrics (e.g DICE score) used to evaluate the segmentation methods. To overcome those challenges, the centerline representation of vascular networks has recently gained interest.

Some recent segmentation approaches propose to integrate the vessel centerline information to build more topology-oriented metrics. (Keshwani et al., 2020) proposed to segment the vascular network from its skeleton by learning a connectivity metric between center-voxels. Besides, (Shit et al., 2021) introduced a novel topology-preserving loss for the training of neural networks, which relies on the centerlines of the predicted segmentation. The information provided by the centerlines allowed the neural network to improve the topology correctness of the segmentations. In this context, we believe that the use of centerlines in the meshing process can offer many advantages for CFD applications.

### 2.2. *Centerline-based meshing*

In this part, 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 - and 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 into explicit methods, where a tetrahedral mesh of the surface is produced, and implicit methods where the surface is represented by implicit functions. Implicit methods employ radial basis functions (Hong et al., 2018), implicit extrusion surfaces (Hong et al., 2020) or local implicit modeling (Kerrien et al., 2017) to reconstruct vascular networks from medical images or from centerlines (Abdellah et al., 2020). If they stand out by their ability to reconstruct complex branching topology, they do not allow as much control on the final mesh as explicit methods (e.g. hexahedral meshing) which makes them less suited for CFD using the finite element method.

In explicit methods, the vessel surface is obtained by sweeping along the centerlines. The quality of the reconstruction depends on the way centerline points are approximated, usually with Bezier or spline functions (Guo et al., 2013; Kociński et al., 2016; Ghaffari et al., 2017), and might be affected by the quality of the input centerline. The details of the approximation methods and the evaluation of the accuracy of the estimation of coordinates and derivatives were not provided in previous studies. For the branching part, various bifurcation models were proposed. In the work of Kociński et al. (2016) and Ghaffari et al. (2015), 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 Kociński et al. (2016) and Bezier segments for Ghaffari et al. (2015). This geometric model facilitates the creation of hexahedral meshing. However, the realism of the bifurcation shape 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.

Han et al. (2015) and Guo et al. (2013) modeled bifurcations using three tubes connecting the inlet and outlet sections. Half of each section of the tube is meshed by sweeping and the hole left in the middle is filled afterward. Their method guarantees the smoothness of the model but results in unnatural-looking bifurcations. Finally, Zakaria et al. (2008) proposed a physiologic model where the bifurcations are represented by two merged tubes. It was validated with regard to both the accuracy of the anatomy and the CFD simulations. It does not rely on the geometrical center of the bifurcation, but on a set of physiolog-ical parameters (apex, apical sections, inlet and outlet sections), and shows more anatomical realism. However, the authors extracted the model parameters from a surface mesh, and they did not suggest a way to extract them from centerlines.

## 2.3. Hexahedral meshing

For applications in CFD with the finite element method, the inside of the surface mesh must be discretized into cells. The commonly used cell shapes include tetrahedral, prismatic and hexahedral. Hexahedral meshes can be further divided into two categories; the structured meshes, where the neighborhood relationships between the cells are defined in the mesh structure

(e.g. regular grid), and the unstructured mesh. In the case of blood vessels, structured and unstructured hexahedral meshing also allows for the creation of flow-oriented cells. Studies of the literature show that both the shape of the cells (tetrahedral, hexahedral) and the type of mesh (structured or unstructured) influence the cost and the stability of the numerical simulation.

Vinchurkar and Longest (2008), De Santis et al. (2010) and Ghaffari et al. (2017) compared the performances of hexahedral and tetrahedral meshes for different models (airways, coronary tree and cerebral arteries) and applications. Those studies demonstrated that hexahedral meshes in general, and more specifically structured hexahedral meshes, converge better for the same accuracy of the result. De Santis et al. (2010) and Ghaffari et al. (2017) reported that 6 times less cells (resp. 10 times) and 14 times (resp. 27 times) less computational time were required. Finally, Vinchurkar and Longest (2008) insisted on the importance of having hexahedral flow-oriented cells for near-wall measurements (e.g particle deposition, wall shear stress). The advantages of hexahedral cells are not limited to CFD; this type of mesh simplifies the boundary layer creation, bridges the gap between representation and physical simulation, and provides a basis for NURBS modeling (Zhang et al., 2007) and isogeometric analysis. Hexahedral meshing, and more specifically structured hexahedral meshing, is however limited by a far more complex generation process than standard tetrahedral meshes (Vinchurkar and Longest, 2008).

In the application to the arterial networks, the main challenge is the generation of the mesh at the bifurcations part. 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 methods were proposed to obtain a robust branch decomposition. De Santis et al. introduced semi-automatic methods, ranging from the manual selection of the most relevant slices of the input surface mesh (De Santis et al. (2010)), user-defined bifurcation coordinate system (De Santis et al., 2011a), to the generation and adjustment of a block-structure representation of the network (De Santis et al., 2011b). Automatic methods are based on Voronoi diagram (Antiga et al., 2002), resolution of the Laplace's equation (Verma et al., 2005), random-walk algorithm (Xiong et al., 2013) or branching templates (Zhang et al., 2007) or parametric models (Ghaffari et al., 2017). The hexahedral meshing can then be created from the decomposition through various techniques; Copper scheme in the work of Antiga et al. (2002), template grid sweeping for Verma et al. (2005), Zhang et al. (2007) and Ghaffari et al. (2017), Bezier spline modeling followed by an iso-parametric transformation of a template mesh (De Santis et al., 2011a), projection and refinement of block-structures (De Santis et al., 2011b), Laplacian-based harmonic functions combined with Catmull-Clark subdivision (Xiong et al., 2013).

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 De Santis et al. (2011a) provided their code through the user-friendly inter-
face PyFormex, which enables to generate hexahedral meshes semi-automatically from a single bifurcation vascular geometry. Only Ghaffari et al. (2017) used 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 or the meshing. They can not be applied directly from realistic centerlines extracted from medical images.

## 2.4. Contributions

In this work, our purpose is to meet the challenges arising from this state-of-the-art with a framework integrating a modeling and a meshing step. The shortcomings of segmentation-based meshing are addressed by developing a method based on centerlines. A parametric model is used to overcome the common defect of centerlines and reconstruct a realistic vessel surface. Finally, an original meshing algorithm is proposed to create high-quality hexahedral meshes suitable for CFD simulations, as opposed to the largely used tetrahedral meshes. Our main contributions in this work are:

* We introduced an original vessel model and its approximation algorithm based on penalized splines, which enable to model both the spatial coordinates and the radius in a single function and offer a good robustness to noise and low-sampling.

* The vessel model is combined with a physiological parametric model of bifurcation proposed by Zakaria et al. (2008) to form a light parametric model of the entire vascular network. A method to extract the parameters of the bifurcation model directly from centerlines is proposed, and the model was generalized to planar n-furcations. If the bifurcation model itself is not new, the use of this type of physiologic bifurcation model - as opposed to geometric bifurcation models - to reconstruct a realistic vascular shape has not been investigated in previous studies.

* We developed a parametric method to create a structured hexahedral volume mesh with flow-oriented cells from the vessel and bifurcation models proposed. It includes relaxation and smoothing steps to improve the quality of the cells without deforming the model shape. This method gives more control over the distribution and density of the cells than the commonly used tetrahedral meshing.

* The model and the mesh are stored in a graph structure which enables to easily and inexpensively edit the topology and geometry of the vascular networks.

The proposed framework is fast and automatic. As it is based on centerlines only, it opens the way to numerical simulation in large cerebral vascular networks. It was evaluated qualitatively and quantitatively against other explicit and implicit centerline-based meshing methods, as well as segmentation-based meshing methods. Finally, several practical applications are presented, including the meshing of a large database of 60 large cerebral networks, pathology modeling, topology and geometry editing, and finally a CFD study comparing a healthy and stenotic middle carotid artery.

## 3. Input data

The input vessel centerlines we consider are composed of a 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 called n-furcation. The centerlines are stored using the swc format or VMTK format of Izzo et al. (2018). In this work, we used two publicly available datasets. The Aneurisk database (Aneurisk-Team, 2012) provides 3D models of the main arteries of the circle of Willis for patients with an aneurism. High-resolution centerlines were extracted from the surface meshes using the VMTK software. The BraVa database (Wright et al., 2013) 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 the ImageJ plugin Neurite Tracer (Longair et al., 2011) 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 centerline were proposed in the literature, based on the approximation of data points by Bezier segments (Ghaffari et al., 2017), regression splines (Kociński et al., 2016), free knot regression splines or local polynomial smoothing (Sangalli et al., 2009b). Only Sangalli et al. (2009b) gives the detail of the implementation of the approximation method and provides a thorough study of the accuracy of their model regarding the spatial coordinates and the derivatives. However, the accuracy of both the first and second derivatives is crucial, as the vessel curvature impacts the hemodynamics (Sangalli et al., 2009a). Moreover, the meshing techniques are often based on the normals of the centerline (Kociński et al., 2016; Ghaffari et al., 2017). It is important to note that the proposed approximation methods (Sangalli et al., 2009b; Kociński et al., 2016; Ghaffari et al., 2017) 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 combining 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, \quad (1)$$

for $u \in [0, 1]$, where $N_{i,p}$ is the ith basis spline function of order p and {$P_0, P_1, ..., Pn-1$} the n control points of the spline.
The shape of the basis splines functions and therefore the part of the spline controlled by each 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 Sangalli et al. (2009b). In the other approach, a relatively large number of control points and a uniform knot vector are used and the smoothness is constrained by a penalty on the second derivatives (Craven and Wahba, 1978; Eilers and Marx, 1996). 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 Eilers and Marx (1996). 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, \dots, P_{n-1}) = \sum_{k=0}^{m-1} |D_k - s(t_k)|^2 + \lambda \sum_{j=2}^{n} (P_j - 2P_{j-1} + P_{j-2})^2, \quad (2)$$

where $t$ is a time parametrization vector that 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 scales that might show different noise levels. For this reason, we propose to approximate them 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)}, \quad (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$ detailed in Section 6.1. A comparison study with other optimization criteria for the smoothing parameter $\lambda$, such as the Bayesian information criterion or cross-validation, detailed in supplementary materials, led to the choice of this criterion.

The linear system is then solved for the data $(t, r)$ where $t$ is the time parametrization of each data point and $r$ their radius value:

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

The value of $\lambda_r$ is also selected by minimizing the Akaike criterion 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><img src="image_35.png" /><figcaption>Fig. 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 as red dots with radius represented by black circles. In (d), the vessel surface defined by the approximated spline is represented.</figcaption></figure>

## 4.2. Bifurcations

In this part, we focus on the modeling of bifurcations.

### 4.2.1. Zakaria's model

Zakaria et al. (2008) proposed a parametric model for non-planar bifurcations. Their model was validated regarding 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 the reconstruction of bifurcations from sparse data. In this model, bifurcations are created by merging two tubes that represent 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 merge. 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. Each circular cross-section C is represented by the three spatial coordinates of its center $P_c$, the radius $r_c$ and the normal vector $\vec{n}_c$. The centerline of each tube is defined by a spline function $spl_{1,2}$. The first segment of the centerline connects the inlet section $C_0$ to the apical section, and the second connects the apical section to the outlet sections. The tangent of the centerline segments matches the normal of the cross-sections that they connect. 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><img src="image_36.png" /><figcaption>Fig. 2. The five cross-sections bifurcation model introduced by Zakaria et al. (2008).</figcaption></figure>

### 4.2.2. Parameter estimation

In order to apply this bifurcation model to our framework, we introduce an algorithm to estimate the parameters of the bifurcations directly from the input centerline data. For all the bifurcations in the centerline data, the inlet data points (in light blue in Figure 3 (a)) are concatenated with each of the outlet data points (resp. in deep blue and green in Figure 3 (a)) to form two individual vessel 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 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$ and $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.

### 4.2.3. Tangent continuity

The full vascular network model is created by assembling the vessels and bifurcations models presented in the previous sections. In order to preserve the continuity of the different parts of the network, the inlet and outlet tangents of the vessels must

<figure><img src="image_37.png" /><figcaption>Fig. 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.</figcaption></figure>

match the normal of the inlet and outlet cross-sections of the bifurcations they connect. For this, we introduce an additional constraint on the endpoints and tangents 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 Piegl and Tiller (2000). In this work, to limit the influence of the end constraints on the approximation of data points, we propose a weaker constraint that fixes the end tangent while the derivative is free. We consider a spline s as defined by Equation 1. Because we work with clamped curves, $s(0) = P_0$ and $s(1) = P_{n-1}$. Moreover, $s'(0)$ (respectively $s'(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:

$$ \left\{ \begin{array}{l} 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{array} \right. \qquad (5) $$

where $\alpha$ and $\beta$, the end tangent magnitude, are additional parameters to optimize. Those constraints guarantee the $G_1$ continuity of the final network model. The details of the system resolution are given in the supplementary materials.

## 5. Structured hexahedral meshing

In this section, we present the meshing algorithm developed to produce a hexahedral mesh with flow-oriented cells from the parametric model described in Section 4.

### 5.1. Bifurcations

#### 5.1.1. Decomposition

The bifurcations are the most challenging parts to mesh with hexahedral elements. The meshing approach proposed relies on a decomposition scheme to split the bifurcation into three geometrical branches; one inlet branch and two outlet branches. This method allows to have the meshing advantages offered by the geometrical bifurcation models of other methods of the state-of-the-art, while keeping the anatomical realism of the physiological bifurcation model used in this work. Figure 4 (b) gives an example of branch splitting using three separation
planes. Antiga et al. Antiga and Steinman (2004) proposed a bifurcation decomposition scheme based on the Voronoi diagram of a surface mesh. The proposed scheme is robust to variations in input geometry and has been successfully used for hexahedral meshing (Antiga and Steinman, 2004). 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. Based on this work, we introduced a decomposition scheme that relies on the spline and bifurcation models described above.

In this decomposition method, three separation planes are defined by a set of five points; the apex point AP, which is 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_{m1}$ and $p_{m2}$, where $p_{m1,2}$ are the projection of the key points $m1,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><img src="image_38.png" /><figcaption>Fig. 4. Geometric decomposition of the bifurcation model. In (b), the end cross-sections are represented in red, and the separation planes in black.</figcaption></figure>

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 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 method enables us to handle large radius differences between the daughter vessels, as 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, the surface mesh of the bifurcation is created using the separation planes defined in the previous section. First, we create an initial mesh grid that connects the end cross-sections to the separation planes with a set of successive sections, as illustrated in Figure 6. Each section of the mesh has a number *N* of nodes where *N* can be any multiple of 4. Figure 5 illustrates the initial mesh creation process. 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 connecting the separation point AP, $SP_1$ and $SP_2$ to both center points $CT_0$ and $CT_1$.

<figure><img src="image_39.png" /><figcaption>Fig. 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}$.</figcaption></figure>

Once the nodes of the end sections and the nodes of the separation half-sections are computed, they are connected to form a surface mesh. The nodes of $C_0$ 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 that is used to control the topology and geometry of the final mesh grid, but it does not necessarily lie on the exact surface of the bifurcation at this stage. 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$ do not constitute a good approximation of the centerline of the three geometric branches. 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.

The properties of the resulting mesh depend 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
<figure><img src="image_40.png" /><figcaption>**Fig. 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.</figcaption></figure>

so that the normal of the end sections is preserved in the out-
put surface mesh. If 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, causing less displacement of
the nodes during projection. Moreover, the preservation of the
normal of the end sections facilitates the inclusion of the bifur-
cation mesh in larger arterial networks; 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 bifurcation meshing method results
in an uneven sampling of the nodes along the trajectories that
can produce faces with heterogeneous size or important skew-
ness. Moreover, a rupture of continuity is observed when the
curves cross the separation between two branches. Those un-
wanted features are corrected by the relaxation of the nodes of
the surface mesh. Mesh smoothing methods are an easy way to
reduce the skewness of faces but it triggers important deforma-
tions of the general shape of the model. To avoid deformations,
Vidal et al. (2015) 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. The nodes are then projected back to the
original surface. To prevent cross-sections from intersecting,
the projection is made in the direction of the vector connecting
the center of the cross-section to the node to project. This pro-
cess 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, mea-
sured by the scaled Jacobian. We observe that while the shape
of the model is preserved, the quality of the faces near the sep-
aration planes is improved, and the grid now smoothly crosses
the separation planes. Based on the average quality of the faces,

we estimated that 5 relaxation iterations give the best results.

<figure><img src="image_41.png" /><figcaption>Fig. 7. Original bifurcation mesh and mesh after 1 and 5 relaxation itera-
tions. The scaled Jacobian is used to measure the quality of the cells be-
tween -1 (poor quality) and 1 (high quality).</figcaption></figure>

**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 an-
gle where the two vessels merge. The curvature in the apex
regions impacts the pressure and velocity fields obtained by nu-
merical simulation, as shown by Haljasmaa et al. (2001). The
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 strug-
gle to generate important local deformations. Zakaria et al.
(2008) proposed to smooth the apex region by projecting the
nodes on a sphere of a 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 our
hexahedral surface mesh, we propose a method to reduce this
complex 3-dimensional problem to a 2-dimensional problem.

<figure><img src="image_42.png" /><figcaption>Fig. 8. Illustration of the apex smoothing pipeline.</figcaption></figure>

Figure 8 illustrates the proposed smoothing method. The curves connecting two nodes of the end sections of the bifurcation are extracted from the 3D mesh (e.g, the curve in red in 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 corresponds to the desired apex radius of curvature is rolled along those 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 the 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><img src="image_43.png" /><figcaption>**Fig. 9.** Apex smoothing with different radius of curvature values. The columnar map encodes the local distance to the original mesh, on the left.</figcaption></figure>

The described smoothing method enables to control the direction of projection and the sampling of the projected nodes, preserving the quality of the faces. As shown in Figure 9, the smoothing is very local and does not affect the shape of the bifurcation 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 Zakaria et al. (2008) 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><img src="image_44.png" /><figcaption>**Fig. 10.** (a) and (b) respectively illustrate the parametric model and the branch decomposition scheme for a trifurcation. (c) shows an example of trifurcation mesh.</figcaption></figure>

## 5.2. Vessels

For the vessels, we adapted the meshing method proposed by Ghaffari et al. (2017) to the parametric model proposed in section 4.1 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 cross-sections per mm) is proportional to the mean radius of the vessel, with a proportional coefficient $d$ which can be set by the user. From each center position, $N$ nodes are radially projected on the model surface to form a circular cross-section. The projection vector is swept along the centerline so that the twisting between sections is minimized. The successive sections are finally connected 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 Ghaffari et al. (2017). From each cross-section of the surface mesh, a structured O-grid pattern is created. This pattern is composed of three different areas; the boundary layers, the intermediary layers, and the 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 handled by combining $n+1$ halves grids. The successive O-grid patterns are connected to form the hexahedral cells of the volume mesh, as shown in Figure 11.

<figure><img src="image_45.png" /><figcaption>**Fig. 11.** Illustration of the O-grid pattern and volume meshing method.</figcaption></figure>

### 5.2.3. Data encoding

In the proposed framework, the results of the modeling and meshing algorithms are encoded in graph structures. Four graphs are created, corresponding to the different steps of the process; data, topology, model, and mesh. The geometric and topological information (e.g centerline data points, model parameters, mesh nodes) are stored in the nodes and edges of the graph. This storage method facilitates data manipulation and editing as it allows the use of graph theory-based algorithms such as depth-first search or neighbor identification. The advantages of this data structure are further demonstrated in Section 7.2.
# 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 are assessed in a comparative study performed on a synthetic dataset of distorted centerlines. Then, our meshing pipeline is compared quantitatively and qualitatively with two concurrent state-of-the-art methods: deep learning-based segmentation and implicit meshing. Finally, we provide additional performance indicators 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 it does not include pathologies but goes 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 high-quality centerline data points obtained were then approximated by a 4-coordinates spline *s* that constitutes the ground truth. The control points were manually added and the accuracy of the fitting of the spatial coordinates, radius and first derivatives was 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, the ground truth centerlines were distorted to mimic defects commonly observed in realistic centerlines; 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}$. The direction of the displacement is normal to the ground truth spline *s* so that unwanted radius noise is not created during the operation. Random radius noise is generated from a zero-centered Gaussian distribution of standard deviation $\sigma_{radius}$ and added to the ground truth radius. In both cases, the applied standard deviation value is proportional to the point radius, as indicated in Table 1, in order to keep similar levels of noise between big and small vessels. Finally, low sampling is obtained by removing data points along the centerline to reach a target density.

<table><caption>Table 1. Parameters used for the distortion of the ground truth centerlines</caption><thead><tr><th>density (mm<sup>-1</sup>)</th><th>2</th><th>4</th><th>10</th><th>16</th><th>20</th></tr></thead><tbody><tr><th>σ<sub>radius</sub>(mm)</th><td>0.01r</td><td>0.05r</td><td>0.1r</td><td>0.3r</td><td>0.5r</td></tr><tr><th>σ<sub>spatial</sub>(mm)</th><td>0.01r</td><td>0.05r</td><td>0.1r</td><td>0.3r</td><td>0.5r</td></tr></tbody></table>

For each density value in Table 1, ten combinations 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 null. Each noise combination is repeated three times to account for the 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 distorted dataset to 600.

### 6.1.2. Approximation methods

To demonstrate the robustness and the accuracy of the approximation strategy used to reconstruct the surface of the vessels presented in Section 4.1, we compared it to other explicit centerline-based meshing methods of the literature (Kociński et al., 2016; Ghaffari et al., 2017). Those methods also rely on splines or Bezier curves to approximate the centerline data points and reconstruct the vessel surface. As most of the authors did not provide the details of the fitting method employed, we implemented four commonly used spline-based approximation methods with incremental complexity 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. We call it global because the 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 (Akaike (1973)) $AIC_1$:

$$AIC_1 = m \times \log(\text{SSE}) + 8(n + p). \quad (6)$$

where *m* is the number of data points, *p* is the degree of the spline, *n* is 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 Eilers and Marx (1996). It uses the same global approach as in GNP, but with a smoothing penalty controlled by a parameter $\lambda \neq 0$ as in Equation 2.

* **Spatial coordinates and Radius Penalized with Akaike criterion (SRP-AIC):** It corresponds to the approximation strategy that we propose in this work. Spatial and radius dimensions are approximated separately with two smoothing parameters $\lambda_s$ and $\lambda_r$. The comparison of our strategy with GP-AIC allows us 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 $λ$ values is another formulation of the Akaike information criterion ($AIC_2$), adapted to penalized splines, as proposed by Eilers and Marx (1996):

$$AIC_2 = m \times \log(SSE/m) + 2tr, \quad (7)$$

where $tr$ is the trace of the matrix $H = N(N'N + λΔ)^{-1}N'$.

Unlike $AIC_1$, it is not employed to choose an optimal number of control points but to select the optimal value for the smoothing parameter $λ$. This criterion was compared to other criteria of the literature; the corrected Akaike information criterion of Hurvich et al. (1998), the Schwarz's Bayesian criterion of Schwarz (1978), the cross-validation criterion and generalized cross-validation criterion of Craven and Wahba (1978). According to this study, given in supplementary material (section 1.3), the Akaike criterion gave 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><img src="image_46.png" /><figcaption>Fig. 12. Matching time parameters by minimum distance projection from $s$ onto $\hat{s}$</figcaption></figure>

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 $RMSE_{radius}$ (respectively $RMSE_{spatial}$) the root mean squared error of the radius (respectively the spatial coordinates). 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 $RMSE_{der_{spatial}}$ and $RMSE_{der_{radius}}$. As curvature is commonly considered in hemodynamic studies, the model performance concerning the centerline curvature is also measured, by the metric $RMSE_{curv}$. Finally, the length of the vessel affects the delay of blood arrival between the inlet and the outlet of the vascular tree in numerical simulations. Therefore, the difference $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 a very high spatial error. The overfitting and noise problems are partially solved by optimizing the number of control points with the 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><img src="image_47.png" /><figcaption>Fig. 13. Mesh resulting from the approximation of distorted data ($density = 1mm^{-1}$, $σ_{radius} = 0.1$) by three of the methods compared in section 6.1.</figcaption></figure>

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 was too smoothed. On the other hand, both the radius and trajectory of the vessel produced with SRP-AIC are closer to the ground truth. In conclusion, the proposed approximation method shows 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.

## 6.2. Comparison with state-of-the-art methods

In this section, meshes obtained with our method are visually and quantitatively compared to meshes produced by state-of-the-art deep learning-based segmentation methods (Tetteh et al., 2020; Livne et al., 2019), as well as a recent implicit centerline-based meshing method (Abdellah et al., 2020).
<table><caption>**Table 2.** Overall evaluation of the approximation methods: mean values of the quality criteria for all the centerlines distorted by radius noise addition. The cells in gray correspond to the lowest error for each metric.</caption><thead><tr><th></th><th>GNP</th><th>GNP-AIC</th><th>GP-AIC</th><th>SRP-AIC</th></tr></thead><tbody><tr><td>RMSE<sub>spatial</sub></td><td>8.462</td><td>0.034</td><td>0.053</td><td>0.029</td></tr><tr><td>RMSE<sub>radius</sub></td><td>17.523</td><td>0.095</td><td>0.042</td><td>0.043</td></tr><tr><td>RMSE<sub>er spatial</sub></td><td>0.218</td><td>0.118</td><td>0.042</td><td>0.009</td></tr><tr><td>RMSE<sub>er radius</sub></td><td>0.391</td><td>0.214</td><td>0.032</td><td>0.032</td></tr><tr><td>RMSE<sub>curv</sub></td><td>1919.428</td><td>190.531</td><td>0.060</td><td>0.035</td></tr><tr><td>L<sub>diff</sub></td><td>718.906</td><td>0.057</td><td>0.207</td><td>0.004</td></tr></tbody></table>

<table><caption>**Table 3.** Overall evaluation of the approximation methods: mean values of the quality criteria for all the centerlines distorted by spatial noise addition. The cells in gray correspond to the lowest error for each metric.</caption><thead><tr><th></th><th>GNP</th><th>GNP-AIC</th><th>GP-AIC</th><th>SRP-AIC</th></tr></thead><tbody><tr><td>RMSE<sub>spatial</sub></td><td>0.511</td><td>0.152</td><td>0.099</td><td>0.096</td></tr><tr><td>RMSE<sub>radius</sub></td><td>0.008</td><td>0.009</td><td>0.018</td><td>0.007</td></tr><tr><td>RMSE<sub>er spatial</sub></td><td>0.314</td><td>0.343</td><td>0.075</td><td>0.076</td></tr><tr><td>RMSE<sub>er radius</sub></td><td>0.015</td><td>0.019</td><td>0.021</td><td>0.013</td></tr><tr><td>RMSE<sub>curv</sub></td><td>1.524</td><td>2.362</td><td>0.085</td><td>0.091</td></tr><tr><td>L<sub>diff</sub></td><td>50.180</td><td>15.071</td><td>0.252</td><td>0.207</td></tr></tbody></table>

### 6.2.1. Comparison pipeline

Centerlines can be extracted either from the grayscale image directly or from a segmented image. In this way, centerline-based meshing can be used either as a substitute or a complement to segmentation. In our comparison study, we investigated both approaches, as illustrated in Figure 14. In what follows, the centerlines manually extracted by experts from MRA images of the BraVa database (Wright et al., 2013) are considered as reference "expert centerlines".

For the first part of our comparison pipeline (in blue in Figure 14), the MRA images are segmented by state-of-the-art segmentation methods (Tetteh et al., 2020; Livne et al., 2019). A surface mesh is produced from the segmented images by the marching cube algorithm and smoothed using a Taubin filter. The very small components of this mesh are removed to keep only the largest connected parts. A set of centerlines, which we call "segmentation-based centerlines", are extracted from the segmentation. In Section 6.2.3, they are quantitatively compared to the expert centerlines to evaluate their accuracy, and thus the topological and geometrical correctness of the segmentation-based mesh.

For the second part of our comparison pipeline (in red in Figure 14), both the expert centerlines and the segmentation-based centerlines are used as input for the centerline-based meshing methods. To improve their quality and to match the input requirements of our meshing method, the segmentation-based centerlines underwent some automatic post-processing before meshing; the small ending segments are cut out, the cycles are removed by computing a maximum spanning tree of the network, and the edges connecting the data points are re-oriented in the flow direction.

Two centerline-based methods are used to create meshes from those post-processed centerlines and the expert centerline; our explicit meshing method and the method of Abdellah

<figure><img src="image_48.png" /><figcaption>**Fig. 14.** Pipeline used to compare segmentation-based meshing and centerline-based meshing.</figcaption></figure>

et al. (2020), which is based on implicit structures and for which we used the Blender plug-in provided by the authors, in metal-balls reconstruction mode. In Section 6.2.4, those centerline-based meshes are visually compared to the segmentation-based meshes produced by deep learning methods Tetteh et al. (2020); Livne et al. (2019).

### 6.2.2. Datasets and segmentation models training

Two databases of whole-brain MRA images for healthy patients are used for this evaluation: The BraVa database gives access to 62 expert centerlines and 50 MRA, and the TubeTK database Bullitt et al. (2005) is composed of 34 MRA and the 34 in-house expert segmentations associated. These databases were chosen because they include the same type of images while offering complementary information (resp. expert centerlines and expert segmentation). For the vessels segmentation, two state-of-the-art methods were implemented, both based on neural networks: DeepVesselNet Tetteh et al. (2020) and U-net Livne et al. (2019). The DeepVesselNet architecture provided by the author was used (https://github.com/giesekow/deepvesselnet), and the U-net neural network was re-implemented. Both models were trained on the expert segmentations of the TubeTK database. 27 images were included in the training set and 7 in the test set.

The loss function used during training was a combination of dice loss and cross-entropy loss. The stochastic gradient descent algorithm was used for the optimization, with a learning rate of 0.01 for U-net and 0.001 for DeepVesselNet. The batch size was set to 5 for U-net and 10 for DeepVesselNet due to memory constraints. These hyperparameters were set empirically by testing a large selection of values for each hyperparameter. The U-net was trained for 200 epochs and DeepVesselNet for 300 epochs.

For this comparison study, we considered the segmentations produced for the 7 images of the test set of TubeTK and the 50 MRA of the BraVa database - segmented with the models trained on TubeTK -. Table 4 summarizes the data type and the number of patients of the different datasets created for this comparison study.