Jonathan Boisvert
email: <jonathan.boisvert@polymtl.ca>

Farida Cheriet

Xavier Pennec

Hubert Labelle

Nicholas Ayache

Principal Deformations Modes of Articulated Models for the Analysis of 3D Spine Deformities

Keywords: Shape Analysis, Articulated Models, Spinal Deformities, Scoliosis, 3D Reconstruction, Surgical Classifications

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Articulated models are collections of objects linked together so that the relative motion of any two parts is rigid (combination of a 3D rotation and a translation). This family of models can be used to describe a wide variety of objects. A popular application is human gait analysis from video sequences. In this type of application, human body motion is characterized by articulations states and about thirty degrees of freedom are usually necessary to obtain convincing results. This number of degrees of freedom (DOF) makes model estimation from video sequences challenging from a computational perspective. However, probabilistic methods can be used to discard unlikely positions and to deal with occlusions [START_REF] Deutscher | Articulated body motion capture by stochastic search[END_REF][START_REF] Green | Quantifying and recognizing human movement patterns from monocular video images -part i: A new framework for modeling human motion[END_REF][START_REF] Sminchisescu | Estimating articulated human motion with covariance scaled sampling[END_REF][START_REF] Ziheng | A bayesian framework for extracting human gait using strong prior knowledge[END_REF].

The human body being composed of 206 bones, it is easily understandable that those models are highly simplified representations of the human skeleton. Only articulations critical to gait description are considered and only these articulations' most important DOF are included in the models. Thus, subtle pathological motions could be inadvertently discarded. For example, most articulated models used for motion capture only attribute one degree of freedom to the knee, whereas many musculoskeletal problems of the knee are characterized by more complex 3D rotations.

This situation contrasts with statistical shapes models currently used to process radiological spine images. Those statistical models are based on point-to-point correspondences [START_REF] Lorenz | Generation of point-based 3d statistical shape models for anatomical objects[END_REF][START_REF] Benameur | 3D/2D registration and segmentation of scoliotic vertebrae using statistical models[END_REF] and many points are sampled on each vertebra, therefore the number of DOF could reach several hundreds. These statistical shape models describe efficiently the local shape of the vertebrae and they can improve the performances of segmentation and registration algorithms. However, they fail to take into account that the spine is a collection of vertebrae linked together by soft tissues (discs, muscles, ligaments, etc.). Capturing the articulated nature of the spine will lead to a statistical shape model that can be interpreted more intuitively and that will uncover more clinically relevant relationships.

Point-to-point correspondences are therefore not the best choice of primitives. A more natural choice would be to use frames (points associated with three orthogonal axes) to deal with articulated anatomical structures. The main reason for this choice is that frames enable a more natural analysis of the relative orientations and positions of the models, whereas point primitives would only properly deal with positions. In this context, articulated models of the spine will be created by associating a frame to each vertebra. The deformations of the spine shape will then be described in terms of rigid transformations applied to those frames.

A variability model of inter-vertebral rigid transformations based on Lie group properties was recently proposed [START_REF] Boisvert | 3D anatomical variability assessment of the scoliotic spine using statistics on Lie groups[END_REF]. However, the inter-vertebral transformations were analyzed individually thus discarding the coupling between different vertebral levels. It is difficult to take into account the coupling between all vertebral levels because of the high dimensionality of the model. In fact, there are 5 lumbar and 12 thoracic vertebrae (excluding cervical vertebrae) and a rigid transformation has 6 DOF, which means that an articulated description of the spine shape is characterized by at least 102 DOF. A clinician can hardly perform the analysis of such large variability model. It is therefore necessary to find a way to extract only the most meaningful modes of variability. These modes will indicate what varies the most from a geometric perspective in a given group, thus helping the physician in the task of analyzing large sets of 3D spine models, which is a necessary but tedious task to define a clinically relevant 3D classification.

In this context, the main contributions of this paper are to propose a method based on Riemannian geometry to perform principal components analysis on articulated spine models and to apply that method to a large database of scoliotic patients in order to construct the first statistical atlas of 3D deformation patterns for adolescent idiopathic scoliosis (a pathology that causes spine deformations). The combination of point-based and frame-based model will also be demonstrated to analyze the relation between the shape of the spine and the shape of individual vertebrae.

The remaining of the manuscript is divided in four sections. The next section review relevant works performed about scoliosis deformities analysis and about dimensionality reduction applied to articulated models. Section 3 then presents the proposed methods and the acquisition materials used for this study. The results obtained on large databases of scoliotic patients are presented in Section 4. Finally, Section 5 summarizes the conclusions and new perspectives that resulted from this study.

Background

3D Spine Models Classifications

Analyzing large databases of 3D spine models is becoming more and more important because there is a consensus among orthopedic surgeons that spinal deformations such as adolescent idiopathic scoliosis (AIS) are three-dimensional. However, there is no consensus yet regarding how to use 3D information to improve patients' care.

Classification is a decisive part of AIS patient assessment because it is used in the preoperative planning of corrective surgeries. King first proposed a classification of frontal thoracic curves, which was used to guide surgeries performed with the Harrington instrumentation [START_REF] King | Selection of fusion levels for posterior instrumentation and fusion in idiopathic scoliosis[END_REF]. This classification did not address sagittal deformations neither lumbar curves. Moreover, its reliability was contested [START_REF] Lenke | Intraobserver and interobserver reliability of the classification of thoracic adolescent idiopathic scoliosis[END_REF][START_REF] Ian | Identifying sources of variability in scoliosis classification using a rule-based automated algorithm[END_REF]. Lenke later proposed another classification [START_REF] Lenke | Adolescent idiopathic scoliosis: a new classification to determine extent of spinal arthrodesis[END_REF] that includes lumbar and sagittal modifiers. Lenke classification is more complete; however combining 2D measures is different from performing a real 3D analysis. Thus, new three-dimensional classifications of the spinal deformations associated with AIS are being investigated.

However, taking into account the three-dimensional nature of the deformations is far from being trivial. Many classifications were recently introduced. Duong et al. [START_REF] Duong | Three-dimensional classification of spinal deformities using fuzzy clustering[END_REF] used a wavelet transform of the vertebrae centroids and a clustering method to agglomerate similar 3D spine shapes. The method is technically elegant but its complexity might not be compatible with clinical practice. Poncet et al. [START_REF] Poncet | Geometric torsion in idiopathic scoliosis: three-dimensional analysis and proposal for a new classification[END_REF] proposed a classification that is based solely on spine torsion. Spine torsion is certainly an important factor to take into account but Duong et al. demonstrated that it is not the only one. Finally, Negrini et al. [START_REF] Negrini | Three-dimensional easy morphological (3-DEMO) classification of scoliosis, part I[END_REF][START_REF] Negrini | The three-dimensional easy morphological (3-DEMO) classification of scoliosis -part III, correlation with clinical classification and parameters[END_REF][START_REF] Negrini | The three-dimensional easy morphological (3-DEMO) classification of scoliosis, part II: repeatability[END_REF] recently proposed a classification based on "quasi-3D" indices (measure performed on multiple 2D projections). The main drawbacks of their approach are the ad hoc character of the indices selected and the two-dimensional nature of the classification.

In summary, analyzing large databases of 3D spine models is an important task but it is also a difficult and time-consuming one. Automatically extracting the most important deformation modes from a set of articulated spine models would greatly help surgeons to reliably and efficiently explore 3D spine models databases, which is a necessary prerequisite to the definition of a clinically relevant 3D classification of scoliotic deformities.

Principal components analysis of the 3D articulated spine models

Dimensionality reduction applied to the spine or to articulated models is not a new idea and methods were proposed in the past. As a part of a method that aim to predict the geometry of the spine based on the geometry of the trunk, Bergeron et al. [START_REF] Bergeron | Prediction of anterior scoliotic spinal curve from trunk surface using support vector regression[END_REF] performed a principal component analysis on the 3D coordinates of vertebrae centers in the frequency domain. The most important limitation of this approach is that vertebrae orientations are not modeled. Moreover, performing the analysis in the frequency domain does not allow a physically grounded intuitive interpretation of the results. Benameur et al. [START_REF] Benameur | 3D/2D registration and segmentation of scoliotic vertebrae using statistical models[END_REF] and Fleute et al. [START_REF] Fleute | Integrated approach for matching statistical shape models with intra-operative 2D and 3D data[END_REF] both proposed principal components analysis (PCA) based spine registration method but only the shapes of the vertebrae were considered not the global spine shape.

