Tobias Heimann
email: tobias.heimann@inria.fr

Hervé Delingette

Model-based Segmentation

This chapter starts with a brief introduction into model-based segmentation, explaining the basic concepts and different approaches. Subsequently, two segmentation approaches are presented in more detail: First, the method of deformable simplex meshes is described, explaining the special properties of the simplex mesh and the formulation of the internal forces. Common choices for image forces are presented, and how to evolve the mesh to adapt to certain structures. Second, the method of point-based statistical shape models is described. The model construction process is explained and the point correspondence problem is treated in detail. Several approaches of how gray level appearance can be modeled are presented, and search algorithms that use this knowledge to segment the modeled structures in new images are described.

Introduction

Automatic medical image segmentation is such a challenging task because it involves highly variable objects that have to be extracted from images of very low quality. Often, lack of contrast or artifacts lead to missing data in the image to be segmented, and the boundary of the object of interest cannot be determined from local information as edges or region homogeneity alone. Medical experts are still able to delineate the object because they know what it is supposed to look like: They have a model of the object in their mind, they have an a-priori knowledge about its shape and appearance. Model-based segmentation methods strive to translate this knowledge into smart algorithms that have a prior knowledge about the structures of interest. Those methods can be qualified as top down and usually consists of two stages. The former initializes the location and appearance of the model. It is based on user input (mouse clicks in an image), the registration of an atlas or machine learning approaches. The latter optimizes the shape and appearance of the model such that it closely matches the ones measured in the images.

As the shape of a model is optimized during the segmentation process, one important property of model-based segmentation lies in the choice of the deformation strategy. There are two basic ways of deforming a shape: by deforming the embedding space of a shape and by modifying their parameters or Degree of Freedom (DOF). For instance, applying a rigid transform (i.e., rotation and translation) lies in the first category while modifying the vertex position of a mesh lies in the second category. Of course, those two mechanisms are often combined but their formalism and approaches are fairly different. Deforming a shape by deforming the embedding space can be qualified as a registration approach whereas modifying the parameters of a shape can be qualified as a deformable model approach. The geometric representation of the models is specifically important in the latter case. In Fig. 12.1, we show a coarse taxonomy of possible geometric representations. Each representation having pros and cons, the choice of a representation is often problem dependent. In Sec. 12.2, we detail the deformable model formulation based on simplex meshes. Among others, common representations are level-sets, triangulated meshes and spline curves or surfaces.

If many segmentations of a specific object are available, there are further possibilities to make model more specific. In addition to knowledge about an average template shape, Statistical Shape Models (SSMs) have also knowledge about the principal modes of variation. The SSMs presented in Sec. 12.3 are based on the Point Distribution Model (PDM). To perform statistical analysis on these shapes, point correspondences have to be known: a challenging problem, for which just recently automatic methods have been proposed [?]. The visual appearance of the object can also be captured by statistics on training samples, using techniques from computer vision and machine learning. Finally, there exist specific algorithms to employ the statistical information to drive a model towards an instance of the object in an image.

To conclude the chapter, we compare the concepts of a deformable simplex mesh and the SSM and give hints which technique to use under which circumstances.

Deformable Simplex Meshes

In a deformable model approach, the shape of a model is optimized in order to match that of a structure of interest in an image. This technique has been pioneered in 1987 by Terzopoulos et al. [START_REF] Terzopoulos | Multiresolution computation of visible-surface representation[END_REF] with the introduction of active contours or snakes [START_REF] Kass | Snakes: active contour models[END_REF]. This has been later generalized to active surfaces [START_REF] Metaxas | Constrained deformable superquadrics and nonrigid motion tracking[END_REF] but one difficulty arises when dealing with three-dimensional (3D) surfaces: the continuous parameterization of surfaces. Indeed, even for surfaces having the same topology as the sphere, it is not straightforward to define a C 1 parametric surface.

In order to avoid this issue, authors have proposed to work with implicit parameterization (e.g., using level-sets function) or with C 0 continuous meshes (e.g., using triangulated or simplex meshes). In the latter case, shape smoothing or regularization cannot be defined in the continuous domain (since the representation is not even C 1 ) but must be defined in a discrete fashion.

A simplex mesh is an original mesh representation that has been introduced in 1992 [START_REF] Delingette | Shape representation and image segmentation using deformable surfaces[END_REF] as a model-based segmentation method for range data and volumetric images. A k-simplex meshes can represent a manifold surface of dimension k (k = 2 for 3D surfaces) and can be defined as a k + 1-cell. We provide below a recursive definition of a cell (Fig. 12.2). Definition 1. A 0-cell is defined as a vertex P and a 1-cell as an edge, i.e., an unordered pair of distinct vertices (P, M ). We then recursively define a p-cell (p ≥ 2) C as the union of c (p -1)-cells, c ∈ N, such that:

1. Each vertex belonging to C also belongs to p distinct (p -1)-cells. 2. A (p -2)-cell belongs to 2 and only 2 (p -1)-cells.

The intersection of 2 (p -1)-cells is empty or is a

(p -2)-cell. 1-Cell= Edge 3-Cell 2-Cell= Face 0-Cell= Vertex Fig. 12.2. Examples of p-cells. A collection of p-cells is given for 0 ≤ p ≤ 3.
An edge is a 0-simplex mesh, a face is a 1-simplex mesh, and a 3-cell a 2-simplex mesh.

4.

A p-cell is simply connected, i.e. that given 2 vertices of that cell, there exists at least one set of edges that connect those 2 vertices.

Two fundamental topological properties result from this definition. First, a k-simplex mesh has a fixed vertex connectivity: every vertex is adjacent to k + 1 distinct vertices. Second, a k-simplex mesh is a topological dual of a k-triangulation. In this duality relationship, a triangle is associated with a vertex of a 2-simplex mesh, a triangulation edge with an edge and a vertex with a face (a 2-cell).

In a deformable 2-simplex mesh for image segmentation, the vertex position P i evolves over time under the influence of internal f int i and external forces f ext i capturing the influence of prior knowledge on shape and appearance, respectively

P t+1 i = P t i + (1 -δ)(P t i -P t-1 i ) + αf int i + βf ext i (12.1)
A description of those forces are the topic of the next two sections.

Internal Forces on Simplex Meshes

Local Internal Forces

There are two types of prior knowledge that can be applied. A weak shape prior consists in assuming that the structure of interest is smooth, which is a reasonable assumption for many anatomical structures (liver, heart, bones...) but not all of them (e.g., infiltrating tumor, brain gray matter). With a strong shape prior, it is assumed that the structure of interest has a typical shape, i.e.that its shape varies to a "small extent" around a reference shape. Simplex meshes are well suited to enforce both types of priors. Indeed, one can define at each vertex of a 2-simplex mesh, discrete geometric quantities that locally control the shape of the mesh. At vertex P i those quantities are the three metric parameters (ǫ 1i , ǫ 2i , ǫ 3i ) and the simplex angle φ i (Fig. 12.3, left). Metric parameters control the relative spacing of vertices through a tangential internal force while the simplex angle controls the local surface smoothness through a normal internal force. By combining those two tangential and normal components of the internal force, it is possible to define weak and strong shape priors in a simple way.

For instance, Fig. 12.3 (right) shows the effect of a specific tangential force that leads to a concentration of vertices at parts of high curvatures. In Fig. 12.4, the tangential component impose metrics parameters equal to 1/3, while the normal component of the internal forces smooths the surface with normal and mean curvature continuity. In some restricted cases, discrete mean curvature on a simplex mesh has been proved [START_REF] Delingette | déformation et Reconnaissance d'objets tridimensionnels a l'aide de maillages simplexes[END_REF] to converge towards its continuous value.

The straightforward definition of forces imposing mean curvature continuity is a major advantage of simplex meshes. Indeed, with such regularizing constraint, no shrinking effect is observed unlike the mean curvature flow [START_REF] Chopp | Computing minimal surfaces via level set curvature flow[END_REF] widely used in the level-sets method [START_REF] Malladi | Shape modeling with front propagation: a level set approach[END_REF]. With those flows minimizing surface areas, it is often necessary to limit the shrinking effect by adding a "balloon force", which may lead to additional problems [START_REF] Cohen | Finite element methods for active contour models and balloons from 2-D to 3-D[END_REF].

