2.9 Illustration of the boundary conditions. The leaflets are in pink, the chordae in black. Yellow points correspond to zero-displacement nodes. . . . . . . . . . . . .

2.11

Constitutive law of mitral valve leaflet tissue in two fiber directions: circumferential (runs parallel to the annulus) and radial (is perpendicular to the circumferential direction), taken from [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF] 3.9 Computation of the sub-mesh by mesh propagation (lines correspond to geodesic paths approximation, green points are the sub-mesh points). . . . . . . . . . . . .

3.10

Comparison of geodesic paths calculated from the distance maps (lines with black markers) with the paths resulting from the algorithm (lines with asterisk markers).

Introduction

The work presented in this thesis focuses on image-based modeling of the mitral valve with the goal of reproducing the subject-specific valve behavior. This study stands at the intersection of the domains of medical imaging, computer vision and biomechanical simulation.

Medical context

Cardiovascular diseases remain one of major causes of death globally according to World Health Organization [WHO]. Often the development of such medical conditions is associated with mitral valve malfunction. Since the anatomical background is essential for the problem understanding, we present first the valve anatomy along with the most frequent valve pathologies.

Mitral valve anatomy

The mitral valve is one of the four cardiac valves. It separates the left atrium from the left ventricle (see figure 1). During the left ventricle relaxation (diastole) the valve opens and blood fills the left ventricle (see figure 1 1

The mitral valve is constituted of the flat anterior and posterior leaflets tethered by the chordae attached by the papillary muscles on the ventricle walls (see figure 2). The leaflets are asymmetric: the anterior leaflet is wider and occupies only one third of the annular circumference, the rest corresponds to the posterior leaflet [START_REF] Carpentier | Carpentier's Reconstructive Valve Surgery: From Valve Analysis to Valve Reconstruction[END_REF]. Chordae are often classified according to their attachment to the leaflets [Mur09, GWF15]: first-order (marginal or free-zone) chordae attach to the free margin of the leaflets, second-order (rough-zone) chordae are connected to the ventricular surface of the leaflets, and third-order or basal chordae originate from the ventricle wall (see figure 3). The two largest second-order chordae on the anterior leaflet side are called strut. Strut chordae are reported to be essential for ventricle shape maintenance [START_REF] Muresian | The clinical anatomy of the mitral valve[END_REF].

The line where the leaflets are attached to the ventricle and atrium walls is called annulus (shown in red in figure 2 right). It defines the hinge for the leaflets movement. In fact the annulus does not exist in the anterior leaflet as it is continuous with the aorto-mitral curtain which extends from the aortic valve [START_REF] Carpentier | Carpentier's Reconstructive Valve Surgery: From Valve Analysis to Valve Reconstruction[END_REF].

The mitral valve leaflets are subjected to the highest pressure since they have the largest surface area of the four heart valves. This may be one of the possible causes of the valve pathologies development [START_REF] Lau | Mitral valve dynamics in structural and fluid-structure interaction models[END_REF].

Mitral valve diseases and treatment

There are three main types of problems associated with the mitral valve: stenosis, prolapse and regurgitation.

Valve stenosis is the narrowing or blockage of the valve, which becomes dome-shaped (see figure 4(a)). It can be caused by rheumatic fever or calcification of the leaflets and annulus Medical treatment of the valve pathologies involves surgical intervention either for the valve replacement or for the valve repair. Repair is preferred as it gives better results in the long term [START_REF] David | Late outcomes of mitral valve repair for floppy valves: Implications for asymptomatic patients[END_REF]vdMC17].