Principal components analysis was also used to process articulated body models. Yacoob and Black [START_REF] Yacoob | Parameterized modeling and recognition of activities[END_REF] applied PCA analysis on motion parameters measured with respect to the image plane in the context of activity recognition. Unfortunately, since the measurements are not intrinsic to the geometry of the articulated object the same viewpoint and camera have to be used, which is unrealistic for clinical studies that involves multiple research centers. Gonzalez et al. [START_REF] Gonzalez | Analysis of human walking based on aspaces[END_REF] and Green et al. [START_REF] Green | Quantifying and recognizing human movement patterns from monocular video images-part ii: applications to biometrics[END_REF] both performed a PCA analysis on aligned time-series of articulated models descriptions in order to recognize human actions. Jiang and Motai [START_REF] Jiang | Learning by observation of robotic tasks using on-line pca-based eigen behavior[END_REF] proposed a similar approach for on-line robot learning. Finally, Al-Zubi and Sommer [START_REF] Al-Zubi | Learning deformations of human arm movement to adapt to environmental constraints[END_REF] performed PCA to learn articulated motions of human arms. These approaches used representations that were either only based on 3D coordinates or only based on an angular description of the articulated body. However, using both positions and orientations would allow a better separation of different physiological phenomena such as pathological deformations and normal growth (that is expected to be dominated by translation). Moreover, these methods compute the mean without taking into account that 3D rotations form a curved manifold. Therefore, the obtained results depend on the arbitrary configuration of the local coordinates systems used to perform angular measurements.

Material and Methods

This section is divided in three subsections. One of the problems with previous attempts to find variability modes from 3D spine models was that their underlying spine geometry representations were based on pointto-point correspondences, which are not adapted to an analysis of orientations. The first subsection (3.1) will thus present an articulated description of the human spine that includes both vertebrae's relative position and orientation as well as the method used to estimate these articulated spine models from radiographs. The chosen articulated description uses rigid transformations to encode the relative position and orientation of the vertebrae. However, rigid transformations naturally belong to a Lie group and conventional statistical methods usually apply only in vector spaces. To alleviate this problem, the second subsection (3.2) will present how to rigorously generalize concepts such as the mean and the covariance to articulated spine models.

Finally, the last subsection (3.3) will introduce a simple mechanism to extract and visualize the most significant deformations modes observed in a database of articulated spine models.

Articulated Spine Models from Radiographs

Articulated models could be estimated from a wide variety of image modalities. Computed tomography or magnetic resonance could be used since the images recorded using those modalities are three-dimensional. However, large portions of the spine deformations are lost when patients are asked to lie down. Moreover, adolescent idiopathic scoliosis patients often present postural problems, which are relevant to the analysis. In Multi-planar radiography is a simple technique where two (or more) calibrated radiographs of a patient are taken to compute the 3D position of anatomical landmarks using a triangulation algorithm. It is one of the few imaging modalities that can be used to infer the three-dimensional anatomy of the spine when the patient is standing up. Furthermore, bi-planar radiography of scoliotic patients is routinely performed at Sainte-Justine hospital (Montreal, Canada).

The anatomical landmarks used in this study are illustrated in Figure 1. The identification of all anatomical landmarks is time consuming and costly because it is performed by a qualified technician. However, it provides a more robust anatomical matching of landmarks, which ensures a more precise 3D reconstruction. Thus, only the anatomical landmarks 1 to 6 were used when only the shape of the spine and not the local shapes of the vertebrae was analyzed. These anatomical landmarks were identified on each vertebra from T1 (first thoracic vertebra) to L5 (last lumbar vertebra) on a posterior-anterior and a lateral radiograph. The 3D coordinates of the landmarks were then computed and the deformation of a high-resolution template using dual kriging yielded a total of 24 reconstructed landmarks. The accuracy of this method was previously established to 2.6 mm [START_REF] Aubin | Morphometric evaluations of personalised 3d reconstructions and geometric models of the human spine[END_REF].

Once the anatomical landmarks presented in Figure 1 are reconstructed in 3D, they are expressed in the local coordinates system of their respective vertebra (see Stokes [START_REF] Stokes | Three-dimensional terminology of spinal deformity. a report presented to the scoliosis research society by the scoliosis research society working group on 3-D terminology of spinal deformity[END_REF] for a formal definition of the local coordinates system of a vertebrae). These local landmarks coordinates p 1 , p 2 , . . . , p n (where n is the number of anatomical landmarks used for each vertebra) will be used to analyze the vertebrae shapes. Each vertebra is then rigidly registered to its first upper neighbor and the resulting rigid transformations are recorded. By doing so, the spine is represented by a set of rigid transformations (see the Figure 2). These inter-vertebral transformations T 1 , T 2 , . . . , T 17 will be used to analyze the spine shape.

The combination of the local landmarks coordinates and the inter-vertebral transformations is especially well adapted to an analysis of the spinal deformities since the inter-vertebral rigid transformations describe the state of the physical links that alter the shape of the whole spine while the local landmarks coordinates can quantify the local vertebral deformations.

Most scoliotic patients are adolescents or pre-adolescents and the spine length of these patients varies considerably. At this point one might be tempted to normalize the spine models to improve comparability. However, the development of many musculoskeletal pathologies, such as adolescent idiopathic scoliosis, is tightly linked with the patient growth process. Thus, normalization could discard valuable information and should be avoided in most applications.

Elements of Articulated Models Statistics

The articulated spine models presented in the last subsection express the spine shape using a vector of intervertebral rigid transformations. This idea is very intuitive and models closely the pathological progression mechanisms of spinal deformities. However, rigid transformations do not naturally belong to a vector space (they belong to a curved Riemannian manifold) and simple operations necessary to conventional multi-variate analysis such as addition and scalar multiplication are not defined on those. A generalization of concepts as simple as the mean and the covariance is therefore needed.

One could be tempted to parameterize the rigid transformations using one of the many available representations (Euler angles and translation vector, for example) and then to compute the mean and covariance on that representation. However, the result would not be intrinsic. In other words, the mean and covariance would be dependent on the viewpoint from which the orientation and position would have been measured.

The distance is a general concept that can be used to perform those generalizations and Riemannian geometry offers a mathematical framework to work with primitives when only a distance is available. Concepts of probability and statistics applied to Riemannian manifolds were recently studied [START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF]. Moreover, a Riemannian framework was also used in the context of statistical shape modeling to perform PGA (principal geodesic analysis) on medial axis representations [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF].

Exponential and Logarithmic Maps

In a complete Riemannian manifold M the shortest smooth curve γ(t) such that γ(0) = x and γ(1) = y (with x, y ∈ M) is called a geodesic and the length of that curve is the distance between x and y. Two important maps can be defined from the geodesics: the exponential map Exp x which maps a vector ∂ x of the tangent plane T x M to the element reached in a unit time by the geodesic that starts at x with an initial tangent vector ∂ x and the logarithmic map Log x which is the inverse function of Exp x . In other words, these two maps enable us to "unfold" the manifold on the tangent plane (which is a vector space) and to project an element of the tangent plane to the manifold (see Figure 3 for an example of Log x and Exp x ).

In this paper, the spine shape is modeled as a set of frames associated to local coordinates systems of vertebrae. The modifications of the spine geometry are then modeled as rigid transformations applied to those frames. In order to compute the principal deformations modes (from equation 7), the exponential and logarithmic maps associated with a distance function on rigid transformations are thus needed.

A rigid transformation is the combination of a rotation and a translation. Defining a suitable distance on the translational part is not difficult since 3D translations belong to a real vector space. However, the choice of a distance function between rotations is more complex. Rotations are usually expressed using a three by three matrix that belongs to the special orthogonal group SO 3 and simple solutions, like the norm of the difference, does not respect the prerequisites conditions for distance functions.

To define a suitable distance between rigid transformations, another representation of the rotations called the rotation vector is needed. This representation is based on the fact that a 3D rotation can be fully described by an axis of rotation supported by a unit vector n and an angle of rotation θ (see figure 4). The rotation vector r is then defined as the product of n and θ.

The conversion from the rotation vector to the rotation matrix is performed using the Rodrigues rotation equation (see [START_REF] Goldstein | Classical Mechanics[END_REF] for details):

R = I + sin(θ)S(n) + (1 -cos(θ))S(n) 2 with S(n) =   0 -n z n y n z 0 -n x -n y n x 0  
And the inverse map (from a rotation matrix to a rotation vector) is given by the following equations:

θ = arccos( T r(R) -1 2 ) and S(n) = R -R T 2 sin(θ) (1)
Using the rotation vector representation, a left-invariant distance (d(T

3 • T 1 , T 3 • T 2 ) = d(T 1 , T 2
)) between two rigid transformations can easily be defined:

d(T 1 , T 2 ) = N ω (T -1 2 • T 1 ) with N ω (T ) 2 = N ω ({r, t}) 2 = ωr 2 + t 2 (2)
Where ω is used to weight the relative effect of rotation and translation, r is the rotation vector and t the translation vector. Because the selected distance function is left-invariant, we have Exp

µ (x) = µ • Exp Id (D(µ) -1 x) and Log µ (T ) = D(µ)Log Id (µ -1 • T ) with D(x) = ∂ ∂y x • y| y=Id .
Furthermore, it can be demonstrated that the exponential and log map associated with the distance of equation 2 are the mappings (up to a scale) between the combination of the translation vector and rotation vector and the combination of the rotation matrix and the translation vector [START_REF] Pennec | Computing the mean of geometric features -application to the mean rotation[END_REF].