To provide a strong shape prior, a specific internal force has been proposed that constrains the metrics parameters and simplex angle to be equal to some user-defined values corresponding to a reference shape. Under the influence of this shape-memory force, a simplex mesh is attracted towards a reference shape up to a translation, rotation and scale. This local representation of shape is also a unique feature of simplex meshes compared to other representations.

A B C F P Fig. 12.3. Metric parameters. Left: point P and its three neighbors on a 2simplex mesh define a circumscribed sphere. The simplex angle in P is related to the height F P of point P with respect to triangle A, B, C such that the angle stays constant when P lies on the circumsphere. The metric parameters (ǫ1, ǫ2, ǫ3) are simply the barycentric coordinates of F in the triangle A, B, C; Middle: cube simplex mesh with metric parameters equal to 1/3; Right: mesh with optimized metric parameters so that vertices are concentrated at parts of high curvature. Fig. 12.4. Regularization of a simplex mesh. Left: simplex mesh with zero mean curvature; Middle: the red part is smoothed with a mean curvature continuity to the four cylinders; Right: normal continuity between the surface and 3D contours can also be enforced.

Image Forces

Image forces drive a deformable simplex mesh towards the apparent boundary of a structure of interest and rely on some prior knowledge about the appearance of that structure. In all cases, the external force is formulated as displacement in the normal direction towards the closest boundary point C i of the structure:

f ext i = ((C i -P i )•n i )n i .
Again one can rely on a weak appearance prior that assumes generic boundary features (such as gradient maxima) or on a strong appearance prior that relies on machine learning trained on segmented dataset. In fact, there is a continuum of priors between those two extremes:

• Contour prior : The closest boundary points are searched among the maxima of the gradient norm in the gradient direction. Additional criteria may be used to select the proper boundary point such as intensity values or the dot product between the surface normal n i and the gradient vector at that point [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF]. • Contour and region prior: The principle is to characterize the inner (or outer) region by a a range of intensity. The boundary is then determined as the voxel located next to this region whose gradient norm is greater than a threshold [START_REF] Montagnat | Cylindrical echocardiographic images segmentation based on 3D deformable models[END_REF]. • Region histogram prior: From a collection of segmented structure, one can build the probability density function that a pixel belongs to the structure from the region histogram (using Parzen windowing for instance). The closest point is then determined as the one which maximize (resp. minimize) the probability for the inner (resp. outer) voxels to belong to the inner region [START_REF] Herbulot | Segmentation of vectorial image features using shape gradients and information Measures[END_REF]. • Intensity profile matching prior: At each mesh vertex, an intensity profile along the normal direction is extracted. The external force is computed by searching along the normal direction the displacement which maximize the correlation with a reference intensity profile. That reference profile may result from a Principal Component Analysis (PCA) as described in Sec. 12.3.4. Several criteria may be used such as the sum of the square differences, the linear correlation or the Mahalanobis distance. • Block matching prior : This approach is similar to the intensity profile approach but based on the matching of image blocks rather than intensity profiles. The main difference is that the search for the closest point is no longer constrained to be along the normal direction but may be in all directions. This is of interest when homologous points can extracted from images (such as points of high curvatures or anatomically relevant features). • Texture classification prior : The image force can be applied after performing some image classification based on some machine learning method. The output of the classification is a probability map that a voxel belongs to the structure of interest. The closest point along the normal direction whose probability is equal to 50 % may then be used to compute the image forces. Linear classifiers, support vector machines [START_REF] Pitiot | Expert knowledge guided segmentation system for brain MRI[END_REF] or even neural networks [START_REF] Pitiot | Texture based MRI segmentation with a two-stage hybrid neural classifier[END_REF] can be trained on segmented datasets to classify the image.

Globally Constrained Deformation

As mentioned previously, the first stage of model-based segmentation consists in positioning the deformable model in the image through various methods, manual or automatic. Due to the convergence to a local minimum of the functional, the second stage corresponding to the optimization of the model parameters may not give good segmentation results if the initialization is not close enough to the final shape. One common strategy to decrease the sensitivity of the optimal shape to the initialization position is to use a coarseto-fine strategy. In this approach, the complexity (in terms of DOF) of the model is slowly increased in order to start the optimization with a small number of local minima and then track the global optimum during the increase of complexity.

A first method to control the DOF consists in using a multi-resolution approach where a coarse mesh is deformed on a given image and replaced with a finer mesh when convergence is reached. This technique is fairly straightforward to implement but requires to create meshes of increasing resolution which is not trivial on unstructured meshes [START_REF] Ciampalini | Multiresolution decimation based on global error[END_REF]. On simplex meshes, the √ 3 subdivision scheme is well suited to globally or even locally refine the mesh [START_REF] Kobbelt | Sqrt(3) subdivision[END_REF]. One limitation of multi-resolution schemes is that the coarser mesh The meshes are driven by the matching of intensity profiles; Upper row: reference CT image scan (left) and head segmented with a simplex mesh (right). Lower row: an intensity profile is stored at each vertex as a reference and this mesh is used to register the head of the same patient on an MRI: mesh position before (left) and after rigid registration (right). An alternative is to use globally constrained deformations [START_REF] Montagnat | Globally constrained deformable models for 3D object reconstruction[END_REF], where global spatial transforms are used in the early stages of segmentation in order to constrain the DOF. Table 12.1 lists some transforms that may be applied with a growing number of DOF. Starting from a rigid transform, it is possible to use only six DOF which greatly constrain the optimization problem.

However, is it often difficult to smoothly increase the number of DOFs, except when using the PCA (Sec. 12.3.3). To this end, the concept of globally constrained deformation has been introduced [START_REF] Montagnat | Globally constrained deformable models for 3D object reconstruction[END_REF] that can tune the complexity of the deformation with a locality parameter λ. More precisely, the motion equation (12.1) is refined by adding a global force f global i = T (P t i ) -P t i computed from the application of a global transform T

P t+1 i = P t i + (1 -δ)(P t i -P t-1 i ) + λ αf int i + βf ext i + (1 -λ)f global i (12.2)
When λ = 0 the simplex mesh deforms according to the application of a global transformation (with few DOFs) whereas when λ = 1, the deformation is only driven by local internal and external forces. When 0 < λ < 1 the deformation is a trade-off between global and local behaviors corresponding to an intermediate number of DOFs. The transformation T is computed such that it minimizes the discrepancy between the vertex position P t i and its closest point on the image C t i . In many cases, it is sufficient to apply global transforms from rigid to affine and then to slowly increase the locality parameter from 0 to 1 or to a smaller value (Fig. 12.6). The concept of globally constrained deformation may also be successfully applied to enforce to some degree the axis-symmetry around a deformable axis (for instance to segment vessels) or to take into account the statistical shape variation around a mean shape.

3D+T Deformable Simplex Meshes

For the analysis of time series of medical images, e.g., gated cardiac imaging, it is important to not only to segment a structure but also to track its deformation over time. A common approach consists in using the result of the segmentation at time t as the initialization of the segmentation at time t + 1.

However this approach is prone to the propagation of errors if one image is especially challenging and tends to bias the result noticeably by underestimating the motion of the structure.

Rather than processing the images one after the other, 3D+T deformable models [START_REF] Montagnat | 4D Deformable models with temporal constraints: application to 4D cardiac image segmentation[END_REF] can be used to process all images of the time series at once (if they can fit in the memory of the computer). In this framework, a family of simplex meshes having the same topology (Fig. 12.7) are deformed on each image of the time series, the mesh deformation at time t being coupled with that of time t -1 and t + 1.

The temporal coupling forces act as additional internal force and can be formulated in a weak or a strong way. The weak prior consists in minimizing the kinetic energy of the system and thus tries to make vertex P t i close to both points P t-1 i and P t+1 i . With a strong motion prior, the shape of the trajectories of each vertex are constrained to be similar to some reference trajectory curves, described with geometric quantities that are invariant to translation, rotation and scale. Examples of segmentation of time-series of echocardiographic images are shown in Fig. 12.8.

Advanced Segmentation Strategies

The segmentation of anatomical structures from medical images is often a difficult problem due to their low contrast with surrounding structures, the inter-subject variability, the presence of pathologies, etc. Problem specific segmentation strategies must be defined in order to achieve a robust and accurate delineation.

Atlas-Based Initialization

