**1. Introduction**

#### **1.1. Aneurysm**

24 Will-be-set-by-IN-TECH

[43] Saylor, P.E. (1988). Leapfrog variants of iterative methods for linear algebraic equations.

[44] Shojima, M.; Oshima, M.; Takagi, K.; Torii, R.; Hayakawa, M.; Katada, K.; Morita, A.; Kirino, T. (2004). Magnitude and Role of Wall Shear Stress on Cerebral Aneurysm. Computational Fluid Dynamic Study of 20 Middle Cerebral Artery Aneurysms. *Stroke*,

[45] Sprengers, M.E.; Schaafsma, J.; van Rooij, W.J.; Sluzewski, M.; Rinkel, G.J.E.; Velthuis, B.K.; van Rijn, J.C.; Majoie, C.B. (2008). Stability of Intracranial Aneurysms Adequately Occluded 6 Months after Coiling: a 3T MR Angiography Multicenter Long-Term Follow-Up Study. *AJNR American Journal of Neuroradiololgy*, Vol. 29, pp. 1768-1774. doi:

[46] van Rooij, W.J. (1998). Endovascular Treatment of Cerebral Aneurysms. PhD Thesis,

[47] van Rooij, W.J.; Sprengers, M.E.; Sluzewski, M.; Beute, G.N. (2007). Intracranial aneurysms that repeatedly reopen over time after coiling: imaging characteristics and treatment outcome. *Neuroradiology*, Vol. 49, pp. 343-349. doi: 10.1007/s00234-006-0200-2 [48] van Rooij, W.J.; Sprengers, M.E.; de Gast, A.N.; Peluso, A.P.P.; Sluzewski, M. (2008). 3D Rotational Angiography: The New Gold Standard in the Detection of Additional Intracranial Aneurysms. *AJNR American Journal of Neuroradiology*, Vol. 29, pp. 976-979.

[49] Vasilyev, O.V.; Kevlahan, N.K.R. (2002) Hybrid wavelet collocation-Brinkman penalization method for complex geometry flows. *International Journal for Numerical*

[50] Verstappen, R.W.C.P. & Veldman, A.E.P. (2003). Symmetry-Preserving Discretization of

[51] Wanke, I.; Dörfler, A.; Forsting, M. (2008). Intracranial Aneurysms. *Intracranial Vascular*

[52] Wiebers, D.O.; Whisnant, J.P.; Forbes, G. et al. (1998). Unruptured intracranial aneurysms – risk of rupture and risks of surgical intervention. International Study of Unruptured Intracranial Aneurysms Investigators. *New England Journal of Medicine*, Vol.

[53] Young, D.F.; Munson, B.R.; Okiishi, T.H. (1997). A brief introduction to fluid mechanics.

*Methods in Fluids*, Vol. 40, No. 3-4, pp. 531-538. doi: 10.1002/fld.307

Turbulent Flow. *Journal of Computational Physics*, Vol. 187, pp. 343-368

*Malformations and Aneurysms*. Springer, ISBN: 978-3-540-32919-0

*Journal of Computational and Applied Mathematics*, Vol. 24, pp. 169-193

Vol. 35, pp. 2500-2505. doi: 10.1161/01.STR.0000144648.89172.0f

10.3174/ajnr.A1181

doi: 10.3174/ajnr.A0964

339, No. 24, pp. 1725-1733

John Wiley & Sons, Inc., ISBN: 0-471-13771-5

University of Utrecht. ISBN: 90-393-1818-2

An aneurysm is an abnormal widening of a blood vessel. As the vessel widens, it also gets thinner and weaker, with an increasing risk of rupture. Aneurysms are essentially found in the aorta, the popliteal artery, mesenteric artery, and cerebral arteries. Intracranial aneurysms are smaller than other types of aneurysm and mostly saccular. Though most patients do not experience rupture, it can lead to a stroke, brain damage and potential death. The mortality rate after rupture is considerably high: the incidence of sudden death was estimated to be 12.4% and death rate ranged from 32% to 67% after the hemorrhage [16]. Each year, over 12,000 people die in the United States due to rupture of intracranial aneurysms [17].