Common repair techniques involve modifying the valve tissues, such as leaflet resection, annuloplasty, chordae shortening, or adding of new chordae [START_REF] Chemtob | Choosing the "best" surgical techniques for mitral valve repair: Lessons from the literature: Chemtob et al[END_REF]. Leaflet resection is performed by excising the prolapsing part of leaflet (see figure 5 Most of these procedures, except for adding artificial chordae, require an open heart surgery where the heart is stopped. This means that the valve competence (ability to prevent the back-flow) cannot be checked before the blood flow is restored. The surgeon has to anticipate the valve behavior resulting from the gesture performed during the intervention. This task is challenging due to high intra-subject variability of mitral valve morphology and complex valve behavior. As a result, the repair outcome is highly dependent on the surgeon's experience [B + 20].

Computer-based modeling could facilitate the pre-operative planning. A computational model reproducing the valve geometry and its mechanical behavior would allow clinicians to simulate medical gestures and their outcomes, and in such a way to choose the most appropriate repair technique. This would improve the quality and the rate of success of treatment procedures. Such computational models would also enable the analysis of the function of healthy and pathological valves, which is important for global scientific understanding. Despite significant progress in mitral valve modeling, no patient-specific framework applicable in clinics has been proposed yet. Patient-specific modeling is defined by the use of the patient medical data such that it can be conveniently exploited to perform the individualized simulation. A variety of studies in the domain contributed to different aspects of the problem and provided insights into a more realistic valve representation and better understanding of its function. Depending on the main focus of the study, we can distinguish several groups of contributions.

The first group includes biomechanical models relying on a simplified valve geometry. In the pioneering works [START_REF] Kunzelman | Annular dilatation increases stress in the mitral valve and delays coaptation: a finite element computer model1presented at the 45th international congress of the european society for cardiovascular surgery, venice, italy, september 1996.1[END_REF][START_REF] Cochran | Effect of papillary muscle position on mitral valve function: Relationship to homografts[END_REF], only post-mortem dissection measurements were exploited to represent the global features of valve components. With the development of imaging modalities, both generic geometry and manually or semi-automatically segmented image-based geometry were used [START_REF] Prot | Finite element analysis of the mitral apparatus: Annulus shape effect and chordal force distribution[END_REF][START_REF] Skallerud | Modeling active muscle contraction in mitral valve leaflets during systole: a first approach[END_REF]. These works attempted to define the influence of the valve components alteration on its functional replication. Although the valve geometry was not fully based on subject data, simulations with such models allowed to show that annulus shape, papillary muscles position and contraction have an important impact on the reproduced valve behavior.

In contrast to these approaches, another category of methods relies on exact valve geometry and focuses on the reproduction of mitral valve function based on biomechanical models that can be very complex. The geometry can be synthesized (not image-based) like in [CMK17, KXW + 18], or segmented from images as in [LOR + 13, TEBI + 17, DKS17, KRD + 19]. The data acquisition for such image-based models was performed in an in vitro environment from one or several excised valves, resulting in clean data with the minimum of artifacts. These works demonstrated their ability to predict the valve biomechanical response, but in a highly controlled context, where the ground truth of mechanical parameters corresponding to the closed valve state was available. Such conditions are difficult to reproduce and do not apply to a clinical setting. In addition, the data went trough an arduous post-processing stage aiming to adapt the mesh manually and ensure the simulation convergence.

Subject-specific models are the closest to clinical application in a longer perspective. In such models, valve geometry representation involves image segmentation using different modalities. Ultra-sound imaging allows to extract the annulus and partly the leaflets, but the chordae are almost not distinguishable. The chordae geometry is completed based on prior anatomical knowledge [SMAC + 11, BSM + 13, SSW + 16, PSO + 17, GSZ + 19]. Computed tomography has higher resolution and makes it possible to distinguish the chordae. However, most of the time, the valve components are segmented manually, including several interactive post-processing stages [WS13, PKM + 17]. This is a tedious task requiring narrow expertise. It cannot be easily reproduced and performed for different patients. Such models allowed to obtain leaflet profiles similar to those resulting from in vivo data, and to get closer to a more realistic simulation.

The ability of a model to faithfully reproduce the valve behavior is essential for its application in repair planning. Preliminary applications of valve models were shown by simulating surgical operations like annuloplasty [THS + 11, SMAC + 11, SSW + 16] and chordae addition [THS + 11, GSZ + 19]. Though the valve function depends on the complex interaction between all the components, there is no method for complete geometry representation of the valve, especially few approaches take into account the actual geometry of the chordae. An accurate geometrical and topological model of chordae plays an important role: whereas some works used a simplified model of the chordae for simulations [KDB + 17], Feng demonstrated in [FQG + 18] that more precise chordal structure could help to achieve more realistic simulations of the valve behavior. Furthermore, extraction of the chordae network with the minimum of human interaction would allow to decrease the manual segmentation bias and thus improve simulation reliability.

In this thesis, we aim to set up a biomechanical model and guide it with real data extracted from computer vision procedures. The global objective of the study is to provide a method for realistic simulation of the valve behavior, which would be a step towards future clinical applications. We use the existing biomechanical model [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF], which will be eventually applied to simulate a quasi-static closed valve state with the goal to ensure that the valve geometry is appropriate. In particular, we aim to provide a representation of the chordae being compatible with the requirements of the biomechanical model. This study was carried out as a part of the associated team Inria CURATIVE3 between MAGRIT team and Harvard Biorobotics Lab. This gave an access to the home-made biomechanical model framework [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF] and acquisition of the data, and included numerous exchanges with experts in the domain.

There are several challenges to overcome in this work. First, it is important to define the imaging modality for data acquisition. A compromise between data resolution and proximity to clinical conditions has to be found. Like in previous studies of Harvard Biorobotics Lab group [THS + 11, HdNH11, HSdNH11], we choose to use micro-CT scans obtained from animal subjects. High spatial resolution of micro-CT allows to distinguish all the valve components, including the chordae. Nevertheless, images may be noisy, suffer from numerous artifacts or be incomplete (detailed in chapter 2). This is at the origin of many difficulties: contrary to usual mechanical applications, which are most of the time guided with a restricted set of accurate boundary conditions, our model will be guided by a lot of uncertain and dense images features. Without preprocessing, regularization or explicit consideration of uncertainty, this usually leads to divergence due to excessive strain level.

As our data acquisition was performed in the closed valve state, the valve components separation is not well determined and requires specific knowledge of the valve anatomy and mechanics. This complicates both automatic segmentation procedure and the validation, since the ground truth segmentation is laborious. Moreover, there is no well established validation protocol for such problems, and we would have to design one.

The choice of the biomechanical model type has to be motivated by a trade-off between precision and usability in terms of reasonable simulation time. The level of biomecanical model simplification allowing to remain in agreement with the effective function of the mitral valve is an open question. Among those proposed in the literature, we can find more precise but timeconsuming fluid structure interaction models, and faster structural models, yet not appropriate for dynamic simulations. As only the static closed state of the valve will be simulated, we have chosen the structural model proposed by our collaborators [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF] due to its relative speed comparing to fluid structure interaction models, fair precision, and availability of a code implementation.

Finally, geometry segmentation and biomechanical modeling are two joint procedures. The biomechanical model itself defines the requirements for the segmentation method. Geometry representation must be realized in terms of primitives used in the biomechanical model and be consistent with mechanical constraints (chordae are straightened in the closed valve state). And vice versa: the information that can be potentially retrieved from the data may limit the complexity of the model that can be considered.

Thesis organization

The rest of this manuscript is organized as follows. In chapter 1 we give an overview of the state of the art of mitral valve models. Chapter 2 covers the description of the dataset and the biomechanical model chosen for the simulation.

In order to incorporate the image-based geometry into the biomechanical model, an approach for automatic chordae geometry extraction is proposed, which includes three consecutive steps:

• coarse segmentation of the chordae with a topology-based method (chapter 3). The method allows to separate the chordae from the rest of the tissues by exploiting the chordae tubular nature.

• chordae tracking based on a local model-fitting approach (chapter 4). The chordae are represented in a form of trees of line segments, which are compatible with the chosen biomechanical model.

• optimization approach for correction of possible anatomical and mechanical inconsistency in the chordal trees (chapter 5). In this procedure, prior knowledge is considered and coupled with image data in order to obtain a more reliable input for further biomechanical simulation.

The extracted chordal tree is a complex structure with numerous branchings in the upper part, which is hard to evaluate. To be able to verify the coherence of the segmented geometry with the ground truth, we propose a validation procedure based on graph-similarity (chapter 6). This allows to quantitatively assess the automatic segmentation quality.

Finally, the chordal geometry segmented with the method proposed in chapters 3, 4, 5 is integrated into the biomechanical model for simulation. The closed valve state is reproduced and compared to the outcome obtained using the ground truth chordae geometry. To evaluate the influence of the chordae network completeness, additionally a pathological case is simulated. Validation criteria are based on the verification of the consistency of obtained simulation outcomes along with mechanical values.

In the last chapter, limitations of our methods and perspectives are discussed.

Chapter 1

State of the art in mitral valve modeling

Various mitral valve (MV) modeling methods proposed in the literature helped to broaden the understanding of the valve function and improve the repair techniques in a certain way. Nevertheless, a range of assumptions and approximations makes them not transferable into clinics yet.

The modeling objective is to be able to reproduce the valve behavior in the most realistic way possible. This process includes representation of the valve geometry and definition of biomechanical model, which combines boundary conditions, mechanical parameters of the valve tissues and resolution methods. Each of these aspects is examined in this chapter.

The majority of works in the literature investigate mechanical simulation rather than geometry extraction. The main contributions in the valve mechanical modeling are presented in section 1.2. In these studies few attention is given to geometry extraction step. The valve components geometry can be generic and simplified or obtained from image data. The image-based geometry is either not complete (not all the valve components can be extracted) or very detailed containing little noise, but requiring controlled acquisition environment and manual image segmentation with elaborated post-processing. We aim to perform patient-specific simulation and hence automatically segment the valve components. An overview of existing MV models classified by geometry extraction method is given in section 1.1.

A challenging point of the modeling is the validation procedure. Its definition depends on focus and setting of the study. The most common validation approach consists in checking qualitatively and quantitatively the coherence of the simulated behavior. For instance, valve competence is ensured for a healthy valve, whereas bulging leaflets and leakage are observed for a pathological case. When biomecanical characteristics are available, the simulated outcomes can be compared to the ground truth measurements obtained in in vitro setting or to the expected values, e.g. checking the stress distribution. Geometrical metrics can be used to compare simulated and ground truth valve profiles. The validation protocol and challenges associated with MV modeling are discussed in sections 1.3 and 1.4.

Mitral valve geometry representation

The MV geometry can be obtained by direct laboratory measurements of a dissected valve or from image data. The first method was widely applied in the early works. The latter is the most common approach nowadays due to the wide use of imaging in cardiac clinical procedures.

Medical imaging modalities used in cardiology include ultrasound (transesophageal echocardiogram), magnetic resonance imaging (MRI) and computed tomography (CT). Depending on the type of imaging procedure different levels of image resolution can be achieved.

Transesophageal echocardiogram (TEE) is widely used during pre-operative periods. The probe with an ultrasound transducer is placed in the patient's esophagus. An example of a TEE image is shown in figure 1.1(a). MRI is performed with a strong magnetic field around the area to be imaged and is used as an alternative to echography in the case of poor acoustic windows [START_REF] Peter | Cardiovascular magnetic resonance imaging for valvular heart disease[END_REF], see an example in figure 1.1(b). Both TEE and MRI are not invasive and allow to observe the dynamics of the valve, but TEE modality has a relatively simple acquisition protocol comparing to MRI. The main drawback of TEE is related to globally poor image resolution. The annulus and the leaflets can be identified partially, as well as PMs in some cases, but the chordal structure is not visible. A better image quality can be achieved with MRI and it can provide complementary information to that obtained with TEE [GSZ + 20].

CT imaging has higher resolution (see figure 1.1(c)). All the anatomical components of the valve can be distinguished. The disadvantage of the CT modality is that patients have to be exposed to a high dose of radiation. Given the above, the valve geometry segmentation problem can be divided in two: extraction of the leaflets and of the chordae. Depending on the modality and the expected geometry precision, the following approaches for the MV geometry acquisition can be applied:

• Entirely based on measurements after dissection.

• Based on measurements after dissection for the chordae and image-based for the leaflets.

• Image-based for both leaflets and chordae.

Geometry entirely based on measurements after dissection

In the early stages of MV modeling, the imaging modalities did not allow to obtain data with a sufficient level of precision. Therefore, the entire valve geometry was based on post-mortem measurements of dissected valves.

One of the first models was proposed by Kunzelman et al. [KCC + 93]. In this work, the valve geometry and material properties were reconstructed from porcine heart measurements. The valve structural elements (annulus, anterior and posterior leaflets, chordae, and papillary muscles (PMs) were modeled with a finite-element method (FEM)-based approach (detailed in section 1.2). Leaflets were represented as a triangulated mesh. The chordal system included different types of chordae (basal and marginal) but its complex topology was very simplified: the chordae origins were limited to two points on PMs and the branching points were not taken into account (see figure 1.2).

Further studies of the group were focused on model improvements such as addition of strut chordae and asymmetric leaflets, and simulation of pathological cases like annular dilatation [START_REF] Kunzelman | Annular dilatation increases stress in the mitral valve and delays coaptation: a finite element computer model1presented at the 45th international congress of the european society for cardiovascular surgery, venice, italy, september 1996.1[END_REF] or ischaemic dysfunction [START_REF] Cochran | Effect of papillary muscle position on mitral valve function: Relationship to homografts[END_REF]. Another approach consists in complete leaflets segmentation. In [BSD + 10, BSM + 13, GSZ + 19] the annulus was segmented manually and the leaflets were extracted in a semi-automatic way from 3D TEE images. The segmentation method proposed by Burlina et al. [BSD + 10, BSM + 13] is based on a dynamic contour approach coupled with a thin tissue detector followed by manual in-the-loop user correction (see figure 1 The models using leaflet geometry extracted from the images allow to better approximate the subject valve morphology. The missing components (chordae) are completed based on the knowledge of their anatomy from previous post-mortem studies. The usage of partially imagebased geometries is explained by the prevalence of echography data, which is more frequent and thus is more available for analysis. However, reproduction of the chordae from the post-mortem measurements implies a range of simplifications: often the origins of all the chordae are limited to two points representing the PMs and the bifurcations are not considered. In addition, such an approach requires a procedure for the chordae attachment to the leaflet surfaces, which is a source of variability of results. The chordal structure is assumed to be similar among the subjects, whereas in reality, high intra-subject variability is observed. The lack of details in the chordal structure could be a cause of the discrepancy between the simulated behavior and image data [GQF + 14]. Hence, image-based extraction of all the MV components, including the chordae, could be a way of obtaining more realistic simulations.

Geometry fully based on images

The imaging modalities allowing to acquire high-resolution data are CT and micro-CT scans. Some research groups use data obtained in conditions closer to clinics: the valve is imaged without being excised. Other groups use in vitro highly controlled setups allowing to simulate the valve movement after the valve was isolated from the rest of the heart structures. In this case, two states of the valve (systolic and diastolic) are imaged. When the valve is open, the leaflets can be easily separated, but the chordae are relaxed and not distinguishable. On the contrary, in the closed state, the chordae are stretched and visible, but the leaflets contact surface is difficult to define.

Geometry obtained in conditions close to clinics

The CT imaging has lower resolution than the micro-CT, but its main advantage is that it can be performed in vivo. In the studies [START_REF] Wang | Finite element modeling of mitral valve dynamic deformation using patient-specific multi-slices computed tomography scans[END_REF] and [FQG + 18] CT scan of a healthy human valve was used (see figure 1.6). Feng et al. have shown that detailed chordae geometry gives better results than the simplified topology used previously: leaflets coapt (seal) properly and less leaflet bulging is observed. Pham et al. experimented on pathological human valves [PKM + 17]. In these works, the geometries of the valve components were obtained manually with the use of software for segmentation and post-processing. Nevertheless, in the case of CT scans, as it was reported by Tautz et al. [TNH + 18] the leaflets are not easily distinguishable. As a solution, a semi-automatic method was proposed. In this work, the valvular axis is estimated automatically with a slice-wise region growing approach. Then leaflet surface is fitted to landmarks that were manually put by the user on 18 planes both on leaflets and on the orifice (see figure 1.7). Comparing to CT, micro-CT scanning allows to perform more detailed manual or semimanual segmentation of the entire valve. The limitation of this modality is high radiation levels and the need to explant the subject heart before the acquisition, which makes it not feasible in vivo.

In the work of Tenenholtz et al. [THS + 11] micro-CT scans were used for the chordae and the leaflets extraction. With this approach, a surgical planning system for MV repair was proposed requiring manual segmentation for both leaflets and chordae. In the works of Hammer and Villard et al. [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF], [VHP + 18] the leaflets and chordae were extracted manually with adhoc software from porcine micro-CT scans.

Geometry extraction proposed in mentioned works is a very time-consuming procedure that cannot be performed easily for large amounts of data.

Geometry obtained from excised valve with an in vitro setup

The segmentation of the valve components is easier when the valve is excised. In this case, the annulus is already identified and the left ventricle walls are absent, which simplifies the leaflet segmentation. In order to be able to image the valve in different states, it is placed in a heart simulator (shown in figure 1.8(a)). The leaflet surface is tracked with fiducial markers uniformly placed on their surface (see figure 1.8(c)).

Such an approach is used in the works of Lee et The valve components obtained in the work of Drach, Khalighi et al. [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF], [START_REF] Amir | Multi-resolution models for the mitral valve image-based simulations[END_REF] are displayed in figure 1.8(b). The segmentation steps included a parametrized representation for the chordae [KDB + 17] as a graph and as a super-quadric surface for the leaflets [KDG + 18]. Regardless of high precision of the obtained geometry, the use of such a technique is limited as it is not reproducible in real-life conditions.

Since the leaflets and chordae are extracted separately in diastolic and systolic states, the chordal connection to the leaflet surface must be done. In some works this is performed manually [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF], [VHP + 18] or with the use of fiducial markers (see figure 1.8(c)) allowing to match both states [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF], [TEBI + 17] and thus infer the connection positions. The latter method may induce mechanical properties alterations due to the markers fixation on the surface of the leaflet. Another commonly used technique consists in connecting the extremities of the chordae to the leaflet surface and adding the branching chordae within a sphere of certain radius (see in figure 1.9). However, the exact way of defining the chordae connection, the radius, and the number of branching segments is not explained. Moreover, the manual procedure may introduce bias.

To sum up, automatic extraction of the valve components remains challenging, especially the chordae segmentation, which is essential for eventual patient-specific applications.

Biomechanical models

A MV biomechanical model includes boundary conditions, constitutive laws describing tissuelevel mechanical behavior, resolution method and discretization method. For each of these elements several approaches could be applied.

Boundary conditions

Imposed boundary conditions may include displacement of the annulus and PMs due to the ventricular contraction. The importance of the PMs position was shown in [START_REF] Cochran | Effect of papillary muscle position on mitral valve function: Relationship to homografts[END_REF]. Annular motion is taken into consideration in most of the works when the valve dynamics is simulated, e.g. [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states. Philosophical transactions of the Royal Society of London[END_REF], [SMAC + 11], [START_REF] Rim | The effect of patient-specific annular motion on dynamic simulation of mitral valve function[END_REF]. In the case when a stable valve state is simulated, points Figure 1.9: Connection of the chordae extremities to the leaflet surface with a fan-like structure proposed in [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF].

corresponding to annulus and PMs are used as zero-displacement boundary conditions [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF]. Chordae attachment points on the leaflets apply tension forces.

Another condition concerns modeling of the interaction between the liquid (blood) and the valve structural elements. Fluid-structure interaction (FSI)-based or "wet" models take this into account. The models where the fluid interaction is not considered are called structural or "dry" [START_REF] Lau | Mitral valve dynamics in structural and fluid-structure interaction models[END_REF].

A comparative study of FSI-based and structural models was done by Lau et al. [START_REF] Lau | Mitral valve dynamics in structural and fluid-structure interaction models[END_REF]. This analysis was performed on the valve with the geometry reconstructed from ex vivo measurements with a range of simplifications including a linear elastic model for the leaflets (detailed later in the section) and the presence of the first-order chordae only. For the comparison, the valve closure and the whole cardiac cycle were simulated with a structural model and two different FSI-based models: with tubular and valvular volume. FSI-based models require high mesh precision. Thus in [START_REF] Lau | Mitral valve dynamics in structural and fluid-structure interaction models[END_REF] sizes of tubular and valvular models were 81920 and 43008 solid brick elements, respectively in comparison to 2688 shell elements in structural model. In must be noted that in structural models often meshes with small precision are used, i.e. 624 shell elements in [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF]. The main difference between structural and FSI-based approaches is that in a structural model the pressure is applied directly to the valve structure as a boundary condition and modeled to be uniformly distributed. On the contrary, in FSI-based method there is an interaction between fluid and structural domains, thus changes in one domain influence the other [START_REF] Lau | Mitral valve dynamics in structural and fluid-structure interaction models[END_REF]. It was shown that the structural model gives similar results of predicted stress to those obtained with the FSI-based models. In the same time, the final configurations are not equivalent. Because of the assumption of uniform pressure distribution in the structural model the leaflets coaptation (contact) zone is displaced towards posterior leaflet (see figure 1.10), which is not the case in reality and is considered in the FSI-based models. The simulation with the structural model is considerably faster (almost two orders of magnitude) than a FSI-based.

Globally, structural models can be used for the simulations of static valve states, although for the whole cardiac cycle the FSI-based models are required as they allow to model the full The choice of boundary conditions depends on the simulation goals. For instance, in the case of the whole cycle simulation the FSI must be considered and annular and PMs motion has to be introduced. But for the closed static state evaluation a structural model could be sufficient and less expensive in computational time.

In most of state-of-the-art studies [HdNH11, BSM + 13, RMK14], in order to model the contact sliding between the valve tissues, an approach introducing contact repulsive forces is applied. As proposed in [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF], first, leaflet elements (triangles in this case) being potentially in collision are identified. Then, the contact force is applied to the colliding leaflet nodes. The contact force direction is opposite to the external force which allows to balance them and avoid the overlap.

Constitutive model

The constitutive models for valve tissue are usually defined separately for the leaflets and the chordae since they express different mechanical behaviors.

There are several common approaches in the literature for the leaflets material modeling. Most of material models are anisotropic to take into account the orientated collagen fibers.

The simplest is linear elastic model: the least accurate but fast and acceptable for the case of small deformations [KRC97, LDSB10]. More complex models are hyperelastic. A constitutive law representing leaflets as highly non-linear material based on biaxial mechanical testing was proposed in [START_REF] Kk | A constitutive law for mitral valve tissue[END_REF]. Two fiber directions are considered: preferential direction (generally corresponding to the circumferential direction in the leaflet) and cross-fiber direction (radial direction in the leaflet) [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF]. An example of stress-strain curves is shown in figure 1.11 for two fiber directions. Another hyperelastic constitutive law proposed by Holzapfel et al. [START_REF] Holzapfel | Determination of layerspecific mechanical properties of human coronary arteries with nonatherosclerotic intimal thickening and related constitutive modeling[END_REF] was adapted in the models proposed in [PHS09, WS13], see the stress-strain relation in two fiber directions in figure 1.12. The chordae can be modeled as an elastic isotropic material like in [START_REF] Lau | Mitral valve dynamics in structural and fluid-structure interaction models[END_REF][START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF] or as incompressible hyperelastic isotropic material, as in [KEC07, PHS09, VAI + 09, RMK14].

Choice of physical model

The valve can be modeled with continuous (Lagrangian, FEM-based) or discrete (mass-spring (MS)-based) approaches. The Lagrangian approach is used by Burlina et al. [BSD + 10], [BSM + 13] being classically applied for cloth simulation. This method is based on minimizing the total energy of the system and finding the mesh displacement taking it to an equilibrium configuration [BSM + 13]. In the majority of studies [KEC07, VAI + 09, WS13, SDL + 19, GSZ + 19] FEM is applied. FEM is a numerical method for solving partial differential equations. In FEM-based models a complex geometry is modeled by a mesh. Thus, valve components are represented by finite elements having the material characteristics of the physical valve [START_REF] Kunzelman | Annular dilatation increases stress in the mitral valve and delays coaptation: a finite element computer model1presented at the 45th international congress of the european society for cardiovascular surgery, venice, italy, september 1996.1[END_REF].

In the following works [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF], [THS + 11] and [PSO + 17] the MS-based approach is exploited. In this case, the valve components are represented as a number of masses being connected by a network of springs [VLS + 13]. The main advantage of MS-based models is the possibility to obtain real-time simulation. Nevertheless, the simulation is less realistic than in the case of FEM-based models, and as reported by Pappalardo [PSO + 17], the stress computation is less reliable.

Figure 1.12: Valve stress-strain behavior for anterior (AML) and posterior (PML) leaflets in circumferential (fiber) and longitude (cross-fiber) direction. Taken from [START_REF] Wang | Finite element modeling of mitral valve dynamic deformation using patient-specific multi-slices computed tomography scans[END_REF].

Influence of meshing

Once segmented, anatomical components of the valve are discretized for their integration into a mechanical model. The leaflets could be represented as shell element triangle mesh with a constant thickness [START_REF] Kunzelman | Annular dilatation increases stress in the mitral valve and delays coaptation: a finite element computer model1presented at the 45th international congress of the european society for cardiovascular surgery, venice, italy, september 1996.1[END_REF], [START_REF] Prot | Finite element analysis of the mitral apparatus: Annulus shape effect and chordal force distribution[END_REF], [START_REF] Rim | The effect of patient-specific annular motion on dynamic simulation of mitral valve function[END_REF], [VHP + 18] or with a thickness map [KDB + 17]. The main advantage of flat structure leaflets is the numerical resolution speed. Another way of leaflets meshing is to use 3D-elements: tetrahedrons [TEBI + 17], [GQF + 14], pentahedrons and hexahedrons [START_REF] Wang | Finite element modeling of mitral valve dynamic deformation using patient-specific multi-slices computed tomography scans[END_REF] or bricks [START_REF] Prot | Nonlinear solid finite element analysis of mitral valves with heterogeneous leaflet layers[END_REF], [START_REF] Skallerud | Modeling active muscle contraction in mitral valve leaflets during systole: a first approach[END_REF]. Such an approach allows to better represent irregular thickness and get higher precision, although with slower simulations.

The chordal elements could be represented as 3D line segments with constant [HdNH11], [VHP + 18] or variable [KRD + 19] cross-section size. Another option is using three-dimensional solid elements (tetrahedrons), similarly to the leaflets [TEBI + 17]. As it was reported by Feng et al. [FQG + 18] 3D line segments representation with local thickness information could be a good approximation, since the simulation results are similar to those obtained with the 3D solid chordae.

Validation protocol

There are several approaches used in the literature in order to validate the results of the simulation of the MV behavior.

One of the most common methods is to compare the valve state obtained with a model to the ground truth. To quantify the results different measurements can be done:

• valve geometry comparison: measure of the closure quality by coaptation surface or coaptation line, distance error (point-to-surface error, point-to-point error).

• mechanical properties comparison: leaflet stress, strain, PMs and chordal tension.

These metrics depend on the available experimental environment. In a controlled setup the comparison can be done to the direct mechanical measurements. In other cases, obtained values 20 1.4. Challenges in computational valve modeling are compared to those from previous studies, like in [GQF + 14], or can be used to asses the valve behavior, for instance by stress distribution in the valve components. It must be noted that there are no methods for chordae similarity evaluation.

Challenges in computational valve modeling

The complex MV anatomy can be represented in a biomechanical model within a range of simplifications: from generic to subject-specific based on image data.

As it was mentioned in works of Gao et al. [GQF + 14] and Feng et al. [FQG + 18], the geometry has considerable influence on the valve behavior, specifically the chordae. For this reason, geometries based on measurements resulting from the valve dissection cannot be used for pre-operative subject-specific planning. However, the chordae geometry is hard to recover without high-resolution imaging modalities.

In the models based on image-based valve geometry, the valve components are extracted manually, which is time-consuming and requires a high level of expertise. Another issue associated with this manual approach is that the results vary significantly depending on the expert. The influence of the segmentation variability is studied by Vellguth [VBT + 19].

One of the main challenges in MV modeling remains incorporation of subject geometry into a mechanical model based on the image data. Providing detailed geometry is not enough for plausible simulation, as it could be not compatible with the model. To increase robustness in this joint process of geometry extraction and mechanical simulation, a promising path to follow is to minimize the human supervision in the geometry segmentation. This means development of solutions for the automatic MV components extraction, particularly of the chordae. This would allow to segment the MV geometry based solely on the imaging, and thereby advance toward subject-based modeling.

Finally, the validation procedure depends on the chosen acquisition environment and the simulation objectives. In a setting close to clinical conditions, a ground truth corresponding to mechanical measurements (e.g. stress in the valve components) cannot be retrieved. However, an image-based validation protocol can be established, such as comparison of the ground truth and simulated valve profiles. The coherence of the simulated valve behavior and stress distribution can be as well used as validation criteria.

Chapter 2

Dataset and biomechanical mitral valve model

In this chapter we give a detailed description of the image dataset (section 2.1) and of the biomechanical model (section 2.2 ). The nature of the data and the biomechanical model choice influence the segmentation procedure. The general requirements for the valve components segmentation methods and the validation strategy are presented in sections 2.3 and 2.4.

Data

As it was mentioned in chapter 1, the precise mitral valve segmentation requires a high resolution imaging modality, such as micro-CT scanning. Currently, in clinical conditions, human data cannot be acquired with such a modality due to the high radiation levels. Thus, a common practice in the field of cardiac valves modeling is the use of animal data, particularly, porcine or ovine valves. This choice is explained by the fair similarity in anatomy to human data and their wider accessibility. In our work we use ex vivo porcine hearts that are explanted prior to micro-CT imaging.

Acquisition protocol

The porcine hearts were scanned using high-resolution CT (X-TEK HMXST225). A protocol was established to ensure high-quality micro-CT scans. This included fixing the heart in a stable position to avoid micro motions and pressurizing it with air (see figure 2.3). Contrary to classical CT, the X-ray is static, and it is the subject that is moving. The acquisition time was optimized to get a fair ratio between precision (number of scanning angles) and time of exposure. This was done to avoid the valve tissues drying which can cause the artifacts. The number of scanning angles was 1600, the exposure time for each image was 500 ms, which gives the total acquisition time ≈ 15 minutes. The X-ray tube power was set to 10.5 W with U = 175 kV, and a 0.5 mm copper filter was used to achieve sufficient contrast and resolution. The result after reconstruction was a 2000 × 2000 × 2000 3D image with isotropic voxel resolution of 0.05 mm, from which a region of interest for the mitral valve was selected. These scans produced high-contrast volumetric images mainly composed of air (Hounsfield density ≈ 1000) and soft tissue (Hounsfield density ∈ [-100, 300]). We must note that in in vivo conditions, the PMs and annulus contract during the cardiac cycle. We do not simulate this movement in our acquisition procedure. Nevertheless, the acquired data is valid as such a simulation corresponds to even more challenging conditions than in vivo. This is due to the fact that the annulus is smaller during the contraction, which facilitates the closure [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF].

Dataset analysis

Overall, the dataset of 26 porcine hearts was obtained. Several examples in the dataset contain acquisition artifacts. Regardless the installed protocol, some scans have artifacts (an example is illustrated in figure 2.4(a)) which are possibly caused by valve tissue shrinking or are due to the micro rotation motions of the subject during the scanning. Another type of artifacts observed in the data is due to the presence of extrinsic objects in the images, like circles in figure 2.4(b) (left) and air tube (right). In in vivo cases such phenomena are not observed, in that respect, 3 data examples have been discarded.

In order to minimize micro movements the heart was fixed in the stable position with a foam rubber put around it. The foam could slightly deform the ventricle by pushing its walls, whereas in real-life conditions the ventricle fills with blood without being contracted by surrounding organs. This may cause abnormalities, such as adhesion of chordae to the ventricle walls (see figure 2.5). In such cases the chordae are indistinguishable from the walls tissue as all the tissues have the same level of contrast (see in figure 2.5(a)). For 18 images the majority of chordae suffer from this flaw, these data were discarded from the dataset.

In some cases the chordae are entangling as shown in figure 2.6. This makes the segmentations more challenging, but this problem does not affect the majority of the chordae.

Even though the valves have been thoroughly cleaned before the imaging, blood clots are difficult to remove and some of them appear in the images (see figure 2.7), which is not specific for in vivo valves. This represents important complications for the segmentation procedure.

The data comes from farm animals that have been intensively breed. Regardless of the high level of similarity between human and porcine hearts, our data comes from animals that may have numerous malformations non-specific for healthy subjects.

To sum up, 5 CT images out of 26 were kept for this work (see the data parameters in appendix A), since the rest of the data is affected by the artifacts that makes the segmentation not feasible (see the proportions in figure 2.8).

FEM-based mitral valve biomechanical model

The valve biomechanical model choice was based on the simulation precision and the biomechanical complexity. Since during peak systole the blood flow influence on the valve is absent, structural models can be used (see section 1.2). Among such models, FEM-based ones provide more realistic valve behavior. We use a FEM-model based on the framework from [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF] except that the mass-spring system has been replaced by a finite element continuum system like in [START_REF] Hammer | Mass-spring model for simulation of heart valve tissue mechanical behavior[END_REF].

In the model valve leaflets are represented as a triangulated surface (in pink in figure 2.9) having a constant thickness of 1 mm. This simplification is made since the leaflets thickness is significantly smaller comparing to other dimensions, so leaflets can be modeled as shell elements. Such an approach is also commonly used in the literature [GQF + 14], [START_REF] Rim | Effect of leaflet-tochordae contact interaction on computational mitral valve evaluation[END_REF]. The thickness is taken into account during the forces calculation.

The chordae are linear elastic rods with constant radius R ch of 0.5 mm (in black in figure 2.9) that only support tension. The chordae are connected to the leaflets and fixed in the initial points representing the PMs. These points together with the leaflet nodes along its annulus are treated as zero-displacement boundary condition.

The nodes in the model gather forces that are due to blood pressure, internal forces, contact forces to avoid auto-collision and tension forces where there is an attachment to the chordae. The sum of forces is used to compute the acceleration according to Newton's second law. The new positions are computed by solving the system of the differential equations.

Pressure forces. The pressure forces result from the pressure applied to the leaflet during peak systole. It corresponds to a pseudo-static state where no fluid mechanics effect is present. It is computed for each node as the weighted sum of all the pressure applied to the faces connected to this node. The direction is normal to the surface (see figure 2.9).

Internal forces. The internal forces are computed with the FEM in the context of large deformation [START_REF] Hammer | Anisotropic mass-spring method accurately simulates mitral valve closure from image-based models[END_REF].

First, the strain tensor E is computed with the nodes displacement using the Green tensor:

E = 0.5(C -I), where C = F F, (2.1)
where I is the identity matrix, C is the right Cauchy-Green deformation tensor and F is the deformation gradient. The deformation gradient F is the derivative of each component of the deformed x vector with respect to each component of the reference X vector (see figure 2.10): Secondly, the stress S is computed from the strain E with the constitutive law using the 2 nd Piola-Kirchhoff stress tensor. The constitutive law used in the model is the same as in [START_REF] Kk | A constitutive law for mitral valve tissue[END_REF] (see figure 2.11) for two valve fiber directions [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF]. It consists in computing Q with equation (2.3) and then the stress tensor terms with equation (2.4) with the following parameter choice: C = 30000 Pa, a 1 = 12, a 2 = 10, a 3 = 0, a 4 = 2, a 5 = 0 and a 6 = 0.

F =    ∂x 1 ∂X 1 ∂x 1 ∂X 2 ∂x 1 ∂X 3 ∂x 2 ∂X 1 ∂x 2 ∂X 2 ∂x 2 ∂X 3 ∂x 3 ∂X 1 ∂x 3 ∂X 2 ∂x 3 ∂X 3 .    (2.2)
Q = exp(a 1 * E 11 * E 11 + a 2 * E 22 * E 22 + 2 * a 3 * E 11 * E 22 + a 4 * E 12 * E 12 + 2 * a 5 * E 11 * E 12 + 2 * a 6 * E 22 * E 12 ).
(2.3)

       S 11 = C/2 * Q * (2 * a 1 * E 11 + 2 * a 3 * E 22 + 2 * a 5 * E 12 ) S 22 = C/2 * Q * (2 * a 2 * E 22 + 2 * a 3 * E 11 + 2 * a 6 * E 12 ) S 12 = C/2 * Q * (2 * a 4 * E 12 + 2 * a 5 * E 11 + 2 * a 6 * E 22 ) (2.4)
Finally, the forces are computed from the strain tensor:

f = -AHB    S 11 S 22 S 12    . (2.5)
where A and H are area and thickness of the triangle and B is is the strain displacement matrix computed from nodal positions in both present and reference configurations. Tension forces. A tension force F γ is applied to the leaflet node γ, where a chord is inserting it. For a chord the stress σ γ and the strain γ are given by equation (2.6) and equation (2.7).

σ γ = F γ /S ch ,
(2.6)

S ch = πR 2 ch ,
where S ch is the section area same for all the chordae given by the constant radius R ch .

γ = δL γ /L γ , (2.7)
Figure 2.11: Constitutive law of mitral valve leaflet tissue in two fiber directions: circumferential (runs parallel to the annulus) and radial (is perpendicular to the circumferential direction), taken from [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF].

where L γ and δL γ are respectively the initial length and the elongation of the chord γ.

Given equation (2.6) and equation (2.7) and by using the linear elasticity hypothesis in the elastic case, the tension force F γ is calculated as follows:

F γ = E ch δL γ /L γ S ch , (2.8)
where E ch is the chordae Young's modulus. Contact forces. The contact force is applied to the nodes located in the potential collision zones in the direction opposite to the external force equal to the sum of pressure forces. The collision detection is based on computing bounding boxes for all the triangles and detecting if they are overlapping. If an overlap is detected, it means that the node X t may intersect with the triangle defined by nodes A t ,B t , C t and by the normal N t . The distance d between the node and the triangle is computed as follows:

d = A t X t • N t .
(2.9)

The contact force F contact is defined as:

f contact =        0 if d h |f ext |exp(-K•d |fext| ) if h > d 0 |f ext | -K • d if d < 0, (2.10)
where h is the collision zone size, K is the contact stiffness, and f ext is the external force applied to the node. Practically we used h = 0.2 mm and K = 0.5 N m -1 . The value of K was chosen to be not too high to counterbalance the pressure force and avoid the leaflets crossing, and, at the same time, not too small, so the leaflets can be properly sealed. The ]h, 0] interval with the exponential force is used to have a sliding contact at the coaptation zone (see figure 2.12). For the interpenetration zone (d < 0) the contact force is linear (see figure 2.12). The collision detection is more robust with parallel surfaces and once the surfaces intersect on a large scale, the contact force is not high enough to counterbalance the surface pressure. Iteration and convergence of the simulation The partial differential equation system is solved using a semi-implicit numerical integration method. The time-step is adaptive and based on step doubling [START_REF] William | Numerical Recipes 3rd Edition: The Art of Scientific Computing[END_REF] with a target position error of 0.01 mm. The stabilized biconjugate gradient method [BBC + 94] is used to solve the sparse linear system.

There are several approaches to define the convergence criteria for the simulation. The choice depends on the parameters of the system. The list of physical (leaflet and chordae mechanical properties, blood pressure) and non-physical (leaflet sliding simulation and numerical scheme) parameters used in the model is given in table 2.1. In the model the time-step is adaptive, which means that when the target position error is reached, the time-step is increased and the simulation continues. In this case, local indicators such as convergence of the velocities per node, and global indicators, like mean/median node velocities can be used. Convergence of velocity values corresponds to the system stabilization. Since in our case the sliding contact between the leaflets is simulated by the repulsive forces, the velocities of the nodes being in the contact zone are likely to oscillate. Therefore, we will analyze the velocities globally, namely mean velocity values. We do not consider median of values since in the closed valve state the majority of the nodes is stable, and such a metric will not be indicative.

In order to illustrate the model function and its convergence properties, we performed tests on the manually extracted valve geometry used in [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF]. The valve closure was simulated (see initial and final states in figure 2.13). The values of the mean node velocities obtained in the simulation are shown in figure 2.14. The convergence of the velocities indicates that the closed state is obtained, which is coherent with the final leaflet profile (see figure 2.15). In the original work [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF] convergence criterion was based on time required for the valve closure in simulation. In our experiments, we chose to apply an approach based on node velocities analysis in order to ensure the simulation stability. More precisely, convergence of z-component of the node velocity will be verified as it corresponds to potential leaflet bulging direction, and of velocity norm as it displays global movement of the nodes.

Segmentation method requirements

A method for the automatic chordae segmentation is the main contribution of this work and it will be presented in the following chapters 3, 4, 5. The approach will be designed in agreement with the biomechanical model that we have chosen. This implies representation of the chordae in a form of connected line segments that are connected to the leaflet surface.

The leaflet segmentation will be performed manually in this work, however a number of conditions must be respected. The imaging was performed only in the closed valve state. In a contrary to the open state (see an example from valve closure simulation with the model detailed in section 2.2 in figure 2.13(a)), the leaflets are sealed and not easily distinguishable in the coaptation zone. The coaptation zone is the region where the anterior and posterior leaflets are in contact as it is shown in figure 2 The annulus (in yellow in figure 2.2) is well visible in the closed state. However, its delineation must be performed carefully, since variations in the annulus position have an impact on the simulation outcome.

Finally, two leaflets are connected as it is displayed in figure 2.16. There are commissures zones (marked in red in figure 2.16), which correspond to the tissue between anterior and posterior leaflets. These zones must be preserved when the leaflets are segmented. ). Annulus is in yellow. The leaflet free margin is in blue.

Validation strategy

The automatically segmented chordae set can be validated with both image-based and biomechanical simulation criteria.

The image-based approach requires extraction of the ground truth chordae in a form corresponding to the automatic segmentation outcome, more precisely, a set of connected line segments obtained via manual pipeline: labeling of the micro-CT images followed by manual delineation of the initial, branching, and ending points for each chord.

The validation of the segmented chordae via the biomechanical simulation must include verification of the simulation stability and the coherence of the valve behavior. The simulation is stable when the globally evaluated node velocities converge to zero without oscillations. The coherence of the valve behavior is checked by the leaflets profile. When the chordae set is complete, the valve must remain sealed, minimum leaflet bulging is observed and no leaks appear. At the same time, the absence (rupture) of chordae must induce leaflet bulging and possibly violate the valve competence.

Conclusion

The available image dataset contains data with numerous artifacts (e.g. blood clots, chordae adhesion to the ventricle walls), which are mostly due to the established post-mortem acquisition protocol. The most affected images were rejected, and the reduced dataset for further experiments was prepared.

Since only the static closed state of the valve is studied, we have chosen a structural FEMbased biomechanical model for simulations. The home-made implementation of the model is available which enables the tuning of the model parameters.

Given the characteristics of the dataset and the biomechanical model, requirements for the segmentation methods were determined, and a validation strategy was proposed. In the following chapters (3, 4 and 5) the full pipeline for the automatic chordae segmentation is presented. The image-based validation of the obtained chordae segmentation results is given in chapter 6. Biomechanical simulation of the valve closed state with the automatically segmented chordae set is detailed in chapter 7.

Chapter 3

Topology-based chordae segmentation

The images in the dataset (see section 2.1) contain various soft tissues including chordae, leaflets, and ventricle walls. All the tissues appear as white structures in the micro-CT scan due to the established experimental procedure (see section 2.1). As it can be seen in the manual segmentation slice illustrated in figure 3.1(a), the transition between chordae and leaflets is rather smooth and not distinct. Therefore, the main challenge consists in separating correctly the chordae (in red in figure 3.1(b)) from the other valve tissues.

The MV chordae have complex topology: cross-section is non-regular and variable in size, numerous branchings are present, closer to the leaflets the chordae get thinner and entangle. As it was demonstrated in [KDB + 17] the chordae follow bifurcation branching pattern meaning that at on point a chord can split only in two directions. Moreover, besides the target tubular structures, the other tissues are present in our data, and the contrast difference is not significant enough to separate them. Given the above, classical methods for vascular vessels segmentation cannot be applied (detailed in section 3.1). In this chapter, we present a method allowing to perform coarse chordae segmentation from the images. The method exploits the tubular nature of the target structures. Since the chordae have elliptic non-regular sections, rather than a purely geometry-based approach, we propose to locally classify structures in images as chordae or non-chordae depending on their topological properties. This requires a pre-processing step resulting in the image contour mesh (detailed in section 3.2). The segmentation method (section 3.3) is based on the local extraction of sub-meshes and the labeling of the sub-mesh points according to their topology shape type.

Vascular vessel segmentation methods

The MV chordae are topologically complex structures with tubular characteristics on the local level. Therefore, classical methods for vascular vessels segmentation were considered first. Among such methods we can distinguish differential filters based on Hessian matrix such as [FNVV98, KMA + 00, MVN06, BB08]. Another group of approaches includes mathematical morphology methods such as RORPO [START_REF] Merveille | Tubular Structure Filtering by Ranking Orientation Responses of Path Operators[END_REF] incorporating both structural-element-based and connectivity-based approaches.

The Frangi filtering [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] and RORPO method [START_REF] Merveille | Tubular Structure Filtering by Ranking Orientation Responses of Path Operators[END_REF] were tested on our data. Examples of the chordae segmented with the two methods after the analysis of the parameter variation are shown in figure 3.2. Even though Frangi method performed better than RORPO, more false detections are observed (see top and bottom parts in the 3D rendering in figure 3.2(c)). As can be seen, none of the considered methods allows extracting the whole chordae network and thus cannot be applied for our problem. Hence, we opt for a topology-based method described hereafter to discard non-chordae parts and apply a geometry-based method for centerline extraction in chapter 4.

Initial data processing

In the first stage, to minimize false detections, images are denoised with the non-local transformdomain filter [START_REF] Maggioni | Nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF], see an example displayed in figure 3.3. Then tissue contours are extracted using the classical 3D edge detection procedure based on the Canny filter [START_REF] Canny | A computational approach to edge detection[END_REF] (see an example in figure 3.4(a)). The resulting set of points contains little noise. A surface mesh can be deduced from Screened Poisson surface reconstruction [START_REF] Michael | Screened poisson surface reconstruction[END_REF] to obtain the connectivity information between the points and points normals as it is shown in figure 3.4(b). The parameters for the data processing steps were tuned in order to obtain the optimal values across the dataset. The Canny filtering tuning included the choice of the lower and the upper bounds of gradient intensity T inf and T sup , and the standard deviation for Gaussian smoothing σ (see table 3.1). The parameters of the Screened Poisson method were set to get a trade-off between mesh precision and complexity.

Mesh labeling method

After the valve structure is represented in form of a mesh, mesh segmentation techniques are used to extract the chordae. Topological constraints are commonly applied for mesh segmentation, but it often boils down to restrict each sub-region to be topologically equivalent to a disk [START_REF] Shamir | A survey on Mesh Segmentation Techniques[END_REF].