The first stage of model-based segmentation consists in finding the location of the structure of interest in the image. This can be done manually by the user but the inter-user variability in providing this information may be very large and may impair the reproducibility of the method. An interesting alternative is to use non-rigid registration [START_REF] Pitiot | Expert knowledge guided segmentation system for brain MRI[END_REF] to find a geometric transformation between a reference image (or an atlas) and the target image. In the reference image, meshes of the structures of interest can then be mapped by this transformation in the space of the target image and finally serve as initial shape for modelbased segmentation. This approach has the advantage of being automatic (and therefore reproducible) and usually provides good model initializations. However, when the structure location varies greatly, it is required to provide a suitable initialization of the non-rigid registration algorithm, which can be achieved by performing some coarse exhaustive search [START_REF] Montagnat | Volumetric medical image segmentation using shape constrained deformable models[END_REF].

Multi-Structure Segmentation

It is often a good strategy to first segment the most contrasted structures and then proceed with more challenging ones (least contrasted). Indeed, already delineated structures can restrict the possible shape of neighboring structures during the model-based segmentation of other structures. Furthermore, when boundaries between structures are hardly visible, it is best to segment both structures at the same time, the shape of the former providing interpenetration constraints [START_REF] Costa | Automatic segmentation of bladder and prostate using coupled 3D deformable models[END_REF] for the latter. Those constraints are based on the computation of distance maps measuring the distance of a given points to the two surface meshes. Fig. 12.9 shows an example of the joint segmentation of bladder and prostate.

Rule-based segmentation

Ideally, all model-based segmentation algorithms should proceed in two stages: initialization followed by an optimization stage. However, most segmentation algorithms are meant to be applied on patients with some pathologies (e.g., tumor, metastasis) or some rare anatomical configurations (e.g., abnormal branching of vessels, already ablated tissue). In such case, it is difficult to build a shape and appearance model that can cope with such variation of cases. Therefore, to improve the robustness of the segmentation, some alert mechanisms should be included in order to detect that the segmentation is not performing well. For instance, for segmentation of hippocampi in brain MRI, a projection on the main PCA mode was used [START_REF] Pitiot | Expert-Knowledge-Guided Segmentation System for Brain MRI[END_REF] to detect the leakage of a deformable model thus signaling that the locality parameter λ should be decreased. Additional rules and meta-rules can be used to test the validity of hypothesis and eventually reset the optimization stage with new parameters.

Geometric Representations for Model-Based Segmentation

The choice of a geometric representation of a model has a significant impact on the performance of a model-based segmentation algorithm. Many components of those models can be implemented independently of the representation, however some implementations may be far easier than others. For instance, PCA of the shape variability may be implemented for both point-set distributions [START_REF] Cootes | Active shape models -their training and application[END_REF] (Sec. 12.3.1) and level-set functions [START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF] but are far simpler for the former.

Simplex Meshes vs. Triangulations

Triangulations are widely used representations for and computational purposes. They can represent manifold and non-manifold surfaces while simplex meshes are only suitable for manifold surfaces. Furthermore, it is straightforward to compute the closest distance between a point to a triangulation because of its planar faces (triangles). On simplex meshes, the faces are curved and therefore computing accurately the closest point is more complex (no closed form solution). However, due to its fixed vertex connectivity, it is much easier to control the spreading of vertices on simplex meshes for instance to concentrate vertices at parts of high curvature. Also one can define on simplex meshes shape memory regularization as well as curvature continuous smoothing without shrinkage. On triangulations, the most common way to smooth a shape is to use the mean-curvature flow [START_REF] Chopp | Flow under curvature: singularity formation, minimal surfaces, and geodesics[END_REF] which entails shrinkage.

Simplex Meshes vs. Parametric Surfaces

Parametric surfaces, such as B-splines, the Non-Uniform Rational B-Spline (NURBS), or Finite Element Model (FEM), can provide C k continuity (k > 0) everywhere while for simplex meshes only C 0 continuity is enforced. Therefore, they describe smoother surfaces which can be important for some applications (e.g., rapid prototyping, physical simulation). Also, differential quantities such as surface normals or curvatures can be described everywhere, which is important for shape analysis (e.g., extraction of feature points). However, using those parametric surfaces brings additional complexity especially when the surface cannot be parameterized by a single planar patch (i.e., surface of planar, cylindrical and toroidal topology). Indeed, even for surfaces of spherical topology, parametric surfaces may be complex to handle for instance with the creation of poles (deformable super-quadrics [START_REF] Metaxas | Constrained deformable superquadrics and nonrigid motion tracking[END_REF]). Simplex meshes as well as triangulated surfaces can cope with surfaces of any topology. Furthermore in many image segmentation problems, only the position, volume and distance of the segmented structure are of interest. Therefore, describing surfaces with C 0 continuity is a good trade-off between complexity and accuracy.

Simplex Meshes vs. Level-Sets

Level-setslevel-set rely on an implicit description of the geometry. They can describe manifolds of co-dimension one (i.e., hyper-surface) of any topology but cannot represent surfaces with borders (except if the border is on the edge of the image) such as a planar patch or a cylinder. Simplex meshes can represent a manifold of any dimension and co-dimension. In level-sets, differential quantities (normal, curvatures) may be computed everywhere by filtering the embedding image while in simplex meshes they are only available at vertices. Level-sets can adapt their topologies easily while it is much more difficult for simplex meshes [START_REF] Delingette | New algorithms for controlling active contours shape and topology[END_REF] or triangulations [START_REF] Mcinerney | T-snakes: topology adaptive snakes[END_REF]. Also, the coupled segmentation of structures is easier to implement with level-sets since the distance function is already computed to update their position. Since they are based on an Eulerian scheme, level-sets are not suited to transport a-priori information from a model to the segmented structure. Also, the regularization of level-sets is usually local and based on the mean curvature flow which entails shrinkage.

Statistical Models of Shape and Appearance

As presented in Sec. 12.2.1, shape priors can be employed as stabilizing forces during the evolution of a deformable model. With strong shape priors, it is possible to enforce smooth shape changes around a defined template. Apart from the smoothness constraint, vertices generally move independently from one another, driven by their individual image forces. For many structures in the human body, it is possible to constrain the possible variations further, creating more specific models while still maintaining a good ability to adapt to individual instances. The higher the specificity of a model, i.e. the fewer non-plausible shapes it produces, the more robust it can segment images in the presence of noise and artifacts.

To determine plausible and non-plausible shapes for a certain structure requires either an artificial parametric model (which has to be hand-crafted each time), or a large amount of examples from which this information can be extracted automatically. The SSM employs the latter approach; it generally uses second order statistics to extract the principal modes of variation from a set of example shapes. These modes can be used instead of the internal forces of deformable models. In addition, it is also possible to employ statistics to generate more specific image forces. As for shape modeling, the required information is extracted from a set of example images. Using both techniques together allows to develop very specific and accurate segmentation algorithms for medical applications. A more detailed discussion of the techniques and approaches is given in the comprehensive review by Heimann and Meinzer [START_REF] Heimann | Statistical shape models for 3D medical image segmentation: a review[END_REF].

Shape Representation

As mentioned in Sec. 12.2.6, there exist several different representations for shapes. In the following, we assume the most common and simplest representation: a cloud of points densely sampled on the surface of each shape. Thus, each training shape can be represented as a single vector x that concatenates coordinates from all k points on the surface

x = (x 1 , y 1 , z 1 , . . . , x k , y k , z k ) T (12.3)
Connecting edges between neighboring points leads to a mesh structure. Although these connections are not required for a Point Distribution Model (PDM) [START_REF] Cootes | Active shape models -their training and application[END_REF], they are essential to define inside, outside, and normal vectors for segmentation applications.

Point Correspondence

To be able to perform statistical analysis on a set of training shapes, the individual points representing each training shape have to correspond over the entire set of examples (Fig. 12.10). This means that if a point is placed on a certain anatomic landmark in one example, it should be placed on the same landmark in all other shapes (and the number of sample points is necessarily equal for all training shapes). Thus, all points do represent a kind of landmark, and the two terms are often used synonymously in shape modeling.

Setting landmarks manually on 3D shapes is an extremely tedious and time-consuming process. Additionally, results are not reproducible, as different experts will label different landmarks, especially in areas without discernible features. For this reason, methods to automatically determine correspondences are of paramount importance to shape modeling, and the quality of extracted correspondences has a large impact on the quality of the resulting shape model. The variety of approaches to determine correspondences over a set of shapes can be classified according to the type of registration used to match one shape to the other. In the following sections, we briefly present the most important approaches.

Mesh-to-Mesh Registration