In order to prevent the rupture, or rerupture, of an aneurysm, several treatments have proved successful: neurosurgical clipping, endovascular coiling and stenting. The aneurysm can be permanently sealed from the normal blood circulation by placing a tiny metal clip across the aneurysm neck. This open surgery requires to perform a craniotomy, which is invasive and associated with risks of complications during or shortly after surgery. In recent years, the development of interventional radiology techniques made it possible for a growing number of patients to be treated with minimally invasive strategies, essentially endovascular coiling. The procedure of coil embolization starts with the insertion of a catheter into the femoral artery, which is then advanced through the arterial system all the way to the location of the intracranial aneurysm. Then the radiologist delivers the filling material (the coils) through a micro-catheter into the aneurysm sac. The presence of coils in the aneurysm reduces blood velocity, and decreases the pressure against the aneurysmal wall, progressively creating a favorable hemodynamic environment for thrombus embolization. Finally, the formation of a blood clot blocks off the aneurysm, thus considerably reducing the risk of rupture. In the case of irregularly-shaped or fusiform aneurysms, or aneurysms with wide necks, stenting of the parent artery can be used in combination with coils. Although endovascular techniques

©2012 Wei et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 Will-be-set-by-IN-TECH 224 Aneurysm

are less invasive than open surgery, they remain complex and risky to perform, requiring an important expertise and careful planning.

#### **1.2. Computer-based medical simulation**

Medical simulation provides a solution to the current need for residency training and procedure planning, by allowing trainees to experience realistic scenarios, and by repeatedly practicing without putting patients at risk. With the ongoing advances in biomechanics, algorithmics, computer graphics, software design and parallelism, computer-based medical simulation is playing an increasingly important role in this area, particularly by providing access to a wide variety of clinical scenarios, patient-specific data, and reduced training cost. Even for experienced physicians, medical simulation has the potential to provide planning and preoperative rehearsal for patient-specific cases. In some cases it can also offer some insight into the procedure outcome. Nevertheless, several fundamental problems remain to be solved for a wide and reliable use of computer-based planning systems in a clinical setup, and in particular for coil embolization. Such a planning system is expected to have a high level of realism and good predictive abilities. In particular we are interested in a prediction of the hemodynamics before and after the procedure, an evaluation of the number of coils required to achieve embolization, and an interactive simulation of coil deployment for rehearsal. This involves the following challenges:


#### **1.3. Previous work**

Generally speaking, there are three main approaches to obtain hemodynamic data. Experimental techniques have been widely used in clinical analysis, but restricted to idealized geometries or surgically created structures in animals. However, a better method is *in vivo* analysis, which can be provided through medical imaging. For example, angiography can provide *in vivo* flow data, but only in 1D; Magnetic Resonance Imaging (MRI) data is 3D, but quite noisy. Finally, the computational approach, which can provide a 3D representation of detailed flow patterns in patient-specific geometry, becomes more and more attractive in this area.

2 Will-be-set-by-IN-TECH

are less invasive than open surgery, they remain complex and risky to perform, requiring an

Medical simulation provides a solution to the current need for residency training and procedure planning, by allowing trainees to experience realistic scenarios, and by repeatedly practicing without putting patients at risk. With the ongoing advances in biomechanics, algorithmics, computer graphics, software design and parallelism, computer-based medical simulation is playing an increasingly important role in this area, particularly by providing access to a wide variety of clinical scenarios, patient-specific data, and reduced training cost. Even for experienced physicians, medical simulation has the potential to provide planning and preoperative rehearsal for patient-specific cases. In some cases it can also offer some insight into the procedure outcome. Nevertheless, several fundamental problems remain to be solved for a wide and reliable use of computer-based planning systems in a clinical setup, and in particular for coil embolization. Such a planning system is expected to have a high level of realism and good predictive abilities. In particular we are interested in a prediction of the hemodynamics before and after the procedure, an evaluation of the number of coils required to achieve embolization, and an interactive simulation of coil deployment for rehearsal. This

• **Geometry modeling**: although 3D imaging of the patient's vasculature is now widely available under various modalities, extracting the actual geometry of the blood vessels is still an issue due to the limitations in spatial resolution and the presence of noises and artifacts. Particularly, accurate and exact geometry extraction of blood vessel is a challenge in the vicinity of intracranial aneurysms due to the small size and complex shape of the

• *In vivo* **data**: for the purpose of realism, patient-specific *in vivo* data is necessary for biomechanical modeling and hemodynamic simulation, such as mechanical parameters of vessel wall, blood flow rate. Acquisition of *in vivo* data is quite challenging, because imaging techniques are not currently capable to provide images with a high resolution

• **Fast computation**: the planning system targets at assistance for rapid decision making or rehearsal in a virtual environment. As such it requires fast or even real-time computation of blood flow and blood-structure interaction, which cannot always be guaranteed by modern computers with certain limitations in memory and frequency. Therefore, the need still remains to increase computing speed by optimizing algorithms or proposing