We propose to extract chordae components as features which are homeomorphic to a tube or to a branching zone, whereas components homeomorphic to a disk are discarded. As demonstrated in [START_REF] Georges | The classification of surfaces with boundaries[END_REF] two compact surfaces are homeomorphic if and only if their boundaries have the same number #CC of connected components (CC). Examples of surfaces with different topological characteristics are shown in figure 3.5(a). Blue ones have one CC (type-1 shapes) and are equivalent to a disk, red ones have 2 CC (type-2 shapes) and are equivalent to a tube, green and black ones have 3 and 4 CC, respectively, and are equivalent to branching zones (type-3 and more shapes). It should be noted that there are no holes in the dataset meshes, so no additional connected components are created. As the chordae section is variable in the form (ellipses with different eccentricity), in practice, we use the knowledge of an upper bound of the maximum chordae radius R max as a parameter for the segmentation. The maximal chordae circumference can be estimated as 2πR max . At each point p of the mesh, we build a sub-mesh centered on p which contains points that are at a geodesic distance from p less than πR max . This allows capturing the whole tube section. Indeed, if p belongs to a chordae, then the sub-mesh should be closed and homeomorphic to a tube or a branching zone. The sub-mesh generation procedure is described in section 3.3.1. The points are labeled depending on the number of CC in a sub-mesh. Figure 3.5(b) gives representative segmentation results on synthetic cases in dimension 1: for a tube, the radius of which is less than R max (top), the red sub-mesh generated from the green point is closed and is marked as chordae. For a point with a smaller local curvature (bottom), the generated sub-mesh in red is not closed and the point is labeled as non-chordae. Figure 3.6 displays the sub-meshes computed for real cases for the point (in red). The point is identified as the one belonging to non-chordae in figure 3.6(a), to a tube in figure 3.6(b), and to a branching zone in figure 3.6(c). The labeling process is detailed in section 3.3.2.

Sub-mesh extraction

A sub-mesh corresponds to the region enclosed by the iso-geodesic distance contour. In our case, the distance is limited by semi-circumference of the largest chord (see figure 3.7(a)). By definition, the geodesic distance is the length of the smallest curve between two points x s and x e on the surface [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF] (see figure 3.7(b)). The geodesic distance from the central point to each point in the set defines a geodesic distance map. An example in figure 3.7(c) corresponds to a discrete case where the color of the map indicates the proximity to the central point, red lines are iso-geodesic distance lines, and yellow lines are examples of geodesic paths [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF]. Considering the size of our data (> 300K points), the computation of a distance map for every point is hardly feasible in practice. Therefore, we propose an approximation scheme based on mesh connectivity. The sub-mesh region grows from the central point through its direct neighborhood (one-ring), in figure 3.8(a), the central point is in green and its direct neighbors are in red.

The procedure goes as follows. Let ps be an initial neighborhood (one-ring) of the central point p c . We limit the construction of the geodesic paths P to k directions defined by the direct neighbors ps (black lines in figure 3.8(a)). At each step of propagation in mesh the one-ring of every point p i ∈ ps is taken, excluding the points of the ps and the central one p c in order to avoid backward propagation (blue points in figure 3.8(b), noted ps i ). The point p i+1 minimizing the Euclidean distance ||p i , p i+1 || is added to the path P k (circled in red in figure 3.8(b)). The union of all the neighborhoods ps i forms new ps, which is an approximation of iso-geodesic line, see red points in figure 3.8(c). The propagation continues till the distance to ps (accumulated distances between the points in paths P k ) reaches πR max . In figure 3.9 initial, intermediate and final stages of mesh propagation are shown, lines are the approximated geodesic paths limiting the sub-mesh points (in green). This approach also takes into account cases when mesh boundaries are present: as soon as a boundary is reached, propagation in this direction stops.

The output of the sub-mesh extraction procedure is a set of points being connected by mesh to p c and laying at a distance smaller than πR max (green points in figure 3.9). The geodesic distance map was calculated for extracted sub-meshes and the geodesic paths obtained from the distance map (line with black markers in figure 3.10) were compared to approximated paths resulting from the mesh (lines with asterisk markers in figure 3.10). The obtained distance values were compared to those calculated from the distance map (see table 3.2). The mesh approximation approach gives less regular paths and thus the distances tend to be overestimated. Nevertheless, taking into consideration the small difference in value (about 10.6 %) and the significant gain in computation time, this approximation is sufficient for the correct mesh

Labeling procedure

After a sub-mesh centered on a point is extracted, its central point is labeled according to the number of sub-mesh connected components. This procedure is repeated for all the points in mesh till all of them are labeled. The segmentation result is in figure 3.11(a), tubes are shown in red, branching zones being type-3 shapes are in green and branchings zones being type-4 and more are in magenta, blue corresponds to non-chordae structures. The separation between tubular structures and non-tubular ones is correct. Branching zones are also well identified. However, this approach involves redundant calculations caused by the fact that sub-meshes centered on points being close to each other overlay. This considerably increases the computational time.

To reduce the computational time, we propose to label all the points of a sub-mesh at once and not only the central one. Such a faster procedure requires fewer iterations, but the overlays of sub-meshes still can occur, and for points in such regions, the labeling may be repeated several times. To minimize the superposition of sub-meshes, a uniformly distributed sample of the points having no label yet is selected (red points in figure 3.12(a)). This assures the spatial separation between sub-meshes corresponding to the sample points (see figure 3.12(b)). The points of each sub-mesh are labeled according to the following scheme: if a point was not labeled yet it is assigned a label according to the #CC in the current sub-mesh. Otherwise, the existing label is overwritten if it is superior to the current one. In other words, the preference is given to a lower label to minimize the false detections. This allows to re-label the points in the overlaying zones. An example is illustrated in figure 3.12(c). First labeled as a branching zone, points enclosed by rectangular are relabeled as a tube afterwords.

Segmentation results with the labels assigned during the faster labeling procedure are shown in figure 3.11(b), the color legend is the same as with the exhaustive approach where submesh is extracted for each point (figure 3.11(a)). When compared, both methods allow to separate correctly chordae and non-chordae structures. It must be noted that the branching zones are better identified with the exhaustive approach (see orange zones in figure 3.11). This can be explained by the section size variation of chordae in the data. To capture the large chordae, the maximum structure size is used as a parameter during mesh segmentation, which means that for smaller chords sub-meshes contain additional structures. On the other hand, with the faster approach, the number of false positives is decreased (see black zones in figure 3.11) since points from the overlaying regions are re-labeled. We performed tests on a Intel Xeon 3.50GHz workstation using the Matlab environment. Given the computation time gain (≈ 0.2 hours against ≈ 21.2 hours) and fair quality of results, we adopt the faster labeling approach for the segmentation.

Results

Results

The method has been tested on 5 micro-CT scans. For each image the maximum chordae size R max was estimated by an informed user (see values in table 3.3).

We have merged all the labels for chordae-type structures. In figure 3.13 the segmented chordae are shown in red, blue color corresponds to non-chordae valve components. Overall, non-chordae structures are correctly identified and detected chordae ending points match up with actual chordae attachment points. We have compared the results of the algorithm with a ground truth chordae mesh segmented manually by the informed user. The ground truth acquisition will be detailed later in chapter 6. In figure 3.14 the results of our segmentation algorithm (contours in red) are superposed with the detection provided by the informed user (labels in blue) on the axial (left), coronal (center) and sagittal (right) views. Two sets of views are provided for one point at half the maximum height of the chordae (first line) and one near the chordae ending close to the leaflets (second line). All the chordae are detected at half-height but some are missing or unduly detected in upper slices. In most of the cases, the chordae contours obtained with the proposed methods match with the manual segmentation. Even if some false detections are present, the results are promising since the main structure of the chordae is correctly detected. The method robustness was validated on 5 images from the dataset.

Conclusion

The chordae segmentation method proposed in this chapter is efficient regardless of the presence of the other tissues in the data.

The chordae mesh is extracted automatically with a procedure requiring only an upper bound of the radius and not an accurate value. This method is flexible since it detects tubes or branch- ings without requiring an explicit model of the section. The robustness of the approach was validated on 5 images. Average computational time is 12 minutes (on a Intel Xeon 3.50GHz workstation using the Matlab environment).

Despite the fact that the majority of the extracted structures correspond to the chordae, some false detections may be present. The latter correspond either to the parts of papillary muscles and leaflets in the vicinity of the chordae (due to the extraction method drawbacks) or to the blood clots and chordae entangling (originating from the data and which are not explicitly considered in the modeling process).

In the next chapter 4 the resulting chordae mesh will be used for precise geometry extraction with a model-fitting approach allowing to correct possible false detections of the topology-based method, and represent the chordae in a form suitable for the biomechanical model.

Chapter 4

Model-based extraction of chordae geometry

In this chapter we present a method allowing to refine the chordae segmentation obtained with a topology-based approach (chapter 3) and to extract the precise chordal geometry in a form that is suitable for the biomechanical valve model detailed in section 2.2. The approach is based on the local fitting of a cylindrical model to the chordae mesh subregions defined in chapter 3. A chord is tracked in a form of a tree of interconnected cylinders and branching zones. The resulting tree is then simplified into a set of connected line segments obtained as cylinder centerlines and approximated skeletons of branching zones.

Chordae as trees

The chordae mesh resulting from the extraction procedure detailed in chapter 3 mainly contains two types of structures: type-2 shapes (having 2 CC of the mesh boundaries) and branching zones (with 3 and more CC of the mesh boundaries). However, this segmentation is coarse since for the sake of computational time the labeling was performed by sub-region and not for each point. This may result in falsely detected type-1 shapes in the areas corresponding to the papillary muscles (see figure 4.1(a)) or to the leaflets (see figure 4.1(b)). For this reason we cannot not fully rely on the labeling obtained in chapter 3 for the precise chordae segmentation and have to re-estimate the topological properties. In addition, in order to be applied in the biomechanical model, the chordae have to be represented in a form of connected line segments, which requires extraction of the chordae centerlines. It is important to note that we do not aim to obtain a precise 1D chordae shape representation, but rather a structure that would be mechanically coherent with the model.

State of the art on centerline extraction methods

There are numerous 3D mesh skeletonization methods in the literature. Methods providing curve skeletons being 1D structures locally centered in a shape [TDS + 16] could be the most appropriate for our task. However, the skeletonization of 3D objects remains challenging, especially when a shape is not well isolated. In our case, the segmented chordae mesh could be connected to Red zones depict the problematic zones were the topology of the structures changes. Branching zones having classical configuration are well detected with the methods (see in circled black zones in figure 4.2). However, in our case, branching zones are large and often have non-typical topology because of the chordae adhesion (such as the one circled in magenta in figure 4.2). This complicates the skeletonization. As we are looking to disconnect chordae one from another and the centerline precision is not essential for our problem, classical methods for skeletonization were judges not applicable in our case.

Our problem is also close to vascular segmentation approaches widely developed in medical imaging. Such approaches allow to obtain tubular structures centerlines. A group of approaches is based on the Hessian matrix, like [START_REF] Zhou | Vascular structure segmentation and bifurcation detection[END_REF] where branching zones are identified using AdaBoost classification or by applying multi-scale vesselness as in [KAR + 15]. Multiple hypothesis template tracking approach proposed in [START_REF] Friman | Multiple hypothesis template tracking of small 3d vessel structures[END_REF] exploits an image-based vessel profile. Other methods are based on tracking tubular vessels with a mobile tracking area, for instance parallelepiped [START_REF] Flasque | Accurate detection of 3d tubular tree structures[END_REF], superellipsoid [TdF + 07], or cylinder [START_REF] Yureidini | Robust ransac-based blood vessel segmentation[END_REF]. RANSAC-based cylinder estimation in [START_REF] Yureidini | Robust ransac-based blood vessel segmentation[END_REF] allows to geometrically constrain the search for subsequent segments and in such a way increase the robustness. Method based on minimal path techniques where a vessel is represented as a 4D curve was proposed in [START_REF] Li | 3d multi-branch tubular surface and centerline extraction with 4d iterative key points[END_REF].

Majority of proposed centerline extraction methods are developed for structures having circular section and branching zones of regular size without large transitions (see examples in figure 4.3), whereas MV chordae have non-circular section and comprise large branching zones with non-regular profiles. Moreover, the user input of seed points is often required for the procedure [TdF + 07, LYC09, YKC12], whereas we aim to develop a procedure with the minimum of human interaction.

As a way to improve the performance of the segmentation for the cases where the structures are irregular and complex, machine learning-based approaches were proposed, for instance [SLF14] using training of a regression for the segmentation. More recently, methods exploiting deep learning techniques were also introduced, such as DeepVesselNet [TEF + 18], where the network is trained on the synthetic data, including the branching zones. Such methods neither can be applied in our case due to the differences in chordae and vessel morphology.

Given that in our data branching zones are large and have non-typical profiles, and the section of the chordae is not circular, we opt for a model-fitting approach applied on a local level without user intervention. A method with prior identification of branching zones and cylinders would allow to increase robustness. Possible inaccuracy of the centerlines is acceptable, since is not the main goal of the procedure, and the centerline will be adjusted after in order to satisfy the mechanical constraints for the simulation.

Chordae centerline extraction strategy

From this point we use chord to refer to the whole structure originating from a single initial point at a papillary muscle. In terms of topology a chord can be seen as a tree with one or several branches. Due to high variability of chordae configurations and sizes it is not possible to define a holistic geometrical model of a chord. Therefore, finding a geometrical pattern that is common for the data is a more suitable approach for our problem.

Given that the chordae are not necessarily straight, it is better to geometrically approximate them on the local level (a sub-region of the chordae mesh). The section of the tubular chordae is rather elliptical, thus its closest geometrical representation is an elliptical cylinder. A chord can be then represented with a series of such cylinders.

Nevertheless, branching zones are difficult to associate with a geometrical primitive. Typi- For the tree extraction, first, a RANSAC-like method is defined in section 4.2. The tree tracking algorithm is detailed in section 4.3.

Cylinder model fitting

Our goal is to associate a chord sub-region (sub-mesh) points with a cylinder. A common approach for such model-fitting tasks is RANSAC (RANdom SAmple Consensus) method [START_REF] Fischler | Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF] and its variants. RANSAC is an iterative non-deterministic method allowing to define a set of data points corresponding to a parametric model. The data contains points that correspond to the model (inliers) and abnormalities (outliers). The fitted model must ensure the highest percentage of inliers, whereas the outliers are rejected. The strong point of this approach is that it can optimally explain the data even when many outliers are present.

The sample size q is the minimum number of data points needed for model estimation (e.g. two points for a line, three point for a circle). At each iteration, a model is fitted to a randomly chosen sample. The consensus set (inliers) and outliers are then defined with respect to the fitted model using a model-specific loss function. As our data is a set of 3D points, the Euclidean distance is used as a loss function: the inliers are the data points closer to the model than the maximum distance . The algorithm stops when a sufficient number of iterations have been performed or when the size of the consensus set is considered sufficient.

The number of iterations depends directly on the rate of inliers present in the data. The more important it is, the less it is necessary to increase the number of iterations. A common way to approximate the number of iterations is to ensure that the confidence γ (probability of taking at least one sample consisting only of inliers) is high enough [START_REF] Zhang | Parameter Estimation Techniques: A Tutorial with Application to Conic Fitting[END_REF]. The inlier proportion is denoted ω. The number of samples m to be tested is recalculated at every iteration as follows:

γ = 1 -(1 -ω q ) m ), m = ln(1 -γ)/ln(1 -ω) q ), (4.1)
where 1 -ω q is the probability that at least one of the q points is an outlier.

In the RANSAC algorithm, the fitted models are classified by a criterion based on the size of the consensus set. The solution corresponding to the largest consensus set is retained. This does not allow to separate hypotheses which would be more precise since all the points being at a distance less than are considered equally close. For our problem, we choose a RANSAC variation MLESAC [TZ00]. In the MLESAC method, the model acceptance criterion is based on the sum of distances from both inliers and outliers to the model, while truncating the distance values exceeding the maximum distance to . We define hereafter a cylinder model and the model acceptance criterion for the MLESAC procedure.

Cylindrical model with elliptical base

A parametric model of a cylinder with an elliptical base is used for the local fitting. The parameters of such a cylinder include 1) cylinder axis direction d and the elliptical base descriptors: 2) center point c, 3) major a and 4) minor b semi-axes, and 5) ellipse rotation angle α (see figure 4.5). Height h is then estimated from the inliers set projected on the cylinder axis as difference between maximum and minimum.

A minimum of 5 points is needed to estimate the elliptical model. However, it was demonstrated in [START_REF] Gander | Least-squares fitting of circles and ellipses[END_REF] that more robust ellipse estimation with the least-squares method is obtained with 8 points (shown in green in figure 4.6). We fix then the minimum sample size for the MLE-SAC procedure to 8, which is enough for a cylinder model fitting as normals of the points are Projected points are displayed in magenta in figure 4.6(b). These points are then used for the ellipse fitting in the basis defined by P and axis d with least-squares method [START_REF] Gander | Least-squares fitting of circles and ellipses[END_REF]. An example of a fitted ellipse is shown in figure 4.7(a).

The center of the ellipse corresponds to the cylinder center c. The centers of lower and upper bases are then defined by the projection of the data points on the cylinder axis and choosing the minimum and maximum, respectively. The cylinder height h is estimated as the distance between the centers of the bases.

Model acceptance criterion

Once an ellipse is fitted to the sample points, the cylinder model can be evaluated with the distance-based criterion. All data points are projected on the ellipse (see figure 4.7(b)) and the distance threshold is used to define inliers (located at distance ≤ , shown in green in figure 4.7(b)) and outliers (located at distance > , shown in red in figure 4.7(b)). Afterward, according to MLESAC, the sum of distances Σ dist of all the points (inliers and outliers) is calculated. The model is accepted if Σ dist is smaller than the best Σ dist obtained so far.

Structures located near the leaflets could be not chordae and topologically equivalent to semitubes (see figure 4.8). To avoid such cases, a criterion of uniformity was added to the model acceptance step. More precisely, the model is accepted only if the data points projected on the cylinder base are distributed relatively uniformly. For that, we estimate the angles between the points projected on the cylinder base. The distribution is considered uniform if the maximum angular difference α 0 between two consecutive points does not exceed a threshold (see figure 4.9). The threshold was fixed to π/4 as it is enough to prune not-complete structures and it allows to be more flexible with the cylinder fitting. A small variation of the estimated cylinder axis is acceptable for the tracking but may result in a not totally uniform distribution of projected points on the cylinder base. Fixing such a threshold allows the algorithm to be more permissive and do not reject satisfactory cylinder models.

An example of a fitted cylinder model is shown in figure 4.10. The proposed model-fitting method is then used for the chordae geometry extraction procedure detailed in the next section.

Tree tracking

We develop a tracking approach to extract a chord which is a series of interconnected cylinders and branching zones. The tracking starts from the lowest data point of the chordae mesh since this corresponds to the zone of anatomical chordae origin, and in these areas, the chordae are larger, thus it is easier to define the consensus set as more points are present. Initial sub-mesh centered on the initial point (in magenta in figure 4.11(a)) is extracted in the same manner as in chapter 3 using the only parameter R max corresponding to the largest chordae size in the data.

As it was mentioned, the labeling obtained with the coarse topology-based segmentation (chapter 3) is not fully reliable, therefore during the tracking, the shape type of the sub-mesh is re-estimated. The tracking proceeds according to the shape type. If the sub-mesh is a type-2 shape, it is a potential chord and model-based cylinder fitting is performed. In the case of type-3 or more shape (potential branching), the tracking is recursively applied in n branching directions. Finally, when a type-1 shape is encountered corresponding to the remaining false detections from the coarse segmentation (chapter 3), the tracking procedure stops since no cylinder model can be fitted. The procedure of subsequent tracking schematically shown in figure 4.12 is repeated until the whole chordae is extracted (see figures 4.11(b) and 4.11(c)). More precisely, the process stops when it is not possible to get an initial point (or points) for the following sub-mesh extraction. We describe hereafter separately the procedures and the stopping conditions for the cases of type-2 and type-3 and more shapes.

Type-2 shapes