Assuming that the original training shapes are available as meshes, a straightforward solution to the point correspondence problem is to register the different surfaces to each other. One training shape is selected as template and matched to all others by a generic surface matching method as Iterative Closest Point (ICP) [START_REF] Besl | A method for registration of 3-D shapes[END_REF] or Softassign Procrustes [START_REF] Rangarajan | The Softassign Procrustes matching algorithm[END_REF]. For each point of the template, these algorithms deliver the closest point in the target mesh after both have been aligned by an optimal similarity transform. However, this solution has a number of draw-backs. Firstly, the selection of the template induces a bias in the process; different templates usually lead to different correspondences. Secondly, possible corresponding points are limited by the resolution of the target mesh: If this mesh is sampled densely enough, the approach can introduce considerable artificial variation in the resulting SSM. Finally, the most serious problem is that the employed similarity transform may be insufficient to match the template to the other training shapes if the training set shows large geometric variations. In this case, the approach leads to invalid correspondences and the resulting statistical model is not representative of the training data.

An alternative approach is to use non-rigid registration to match the training surfaces. For this, a limited number of matching landmarks is determined in advance, either by manual labeling or by automatic feature matching. These sparse point sets can then be used to initialize a non-rigid registration as the thin-plate-spline deformation, which produces exact matches for the predefined landmarks. Correspondences for the remaining points can be determined simply by Euclidean distance. However, to avoid topological problems it is safer to regularize this matching, e.g., by mesh relaxation [START_REF] Lorenz | Generation of point-based 3D statistical shape models for anatomical objects[END_REF].

Mesh-to-Volume Registration

Most training shapes for medical imaging applications originate from volumetric image modalities as CT or MRI. In these cases, it is possible to acquire point correspondences by registering a deformable mesh to each volumetric training image. The basic idea is that the positions to which vertices of the deformable mesh converge are corresponding over the entire training set. Obviously, the deformable mesh has to adapt accurately to each shape, which is difficult when using one of the typical image forces (Sec. 12.2.2) on intensity images. Therefore, the deformable mesh is mostly used with the segmented, binary volumes, which provide a hard feature for the image force. The geometry for the deformable mesh can be extracted from an arbitrary instance of the training set. However, this procedure can introduce a bias in the corresponding detection. A technique to minimize this bias has been proposed in [START_REF] Zhao | A novel framework for automated 3D PDM construction using deformable models[END_REF], where intermediate templates are generated for difficult cases.

Volume-to-Volume Registration

Similar to registering a deformable mesh to all training volumes, it is also possible to employ a volumetric atlas for that purpose. The atlas can either consist of a single volume or can be built from a collection of images. After registering the atlas to an instance of the training set, the resulting deformation field is used to propagate landmarks on the atlas to the respective training shape. As with mesh-to-volume registration, the atlas registration is usually conducted on the segmented volumes to increase accuracy. To provide the required flexibility, non-rigid transforms are essential for the registration; popular methods include B-splines or thin-plate-splines deformations. As usual, a multi-resolution approach increases robustness of the matching.

Parameterization-to-Parameterization Registration

A parameterization is a bijective mapping between a mesh and an appropriate base domain [START_REF] Floater | Surface Parameterization: a Tutorial and Survey[END_REF]. For closed 2D contours, this base domain is the circle, and determining correspondences by parameterization is equivalent to matching points by relative arc-length. Depending on the topology of training shapes the respective base domains can change, which makes registration between them impossible. This means that a prerequisite for using this approach is that all training shapes have to feature the same topology. Please note that this is also true for the previously presented registration approaches, with the difference that this constraint is not verified there. For 3D shapes, most parameterization approaches are geared to s, i.e. 2-manifolds without holes or self-intersections. The corresponding base-domain is the sphere. When using unconstrained parameterizations as spherical harmonics [START_REF] Kelemen | Elastic model-based segmentation of 3-D neuroradiological data sets[END_REF], the obtained correspondences are often arbitrary. Therefore, a common approach is to determine a limited set of of matching landmarks in advance (as for the mesh-to-mesh registration) and employ to match the remaining regions of the shape.

Population-Based Optimization

In all the correspondence approaches presented above, matching points are determined in a pair-wise comparison or without a comparison at all (for the parameterization-based approaches). Population-based correspondence methods analyze the complete set of training shapes at once and select corresponding points based on a global cost function. Starting from initially arbitrary landmarks, an optimization algorithm iteratively modifies correspondences until the cost function converges at its minimum. Although this procedure sounds straight-forward and should lead to optimal results, there are a number of challenges involved in the design.

First, as it is difficult to define what good correspondences are, coming up with a practically usable and theoretically sound cost function is not evident. Practically all current approaches rely on a form of Minimum Description Length (MDL) as cost function, which was pioneered for correspondence optimization by Davies [START_REF] Davies | 3D statistical shape models using direct optimization of description length[END_REF]. MDL implements the principle of Occam's razor that a simpler model is better than a complex one. Translated to SSMs, this means that a model is better when it features fewer and more compact modes of variation. As the cost function for the full MDL is computationally very expansive, a simplified version [START_REF] Thodberg | Minimum description length shape and appearance models[END_REF] is often used. Here, costs F are directly calculated from the eigenvalues λ of the shape model

F = m L m with L m = 1 + log(λ m /λ cut ) for λ m ≥ λ cut λ m /λ cut for λ m < λ cut (12.4)
where λ cut describes the amount of noise in the dataset.

A second challenge is to adjust the landmark positions in the course of optimization. The difficulty lies in the fact that the optimization problem is heavily constrained, as individual landmarks cannot be moved independently from their neighbors (to avoid flipping triangles in the model mesh). Additionally, all 3D landmarks have to remain on the 2D manifolds of the surfaces of the training shapes. This is either accomplished by employing parameterizations to modify landmark positions [START_REF] Davies | 3D statistical shape models using direct optimization of description length[END_REF][START_REF] Heimann | 3D Active shape models using gradient descent optimization of description length[END_REF] or by particle approaches [START_REF] Cates | Shape modeling and analysis with entropy-based particle systems[END_REF].

Construction of Statistical Shape Models

Once point correspondences are known, each training shape i can be represented as coordinate vector x i in (12.3). Before applying the dimensionality reduction that delivers the predominant modes of variation, the individual shapes first have to be aligned into the same coordinate system.

Alignment

The commonly used procedure to align a set of training shapes is the Generalized Procrustes Analysis (GPA) [START_REF] Goodall | Procrustes methods in the statistical analysis of shape[END_REF]. The underlying Procrustes analysis determines the optimal similarity transform to minimize the sum of squared distances between two corresponding point sets. GPA employs this method iteratively to match a whole set of shapes to their unknown mean. An additional step for constructing a shape model is a rescaling of the thus aligned shapes. This is required since the minimization of squared distances introduces non-linearities in the shape space, which are difficult to capture by commonly used techniques. A straight-forward solution is to scale each training shape i with 1/(x i • x), where x is the mean as determined by GPA.

Dimensionality Reduction

Statistical shape modeling is essentially a problem of dimensionality reduction:

The input data of s training shapes x i (with the dimension of x i generally in the thousands) should be reduced to a limited number of modes of variation. This task is usually solved by Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF]. The first step is to calculate the mean shape x of the training set

x = 1 s s i=1 x i (12.5)
Then, the covariance matrix S of the input data is computed by:

S = 1 s -1 s i=1 (x i -x)(x i -x) T (12.6)
An eigendecomposition of S yields the eigenvectors φ m (representing the principal modes of variation) and the corresponding eigenvalues λ m (indicating the variance per mode). Sorting all modes from largest to smallest variance, the first c modes are employed to model the distribution, while the remaining modes are discarded. Thus, all valid shapes can be approximated by a linear combination

x = x + c m=1 b m φ m (12.7)
with b as the vector of shape parameters. To constrain variation of the model, b is usually bound to certain limits, either by constraining |b m | < 3λ m for each mode individually or by constraining b globally to lie within an appropriate ellipsoid. Typically, c is chosen so that the accumulated variance c m=1 λ m reaches a certain percentage (e.g. 95 %) of the total variance. Fig. 12.11 shows the first three principal modes of variation as extracted from a collection of liver shapes.

Modeling Object Appearance

In Sec. 12.2.2, we already described possible image forces for deformable models. As we assume to have training data available for statistical modeling, we now focus on strong appearance priors.

Boundary-Based Features