• **Coil modeling**: modeling intracranial coils, which are thin platinum wires, remains challenging, because mechanical parameters are not provided by device manufacturers. As coils are designed to conform to the aneurysm wall, it is also important to compute the

Generally speaking, there are three main approaches to obtain hemodynamic data. Experimental techniques have been widely used in clinical analysis, but restricted to idealized

multiple contacts (or self-contacts) between the coils and the aneurysm.

important expertise and careful planning.

involves the following challenges:

surrounding vascular network.

either in space or in time.

alternative numerical methods.

**1.3. Previous work**

**1.2. Computer-based medical simulation**

To obtain patient-specific geometry, excellent review papers [21] reported on the vast literature that addressed blood vessel segmentation. The diversity in the methods reflected the variety of contexts in which the question of blood vessel segmentation araised, requiring us to be more specific about our expectancies. First, the vessel surface model should be smooth enough for its use in the simulation. Implicit surface representations, where the surface was defined as the zero-level set of a known function *f* , were arguably well suited [5]. C1-continuous models (with C0-continuous normal) allowed for much smoother sliding contacts. Unwanted friction might occur with polyhedral surfaces [26] or level-sets such as vesselness criteria [11, 34] defined on a discrete grid. Implicit surfaces also offered an improved collision management over parametric surfaces [10]. Indeed, the implicit function value at a point telled whether this point was inside or outside the surface, detecting a collision in the latter case. Furthermore, the implicit function gradient gave a natural direction for the contact force used to handle this collision. In many cases, segmentation or vessel enhancement aimed at improving vessel visibility. When quantitative assessment was required, most previous works grounded on the same seminal idea as Masutani [22], which tracked the vessel centerline by fitting local vessel models: graphics primitives [43], convolution kernel [33], or cylinder templates [12]. Both the centerline location and radius estimation might be correct, but such models were unable to accurately cope with irregularities on the vessel surface, especially when the vessel section was not circular (or elliptic). Also, they did not correctly handle pathologies such as aneurysms. Various methods were proposed to improve the vessel delineation in cross-sections. Ray-casting methods were of particular interest, as they were able to extract candidate points at the vessel surface. Indeed, they opened the path to using a wealth of techniques developed by the graphics community to fit an implicit surface to a set of scattered surface points. Radial Basis Functions (RBF) had a promising potential, especially in their variational formulation [37] with the capacity to produce compact models. Closer to our work, Schumann [35] used Multi-level Partition of Unity (MPU) implicits to get a locally defined model. However, RBF or MPU implicit gradients gave appropriate contact directions close to the surface but could mislead contact forces elsewhere. We propose in Section 2 a new and efficient algorithm to model local blood vessels with blobby models [28].

Most computational methods for blood flow simulation relied on the science of Computational Fluid Dynamcis (CFD), and approximated blood flow as a continuous incompressible Newtonian fluid, described by the unsteady incompressible Navier-Stokes equations [23]. For instance, Groden et al. constructed a simple geometrical model by only straight cylinders and spheres to approximate an actual aneurysm, and simulated the flow by solving the Navier-Stokes equations [14]. The geometry model they used could not accurately describe the real patient's case, therefore, had little use in surgery planning for a specific patient. In the last ten years, the successful effort in the research of combining image processing and CFD made it possible to compute patient-specific hemodynamic information in this community. Kakalis et al. employed patient-specific data to get more realistic flow patterns [19]. However, both of their methods, as well as most similar studies, relied on the CFD commercial software to simulate the flow, and the computational times (dozens of hours in general) were incompatible with interactive simulation or even clinical practice. In order to reach fast computation, several template-based methods were designed [24] [25]. These methods pre-computed hemodynamics on a set of similar template geometries, and then interpolated on a patient-specific geometry instead of fluid simulation. However, they required to set up a pre-computed database, and were only tested on simple artery structures. But for the brain vessels around intracranial aneurysm, the network and shape are much more complicated, so the high accuracy cannot be expected.

More methods for fast computation were proposed in the field of computer graphics, essentially required the results to be visually convincing, but not physically accurate. The stable fluids approach [39] was a significant milestone, as it brought in fluid advection and the Helmholtz-Hodge decomposition to ensure the mass conservation law. However, this approach relied on discretization of Eulerian space using a regular grid, thus making it inappropriate for complex geometry with irregular boundaries, as it is usually the case in anatomical structures. Recently, the Discrete Exterior Calculus (DEC) method, based on unstructured mesh for incompressible fluid simulation, was proposed in [9]. It presented several benefits in terms of stability and computational efficiency. However, the context, in which this method was applied, did not aimed at the practical use in medical simulation. We further assess the stability, accuracy and computational time of this method, more importantly, improve the method for medical applications in (near) real-time.