The sub-meshes being type-2 shapes are likely to be the chordae parts. In this case, a cylinder is fitted to the sub-mesh points, as shown in figure 4.11(a) using the method defined in section 4.2. Since the lower and upper bases may be inappropriate as we do not verify the density along the axis, we introduce additional checks. The fitted model is then used to verify the initial estimation of the topological properties of the sub-mesh. This is aimed to avoid missing branching zones when the sub-mesh is not covering the whole structure (as shown in blue in figure 4.13(a)). The new sub-mesh (in red in figure 4.13(b) is defined as a connected component limited by the cylinder bases (see two planes in blue in figure 4.13(b)). In the case when the newly extracted sub-mesh is not type-2 shape, it it processed accordingly (as type-1 or type-3 and more shape). For an example shown in figure 4.13, a 3-type shape was identified after the verification. This means that tracking must be divided in several directions. In the case when the newly extracted sub-mesh is still a type-2 shape, the tracking pursuits in one direction defined by the fitted cylinder.

The next tracking step requires an initial point for the subsequent sub-mesh extraction. It is chosen by projecting the points of the boundaries of the previously examined sub-mesh (in green in figure 4.11(a)) and taking the one being the lowest by z-axis (in yellow in figure 4.11(a). This allows avoiding holes in a tracked structure.

To provide coherent chordal trees, where cylinders sections do not change abruptly in size, as well as cylinder directions do not deviate, the constraints for a currently fitted cylinder are added. The currently fitted cylinder section must not vary more than 1.5 times from the section of its preceding cylinder and the angular difference between the axes of currently fitted cylinder and preceding one must be inferior to 15

• . The threshold value helps to speed up the fitting procedure and does not influence the final outcome. Thereby, processing of a type-2 shape for which a model was fitted, results in an initial point for the next sub-mesh estimation and the model parameters, which are used as constraints for the subsequent cylinder fitting. When no model was adjusted the chord tracking ends.

Type-3 and more shapes

The type-3 and more shapes are the potential branching zones. Before extracting the directions for the tracking, these structures must be validated. This is explained by the presence of numerous artifacts in the data such as blood clots and chordae adhesion (see chapter 2).

This may result in type-3 and more shapes with untypically large boundary connected components as shown in figure 4.14(a), or inverted structures as displayed in figure 4.14(b). As such sub-meshes would disturb the oriented tree structure (from bottom to top) and lead to incorrect tracking, the tracking procedure is stopped when they are identified.

The structures with untypically large boundary connected components are detected as those having the length of CC of the mesh boundary superior to the doubled length of the boundary CC corresponding to a chord with R = R max .

Badly oriented structures (see an example in figure 4.14(b)) are detected by analyzing the segments that connect the gravity centers of the boundary connected components (shown in green). The idea is to detect the gravity center G o where the two segments remaining the most inside in the sub-mesh (see in blue in figures 4.14(b) and 4.14(c)) originate. If the structure is correctly orientated, the mesh boundary with G o contains as well the point from which the sub-mesh was initialized (see in yellow in figure 4.14(c)). The opposite means that the structure is badly orientated (see G o and the sub-mesh origin point in figure 4.14(b)).

The type-3 and more sub-mesh boundaries are used to define the directions for the subsequent tree construction. The sub-mesh boundary containing the sub-mesh initial point is not taken into account. The remaining sub-mesh boundaries are checked one by one as shown with numbers in figure 4. 15(a). The tree is constructed following each of the defined directions until the fitting is possible (see figures 4.15(b) and 4.15(c)).

As the difference in the chordae section may be significant at the chordae origin and at the chordae endings, we adopt a resizing scheme at the branching zone processing step. The radius R b of the sub-mesh boundary is roughly estimated as the length of boundary CC. The radius for the subsequent sub-mesh extraction is reduced by two if the estimated radius R b ≤ Rmax 2 . This strategy allows being more flexible in the zones where chordae are splitting into numerous and small segments closer to the leaflets.

The described tracking procedure is repeated until all the data points are examined and results in a set of trees each corresponding to a separate chord.

Trees extraction results

The proposed tree extraction approach was tested on 5 CT images. The following choice of parameters was done. The maximum distance to the model was fixed to 0.2 px, which was

Trees extraction results

validated by tests performed on the data. The maximum number of iterations N iter was 2000 with confidence of 99%.

The MLESAC procedure includes the adjustment of the remaining iterations according to the inliers proportion for the fitted model. When the maximum number of iterations is high enough, the robustness of the results is ensured. But the procedure may be time-expensive close to the leaflets as more hypotheses are checked. To validate the choice of the iteration upper bound, we performed tests while increasing the N iter . The results for different values of N iter (2000, 3000, 5000) are shown in figure 4.16. The trees extracted with three different values of upper bound for the iteration number do not significantly vary. The exceptions are the regions (circled in red in figure 4.16) containing the structures such as blood clots or the leaflet parts being not complete tube (see in figure 4.17). These parts do not fit the model and give not reliable results.

Reproducibility of the results

Since the MLESAC procedure is non-deterministic, it induces some variance in the parameters of the cylinder fitted to a sub-mesh. This leads to the variability in the sub-meshes initial points and consequently the variation in the obtained results. In order to evaluate the influence of the Figure 4.17: Leaflet part being a not complete tube. Fitted model is in red. Chordae mesh is in pink.

MLESAC procedure on the results, we performed 5 test iterations of the method for each data. We show the results for one image (data 3) in figure 4. 19 and the rest of the results are given in appendix B. It appears from that study that the results are reproducible. The variation between the iterations is observed only in the uncertain zones (circled in red in figure 4. 19) that represent complications for the tracking.

Tree structure post-processing

The representation of the chordae branches as trees is an intermediate step. With respect to the biomechanical model (chapter 2), the chordae should have only interconnected line segments, whereas tree elements are subsequent cylinders and branching zones (see figure 4.4(a)). For this reason, cylinders are simplified to their centerlines. Branching zones are represented with skeleton approximation as the segments connecting the gravity centers of the boundary connected components (see figure 4.4(b)).

The connection of the subsequent cylinder centerlines and bifurcation skeletons boils into the presence of intermediate connection points (shown in blue in figure 4.20(a)), which is not an effective way to represent the chordae in the model. The chordae segments must be rectified. The intermediate connection points are omitted and only the initial and the terminal points along with the branching points are kept as displayed in figure 4. 20

(b).

The obtained representation of the chordal geometry (see an example in figure 4.20(b)) mostly agrees with the chordae boundaries in the image. Nevertheless, the chordae segments extracted from the trees may contain discontinuities, more precisely gaps between the segments belonging to the same chordal branch. Three examples in figure 4.21 illustrate the major causes of gaps in the trees. The example in figure 4.21(a) is the result of blood clots presence. The latter are nontubular irregular structures of a large size where the tracking procedure stops (see section 4.3). Another cause of discontinuities is the low image contrast which may result in an incomplete chordae mesh. An example is displayed in figure 4.21(b). Finally, discontinuities can occur in the proximity of large branching zones as shown in figure 4.21(c). In this case, the chord tracking stops due to the insufficient size of the locally extracted sub-mesh which is categorized as nontubular structure. The connection between such gaps is required because the beginning of the connected segments group denotes the papillary muscle. Hence, the group of segments starting after the gaps would misrepresent the chordae anatomy, as in reality there is no chordae origin at this point.

To achieve the connection between the discontinuities, we propose a merging procedure based on the angular proximity of the tree segment directions. The possibility of merging is considered for the trees whose initial segments have their starting points above the endings of the other trees. The spatially closest segments are defined for such trees and the collinearity constraints are checked afterward. More precisely, three directions are exploited in the process, the direction -→ s 1 of the spatially closest segment (in magenta in figure 4.22), the vector -→ s 2 connecting the ending point of the spatially closest segment and the starting point of the initial segment of the tree to be merged (in green in figure 4.22), and the direction -→ s 3 of the initial segment of the tree to be merged. Here we take into account the trees starting with a branching point like shown in figure 4.22(a). In this case, the initial segment direction is given by the mean direction of branching segments (in blue in figure 4.22(a)). Otherwise, the direction -→ s3 matches the initial vector orientation (in blue in figure 4.22(b)). For the collinearity constraints, two vectors are collinear if their cross product is a zero vector. In our case, the merging of the two segments is authorized if the following conditions are satisfied:

- → s1 × - → s2 ≤ T, - → s2 × - → s3 ≤ T, (4.3)
where T is a threshold value for the cross-product. The threshold value T is chosen experimentally in way that it is not over permissive but at the same time the angular differences due to the distance are tolerated. In practice, T is fixed to 0.3, which corresponds to 17.5

• angular distance.

Checking two cross-products defines the alignment of points of two segments s 1 and s 3 . It is necessary to avoid the connection of the badly orientated higher-lying segment whose starting point can be aligned with the lower positioned segment direction.

Using the proposed strategy the discontinuity flaw is removed in majority of cases (circled blue zones in figure 4.23).

Conclusion

The proposed approach is an effective way of unsupervised chordae geometry extraction, which only requires prior knowledge on the maximum chordae size present in the data.

The method consists in the iterative tracking and local fitting of a cylinder model to the chordae. Adopting such a technique enables the pruning of possible false detections resulting from the previous segmentation step (chapter 3) and allows to represent the chordae in a form compatible with the biomechanical model. More precisely, obtained chordae are trees of connected cylinders and branching zones, which are simplified in trees of connected line segments.

The efficiency of the method was demonstrated on 5 micro-CT images. The reproducibility of the results is limited in the uncertain zones of the chordae mesh.

The segmented chordal trees visually agree with the chordae mesh, yet anatomical and mechanical nonconformities can be observed in some cases. This problem will be addressed in the next chapter 5 before the chordae segmentation validation. Chapter 5

Adjustment of the chordae geometry with respect to the biomechanical model requirements

The chordal trees resulting from the method presented in chapter 4 represent the MV chordae in a form of connected line segments. The tree starting points are located in the zones where the chordae attach to the papillary muscles, the terminal points are in the proximity of the leaflets. The points interconnecting the segments are the branching nodes. However, the segmented trees are not directly suitable for the biomechanical model (chapter 2) and do not completely agree with the chordae anatomical characteristics (chapter 3). Therefore, adjustments have to be done to satisfy these requirements. Firstly, the coherence of the segmented chordal trees with the biomechanical model must be ensured. In our case, the initial configuration corresponds to the closed valve state, therefore the chordae have to be in tension. The chordae having slack (schematically shown in red in figure 5.1) may cause unrealistic simulation outcomes and convergence instability once the pressure is applied (indicated by the arrows in figure 5.1) such as leaflet bulging (see the leaflets profile in blue in figure 5.1(b)) or leaking. These configurations are caused by the fact that locally fitted cylinders centerlines and branching zones extremities points are connected one by one (see chapter 4). As shown in figure 5.2, the chordae obtained in this way (see figure 5.2(a)) may be loose and thus not compatible with the biomechanical model. The chordae slack is a global characteristics that can be properly measured only after the simulation. Moreover, it also depends on the leaflet elasticity. To address the task, we will focus on finding an equivalent of the chordal set with the minimum of slack being in agreement with the image data. An example of desired configuration is shown in magenta in figure 5.2(b)). The second adjustment of the chordal trees is aimed to correct abnormal topologies in branching zones. The mitral valve chordae branchings should follow the bifurcation pattern, i.e. have exactly two outgoing segments (see chapter 3), therefore three and more segments cannot start at one branching node. Such configurations must be revised in order to accurately represent the chordae geometry.

In this chapter we introduce a method to correct the faulty topology cases and then to straighten the chordal trees considering the anatomical constraints along with the image data (chordae limits). The approach is based on an optimization procedure for the branching nodes rearrangement (presented in section 5.1) with different sets of objectives. For the topology correction an image-driven objective function is proposed (section 5.2) while the straightening of the chordae is based on a trade-off between trees straightness and image-driven term (section 5.3).

Optimization procedure for chordal trees adjustment

The anatomical and mechanical coherence of the chordal trees can be achieved by modifying the branching node positions while considering the chordae limits and reducing the slack. In the case of incorrect topology, branching nodes must be added to respect the valve anatomical characteristics. An example of a chord with a not-satisfying geometry resulting from the method detailed in chapter 4 is displayed in figure 5.3(a). Examples of the desired outcomes of the adjustment procedure are shown in blue and purple in figure 5.3(b). For both examples the topology is changed to respect the bifurcation pattern of the MV chordae (a node is added) and the branching nodes are displaced. The definition of an appropriate modification is challenging as both proposed configurations of the chordal structure in figure 5.3(b) have correct topology and are more straightened while staying within the chordae limits (shown in black). For this reason, the chordae geometry adjustment procedure is formulated as an optimization problem, where the solution is a set of branching nodes. For the case of branching nodes position adjustment, the slack is minimized while allowing to maintain the segments inside the chordae limits. On the contrary, for the faulty topology case, only image constrains are applied. Slack is measured globally for the whole chord, while the topology adjustment is a local problem since an isolated branching zone is considered. To describe the geometry of a chord two angles are defined: the angle α between branching segments (in red in figure 5.4) called branching angle and the angle θ between the base segment (in green in figure 5.4) and the bisector of the branching segments (in gray in figure 5.4) called straightness angle. Starting from the observation that chordae straightening reduces slack, we have chosen to quantify the slack with a cost function C straight defined as the sum of straightness angles. Similar properties were observed in the study of the chordae geometrical characteristics [KDB + 17] and had shown that the base segment of the branch generally lies in the plane defined by the branching segments.

The optimization process will be driven by various constraints. Firstly, each chord segment must remain inside the chordae limits in the image. Two additional constraints are the threshold values of the straightness and the branching angles that must stay in the intervals [0 Three types of constraints influence the choice of the tree: mechanical, image-based and anatomical ones. Mechanical constraints aim at recovering a tree with the minimum slack while images-based constraints ensure that the optimized tree is inside the chordae mesh. Finally, known anatomical constraints on the chordae which mainly consist in linear constraints with limit angles θ and α for the branching structure are incorporated.

Except for the anatomical constraints which are linear, building a framework which integrates all these constraints is difficult. Formulating the problem as a linear combination of three energy terms is not appropriate since the measurements do not concern the same physical quantity: pixels distances for the image term, angle for the straightness term. Another possibility would be to perform constrained optimization, one of the terms being optimized under constraints provided by the two others. However, considering straightness as a constraint is not possible since it is not straightforward to relate a stable state of minimal slack with limits on C straight .

For these reasons, we opt for a constrained framework which minimizes C straight under imagebased and anatomical constraints.

We define a binary function F inIm (p) indicating whether a point of the chordae segment is within the chordae limits, where p is a 3D point belonging to the segment s j (b), j ∈ [1, m]. n straightness angles θ i (b) and n branching angles α i (b), i ∈ [1, n], equally depend on the given set of branching points {b}. The straightness term C straight (b) is defined as follows:

C straight (b) = n i=1 θ i (b), b ∈ R 3 × n.
(5.1)

We have two groups of criteria: straightness cost to quantify slack and image-based term. Whereas the image-based criterion can be defined as a constraint with a binary function F inIm (p), for the straightness term C straight the threshold values are not known. Therefore, we choose to minimize the function C straight under the image-based and angles constraints. It must be noted that the minimum of the function C straight can be easily reached when all the branching points are displaced to the chord origin. To avoid this, we add a constraint to limit the space for branching points b coordinates by the neighborhood of their initial positions denoted N b .

The constrained optimization problem is then defined as:

min b C straight (b), 0 ≤ θ i (b) ≤ 30 • , ∀i ∈ [1, n], 0 ≤ α i (b) ≤ 90 • , ∀i ∈ [1, n], F inIm (p) = 1, ∀p ∈ s(b), b ∈ N b , b ∈ R 3 × n.
(5.

2)

The solution of the optimization problem is a set of branching nodes {b} satisfying the constrains (equation (5.2)) and minimizing the cost C straight (equation (5.1)). However, a set of branching node coordinates belonging to the neighborhood N b , minimizing the sum of the straightness angles and having all the segments completely inside the chordae limits does not necessarily exist. An example of such a chord for which a solution under the constraints cannot be found is displayed in figure 5.5(a). Displacement of the branching nodes limited by the neighborhood N b and bringing all segments inside the chordae (in magenta in figure 5.5(a)) does not allow to reach a sum of the straightness angles being lower than the initial one. In other cases, the solution may be not satisfactory in terms of proximity of the segments to the chordae centerlines, like the one shown in figure 5.5(b). The solution corresponding to the green segments allows to minimize the slack, but at the same time it is close to the border, since only the belonging to the chordae inner volume is required.

The binary image constraint is too strict. In order to make it more permissive an imagebased term C im expressing the proximity to the image skeleton must be taken into account. As a solution we apply a multi-objective optimization procedure with two objective functions: the straightness cost function C straight and the image-based cost function C im indicating the distance of the segments to the chordae center. Since there is no single solution simultaneously minimizing both functions, we need to order the set of objectives. The concept of Pareto dominance is widely used to compare the objective vectors [START_REF] Ngatchou | Pareto multi objective optimization[END_REF]. We rank the two functions by the order of importance and then apply a lexicographic method for the minimization. As reduced slack is essential for the physical behavior of the model, we choose the angle sum cost function C straight as the first considered criterion while enforcing the proximity to chordae center is considered at the second step. Each time a local minimum of C straight is reached the image-based cost C im is evaluated, and if a local minimum of C im is reached, the current solution is updated. We formulate the optimization problem as such:

1) min b C straight (b) = n i=1 θ i (b), 2) min b C im (b) = m j=1 C s (s j (b)), 0 ≤ θ i (b) ≤ 30 • , ∀i ∈ [1, n], 0 ≤ α i (b) ≤ 90 • , ∀i ∈ [1, n], b ∈ N b , b ∈ R 3 × n, (5.3)
where C s is an image-data term for chord segments, which is defined later. The constraints in this case include branching and straightness angles limits and the neighborhood N b limiting the branching points possible coordinates.

In the optimization procedure for the case of topology correction only the image-based cost function is minimized.

Image-based cost function C im

The image-based cost function C im for a given set {b} of branching nodes is defined based on a distance map D indicating how far a segment point is located from the center of the image structure (see an example in figure 5.6(a)). To simplify the calculation, each segment is discretized, resulting in a set of l points p

s(b) t , t ∈ [1, l],
shown in red in figure 5.11. In this case, the cost function of a segment C s is the sum of corresponding values of the distance map for its points:

C s (s(b)) = l t=1 D(p s(b) t ).
(5.4)

The values of the distance map in the 3D case correspond to the distance from a given voxel to the skeleton of the volume of the detected chordae. The skeleton itself is a geometrically simplified representation of the shape. Therefore, given a binary image, the classical methods for the skeletonization based on thinning operation ([LKC94]) can be applied. Whereas skeletonization of the whole volume would not permit to correctly extract the form of the chordae. We consider only the part of the image containing chordae, i.e. the chordal trees resulting from the tracking procedure based on cylindrical model fitting (chapter 4). In this case, the representation in a form of connected cylinders and branching zones is used as this allows to skeletonize the tree elements one by one. For the tree segments corresponding to the cylinders the skeletons are defined as cylinder center lines (see an example of corresponding meshes in blue and fitted cylinders in magenta in figure 5.7). The yellow points correspond to the voxels of the obtained binary volume. For the branching zones (see the corresponding sub-meshes in red in figure 5.7) the sub-mesh voxelization is performed and the skeletons are extracted using the classic techniques for binary volumes ( [START_REF] Lee | Building skeleton models via 3-d medial surface axis thinning algorithms[END_REF]). As it can be seen in the close-up in figure 5.7, the skeleton of the bifurcation zone (in yellow) adequately represents the shape and does not require any post-processing such as deletion of the short branchings on the extremities. The distance map to the chordae skeleton is then calculated. However, if the chordae limits are not considered (shown in red in figure 5.6(a)), the segments going out of the limits are not enough penalized by the distance value difference. For this reason a cost map allowing to evaluate the proximity of the segments to the chordae skeleton and thus penalize the segment crossing the chordae limits is built based on the distance map values.

A gap is introduced between the distance map values corresponding to the interior and the exterior of the chordae. At the same time, for the interior, the difference between distance values in the center of the structure and on the border is enhanced. To define the chordae interior we use the binary volume V of the mesh (in cyan in figure 5.9) resulting from the topology-based segmentation method (chapter 3). This allows to avoid the discontinuities in the cost map due to the holes that are present in the chordae skeleton (see an example in figure 5.7, the whole corresponds to the zone where no cylinder of branching zone was identified).

Considering the presence of configurations where the segments cannot remain entirely inside the chordae limits (figure 5.5(a)), an additional buffer volume V b between exterior and interior is introduced. Thus, segments slightly going out of the chordae borders are not over penalized. The buffer volume is obtained by dilation of the volume V with a sphere of size equal to 0.2 • R max and subtraction of the volume V (see magenta zones in figure 5.9). The size of buffer volume is chosen by tests performed on the available dataset. For the example shown in figure 5.9 R max = 10 px, so the buffer zone size is equal to 2 px.

A cost map allowing to evaluate the proximity of the segments to the chordae skeleton is built based on the distance map values depending on the zone (inner, buffer or outer volume). To favor the proximity to the chordae center the distance values are taken to power 2 for the voxels being inside the chordae limits. A power function with a higher factor increase the influence of the holes in the skeleton as shown in figure 5.8 (scales are different for the two cases). For the points being in the buffer volume, the distance values are taken to the power 2.5. This authorizes the segments to go out slightly of the chordae limits while still prioritizing the position inside the structure. Finally, for the chordae exterior, the cost map values are homogenized by assigning a constant value dependent on the size of the largest chordae in data equal to R 3 max . This value is high enough to penalize the undesired configurations. Thus, the cost map value C(p i,j,k ) for a volume point p i,j,k with coordinates i, j, k is defined as follows:

C(p i,j,k ) =      D(p i,j,k ) 2 if p i,j,k ∈ V D(p i,j,k ) 2.5 if p i,j,k ∈ V b R 3 max otherwise , (5.5)
where D(p i,j,k ) is the distance from a point p i,j,k to the skeleton. The profiles of the cost map obtained based on the distance map values are shown in logarithmic scale in figure 5.10, for an example with R max = 10 px and for a chord of size R = R max . In this case, the lower and upper limits of the buffer zone are D lim = 10 px and D lim + 0.2 • R max = 12 px, respectively. A 2D slice of the calculated distance map C is displayed in figure 5.6(b).

The estimation of the image-based cost function C im given a set {b} of n branching nodes must be done for each segment s j of the chord, j ∈ [1, m]. As before each segment is discretized with regular spacing, resulting in a set of l points p ) for the point p s j (b) t is obtained by rounding its coordinates (an example for a point in figure 5.11, the result of rounding is shown in green). We do not perform interpolation in order to estimate the cost map value due to the fact that the experiments had shown that rounding gives equivalent results and is considerably faster in terms of computation. The image-based cost function value is calculated as a sum over the costs of discrete points of all chord segments: ).

C im (b) =
(5.6)

Optimization method choice The defined cost functions are non-convex with multiple local minima. Taking this into account, the convergence of classic iterative optimization methods such as [START_REF] Lagarias | Convergence properties of the nelder-mead simplex method in low dimensions[END_REF] is expensive in time. Furthermore, the found result is not guaranteed to be optimal as the initial guess may be far from the desired solution.

Therefore, we opt for an exhaustive procedure with a discrete set of candidate solutions. In our case a candidate solution b ∈ R 3 × n is the coordinates of n branching nodes. The set of candidate solutions is built as follows. First, the list of possible candidate points L i is defined for each branching point, i ∈ [1, n]. By possible candidate points for a branching node we mean those being located inside the chordae limits and in proximity to the branching node initial coordinates (see blue points in figure 5.12). This can be assured by the extraction of the points enclosed by the sub-mesh centered at the initial branching nodes. The procedure is detailed in the following sections.

Afterwards, the candidate solutions are defined by exploring all combinations of the points from the L i i ∈ [1, n]. In the general case, the number of candidate solutions is The elements of the obtained candidate solution set are iteratively evaluated with the defined objective functions.

N sol = |L 1 | • |L i | • ... • |L n |,
For the topology correction case the solution is the set of branching nodes minimizing only the image-based cost function C im defined in equation (5.3). The detailed procedure of chordae topology optimization is presented in the next section.

For the branching nodes position adjustment two minimal values C 0 straight and C 0 im are found by exhaustive search while checking all the combinations of possible bifurcation node positions. Each time a local minimum of the straightness is reached, the value of C im is evaluated. If this value is less than C 0 im , then C 0 straight and C 0 im are updated. The choice of the sets L i considered for discrete optimization is detailed in section 5.3.

Optimization of chordae topology

The chordal trees follow the binary branching pattern, i.e. all branching points are bifurcations, which excludes presence of three and more chordae segments sharing the same initial point (chapter 3). At the chordal trees extraction step, some branching zones are too large and enclose several bifurcations being located close one to another as displayed in blue the zones in figure 5. 13. This results in a distorted topology in which a branching point is an origin for three and more chord segments (in red).

In order to correct the chordae topology we need to generate a set of compatible topology configurations and, afterwards, choose the best one in terms of the proximity to the chordae center. Consequently, the optimization process that gives the coordinates of branching nodes as Figure 5.13: Example of a chord (data 1) with faulty topology in two orthogonal views. The terminal points are in yellow, the initial branching node is in cyan, the branching nodes after the correction are in blue. Circled zones correspond to the overlapping bifurcations. The chordae mesh is in pink. solution is repeated for every feasible topology configuration of the chord.

Generation of the topological configurations set

Let us explain the generation of the set of feasible topological configurations for the simplest case when a chord splits into three segments at one point, schematic representation is shown in figure 5.14(a) left. In order to correct this topology, according to the branching rules, we need to add one bifurcation point. A possible solution is given in green in figure 5.14(a) right. Thus, at every branching point, the chord splits only into two segments. We illustrate as well the cases when a chord splits into four and five segments, schematically shown in figures 5.14(b) and 5.14(c). Two and three bifurcation points must be added, respectively, in order to correct the topology (see possible solutions in figures 5.14(b) and 5.14(c)). In the general case, the number of bifurcation points needed to result in a given number of ending points N e is defined as N e -1. Since the initial bifurcation point is already present in the data, the number of bifurcations to be added is

N bnew = N e -2.
The way of connecting the exiting bifurcation point (in yellow in figure 5.14) with added ones (in green) and ending points (in gray) is not unique as it can be seen in figures 5.14(b) and 5.14(c). This plurality of connection ways forms the limited set of possible topological configurations that should be examined. In fact, the number of ending nodes N e is sufficient to define the set of topological configurations of a chord, each configuration being a strictly binary tree (nodes have exactly two children). There are two types of nodes in such binary trees: bifurcations allowing further division and ending points being the tree leaves. To define the set of possibilities we use the following scheme of possible divisions starting from initial bifurcation node (in cyan in figure 5.15). The bifurcation must have two child nodes. Possible pairs of child nodes are bifurcation-bifurcation (level 1 left in the scheme in figure 5.15), bifurcation-ending point (level 1 center) and two ending points (level 1 right). This corresponds to the first level of branching. Further split is possible only in cases where a bifurcation node is present (see second and third levels in the scheme in figure 5. 15).

The number of ending nodes N e defines the number of bifurcations to be added N bnew . To retrieve a feasible topological configuration it is enough to perform the traversal of the scheme of possible chordae splits till the number of needed bifurcation nodes is reached. In order to obtain all the topological configurations, the scheme should be traversed in predefined order until no feasible configuration is found. For that, the traversal order is fixed from top to bottom and from left to right.

To demonstrate the procedure, we consider the example shown in figure 5.14(c), for which N bnew = 5 -2 = 3 bifurcations must be added. Therefore, at the first level of branching two possibilities are acceptable: two bifurcations or one bifurcation and one ending point (see figure 5.16). The split into two ending points is rejected because in this case no bifurcations could be added at all. As the left direction is checked first, at the second level the initial bifurcation splits into two bifurcation nodes (blue nodes in figure 5.16), which leaves the possibility to add only one bifurcation node at the second branching level. The only solution satisfying this constraint is the split of one of the bifurcation nodes into one bifurcation and one ending point and the other one into two ending points (second branching level in figure 5.16). The limit of the needed bifurcation additions is reached, so the only remaining bifurcation has two ending points as child nodes (third branching level in figure 5.16). The resulting traversal path indicated by blue arrows in figure 5.16 corresponds to the highlighted topological configuration in figure 5.14(c). The leaves of the retrieved solution tree are marked as those not to be considered for the following traversals (crossed nodes in figure 5.16). The rest of topology configurations is Figure 5.15: Scheme of possible branching in the chordae: cyan is the initial bifurcation node, gray corresponds to the bifurcations to be added, red corresponds to the ending points. Blue dashed lines correspond to possible branching type, black lines unfold the two nodes in the possible branching type. Rectangular zones limit the levels of branching. extracted in the same manner while checking only the unmarked nodes in the traversal order defined previously.

Once the set of topological configurations is defined, it is expanded by permuting the ending nodes position. An example displayed in figure 5.17 shows possible permutation for the chord with three ending nodes. The numbers stand for identifiers of the ending nodes. In practice, the initial bifurcation node and the added ones have the same coordinates at the beginning of the optimization process, so the gap is introduced in figure 5.17 for the illustration purposes.

Candidate solution set for the topology optimization

The optimization problem defined in section 5.1 is solved for each feasible topology configuration. Since the branching nodes are located at the same coordinate, at the initialization of the optimization process, the list of the candidate points L 1 is defined only for the initial branching point and is duplicated for the others. The candidate points are chosen from a volume enclosed by a sub-mesh centered at the initial bifurcation point. The sub-mesh is obtained by extracting the zone limited by two planes as shown in light blue in figure 5. 18(a). The planes are given by the branch direction (in gray in figure 5.18(a)) and two points (in blue in figure 5.18(a)) situated at a distance large enough to capture possible solutions and small enough so that points that are too distanced from the split area are not taken into account. In practice, the upper limit is defined by the half of projection of the shortest branching segment (in magenta in figure 5.18(a)) on the branch direction. The lower plane is defined by a point located at a distance of half of the base segment length (in green in figure 5.18(a)) from the bifurcation point. Both limits allowing to capture a zone of adequate size were determined empirically.

The sub-mesh inner points are quite dense depending on the image resolution, shown in blue in figure 5.18(b), which results in a high total number of hypothesis to explore equal to |L 1 | n for n bifurcation nodes. For the sake of feasible computation time, the size of L 1 must be reduced. This can be achieved by subsampling the points space with a cubic grid and taking the centroids of the grid elements. The grid size is chosen depending on the largest chordae in the data, we fix it to 1 2 • R max (see figure 5.18(b)). The candidate points after the decimation are displayed in green.

For every topology configuration, the set of the coordinates of the bifurcation nodes is obtained with the optimization procedure defined in equation (5.3) (only image-based term is minimized). Afterwards, the topology configurations are compared by the value of the image cost. An example, shown in figure 5.19, illustrates the optimization results for all three possible configurations for a given chord. The rightmost configuration allows to minimize image-based cost function, therefore it is retained as the best-fitting topology configuration. This solution corresponds to the optimization results shown in figure 5.13 in green. The topology obtained with the proposed approach gives a more realistic representation of chordae according to the image. In general, the topological configurations obtained with the proposed method match the desired outcome, see the examples of the chordae with corrected topology in figures 5.20(a) and 5.20(b). Nevertheless, there are ambiguous cases where the desired modification is hard to define (20 % of faulty topology cases). Examples of such chords are displayed in figures 5.21(a) and 5.21(b). For both cases the bifurcation zones are located too close to one another, which makes it difficult to identify the splitting order. And yet the obtained solutions are perceptually correct as they are coherent with the image information and chordae anatomical characteristics.