Exp Id (T ) = R(r) t and Log Id (T ) = ωr(R) t
The Riemannian structure of the anatomical landmarks used to quantify the local shape of the vertebrae is much more simple. In fact, anatomical landmarks simply belong to 3 , which is a vector space. Thus, with respect to the Euclidian distance, the anatomical landmarks belong to a flat Riemannian manifold. The Exp and Log maps are then as simple as translations of the original landmarks: Exp p (q) = p + q and Log p (r) = r -p with p, q, r ∈ 3

Fréchet Mean

For a given distance, the generalization of the usual mean can be obtained by defining the mean as the element µ of a manifold M that minimizes the sum of the distances with a set of elements x 0...N of the same manifold M:

µ = arg min x∈M N i=0 d(x, x i ) 2
This generalization of the mean is called the Fréchet mean. Since it is defined using a minimization, it is difficult to compute it directly from the definition. However, it can be computed using a gradient descent performed on the summation. The following iterative equation summarizes this operation:

µ n+1 = Exp µn ( 1 N N i=0 Log µn (x i )) (3)

Generalized Covariance

The variance (as it is usually defined on real vector spaces) is the expectation of the L 2 norm of the difference between the mean and the measures. An intuitive generalization of the variance on Riemannian manifolds is thus given by the expectation of a squared distance:

σ 2 = 1 N N i=0 d(µ, x i ) 2 (4)
To create statistical shape models it is necessary to have a directional dispersion measure since the anatomical variability of the spine is anisotropic [START_REF] Boisvert | 3D anatomical variability assessment of the scoliotic spine using statistics on Lie groups[END_REF]. The covariance is usually defined as the expectation of the matricial product of the vectors from the mean to the elements on which the covariance is computed. Thus, a similar definition for Riemannian manifolds would be to compute the expectation in the tangent plane of the mean using the log map:

Σ = 1 N N i=0 Log µ (x)Log µ (x) T (5)

Multivariate Case

The Fréchet mean and the generalized covariance make it possible to study the centrality and dispersion of one primitive belonging to a Riemannian manifold. However, to build complete statistical shape models, it would be most desirable to study multiple primitives altogether. Therefore, a generalized cross-covariance Σ f g is needed.

Σ f g = 1 N N i=0 Log µ f (f i )Log µg (g i ) T
A natural extension is to create a multivariate vector f = [f 1 , f 2 , f 3 , . . . , f k ] T where each element corresponds to a part of a model made of several primitives. In our case, these primitives are the inter-vertebral rigid transformations and the local anatomical landmarks described in Figure 1. The mean and the covariance of this multivariate vector will thus be:

µ =      µ 1 µ 2 . . . µ k      and Σ =      Σ f 1 f 1 Σ f 1 f 2 . . . Σ f 1 f k Σ f 2 f 1 Σ f 2 f 2 . . . Σ f 2 f k . . . . . . . . . Σ f k f 1 Σ f k f 2 . . . Σ f k f k      (6)
This is very similar to the conventional multivariate mean and covariance except that the Fréchet mean and the generalized cross-covariance are used in the computations.

Principal Deformations

The equation 6 allows us to compute a statistical shape model for a group of models made of several primitives. However, the different primitives will most likely be correlated which makes the variability analysis very difficult. Furthermore, the dimensionality of the model is also a concern and we would like to select only a few important uncorrelated components.

Unlike the manifold itself, the tangent plane is a vector space and its basis could be changed using a simple linear transformation. Thus, we seek an orthonormal matrix A (AA T = I) to linearly transform the tangent plane ( Log µ (g) = ALog µ (f ) ) such as the generalized covariance in the transformed tangent space is a diagonal matrix (Σ gg = diag (λ 1 , λ 2 , . . . , λ k )). The covariance matrices of the transformed tangent space and of the original tangent space are connected by the following equation:

Σ gg = diag (λ 1 , λ 2 , . . . , λ k ) = AΣ f f A T If A is rewritten as A = [a 1 ,
a 2 , . . . , a k ] T , then it is easy to show that:

[λ 1 a 1 , λ 2 a 2 , . . . , λ k a k ] = [Σ f f a 1 , Σ f f a 2 , . . . , Σ f f a k ] (7)
The line vectors of the matrix A are therefore the eigenvectors of the original covariance matrix and the elements of the covariance matrix in the transformed space are the eigenvalues of the original covariance. This is the exact same procedure that is used to perform PCA in real vector spaces. Like for real vector spaces, the variance is left unchanged since σ 2 = Tr(Σ f f ) = Tr(Σ gg ) and the cumulative fraction of the variance explained by the first n components is:

p = 1 σ 2 i=1...n λ i
A shape model can be re-created from coordinates of the transformed tangent space simply by going back to the original tangent space and projecting the model on the manifold using the exponential map. So if α i is the coordinate associated with the i th principal component, the following equation can be used to re-create a shape model:

S = Exp µ ( k i=1 α i a i )

Results and Discussion

This section presents and discusses the results obtained using the method previously described. Principal spine deformations modes are presented in three different contexts. First of all, the principal spine deformation modes of a group of patients afflicted by AIS is presented. These results indicated that the principal deformation modes can efficiently extract the most important modes of variability and that those modes have a strong clinical meaning. Then, a comparison of the principal deformations modes of two groups of scoliotic patients classified according to the Lenke classification [START_REF] Lenke | Adolescent idiopathic scoliosis: a new classification to determine extent of spinal arthrodesis[END_REF] is introduced. This experiment demonstrates that the principal deformation modes can extract clinically relevant patterns that could easily be missed by a direct visualization of the spine models. Finally, we present results obtained from a group of spine models with 14 landmarks by vertebra where we show that global deformations of the spine are correlated with vertebral deformations.

Spine Shape Analysis

The method described in the previous sections was applied to a group of 307 scoliotic patients. The patients selected for this study had not been treated with any kind of orthopedic treatment when the radiographs were taken. Therefore, the inter-patients variability observed was mainly caused by anatomical differences and not by any treatment. A total of 6 anatomical landmarks were identified manually by a qualified technician, this number of landmarks is sufficient to analyze the shape of the spine, but is not enough to properly analyze vertebrae shapes. Thus, only the inter-vertebral rigid transformations were used for this experiment.

The only additional parameter that needs to be selected is ω, which balances the relative importance of the translation and rotation (see Equation 2). In some applications, rotations are more important than translations (or the converse) and this parameter may be used to bias the results toward the most important part of the rigid transformations. In our case, however, we did not want to introduce such bias. Thus, ω was selected so that the rotations and the translations had the same contribution to the global variance. This criteria lead to a numerical value of ω = 20.

To illustrate the different deformation modes retrieved using the proposed method, four models were reconstructed for each of the first five principal deformation modes. Those models were reconstructed by setting

α k to -3 √ λ k , - √ λ k , √
λ k and 3 √ λ k for k = 1 . . . 5 while all others components (α i with i = k) were set to zero (see figures 5 and 6). In theory, a total of 102 principal deformation modes could be extracted. However, 102 figures would certainly overload this paper. The first five modes (in term of variance explained) illustrate well the power of the proposed method since they can be understood without an extensive medical background. The interpretation of less important deformations modes is subtler. Furthermore, the 3D spine models in our databases are corrupted by a certain amount of acquisition noise. Thus, principal deformations modes associated with small variances are more likely to be dominated by noise.

A visual inspection revealed that the first five principal deformation modes have clinical meanings. The first mode is associated with the patient growth because it is mainly characterized by an elongation of the spine and also includes a mild thoracic curve. The second principal deformation mode could be described as a double thoraco-lumbar curve, because there are two curves: one in the thoracic segment (upper spine) and another in the lumbar segment (lower spine). The third principal mode of deformation is a simple thoracic curve (the apex of the curve is in the thoracic spine), but it is longer than the thoracic curve observed in the first principal component. It is also interesting to note that, in addition to the curves visible on the posterior-anterior view, the second and third principal deformation modes are also associated with the development of a kyphosis (back hump) on the lateral view. The fourth component is a lumbar lordosis (lateral curve of the lumbar spine). Finally, the fifth component is related to the patient's frontal balance.

Those curve patterns are routinely used in different clinical classifications of scoliosis (used to plan surgeries). For instance, the reconstructions built from the first principal deformation mode would be classified using King's classification [START_REF] King | Selection of fusion levels for posterior instrumentation and fusion in idiopathic scoliosis[END_REF] as a type II or III (depending on which reconstruction is evaluated), the second deformation mode would be associated to King's type I or III and the third principal deformation could be associated to King's type IV.

Previously those patterns were derived from surgeons' intuition using 2D images and clinical indices, whereas it is now possible to automatically compute those patterns from statistics based only on 3D geometries.

Furthermore, the cumulative variance explained by an increasing number of principal deformations modes (illustrated in Figure 7) shows the capacity of the proposed method to reduce the number of dimensions of the model while keeping a large part of the original variance. This property could be used in registration or reconstruction algorithms. Also, since about 50 % of the variance is explained by the first five principal deformations modes, it appears that a careful analysis of more principal deformation modes by trained medical experts might reveal more clinically pertinent knowledge in the future.