Previous work in the area of real-time or near real-time simulation for interventional radiology mainly focused on training, for instance, [29], [15], [1], [7] and [5], which proposed approaches for modeling either catheter deformation or more general catheter navigation in vascular networks. Real-time simulation of coil embolization was investigated more recently. However, additional difficulties need to be faced: the deformation of the coil is more complex, and the coil is self-colliding during embolization while it also has frictional contact with the vessel wall. A FEM model, based on beam element allowed to adequately capture the deformation of the coil [6]. This model was validated, and an inverse problem was solved to capture the correct rest-shape configuration. Moreover, a solution to contact and self-collision was proposed, which was based on the resolution of Signorini's law and Coulomb friction [4]. A Gaussian deformation, close to the Boussinesq approximation model, allowed to take into account the very small deformation of the aneurysm wall. This approach is developed in more details in Section 3. Recently, a deformation model for the coil, based on inextensible elastic rods was proposed [38], in combination with a geometric approach for coil embolization, based on path-planning [27].

#### **1.4. Our contributions**

In our work, we aim at real-time simulations of the blood flow and blood-structure interaction during aneurysm coil embolization in the patient-specific geometry. While the existing studies of aneurysm embolization only concerned the impact of the deployed coil(s) on the blood flow, we also present an effective way to compute the reverse effect of the blood flow on the coil during the medical procedure, and provide higher reality for the simulation of placing coils into the aneurysm compared to the case without considering the interactive force from blood flow.

First of all, a smooth and accurate model of the vessel wall around aneurysm is required. We propose to use a blobby model whose parameterization is constrained so that the implicit function varies with the distance to the vessel surface. Thereafter, a stable blob selection and subdivision process are designed. Also, an original energy formulation is given under closed-form, allowing for an efficient minimization.

Secondly, we present a coil model based on the serially-linked beams that can handle large geometric deformations. Our model can also handle various shape memories in order to simulate different types of coils (like helical, bird-cage or 3D coils). This model is computed efficiently by a backward Euler scheme combined with a linear solver that takes into account the particularities of our model.

Thirdly, in order to achieve accurate and near real-time blood flow simulation, we introduce the DEC method into hemodynamic simulation for the first time. Compared to the DEC method initially applied in the field of computer graphics, we improve the numerical stability by using more advanced backtracking schemes, and more importantly by optimizing quality of the mesh used in the computation. Moreover, a detailed analysis of the results and comparison with a reference software are performed to understand the stability and accuracy of the method, as well as the factors affecting these two aspects.

Finally, based on the DEC method, we propose a framework for real-time simulation of the interventional radiology procedure (Figure 1). We reconstruct 3D models of vascular structures, and propose a real-time finite element approach for computing the behavior of flexible medical devices, such as coils, catheters and guide wires. Then we present the methods for computing contacts between virtual medical devices and virtual blood vessels. Furthermore, we model bilateral interactions between blood flow and deformable medical devices for real-time simulation of coil embolization. Our simulated results of blood-coil interaction show that our approach permits to describe the influence between coils and blood flow during coil embolization, and that an optimal trade-off between accuracy and computation time can be obtained.

#### **2. Geometry modeling**

4 Will-be-set-by-IN-TECH

the CFD commercial software to simulate the flow, and the computational times (dozens of hours in general) were incompatible with interactive simulation or even clinical practice. In order to reach fast computation, several template-based methods were designed [24] [25]. These methods pre-computed hemodynamics on a set of similar template geometries, and then interpolated on a patient-specific geometry instead of fluid simulation. However, they required to set up a pre-computed database, and were only tested on simple artery structures. But for the brain vessels around intracranial aneurysm, the network and shape are much more

More methods for fast computation were proposed in the field of computer graphics, essentially required the results to be visually convincing, but not physically accurate. The stable fluids approach [39] was a significant milestone, as it brought in fluid advection and the Helmholtz-Hodge decomposition to ensure the mass conservation law. However, this approach relied on discretization of Eulerian space using a regular grid, thus making it inappropriate for complex geometry with irregular boundaries, as it is usually the case in anatomical structures. Recently, the Discrete Exterior Calculus (DEC) method, based on unstructured mesh for incompressible fluid simulation, was proposed in [9]. It presented several benefits in terms of stability and computational efficiency. However, the context, in which this method was applied, did not aimed at the practical use in medical simulation. We further assess the stability, accuracy and computational time of this method, more importantly,