Once all the false topology cases are corrected, the position of the bifurcation nodes is adjusted globally (per chord) in order to minimize the slack by the optimization procedure detailed in the next section.

Optimization of chordae branching nodes position

The position of the branching nodes is adjusted by solving the optimization problem defined in equation (5.3). The approach is similar to the faulty topology optimization. However, for the nodes position optimization we propose a multi-scale strategy comprising two steps. At the first step the candidate solution list is constructed on the coarsely decimated points of the sub-meshes centered at the bifurcation nodes. The second step is a refined search around the solutions resulting from the first step within a smaller zone with higher precision. Both steps are detailed hereafter.

Unlike the topology optimization case, the bifurcation points are quite far apart, therefore the list of candidate points must be obtained for each branching node. We also can reduce the size of the sub-mesh enclosing the candidate points, because the branching zones are not overlapping, like in the topology optimization case. Thus, the higher plane limiting the sub-mesh (in light blue in figure 5.22(a)) is located at a distance of quarter of projection of the smallest branching segment (in magenta in figure 5.22(a)) on the the branch direction (in gray in figure 5.22(a)). The lower plane is at a distance of half of the base segment (in green in figure 5.22(a)).

Like in the topology optimization case, points enclosed in the obtained sub-mesh (in blue in figure 5.22(b)) are too dense. To reduce the number of candidate points they are projected onto the branch direction and only those being at the distance closer than 1 2 R max are kept (in gray in figure 5.22(b)). This pruning step is performed since the desired solutions are located closely to the branching nodes and thus are not omitted. Afterwards, a cubic grid of size equal to 1 2 •R max (same as for the topology optimization) is used for the points downsampling. The points enclosed by a cubic cell are replaced with their centroids resulting in a smaller set of candidate points shown in green in figure 5.22(b). Combinations of the elements of the candidate points list obtained for each bifurcation node define the set of candidate solutions to be explored.

The candidate solutions are evaluated iteratively with the cost functions: firstly with the straightness C straight and secondly with the image-based C im . The one minimizing the objective functions as defined in section 5.1 is retained.

An example of the result obtained after the first step of the optimization (coarse solution) is displayed in figure 5.23(a). When compared to the initial configuration in black the solution in green is more straightened and agrees better with the chordae limits. The value of C straight is minimized to 4.5 comparing to the initial 6.3 and the value of C im decreased to 1495 comparing to the initial 4186.

However, the size of the grid used for the decimation is quite large, which may cause missing the optimal solution. Thus, a second step consisting in the examination of the points around the coarse solution with a finer grid is performed. The candidate points for the second refining step are pruned by distance points enclosed by a sub-mesh (see gray points in figure 5.24) limited by a cube centered at the first-step solution (green point in the center). The cube size is equal to the grid size used at the first step 1 2 R max . Gray points in figure 5.23(a) correspond to the candidate points for the second optimization step. As at the first step, the candidate points are downsampled with a cubic grid of the size equal to 1 4 R max (see cells in magenta in figure 5.24). Centroids of the grid elements (in magenta in figures 5.24 and 5.23(a)) are finally used for the second step of the optimization.

An example of the optimization result is shown in figure 5.23(c) in magenta. The improvements are visible comparing to the first optimization procedure step: the solution is more straightened (C straight is minimized to 1.96 from 4.6) and it is closer to the chordae skeleton (C im is minimized to 1295 from 1495). This allows us to consider the multi-scale procedure as an effective way of finding the position of branching nodes resulting in chordae configurations that are less slack and more centered in the chordae.

Results

The chordal trees obtained after adjustment with the proposed method for data 1 are displayed in figure 5.26. The results are shown only for the chordae subjected to the optimization, i.e. having bifurcation nodes. For most chords, the optimization outcome meets the requirements: only the bifurcations are present in the data (topology corrected in the chord in figure 5.26(a)), and the chordae are straightened and are closer to the skeleton.

The results obtained for the other data examples (data 2, data 3, data 4 and data 5) are shown in figures 5.27, 5.28, 5.29 and 5.30, respectively (see video link for 3D visualization). Similarly, significant improvement is observed for the majority of the chordae: see figures 5.27(a) and 5.27(c)-5.27(g) in data 2, figures 5.28(b)-5.28(m) in data 3, figures 5.29(a)-5.29(g), 5.29(i)-5.29(s) in data 4, figures 5.30(b), 5.30(c), 5.30(e)-5.30(h) in data 5. For these chordae the topology is correct and the chordae slack is reduced while the segments are put closer to the chordae center. The absence of a better solution for the chords in figures 5.27(b) and 5.30(a) is explained by the insufficient length of the base segment extracted in the previous stage (chapter 4). The configurations of the chordae shown in figures 5.28(a) in data 3 and 5.29(h) in data 4 derive from the incorrect detection of the chord base resulting as well from the trees tracking with modelfitting approach (chapter 4). The slack chordae in 5.30(d) derive from the topology optimization step. The topology was corrected, however the discretization for the branching nodes hypothesis generation was not fine enough, which resulted in a solution being too far from a desired one. For some chordae no solutions better than the initial ones were found with the branching node position optimization. These chordae are displayed in figure 5.31. It can be seen that the tree segments for these cases are sufficiently straightened before the optimization. For this reason a better agreement with the image was not achieved, since such solutions would not allow to minimize the initial straightness cost.

Given that the larger chordae branching have more significant impact on the global behavior of the chordal system, we experimented with another strategy where prioritizing the slack minimization in such branching can be achieved by weighting the straightness cost C straight . Each angle θ i (b) could be weighted by the relative size of the i-th branching zone S br i /S (diagonal of the bounding box), where S is the diagonal of the whole chord bounding box in the initial state, see an example in figure 5.25. The straightness cost C straight defined in equation (5.1) could be replaced by:

C straight (b) = n i=1 θ i (b) • w i , w i = S br i /S, (5.7)
where w i is the weight of the straightness angle θ i corresponding to a i-th branching. Nevertheless, such a modification does not allow to obtain better results in terms of slack reducing and image-cost value minimization at the same time. The results obtained for the chord in figure 5.32(a) illustrate the better correspondence to the image and more slack for the unweighted straightness angle sum case (chord in figure 5.32(a)). Oppositely, more straight- ened configuration in worse agreement with the chordae limits is obtained for the unweighted straightness angle sum case for the chord in figure 5.32(b). For this reason, the unweighted sum approach is kept as the final solution.

Robustness to the parameters

The only parameter used in the optimization process is the maximum chordae size R max . In order to verify the robustness of the procedure to the R max variation, we perform the tests with values R max ± 20 % on 5 CT images. No significant difference is observed for the optimization results obtained with variation of maximum chordae size value, examples are displayed in figure 5.33. In order to globally assess the parameter variation influence, the Euclidean distance was calculated between the bifurcation nodes coordinates obtained with the initial R max and with smaller value 0.8R max and, similarly, for the bigger value 1.2R max . The distances were normalized by the initial R max . The boxplots of normalized distance values are shown for two cases in figures 5.34 per data. The median of the normalized distances is close to 0.5R max , except for data 4. This can be explained by the fact that in this example the branching density is more important, thus the variation of chordae cross-sections within the data is higher. We can therefore conclude that our optimization process is not very sensitive to the change of the parameter R max .

Conclusion

An effective approach for the adjustment of the chordae geometry representation in form of trees was proposed. Both anatomical knowledge (topological characteristics, angles in the chordae branching) and image data (chordae limits) are considered in the proposed optimization procedure. The method allows to make the chordal structure consistent with the requirements of the biomechanical model by reducing the slack and to correct the faulty topologies while remaining in agreement with image information. We believe that the application of this optimization step is essential for further simulations with the biomechanical model (see chapter 2) since this results in a more realistic chordae representation and thus provides better base for the simulation outcomes.

The optimization of the chordae set geometry is performed automatically with no manual intervention. The only parameter used is the maximum chordae size R max . The robustness of the method was demonstrated for variation of the parameter values in the interval of ±20%.

For the majority of chordae, there is a significant decrease in the image and straightness cost values and the improvements are visually perceived after the application of the proposed optimization procedure. An image-based procedure for quantitative evaluation of the results will be presented in the next chapter. were extracted using a classical thinning algorithm [START_REF] Lee | Building skeleton models via 3-d medial surface axis thinning algorithms[END_REF]. This was possible because the segmented meshes are disconnected and watertight and do not contain other structures than chordae, contrary to meshes in our automatic procedure (see chapter 4). Centerline extraction is applied to facilitate the final segmentation step, at which the expert manually selected terminal points (shown in yellow in figure 6.1(b)) and bifurcation points (shown in magenta in figure 6.1(b)). Once connected, the selected points form a chordal tree (in black in figure 6.1(b)). The obtained GT is not a classical manual segmentation: the informed user performing the segmentation was aware of the data acquisition conditions (the chordae are straightened during the valve closure) and the application in the biomechanical simulation. This influenced the choice of the bifurcation nodes position in a way that resulting chordal segments have as less slack as possible. However, the manual segmentation depends on the informed user experience and on the variability in the chordae shapes. This may result in differences in nodal positions, while the chordae topology remains more robust.

a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) (q) (r) (s)

Graph similarity metric

Our goal is to design a metric that would allow to find the actual match for a given GT graph in the set of automatically segmented chordae. If the correspondence exists, the quality of the segmentation of the matched chordae is measured with the metric.

The number of methods for graph matching studied in the literature is vast. A group of approaches such as [START_REF] Leordeanu | A spectral technique for correspondence problems using pairwise constraints[END_REF] and [START_REF] Cour | Balanced graph matching[END_REF] are based on the maximization of similarity between two sets of elements with pairwise constraints. Other methods are based on GED as the one proposed in [START_REF] Abu-Aisheh | An Exact Graph Edit Distance Algorithm for Solving Pattern Recognition Problems[END_REF]. The methods have been studied and tested on our data. We will see that the latter fits better our context and therefore is applied later on.

Let us introduce the notations first. The graph of a given ground truth chord is denoted as G 1 and a graph from the automatically segmented set as G 2 . Each graph is represented with a list of vertices V and a list of edges E, thereby

|V 1 | = n 1 , |E 1 | = m 1 and |V 2 | = n 2 , |E 2 | = m 2 .

Similarity maximization methods

The methods based on similarity maximization [START_REF] Leordeanu | A spectral technique for correspondence problems using pairwise constraints[END_REF], [START_REF] Cour | Balanced graph matching[END_REF] allow to find such a cluster C of correspondences between elements of two sets that maximizes the inter-cluster score:

S = a,b∈C M (a, b), (6.1)
where a = (i, i ) and b = (j, j ) are the correspondences between elements i, j of the first set with the elements i , j of the second set. M is the affinity matrix that shows how well a pairwise relationship is preserved for each pair of correspondences. When applied to our case, the two sets of elements are the nodes of two graphs, therefore i, j ∈ V 1 and i , j ∈ V 2 . The pairwise relationship is the similarity between the edges ij and i j . The number of possible pairs is n 1 × n 2 , and the size of M is (n 1 × n 2 , n 1 × n 2 ). We define the affinity measure M (a, b) based on the length difference ∆l and the distance between the midpoints ∆d of the edges i -j and i -j corresponding to the pairs a and b as follows:

M (a, b) = 1 2 (e -∆l + e -∆d ). (6.2)
The values of M belong to ∈ [0, 1], the higher values of M correspond to more plausible matches.

The tests performed on our dataset show that similar graphs can be correctly matched as shown in figure 6.2(b). However, it is difficult to define the affinity measure for ambiguous cases where a significant part of segments is missing, as the one shown in figure 6.2(a). In this case, the vertices are wrongly matched and the edge connectivity is not preserved. Another drawback of the method is that it only allows to get the list of matches but not to quantify the quality of the graph matching. In this regard, an approach based on graph edit distance (GED) is more suitable for the problem since it allows to measure how far the matched graphs are one from another.

Graph edit distance method

We used the GED algorithm proposed in [START_REF] Abu-Aisheh | An Exact Graph Edit Distance Algorithm for Solving Pattern Recognition Problems[END_REF] and adapted it to our problem, as it was done with the similarity maximization methods. The value of GED indicates the level of proximity of two graphs. It is calculated as the sum of the costs of the modification operations required to align the two graphs. The operations include substitution, insertion and deletion of vertices and edges.

Given these two graphs G 1 and G 2 , the GED is defined as in [START_REF] Abu-Aisheh | An Exact Graph Edit Distance Algorithm for Solving Pattern Recognition Problems[END_REF]: where c(e i ) is the cost of the edit operation e i , i ∈ [1, k], k is the number of operations performed to reach the graphs alignment and γ(G 1 , G 2 ) is the set of edit paths allowing to align G 1 with G 2 . The cost of insertion and deletion operations for both vertices and edges is constant and fixed to 1 for easier interpretation. In order to define the substitution cost for the vertices of two given graphs G 1 and G 2 , the possibility of matching the i th vertex of the graph G 1 (i ∈ [1, n 1 ]) with j th vertex of the graph G 2 (j ∈ [1, n 2 ]) is first determined: if the alignment is possible the substitution cost is 0, otherwise it is 1. In our case, vertices are 3D points, so the Euclidean distance is chosen as the matching criterion. For the distance threshold value, we have to take into account possible imprecision in the location of the tree terminal nodes. This is illustrated in figure 6.3, where the displacement of the bifurcation node (orange to blue) is more significant for the larger chordae (on the left). Therefore, the alignment of vertices is based on a threshold fixed to 1 4 • S, where S is the bounding box diagonal of the chordal tree (see in figure 6.3). This allows to take the chordae size into account and to be more permissive for the longer chordae where the optimization procedure may induce larger displacements, while not accepting the correspondences between distant vertices.

GED(G 1 , G 2 ) = min e 1 ,...,e k ∈γ(G 1 ,G 2 ) k i=1 c(e i ), (6.3)
To define similarity of edges both spatial and orientation proximity are considered. The spatial proximity of two edges D es is measured by the mean of the mutually dropped perpendiculars h 12 and h 21 from the middle points (as shown in figure 6.4). The angular difference between the edges directions indicates the orientation proximity. In our procedure, the k th edge of the graph G 1 (k ∈ [1, m 1 ]) can be matched with l th edge of the graph G 2 (l ∈ [1, m 2 ]) if the distance D es between edges is less that 3 4 • R max (R max is the largest chordae in the data), and the angle difference between edges does not exceed 15 • . The thresholds were estimated based on anatomical knowledge ([KDB + 17]), tested on one image (data 1) and validated on the remaining examples in the dataset.

The graphs vary in complexity. This means that graphs with more branchings have higher The results obtained by the model-fitting approach (chapter 4) may contain false detections. In this respect, we match the GT graphs with the automatically segmented ones using the proposed metric. The chosen correspondence is the one having the minimum GED norm as it Since the GED-based metric has more advantages comparing to the similarity minimization methods, we apply it for the validation of the segmentation results presented in the next section.

Validation results

The segmentation method was tested on 5 micro-CT images using the GT provided by the first informed user. Some chordae are ambiguous for the manual segmentation due to entangling between them (see an example in figure 6.6(a)), blood clots (see an example in figure 6.6(b)), adhesion to the ventricle wall (see an example in figure 6.6(c)) and the image artifacts. In some cases, the topology could not be identified by the informed user with high certainty. Such examples were not taken into account for the validation. The number of total chordae identified and rejected due to these issues is given in table 6.1. We must note that in data 2 and data 5 not all the chordae were segmented because of a high level of ambiguity.

A total set of 74 chordae was examined. The distribution of the GED norm values is shown in figure 6.8. For 40 chordae (54.1%), the actual match was found. Examples of such chordae are shown in figure 6.9. Low values of GED norm (less than 0.4) correspond to slight inaccuracies such as shown in figures 6.10(a) and 6.10(b). Higher values of GED norm (between 0.4 and 0.6) indicate the absence of significant components in the extracted chordae (see an example in figure 6.10(c)). Nevertheless, in all mentioned cases (GED norm ≤ 0.6) the GT chordae were matched with the actual correspondences in the segmentation results which is 85.2% of tested chordae.

For the correctly matched GT chordae the image-driven C im and straightness costs C straight (see chapter 5) were evaluated in order to be compared with the corresponding costs of optimized chordae. The difference in values of image-driven costs C Res im -C GT im and straightness costs C Res straight -C GT im between the segmented optimized chordae and the GT chordae is displayed in figures 6.7(a) and 6.7(b) (y-axis) per chord (x-axis). The blue part in figure 6.7 corresponds to lower values for the automatic segmentation, the red to higher values for the manual segmenta- tion. For 92.3% of the cases the optimization allows to obtain lower image-driven costs, and for 89.7% of chordae the straightness term is decreased. The chordae set optimized with the method detailed in chapter 5 is more reliable by both image-driven and straightness costs. This shows the coherence of the proposed validation metric and the efficiency of the proposed method of the chordae optimization.

The interval [0.6-1] of GED norm corresponds to false detections (8 in total) and non-detections (GED norm is 1 for 3 cases). This is explained by the absence of real correspondences in the segmented chordae resulting either from non-detection at the image contours extraction step 10 out of 11 chordae) or due to the false rejection by the model-fitting algorithm (1 out of 11) in case of entangled chordae as displayed in figure 6.6(a).

The proposed metric is discriminative as it allows to find one-to-one correspondence between GT and automatically segmented chordae. Moreover, normalized GED values can be used to evaluate the quality of the segmentation and classify the cases of segmentation.

We note that the proposed metric is designed to distinguish the graphs primarily by the topology. This may lead to the over-penalization of the topology inaccuracy. Examples of such cases are shown in figure 6.10(d) and figure 6.10(e). For both chordae the geometry was globally correctly detected, even if the GED norm values are quite high (0.7 and 0.57, respectively). We discuss this later in the limitations of the method.

Conclusion

The quality of the chordae extraction has been evaluated with a metric based on graph-edit distances. By the means of proposed validation procedure, the completeness of the extracted chord has been assessed along with its topology. The efficiency of our chordae extraction approach has been evaluated with the metric on 5 micro-CT images: 85.2% of tested chordae were correctly extracted with our automatic approach.

A limitation of the method is the over-penalization of the topology inaccuracy. This issue could be addressed by giving more importance to the segments being longer and closer to the tree root.

Comparing the segmentation obtained with our method with the manually extracted chordae is not sufficient as the correspondence does not guaranty that the chordae would maintain the valve competence. A biomechanical simulation is a way to verify if the valve behavior is coherent, which is addressed in the next chapter 7.

Chapter 7

Biomechanical simulation of the valve closed state

The objective of the study is to reproduce the valve behavior at peak systole via biomechanical simulation. The chordae geometry representation resulting from the proposed automatic segmentation procedure is used in the biomechanical model to simulate the closed state of the valve.

In this chapter, stability and coherence of obtained simulation outcomes is checked. We compare the outcomes corresponding to the automatic and manual chordae segmentations while verifying that the proposed segmentation procedure does not induce leaflet anomalies such as bulging or leaking. The consistency of the reproduced valve behavior is also analyzed for a modeled pathological case (ruptured chordae). As quantitative criteria we use the mechanical stress distribution (von Mises criterion) and the bulging volume.

Valve geometry integration into a biomechanical model

For the simulation we use the biomechanical mitral valve model detailed in chapter 2. The method for the automatic segmentation of the mitral valve components was proposed only for the chordae with the methods proposed in previous chapters. Therefore, prior to the simulation, the leaflets surface was obtained via manual segmentation and afterwards, the chordae were connected to them.

Leaflet segmentation

The segmentation of the leaflet surface is a challenging task since the images are acquired in the closed valve state. As it was mentioned in chapter 2, the definition of the coaptation zone is difficult because the leaflets are sealed and there is no gap between them as it can be seen in figure 7.1(a). Moreover, numerous folds are present on the leaflet surface (see figure 2.15(b)).

To address this task, the 3D leaflet contours were extracted by the informed user in ITK Snap software [YPCH + 06]. Knowing that the leaflets must be separated when used in the biomechanical model, a gap (in yellow in figure 7.1) across the coaptation line was introduced by the informed user allowing to detach the anterior (in red in figure 7.1) and posterior (in green in figure 7.1) leaflets while preserving the tissue near the commissures (see section 2.3). After that, the anterior and posterior leaflet labels were merged and transformed into a mesh with the Marching cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] as shown in figure 7.2. The segmentation procedure is quite time consuming as the level of uncertainty of the coaptation zone position is high and thus the segmentation requires a lot of adjustments. In the biomechanical model described in chapter 2, the leaflet is modeled as a surface with constant thickness. Therefore, only the upper surface is used in the simulation. It was manually isolated using MeshLab software [CCC + 08]. High data resolution and the application of the Marching cube algorithm result in excessive number of vertices (≈ 90k) in the mesh, which can slow down the simulation and is not necessary for the model precision as it was demonstrated in [START_REF] Hammer | Simulating heart valve mechanical behavior for planning surgical repair[END_REF]. In order to reduce the number of vertices, the mesh was decimated using the Vorpaline software [START_REF] Levy | Variational anisotropic surface meshing with voronoi parallel linear enumeration[END_REF]. The final leaflet mesh contains 591 vertices and is used in the biomechanical model afterwards (see figure 7.2(b)).

Chordae connection to the leaflets

There is a gap between the leaflet surface and the chordal trees when they are integrated into the model. The gap corresponds to the zone where the chordae connect to the leaflets (see lower part of the leaflet mesh in figure 7.2(a)). In this area the chordae intermingle with the very leaflet substance [START_REF] Muresian | The clinical anatomy of the mitral valve[END_REF] (see in figure 7.3(a)) and smoothly transition to the leaflet surface. This sets a limit for the model-based chordae segmentation method (chapter 4) allowing to extract the chordal structure until certain level, which is still far from the upper leaflet surface (see the gap in figure 7.4).