The variability observed from our dataset is predominantly associated with anatomical variability, but it also includes variability caused by other factors such as patients' postures and landmark reconstruction errors. Therefore, the principal deformations modes associated with small eigenvalues are likely to be more closely related to acquisition noise than to anatomical variations between individuals.

Surgical Classifications Analysis

Classifications of spinal deformities are used for several purposes : as a training tool for new residents, as mean of communication between practitioners, as a guideline for the selection of a corrective treatment, etc. Current classifications used for adolescent idiopathic scoliosis are based on 2D measures performed on radiographs. However, the deformities are three-dimensional and therefore it is very likely that clinically relevant sub-groups could be identified by a 3D analysis.

Two groups of untreated AIS patients with two different curve types (according to the Lenke classification [START_REF] Lenke | Adolescent idiopathic scoliosis: a new classification to determine extent of spinal arthrodesis[END_REF]) were identified from the Sainte-Justine Hospital's archive: 86 patients with type I curves and 47 patients with type V curves. Type I curves are characterized by a main thoracic curve while type V curves are thoracolumbar or lumbar curves. A qualified physician performed the classification and six anatomical landmarks were identified on each radiograph. Only the inter-vertebral rigid transformations were used to compute the principal deformation modes (not the anatomical landmarks). The first three principal modes of deformations (with ω = 20) computed from the two groups of patients are presented in Figure 8 and Figure 9.

Figure 8 shows that the most important mode of deformation for scoliotic spines of type I is an elongation of the spine combined with the development of a right thoracic curve and little curvature modification visible from the lateral view. Type I curves are defined as main thoracic curvatures and right curves are more prevalent then the left among scoliotic patients. Moreover, smaller patients are usually younger and scoliosis progress during growth. Thus, this deformation mode is compatible with the progression of the disease for patient that suffers from a curve of the first type. The second mode of deformation is associated with a double curve and strong sagittal curvature of the lower spine (lordosis) caused largely by a rotation of L5 (the lowest lumbar vertebra). This L5 rotation could indicate a deformation of a structure below the spine, thus it would be very interesting to also analyze the position and orientation of the pelvis. Finally, the third mode exhibit the development of a right thoracic curve, but a curve with a lower apex (close to T9-T10) than the right thoracic curve present in the first deformation mode.

The deformations modes presented in Figure 9 were computed in the database of scoliotic models with curves of the fifth type. The first mode of deformation in this case is characterized by an elongation of the spine associated with the creation of a left thoracolumbar curve. Lenke defined curves of the fifth type as thoracolumbar or lumbar curves. This first deformation mode is therefore similar to the first deformation mode of type 1 curves as they both depict the worsening of the 2D curve they were aimed to represent and are combined with an elongation of the spine. The second mode of deformation is thoracolumbar curve going from left to right with an increasing kyphosis. Finally, the third mode is the transformation of a double curve with very little lordosis and kyphosis to a small thoracolumbar curve with large lordosis and kyphosis.

These two figures indicate that the principal deformations modes can efficiently capture the clinically in- tended deformation mode in a given category of a surgical classification system (main thoracic for curves of the first type and thoracolumbar/lumbar for curves of the fifth type). It is also interesting to note that other important 3D curve patterns could also be identified. Those minor patterns would most probably be missed by a direct examination of the 3D spine models, which demonstrates the usefulness of the method to analyze rapidly and efficiently large databases of articulated models.

Combined Spine Shape and Local Vertebrae Shape Analysis

Surgical classifications of the spine are more concerned with the shape of the whole spine and less with the shape of individual vertebrae. This small interest is due to the fact that little can be done to correct vertebral deformations, especially if the patient's growth is over. However, there are several applications where a combined analysis of both spine and vertebrae shapes is very important. First of all, image-guided spine surgeries are based on efficient and reliable registration methods. Those methods optimize a similarity measure between a spine model and spine images by deforming the spine model. Thus, registration methods would benefit from the correlation that exists between spine and vertebrae deformations. Moreover, vertebral deformations and spinal deformations in AIS patients are probably physiologically linked [START_REF] Stokes | Mechanical modulation of vertebral body growth. implications for scoliosisprogression[END_REF]. Thus analyzing vertebral deformations in the context of their underlying spinal deformation will lead to a better understanding of scoliosis progression.

Vertebrae shapes were already studied on dry human specimens [START_REF] Parent | Morphometric analysis of anatomic scoliotic specimens[END_REF], their deformation patterns were documented on scoliotic chickens [START_REF] Coillard | Vertebral deformities and scoliosis[END_REF] and asymmetric vertebral growth was experimented on rats [START_REF] Stokes | Mechanical modulation of vertebral body growth. implications for scoliosisprogression[END_REF]. However, our method offers the possibility of analyzing real patients and to visualize the relationships between the spinal deformities and vertebral deformities.

For this experiment 117 patients were selected and 14 anatomical landmarks for each vertebrae were digitized from the patients' radiographs. The principal deformation modes were then extracted using both the intervertebral transformations and the anatomical landmarks. Figure 10 illustrates that large deformations of the spine can be connected with vertebral deformations. In this particular case, the development of a thoracic curve (the camera does not move) is associated with an asymmetrical deformation of the vertebral bodies near the curve apex.

Although, vertebral deformations appears to be connected with spinal deformations our method cannot conclude if there is a causality link (or what this link is). The relationships observed between spinal curves types and vertebral deformations could however be used to validate conclusions obtained from animal studies or from biomechanical simulations.

Conclusion

A method to extract the principal modes of deformation from articulated models was described. The method consists in performing a principal component analysis in the tangent space of a Riemannian manifold. This manifold is composed of rigid transformations that describe the spine shape and also includes anatomical landmarks. We applied this method to databases of scoliotic patients reconstructed in 3D using biplanar radiographs.

Results suggest that PCA applied to a suitable representation of the spine, namely a set of rigid transformations, leads to an algorithm that can expose natural modes of deformation of the spine. Clinically relevant patterns of deformations were extracted and dimensionality reduction was successfully achieved. Patterns defined in existing surgical classifications were automatically identified from an unclassified group of scoliotic patients. The relationships that exists between clinically used 2D classifications and principal deformations modes were also assess qualitatively. Quantitative analysis will be performed in the future to demonstrate the qualities of the proposed method to the medical community.

Moreover, the application of the proposed method in the context of 3D classification of the scoliotic deformations was demonstrated by comparing the principal deformation modes of two groups of the Lenke classification. This experiment lead to the identification of deformations patterns that would have been very difficult to identify using another approach. An analysis of different patterns extracted from all Lenke's classes using the proposed method would thus be interesting to realize in the future.

One important characteristic of the proposed method behind these promising results is the absence an arbitrary selection of clinical indices. All the information needed to describe the 3D spine models is used, but it is regrouped into uncorrelated modes and sorted in descending order of explained variance. The proposed method is therefore a promising analytical tool and it will certainly contribute to the definition of new 3D classifications.

We also showed that the vertebral deformations could be linked to the deformation of the whole spine. This experiment is compatible with the "vicious cycle" described by Stokes et al. [START_REF] Stokes | Mechanical modulation of vertebral body growth. implications for scoliosisprogression[END_REF] that could explain a progression mechanism of scoliosis. However, describing vertebrae shapes with anatomical landmarks might be simplistic and alternate representation such as spherical harmonics or medial axis representations could also be used.

A very promising perspective is the longitudinal analysis of patients. Multiple acquisition of the patients' spine geometry could be used to extract progression modes instead of deformation modes. This could be a significant breakthrough since patients belonging to more severe progression modes could be followed more closely by their physician and patients with lower progression risks could be treated less aggressively.

Finally, one of the reasons to perform dimensionality reduction on statistical shape models is to reduce the number of DOF that needs to be optimized during model registration. The proposed method will therefore be integrated to a spine registration algorithm in the future. It might also be useful for the integration of a large number of rigid structures in non-rigid registration procedures [START_REF] Little | Deformations incorporating rigid structures[END_REF] of the whole human torso.

Figure 1 :

 1 Figure 1: Anatomical landmarks used to digitize the spine geometry and to analyze local vertebra shape. Top: lateral view. Bottom: posterior view.

Figure 2 :

 2 Figure 2: Frames and transformations used to express the spine as an articulated model.

Figure 3 :

 3 Figure 3: Graphical representation of the Log X and Exp X maps on a three-dimensional sphere.

Figure 4 :

 4 Figure 4: Rigid transformation expressed by an axis of rotation n, an angle of rotation θ and a translation vector t

Figure 5 :

 5 Figure 5: First, second and third principal deformation modes of an AIS patient database. Posterior-anterior view (a) and lateral view (b) of the first principal deformation mode. Posterior-anterior view (c) and lateral view (d) of the second principal deformation mode. Posterior-anterior view (e) and lateral view (f) of the third principal deformation mode. Spine models were rendered for -3, -1, 1 and 3 times the standard deviation explained by the corresponding deformation mode.