Previous work in the area of real-time or near real-time simulation for interventional radiology mainly focused on training, for instance, [29], [15], [1], [7] and [5], which proposed approaches for modeling either catheter deformation or more general catheter navigation in vascular networks. Real-time simulation of coil embolization was investigated more recently. However, additional difficulties need to be faced: the deformation of the coil is more complex, and the coil is self-colliding during embolization while it also has frictional contact with the vessel wall. A FEM model, based on beam element allowed to adequately capture the deformation of the coil [6]. This model was validated, and an inverse problem was solved to capture the correct rest-shape configuration. Moreover, a solution to contact and self-collision was proposed, which was based on the resolution of Signorini's law and Coulomb friction [4]. A Gaussian deformation, close to the Boussinesq approximation model, allowed to take into account the very small deformation of the aneurysm wall. This approach is developed in more details in Section 3. Recently, a deformation model for the coil, based on inextensible elastic rods was proposed [38], in combination with a geometric approach for coil embolization,

In our work, we aim at real-time simulations of the blood flow and blood-structure interaction during aneurysm coil embolization in the patient-specific geometry. While the existing studies of aneurysm embolization only concerned the impact of the deployed coil(s) on the blood flow, we also present an effective way to compute the reverse effect of the blood flow on the coil during the medical procedure, and provide higher reality for the simulation of placing coils into the aneurysm compared to the case without considering the interactive force from blood

complicated, so the high accuracy cannot be expected.

improve the method for medical applications in (near) real-time.

based on path-planning [27].

**1.4. Our contributions**

flow.

In this section, we present our work on vascular modeling using implicit representations. We show how such models can be obtained from actual patient data (3D rotational angiography, 3DRA).

#### **2.1. Local implicit modeling**

#### *2.1.1. Implicit formulation as a blobby model*

An implicit isosurface generated by a point-set skeleton is expressed as the zero-level set of a function *f* , a sum of implicit spheres:

$$f(\mathbf{X}; p) = T - \sum\_{j=1}^{N\_b} \alpha\_j \phi \left( \frac{|\mathbf{X} - \mathbf{C}\_j|}{\rho\_j} \right) \mathbf{A}$$

(a) From patient-specific data (left), we reconstruct vessel surface model (middle), and then generate tetrahedral mesh in the region of aneurysm (right).

(b) Given the boundary conditions, we simulate the blood flow velocity, from which we compute the interactive force applied on the coil. Until there is a significant increase of coil density, we update the velocity field.

**Figure 1.** The framework of real-time simulation of coil embolization: (a) mesh generation; (b) simulation of bilateral interaction of blood and medical devices.

where *T* is the isosurface threshold, {*αj*} are positive weights, and {*Cj*} is the point set skeleton. Each implicit sphere #*j* is defined by a symmetric spherical function, centered on *Cj*, of width *<sup>ρ</sup>j*. The local field function, or kernel, is a function *<sup>φ</sup>* : **<sup>R</sup>** <sup>→</sup> **<sup>R</sup>**+, rapidly decreasing to 0 at infinity. Thus, model locality is ensured: an implicit sub-model can be extracted by merely selecting neighboring blobs. For example, all results presented below are produced using the 'Cauchy' kernel [36]:

$$\phi(\mathbf{x}) = (1 + \mathbf{x}^2/\mathbf{5})^{-2}\mathbf{y}$$

where the dividing factor 5 normalizes the kernel such that *φ*��(1) = 0.

Such objets are called differently depending on the kernel used [36]. Our method is not kernel-dependent, and was successfully used with the computationally less efficient Gaussian kernel. Muraki [28] was the first to use this type of model in the context of object reconstruction. Following this seminal work, we will use the terms *blob* for an implicit sphere, and *blobby models* as a generic name for implicit models.

#### 228 Aneurysm A (Near) Real-Time Simulation Method of Aneurysm Coil Embolization <sup>7</sup> A (Near) Real-Time Simulation Method of Aneurysm Coil Embolization 229

**Figure 2.** Redundancy in implicit function parameterization: (left): 2D shape to model; middle: implicit modeling using a uniform weight of 1 (*α<sup>j</sup>* = 1); right: implicit modeling using weights equal to radii (*α<sup>j</sup>* = *ρj*). The implicit function looks more like a distance function in the latter case.