There are several possible approaches for the chordae connection to the leaflet surface. In the models where the chordae geometry is not image-based, such as [START_REF] Prot | Finite element analysis of the mitral apparatus: Annulus shape effect and chordal force distribution[END_REF], the chordal network is often represented with fan-shaped structures originating at papillary muscles as displayed in figure 7.5(a). In the models that use chordae geometry segmented from images acquired in a controlled environment (e.g. in [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF]), the most common method consists in connecting the chordae extremities to the leaflet surface within a sphere of fixed radius (see in red in . We use an image-based approach for the chordae geometry representation, and hence consider two methods for chordae prolongation: fan-shaped and one-point connection. Fan-shaped connection approach makes the system more constrained and could result in an artificially more stable simulation outcomes depending on the fan radius. On the contrary to the study [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF], we extracted the chordae automatically. The chordae are located at different distances from the leaflets, therefore applying the constant radius approach is not suitable for our data. The chordae connection to the leaflets is a complex task, solution for which is neither straightforward even for an experienced user.

Before describing our strategy for the chordae connection to the leaflets, we first introduce a global constraint C orient between the general direction of a chord and the orientation of a local connection point that we use to discard erroneous connections.

The chordae connection must be performed considering the general chordae orientation, because the extremities of a chord have to connect to the same leaflet (see figure 7.6). Otherwise, they will prevent the valve from opening. Only the chordae near the commissures have a cleft shape and can be attached to both leaflets (see figure 7.3(b)). For each chord, the general direction is defined (in cyan in figure 7.7) as a vector connecting its origin point and the mean position of its terminal nodes. Given that the surfaces of the anterior and posterior leaflets face each other in the coaptation zone, we use the dot product sign in order to define the position of the chordae with respect to the leaflet node. More precisely, the dot product between the normal n c to the point where the mean chordae direction intersects the leaflet (in cyan in figure 7.8(b)) and the normal n con (in blue in figure 7.8(b)) to the point where the chordae is attached must be positive. In this case, the constraint C orient is satisfied. It must be noted that this process is sensible to the segmented chordae orientation, which may result in incorrect attachment close to the coaptation surface. We consider two approaches: the connection by a fan-shaped structure to several leaflet points and the one-point connection to the leaflet. To illustrate both approaches the terminal chord point is shown in green, the leaflet point being the nearest to it is displayed in magenta and denoted p nr . The supporting segment of the terminal point determines the vector for the chordae prolongation. The intersection point of the leaflet surface with the supporting vector is shown in blue and denoted p int .

Fan-shaped structure connection Using a fan-shaped structure is a way to represent the missing chordae in the zone where the chordae are connected to the leaflets. In order to define the area of connection two points are used: the intersection point p int and the nearest point p nr . According to the position of these points, we separate two cases for the chordae connection.

In the first case, the terminal point (see figure 7.8(a)) is located close to the leaflet surface, and p int and p nr are more or less in the same point. The connection is performed to points being enclosed by a sphere of a radius d c as shown in dash black in figure 7.8(a). The radius d c is In the second case, the intersection point p int and the nearest leaflet point p nr (see figure 7.9(a)) are located far from each other. The connection is then done to the leaflet points (see in dash black in figure 7.9(b)) located in a cone defined by the bisector between the vectors connecting the terminal point with p int and p nr .

In both cases the dot product sign between the normal n c (defined by the mean chordae direction, shown in cyan in figures 7.8(a) and 7.9(a)) and the normal to a possible connection point n con (displayed in blue in figures 7.8(b) and 7.9(b)) is verified. The connection is performed if the dot product is positive (see in red in figures 7.8(b) and 7.9(b)).

One-point connection

The second type of connection to one leaflet point is performed as follows. Each terminal point (see figure 7.10) is projected following the direction given by its supporting segment on the surface of the leaflet. If the intersection with the leaflets surface exists and the constraint C orient is satisfied, the terminal point is replaced by the surface point as shown in figure 7.10(a). In the case when the intersection with the leaflet cannot be found, nearest points are considered: they are checked one by one, and the connection is performed to the point that satisfies the constraint C orient (dot product sign is positive as shown in figure 7.11(b) in magenta).

Biomechanical simulation

The chordae connection obtained with the two proposed methods is shown in figure 7.11. As the choice of the best suited technique is not straightforward, tests via the biomechanical simulation were conducted. The goal is to find the chordae connection type minimizing the simulation instability without adding unreasonable number of supplementary chordae that will over constrain the system. We propose to analyze the convergence curves of the node velocities (detailed in chapter 2). In order to evaluate the whole system, the mean of z-component of the nodes velocity and the mean of the nodes velocity norm are checked. The simulation is considered stable when few oscillations appear.

With the fan-shaped and the one-point connection the node velocity components have slightly different convergence speed. The fan-shaped connection converges to lower z-component velocity values: 0.09•10 -3 against 1.44•10 -3 (see figure 7.12). This is due to the high influence of this type of connection on the model, as the more segments are used for the connection the more the system is constrained. This reduces the degrees of liberty of the system and makes it less realistic. The fan-shaped connection causes higher instability as higher velocity norm values are observed (see figure 7.12) which indicates more movement in x-and y-direction. This can be explained by the slack added with this type of attachment in typical cases as the one shown in figure 7. 13(a).

In view of these considerations, the one-point connection is used for further simulations because we privilege simulation stability.

The extracted chordae set is validated via biomechanical simulations. First, the coherence of the valve behavior with the expected outcome is verified for the complete chordal set case.

In such conditions, the valve should remain sealed and the leaflet bulging has to be minimal. The influence of the chordae optimization is studied as well. Secondly, a pathological valve case is modeled and simulated. The validity of the chordae segmentation may be evaluated by appearance of leaks and leaflet bulging.

Validation with full chordae set

Three chordae configurations have been tested for the experiments: the GT chordae resulting from the manual segmentation (see figure 7.14(a)), the chordae obtained from our complete pipeline (see figure 7.14(b)) and the same configuration but with enforced slack. The slack configuration was added since it is difficult to compare the configurations with little differences in the attachment points between manually segmented chordae and the optimized results. Therefore, to analyze the influence of slack on the simulation outcome, we decided to induce it into the optimized chordae by uniformly elongating them by 10% in the biomechanical model. It has to be noted that some chordae were not detected with our method (the four chordae shown in green in figure 7.15). As it can be seen from the points of their attachment to the leaflets, these are second-order chordae whose impact on the valve behavior should be restricted (see introduction section 1 for the explanation).

Leaflet pre-strain

The resting state of the leaflet is the open valve state, so the stresses and the strains are null. The open state corresponds to the green point on the plot of the constitutive law (equation (2.4)) displayed in figure 7.16. At this state, the leaflet is very extensible as it can be seen with the curve slopes. However, the simulation is not starting from this stress state as the leaflet is in a closed configuration and its behavior is locally becoming stiffer. A stress state S pr should be estimated to be used in the simulation as a way to "pre-strain" the leaflets.

To roughly estimate this initial stress state, we use the Laplace law which describes the pressure-volume relationships of spheres (equation (7.1)) where P is the pressure, h is the sphere thickness and r is the sphere radius and T is the tension. It is classically used to understand physiological features in cardiac analysis [START_REF] Thomas | The law of laplace its limitations as a relation for diastolic pressure, volume, or wall stress of the left ventricle[END_REF]. It states that the tension within the wall of a sphere filled to a particular pressure depends on the thickness of the sphere.

P = 2 • h • T r (7.1)
The mitral valve shape at peak systole could be estimated by a thin half-sphere. During this stage, there is also a uniform blood pressure, so it fits well the Laplace law conditions. The tension term T corresponds to the stress state S pr . Equation (7.1) allows to estimate S pr as detailed in equation (7.2).

S pr = P • r 2 • h . ( 7.2)
With P = 100 mmHg, r = 1 cm and h = 1 mm, it gives S pr = 66.5 KP a.

Once the initial stresses of our simulation are estimated, the initial strains for both anisotropic directions E 11

pr and E 22 pr are found by using the constitutive law (blue and orange points on the plot displayed in figure 7.16). The resulting values are E 11 pr = 0.1974 and E 22 pr = 0.2098. In our biomechanical simulation, the strains E comp are computed from the vertex displacements, then the stresses S comp are computed with the constitutive law and finally the internal forces are computed from the stresses. In order to take these initial strains into account, the strain states E curSt are computed as follows:

E 11 curSt = E 11 pr + E 11 comp , E 22 curSt = E 22 pr + E 22 comp . (7.
3)

The stresses are then computed from E curSt and not from E comp as initially.

After experimenting, we noticed that the addition of the initial strains is important for the simulation to ensure the valve closure. Otherwise, we observed significant leaflet bulging, since without the pre-strain the leaflets are too distensible. This is illustrated in figure 7.17 where the profiles of the leaflets with and without pre-strain conditions are compared at the end of the simulation (case of the automatically extracted chordae set).

Quantitative metrics

In order to quantify the results we compare the bulging volume. By the bulging volume we mean the volume between the leaflet surface at the beginning and at the end of the simulation (see two surfaces in figure 7.18). The volume is computed per face as volume of four tetrahedrons that can be constructed on the face vertices in the initial (

V 1 , V 2 , V 3 ) and final (V 1 , V 2 , V 3 ) states (see figure 7.19): 1 6 (A -D) • ((B -D) ∧ (C -D)), (7.4) where A, B, C, D ∈ {V 1 , V 2 , V 3 , V 1 , V 2 , V 3
} are the vertices of the tetrahedron. Another metric used for comparison of the simulation outcomes is the von Mises stress distribution σ V on the leaflets. The von Mises stress criterion σ V is a measure of the global stress acting on a material. Since the stress is biaxial, the σ V is calculated from the Cauchy stress Simulation results First, the leaflet profiles obtained in the simulation with the three chordae configurations were compared (see figure 7.20) in the cut plane shown in figure 7.17(a). As expected, the configuration with induced slack results in a more significant leaflet bulging, comparing to the cases of manually segmented chordae and the chordae obtained with our method, which have similar leaflet profiles.

The corresponding bulging volumes were measured: 0.385 cm 3 for the GT chordae, 0.411 cm 3 for the automatically segmented chordae and 0.803 cm 3 for the chordae with induced slack. Even if porcine and human anatomic measurements are not fully comparable, we obtained volumes in the range of the values found in the literature, for example, in the study of MV morphology in children and young adults with Marfan syndrome [JHG + 18]. As is customary in cardiology, the values are normalized by body square area (1.4 m 2 ) and are 1.07 cm 3 /m 2 for patients with MV regurgitation, 0.13 cm 3 /m 2 for patients with MV prolapse. The corresponding not normalized values are 1.498 cm 3 and 0.182 cm 3 , respectively.

The map of bulging volume differences calculated per face (not including the coaptation zone) between the optimized chordae set and the manually segmented chordae is shown in figure 7.21(a), while the map of differences between the chordae set with induced slack and the manually segmented chordae is in figure 7.21(b). It can be seen that considerably more leaflet bulging is observed for the slack chordae. The zone were the difference between the volumes obtained with the automatically segmented chordae and the manually segmented set corresponds to the location of non-detected chordae i.e. false negatives of the segmentation algorithm (see in green in figure 7.15). However, as these chordae are second-order, their influence is limited. This explains the small differences in the bulging volumes and in the leaflet profiles regardless of their absence in the automatic segmented chordae set.

We must note that the bulging volume in the case of complete chordal network for both manually and automatically segmented sets is quite high. This may result from inaccuracies and errors in the segmented chordal sets and from inappropriate constitutive law that does not correctly describe the real valve material properties.

Finally, the stress distributions obtained for the three cases were analyzed (see figure 7.22). The results are quite similar to those found in the literature and are in similar ranges of values: maximum of 3.5 •10 5 Pa in [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF], maximum of 1 •10 6 Pa in [TJE + 16]. The difference map between von Mises stress obtained with the manually segmented chordae and the optimized chordae set in displayed in figure 7.23(a) and with the chordae set with induced slack in figure7.23(b). Higher stress value differences are observed for the chordae configuration with slack, while for the automatically segmented set the differences are close to zero except for the zones of the chordae attachment (see figure 7.15(a)).

The analysis of the simulation outcomes shows that the chordae architecture resulting from our algorithm can give results similar to experienced users while providing an equivalent biomechanical simulation. This means that the chordae obtained with our method allow to spread the stress similarly to those segmented manually, which is not the case for the slack chordae. In this case, lower stress values and more important leaflet bulging is observed, which proves that straightening of the chordae can improve the simulation outcome.

The simulation convergence was analyzed for the case of the chordae extracted automatically and optimized with the procedure proposed in chapter 5. The curves of mean node velocity z-component and mean node velocity norm are displayed in figure 7.24. At a coarse scale, the curves converge to small but non-zero values (1.19•10 -3 and 1.57•10 -3 ). However, at a smaller scale, oscillations prevent from adopting a criterion for automatically stopping the simulation.

Validation on pathological case

The simulation of a mimicked pathological case would allow us to verify that the behavior of the valve is coherent and the leaflet will not remain sealed. In order to model a pathological valve with ruptured chordae, some chordae were deleted from the sets obtained via manual and automatic segmentations as shown in red in figure 7.25. We expect that applying a constant pressure will lead to valve leakage.

We perform the experiment with the GT chordae and the chordal set obtained with the proposed segmentation pipeline. The leaflets profiles correspond to the expected outcome. As is shown in figure 7.26, the chordae deletion leads to a leak in the area where chordae were ruptured in both cases (in red in figure 7.26) while the leaflets remain sealed in the other zones. This proves the coherence of the simulated valve behavior when the chordae are extracted automatically with the proposed method.

Similarly to the simulations with the complete chordal set, in order to quantify the results, we compare the distribution of the von Mises stress on the leaflets. The stress distributions show that the simulations with manually and automatically extracted chordae lead to similar final valve configurations (see figure 7.27). Moreover, higher stresses are observed in the areas close to the zones corresponding to the leaflet leakage areas (circled black zones) when compared to the results in figures 7.22(a) and 7.22(b) for the GT and automatically segmented chordae, respectively.

Conclusion

The extracted chordal set was validated via the biomechanical simulation of the valve closed state. The application was shown only for one subject and the results of the simulation application are mainly qualitative. Regardless of that, the leaflet state at peak systole is well reproduced in the challenging context of real data.

The simulation stability was verified by analysis of the node velocities and showed the convergence towards small values close to zero.

The coherence of the segmentation was also demonstrated within the simulation of a pathological case where the chordal set is incomplete. The similarity of the simulation outcome obtained with manually segmented chordae to the one observed when the chordae are extracted with our algorithm straightens our belief that the proposed chordae extraction method can be applied in the computational valve modeling. However, bulging is observed even in the case of the full chordal network for both manually and automatically segmented chordae. We suppose it is due to the possible detection errors in both chordae sets or to the leaflet constitutive law which is not appropriate for the valve tissue in our experiments.

The proposed validation approach is limited by the ambiguity in the leaflet segmentation since the image acquisition performed only in the closed state. Usage of temporal data would allow to get more precise leaflet geometry and enable a more realistic simulation.

Conclusion

In this thesis we aimed to develop methods for realistic simulation of the mitral valve behavior using image-based geometry. To achieve this, an approach linking the biomechanical model and medical images was proposed, including methods for valve components segmentation and their integration into the model. We focused on automatic chordae extraction as it would improve the simulation reproducibility on various data examples, and thus help to advance toward patientspecific clinical applications.

Contributions

Automatic method for chordae segmentation

A multi-stage approach for automatic chordae segmentation from medical images was proposed. Given that the study implies biomechanical simulation, the method was designed in connection to the chosen mechanical model.

Topology-based segmentation: the first segmentation step is an approach exploiting topological properties of the structures, which provides coarse but satisfactory data labeling without any model of the chordae cross-section and guided only by one parameter (the largest chordae size).

The labeling performs better for data having small variations in chordae cross-sections. Large cross-section discrepancy may lead to the detection of disproportionate branching zones comprising several bifurcations of thinner chordae. However, this limitation is corrected at the following segmentation step.

Model-based chordae extraction: the second stage describes a tracking procedure based on model-fitting. We approximate a chord with an elliptical cylinder on the local level and process branching zones separately. The resulting chord structure is an oriented tree of connected line segments, where the discontinuities are corrected relying on segment collinearity.

The tree extraction method proves to be reproducible except for the zones containing untypical erroneous structures. The lower parts of the chordae are detected more accurately with the method. In the areas close to the leaflets the segmentation may be incorrect as the chordae are thin and entangled, and thus less distinguishable and not geometrically equivalent to elliptical cylinders.

The topological errors in the extracted trees are mostly due to the excessive size of the submesh with respect to the size of the branching zone. To address this issue, a multi-scale approach could be considered, which would adapt locally by enlarging or reducing a tracked sub-region. This would enable fitting of smaller cylinders and better locating branching zones, and hence refine the resulting structure.

Chordae geometry optimization: at the final step of the chordae extraction, we addressed the compatibility between the segmented chordal network and the biomechanical model. An optimization procedure with anatomical and mechanical constraints joined with image-based chordae limits effectively corrects the majority of inappropriate chordae configurations, such as faulty topology cases and slack chordae.

The generation of branching nodes hypotheses requires the tuning of parameters dependent on the largest chordae size. The robustness of the method to the parameters was shown. For data where chordae are similar in cross-section, the robustness of the approach is higher. The strategy in perspective is to adapt to the search zone for the candidate points choice.

For the sake of computation time, the optimization method uses a discrete set of candidate points for branching nodes. Another solution could be to implement parallel hypotheses processing in order to speed up the procedure and to test more hypotheses.

The slack reduction is performed for each chord independently. Though this strategy has proven to improve the quality of the prediction, only complete biomechanical simulation with the full set of chordae will allow us to properly evaluate the slack.

Graph-based validation procedure

We have designed a graph-based approach for chordae segmentation validation in the context of further biomechanical applications. The method uses graph similarity between automatically extracted and ground truth chordal trees. The approach allows to assess the segmentation quality and distinguish false detection cases.

As the chordae graphs vary in complexity, we proposed a normalized metric to compare them, and not a simple threshold-based criterion. Even if the global structure is extracted correctly, the method limitation is the over-penalization of topology inaccuracies. One of the possibilities to overcome this problem could be to weight the chordae segments depending on their length and proximity to the chordae origin.

The method was tested with the ground truth provided by one informed user. For the sake of generality, more than one user would be desirable.

Biomechanical simulation with subject-specific geometry

The results of chordae segmentation and the manually segmented leaflet geometry were integrated into the biomechanical model. A method for the chordae extension till the leaflet surface was based on a simulation stability analysis.

The closed valve state was simulated with complete and incomplete chordae sets as a way to model healthy and pathological cases. The coherence of the simulation was demonstrated by comparing the outcomes obtained with the chordae segmented with our automatic methods and the ground truth chordae. The comparison was done regarding general leaflet profile (valve competence), bulging volume, and stress distribution. The similarity of the obtained results for automatic and manual segmentations indicates that the proposed chordae extraction approach can provide a reliable input for the biomechanical model.

The simulations demonstrate the possibility to guide a mechanical model with a valve geometry extracted with the minimum of human interaction. These results are promising, however a deeper analysis has to be done. More precisely, simulation convergence criteria are not straightforward to define, which is reflected in the long time required for the simulation. Adding user interactivity may be a possible interim solution for the simulation stopping. Among the numerous parameters involved in the model, some of them are linked to unknown (or hardly measurable) specific mechanical properties and cannot be well controlled. A better grasp of parameters is one of the future challenges of this work.

The way in which the chordae connect to the leaflets may strongly impact the simulated valve behavior, therefore more robust connection methods are worth being explored.

Perspectives

Short term perspectives and software engineering

We aimed to design methods for mitral valve modeling to be as automatic as possible. However, a certain level of interaction could help to eliminate the ambiguity at some points of the segmentation algorithms, for instance, with seed initialization points near complex configurations and with badly positioned connections of the chordae to the leaflets. To this aim, improved tools for 3D interaction with the segmentation will be helpful. Even if an experimented user is not always able to perform the whole segmentation process better than the algorithm, he may, however, improve the automatic segmentation by local modifications.

In the perspective of assessing one simulation or to compare several simulations, various criteria have been used: leaflet bulging, strain map, mean velocity, local velocity. We think that a graphical interface with these features could help the clinician to understand the dynamic behavior of the mitral valve and to stop the simulation earlier.

Improved data acquisition

Even if the available dataset still can be used for further investigations along with the prepared manual segmentation, obtaining more accurate ground truth would improve the validation procedure.

The dissociation of geometry and mechanics could help to study the parameters influence. To disconnect the data uncertainty and the mechanical parameters, we plan to acquire the motion of a valve phantom with high-rate stereo-vision cameras. With such a setup, characteristics of the leaflet are perfectly known and the stereo-vision would allow us to acquire ground truth temporal reconstruction of the leaflet.

As another practical application for the future, an image modality being closer to the clinical setting may be considered, for instance temporal MRI. The progress of MRI would probably help to segment the leaflets and the chordae from images corresponding to different valve states.

Improving the modeling of the chordae

The results obtained with the proposed model-based chordae segmentation method are rather good on the local level: branching zones are identified and chordae sub-regions are well approx-imated with cylinder models. Defining whether obtained geometry precision is sufficient or if a finer segmentation is required could be tested on the global functional level. Providing a coherent global chordal structure is challenging. The current tracking procedure does not allow to reconsider errors occurring earlier. Especially segments that may be erroneously attributed to a chord cannot be re-affected when this chord is considered. An opportunity to explore is enforcing topological constraints (tree topology type) on the chordae during tracking procedure. Another possibility is to process the chordae in parallel, which would allow to correct the ambiguous cases once they are encountered. Such an approach may also require user interaction for chordae origins identification.

Our experiments have proven that our chordae extraction method achieves good performance except for the upper region of the trees where many entangled and small segments are present. A possible way to solve this problem could be to adopt a heterogeneous model of the tree where a model-based cylinder detection is kept in the bottom part of the tree whereas a representation of the region around the leaflet could be modeled as a thick material with appropriate elastic properties. The fundamental question that we want to investigate in future studies is: could we have a simplified model of the valve which gives equivalent functional results?

Real time simulation and outliers handling

An important bottleneck of the method is the time needed to perform mechanical simulations. It is not so much because of the time required by one simulation but because it prevents us from using robust strategies at the global mechanical level. Indeed, hypotheses on branching points are currently generated and solved locally. If simulation time was fast, we could imagine to keep trace of all these hypotheses and generate a set of tree configurations that are ranked based on physical simulation. However, this requires fast ways to perform the simulation. In the short term, resorting to SOFA [FDD + 12, VHP + 18], which provides fast but approximated physical behavior could be a possibility. In the longer term, another alternative could be to investigate new and promising methods based on convolutional neural networks that are able to mimic physically-based simulation while requiring much less computation time [LLMS18, MTH + 20]. This can open a way towards effective global methods for improved pruning of tree hypotheses.

Discussion on the choice of the biomechanical model

We have used a FEM-based structural model being relatively fast in comparison to FSI models and allowing quite precise simulations. Although the code is available and the mechanical parameters can be tuned, it is difficult to define their impact on the simulation. Finding an approach to study the influence of each isolated parameter is a path to investigate. This would help to identify the most influencing ones and keep a reduced set of parameters resulting in a more compact model.

Several clinical gestures like leaflet resection of a bulging valve or annuloplasty could be simulated with the proposed framework with a structural model. This is possible due to the fact that we studied the static closed valve state. Increasing the model complexity and opting for a FSI model may enable the simulation of other procedures, such as leak repair. This cannot be carried out with the applied "dry" model where blood flow is not considered. On the other hand, even if the simulated intermediate states of valve motion are not physically valid, for the cases where the objective is to reproduce the final static state, a structural model still can be a suitable option. It is worth mentioning that FSI-based models convergence is highly sensitive to mesh size and regularity. As the segmented mitral valve contains irregularities, it is not clear if our current models could adapt to FSI methods with a simple mesh densification or if FSI methods require deeper changes in segmentation to meet the needed accuracy and boundaries regularity. The choice of a mechanical model type remains complex and should be motivated by a trade-off between speed and simulation fidelity required for the intended application.

Résumé

Dans cette thèse nous avons proposé des méthodes pour la simulation réaliste du comportement de la valve mitrale en inférant sa géométrie à partir d'images. Afin de réaliser cet objectif, la segmentation des composants de la valve a été guidée par la nature du modèle biomécanique choisi.

Des méthodes pour la segmentation des composants de la valve et leurs intégration dans le modèle ont été proposées. Nous nous sommes concentrés sur l'extraction automatique des cordages, car cela permettrait d'améliorer la reproductibilité de la simulation et ainsi d'avancer vers des applications cliniques basées patient.

La procédure proposée pour l'extraction automatique des cordages comporte trois étapes. Premièrement, la segmentation des cordages utilise les propriétés topologiques des structures. Ensuite, la représentation de la géométrie sous la forme d'arbres de segments de droite connectés est extraite en utilisant une approche de type suivi couplée avec l'ajustement local d'un modèle cylindrique. Afin de rendre la géométrie obtenue conforme au modèle biomécanique, une étape finale d'optimisation de la structure des cordages est effectuée. Elle vise notamment à réduire le mou des cordages et de corriger certaines fausses configurations topologiques en prenant en compte des considérations mécaniques et anatomiques ainsi que les données d'image.

Nous avons proposé une métrique de validation des résultats de la segmentation basée sur la similarité de graphes. Elle permet de définir des correspondances entre les cordages obtenus avec notre méthode automatique et ceux provenant de la segmentation manuelle, ainsi que d'évaluer la qualité de la segmentation.

Nous avons ensuite utilisé la géométrie des cordages extraite pour la simulation biomécanique de l'état fermé de la valve. La cohérence du comportement reproduit a été vérifié en comparant les résultats obtenus avec les cordages issus de notre procédure automatique et les cordages de la vérité terrain. Le volume du gonflement des feuillets et la distribution de la contrainte mécanique ont été utilisés comme critère d'évaluation. La conformité de la simulation a été également analysée sur un cas pathologique modélisé avec des cordages rompus.

Mots-clés: valve mitrale, modélisation basée image, segmentation, simulation biomécanique.

Abstract

In this thesis we aimed to develop methods for the realistic simulation of the mitral valve behavior using image-based geometry. To achieve this goal, the segmentation of valve components was guided by the nature of the chosen biomechanical model.

Methods for the valve components segmentation and their integration into the model have been proposed. We focused on automatic chordae extraction as it would improve the simulation reproducibility on various data examples, and thus help to advance towards patient-based clinical applications.

  Introduction 1 Medical context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Thesis organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 State of the art in mitral valve modeling 1.1 Mitral valve geometry representation . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Geometry entirely based on measurements after dissection . . . . . . . . . 1.1.2 Geometry based on measurements after dissection for the chordae and image-based for the leaflets . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Geometry fully based on images . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Biomechanical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Validation protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Challenges in computational valve modeling . . . . . . . . . . . . . . . . . . . . . 2 Dataset and biomechanical mitral valve model 2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Acquisition protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Dataset analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 FEM-based mitral valve biomechanical model . . . . . . . . . . . . . . . . . . . . 2.3 Segmentation method requirements . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Validation strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Topology-based chordae segmentation 3.1 Vascular vessel segmentation methods . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Initial data processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Mesh labeling method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii 1 Cardiac cycle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Mitral valve anatomy. Open (left) and closed (right) states (taken from [CAF10]). 3 Chordae classification (taken from [Mur09]). Anterior leaflet is denoted as AoL, posterior leaflet is denoted as ML. . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Mitral valve pathologies (taken from [Tur04]). . . . . . . . . . . . . . . . . . . . . 5 Examples of mitral valve repair procedures (taken from [CWMG19]). . . . . . . . 1.1 MV imaging using different modalities. . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Valve model based on geometry from post-mortem surgical studies proposed in[KRC97], [CK98]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Valve model based on both image data and post-mortem measurements proposed in [PHS09]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Manually segmented leaflet geometries. . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Automatic leaflet segmentation approaches. . . . . . . . . . . . . . . . . . . . . . 1.6 Valve model model proposed in [WS13]. . . . . . . . . . . . . . . . . . . . . . . . 1.7 Leaflet segmentation process proposed in [TNH + 18]. . . . . . . . . . . . . . . . . 1.8 Valve geometry acquisition with an in vitro setup proposed in [Kha18]. . . . . . . 1.9 Connection of the chordae extremities to the leaflet surface with a fan-like structure proposed in [DKS17]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 Comparison of two FSI-based models (tubular-shape and ventricular-shape fluid volumes) with a structural model by leaflets coaptation zone profile in 2D in [LDSB10]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Valve stress-strain behavior in fiber and cross-fibre directions proposed in [KEC07]. 1.12 Valve stress-strain behavior for anterior (AML) and posterior (PML) leaflets in circumferential (fiber) and longitude (cross-fiber) direction. Taken from [WS13]. . 2.1 Ex vivo valve open state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Ex vivo valve closed state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Data acquisition : micro-CT device with pressurized air tube going inside and its screen (in red and zoomed on the right) displaying an X-ray of a porcine heart. . 2.4 Examples of the artifacts present in the data. . . . . . . . . . . . . . . . . . . . . 2.8 Proportion of data exploitable data in the dataset. . . . . . . . . . . . . . . . . . 2.5 Example of chordae being close to the ventricle wall. . . . . . . . . . . . . . . . . 2.6 Example of entangling chordae in the data (data 2). . . . . . . . . . . . . . . . . vii List of Figures 2.7 Example of a blood clot in the data (data 1). . . . . . . . . . . . . . . . . . . . . 2.10 Illustration of a mesh triangle deformation. . . . . . . . . . . . . . . . . . . . . .

3. 4

 4 Contours and mesh extraction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Topological characteristics of non-chordae, cylinder and bifurcation surfaces. . . . 3.6 Different types of topology. Locally extracted sub-mesh is shown in green. Central point of the sub-mesh is in red. The data mesh points are in blue. . . . . . . . . 3.7 Sub-mesh size definition based on geodesic distance. . . . . . . . . . . . . . . . . 3.8 Mesh propagation scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 11

 11 Segmentation results (data 1). Non-chordae structures are in blue, chordae are shown in red, branching zones in green and magenta. False detections are in black rectangles. Differently labeled with two methods branching zones are in orange rectangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Labeling process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Results of the chordae segmentation on 5 micro-CT images. Blue corresponds to the whole boundary mesh, the segmented chordae are in red. . . . . . . . . . . . 3.14 Superposition of the segmentation results (red contour) with the ground truth (blue inside) shown in axial, coronal and sagittal views (data 1). Yellow lines show correspondence between the slices. First line: for a point at half maximum height of the chordae. Second line: for a point near the chordae ending close to the leaflets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .viii 4.1 Examples of incorrectly detected type-1 shapes (circled in black) with a topologybased segmentation method. Chordae are in red, non-chordae structures are in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Examples of the results obtained with the skeletonization methods proposed in [CTO + 10] (ROSA 2.0) and in [TAOZ12] (part of CGAL). Black zones correspond to correct detections. Red and magenta zones indicate regions where the outcomes are not satisfactory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Examples of branchings profiles typically considered in the vascular structures segmentation methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Tree transformation. Terminal points are in yellow, branching points in magenta, intermediate points in gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Parametric cylinder model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Model fitting procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Ellipse fitting procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Example of a semi-tube structure shown in two views (circled in black). Chordae are in red, non-chordae structures are in blue. . . . . . . . . . . . . . . . . . . . . 4.9 Distribution uniformity check. Data points projected on the cylinder base are in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Example of a fitted cylinder model (in light blue) to data points (in blue). . . . . 4.11 Tree tracking procedure. Chordae mesh points are in gray, initial tracking point is in magenta, locally extracted sub-meshes are in pink, fitted cylinders are in blue, initial points for the sub-mesh extraction are in yellow, type-2 shape submesh boundary points are in green, branching zone mesh is in red, branching zone sub-mesh boundary points are in orange. See video link for 3D visualization. . . . 4.12 Chordal tree tracking procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Detection of a missed branching zone. . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Non-valid cases of a potential branching zone mesh (in red). Chordae mesh is in pink. Initial sub-mesh point in in yellow. Gravity centers of the mesh boundary connected components are in green. . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Tree tracking at a branching zone. Chordae mesh points are in gray, initial points for sub-mesh extraction are in yellow, locally extracted sub-meshes are in pink, fitted cylinders are in blue, branching zone sub-mesh is in red, branching zone sub-mesh boundary points are in red. . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Variation of the maximum number of the iterations. Each chord is in separate color. The problematic zones are circled in red. . . . . . . . . . . . . . . . . . . . 4.17 Leaflet part being a not complete tube. Fitted model is in red. Chordae mesh is in pink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Results of the tree extraction. Each chord is in separate color. . . . . . . . . . . 4.19 Results of the tree extraction for data 3 obtained with 5 iteration of the method. Each chord is in separate color. The problematic zones are circled in red. . . . . ix 4.20 Representation of the chordae in a form of connected line segments (data 1). Chordae segments are shown in black, chordae mesh is in gray. Starting points of chordal trees are in red, branching points are in magenta and terminal points are in yellow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Example of discontinuities in the trees. Magenta and blue zones correspond to the tree elements preceding and following the discontinuities, respectively. Red segments are the connected centerlines of the detected cylinders, green dashed lines are the missing connecting segments. . . . . . . . . . . . . . . . . . . . . . 4.22 Collinearity-based merging of segments. Segments to be merged are in blue, their starting points are in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Trees (in green) with corrected discontinuities (data 1). Initial trees are in black. Chordae starting points are in red, branching points are in magenta, terminal points are in yellow. Chordae mesh is in pink. Corrected discontinuities are circled in blue, not corrected case is circled in red. . . . . . . . . . . . . . . . . . 5.1 Example of the initial (a) and the final (b) states of the simulation in the case of presence of slack chordae. Chordae are in black, the slack ones are in red, the leaflets are shown in blue. Arrows depict the pressure applied to the leaflets, terminal and zero-displacement annulus points are in yellow. . . . . . . . . . . . 5.2 Compatibility of the chordae with the biomechanical model. Chordae mesh is in pink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Example of a chord with false topology and badly positioned branching node (data 1): the initial branching point is in magenta, terminal points are in yellow. Chordae mesh is in pink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Angle constraints illustration on a chord shown in two orthogonal views: α is the angle between branching segments (in red), θ is the angle between the base segment (in green) and the branching segments bisector (in dashed gray), the terminal points are in yellow, the branching node is in magenta. The chordae mesh is in pink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Examples of the results obtained with straightness term minimization under the image-driven, angle and branching nodes coordinates constraints (data 3): the initial branching nodes are in magenta, the terminal points are in yellow. The chordae mesh is in pink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Cost map calculation. The chordae limits are in red. . . . . . . . . . . . . . . . . 5.7 Cost map calculation from the skeleton. Gray zone corresponds to the plane P of the 2D slice. Example from data 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Examples of cost maps obtained with different power function factors on synthetic case with a hole in skeleton (shown in red). . . . . . . . . . . . . . . . . . . . . . 5.9 Binary volumes used for the cost map calculation. The chordae limits are in black, the skeleton points are in yellow. Example on data 2. . . . . . . . . . . . . . . . . 5.10 Functions applied to the distance map values for the cost map calculation. The red segment corresponds to the inside of the chordae, the blue one to the buffer zone and the green one to the outer volume. The values on the y-axis are given on the logarithmic scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 5.11 Example of the discretization of chordae segments (red points) and the corresponding cost map value for a point estimated by rounding (in green in the caption). . 5.12 Example of candidate points set (in blue) for branching nodes of the chord (in magenta). Terminal points are in yellow. . . . . . . . . . . . . . . . . . . . . . . . 5.13 Example of a chord (data 1) with faulty topology in two orthogonal views. The terminal points are in yellow, the initial branching node is in cyan, the branching nodes after the correction are in blue. Circled zones correspond to the overlapping bifurcations. The chordae mesh is in pink. . . . . . . . . . . . . . . . . . . . . . 5.14 Schematic examples of the correction of the false topology configurations (leftmost). Initial branching node is shown in yellow, added branching nodes are in green, ending points are in gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Scheme of possible branching in the chordae: cyan is the initial bifurcation node, gray corresponds to the bifurcations to be added, red corresponds to the ending points. Blue dashed lines correspond to possible branching type, black lines unfold the two nodes in the possible branching type. Rectangular zones limit the levels of branching. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Example of the solution scheme traversal: cyan is the initial bifurcation node, blue are the nodes added to the tree, gray are the bifurcation nodes, red are the ending points. Blue arrows correspond to the split direction, black lines unfold the two nodes in the possible branching type. Crossed nodes correspond to the tree leaves. Rectangular zones limit the levels of branching. . . . . . . . . . . . . . . . . . . 5.17 Schematic example of ending nodes (in gray) permutation for a given topology configuration (left-most). The initial bifurcation node is in yellow and the added bifurcation node is in green. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Candidate point selection for the topology optimization. . . . . . . . . . . . . . . 5.19 Example of optimization results for three possible topological configurations of a chord (data 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.20 Examples of correct topology optimization results for two chords (a, data 4) and (b, data 2) shown in two orthogonal views (left and right). The initial branching nodes are in cyan, the bifurcation nodes after the optimization are in blue and the terminal points are in yellow. The chordae mesh is in pink. . . . . . . . . . . . . 5.21 Examples of ambiguous topology cases for two chords (a, data 4) and (b, data 4) shown in two orthogonal views (left and right). The initial branching nodes are in cyan, the bifurcation nodes after the optimization are in blue and the terminal points are in yellow. The chordae mesh is in pink. . . . . . . . . . . . . . . . . . 5.22 Candidate points selection for branching nodes position optimization. . . . . . . 5.23 Example of nodes position optimization. The chordae mesh is in pink. . . . . . . 5.24 Decimation of the candidate points with a cubic grid. . . . . . . . . . . . . . . . . 5.25 Illustration of the size estimation of the whole chord (blue line) and of one of the branching zones (red line). Branching points are in yellow, terminal points are in gray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.26 Results of the chordal trees adjustment. Data 1. . . . . . . . . . . . . . . . . . . 5.27 Results of the chordal trees adjustment. Data 2. . . . . . . . . . . . . . . . . . . xi 5.28 Results of the chordal trees adjustment. Data 3. . . . . . . . . . . . . . . . . . . 5.29 Results of the chordal trees adjustment. Data 4. . . . . . . . . . . . . . . . . . . 5.30 Results of the chordal trees adjustment. Data 5. . . . . . . . . . . . . . . . . . . 5.31 Examples of the chordal trees with no change resulting from optimization procedure. 5.32 Example of the nodes position optimization results with weighted angle sum for two chordae (data 2) shown in two orthogonal views (left and right). The terminal nodes are in yellow. The chordae mesh is in pink. . . . . . . . . . . . . . . . . . 5.33 Results of the chordae optimization with R max variation. . . . . . . . . . . . . . 5.34 Boxplots of normalized distance between the bifurcation nodes coordinates with R max variation per data example. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Ground truth extraction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Examples of results of graph matching with the similarity maximization graph matching method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Choice of nodes proximity threshold. Initial chord is in blue, chord after optimization is in orange. S is the bounding box diagonal. The terminal tree points are in yellow, chordae mesh is in pink . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Calculation of the distance between two edges. Edge middle points are in gray. .6.5 GED max depending on the graph complexity. Node costs are in black, edge costs are in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Examples of the ambiguous for the segmentation chordae. . . . . . . . . . . . . . 6.7 Comparison of the cost function values obtained for the GT and segmented chordae after optimization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8 Distribution of normalized GED values. . . . . . . . . . . . . . . . . . . . . . . . 6.9 Examples of exact match of automatically segmented chord (in blue) with corresponding GT chord (in red). The matches between vertices are shown in green. . 6.10 Examples of correctly matched chordae with some inaccuracies: (a) Slight inaccuracy (GED norm = 0.17), data 3. (b) More pronounced inaccuracy (GED norm = 0.33), data 1. (c) Important parts missing (GED norm = 0.46), data 5. (d) Inaccurate topology (GED norm = 0.7), data 4. (e) Inaccurate topology (GED norm = 0.57), data 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Leaflet segmentation in ITK snap software. . . . . . . . . . . . . . . . . . . . . . 7.2 Extraction of the leaflet surface. The coaptation line is in yellow. . . . . . . . . . 7.3 Excised porcine mitral valve. Annulus is in yellow. The leaflet free margin is in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Mitral valve geometry in the biomechanical model (in two orthogonal views). The zero-displacement points (annulus and starting points of the chordae) are in yellow, terminal chordae points are in red. Leaflets are in pink. . . . . . . . . . . . . . . 7.5 Types of chordae connection to the leaflets. . . . . . . . . . . . . . . . . . . . . . 7.6 Unfold porcine mitral valve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.7 Definition of the mean direction (cyan) of a chord. . . . . . . . . . . . . . . . . . 7.8 Fan-shaped connection for terminal points being close to the leaflet. . . . . . . . 7.9 Fan-shaped connection for terminal points being far to the leaflet. . . . . . . . . . xii 7.10 One-point connection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.11 Connection of the chordae to the leaflets.

  (a)). When the ventricle contracts during systole, blood flow is ejected from the left ventricle to the aorta. The mitral valve closes and prevents the reverse flow into the left atrium (see figure 1(b)).

Figure 1 :

 1 Figure 1: Cardiac cycle.

Figure 2 :

 2 Figure 2: Mitral valve anatomy. Open (left) and closed (right) states (taken from [CAF10]).

Figure 3 :

 3 Figure 3: Chordae classification (taken from [Mur09]). Anterior leaflet is denoted as AoL, posterior leaflet is denoted as ML.

Figure 4 :

 4 Figure 4: Mitral valve pathologies (taken from [Tur04]).

  (a)), which allows to restore the valve function. As an alternative to resection, artificial chordae are used. They are introduced trough the papillary muscles and are adjustable during the intervention as shown in figure 5(c), or are applied as chordal loops based on prior measurements (see figure 5(d)) [CWMG19]. Annuloplasty is performed to prevent further annular dilatation by placing an annuloplasty device around the valve (see figure 5(b)).

  (a) Resection. (b) Annuloplasty. (c) Added adjustable chordae. (d) Added chordal loops.

