<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 topological 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) do 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 trifurcations 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

the arterial system and reconnect or separate vessels. The topological 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 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) offers 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 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 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
<table><caption>Table 5. topological 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 rowspan="2">BraVa</td><td>DeepVesselNet</td><td>43</td><td>369</td><td>451</td><td>59</td><td>26</td><td>76</td><td>15</td><td>0</td><td>5</td></tr><tr><td>Expert</td><td>0</td><td>205</td><td>1</td><td>205</td><td>0</td><td>102</td><td>1</td><td>1</td><td>4</td></tr><tr><td rowspan="2">BraVa</td><td>Unet</td><td>38</td><td>504</td><td>380</td><td>200</td><td>44</td><td>136</td><td>30</td><td>0</td><td>6</td></tr><tr><td>DeepVesselNet</td><td>52</td><td>552</td><td>557</td><td>212</td><td>46</td><td>140</td><td>28</td><td>0</td><td>6</td></tr><tr><td rowspan="2">TubeTK</td><td>Expert</td><td>12</td><td>551</td><td>26</td><td>508</td><td>98</td><td>230</td><td>48</td><td>0</td><td>6</td></tr><tr><td>Unet</td><td>34</td><td>626</td><td>300</td><td>446</td><td>85</td><td>215</td><td>46</td><td>0</td><td>6</td></tr></tbody></table>

<figure><img src="image_49.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>

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><img src="image_50.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>

Figure 17, together with an example branch.

As shown in 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 results 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.

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

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.
<figure><img src="image_51.png" /><figcaption>**Fig. 17.** Distribution of the scaled Jacobian values of the mesh cells. Histogram (b) represents the quality of bifurcation cells and histogram (c) the quality of the vessel cells. Image (a) illustrates the location of high and low-quality cells within a mesh.</figcaption></figure>

## 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 post-processing. If the user input data is a surface mesh, we propose the following alternative use of our meshing framework:

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.

<figure><img src="image_52.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>

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.

<figure><img src="image_53.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>

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. **The graph-based storage proposed (Section 5.2) facilitates the identification and modification of the data points of a branch of interest. Besides, many parameters of the model (e.g bifurcation cross-sections**
and apex smoothing, vessel smoothing) and mesh (e.g cell longitudinal and circumferential density, boundary layer) can be adjusted. As the bifurcations are modeled by two merging vessels, one branch can be removed without affecting the trajectory of the other branch, as illustrated on the right in Figure 20. Those modifications can be performed inexpensively by a local re-computation of the model and mesh parts.

<figure><img src="image_54.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>

## 7.3. Large cerebral arterial network meshing

To demonstrate further the applications of our method, we applied it to 60 patients of the BraVa dataset. The meshes produced for 4 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 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 oriented bifurcation model, it failed to correctly represent the

geometry. Figure 22 (a), 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 (a), even challenging topologies (e.g short connecting segments between bifurcations) can be successfully meshed with hexahedral elements. An image of all the meshes of the database, with failure areas highlighted, is given in supplementary materials, section 1.3.

<figure><img src="image_55.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>

## 7.4. CFD simulation

In this section, we demonstrate the applicability of the meshing method proposed for CFD simulations. Firstly, we com-
<figure><img src="image_56.png" /><figcaption>Fig. 21. Top view of 4 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>

pare the hexahedral meshes produced by our method to the commonly used tetrahedral meshes in terms of computational cost, convergence, and accuracy of the results. This comparison was conducted in a straight tube model, as this experimental setting allows a comparison to the analytical form of the sectional velocity profile given by the Poiseuille equation. The tube diameter was set to mimic a middle cerebral carotid artery (D = 2.5mm), and the tube length to guarantee that the flow is fully developed (L = 200mm). Five tetrahedral volume meshes of increasing cell density -from coarse to fine- were created by the software TetGen® (Hang, 2015), a state-of-the-art tetrahedral meshing software often used to produce the volume mesh in blood flow studies (Taebi et al., 2020; Shad et al., 2021). In the same way, five hexahedral volume meshes of increasing cell density were created by our method. For the CFD simulations, the fluid properties were selected to mimic blood, with a density $ρ = 1053kg.m^{-3}$, and a dynamic viscosity $μ = 0.0035kg.m^{-1}.s^{-1}$. The flow was assumed laminar -justified by a Reynold number of 150.4-. The inlet boundary condition was set to a fixed velocity $U = 0.2m.s^{-1}$, and the outlet boundary condition to zero pressure. The residuals value for convergence was fixed to $10^{-6}$. The CFD simulations were run using ANSYS Fluent (ANSYS Inc., USA).