#### *2.1.2. Setting the weights*

6 Will-be-set-by-IN-TECH

(a) From patient-specific data (left), we reconstruct vessel surface model (middle), and then

(b) Given the boundary conditions, we simulate the blood flow velocity, from which we compute the interactive force

where *T* is the isosurface threshold, {*αj*} are positive weights, and {*Cj*} is the point set skeleton. Each implicit sphere #*j* is defined by a symmetric spherical function, centered on *Cj*, of width *<sup>ρ</sup>j*. The local field function, or kernel, is a function *<sup>φ</sup>* : **<sup>R</sup>** <sup>→</sup> **<sup>R</sup>**+, rapidly decreasing to 0 at infinity. Thus, model locality is ensured: an implicit sub-model can be extracted by merely selecting neighboring blobs. For example, all results presented below are produced

*φ*(*x*)=(1 + *x*2/5)<sup>−</sup>2,

Such objets are called differently depending on the kernel used [36]. Our method is not kernel-dependent, and was successfully used with the computationally less efficient Gaussian kernel. Muraki [28] was the first to use this type of model in the context of object reconstruction. Following this seminal work, we will use the terms *blob* for an implicit sphere,

applied on the coil. Until there is a significant increase of coil density, we update the velocity field.

where the dividing factor 5 normalizes the kernel such that *φ*��(1) = 0.

simulation of bilateral interaction of blood and medical devices.

and *blobby models* as a generic name for implicit models.

using the 'Cauchy' kernel [36]:

**Figure 1.** The framework of real-time simulation of coil embolization: (a) mesh generation; (b)

**- -**

generate tetrahedral mesh in the region of aneurysm (right).

In the above formulation, there is a redundancy between the weight {*αj*} and the radii {*ρj*}. Figure 2 gives an example of a 2D shape (left image) with two different implicit representations (*z* = *f*(*x*, *y*)) only parameterized by the center locations and the radii: for the image in the center, all the weights {*αj*} were set to 1 ; while on the right we had *α<sup>j</sup>* = *ρj*, ∀*j*. This latter case is particularly interesting when we consider circles of different radii: there is a monotonously increasing relation between the distance function and the implicit function *f* .

In our particular simulation context, in order to help predict collisions, and have the function give a valid contact force direction, the algebraic value *f*(*X*; *p*) at point *X* should relate monotonously to the geometric distance of *X* to the surface.

As a consequence, we set *α<sup>j</sup>* = *ρj*. Thereby, as demonstrated in Figure 2, the function gradient gives a valid contact direction anywhere in space (consider the crease in the middle of the shape in the middle image). Meanwhile, redundancy in the parameters of *f* is also dismissed. This choice for parameterization will also prove useful later to efficiently select the blob to be subdivided.

#### *2.1.3. Energy formulation*

Fitting a surface to *<sup>N</sup>* points {*Pi*}1≤*i*≤*Np* can be written as an energy minimization problem [28, 42]. We propose to combine three energy terms:

$$\mathcal{E} = \mathcal{E}\_d + \mathfrak{a}\mathcal{E}\_c + \beta\mathcal{E}\_a$$

, where (*α*, *<sup>β</sup>*) <sup>∈</sup> **<sup>R</sup>**+<sup>2</sup> , and:

1.

$$\mathcal{E}\_d = \frac{1}{N\_p} \sum\_{i} f(P\_{i\prime}; p)^2$$

, which translates the algebraic relation between data points and the zero-level set. It gives a raw expression of the approximation problem.

2.

$$\mathcal{E}\_{\mathbb{C}} = \frac{1}{(N\_b(N\_b - 1))} \sum\_{j \neq k} \left( \frac{s \sqrt{\rho\_j \rho\_k}}{|\mathbb{C}\_k - \mathbb{C}\_j|} \right)^{12} - 2 \left( \frac{s \sqrt{\rho\_j \rho\_k}}{|\mathbb{C}\_k - \mathbb{C}\_j|} \right)^{6}$$

, which is Lennard-Jones energy. Each term is minimal (with value -1) for |*Cj* − *Ck*| = *s* √*ρjρk*, being repulsive for blobs closer than this distance, and attractive for blobs further away. It imposes some cohesion between neighboring blobs to avoid leakage where data points are missing, while preventing blobs from accumulating within the model.

3.

$$\mathcal{E}\_a = \frac{1}{N\_p} \sum\_i \kappa(P\_i)^2$$