The most popular appearance models for SSMs are related to the abovementioned intensity profile. At each landmark of the shape model, a profile is sampled in normal direction and compared to a model built from training data at the same landmark. This implies that one shape model commonly features a large number of appearance models. The employed profiles are not limited to plain intensity alone; gradient profiles and normalized intensity or gradient profiles are also popular choices. Different features can also be combined in the same model. Each appearance model is built similar to the shape model, i.e. mean profile and covariance matrix are generated. Using the covariance matrix, the Mahalanobis distance to profiles sampled during the search can be determined. For Gaussian distributions, this represents a robust measure of the goodness of fit [START_REF] Cootes | Using grey-level models to improve active shape model search[END_REF].

In cases where the distribution of profiles is essentially non-Gaussian (e.g., because the modeled object is adjacent to several tissues of different appear- ance), non-linear appearance model can deliver better results. A straightforward method to employ a non-linear model is using a k-Nearest-Neighbor (k-NN) classifier [START_REF] De Bruijne | Adapting active shape models for 3D segmentation of tubular structures in medical images[END_REF]. To train it, both true and false examples of boundary profiles must be used. The required false examples can be generated by shifting profiles from their known positions towards the inside and outside of the Object of Interest (OOI).

Region-Based Features

As mentioned in Sec. 12.2.2, region histogram priors can be used to evolve a deformable model in order to match a certain distribution of intensities.

For statistical models, it is straight-forward to model this distribution from training data. However, there are a number of much more specific features that can be used in this case. The most consequent approach for learning object appearance uses the entire texture of the OOI [START_REF] Cootes | Active appearance models[END_REF]. To combine textures from all training images, the respective objects are morphed to the same shape (the mean shape of the SSM). Subsequently, the mean texture and its covariance matrix can be calculated as usual. A drawback of this approach is the high memory footprint of the appearance model, especially for 3D models where a volumetric texture is used. One solution is to model the texture only in certain parts of the object, e.g. around the boundary. Modeling local regions inside and outside of the object of interest opens up more possibilities regarding the search algorithm for model evolution (Sec. 12.3.5), as the goodness of fit can be determined locally (like when using profiles). In addition to the above mentioned techniques of histogram matching and texture matching, local appearance can also be modeled by multi-scale derivatives, Gabor filterbanks, wavelets, etc. In comparison to boundary-based features, region-based features are often more robust as they use more image data to determine the goodness of fit. At the same time and for the same reason, they are also computationally more demanding.

Local Search Algorithms

As mentioned in Sec. 12.2.5, the first step to employ a deformable model for segmentation is to initialize it for the respective image. This initialization places the model close to the object of interest (e.g., by atlas registration), so that the model can iteratively evolve towards the correct contour. In the following, we present several techniques to conduct this subsequent evolution (the local search) for statistical models of shape and appearance.

Active Shape Model

The Active Shape Model (ASM) was introduced by Cootes et al. [START_REF] Cootes | Active shape models -their training and application[END_REF]. It is one of the most popular model-based approaches for segmentation of medical images. An instance y of the model in an image is defined by a similarity transform T and the shape parameter vector b y = T (x + Φb) (12.8) where Φ = (φ 1 . . . φ c ) is the matrix of eigenvectors (12.7). To evolve the shape over time, the optimal displacement dy p is determined for each landmark point separately by querying the goodness of fit for the respective appearance model at various positions along the . The update then consists of two steps: First, transform T is updated by the result of the Procrustes match of the model to y + dy p . This leads to new residual displacements dy s . Second, shape parameters b are updated by transforming dy s into model space and then applying a projection into shape space

db = Φ T T -1 (dy s ) (12.9)
where T is equal to T without the translational part. The resulting adjustments db are added to b. Usually, b is restricted to certain limits to maintain a valid shape (Sec. 12.3.3). The presented two steps are repeated iteratively until the model converges at the best local fit in the image.

To increase , especially against incorrect initialization, the search is commonly conducted in a multi-resolution fashion. Contrary to the multiresolution approach presented for simplex meshes (Sec. 12.2.3), the DOF usually remains unchanged. Instead, appearance models are constructed for a number of down-sampled versions of all training images, and the search is started on a down-sampled version of the original image. Here the search radius is larger, since the goodness of fit is evaluated at points at a greater distance from the original contour. When the model evolution converges on the rough scale, appearance model and search radius are switched to the next (more detailed) versions, until the original resolution is reached.

Another approach to increase robustness is to decrease the effect of outliers, i.e. wrongly detected optimal displacements for specific landmarks. There are two basic techniques for this problem. Firstly, outliers can be detected (e.g. by comparing suggested landmark positions to collected statistics) and corrected before the shape is updated. Secondly, the update of transform and shape parameters can be weighted and each suggested landmark position receives a reliability weighting in each step.

Active Appearance Model

The Active Appearance Model (AAM) [START_REF] Cootes | Active appearance models[END_REF] belongs to the class of generative models, i.e., they can generate realistic images of the modeled data. This is accomplished by storing a complete texture model, consisting of the mean and the principal modes of variation, in addition to the shape model. However, AAMs are more than a shape model with region-based features, they employ a specific search method which is completely different from the ASM approach. Moreover, shape and appearance parameters are combined into one linear system, where shape x and appearance g are described by a common parameter vector c

x = x + Φ s W s Q s c g = ḡ + Φ g Q g c (12.10)
Here, Φ s and Φ g are the independent eigenvector matrices of shape and appearance model, respectively, and W s is a diagonal weight matrix for the shape parameters. Q = Qs Qg is the eigenvector matrix of the combined shape and appearance parameters. It is the result of a PCA on the independent parameters b = Wsbs bg . An instance of the model in an image is defined by a similarity transformation T and the combined shape-appearance parameters c. In the following, we name the latter ones simply parameters p. To evaluate the goodness of fit, the image texture is warped to the mean shape and normalized, resulting in g s . With the modeled appearance g m = g from (12.10), the residuals are given by r(p) = g sg m , and the error by E = r 2 .

The key idea of AAM search is to assume a constant relationship between texture residuals r(p) and parameter updates dp over the entire search dp = -Rr(p) (12.11) The success of this optimization scheme largely depends on the derivative matrix R. In the first presentation of AAMs, R was computed using multivariate linear regression on a large number of simulated disturbances of the training images. Later, regression was replaced by numeric differentiation, claimed to be both faster and more reliable.

The main challenge when employing AAMs for medical image segmentation is the enormous amount of data that has to be captured by the model. For 3D models, appearance is generally modeled as a 3D texture, which quickly leads to prohibitively large equation systems. Therefore, texture resolution generally has to be scaled down. It it noteworthy that the AAM search procedure is not limited to region-based features. Indeed, the general idea of determining parameter updates from residuals can also be used with alternative appearance features as described in Sec. 12.3.4.

Conclusion

Model-based approaches to segmentation are arguably the most robust methods when image data is noisy or includes artifact. Therefore, they are a prime choice for many applications in medical image analysis. The deformable simplex meshes and point-based statistical models presented in this chapter represent state-of-the-art techniques for this purpose.

Comparing both methods, an obvious difference is that building a shape model requires a large collection of training images (typically 20-50 as a minimum for 3D models), while a simplex mesh with strong shape prior can be based on a single example. Another point is that the additional constraints of SSMs generally result in a higher robustness, but at the same time limit the accuracy of the final result. Statistical models are therefore often applied as part of a coarse-to-fine strategy that features an additional refinement step afterwards. As already hinted to in Sec. 12.2.3, it is also possible to combine both techniques in one framework. That way, the strict shape constraints can gradually be lifted and the final free deformation ensures an optimal fit to the data.

List of Acronyms

Fig. 12 . 1 .

 121 Fig. 12.1. Taxonomy of deformable models. The classification scheme of deformable models for medical image segmentation is based on their geometric representation.

Fig. 12

 12 Fig. 12.5. Deformable simplex meshes.The meshes are driven by the matching of intensity profiles; Upper row: reference CT image scan (left) and head segmented with a simplex mesh (right). Lower row: an intensity profile is stored at each vertex as a reference and this mesh is used to register the head of the same patient on an MRI: mesh position before (left) and after rigid registration (right).

Fig. 12

 12 Fig. 12.6. Coarse-to-fine strategy of 3D mesh reconstruction. The coarse to fine strategy runs from clouds of points: the different stages of the deformations with increasing deformation complexity are shown ending with a globally constrained deformation (λ = 0.3 in panel f).