As shown in Figure 23, 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 (as shown in Figure 23), reducing the computational cost even more. These results are consistent with the conclusions given in the works of Vinchurkar and Longest (2008), De Santis et al. (2010), and Ghaffari et al. (2017), demonstrating the advantages of hexahedral meshes over tetrahedral meshes for CFD simulations. More details on the methods and results can be found in Supplementary materials, Section 3.2.1. This experiment was reproduced in a realistic bifurcation case, giving similar results, reported in Section 3.2.2 of the supplementary materials.

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

Secondly, we applied our method to reconstruct a patient-specific mesh of the middle carotid artery (MCA) and downstream vessel in a case where the segmentation failed to produce a valid mesh. Using our framework, fluid extensions were automatically added to the inlet and outlet of the mesh and a stenosis with a user-defined shape was automatically added to the MCA. The stenosis was designed to induce a reduction of 50% of the vessel diameter. Images of the meshes are provided in Supplementary materials, Section 3.3.

Simple boundary conditions were used for this application case, as our goal is not to provide an analysis of this case study but simply to demonstrate the applicability of our method for the study of cerebrovascular pathologies by CFD. Blood is considered a Newtonian fluid ($ρ = 1053kg.m^{-3}$, $μ = 0.0035kg.m^{-1}.s^{-1}$), and the flow is assumed steady and laminar. The inlet boundary condition is set to a velocity of $0.2m.s^{-1}$ (Blackshear et al., 1980), and the outlet is set to zero pressure. The simulation converged with a residual value of $10^{-6}$ in 50
<figure><img src="image_58.png" /><figcaption>Fig. 24. 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>

iterations (resp. 78 iterations) and 4 minutes (resp. 6 minutes) for the healthy case and the pathological case respectively. The velocity streamlines and the velocity fields computed by CFD for the healthy and pathologic cases are presented in Figure 24. With this experiment, we showed the potential of our meshing method to easily design and conduct blood flow studies by CFD. The editing flexibility of our framework allows to study the effect of hemodynamic pathologies or topological changes compared to a reference geometry like in Figure 24. The advantages of our framework are not limited to the modeling and meshing steps, as it also facilitates the analysis of the results (e.g extraction of cross sections, extraction of velocity values along the centerline).

## 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 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 limits its use to cerebral vasculature. 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 of the smoothing method. Founding is acknowledged through AURA region (SIMAVC project) and ANR 20-CE45-0011 (PreSpin project).
## References

Abdallah, 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., Ingebrichtsen, 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 Graphics* 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.

Blackshear, W., Phillips, D., Chikos, P., Harley, J., Thiele, B., Strandness Jr, D., 1980. Carotid artery velocity patterns in normal and stenotic vessels. *Stroke* 11, 67–71.

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 study. *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.

Hang, S., 2015. Tetgen, a delaunay-based quality tetrahedral mesh generator. *ACM Trans. Math. Softw* 41, 11.

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.

Hodneland, E., Hanson, E., Sævareid, O., Nævdal, G., Lundervold, A., Soltészová, V., Munthe-Kaas, A.Z., Deistung, A., Reichenbach, J.R., Nordbotten, J.M., 2019. A new framework for assessing subject-specific whole brain circulation and perfusion using mri-based measurements and a multi-scale continuous flow model. *PLoS computational biology* 15, e1007073.

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.

Shad, R., Kaiser, A.D., Kong, S., Fong, R., Quach, N., Bowles, C., Kasinpila, P., Shudo, Y., Teuteberg, J., Woo, Y.J., et al., 2021. Patient-specific computational fluid dynamics reveal localized flow patterns predictive of post-left ventricular assist device aortic incompetence. Circulation: Heart Failure 14, e008034.

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.

Taebi, A., Pillai, R.M., S. Roudsari, B., Vu, C.T., Roncali, E., 2020. Computational modeling of the liver arterial blood flow for microsphere therapy: effect of boundary conditions. Bioengineering 7, 64.

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.
Click here to access/download

**Supplementary Material for on-line publication only**
*compressed_revised_supplementary_material.pdf*
<figure><img src="image_59.png" /></figure>

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

<figure><img src="image_61.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, ..., 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_62.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_63.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_64.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:

$$ \begin{cases} P_0 = S_0 \\ P_{n-1} = S_{n-1} \\ P_1 = P_0 + \alpha T_0 \\ P_{n-1} = P_{n-2} + \beta T_{n-1}, \end{cases} \quad (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_65.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_66.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_67.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_68.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_69.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_70.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_71.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_72.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(\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 $λ$. 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_73.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_74.png" /><figcaption>Fig. 13. Mesh resulting from the approximation of distorted data ($\text{density} = 1\text{mm}^{-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-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_75.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.
<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 topological 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) do 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 trifurcations 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

the arterial system and reconnect or separate vessels. The topological 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 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) offers 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 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 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
<table><caption>Table 5. topological 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 rowspan="2">BraVa</td><td>DeepVesselNet</td><td>43</td><td>369</td><td>451</td><td>59</td><td>26</td><td>76</td><td>15</td><td>0</td><td>5</td></tr><tr><td>Expert</td><td>0</td><td>205</td><td>1</td><td>205</td><td>0</td><td>102</td><td>1</td><td>1</td><td>4</td></tr><tr><td rowspan="2">BraVa</td><td>Unet</td><td>38</td><td>504</td><td>380</td><td>200</td><td>44</td><td>136</td><td>30</td><td>0</td><td>6</td></tr><tr><td>DeepVesselNet</td><td>52</td><td>552</td><td>557</td><td>212</td><td>46</td><td>140</td><td>28</td><td>0</td><td>6</td></tr><tr><td rowspan="2">TubeTK</td><td>Expert</td><td>12</td><td>551</td><td>26</td><td>508</td><td>98</td><td>230</td><td>48</td><td>0</td><td>6</td></tr><tr><td>Unet</td><td>34</td><td>626</td><td>300</td><td>446</td><td>85</td><td>215</td><td>46</td><td>0</td><td>6</td></tr></tbody></table>

<figure><img src="image_76.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>

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><img src="image_77.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>

Figure 17, together with an example branch.

As shown in 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 results 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.

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

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.
<figure><img src="image_78.png" /><figcaption>**Fig. 17.** Distribution of the scaled Jacobian values of the mesh cells. Histogram (b) represents the quality of bifurcation cells and histogram (c) the quality of the vessel cells. Image (a) illustrates the location of high and low-quality cells within a mesh.</figcaption></figure>

## 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 post-processing. If the user input data is a surface mesh, we propose the following alternative use of our meshing framework:

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.

<figure><img src="image_79.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>

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.

<figure><img src="image_80.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>

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. The graph-based storage proposed (Section 5.2) facilitates the identification and modification of the data points of a branch of interest. Besides, many parameters of the model (e.g bifurcation cross-sections
and apex smoothing, vessel smoothing) and mesh (e.g cell longitudinal and circumferential density, boundary layer) can be adjusted. As the bifurcations are modeled by two merging vessels, one branch can be removed without affecting the trajectory of the other branch, as illustrated on the right in Figure 20. Those modifications can be performed inexpensively by a local re-computation of the model and mesh parts.

<figure><img src="image_81.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>

### 7.3. Large cerebral arterial network meshing

To demonstrate further the applications of our method, we applied it to 60 patients of the BraVa dataset. The meshes produced for 4 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 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 oriented bifurcation model, it failed to correctly represent the

geometry. Figure 22 (a), 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 (a), even challenging topologies (e.g short connecting segments between bifurcations) can be successfully meshed with hexahedral elements. An image of all the meshes of the database, with failure areas highlighted, is given in supplementary materials, section 1.3.

<figure><img src="image_82.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>

### 7.4. CFD simulation

In this section, we demonstrate the applicability of the meshing method proposed for CFD simulations. Firstly, we com-
<figure><img src="image_83.png" /><figcaption>Fig. 21. Top view of 4 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>