Figure 6 :

 6 Figure 6: Fourth and fifth principal deformation modes of an AIS patient database. Posterior-anterior view (a) and lateral view (b) of the fourth principal deformation mode. Posterior-anterior view (c) and lateral view (d) of the fifth principal deformation mode. Spine models were rendered for -3, -1, 1 and 3 times the standard deviation explained by the corresponding deformation mode.

Figure 7 :

 7 Figure 7: Fraction of the variance explained by the principal deformation modes

Figure 8 :

 8 Figure 8: First, second and third principal deformation modes of AIS patient with type I curves. Posterioranterior view (a) and lateral view (b) of the first principal deformation mode. Posterior-anterior view (c) and lateral view (d) of the second principal deformation mode. Posterior-anterior view (e) and lateral view (f) of the third principal deformation mode. Spine models were rendered for -3, -1, 1 and 3 times the standard deviation explained by the corresponding deformation mode.

Figure 9 :

 9 Figure 9: First, second and third principal deformation modes of AIS patient with type V curves. Posterioranterior view (a) and lateral view (b) of the first principal deformation mode. Posterior-anterior view (c) and lateral view (d) of the second principal deformation mode. Posterior-anterior view (e) and lateral view (f) of the third principal deformation mode. Spine models were rendered for -3, -1, 1 and 3 times the standard deviation explained by the corresponding deformation mode.

Figure 10 :

 10 Figure 10: Asymmetrical deformation of a vertebral body located near the apex of a spinal curve (a phenomena also called vertebra wedging) observe in the third principal deformation modes of a group of scoliotic patients which had 14 anatomical landmarks reconstructed on each of their vertebrae. Model reconstructions (with a fixed camera) for -3 √ λ 3 (a), -√ λ 3 (b), √ λ 3 (c) and 3 √ λ 3 (d)

J. Boisvert et al. / Electronic Letters on Computer Vision and Image Analysis 7(4):[START_REF] Green | Quantifying and recognizing human movement patterns from monocular video images-part ii: applications to biometrics[END_REF][START_REF] Green | Quantifying and recognizing human movement patterns from monocular video images -part i: A new framework for modeling human motion[END_REF][START_REF] Jiang | Learning by observation of robotic tasks using on-line pca-based eigen behavior[END_REF][START_REF] King | Selection of fusion levels for posterior instrumentation and fusion in idiopathic scoliosis[END_REF][START_REF] Lenke | Intraobserver and interobserver reliability of the classification of thoracic adolescent idiopathic scoliosis[END_REF][START_REF] Lenke | Adolescent idiopathic scoliosis: a new classification to determine extent of spinal arthrodesis[END_REF][START_REF] Little | Deformations incorporating rigid structures[END_REF][START_REF] Lorenz | Generation of point-based 3d statistical shape models for anatomical objects[END_REF][START_REF] Negrini | The three-dimensional easy morphological (3-DEMO) classification of scoliosis, part II: repeatability[END_REF][START_REF] Negrini | The three-dimensional easy morphological (3-DEMO) classification of scoliosis -part III, correlation with clinical classification and parameters[END_REF][START_REF] Negrini | Three-dimensional easy morphological (3-DEMO) classification of scoliosis, part I[END_REF][START_REF] Parent | Morphometric analysis of anatomic scoliotic specimens[END_REF][START_REF] Pennec | Computing the mean of geometric features -application to the mean rotation[END_REF][START_REF] Pennec | Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements[END_REF][START_REF] Poncet | Geometric torsion in idiopathic scoliosis: three-dimensional analysis and proposal for a new classification[END_REF][START_REF] Sminchisescu | Estimating articulated human motion with covariance scaled sampling[END_REF][START_REF] Stokes | Three-dimensional terminology of spinal deformity. a report presented to the scoliosis research society by the scoliosis research society working group on 3-D terminology of spinal deformity[END_REF][START_REF] Stokes | Mechanical modulation of vertebral body growth. implications for scoliosisprogression[END_REF][START_REF] Ian | Identifying sources of variability in scoliosis classification using a rule-based automated algorithm[END_REF] 2008

Acknowledgement

This work was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Québec's Technology and Nature Research Funds (Fonds de Recherche sur la Nature et les Technologies de Québec) and the Canadian Institutes of Health Research (CIHR).


1099

In recent years there has been a significant development in medical robotics from basic research to product development, evaluation, and feasibility studies. While earlier research in medical robotics was targeted towards the development of computer-integrated systems for surgical interventions, over the years it has branched out to cover a wider variety of areas including image-guided therapy, rehabilitation robotics, and cellular-scale manipulation.

The goal of this special issue was to bring together papers in the wide area of medical robotics and share them with the robotics community at large. While clinical work is highly relevant to this area of research, this special issue has focused primarily on papers with a strong methodological component and clinical evaluation was encouraged primarily as a validation tool rather than a major component of the paper. The special issue encouraged submission in a wide range of areas including: robot-assisted procedures, smart instrumented tools for surgery, haptic feedback in medical robotics, rehabilitation robotics, interventional therapy, image-guided procedures, medical imaging for robotic interventions, cell manipulation, cellular-scale interventions, surgical simulation, soft-tissue modeling, and telesurgery.

The Call for Papers attracted a record number of submissions, 51 in total, encompassing most of the areas mentioned above. After a thorough review process including major and minor revisions of several manuscripts, 19 papers were accepted for publication and are published in two parts. These two issues cover several topic areas as discussed below.

Jagadeesan et al. discuss algorithms for robot-assisted catheter insertions with the eventual goal of minimizing radiation exposure for the patient and potentially improve the outcome of the procedure while Trejos et al. Tavakoli and Howe present their work on the effect of link and joint flexibility of the slave robot on the stiffness of the slave robot in a teleoperated system and its effect on the transparency of teleoperation while Yuen et al. discuss robotic motion compensation for beating heart intracardiac surgery using a one degree-of-freedom motion compensation system. Krupa et al. discuss their work on real-time motion stabilization with B-mode ultrasound using image speckle information and visual servoing. Hu et al. present their work on the design, development, and evaluation of a camera system that is capable of pan, tilt, and zoom with lighting for visualization of organs in minimally invasive surgery.

In the area of rehabilitation robotics, Saglia et al. discuss the development of a high-performance redundantly actuated parallel mechanism for ankle rehabilitation while Ishii et al. present their work on the development and evaluation of an oral rehabilitation robot that provides maxillofacial massage to patients with oral disorders such as temporomandibular joint (TMJ) disorders and dry mouth. Wolf et al. discuss their research on the development of a robotic neural interface for autonomous positioning of extracellular recording electrodes, which can be used to develop neuroprosthetic devices as well as carry out neuroscientific research.

In the area of soft-tissue modeling, Ciarletta et al. present their work on the development of the hyperelastic model for anisotropic fiber reinforcements within intestinal walls while Gosline and Hayward present their work on the haptic synthesis of viscoelastic tissue properties using programmable eddy current brakes. Negeotte et al. discuss feasible paths for the needle through the tissue for suturing task to minimize tissue deformation using kinematic analysis and geometric modeling of the stitching task.

Finally, we would like to thank all of the reviewers who helped to review the manuscripts and the revised versions. We would also like to thank Professor John Hollerbach, Editor of The International Journal of Robotics Research, and Jennet Batten, Managing Editor of The International Journal of Robotics Research, for their support throughout the review process to bring this special issue to the robotics research community.

  discuss the development of a robot-assisted tactile sensing system for tumor localization using a commercially available sensor. Simaan et al. discuss the design and development of a teleoperated system for minimally invasive surgery of the throat while Piccin et al. present their work on the development of a teleoperated needle insertion device with force-feedback capability for percutaneous procedures. Martel et al. present their work in the area of nanorobotics under magnetic resonance imaging (MRI) guidance for interventions in target regions where catheterization is challenging. Lum et al. present their work on the development of RAVEN, a telesurgical system, while Berkelman et al. discuss the development of a compact and modular teleoperated robotic system for laparoscopic surgery. Venture et al. present their work on in-vivo estimation of joint dynamics of the human limbs with applications to diagnosis of neuromuscular disease while Westpahl et al. present their research on robot-assisted long bone fracture reduction.

at INRIA Sophia Documentation on July 26, 2013 ijr.sagepub.com Downloaded from

Guest Editors: Jaydev P. Desai University of Maryland, College Park, MD, USA Nicholas Ayache INRIA, Sophia-Antipolis Méditerranée, France


Ender Konukoglu

Olivier Clatz

Hervé Delingette

Nicholas Ayache

Personalization of Reaction-Diffusion Tumor Growth Models in MR images: Application to Brain Gliomas Characterization and Radiotherapy Planning