Fig. 12 .

 12 Fig. 12.7. 3D+T simplex mesh. The array of meshes has the same topology but different geometries.

Fig. 12 . 8 .

 128 Fig. 12.8. Segmentation of the left ventricle endocardium. Segmentation from 3D+T echo-cardiographic images by 3D+T simplex meshes.

Fig. 12 . 9 .

 129 Fig. 12.9. Joint segmentation of bladder and prostate. The bladder (red) and prostate (green) are shown in a CT image. Left: no inter-penetration constraints are used which causes a leakage; Middle: with inter-penetration constraint and some statistical shape constraint on the prostate, the segmentation is improved despite the low contrast between structures; Right: 3D view of the two organs.

Fig. 12 .

 12 Fig. 12.10. Correspondences for a set of liver shapes. Corresponding areas are marked by colors and a sparse coordinate grid.

Fig. 12 . 11 .

 1211 Fig. 12.11. Statistical shape model of the liver. Panels (a,b), (c,d), and (f,g) show the first three modes of variation for a SSM of the liver with mean shape (d).

Table 12 .

 12 1. Geometric transforms and constrains. The table shows geometric ztransforms that can constrain the deformation of a simplex mesh.

	Transform	Number of DOF	Transformation type
	Rigid	6	global
	Similarity	7	global
	Affine	12	global
	B-spline	3 × number of control points	restricted to a bounding box
	Axis-symmetric	3 × number of points on the central line	restricted to the mesh
	PCA	6 + number of modes	restricted to the mesh
	may still have too many DOF to prevent falling into a local mimimum of the
	functional.

Model-based Segmentation


Tobias Heimann
email: tobias.heimann@inria.fr

François Chung

Hans Lamecker

Hervé Delingette

Subject-Specific Ligament Models: Towards Real-Time Simulation of the Knee Joint

Keywords:

come

Introduction

The human knee joint has been extensively studied in the field of biomechanics, both due to its complex anatomy and its high clinical relevance. It has to cope with high loads and large displacements during many daily activities. The knee ligaments, notably the two cruciate ligaments and the two collateral ligaments, play a key role in ensuring stability of the joint [START_REF] Woo | Biomechanics of knee ligaments[END_REF]. They are also at increased risk during high impact scenarios as occurring e.g. in many sports activities. In case of ligament rupture, kinematics of the knee are severely compromised, and sophisticated surgical interventions are required to restore its function. Virtual simulation of these surgeries for subject-specific cases combined with a subsequent analysis of resulting knee kinematics seem a promising way to increase success rates of interventions [START_REF] Chen | Knee surgery assistance: patient model construction, motion simulation, and biomechanical visualization[END_REF][START_REF] Ramaniraka | Biomechanical evaluation of intra-articular and extra-articular procedures in anterior cruciate ligament reconstruction: a finite element analysis[END_REF].

The finite element method (FEM) has successfully been used to analyze a variety of problems regarding the knee joint [START_REF] Ramaniraka | Biomechanical evaluation of intra-articular and extra-articular procedures in anterior cruciate ligament reconstruction: a finite element analysis[END_REF][START_REF] Ellis | Medial collateral ligament insertion site and contact forces in the ACL-deficient knee[END_REF][START_REF] Fernandez | An anatomically based patient-specific finite element model of patella articulation: towards a diagnostic tool[END_REF]. To achieve realistic simulation of ligaments, sophisticated material models have been presented [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF][START_REF] Pioletti | Viscoelastic constitutive law in large deformations: application to human knee ligaments and tendons[END_REF]. These models implement properties such as the transverse isotropy and non-linear stressstrain curve due to internal collagen fibers, and partly also viscoelastic behavior.

In most works related to real-time simulation, however, material properties of ligaments are highly simplified to maintain the required frame-rate. Often, simple isotropic linear elasticity is used [START_REF] Wu | An improved scheme of an interactive finite element model for 3D soft-tissue cutting and deformation[END_REF]. Picinbono et al. [START_REF] Picinbono | Non-linear anisotropic elasticity for real-time surgery simulation[END_REF] presented an optimized FE model for tetrahedral elements, which also allowed transversely isotropic materials. However, the stress-strain relationship was still linear. A fast non-linear anisotropic material model was presented by Teran et al. [START_REF] Teran | Creating and simulating skeletal muscle from the visible human data set[END_REF], who use the finite volume method to simulate muscles and tendons. This approach is equivalent to FEM for constant strain linear basis functions. Recently, Joldes et al. [START_REF] Joldes | Suite of finite element algorithms for accurate computation of soft tissue deformation for surgical simulation[END_REF] developed a highly efficient framework for computation of soft tissue deformation, which is based on Total Lagrangian and explicit time integration. This framework was also implemented on GPU and used for viscoelastic nonlinear soft tissue simulation by Taylor et al. [START_REF] Taylor | On modelling of anisotropic viscoelasticity for soft tissue simulation: Numerical solution and GPU execution[END_REF].

In this work, we present a material model that incorporates the ligamentspecific properties of transverse isotropy and non-linear stress-strain curve and that is suitable for both explicit and implicit solvers. The model is aimed at real-time simulation and reaches high efficiency due to the use of linear tetrahedral elements, which leads to an optimization of a large part of the involved computations. We also address the problem of generating subject-specific models, which is neglected in many current publications: with transversely isotropic materials for ligaments, the main fiber orientation has to be known for each finite element in the simulation. Techniques as X-ray diffraction or optical diffraction that can be used for experiments with cadavers are not possible in vivo. To employ FEM simulations of the knee for surgery planning, fiber orientations have to be extracted from clinical imaging modalities as MRI or DT-MRI.

The remainder of this paper is organized as follows: in Sec. 2, we present our optimized FE model for ligaments. Subsequently, the method to determine local fiber orientation from segmented geometry is presented in Sec. 3. In Sec. 4, image data and rheological parameters to generate a subject-specific model of the medial collateral ligament (MCL) are described. In Sec. 5, we verify our FE model by comparing it to the established FEBio software 1 and we present our model of the MCL. We close with a discussion of results and future work.

Transversely isotropic hyperelasticity for tetrahedrons

In this section, we present the strain energy for the employed material model and derive its stress tensor (first derivative) and elasticity tensor (second derivative). The latter one is required if the model is used for implicit or quasi-static analysis. Throughout this work, we follow the concepts presented by Delingette [START_REF] Delingette | Triangular springs for modeling nonlinear membranes[END_REF] for triangular finite elements. In particular, the elasticity tensor for each node is assembled by accumulating terms from all adjacent edges, which allows to exploit symmetry to reduce computations.

Properties of tetrahedral elements

The tetrahedral element with linear basis function features a constant strain, which reduces integration to a multiplication by volume. Moreover, the deformation tensor can be calculated efficiently based on the current node positions Q i and the shape vectors of the tetrahedron. Shape vectors D i are defined by the rest configuration of the tetrahedron (with node positions P i ) as:

D i = s V 0 ((P i⊕2 -P i⊕1 ) ∧ (P i⊕3 -P i⊕1 )) (1)
where ∧ denotes the cross product, V 0 is the volume in rest configuration and s = 1 for i = 2, 4, else s = -1. Shape vector D i stands perpendicular to the face opposite node i and points inward. For the many interesting properties of shape vectors, we refer the reader to [START_REF] Delingette | Triangular springs for modeling nonlinear membranes[END_REF]. Here, we note that the deformation gradient tensor ∇Φ can be expressed as:

∇Φ = 4 i=1 Q i ⊗ D i (2)
Consequently, the right Cauchy-Green deformation tensor C is:

C = ∇Φ T ∇Φ = 4 i=1 4 j=1 (Q i • Q j )(D i ⊗ D j ) (3)
The Jacobian J of the deformation gradient tensor ∇Φ can also be expressed as the ratio of volumes in deformed and in rest configuration: J = V /V 0 . We sum up that the signed volume of a tetrahedron is given by:

V = 1 6 ((Q 4 -Q 1 ) • ((Q 2 -Q 1 ) ∧ (Q 3 -Q 1 ))) (4)
where the order of nodes is interchangably as long as the sign is kept. As described in [START_REF] Delingette | Triangular springs for modeling nonlinear membranes[END_REF], forces on node i can be determined directly by F i = -V ∇Φ S D i , where S is the second Piola-Kirchhoff tensor used for large deformations. Using the weak formulation for finite elements [START_REF] Bower | Applied Mechanics of Solids[END_REF], the required derivative ∂S/∂C results in very complex expressions (see e.g. [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF]). Therefore, we employ the Rayleigh-Ritz approach and derive the required expressions directly from the discretized strain energy, using individual node positions Q i .