pare the hexahedral meshes produced by our method to the commonly used tetrahedral meshes in terms of computational cost, convergence, and accuracy of the results. This comparison was conducted in a straight tube model, as this experimental setting allows a comparison to the analytical form of the sectional velocity profile given by the Poiseuille equation. The tube diameter was set to mimic a middle cerebral carotid artery (D = 2.5mm), and the tube length to guarantee that the flow is fully developed (L = 200mm). Five tetrahedral volume meshes of increasing cell density -from coarse to fine- were created by the software TetGen® (Hang, 2015), a state-of-the-art tetrahedral meshing software often used to produce the volume mesh in blood flow studies (Taebi et al., 2020; Shad et al., 2021). In the same way, five hexahedral volume meshes of increasing cell density were created by our method. For the CFD simulations, the fluid properties were selected to mimic blood, with a density $ρ = 1053kg.m^{-3}$, and a dynamic viscosity $μ = 0.0035kg.m^{-1}.s^{-1}$. The flow was assumed laminar -justified by a Reynold number of 150.4-. The inlet boundary condition was set to a fixed velocity $U = 0.2m.s^{-1}$, and the outlet boundary condition to zero pressure. The residuals value for convergence was fixed to $10^{-6}$. The CFD simulations were run using ANSYS Fluent (ANSYS Inc., USA).

As shown in Figure 23, 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 (as shown in Figure 23), reducing the computational cost even more. These results are consistent with the conclusions given in the works of Vinchurkar and Longest (2008), De Santis et al. (2010), and Ghaffari et al. (2017), demonstrating the advantages of hexahedral meshes over tetrahedral meshes for CFD simulations. More details on the methods and results can be found in Supplementary materials, Section 3.2.1. This experiment was reproduced in a realistic bifurcation case, giving similar results, reported in Section 3.2.2 of the supplementary materials.

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

Secondly, we applied our method to reconstruct a patient-specific mesh of the middle carotid artery (MCA) and downstream vessel in a case where the segmentation failed to produce a valid mesh. Using our framework, fluid extensions were automatically added to the inlet and outlet of the mesh and a stenosis with a user-defined shape was automatically added to the MCA. The stenosis was designed to induce a reduction of 50% of the vessel diameter. Images of the meshes are provided in Supplementary materials, Section 3.3.

Simple boundary conditions were used for this application case, as our goal is not to provide an analysis of this case study but simply to demonstrate the applicability of our method for the study of cerebrovascular pathologies by CFD. Blood is considered a Newtonian fluid ($ρ = 1053kg.m^{-3}$, $μ = 0.0035kg.m^{-1}.s^{-1}$), and the flow is assumed steady and laminar. The inlet boundary condition is set to a velocity of $0.2m.s^{-1}$ (Blackshear et al., 1980), and the outlet is set to zero pressure. The simulation converged with a residual value of $10^{-6}$ in 50
<figure><img src="image_85.png" /><figcaption>Fig. 24. 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>

iterations (resp. 78 iterations) and 4 minutes (resp. 6 minutes) for the healthy case and the pathological case respectively. The velocity streamlines and the velocity fields computed by CFD for the healthy and pathologic cases are presented in Figure 24. With this experiment, we showed the potential of our meshing method to easily design and conduct blood flow studies by CFD. The editing flexibility of our framework allows to study the effect of hemodynamic pathologies or topological changes compared to a reference geometry like in Figure 24. The advantages of our framework are not limited to the modeling and meshing steps, as it also facilitates the analysis of the results (e.g extraction of cross sections, extraction of velocity values along the centerline).

## 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 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 limits its use to cerebral vasculature. 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 of the 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., Ingebrichts- sen, 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.

Blackshear, W., Phillips, D., Chikos, P., Harley, J., Thiele, B., Strandness Jr, D., 1980. Carotid artery velocity patterns in normal and stenotic vessels. Stroke 11, 67–71.

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.

Hang, S., 2015. Tetgen, a delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Softw 41, 11.

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.

Hodneland, E., Hanson, E., Sævareid, O., Nævdal, G., Lundervold, A., Soltészová, V., Munthe-Kaas, A.Z., Deistung, A., Reichenbach, J.R., Nordbotten, J.M., 2019. A new framework for assessing subject-specific whole brain circulation and perfusion using mri-based measurements and a multi-scale continuous flow model. PLoS computational biology 15, e1007073.

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.

Shad, R., Kaiser, A.D., Kong, S., Fong, R., Quach, N., Bowles, C., Kasinpila, P., Shudo, Y., Teuteberg, J., Woo, Y.J., et al., 2021. Patient-specific computational fluid dynamics reveal localized flow patterns predictive of post-left ventricular assist device aortic incompetence. Circulation: Heart Failure 14, e008034.

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.

Taebi, A., Pillai, R.M., S. Roudsari, B., Vu, C.T., Roncali, E., 2020. Computational modeling of the liver arterial blood flow for microsphere therapy: effect of boundary conditions. Bioengineering 7, 64.

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.