Figure 5 :

 5 Figure 5: Examples of mitral valve repair procedures (taken from [CWMG19]).

  (a) 3D TEE (taken from [PHS09]). (b) MRI (taken from [CMO09]). (c) micro-CT (taken from [VHP + 18]).

Figure 1

 1 Figure 1.1: MV imaging using different modalities.

Figure 1 . 2 :

 12 Figure 1.2: Valve model based on geometry from post-mortem surgical studies proposed in [KRC97], [CK98].

1. 1 . 2

 12 Geometry based on measurements after dissection for the chordae and image-based for the leaflets With the development of image acquisition techniques, studies focused on subject-specific models, which implied extraction of the valve geometry directly from images. Modalities that are most frequently used in clinics are echography and MRI. The chordal structure is not distinguishable with such images and hence is reconstructed based on the post-mortem measurements of dissected valves. Leaflet extraction may be performed using various approaches. One of the methods is to reconstruct the leaflet geometry from the segmented annulus. The approach proposed by Prot et al. [PHS09] consisted in the annulus plane extracting from 3D TEE images. Then the leaflets surface was reconstructed from the ex vivo measurements (see figure 1.3(a)). The same technique was used in subsequent studies of Skallerud et al. [SPN11]. Studies of Votta et al. [VAI + 09] and Stevanella et al. [SMAC + 11] exploited a similar method to represent the leaflets.

( a )

 a Valve model. (b) Ex vivo measurements.

Figure 1 . 3 :

 13 Figure 1.3: Valve model based on both image data and post-mortem measurements proposed in [PHS09].

  .4(a)). In the work of Gaidulis et al. [GSZ + 19] the leaflets and the annulus were delineated manually per slice (see figure 1.4(b)). Similarly, in the models proposed by Pappalardo et al. [PSO + 17] and Gao et al. [GQF + 14] the leaflets were extracted manually from MRI (see figures 1.4(c) and 1.4(d)). Schneider et al. proposed a method for automatic leaflet extraction based on graph cuts [SPV + 10] (see figure 1.5(a)). Stoll et al. [ELAM + 15], [SSW + 16] developed an automatic pipeline for patient-specific modeling using 4D TEE data. Leaflets were segmented automatically with a thin tissue detector and k-means clustering (suppressing of the components thicker than the leaflet supposed to be) allowing as well the separation between exterior and interior surfaces (see figure 1.5(b)). Promising results of automatic leaflet segmentation from the echography were obtained by Costa et al. [CMS + 19]. The method is based on convolutional neural networks with U-net architecture allowing the identification of the leaflets. As it was noted by the authors, due to the significant number of false positives the segmentation requires a post-processing step.

  from 3D TEE proposed in [GSZ + 19]. (c) Manual leaflet segmentation from MRI proposed in [PSO + 17]. (d) Manual leaflet segmentation from MRI proposed in [GQF + 14].

Figure 1 . 4 :

 14 Figure 1.4: Manually segmented leaflet geometries.

  (a) Automatic annulus extraction proposed in [SPV + 10] (detected annulus is in blue, ground truth is in red). (b) Leaflet segmentation proposed in [SSW + 16].

Figure 1 . 5 :

 15 Figure 1.5: Automatic leaflet segmentation approaches.

Figure 1 . 6 :

 16 Figure 1.6: Valve model model proposed in [WS13].

Figure 1 . 7 :

 17 Figure 1.7: Leaflet segmentation process proposed in [TNH + 18].

  (a) Heart simulator. (b) Resulting valve segmentation.(c) Fiducial markers placement.

Figure 1 . 8 :

 18 Figure 1.8: Valve geometry acquisition with an in vitro setup proposed in [Kha18].

  dynamics [LDSB10], [GQF + 17].

Figure 1 . 10 :

 110 Figure 1.10: Comparison of two FSI-based models (tubular-shape and ventricular-shape fluid volumes) with a structural model by leaflets coaptation zone profile in 2D in [LDSB10].

Figure 1 .

 1 Figure 1.11: Valve stress-strain behavior in fiber and cross-fibre directions proposed in [KEC07].

Figure 2 . 1 :

 21 Figure 2.1: Ex vivo valve open state.

  (a) Valve photo. Annulus is shown in yellow. (b) Corresponding micro-CT image. Leaflets area is circled in red.

Figure 2 . 2 :

 22 Figure 2.2: Ex vivo valve closed state.

Figure 2 . 3 :

 23 Figure 2.3: Data acquisition : micro-CT device with pressurized air tube going inside and its screen (in red and zoomed on the right) displaying an X-ray of a porcine heart.

  (a) Micro-motion artifacts (circle red zones). (b) Artifacts due to the extrinsic objects: circles (left), air tube (right).

Figure 2 . 4 :

 24 Figure 2.4: Examples of the artifacts present in the data.

Figure 2 . 8 :

 28 Figure 2.8: Proportion of data exploitable data in the dataset.

  (a) Close-up in the micro-CT image. (b)Marked-up valve components: ventricle wall is in blue, chord is in red.

Figure 2

 2 Figure 2.5: Example of chordae being close to the ventricle wall.

  (a) Close-up in the micro-CT image. (b) Marked-up components: separated chords are in blue and yellow.

Figure 2

 2 Figure 2.6: Example of entangling chordae in the data (data 2).

Figure 2

 2 Figure 2.7: Example of a blood clot in the data (data 1).

Figure 2 . 10 :

 210 Figure 2.10: Illustration of a mesh triangle deformation.

Figure 2 . 9 :

 29 Figure 2.9: Illustration of the boundary conditions. The leaflets are in pink, the chordae in black. Yellow points correspond to zero-displacement nodes.

Figure 2 .

 2 Figure 2.12: Profile of the contact force depending on the distance d.

  (a) Initial (open) state. (b) Final (closed) state.

Figure 2 . 13 :

 213 Figure 2.13: Example of the outcome of the valve closure simulation obtained with the biomechanical model. The leaflets are in pink, the zero-displacement points are in yellow.

  (a) Mean z-component of velocity. (b) Mean velocity norm.

Figure 2 . 14 :

 214 Figure 2.14: Mean node velocity convergence curves.

  .13(b). In such conditions, the definition of the coaptation line is difficult (shown in blue in figure 2.15(a)) because of the numerous folds of the leaflet surface (shown in red in figure 2.15(b)). These folds are formed either due to the tethering by chords in the point of connection or by folding of the leaflet tissue. Such a configuration considerably complicates the segmentation and requires a high level of expertise.

  (a) Top view. (b) Bottom view.

Figure 2 . 15 :

 215 Figure 2.15: Leaflet surface of the closed valve.

Figure 2 .

 2 Figure 2.16: Excised porcine mitral valve (taken from [Ham11]). Annulus is in yellow. The leaflet free margin is in blue.

  (a) 2D slice of the micro-CT scan in sagittal view. (b) Labeled 2D slice of the micro-CT scan in sagittal view.

Figure 3

 3 Figure 3.1: Soft tissues present in the micro-CT scans.

  (a) Example of Frangi filtering results. 2D slice in the axial view. (b) Example of RORPO method result.2D slice in the axial view. (c) Example of Frangi method result. 3D rendering. (d) Example of RORPO method result. 3D rendering.

Figure 3

 3 Figure 3.2: Test of the methods for tubular structures extraction (data 1). Extracted chordae are in red.

  (a) 2D slice in axial view of the initial micro-CT image (data 2). (b) 2D slice in axial view of the micro-CT image after the nonlocal filtering (data 2).

Figure 3

 3 Figure 3.3: Example of non-local filtering.

  (a) Contours extracted with Canny filter in a form of a point cloud (data 4). (b) Meshed corresponding to the image contours (data 4).

Figure 3

 3 Figure 3.4: Contours and mesh extraction.

  (a) Surfaces with different topological characteristics. (b) Segmentation on synthetic data. Top is marked as chordae. Bottom corresponds to non-chordae.

Figure 3

 3 Figure 3.5: Topological characteristics of non-chordae, cylinder and bifurcation surfaces.

  (a) Non-chordae structure. (b) Tubular structure.(c) Branching.

Figure 3 . 6 :

 36 Figure 3.6: Different types of topology. Locally extracted sub-mesh is shown in green. Central point of the sub-mesh is in red. The data mesh points are in blue.

  Geodesic path between two points on a surface. (c) Geodesic distance map: iso-geodesic distance lines are in red, geodesic paths examples are in yellow (taken from [Pey09]).

Figure 3

 3 Figure 3.7: Sub-mesh size definition based on geodesic distance.

  (a) Central point (green) and its one-ring (red). (b) Choice of the closest point (circled in red) to the central one from the set of neighboring points (in blue). (c) Approximated iso-geodesic line (red points) at the current propagation step, black arrows indicate the propagation paths directions.

Figure 3

 3 Figure 3.8: Mesh propagation scheme.

  (a) Initial step of the algorithm. (b) Intermediate step of the algorithm. (c) Final step of the algorithm.

Figure 3

 3 Figure3.9: Computation of the sub-mesh by mesh propagation (lines correspond to geodesic paths approximation, green points are the sub-mesh points).

Figure 3 . 10 :

 310 Figure 3.10: Comparison of geodesic paths calculated from the distance maps (lines with black markers) with the paths resulting from the algorithm (lines with asterisk markers).

( a )

 a Exhaustive approach (each point is labeled). (b) Faster approach (all points of a sub-mesh are labeled at once).

Figure 3 .

 3 Figure 3.11: Segmentation results (data 1). Non-chordae structures are in blue, chordae are shown in red, branching zones in green and magenta. False detections are in black rectangles. Differently labeled with two methods branching zones are in orange rectangles.

  (a) Uniformly distributed sample of central points for sub-mesh extraction (red points). Data points are in blue. (b) Labeled sub-meshes corresponding to the points of the uniformly distributed sample. Data points are in blue. (c) Re-labeling of the overlaying sub-meshes (central points are in blue, the overlaying zone is in rectangular).

Figure 3 .

 3 Figure 3.12: Labeling process.

Figure 3 . 13 :

 313 Figure 3.13: Results of the chordae segmentation on 5 micro-CT images. Blue corresponds to the whole boundary mesh, the segmented chordae are in red.

  (a) Axial view. (b) Coronal view. (c) Sagittal view. (d) Axial view. (e) Coronal view. (f) Sagittal view.

Figure 3 . 14 :

 314 Figure 3.14: Superposition of the segmentation results (red contour) with the ground truth (blue inside) shown in axial, coronal and sagittal views (data 1). Yellow lines show correspondence between the slices. First line: for a point at half maximum height of the chordae. Second line: for a point near the chordae ending close to the leaflets.

Figure 4 . 1 :

 41 Figure 4.1: Examples of incorrectly detected type-1 shapes (circled in black) with a topologybased segmentation method. Chordae are in red, non-chordae structures are in blue.

Figure 4

 4 Figure 4.2: Examples of the results obtained with the skeletonization methods proposed in [CTO + 10] (ROSA 2.0) and in [TAOZ12] (part of CGAL). Black zones correspond to correct detections. Red and magenta zones indicate regions where the outcomes are not satisfactory.

  Taken from[START_REF] Li | 3d multi-branch tubular surface and centerline extraction with 4d iterative key points[END_REF].

Figure 4 . 3 :

 43 Figure 4.3: Examples of branchings profiles typically considered in the vascular structures segmentation methods.

  Resulting tree of line segments after the omission of intermediate points.

Figure 4 . 4 :

 44 Figure 4.4: Tree transformation. Terminal points are in yellow, branching points in magenta, intermediate points in gray.

Figure 4

 4 Figure 4.5: Parametric cylinder model.

Figure 4 . 6 :

 46 Figure 4.6: Model fitting procedure.

  (a) Ellipse (in blue) fitted to the sample points projection (in magenta). (b) Data points projected on the fitted ellipse. Inliers are shown in green, outliers are in red.

Figure 4 . 7 :

 47 Figure 4.7: Ellipse fitting procedure.

Figure 4 . 8 :

 48 Figure 4.8: Example of a semi-tube structure shown in two views (circled in black). Chordae are in red, non-chordae structures are in blue.

Figure 4

 4 Figure 4.9: Distribution uniformity check. Data points projected on the cylinder base are in blue.

Figure 4 . 10 :

 410 Figure 4.10: Example of a fitted cylinder model (in light blue) to data points (in blue).

  (a) Initial step. (b) Series of fitted cylinders.(c) Resulting tree.

Figure 4 .

 4 Figure 4.11: Tree tracking procedure. Chordae mesh points are in gray, initial tracking point is in magenta, locally extracted sub-meshes are in pink, fitted cylinders are in blue, initial points for the sub-mesh extraction are in yellow, type-2 shape sub-mesh boundary points are in green, branching zone mesh is in red, branching zone sub-mesh boundary points are in orange. See video link for 3D visualization.