Strain energy

We base our strain energy W on the terms proposed by Weiss et al. [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF] for incompressible, transversely isotropic materials. However, instead of using an Augmented Lagrangian method to ensure incompressibility, we employ a simple penalty approach based on the bulk modulus of the material. The energy is composed of three different terms:

W = F G (I 1 ) + F C (I 4 ) + F V (J) (5)
F G describes the contribution of the isotropic ground substance matrix, F C describes the contribution from the collagen fibers, and F V is the penalty term for volume changes. The underlying invariants are defined as:

I 1 =J -2 3 tr C I 4 =J -2 3 aCa = λ 2 (6) J =V /V 0
Here, I 1 is a standard invariant of C for nearly incompressible materials [START_REF] Bower | Applied Mechanics of Solids[END_REF]. I 4 arises from the anisotropy of the material, using main fiber orientation a [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF]. λ = l/l 0 is the stretch along this direction. Note that in our formulation, I 1 and I 4 remain constant under a pure volume change. Thus, F G and F C represent a pure deviatoric response and F V a pure dilational response.

Derivations for the Jacobian

Starting with the formula for the volume of a tetrahedron (4), the first derivative of J w.r.t. point Q i can be determined as:

∂J ∂Q i = 1 6V 0 (b ∧ c) T (7)
Here, b = Q j -Q l and c = Q k -Q l are edge vectors of the triangle opposing Q i . Thus, indices j, k, l depend directly on i and can be stored in a simple table.

Using edge vector c, the second derivative of J is:

∂ 2 J ∂Q j ∂Q i = 1 -δ ij 6V 0   0 -c z c y c z 0 -c x -c y c x 0   (8)
Here, δ ij is the Kronecker delta and matrix elements c are components of the edge opposing points Q i and Q j .

Derivations for isotropic ground substance

Following Weiss et al. [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF], we use a Neo-Hookean material with shear modulus µ to model the response of the isotropic ground substance matrix:

F G = µ 2 (J -2 3 tr C -3) (9)
Stress tensor The individual terms required for the first derivative are:

∂J -2 3 ∂Q i = - 1 9V 0 J -5 3 (b ∧ c) T , ∂tr C ∂Q i = 2 4 n=1 (D n • D i )Q T n ( 10
)
The first derivative can then be calculated by:

∂F G ∂Q i = µ 2 J -2 3 ∂tr C ∂Q i + ∂J -2 3 ∂Q i tr C (11)
Elasticity tensor Again, we start with the individual terms of F G :

∂ 2 J -2 3 ∂Q j ∂Q i = 10 9 J -8 3 ∂J ∂Q j T ∂J ∂Q i - 2 3 J -5 3 ∂ 2 J ∂Q j ∂Q i , ∂ 2 tr C ∂Q j ∂Q i = 2(D i • D j )1
(12) Using these results, the second derivative is:

∂ 2 F G ∂Q j ∂Q i = µ 2   ∂ 2 J -2 3 ∂Q j ∂Q i tr C + ∂tr C ∂Q j T ∂J -2 3 ∂Q i + ∂J -2 3 ∂Q j T ∂tr C ∂Q i + J -2 3 ∂ 2 tr C ∂Q j ∂Q i   (13)
2.5 Derivations for collagen fiber family F C is a scalar function governing the non-linear behavior of the material along the main fiber orientation. We employ the formulation by Weiss et al. [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF] and explicitly model the toe region of the stress-strain curve, which represents the successive recruitment of collagen fibers in ligaments under lower strains λ < λ * .

∂F C ∂λ =0 , λ <1 ∂F C ∂λ =C 3 (exp(C 4 (λ -1)) -1) , λ <λ * (14)
∂F C ∂λ =C 5 λ + C 6 , λ ≥λ *
Stress tensor First, we note that:

∂F C ∂Q i = ∂F C ∂λ ∂λ ∂I 4 ∂I 4 ∂Q i = ∂F C ∂λ 1 2λ ∂I 4 ∂Q i (15)
This leaves us with finding the derivative for I 4 . As we can use the result for ∂J -2 3 /∂Q i from (10), we continue with:

∂aCa ∂Q i = 2(a • D i ) 4 n=1 (a • D n )Q T n (16)
We now know all the individual terms required for the derivative of I 4 :

∂I 4 ∂Q i = J -2/3 ∂aCa ∂Q i + ∂J -2 3 ∂Q i aCa (17)
Using this in (15) delivers the stress tensor for F C .

Elasticity tensor

The second derivatives of the stress-strain curve are:

∂ 2 F C ∂λ 2 =0 , λ <1 ∂ 2 F C ∂λ 2 =C 3 C 4 exp(C 4 (λ -1)) , λ <λ * (18)
∂ 2 F C ∂λ 2 =C 5 , λ ≥λ *
We continue from [START_REF] Weiss | Computational modeling of ligament mechanics[END_REF], which gives us:

∂ 2 F C ∂Q j ∂Q i = ∂ 2 F C ∂λ 2 ∂λ ∂I 4 2 + ∂F C ∂λ ∂ 2 λ ∂I 2 4 ∂I 4 ∂Q j T ∂I 4 ∂Q i + ∂F C ∂λ ∂λ ∂I 4 ∂ 2 I 4 ∂Q j ∂Q i (19) The second derivative of λ w.r.t. I 4 is simply ∂ 2 λ ∂I 2 4 = -1 4 √ I4
3 . Next, we require the second derivative of I 4 :

∂ 2 I 4 ∂Q j ∂Q i = ∂ 2 J -2 3 ∂Q j ∂Q i aCa+ ∂aCa ∂Q j T ∂J -2 3 ∂Q i + ∂J -2 3 ∂Q j T ∂aCa ∂Q i +J -2 3 ∂ 2 aCa ∂Q j ∂Q i (20)
The only new term in this equation is the second derivative of aCa, which is:

∂ 2 aCa ∂Q j ∂Q i = 2(a • D i )(a • D j )1 (21)
Now all terms required for ( 19) are known and we can calculate the elasticity tensor for F C .

Derivations for volume conservation

The penalty term controlling response of an element to volume changes is the same as for a nearly incompressible Neo-Hookean model [START_REF] Bower | Applied Mechanics of Solids[END_REF]:

F V = K 2 (J -1) 2 (22)
with K as the bulk modulus of the material.

Stress tensor

The derivation of the first derivative is straight-forward:

∂F V ∂Q i = K(J -1) ∂J ∂Q i ( 23
)
Elasticity tensor The second derivative of the penalty term is:

∂ 2 F V ∂Q j ∂Q i = K ∂J ∂Q j T ∂J ∂Q i + (J -1) ∂ 2 J ∂Q j ∂Q i (24)

In-situ stress

When joints as the knee are in neutral position, the involved ligaments are usually still under tension [START_REF] Weiss | Computational modeling of ligament mechanics[END_REF]. Since this tension stabilizes the joint even when no muscles are active, the so-called in-situ stress plays an important role in joint kinematics. To achieve realistic simulation results, the employed FE model should thus include pre-stressed materials. In our model, we incorporate in-situ stress by increasing the magnitude of direction vectors in the fiber orientation field. By default, fiber orientation a is represented by a unit vector. If we use a modified orientation vector a ′ = (1 + c)a, the material is under a strain of c even for an identity deformation gradient tensor.

Estimating subject-specifc fiber orientations

We assume that the subject-specific geometry of the ligament to be simulated is given by a high-resolution mesh. The attachment areas to the two bones the ligament connects are specified by two sets of vertices of this mesh, A 1 and A 2 . The first step in our approach consists of finding the centerline between these two areas. This line provides a solid estimate of the main fiber orientation.

Centerline extraction

Centerline extraction is a common problem in medical image analysis, employed e.g. for vessel tracking or for generating fly-through paths for virtual endoscopy. For a brief overview of existing approaches see Bouix et al. [START_REF] Bouix | Flux driven automatic centerline extraction[END_REF]. As we can calculate the line start and end positions as the average positions of A 1 and A 2 , we opt for a method based on minimal cost paths. Among all possible paths P (s) connecting both points, we search the one that minimizes a specific metric C. We include two components in C: first, the path should be as short as possible, and second, it should run at maximum distance from the mesh surface (i.e. as close to the medial sheet of the shape as possible). Both components are weighted by α ∈ [0, 1] such that α = 0 results in the shortest possible path (a straight line) and α = 1 in a path around the geometric center. This leads to the following metric:

C(P ) = 1 0 ∂P ∂s α d 2 (P (s)) + (1 -α) ds (25)
where d(x) is the closest distance to the mesh surface from point x.

The search for the minimal cost path is set in a discrete setting. First, the surface mesh is rendered to an image volume and the corresponding distance transformation is calculated. Then, a graph is constructed with all object voxels as nodes and the 26-neighborhood for each voxel as edges. A Dijkstra search finally delivers the minimum cost path that serves as centerline for the ligament.

Local fiber orientation

To calculate fiber orientations for all positions inside the segmented ligament mesh, we use the gradient of a diffusion field between centerline and mesh boundary. The diffusion field D is constructed such that all voxels on the centerline have a value of 0 and all voxels on the mesh boundary a value of 1. Generation of this field is a classical Dirichlet problem. In our current implementation we use the following geometric method: for each point x within the mesh, the closest point on the centerline C(x) is determined. A ray is cast from C(x) through x until it hits the mesh boundary at B(x). We then set

D(x) = |x -C(x)|/|B(x) -C(x)|.
We assume that fibers are oriented parallel to the boundary for positions on the boundary and parallel to the centerline for positions on the centerline. In between, fiber directions are interpolated so that they run perpendicular to the gradient ∇D, i.e. on an isosurface of the diffusion field. To determine the orientation vector a(x), the direction of the centerline at C(x) is projected onto the isosurface at x and normalized. Finally, orientation vectors for all x inside a specific tetrahedral element are averaged to yield a for this element.

Material for subject-specific ligament model

As proof of concept for the presented methods, we simulate the medial collateral ligament (MCL) in a subject-specific knee joint.

Image data

MRI data was used as source for the geometry of the simulation. An institutional medical-ethics committee approved the image acquisition and the subject gave his written consent. The employed dataset covers the right knee of a healthy male subject (age: 38) in full extension. It was acquired using a 1.5T Siemens MRI device with Flair Sagittal 3D T1 protocol and 256 × 256 matrix. Image resolution was 0.5 × 0.5 mm with a slice thickness of 1.0 mm.

Rheological parameters

While it is possible to extract subject-specific geometry for a living subject, we can think of no way to determine subject-specific material properties. For our simulation of the MCL, we thus employ the rheological parameters presented for this structure by Weiss et al. [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF]: µ=27.7 MPa, λ * =1.03, C 3 =2.07 MPa, C 4 =61.44, C 5 =640.7 MPa. In-situ strain was set to 2.9%. To achieve near incompressible behavior, the bulk modulus was set to K = 1000.

Experiments and results

Verification

To verify our implementation of the material model, we conducted a number of tests on a synthetic cube shape of 1mm edge length. The employed mesh consisted of 4316 nodes and 20737 tetrahedral elements. Material parameters were set to the values given in Sec. 4.2. Constraints were used to fix positions of all nodes on the lower side of the cube. To all nodes on the upper side of the cube, we applied constant forces pointing upwards (tensile strain) or downwards (compressive strain). An implicit solver with velocity damping was run until convergence and the resulting deformed mesh was exported.

The scenario was then simulated with the FEBio software, using the same mesh and identical boundary conditions. The material model in FEBio was set to transversly isotropic Mooney-Rivlin with the same parameters as above. Volume conservation was enforced using the Augmented Lagrangian approach. The solution was obtained for the quasi-static case with the default BFGS solver.

We compared our results to the FEBio reference solution for a variety of cases. All positioning errors are expressed in percent of the edge length of the cube in Table 1. Results of experiments on the synthetic cube with fibers oriented parallel to tensile forces. For each load, the resulting average fiber strain and volume change is given. From comparison to the FEBio solution, we determined the average positioning error for all nodes on the top surface of the cube as well as the maximum error. 1 lists results for tensile forces typically encountered in ligaments under physiological conditions. In Fig. 1, the corresponding node positioning errors are visualized. Table 2 lists results for cases that usually do not occur in ligaments (as loading perpendicular to fiber orientation) and also compares our results to the FEBio solution without Augmented Lagrangian. Figure 2 again visualizes the corresponding node positioning errors.

Load

MCL simulation

Geometry of femur, tibia and MCL was segmented manually from the MRI dataset described in Sec. 4. For the centerline extraction, we set α = 0.1. Figure 3 shows the segmented ligament mesh and the resulting centerline. From the surface, we generated a volumetric mesh consisting of 1193 nodes and 5911 tetrahedra. Using the presented rheological parameters for the ligament and rigid bones, a dynamic simulation was set up with the SOFA software [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF]. An Euler implicit scheme was used for time integration, and a modified conjugate gradient solver [START_REF] Baraff | Large steps in cloth simulation[END_REF] was employed to solve the system of equations. While the tibia was fixed, the femur was moved manually to different flexion angles. Figure 4 shows the resulting strain energy within the MCL for these cases. On an Intel Core 2 Duo notebook with 2.2GHz CPU the obtained framerate was 15 frames/second (including simulation and visualization).

Discussion

The comparison of our FE model to FEBio results in average node positioning errors < 1% for the intended application, which we regard as sufficient accuracy. As a comparison to FEBio without Augmented Lagrangian reveals, the discrepancies are mainly due to the different approach to volume conservation. Our method to estimate fiber orientations requires only minimal interaction for the definition of attachment areas (which are required in any case) and delivers mostly convincing results. However, some problems occur at the attachment areas, where orientation cannot be determined reliably when the centerline does not start at the mesh boundary. While the obtained framerate is not yet suitable for real-time simulation, it is nonetheless much faster than classic FE solvers. Indeed, the largest challenge we see for future virtual surgeries is not the performance of the FE model, but the generation of accurate subject-specific geometries: ligaments are hardly visible in standard MRI data, and small variations in the segmentation can affect simulation results considerably. As a next step, we plan to validate the presented methods against experimental studies on real ligament tissue.

Fig. 1 .

 1 Fig. 1. Node positioning errors to FEBio reference solution for tensile loading, colorcoded from 0 to 3%. From left to right, applied tensile forces per node are 0.05N, 0.1N, and 0.2N, respectively. Fibers are oriented parallel to forces.

Fig. 2 .

 2 Fig. 2. Node positioning errors for alternative test cases. Left: Tensile forces of 0.05N/node, but with fibers oriented perpendicular to forces. Center: Compressing forces of 0.05N/node. Right: Same scenario as in Fig. 1 center, but compared to FEBio without Augmented Langrangian. Note the different error scales (up to 6% for left and center, up to 1% for right).

Fig. 3 .

 3 Fig. 3. Segmented ligament mesh together with the original MRI data (left) and the extracted centerline connecting the attachment areas of femur and tibia (right).

Fig. 4 .

 4 Fig. 4. Simulation of strain in the MCL at 0, 45, and 90 degrees of knee flexion. Strain energy is color-coded from green to red, the lines indicate local fiber orientation.

Table 2 .

 2 /node Total load Avg. fiber strain Vol. change Avg. error Max. error Additional experiments on the synthetic cube. The first two cases represent unusual load scenarios, the last two (marked with * ) are compared to the FEBio solution without Augmented Lagrangian.

	0.02 N	5.8 N	+1.6%	+0.2%	0.1%	0.4%
	0.05 N	14.5 N	+3.0%	+0.5%	0.3%	1.0%
	0.10 N	28.9 N	+4.9%	+1.1%	0.5%	1.9%
	0.20 N	57.8 N	+8.8%	+2.2%	0.9%	3.4%
			Set-up Vol. change Avg. error Max. error
	Fibers perpendicular to tensile load (0.05 N)	+0.7%	3.9%	5.8%
	Fibers parallel to compressive load (0.05 N)	-0.5%	2.4%	5.6%
	Fibers parallel to compressive load (0.10 N) *	+1.1%	0.3%	0.6%
	Fibers perpendicular to tensile load (0.05 N) *	+0.7%	0.4%	0.5%
	rest configuration. Table

FEBio is available at the website of the Musculoskeletal Research Laboratories of the University of Utah: http://mrl.sci.utah.edu/software.php

Acknowledgment This work is supported by the EU Marie Curie project 3D Anatomical Human (MRTN-CT-2006-035763).