*κ*(*P*) is the mean curvature. It can be computed in a closed form at any point in space from the implicit formulation [13]

$$\kappa(P) = \frac{\nabla f^t H\_f \nabla f - |\nabla f|^2 trace(H\_f)}{2|\nabla f|^3}$$

, where ∇ *f* is the implicit function gradient and *Hf* its Hessian matrix, both computed at point *P*. This energy smoothes the surface according to the minimal area criterion. In particular, the wavy effect that could stem from modeling a tubular shape with implicit spheres, is reduced.

Behind the rather classical form given above for the energy terms, it is important to notice that the whole energy is known under a closed-form expression. As a consequence, closed-form expressions were derived for its gradients with respect to the blobby model parameters {*ρj*} and {*Cj*}.

#### *2.1.4. Selection-subdivision*

The blob subdivision procedure proposed in the seminal work [28] was exhaustive and time consuming. A blob selection mechanism was added in [42], measuring the contribution of each blob to E*<sup>d</sup>* within a user-defined window, and choosing the main contributor. We experimentally noted that this technique was prone to favor small blobs, thus focusing on details, before dealing with areas roughly approximated by one large blob. This behavior is caused by this selection mechanism using the algebraic distance to the implicit surface. Our criterion relies upon the geometric distance approximation proposed by Taubin [40]:

$$d\_T(P, f\_0) = \frac{|f(P; p)|}{|\nabla f(P; p)|},$$

where *P* is a point and *f*<sup>0</sup> is the 0-level set of implicit function *f* , parameterized by *p*. As a consequence, the point *Pi*<sup>∗</sup> farthest to the surface is such that:

$$\dot{\iota}^\* = \arg\max\_{1 \le i \le N\_\mathcal{P}} d\_T(P\_{i\prime} f\_0).$$

The blob #*j* <sup>∗</sup> whose isosurface is the closest to *Pi*<sup>∗</sup> is selected (according to Taubin's distance *dT*). Note that this criterion is valid in large area because we set *α<sup>j</sup>* = *ρ<sup>j</sup>* in the definition of *f* . The subdivision step then replaces this blob with two new ones. Their width *ρ*� *<sup>j</sup>*<sup>∗</sup> is chosen such that two blobs, centered on *Cj*<sup>∗</sup> , of width *ρ*� *<sup>j</sup>*<sup>∗</sup> would have the same isosurface as one blob centered on *Cj*<sup>∗</sup> , with width *ρj*<sup>∗</sup> (the formula depends on the kernel). The first new blob is centered on *Cj*<sup>∗</sup> , while the second is translated by *ρj*∗/10 towards *Pi*<sup>∗</sup> .

(a) 2D case: initialization with 4 blobs (b) 2D case: final result (c) 3D case: initialization with one blob (d) 3D case: final result

**Figure 3.** Examples of implicit modeling. In the 2D case (two leftmost images), data points are in black, blob centers are red crosses, and the final implicit curve is in green. In the 3D case (two rightmost images), data points are in red and the implicit surface is in white.

#### *2.1.5. Optimization*

8 Will-be-set-by-IN-TECH

points are missing, while preventing blobs from accumulating within the model.

<sup>E</sup>*<sup>a</sup>* <sup>=</sup> <sup>1</sup> *Np* ∑ *i*

*κ*(*P*) =

*s*

the implicit formulation [13]

spheres, is reduced.

*2.1.4. Selection-subdivision*

and {*Cj*}.

The blob #*j*

3.

, which is Lennard-Jones energy. Each term is minimal (with value -1) for |*Cj* − *Ck*| =

√*ρjρk*, being repulsive for blobs closer than this distance, and attractive for blobs further away. It imposes some cohesion between neighboring blobs to avoid leakage where data

*κ*(*P*) is the mean curvature. It can be computed in a closed form at any point in space from

<sup>∇</sup> *<sup>f</sup> tHf* <sup>∇</sup> *<sup>f</sup>* − |∇ *<sup>f</sup>* <sup>|</sup>

, where ∇ *f* is the implicit function gradient and *Hf* its Hessian matrix, both computed at point *P*. This energy smoothes the surface according to the minimal area criterion. In particular, the wavy effect that could stem from modeling a tubular shape with implicit

Behind the rather classical form given above for the energy terms, it is important to notice that the whole energy is known under a closed-form expression. As a consequence, closed-form expressions were derived for its gradients with respect to the blobby model parameters {*ρj*}