Figure 4 .

 4 Figure 4.12: Chordal tree tracking procedure.

  (a) Sub-mesh is in blue. Fitted cylinder is in dashed black.(b) Sub-mesh (in red) enclosed between fitted cylinder bases (in blue).

Figure 4 . 13 :

 413 Figure 4.13: Detection of a missed branching zone.

  Example of a correctly orientated submesh.

Figure 4 .

 4 Figure 4.14: Non-valid cases of a potential branching zone mesh (in red). Chordae mesh is in pink. Initial sub-mesh point in in yellow. Gravity centers of the mesh boundary connected components are in green.

Figure 4 . 15 :

 415 Figure 4.15: Tree tracking at a branching zone. Chordae mesh points are in gray, initial points for sub-mesh extraction are in yellow, locally extracted sub-meshes are in pink, fitted cylinders are in blue, branching zone sub-mesh is in red, branching zone sub-mesh boundary points are in red.

  (a) Data 1. N umiter = 2000. (b) Data 1. N umiter = 3000. (c) Data 1. N umiter = 5000.

Figure 4 .

 4 Figure 4.16: Variation of the maximum number of the iterations. Each chord is in separate color. The problematic zones are circled in red.

Figure 4 . 18 :

 418 Figure 4.18: Results of the tree extraction. Each chord is in separate color.

  (a) Tree extraction #1. (b) Tree extraction #2. (c) Tree extraction #3. (d) Tree extraction #4. (e) Tree extraction #5.

Figure 4 . 19 :

 419 Figure 4.19: Results of the tree extraction for data 3 obtained with 5 iteration of the method. Each chord is in separate color. The problematic zones are circled in red. 65

( a )

 a Chordae as connected line segments with the intermediate points shown in blue in the close-up. (b) Chordae as connected line segments.

Figure 4 . 20 :

 420 Figure 4.20: Representation of the chordae in a form of connected line segments (data 1). Chordae segments are shown in black, chordae mesh is in gray. Starting points of chordal trees are in red, branching points are in magenta and terminal points are in yellow.

Figure 4 .

 4 Figure 4.21: Example of discontinuities in the trees. Magenta and blue zones correspond to the tree elements preceding and following the discontinuities, respectively. Red segments are the connected centerlines of the detected cylinders, green dashed lines are the missing connecting segments.

  (a) Starting point is a bifurcation. (b) Starting point is not a bifurcation.

Figure 4 .

 4 Figure 4.22: Collinearity-based merging of segments. Segments to be merged are in blue, their starting points are in red.

Figure 4 .

 4 Figure 4.23: Trees (in green) with corrected discontinuities (data 1). Initial trees are in black. Chordae starting points are in red, branching points are in magenta, terminal points are in yellow. Chordae mesh is in pink. Corrected discontinuities are circled in blue, not corrected case is circled in red.

Figure 5 . 1 :

 51 Figure 5.1: Example of the initial (a) and the final (b) states of the simulation in the case of presence of slack chordae. Chordae are in black, the slack ones are in red, the leaflets are shown in blue. Arrows depict the pressure applied to the leaflets, terminal and zero-displacement annulus points are in yellow.

  (a) Diagram of an outcome of the geometry extraction method. Fitted cylinders are shown in blue, branching zone is in red. (b) Diagram of a solution more suitable for the mechanical model (in magenta) comparing to the outcome of the geometry extraction method (in black).

Figure 5 . 2 :

 52 Figure 5.2: Compatibility of the chordae with the biomechanical model. Chordae mesh is in pink.

  (a) Segments before the correction. (b) Two examples of possible corrected configurations (blue and purple).

Figure 5 . 3 :

 53 Figure 5.3: Example of a chord with false topology and badly positioned branching node (data 1): the initial branching point is in magenta, terminal points are in yellow. Chordae mesh is in pink.

  • ] respectively. These values come from the literature data obtained for ovine MV chordae [KDB+ 17] and are consistent with our data.

Figure 5 . 4 :

 54 Figure5.4: Angle constraints illustration on a chord shown in two orthogonal views: α is the angle between branching segments (in red), θ is the angle between the base segment (in green) and the branching segments bisector (in dashed gray), the terminal points are in yellow, the branching node is in magenta. The chordae mesh is in pink.

  Non-satisfactory optimization results obtained with respect to the constraints.

Figure 5 . 5 :

 55 Figure 5.5: Examples of the results obtained with straightness term minimization under the image-driven, angle and branching nodes coordinates constraints (data 3): the initial branching nodes are in magenta, the terminal points are in yellow. The chordae mesh is in pink.

  (a) 2D slice of distance map D with values indicating the distance to the chordae skeleton. (b) 2D slice of cost map C with values used for the cost function evaluation.

Figure 5 . 6 :

 56 Figure 5.6: Cost map calculation. The chordae limits are in red.

  (a) Skeletonization of a chordal tree. Meshes corresponding to the tubular parts are in blue, cylinders fitted to the sub-meshes are in magenta. Meshes corresponding to bifurcation zones shown in red. The points of the skeleton are in yellow. The chordae mesh is in pink. (b) 2D slice of the cost map in the plane P. The chordae limits are in black.

Figure 5 . 7 :

 57 Figure 5.7: Cost map calculation from the skeleton. Gray zone corresponds to the plane P of the 2D slice. Example from data 1.

  (a) Power function factor 2 for the inner volume. (b) Power function factor 3 for the inner volume.

Figure 5

 5 Figure 5.8: Examples of cost maps obtained with different power function factors on synthetic case with a hole in skeleton (shown in red).

  s j (b) t , t ∈ [1, l] (see in red in figure 5.11). The corresponding value of the cost map C(p s j (b) t

Figure 5

 5 Figure 5.9: Binary volumes used for the cost map calculation. The chordae limits are in black, the skeleton points are in yellow. Example on data 2.

Figure 5 . 10 :

 510 Figure 5.10: Functions applied to the distance map values for the cost map calculation. The red segment corresponds to the inside of the chordae, the blue one to the buffer zone and the green one to the outer volume. The values on the y-axis are given on the logarithmic scale.

  where |L i | is the size of the list of possible coordinates for the i-th branching node.

Figure 5 .

 5 Figure 5.11: Example of the discretization of chordae segments (red points) and the corresponding cost map value for a point estimated by rounding (in green in the caption).

Figure 5 .

 5 Figure 5.12: Example of candidate points set (in blue) for branching nodes of the chord (in magenta). Terminal points are in yellow.

  (a) Possible configurations for tree ending nodes. (b) Possible configurations for four ending nodes. (c) Possible configurations for five ending nodes.

Figure 5 . 14 :

 514 Figure 5.14: Schematic examples of the correction of the false topology configurations (left-most). Initial branching node is shown in yellow, added branching nodes are in green, ending points are in gray.

Figure 5 .

 5 Figure 5.16: Example of the solution scheme traversal: cyan is the initial bifurcation node, blue are the nodes added to the tree, gray are the bifurcation nodes, red are the ending points.Blue arrows correspond to the split direction, black lines unfold the two nodes in the possible branching type. Crossed nodes correspond to the tree leaves. Rectangular zones limit the levels of branching.

Figure 5 . 17 :

 517 Figure 5.17: Schematic example of ending nodes (in gray) permutation for a given topology configuration (left-most). The initial bifurcation node is in yellow and the added bifurcation node is in green.

  (a) Sub-mesh definition: the bifurcation point to be multiplied is in cyan, points for the limiting planes definition are in blue, terminal points are in yellow.(b) Decimation of the candidate points set: before is in blue, after is in green.

Figure 5 . 18 :

 518 Figure 5.18: Candidate point selection for the topology optimization.

Figure 5 . 19 :

 519 Figure 5.19: Example of optimization results for three possible topological configurations of a chord (data 1).

Figure 5 . 20 :

 520 Figure 5.20: Examples of correct topology optimization results for two chords (a, data 4) and (b, data 2) shown in two orthogonal views (left and right). The initial branching nodes are in cyan, the bifurcation nodes after the optimization are in blue and the terminal points are in yellow. The chordae mesh is in pink.

Figure 5

 5 Figure 5.21: Examples of ambiguous topology cases for two chords (a, data 4) and (b, data 4) shown in two orthogonal views (left and right). The initial branching nodes are in cyan, the bifurcation nodes after the optimization are in blue and the terminal points are in yellow. The chordae mesh is in pink.

  (a) Sub-mesh definition: the bifurcation point is in cyan, points for the limiting planes definition are in blue, the terminal points are in yellow.(b) Decimation of the candidate points set: all the points of the sub-mesh are in blue, points limited by 1 2 Rmax are in gray, points after the downsampling are in green.

Figure 5 .

 5 Figure 5.22: Candidate points selection for branching nodes position optimization.

  (a) Result of the first optimization step. (b) Choice of the candidate points for the second optimization step. Sub-mesh points limited by 1 2 Rmax are in gray, candidate points after downsampling are in magenta. (c) Results of both optimization steps.

Figure 5 .

 5 Figure 5.23: Example of nodes position optimization. The chordae mesh is in pink.

Figure 5 .

 5 Figure 5.24: Decimation of the candidate points with a cubic grid.

Figure 5 .

 5 Figure 5.25: Illustration of the size estimation of the whole chord (blue line) and of one of the branching zones (red line). Branching points are in yellow, terminal points are in gray.

Figure 5

 5 Figure 5.26: Results of the chordal trees adjustment. Data 1. 91

Figure 5 Figure 5

 55 Figure 5.27: Results of the chordal trees adjustment. Data 2.

(

Figure 5 .

 5 Figure 5.29: Results of the chordal trees adjustment. Data 4.

Figure 5

 5 Figure 5.30: Results of the chordal trees adjustment. Data 5.

Figure 5 .

 5 Figure 5.31: Examples of the chordal trees with no change resulting from optimization procedure.

Figure 5

 5 Figure 5.32: Example of the nodes position optimization results with weighted angle sum for two chordae (data 2) shown in two orthogonal views (left and right). The terminal nodes are in yellow. The chordae mesh is in pink.

Figure 5 .

 5 Figure 5.33: Results of the chordae optimization with R max variation.

( a )

 a Maximal chordae size = 0.8Rmax. (b) Maximal chordae size = 1.2Rmax.

Figure 5 .

 5 Figure 5.34: Boxplots of normalized distance between the bifurcation nodes coordinates with R max variation per data example.

( a )

 a Chordae segmented in 3D Slicer (a separate color per chord), data 1. (b) Example of the extracted chordal tree (in black) from the segmented surface (in pink). The centerline points are in blue, the terminal points are in yellow, the bifurcation point is in magenta.

Figure 6 . 1 :

 61 Figure 6.1: Ground truth extraction.

  (a) False matching. (b) Correct matching.

Figure 6 . 2 :

 62 Figure 6.2: Examples of results of graph matching with the similarity maximization graph matching method.

Figure 6 . 3 :

 63 Figure 6.3: Choice of nodes proximity threshold. Initial chord is in blue, chord after optimization is in orange. S is the bounding box diagonal. The terminal tree points are in yellow, chordae mesh is in pink Figure 6.4: Calculation of the distance between two edges. Edge middle points are in gray.

Figure 6 . 5 :

 65 Figure 6.5: GED max depending on the graph complexity. Node costs are in black, edge costs are in red.

  (a) Entangling chordae (separate color per chord). (b) Blood clot (in red) on a chord (in blue).

  (c) Chord (in red) adhesion to the ventricle wall (blue).

Figure 6

 6 Figure 6.6: Examples of the ambiguous for the segmentation chordae.

( a )

 a Difference in values of Cim. (b) Difference in values of C straight .

Figure 6 . 7 :

 67 Figure 6.7: Comparison of the cost function values obtained for the GT and segmented chordae after optimization.

Figure 6

 6 Figure 6.8: Distribution of normalized GED values.

Figure 6 . 9 :

 69 Figure 6.9: Examples of exact match of automatically segmented chord (in blue) with corresponding GT chord (in red). The matches between vertices are shown in green.

Figure 6 . 10 :

 610 Figure 6.10: Examples of correctly matched chordae with some inaccuracies: (a) Slight inaccuracy (GED norm = 0.17), data 3. (b) More pronounced inaccuracy (GED norm = 0.33), data 1. (c) Important parts missing (GED norm = 0.46), data 5. (d) Inaccurate topology (GED norm = 0.7), data 4. (e) Inaccurate topology (GED norm = 0.57), data 1.

  (a) Axial 2D slice of the image (data 1). (b) Segmented axial 2D slice of the image (data 1).

Figure 7 . 1 :

 71 Figure 7.1: Leaflet segmentation in ITK snap software.

  (a) 3D mesh of the segmented leaflets. (b) Separated upper surface of the leaflet after decimation.

Figure 7

 7 Figure 7.2: Extraction of the leaflet surface. The coaptation line is in yellow.

  (a) Example of entangled chordae close to the leaflet surface. (b) Example of a cleft chord near the commissure.

Figure 7

 7 Figure 7.3: Excised porcine mitral valve. Annulus is in yellow. The leaflet free margin is in blue.

Figure 7 . 4 :

 74 Figure 7.4: Mitral valve geometry in the biomechanical model (in two orthogonal views). The zero-displacement points (annulus and starting points of the chordae) are in yellow, terminal chordae points are in red. Leaflets are in pink.

  (a) Representation of chordae with fan-shaped structures starting at papillary muscles (from [PHS09]). (b) Connection of the chordae extremities to the leaflet surface with a fan-shaped structure (from [DKS17]).

Figure 7

 7 Figure 7.5: Types of chordae connection to the leaflets.

Figure 7

 7 Figure 7.6: Unfold porcine mitral valve.

Figure 7 . 7 :

 77 Figure 7.7: Definition of the mean direction (cyan) of a chord.

  (a) Selection of the candidate points for the connection (limited by a sphere displayed in cyan). Connections are displayed in dash black. (b) Connecting segments (displayed in red).

Figure 7

 7 Figure 7.8: Fan-shaped connection for terminal points being close to the leaflet.

  (a) Selection of the candidate points for the connection (limited by a cone displayed in cyan). Connections are displayed in dash black. (b) Connecting segments (displayed in red).

Figure 7

 7 Figure7.9: Fan-shaped connection for terminal points being far to the leaflet.

  (a) Connection to the intersection point. (b) Connection to the nearest leaflet point.

Figure 7 .

 7 Figure 7.10: One-point connection.

  (a) Fan-shaped connection (b) One-point connection.

Figure 7 . 13 :

 713 Figure 7.13: Choice of the connection type. Ambiguous cases. Connections are shown in red, the terminal segment point is in yellow. The leaflets are in pink.

  (a) Manual segmentation.(b) Our segmentation method with optimization.

Figure 7 . 14 :

 714 Figure 7.14: Chordae configurations used for the simulation. Connection points are in blue. Zero-displacement points are in yellow. The leaflets are in pink.

Figure 7 . 15 :

 715 Figure 7.15: Manually segmented chordae with the chordae resulting from our segmentation procedure (shown in two views). Initial chordae points are in red, bifurcation points in magenta, terminal points in blue. The chordae not detected by the automatic segmentation approach and their terminal points are in green. The annulus points are in yellow, the leaflets are in pink.

Figure 7 .

 7 Figure 7.16: Estimation of the initial leaflet strains.

  (a) Cut plane.(b) Profiles corresponding to the leaflets in pre-strained and not pre-strained conditions.

Figure 7 . 17 :

 717 Figure 7.17: Comparison of leaflet profiles with and without pre-strain shown in the cut plane.

Figure 7 . 18 :

 718 Figure 7.18: Bulging volume calculation.

Figure 7 . 19 :

 719 Figure 7.19: Computation of volume between two faces as a sum of volumes of four tetrahedrons (each in separate color). The lower and upper faces correspond to the initial and final leaflets states.

Figure 7 . 20 :

 720 Figure 7.20: Profiles corresponding to the simulation outcomes obtained with the three chordae sets.

  (a) Our method minus manual segmentation. (b) Induced slack case minus manual segmentation.

Figure 7 .

 7 Figure 7.21: Map of bulging volume differences.

Figure 7 .

 7 Figure 7.22: Von Mises stress distribution for the cases of the three chordae sets.

  (a) Manual segmentation -our segmentation method with optimization.(b) Manual segmentation -our segmentation method with induced slack.

Figure 7 .

 7 Figure 7.23: Difference of Von Mises stress distribution compared to the manually extracted chordae case.

  (a) Mean z-component of velocity. (b) Mean velocity norm.

Figure 7 .

 7 Figure 7.24: Node velocity convergence curves obtained for the optimized chordae set extracted with the proposed segmentation method.

  (a) Manual segmentation.(b) Our segmentation method with optimization.

Figure 7 .

 7 Figure 7.25: Initial configuration of the pathological case simulation for the two chordae sets. Deleted chordae are in red.

( a )

 a Manual segmentation. (b) Our segmentation method with optimization.

Figure 7 .

 7 Figure 7.26: Final configuration of the pathological case simulation for the two chordal sets. Zones where the leaflets are not sealed are circled in red.

  (a) Manual segmentation (lateral view). (b) Our segmentation method with optimization (lateral view). (c) Manual segmentation (top view). (d) Our segmentation method with optimization (top view).

Figure 7 .

 7 Figure 7.27: Von Mises stress distribution obtained with the incomplete chordae sets shown in two views (lateral and top). Zones where the leaflets are not sealed are encircled in red.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Profile of the contact force depending on the distance d. . . . . . . . . . . . . . . 2.13 Example of the outcome of the valve closure simulation obtained with the biomechanical model. The leaflets are in pink, the zero-displacement points are in yellow. 2.14 Mean node velocity convergence curves. . . . . . . . . . . . . . . . . . . . . . . .

2.

15

Leaflet surface of the closed valve. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Excised porcine mitral valve (taken from [Ham11]). Annulus is in yellow. The leaflet free margin is in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Soft tissues present in the micro-CT scans. . . . . . . . . . . . . . . . . . . . . . 3.2 Test of the methods for tubular structures extraction (data 1). Extracted chordae are in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Example of non-local filtering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Connection points are in blue. Zerodisplacement points are in yellow. The leaflets are in pink. . . . . . . . . . . . . . 7.12 Node velocity convergence curves for two types of chordae connection. . . . . . . 7.13 Choice of the connection type. Ambiguous cases. Connections are shown in red, the terminal segment point is in yellow. The leaflets are in pink. . . . . . . . . . 7.14 Chordae configurations used for the simulation. Connection points are in blue. Zero-displacement points are in yellow. The leaflets are in pink. . . . . . . . . . . 7.15 Manually segmented chordae with the chordae resulting from our segmentation procedure (shown in two views). Initial chordae points are in red, bifurcation points in magenta, terminal points in blue. The chordae not detected by the automatic segmentation approach and their terminal points are in green. The annulus points are in yellow, the leaflets are in pink. . . . . . . . . . . . . . . . . 7.16 Estimation of the initial leaflet strains. . . . . . . . . . . . . . . . . . . . . . . . . 7.17 Comparison of leaflet profiles with and without pre-strain shown in the cut plane.7.18 Bulging volume calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.19 Computation of volume between two faces as a sum of volumes of four tetrahedrons .21 Map of bulging volume differences. . . . . . . . . . . . . . . . . . . . . . . . . . . 7.22 Von Mises stress distribution for the cases of the three chordae sets. . . . . . . . Results of the tree extraction for data 1 obtained with 5 iteration of the method. Each chord is in separate color. . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Results of the tree extraction for data 2 obtained with 5 iteration of the method. Each chord is in separate color. . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Results of the tree extraction for data 4 obtained with 5 iteration of the method. Each chord is in separate color. . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Results of the tree extraction for data 5 obtained with 5 iteration of the method. Model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Canny filtering parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	List of Figures
	List of Tables
	2.1
	xiv	xiii

(each in separate color). The lower and upper faces correspond to the initial and final leaflets states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.20 Profiles corresponding to the simulation outcomes obtained with the three chordae sets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77.23 Difference of Von Mises stress distribution compared to the manually extracted chordae case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.24 Node velocity convergence curves obtained for the optimized chordae set extracted with the proposed segmentation method. . . . . . . . . . . . . . . . . . . . . . . . 7.25 Initial configuration of the pathological case simulation for the two chordae sets. Deleted chordae are in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.26 Final configuration of the pathological case simulation for the two chordal sets. Zones where the leaflets are not sealed are circled in red. . . . . . . . . . . . . . . 7.27 Von Mises stress distribution obtained with the incomplete chordae sets shown in two views (lateral and top). Zones where the leaflets are not sealed are encircled in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Each chord is in separate color. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Geodesic distance values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 R max values used for the segmentation. . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Ground truth chordae discarded for the validation. . . . . . . . . . . . . . . . . . xv

Toma et al. [TBP + 17, TEBI + 17]. It permits to obtain high resolution data with a low level of noise. The segmentation is obtained manually with the use of commercial software (ScanIp, Zbrush).

	al. ([LOR + 13]), Drach, Khalighi et al.
	[DKS17], [Kha18],

  To be closer to clinical conditions, we do not use an in vitro setup where the excised valve movement is reproduced like in [LOR + 13],[START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF], [KDG + 18], therefore the ventricle walls are present in the images. This makes the leaflets indistinguishable from the ventricle walls in the open state (see figure2.1(b)), hence, images corresponding only to the closed state during systole peak are used (see figure2.2). In clinical conditions, the valve competence after the repair is checked by injecting of liquid under quasi-static pressure. This justifies the usage of only closed state for the

study. To simulate closed valve state, 100 mmHg air pressure was delivered through the aorta, which corresponds to normal peak systolic pressure. The two states of the valve are displayed in figures 2.1(a) and 2.2(a).

Table 2

 2

.1: Model parameters.

Table 3 .

 3 2: Geodesic distance values Geodesic, px Approximation, px Relative difference, %

	13.8427	16.0128	15.7
	13.5752	14.7587	8.71
	13.0107	13.3948	2.95
	13.8427	14.9095	7.7
	13.575	16.0479	18.2

Table 6 .

 6 1: Ground truth chordae discarded for the validation.

	Data ID Total segmented Total discarded	Issues
	1	20	2	blood clots, entangling chordae
	2	13	6	blood clots,entangling chordae, adhesion to ventricle
	3	32	13	ambiguous topology, image artifacts, entangling chordae
	4	30	5	ambiguous topology, image artifacts
	5	10	5	blood clots, entangling chordae, ambiguous topology

https://team.inria.fr/curative/research/

(a) Tree extraction #1. (b) Tree extraction #2. (c) Tree extraction #3. (d) Tree extraction #4. (e) Tree extraction #5.Figure B.4: Results of the tree extraction for data 5 obtained with 5 iteration of the method. Each chord is in separate color.

Chapter 6

Validation of the segmented chordae set based on graph similarity

The chordal trees extracted with the method presented in chapter 4 were modified with the optimization-based approach proposed in chapter 5. In this chapter we present a validation procedure aimed at quantitatively evaluating the quality of the obtained chordae segmentation.

The method exploits a graph-like representation of the automatically segmented chordal trees which is compared to the ground truth structures. We propose a metric based on graph similarity allowing to define the correspondence between the ground truth and automatically segmented chordae while measuring the accuracy of the segmentation.

Validation ground truth

The classical approach in image segmentation validation is the comparison of the obtained results with a ground truth (GT). In our case, the results of the chordae extraction are not the labeled images but the geometrical representation of the chordae which is necessary for the biomechanical simulation. More precisely, the obtained chordal trees can be seen as directed graphs. Since the final goal is the application in the biomechanical simulation, it is more reasonable to obtain directly the GT in form of graphs via manual segmentation.

The manual segmentation of mitral valve chordae is a narrow field of expertise, therefore it is hard to find an expert providing satisfactory results. After unsuccessful tests with the results provided by a clinician, we decided to perform the GT extraction ourselves by two informed users. The first informed user provided the segmentation of 5 CT scans. The segmentation of 1 image was obtained by the second user (less experienced than the first one). It was noticed that the computed solution was better at several places of the image than the manual delineation obtained by the second user, as it gives a more complete result by the number of detected branching zones and segments. Given the above, we could not compare the results using this ground truth and exploited only the manual segmentation provided by the first informed user.

The segmentation process included several stages. First, the chordae were manually delineated using the software 3D Slicer [FBKC + 12]. This tedious segmentation process resulted in a set of voxels transformed into meshes (see figure 6.1(a)), each mesh corresponding to an isolated chord. Afterwards, chord by chord, the centerlines of meshes (see blue points in figure 6.1(b)) Appendix A

Data parameters

Tree extraction results

The influence of the MLESAC procedure on the results of tree extraction is evaluated with 5 test iterations of the method for each data. The procedure proposed for the automatic extraction of the chordae consists of three steps. First, the segmentation of the chordae is obtained with the method based on the topological properties of the structures. Then, the representation of the geometry in a form of trees of connected line segments is extracted using the tracking approach combined with the local fitting of a cylindrical model. In order to make the geometry resulting from these methods consistent with the biomechanical model, the final stage of the chordae structure optimization is applied. In particular, the optimization makes it possible to reduce slack in the chordae and to correct false topological configurations by taking into account mechanical and anatomical considerations associated with image data.

We have proposed a validation metric for the segmentation results based on the graph similarity. With this metric correspondences between the chordae obtained with our automatic method and those resulting from the manual segmentation are defined, as well as the quality of the segmentation is assessed.

We then demonstrated the application of the extracted chordae geometry in the biomechanical simulation of the closed valve state. The consistency of the reproduced behavior was verified by comparing the results obtained with the chordae from our automatic procedure and the ground truth chordae. The bulging volume and the mechanical stress distribution have been used as quantitative criteria. The coherence of the simulation was also analyzed on a modeled pathological case with ruptured chordae.

Keywords: mitral valve, image-based modeling, segmentation, biomechanical simulation.