Reaction-diffusion based tumor growth models have been widely used in the literature for modeling the growth of brain gliomas. Lately, recent models have started integrating medical images, specifically anatomical and diffusion images, in their formulation. On the other hand, the adaptation of the general model to the specific patient cases has not been studied thoroughly yet. In this chapter we address this adaptation. This chapter is a short summary of the articles (Konukoglu 2009a), (Konukoglu 2009b) and the thesis (Konukoglu 2009c) that we have submitted recently. In the first part, we describe a parameter estimation method for reaction-diffusion tumor growth models using time series of medical (Magnetic Resonance) images. This method estimates the patient specific parameters of the model using the images of the patient taken at different successive time instances. In the second part of the chapter we focus on an application of the personalized models aimed to improve the tumor targeting in radiation therapy. Specifically we address the problem of limited visualization of medical images. We describe a method for extrapolating the invisible infiltration margins of gliomas in the MR images and the usage of these margins in constructing irradiation margins taking into account the growth dynamics of the tumor. Finally for both parts we show preliminary results demonstrating the power and the potential benefits of the personalization.

INTRODUCTION

Gliomas are neoplasms of glial cells which support and nourish the brain. These tumors have varying histopathological features and biological behavior showing different aggressiveness levels, from benign -grade I-to malignant -grade IV (Glioblastoma Multiforme). There has been a vast amount of research in mathematical modeling to describe the growth dynamics of these tumors (Byrne 2006, Cristini 2003, Frieboes 2007, Patel 2001, Stamatakos 2006, Zhang 2007). Lately, a specific type of macroscopic models, the reaction-diffusion models, received considerable attention from the literature in the attempt to link glioma growth models to medical images (Tracqui 1995, integrate information coming from medical images, specifically through anatomical and diffusion images, in their formulation. This integration is crucial for the transfer of mathematical models to the clinical applications since medical images are conventionally used for diagnosis and patient follow-up in the clinical routine. One of the biggest challenges in this transfer is the automatic adaptation of mathematical models to the patient based on images. In this chapter, first we address the problem of adapting the recent reaction-diffusion models to specific patient cases using time series of medical (Magnetic Resonance -MR) images (Konukoglu 2009a). Following this, we address the question of retrieving relevant information for radiotherapy planning from the personalized reaction-diffusion models (Konukoglu 2009c).

Reaction-Diffusion Models

Reaction-diffusion models describe the evolution of gliomas via proliferation of tumor cells and infiltration of the surrounding healthy tissue. The building block of these models is the reactiondiffusion type partial differential equations of this form:

    Ω n D D(x)            u u u u t u , 1  (1)
where u represents the tumor cell density, D is the local diffusion tensor (i.e. symmetric positive definite 3x3 matrix), ρ is the proliferation rate, Ω is the brain domain and ∂Ω represents the boundaries of the brain. The two terms on the right hand side correspond to the two phenomena described by the model: the diffusion term

  u    D
models the migration of tumor cells within the brain tissue and the reaction term ρu(1-u) models the proliferation of tumor cells.

Different models proposed in the literature mostly differ by the construction of the D tensor and the form of the proliferation term (Tracqui 1995, Swanson 2000, Giese 1996, Clatz 2005, Jbabdi 2005, Hogea 2007). In this chapter we concentrate on the personalization of the growth models proposed in (Clatz 2005) and (Jbabdi 2005). However, the methods explained in this chapter are independent from the exact construction of the tensor and can be applied to both of these formulations as well as to other reaction-diffusion models. The two aforementioned models are based on Equation (1) and they construct the diffusion tensor D as follows:

        matter white , , matter gray , x D x I D(x) water w g d f d (2)
The choice of the function f is defined differently in the two models:

        2005) (Jbabdi , , , ) ( diag ) 2005 Clatz ( , 1 T water V(x) V(x) D g g w w d d d x e d f
where d w is the diffusion rate in the white matter (Clatz 2005) or along the white matter fiber tracts (Jbabdi 2005), d g is the diffusion rate in the gray matter, D water =VΛV T is the diffusion tensor of water molecules at point x with e 1 being its highest eigenvalue.

Image-Guided Model Personalization

Once the mathematical model describing the general dynamics of the tumor growth process is created, the personalization is defined as the estimation of the model parameters based on the observations. The image-guided personalization focuses on determining the reaction-diffusion parameters, the diffusion tensor D (d w and d g ) and the proliferation rate ρ, based on the observed evolution of the tumor in the time series medical images. The difficulty in this estimation is due to the sparsity of the available information. The reaction-diffusion models describe the temporal evolution of tumor cell densities while in the images we only observe the evolution of the tumor delineation which is assumed to correspond to an iso-density contour (Burger 1988) as shown in Figure 1. Therefore, reaction-diffusion models are not directly applicable for the estimation problem.

[FIGURE 1 TO BE PLACED HEREimages, visibility problem] The problem of parameter estimation in the context of tumor growth models is a rather unexplored problem. A first attempt was made by Tracqui et al. in (Tracqui 1995) where they optimized the parameters of their model by comparing the area of the tumor observed in CT images at different times and the area of the simulated tumor. The drawback of this approach was to use tumor cell densities requiring an initialization of the density distribution throughout the brain while these densities are not observable in the images ( again, we observe tumor delineations). More recently, in (Hogea 2008), Hogea et al. have optimized their parameters using two different methods: by comparing locations of some manually placed landmarks with the model generated ones and by comparing tumor cell densities extracted from the images and generated by the model.

In addition to the parameters of the reaction-diffusion model they optimize the parameters of their mechanical model as well. The use of tumor cell densities in the optimization process has the same problems as the previously mentioned method, they are not observable. The landmarks on the other hand, seems to be a promising approach, however, the parameter values obtained in this approach depend on the assumed coupling between the mechanical and the patho-physiological model. The resulting parameters are not purely inherent to the tumor. Recently Swanson et al. in (Swanson 2008) proposed a parameter estimation method for the diffusion process in petri-dish experiments, which is consistent with the observed information in the images as it uses the tumor boundaries rather than tumor cell densities. They have derived analytical approximations for the evolution of the tumor delineation for 2 dimensional circular growth. Using the formulation for the tumor delineation they have estimated the diffusion coefficient for the petri-dish experiments. The difficulty one would encounter if one wants to apply this method to medical images is that the method assumes radial symmetric growth which is not the case in the brain (in-vivo). Moreover, the existence of a reaction term results in a different evolution than pure diffusion.

In this chapter, we explain and analyze a parameter estimation method for reaction-diffusion based tumor growth models using time series of MR images. The method is based on the evolution of the tumor delineation rather than tumor cell densities and in this respect it is consistent with the observations in the images. It also takes into account tissue inhomogeneities, fiber structures and the real geometry of both the patient's brain and the tumor while remaining consistent with the image information. We also present preliminary promising results on 2 sets of patient data.

Infiltration Extent of Gliomas and Constructing Model based Irradiation Margins

For the diagnosis and the therapy of gliomas, clinicians rely on medical images, such as Magnetic Resonance (MR) and Computed Tomography (CT) images, which show the mass part of the tumor. Current imaging techniques are not able to expose the low density infiltration (Tovi 1994, Tracqui 1995, Swanson 2004) posing a problem for the experts in outlining the whole tumor and in delineating its extent ( see Figure 1 ). In radiotherapy, this problem of visualizing low density infiltration is addressed by outlining a constant margin of 2 cm around the visible tumor boundary for irradiation and assuming the whole tumor infiltration is contained within this region (Seither 1995, Kantor 2001). This approach does not take into account the infiltration dynamics of gliomas, particularly the higher motility of tumor cells in white matter compared to gray matter (Giese 1996).

As a result, the irradiation region ignoring these dynamics might not reach the full extent of the tumor infiltration in white matter and irradiate healthy gray matter. The importance of this problem was shown by Swanson et al. in (Swanson 2002) where they compare the visible part of the tumor with the extent of the invisible infiltration for virtual tumors grown by reaction-diffusion models.

Personalized tumor growth models can offer solutions to this problem by integrating clinical information and theoretical knowledge about tumor cell dynamics. In this chapter we describe a new formulation which aims to solve the problem of estimating tumor cell density distribution beyond the visible boundary of gliomas in the images. Starting from the tumor delineation either found by a segmentation algorithm or manually drawn by an expert it produces a map of possible tumor infiltration. It uses the anatomical MR images and diffusion tensor images (DTI) to suggest irradiation margins taking into account the growth dynamics. We then use the extrapolated infiltration extents to create variable irradiation margins. The potential benefits of such margins in targeting tumor cells are also described.

METHODS

Parameter Estimation and Model Personalization

The reaction-diffusion model given in Equations ( 1) and (2) describes the temporal evolution of local tumor cell densities. In order to solve the parameter estimation problem we need a formulation consistent with the image: it should model the evolution of the tumor delineation rather than the evolution of the cell densities. In this section, first we present the traveling time formulation for tumor delineation which captures the same growth dynamics as the reactiondiffusion models while remaining consistent with the image information. Then we formulate the parameter estimation problem as an optimization problem based on the MR images using the traveling time formulation.

Traveling Time Formulation for Tumor Delineation

The asymptotic properties of the reaction-diffusion equations under certain conditions allow us to construct a formulation for the tumor delineation as described above, the traveling time formulation. Reaction-diffusion equations and their asymptotic properties have been well studied in the literature (Aronson 1978, Ebert 2000). The property on which we base our derivation is the existence of traveling wave solutions of the reaction-diffusion equations and the convergence of different initial conditions to these solutions in time.

The constant coefficient case of Equation ( 1) admits a traveling wave solution in the infinite cylinder given as

) ( ) ( ) , (  u vt u t x u     x n ,
where v is the speed of the wave (front of the u distribution), n is the direction of motion of the wave and

vt    x n 
is the moving frame of the traveling wave. There are two important characteristics of the reaction-diffusion equations and the traveling wave solutions that are very useful for our derivation:

• For the traveling wave solution, all iso-density contours of the distribution u move with the same speed v.

• Any compact support initial condition u(x,0) of the reaction-diffusion equation ( 1) converges to the traveling wave solution in time ( see Figure 2 ). Therefore, the traveling wave solution serves as a good approximation for a certain class of solutions of the reaction-diffusion equation.

Based on these two characteristics we can model the evolution of the tumor delineation through the speed of the traveling wave.

[PLACE FIGURE 2convergence of the reaction-diffusion initializations to the traveling wave]

The asymptotic speed v of the traveling wave solution is given as a function of the model parameters

Dn n'  2  v , (3)
where this speed is defined in the infinite cylinder under the constant coefficient assumption. In the case of the tumor modeling the shape of the tumor is arbitrary and the model coefficients can be spatially and temporally varying. Therefore, the above assumptions should be relaxed. Under the assumptions that the tumor delineation is planar and coefficients are constant within each image voxel, we can write a preliminary traveling time formulation for the tumor delineation using v:

1 ' 2    T T D  , (4)     x x 0 ) ( T T , (5)
where T(x) is an implicit time function representing the time when the tumor delineation passes through the point x. Equation 5 represents the Dirichlet type boundary condition stating that the tumor delineation reaches the surface Γ in T 0 time since its emergence. The value of T 0 is the absolute time value and it is not known from the images. However, we do not need to know this value in order to evolve the tumor delineation with the model given in Equation ( 4). The formulation given by Equations ( 4) and ( 5) is in a sense a first order approximation. It does not take into account the convergence rate of the initial distribution u to the traveling wave solution -see Figure 2 -and the effect of the curvature of the tumor delineation on its speed. These effects can be included, leading to the final traveling time formulation for the tumor delineation as follows

(Konukoglu 2009a)     1 ' 1 3 . 0 2 3 4 3 . 0 / | |                  T T e T T eff D      (6)           x x D D 0 ) ( , ' T T T T T eff  , where the term       T T 2 / 3 4 
is the effect of the tumor rate of convergence and the term with κ eff is the effect of the tumor delineation's curvature on its evolution. The traveling time formulation given in Equation ( 6) describes the evolution of the tumor delineation in the MR images based on the same growth dynamics as the reaction-diffusion models. In the formulation given in Equation ( 6) we notice the T dependence of the equation. This means that the value of T 0 becomes important for the simulation. As we have noted this value is cannot be known in the clinical routine this seems to pose a problem. However, we solve this problem by treating this value as another model parameter to be optimized for.

In Figure 3 we show an example evolution simulated using the traveling time formulation to show that it captures the same growth dynamics as the reaction-diffusion model given in Equations ( 1) and ( 2). We compare the evolution of a synthetic tumor delineation simulated by the reactiondiffusion model (in white) and by the traveling time formulation (in black). In the case of the reaction-diffusion model the synthetic delineation is obtained by thresholding the tumor cell density distribution, which would not be available in patient cases.

[PLACE FIGURE 3 In addition to these, in the previous section we added another parameter T 0 in the traveling time formulation as a result of integrating the convergence characteristics of the reaction-diffusion solutions. This results in 4 parameters to estimate for. Our aim is to optimize these parameters such that the evolution simulated using the traveling time formulation best matches the real evolution observed in the time series of MR images.

As the first step we define a discrepancy measure between the simulated tumor evolution and the observed one for a given set of parameters. Minimizing this discrepancy then would provide us the optimum model parameters. Our strategy in defining this discrepancy is to use the symmetric surface distances between the simulated and the real delineations and also add a criterion regarding the initial size of the tumor in the first acquired image. The resulting criterion becomes a quadratic measure as follows:

           1 1 2 0 min min 0 ) , dist( , , , N i i g w T T v T d d C  (7)   0 0 0 ) ( with ) ( | ˆ         x x x x T T t T T i i (8) ) ( , 1 . 0 min min min max min max x )n D(x ' n T T v    (9)
where N is the number of images, Γ i is the delineation enclosing the tumor in the i th image, dist() is the symmetric distance between two surfaces, Δt i is the time difference between the i th and the first MR image, x min is the point where the simulated T reaches its minimum and n max is the principal eigenvector of D(x min ). In Equation ( 7 We define the parameter estimation problem as the minimization of the error measure defined in Equation ( 7). Different optimization algorithms can be used for this purpose. Since the gradients of C are not trivial to find analytically we prefer to use a general algorithm which approximates the gradient directions. In this work we use the algorithm proposed by Powell in (Powell 2001) which finds the gradient directions by fitting quadratic surfaces to the underlying minimization surface. Computation times depend on the size of the tumor and the number of images. As an example for a high grade glioma with a time series of 4 images of size (256x256x53) the minimization takes around 50 minutes on a 2.4 GHz Intel Pentium machine with 1 Gb of RAM.

In Figure 4 we provide an example showing the inputs and the outputs of the parameter estimation method.

[PLACE FIGURE 4inputs and the outputs of the parameter estimation method]

Extrapolating Infiltration Extent and Adaptive Irradiation Margins

In the previous section we presented a method to personalize the reaction-diffusion type tumor growth models using the MR images of the patient. Here we assume that we have the personalized model and we present a method to extrapolate the cell density distribution of the tumor beyond its delineation in the image. This method produces a map of possible tumor infiltration that is not visible in the image. Then using this map we describe a method for constructing variable irradiation margins taking into account the growth dynamics of the tumor.

Extrapolating Invisible Infiltration Extents of Gliomas

In order to formulate the extrapolation method first we need an imaging model for the gliomas. We assume that this imaging process can be modeled using a Heaviside function as done in (Tracqui 1995) and in (Swanson 2000). The imaging function Im is given as

      0 0 0 1 )) , ( Im( u u if u u if t u x , (10)
we define 1 as enhanced, 0 as non-enhanced regions in the image and u 0 is the tumor cell density threshold of the imaging modality. This simplistic model assumes that the delineations in the images correspond to an iso-density contour of the tumor cell density distribution. A more general imaging function can be used and similar ideas presented below would be applicable. Based on this definition the problem of extrapolating cell density distribution of a tumor beyond its delineation in the image is defined as

  0 ) Im( | ) , ( ) ( ~0     x x x x x T u u , (11)
where u ~approximates the actual tumor distribution at a time instant T 0 . Unlike the forward modeling of the tumor growth, the construction of the approximation u ~is a static problem and it does not involve the time evolution of the tumor. Moreover, in the clinical situations, the value of infiltration that is not visible in the images. We can use this to construct irradiation margins adapted to this map. Our strategy is to first construct the conventional constant margin irradiation region M c and then create the adapted margins M v by molding M c such that M v takes into account the infiltration extent constructed as explained in the previous section. We refer the reader to (Konukoglu 2009c) for the details of this construction. In Figure 6 we demonstrate the two margins, the constant and the variable ones. For a virtually grown tumor we show the visible part of the tumor, its invisible infiltration and the constant and the variable irradiation margins constructed based on its visible part.

[PLACE FIGURE 6variable irradiation margin]

RESULTS

Parameter Estimation and Model Personalization

The evaluation of parameter estimation for tumor growth models using real patient images is not easy because we do not have access to the real values of the parameters. The real values could be found using microscopic in-vivo analysis however, up to the best of our knowledge such a study has not been performed yet. In this work we perform an indirect evaluation for the proposed parameter estimation method using patient images. For a given patient dataset, we estimate the parameters using all but the image taken at the last time point. Then using the estimated parameters we simulate the evolution of the tumor delineation for the time elapsed between the 2 last images.

We then compare the evolution predicted using the estimated parameters and the travelling time formulation with the one observed in the last image. The correlation between the prediction and the observed delineation provides us with a qualitative evaluation of the estimated parameters.

In the parameter estimation process it is shown in (Konukoglu 2009a) that estimating all the parameters of the model independently using the evolution of the tumor delineations in the images results in non-unique solutions. Therefore in this study we fix the proliferation rates ρ of the tumors to the values suggested in the literature. We only estimate for the diffusion rates and the initial time estimates. In fixing the value of ρ we assume that this value can be estimated based on the biopsy results and microscopic analysis using the relation between the mitotic index (MI) and the labelling index (LI) (Johannessen 2006).

As a preliminary step we use two patient datasets which include anatomical and diffusion tensor MR images. The dataset for the first patient, who suffers from a high grade glioma We observe in Figure 7(c) that the prediction of the tumor delineation is in very good agreement with the final state of the tumor. In the case of the low grade tumor shown in Figure 8, the correlation between the predicted tumor delineation and the final state of the tumor is in line with our previous arguments. We observe that the slow evolution of the tumor is well captured by the estimated parameters. For the proliferation rate we pick a lower value than the one in the previous case since it is a lower grade tumor.

Extrapolating Infiltration Extent and Adaptive Irradiation Margins

In this section we assess the quality of the extrapolation method and the method for constructing the variable, adaptive irradiation margins. In order to understand the potential benefits of extrapolating the infiltration extents of gliomas and constructing irradiation margins adapted to these results we perform experiments on synthetic images. We create a dataset of virtually grown tumors using the reaction-diffusion models where the tumors are grown in the MR images of a healthy subject (1x1x2.6 mm 3 resolution). Different tumors in the dataset are created using different model parameters and different seed points in the brain. For each tumor we construct synthetic images using the Im() function given in Equation 10 with u 0 =0.4, value consistent with (Tracqui 1995). We assume that the detection and the first image acquisition take place when the average diameter of the visible tumor reaches 1.5 cm. After the detection for each tumor we construct a synthetic image every 50 days for 1 year (8 images in total). An example of the created synthetic image with the synthetic tumor is shown in Figure 6. For the tumors in the dataset we construct the constant and the variable irradiation margins based on their visible parts in the images. Since for the virtual tumors, the cell density at every location is knowneven if the image does not show tumor at that pointwe carry out a quantitative comparison. We geometrically compare these margins based on two criteria:

1. R: number of tumor cells not targeted and 2. Vol: volume of healthy tissue targeted by the irradiation margin.

In order to compute these values we follow the values given in (Tracqui 1995) we assume that a voxel of 1x1x2.6 mm 3 can hold a maximum of 9.1x10 4 tumor cells. At the time of detection of high grade gliomas, isolated tumor cells can be found in any region in the brain. Therefore, there is no completely healthy brain. In order to compute the Vol value we need to define "healthy tissue". In this work we define a voxel to be healthy if there are on the average less than 1 tumor cell in it. In Figure 9 for a synthetic tumor we show the comparison between the constant and the variable irradiation margins. We plot the R vs. time and Vol vs. time graphs. We observe that the variable irradiation margins adapted to the extrapolated infiltration extents targets more tumor cells and less healthy tissue.

[PLACE FIGURE 9 The experiments presented above use synthetic images and virtually grown tumors.

Validating the presented methods on real images would require the knowledge of tumor cell density distributions throughout the brain. This could be obtained with post-mortem analysis or with animal models. This would be the topic of a future study.

CONCLUSION

In this chapter we presented methods to personalize reaction-diffusion type tumor growth models in MR images. Our aim was to bridge the gap between the mathematical tumor growth models and the clinical applications. In this spirit, first we focused on a method for personalizing the generic reaction-diffusion growth models based on time series of MR images. The method described the evolution of the tumor delineation in the images based on the growth dynamics captured by the reaction-diffusion models. Based on this formulation it optimized for the model parameters such that the resulting simulations best matches the observed evolution of the tumor boundaries in the MR images. The results shown in Section 3 demonstrated promising preliminary findings.

As a second part, we presented an application of personalized tumor growth models for radiotherapy to demonstrate clinical relevance of the personalization process. We constructed the (Tovi 1994) Tovi., M., Hartman, M., Lilja, A., andEricsson, A. 1994. MR imaging of cerebral gliomas. Tissue component analysis in correlation with histopathology of whole-brain specimens. Acta Radiologica 35:495-505. (Tracqui 1995) Tracqui, P., Cruywagen, G., Woodward, D., et al. 1995

Figures & Tables

  comparison of the traveling time with reaction-diffusion] 2.1.2 The Parameter Estimation Formulation Since we have linked the reaction-diffusion model and the evolution of the tumor delineation through traveling time formulation, we can formulate the parameter estimation problem using this link. In the reaction-diffusion model given in Equations (1) and (2) there are three different parameters: d w , d g and ρ.

  ) the summation measures the distance between the simulated and the real observations. This is obtained by running the traveling time formulation outwards starting from the delineation of the tumor obtained in the first image and comparing the corresponding simulated and observed delineations for the other images. The other term in the account the size of the visible tumor in the first image and quantifies the coherency of this size with the model parameters. We compute this value by simulating the traveling time evolution within the tumor delineation in the first image and compare the minimum T value obtained with the T 0 parameter. If all the parameters of the model are coherent with the observations then this term should equal to zero.

(

  Glioblastoma Multiforme), includes T1-post gadolinium MR images (with the resolution of 0.5x0.5x6.5 mm3) at three successive different time points and diffusion tensor MR image (with the resolution of 2.5x2.5x2.5 mm3) taken at the second time point. The second patient suffers from a low grade glioma (second grade astrocytoma) and the dataset for this patient includes T2 flair MR images (with the resolution of 0.5x0.5x6.5 mm3) at 5 successive time points and the DT-MRI image (with the resolution of 2.5x2.5x2.5 mm3) taken at the first time point. For both cases the tumor boundaries were manually delineated by an expert in each image separately.[PLACE FIGURE 7high grade patient] [PLACE FIGURE 8low grade patient]The images used to estimate parameters, the estimated parameters and the predicted evolution of the tumor delineations along with the real delineations are given in Figures7 and 8. In the images in both Figures, first we show the anatomical images at the time of detection and the intermediate images used in the parameter estimation. On the intermediate images we also plot the manual delineations for the underlying image (white contour) and the simulated evolution of the tumor delineation with the estimated parameters (dark contour) obtained in the course of estimation.Following this we start from the last image (in time) used in the parameter estimation and predict the evolution of the tumor delineation until the acquisition of the final image (which was not used in the estimation). In the corresponding images we show the anatomical MR image taken at the last time point showing the final state of the tumor along the tumor delineation predicted using the estimated parameters drawn as the dark contour. In the accompanying tables we provide the values of the estimated parameters.

  showing the R vs. time and Vol vs. time graph for an example tumor. Give the parameter values in the caption.]

Fig3:Fig4:Fig5:Fig6:Fig7:Fig8: 20 Fig9:

 20 Fig1: (Left) T2-weighted and (Middle) T1-weighted post-gad MR images showing a high grade glioma. The bulk tumor and the infiltrated edema are the enhanced regions in the left and the middle images respectively. (Right) Hypothetical distribution of tumor cell density is a smooth transition starting from the bulk tumor extending beyond the edema. Fig2: Any initial condition with compact support when evolved with reaction-diffusion equations (in the infinite cylinder and constant coefficients) converges to the traveling wave solution. (Left) A Heaviside initialization is evolved in time. (Middle)We plot the distribution at different times on the moving frame. (Right) The traveling wave solution has an asymptotic speed shown in dashed curve however, when we observe the speed of the tumor front (u=0.5 iso-density contour) we notice the low rate of convergence (solid curve). An approximation of the tumor front's speed including time convergence and curvature we get a better approximation (point-dashed curve). Fig3: Comparison between the reaction-diffusion model and the Eikonal Approximation: The temporal evolution of the iso-density contour is demonstrated for a synthetic tumor. Contours are shown for days 400, 600, 800, 1000 and 1200 from the innermost to outermost respectively. The synthetic tumor is virtually grown using the reaction-diffusion model. White contours are obtained by thresholding the tumor cell densities at u = 0.4 for the respective day values (400-600-800-1000-1200). Then in order to simulate the evolution of the iso-density contour (assumed to correspond to tumor delineation in real images) starting from day=400, without the knowledge of the tumor cell density distribution we use the traveling time formulation. Black curves are the contours we obtain at days 600 (2nd innermost) to 1200 (outermost). The tumors were grown in the images of a healthy subject for whom we also have the DT-MRIs. Parameters: (dw = 0.25 mm2/day,dg = 0.01 mm2/day,ρ = 0.012

  Figure 1

*This work was partially supported by Microsoft Research Cambridge

is not available.

The ability to personalize reaction-diffusion models give us the opportunity to construct u for patient images. The asymptotic properties of the reaction-diffusion equations explained in Section 2.1 allow us to write the following formulation

In this equation Γ is the tumor delineation in the image either found by a segmentation algorithm or drawn manually. For the derivation of this formulation we ask the reviewer to refer to (Konukoglu 2009c). Equation ( 12

Adaptive Irradiation Margins

In conventional radiotherapy the irradiation margins are constructed based on the tumor geometry visible in the medical images. The margin takes into account the enhanced area of the tumor in the image plus a constant margin around the delineation to deal with the low cell density infiltration (Kantor 2001). This approach assumes that the tumor cells would diffuse within the brain tissue homogeneously around the visible tumor. It does not take into account the growth dynamics of the tumor especially the differential motility (Giese 1996).

The extrapolation method presented in the previous section gives us a map of possible tumor possible tumor infiltration map of gliomas starting from their visible delineation in the MR images.

Using this map we proposed variable irradiation margins which are adapted to the modeled growth dynamics of the observed tumor. We also showed potential benefits of such variable irradiation margins compared to constant ones.

Although a preliminary validation was performed a thorough validation remains to be done on animal models and post-mortem analysis.

ACKNOWLEDGEMENTS

We