The blob subdivision procedure proposed in the seminal work [28] was exhaustive and time consuming. A blob selection mechanism was added in [42], measuring the contribution of each blob to E*<sup>d</sup>* within a user-defined window, and choosing the main contributor. We experimentally noted that this technique was prone to favor small blobs, thus focusing on details, before dealing with areas roughly approximated by one large blob. This behavior is caused by this selection mechanism using the algebraic distance to the implicit surface. Our

criterion relies upon the geometric distance approximation proposed by Taubin [40]:

<sup>∗</sup> <sup>=</sup> arg max <sup>1</sup>≤*i*≤*Np*

*f* . The subdivision step then replaces this blob with two new ones. Their width *ρ*�

consequence, the point *Pi*<sup>∗</sup> farthest to the surface is such that:

such that two blobs, centered on *Cj*<sup>∗</sup> , of width *ρ*�

*i*

centered on *Cj*<sup>∗</sup> , while the second is translated by *ρj*∗/10 towards *Pi*<sup>∗</sup> .

*dT*(*P*, *<sup>f</sup>*0) = <sup>|</sup> *<sup>f</sup>*(*P*; *<sup>p</sup>*)<sup>|</sup>

where *P* is a point and *f*<sup>0</sup> is the 0-level set of implicit function *f* , parameterized by *p*. As a

*dT*). Note that this criterion is valid in large area because we set *α<sup>j</sup>* = *ρ<sup>j</sup>* in the definition of

centered on *Cj*<sup>∗</sup> , with width *ρj*<sup>∗</sup> (the formula depends on the kernel). The first new blob is


*dT*(*Pi*, *f*0).

<sup>∗</sup> whose isosurface is the closest to *Pi*<sup>∗</sup> is selected (according to Taubin's distance

,

*<sup>j</sup>*<sup>∗</sup> is chosen

*<sup>j</sup>*<sup>∗</sup> would have the same isosurface as one blob

*κ*(*Pi*)<sup>2</sup>

<sup>2</sup>|∇ *<sup>f</sup>* |<sup>3</sup>

<sup>2</sup>*trace*(*Hf*)

Such a gradual subdivision procedure may lead to a dramatic increase in the number of blobs, and hence the size of the optimization problem. The locality of the kernel *φ* allows us to focus the optimization onto the newly created pair of blobs. More exactly, only the new blob that is slightly misplaced is optimized, the other blobs remaining constant. The energy is minimized using Polak-Ribiere conjugate gradient (PR) algorithm, taking advantage of the closed-form expressions of both the energy and its gradients. A single minimization loop consists in one PR minimization over the center (3 variables), followed by one on the width (1 variable). In practice, a maximum of 5 loops proved sufficient.

#### **2.2. Overall modeling algorithm**

This fitting procedure proved to be very robust to initialization. Figure 3 gives two examples, one in 2D and the other in 3D, of typical initializations and results obtained.

The 2D case (Figure 3(a) and 3(b)) mimics the neighborhood of an aneurysm. Starting from 4 blobs dispatched on the aneurysm and parent vessel (Figure 3(a)), the final result on Figure 3(b) shows 20 blobs whose centers naturally span the skeletal line of the shape while providing a close fit (green line). The 3D case (Figure 3(c) and 3(d)) is an actual bilobated aneurysm. Starting from one blob located at the entrance of the aneurysm (Figure 3(c)), the procedure ends up on Figure 3(d) with 100 blobs closely fitting all bumps and other details in the shape.

Previous works [44] demonstrated that good candidate points at the vessel surface could be obtained by casting rays from center points inside the vessel, and keeping only the locations of minimal directional gradient (blood vessels are bright onto a darker background in angiographic images). The same strategy was followed, except that rays were thrown in directions regularly spaced on the unit 3D sphere. A user-defined threshold was applied to the gradient amplitude to remove extraneous points detected on small intensity variations that could be captured in large areas, such as the aneurysm sac or when rays were cast along the vessel direction.

As a consequence, modeling the vicinity of an aneurysm takes three steps:

#### 10 Will-be-set-by-IN-TECH 232 Aneurysm


Once we get a smooth surface model, meshing the domain is done using the interleaved optimization algorithm based on Delaunay refinement and Lloyd optimization [41], and implemented using the CGAL library 1. Each refinement step acts on the size of elements, while each optimization step acts on the shape of elements. Using this method, we can also define a size field to obtain anisotropic and non-uniform mesh. We control the fidelity of the generated mesh to the previously obtained surface model, by specifying the distance tolerance between the two surfaces.
