A Statistical Shape Model of the Heart and its Application to Model-based Segmentation Sebastian Ordasa , Estanislao Oubela , Rub´en Letab , Francesc Carrerasb and Alejandro F. Frangia a Computational

Imaging Lab., Pompeu Fabra University, Barcelona, Spain b Hospital Sant Pau, Barcelona, Spain ABSTRACT

In the present paper we describe the automatic construction of a statistical shape model of the whole heart built from a training set of 100 Multi-Slice Computed Tomography (MSCT) studies of pathologic and asymptomatic patients, including 15 (temporal) cardiac phases each. With these data sets we were able to build a compact and representative shape model of both inter-subject and temporal variability. A practical limitation in building statistical shape models, and in particular point distribution models (PDM), is the manual delineation of the training set. A key advantage of the proposed method is to overcome this limitation by not requiring them. Another one is the use of MSCT images, which thanks to their excellent anatomical depiction, have allowed for a realistic heart representation, including the four chambers and connected vasculature. The generalization ability of the shape model permits its deformation to unseen anatomies with an acceptable accuracy. Moreover, its compactness allows for having a reduced set of parameters to describe the modeled population. By varying these parameters, the statistical model can generate a set of valid examples. This is especially useful for the generation of synthetic populations of cardiac shapes, that may correspond e.g. to healthy or diseased cases. Finally, an illustrative example of the use of the constructed shape model for cardiac segmentation is provided. Keywords: statistical shape model, active shape model, model-based segmentation, computational heart anatomy, heart atlas

1. INTRODUCTION Computational or in silico models of the heart will offer a solid underpinning for research and clinical practise in tomorrow’s cardiology. Their importance stems not only from being a powerful tool to integrate theoretical principles with empirical observations, but to study the function of the heart in terms of anatomy, hemodynamics, tissue mechanics, histology, cellular function, and metabolism as well.1, 2 Because of their complexity, the development of a computational model of the heart exceeds the limits of a single discipline. Therefore, important and very active research is currently being conducted by integrating multi-disciplinary efforts. This paper presents a statistical shape model of the whole heart aimed at providing domain knowledge to different applications, ranging from cardiac image analysis to patient-specific multi-physics modeling and simulation. The shape model was built from a training set comprising 100 Multi-Slice Computed Tomography (MSCT) scans of pathologic and asymptomatic patients tailored to coronary artery disease (CAD) assessment (CT Angiography, CTA). Each scan includes 15 temporal phases. The resulting 1500 image volumes exhibit a large amount of both inter-subject and inter-phase (dynamic) variations. One of the devised applications of the constructed shape model is to provide the basis for the application of statistical constraints for dynamic and three-dimensional (3-D) algorithms to segment the heart in the different imaging modalities nowadays available (MSCT, MRI, 3-D US, SPECT and PET). These technologies have different fields of view, as well as varying spatial and temporal resolution. In order to tailor the shape model to the application at hand, a number of combinations of subparts can readily be defined. For instance, while a complete shape model of the heart can be used for functional analysis in CTA, a simplified configuration with Correspondence to Alejandro F. Frangi: [email protected], CILab, Pompeu Fabra University, Pg. de Circumval.laci´ o 8, E080003, Barcelona, Spain. Telephone: +34 93-542-1451

Figure 1. Shape model construction algorithm. Tg : global transformation, Tl : local transformation, RS: Reference Sample, RCS: Reference Coordinate System, NCS: Natural Coordinate System, PCA: Principal Component Analysis

the two ventricles would be appropriate for MRI, a single-ventricle model for 3-D US and SPECT, or a combined left and right atrial model for intra-operative guidance in electrophysiology. The challenge is thus to define each subpart of the model in a convenient way in order to allow for simplified but still valid model configurations.

2. METHOD The rationale behind the construction of the statistical shape model is to warp all training images onto an average heart anatomy on top of which all cardiac structures that are needed in the model can be specified. Subsequently, these structures are transferred (or propagated) to each sample by inverting those warping transformations.

2.1. Building Procedure A first challenging issue is to construct the average anatomy. If an individual sample is used, the warping transformations would be biased to the chosen sample and to similar samples. This would imply having good matchings towards those, but deficient for others. Therefore, a more complex mechanism must be employed. In our work the construction of an average anatomy, in the sense of having the smallest average residual transformation to all the training images, involves an iterative procedure similar to the one proposed by Frangi et al.3 The deformations required to establish dense correspondences between the average anatomy and the training samples were successively obtained with the use of rigid4 and nonrigid volumetric image registration.5, 6 The nonrigid

Figure 2. Atlas building procedure. (a) visual comparison between the reference sample (RS) and the average heart in the reference coordinate system (RCS). (b) Average heart representation in the Reference Coordinate System (RCS) and the Natural Coordinate System (NCS) (after the first iteration of its building procedure). The vertical and longitudinal lines in the center of the image help to depict the effect of the transformation between the two images (see Sec. 2.1 for details)

registration algorithm based on free-form deformations (FFDs) and normalized mutual information (NMI) that in the work of Rueckert et al.7 was used to build statistical deformation models of the brain in MRI, could also cope with the large inter-subject and temporal variability of our population of hearts in MSCT. A practical limitation in building statistical shape models, and in particular 3-D PDMs like in the aforementioned work,3 is the manual delineation of the training set. A key advantage of the proposed method is to overcome this limitation by not requiring them. Again, like in the work of Rueckert et al.,7 it was possible to operate on the intensity images directly. Figure 1 offers a diagram of the complete algorithm, which is composed of the following 10 steps: 1: Choose any sample from the training set. This will be denoted Reference Sample (RS) and constitutes the first average heart anatomy. 2: All samples are rigidly and nonrigidly aligned to RS. The rigidly aligned versions of the samples are used in the following steps of the algorithm. 3: A new heart is constructed by averaging the nonrigidly registered images. At this point, the average heart is said to be in a Reference Coordinate System (RCS), i.e. in the coordinates system of RS, and quite biased to it as can be seen in Figure 2-(a). 4: A new set of transformations is recalculated from the RCS average heart to each rigidly aligned sample. 5: The obtained transformations are inverted and averaged, and the resulting transformation is applied to the average heart in RCS, yielding the so-called Natural Coordinate System (NCS) representation,8 which is the anatomy requiring the least residual nonrigid deformation to explain the anatomical variability across all individual samples. Figure 2-(b) compares the two average hearts. 6: Under the assumption of perfect nonrigid registrations, the NCS average heart as yet obtained is independent to the RS election. However, to make the average representation less sensitive to the initial RS election, a new set of nonrigid registrations is recalculated, but this time from the NCS coordinate system to each rigidly aligned sample. Again, the obtained transformations are inverted, averaged, and applied to the heart in the NCS system. 7: The previous step is repeated a number of times until the NCS average heart does not change. 8: The NCS heart is manually segmented and labeled. Each subpart is triangulated separately as explained in the next section. 9: The transformations corresponding to the final iteration are propagated to the rigidly aligned samples. 10: Finally, the statistical shape model is built by applying PCA on the set of propagated shapes. In synthesis, any autolandmarked sample shape (S i ) of the training set can be obtained by S i = Tg−i (Tl−i (A)) where Tg−i and Tl−i are its global and local inverse transformations, and A is the landmarked average heart.

Figure 3. Shape model subparts. From left to right: LV myocardium, RV myocardium, aorta, left atrium, right atrium and fibrous skeleton. Table 1. Model Subparts Subpart 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Main subpart LV Myocardium LV Myocardium RV Myocardium RV Myocardium RV Myocardium LA Myocardium LA Myocardium LA Myocardium LA Myocardium LA Myocardium LA Myocardium LA Myocardium RA Myocardium RA Myocardium RA Myocardium RA Myocardium RA Myocardium RA Myocardium Aorta Fibrous skeleton

Acronym LVEN DO LVEP I RVEN DO RVEP I RV OT LAEN DO LAEP I LAP V LP P V RAP V RP P V LAA RAEN DO RAEP I SV C IV C CS RAA AAo FS

Name LV endocardium LV epicardium RV endocardium RV epicardium RV outflow tract LA endocardium LA epicardium Left Anterior Pulmonary Vein Left Posterior Pulmonary Vein Right Anterior Pulmonary Vein Right Posterior Pulmonary Vein LA Appendage RA endocardium RA epicardium Superior Vena Cava Inferior Vena Cava Coronary Sinus RA Appendage Ascending Aorta Fibrous skeleton

Number of Nodes (total: 30964) 4500 4500 3500 3500 1500 2198 2386 234 243 277 222 379 2250 2350 649 324 367 487 3650 2130

2.2. Labeling After convergence of the iterative procedure, the average heart was manually segmented to generate 6 main constitutive parts (see Figure 3). Each of these regions was triangulated using AMIRA 3.0 (Mercury Computer Systems Inc., USA) and further optimized with the hierarchical mesh refinement method available in NetGen 4.4 (Johannes Kepler University, Linz, Austria) for volumetric mesh generation. By subdividing each main part, the final multi-surface representation results in 20 subparts as summarized in Table 1. The pulmonary veins were considered up to approximately 2 cm from their ostium in the left atrium. First segments of superior and inferior vena cava, and coronary sinus, were set in a similar fashion. The RV outflow tract was defined up to the pulmonary valve. The fibrous skeleton includes part of the interatrial septum, the little “membranous spot” at the top of the interventricular septum, and the mitral, tricuspid and aortic fibrous rings. The atrial myocardium, thicker in the left atrium than in the right atrium, was assigned a regular thickness of approximately 2 mm and 3 mm, respectively. In the current implementation of the model, the papillary muscles were excluded. By combining the different subparts of the model, many possible configurations can be assembled. In this work we built 5 model configurations, namely full-heart, 4-chamber, 2-chamber, 1-chamber and lv-la.

2.3. Point Distribution Model (PDM) Consider the set X = {xi ; i = 1 · · · n} of n propagated shapes obtained from the definition of the previously described cardiac components on top of the NCS average heart. Each shape is described by the concatenation of m 3-D landmarks pj = (p1j , p2j , p3j ); j = 1 · · · m. Note that only the points and not their triangulation is considered. To construct the PDM, Principal Component Analysis (PCA) is used to obtain the following

Figure 4. Principal modes of variation of the full-heart model √ for a single individual. Each mode mi is generated by moving the corresponding bi parameter to ±3SD with SD = λi . The other parameters in the b vector are set to zero.

Pn linear model x = x ˆ + Φb, where x ˆ = n1 i=1 xi is the average landmark vector, b is the shape parameter vector of P the model, and Φ is a matrix whose columns are the principal components of the covariance matrix n 1 S = n−1 ˆ)(xi − x ˆ)T . The principal components of S are the eigenvectors, φi , with corresponding i=1 (xi − x eigenvalues, λi (sorted so that λi ≥ λi+1 ). If Φ contains only the first t < min{m, n} eigenvectors corresponding to the largest non-zero eigenvalues, we can approximate any shape of the training set, x, using the aforementioned linear model, where Φ = (φ1 |φ2 | · · · |φt ) and b is a t dimensional vector given by b = ΦT (x − x ˆ). As a prerequisite for a correct shape description, the shapes have to be normalized with respect to a reference coordinate frame. Procrustes alignment9 was used for this purpose. The vector b defines the shape parameters of the PDM. By varying its components different instances of the shape class under analysis can be generated. Assuming that the cloud of landmark vectors follows a multi-dimensional Gaussian distribution, the variance of the i-th √ parameter, bi , across the training set is given by λi . By applying limits to the variation of bi , |bi | ≤ ±β λi , with β typically 2-3, it can be ensured that a generated shape is similar to the training shapes. Figure 4 shows the first four principal modes of variation of the full-heart representation for a single individual (15 samples, i.e. one sample per phase). The main characteristic variations consist in relative size differences between subparts, twisting, blending, and their combinations. The 1st mode describes the vertical displacement of the base and the 2nd mode the expansion and contraction of the ventricles and atria. Less vigorous modes depict shape variations of the RV, RVOT and connected vasculature of atria.

3. EXPERIMENTS 3.1. CTA Training Set The training set used to construct the shape model comprises 100 individuals (age 58 ± 8 years, 56% men) in whom MSCT was performed for non-invasive evaluation of the coronary arteries (CTA studies). The population included 60 subjects who had normal coronary arteries, 20 patients with significant coronary artery disease (CAD) without a history of previous myocardial infarction (MI) and 20 patients with CAD and a history of MI. Each scan includes 15 temporal phases with a voxel dimension of 0.4 x 0.4 x 2.0 mm3 . Imaging was performed with a 64-row detector Toshiba Aquilion 64 system (Toshiba America Medical Systems, USA), operating at Clinica Creu Blanca (Barcelona, Spain). Between 80 and 100 ml of contrast material (Xenetic 350) at an injection rate of 5 ml/s was used. Rotation time ranged from 400 to 500 ms depending on heart rate, and tube voltage was 120 kV at 400 - 430 mA.

Shape models for different groups of samples can easily be constructed. For the experiments we built both, shape models describing the inter-subject variability at different time frames (with 100 samples each) as well as a more general one combining the inter-subject and temporal variability in a single model (with 1500 samples). We will refer to them as single-phase and multi-phase, respectively.

3.2. Computing Issues The algorithm was run on a SGI Altix 350 system with SGI Suse Linux 9.2. The system contains 28 processors and 30 GB of shared memory. Each processor is a 64 bit Intel Itanium at 1.5 GHz. All the registrations (rigid and nonrigid) were obtained with the Image Registration Toolkit∗ , which is an implementation of the previously referenced methods.4–6 The registration tasks were distributed among the processors using Perl scripting language. With this facility, the construction of a set of rigid followed by nonrigid registrations took for the total of 1500 training samples, approximately 45 hours The building procedure has shown to be convergent. After four iterations, the agreement was excellent (NMI = 4.0) and hardly any visual difference between images in successive iterations could be identified. The average heart of the fourth iteration was used in this work. Therefore, adding the time for the 4 iterations, the total processing time was 180 hours (7.5 days).

3.3. Statistical Analysis To explore the statistical behavior of the constructed shape model we have assessed its compactness capacity and generalization ability as usually assessed in shape analysis studies.10 We found that these characteristics were quite similar among single-phase models and marginally dependent on the number of nodes (as long as the global shape is reasonably preserved). Therefore, only the analysis of the multi-phase and a single-phase (end diastole) models is reported. 3.3.1. Compactness Capacity A compact model is one that requires as few parameters as possible for the generation of a valid instance of the modeled object. The compactness capacity of the shape model, denoted C(M ) is measured as the cumulated variance (absolute or percentile with respect to the total shape variance) for the first M modes, PM C(M ) = m=1 λm . Figure 5-(a) compares C(M ) between the single-phase and multi-phase models. The curves differ due to the different number of samples (1500/100 = 15 times) and hence, the total variance in each shape model. For instance, to explain 98% of shape variance, 49 and 73 modes were necessary for the single-phase and multi-phase models, respectively. On the other hand, the more subparts are incorporated in the shape model, the larger the required number of modes to explain a given percentage of the total variance. Figure 5-(b) compares C(M ) for the different model configurations. 3.3.2. Generalization Ability The property of generalization of a shape model measures its ability to represent unseen instances of the object class modeled. It is explored by evaluating the ability to reconstruct a sample not included in the training set (leave-one-out test), increasingly employing more modes. The reported reconstruction error, denoted G(M ), is calculated as the point-to-point (p2p) distance averaged over all shapes. Figure 5-(c) shows the generalization ability of the single-phase model. The behavior of other phases and mesh densities was very similar. Note that G(M ), although adequately low (below 2 mm), does not stabilize with larger number of modes. This suggests that 100 subjects is not yet large enough to completely model the shape variation in the corresponding cardiac phase. Therefore, a larger training set could still improve the generalization ability of the model. It is important to point out that in the leave-one-out experiments, none of the samples corresponding to the “left out” individual were included in the corresponding “left in” training set. In this way, the sample to be reconstructed was completely unknown to the shape model. ∗

http://wwwhomes.doc.ic.ac.uk/∼dr/software/

Figure 5. Statistical analysis of the shape model. (a) compaction capacity between different populations (single-phase and multi-phase). (b) compaction capacity between different model representations (full-heart, two-chamber, one-chamber and lv-la). (c) generalization ability (error ± standard deviation).

Table 2. Point Correspondence Accuracy: intra-observer and inter-observer variability (distance in millimeter) of manual landmark definition on samples, onto the average heart, and after propagation to each image sample. Lmk #

1 2 3 4 5 6

Variability of landmarks manually defined on samples inter intra intra (Obs.1 vs. (Obs.1) (Obs.2) Obs.2) 2.38 ± 1.45 1.56 ± 1.05 2.19 ± 1.04 3.97 ± 0.99 3.10 ± 1.32 6.26 ± 2.32 2.98 ± 1.69 5.55 ± 3.13 3.79 ± 1.15 2.62 ± 1.14 3.38 ± 2.71 2.95 ± 0.89 1.31 ± 1.62 1.05 ± 0.71 1.29 ± 0.56 1.13 ± 0.84 0.83 ± 0.56 0.97 ± 0.54

Variability of landmarks manually defined on the average heart inter intra intra (Obs.1 vs. (Obs.1) (Obs.2) Obs.2) 1.81 0.61 1.16 ± 0.29 1.71 3.48 4.74 ± 1.23 1.92 1.60 1.35 ± 0.27 1.33 0.36 2.30 ± 0.23 2.08 1.63 1.75 ± 0.81 0.94 0.62 1.54 ± 0.43

Variability of landmarks after propagation inter intra intra (Obs.1 vs. (Obs.1) (Obs.2) Obs.2) 1.70 ± 0.25 0.63 ± 0.10 1.08 ± 0.21 1.57 ± 0.31 3.12 ± 0.96 4.81 ± 0.66 1.77 ± 0.14 1.60 ± 0.07 1.20 ± 0.10 1.32 ± 0.20 0.27 ± 0.06 2.20 ± 0.21 2.19 ± 0.37 1.59 ± 0.25 1.73 ± 0.28 1.00 ± 0.12 0.65 ± 0.13 1.63 ± 0.25

3.4. Point Correspondence Accuracy To evaluate the quality (accuracy and reproducibility) of the point correspondence achieved with the method, two observers were asked to individualize in two independent sessions (separated by 1 week), 6 anatomical landmarks in a total of 20 images of the training set. The selected landmarks were: (1) aortic valve centroid, (2) centroid of the tricuspid valve, (3) centroid of the mitral valve, (4) LV endocardial apex, (5) left coronary root and (6) right coronary root. In each landmarking session, the observers were also asked to identify the same landmarks in the average heart anatomy. Table 2 presents the results. Note from the second column the intrinsic inaccuracy of manual landmark positioning on the samples. Among the explored anatomical landmarks, the centroid of the aortic valve (lmk #1) had the lowest variability, whereas the centroid of the tricuspid valve (lmk #2) was the most difficult to localize due to the little contrast media in the RV, typical of CTA studies. From the third column, the variability of manual landmarking on the average heart was very low because all landmarks were easily identifiable on such a “smoothed” anatomy. Subsequently, the set of landmarks on the average heart was propagated to each subject using the deformation field described in Section 2.1. The fourth column shows the variability of the resulting positions, which has remained low. Table 3 shows the distance between the mean position of each landmark as the result of the manual and automatic approaches. For example, for the case of the aortic valve centroid, the variability for manual definition on the average heart was 1.81 mm (obs.1 intra), 0.81 mm (obs.2 intra) and 1.16 ± 0.29 mm (inter). These errors remained low after propagation: 1.70 ± 0.25 mm (obs.1 intra), 0.63 ± 0.10 mm (obs.2 intra), and 1.08 ± 0.21 mm (inter). For this same landmark, the average distance between the average position by the manual and automatic approaches was 1.78 ± 1.34 mm (in average, for the two observers). This variability is less than both intra- (2.38 ± 1.45 mm, 1.56 ± 1.05 mm) and inter- (2.19 ± 1.04 mm) observer variability of manual definition on the samples.

Table 3. Point Correspondence Accuracy: distance in millimeters between the average landmark position resulting from the manual and automatic approaches. Lmk #

Anatomical Location

1 2 3 4 5 6

aortic valve centroid tricuspid valve centroid mitral valve centroid LV endocardial apex left coronary root right coronary root

distance between average landmark position of: Auto Obs.1 - Manual Obs.1 Auto Obs.2 - Manual Obs.2 Auto Obs.1,2 - Manual Obs.1,2 1.88 ± 1.21 2.16 ± 1.56 1.78 ± 1.34 6.75 ± 2.98 7.33 ± 3.51 6.83 ± 3.16 3.72 ± 1.84 3.64 ± 0.96 3.31 ± 1.22 3.68 ± 1.25 2.76 ± 1.84 2.72 ± 1.80 3.86 ± 2.76 3.64 ± 2.74 3.66 ± 2.66 3.90 ± 2.45 4.18 ± 2.14 3.97 ± 2.33

4. MODEL-BASED SEGMENTATION In this section we describe the inclusion of the constructed shape models in the popular Active Shape Model (ASM) algorithm pioneered by Cootes et al.11 and further extended by several researchers. An ASM is composed of three parts: (1) a shape model like the ones constructed, (2) an appearance model that describes texture patterns around the object’s boundaries in the different imaging modalities and (3) a matching algorithm that drives the shape model deformation towards the target in an iterative manner.

4.1. Appearance Model Since organ boundaries are related to gray-level transitions, a common feature function used is the gradient operator. In the seminal work of 2-D ASMs11 the image structure around each landmark was represented from first derivative profiles perpendicular to the object’s border, sampled with bi-linear interpolation. Taking k positions at each size, the total length is 2k + 1. The effect of global intensity changes is reduced by normalizing the patches such that the sum of the absolute intensity values equals 1. The resulting normalized Gray-level with gik as the k − th component in the i − th profile. Denoting Profile (nGLP) takes the form: gˆi = PL gi k=1

|gik |

these normalized derivative profiles as gˆ1 , ..., gˆs , the mean profile g and the covariance matrix Σg are computed for each landmark. These two constitute the nGLP appearance model for that landmark. The Mahalanobis distance f (ˆ gi ) = (ˆ gi − g)Σ−1 gi − g) between a new profile gˆi sampled on a new image and g (ˆ the corresponding mean profile g in the learning set is used as cost function to minimize in order to displace each landmark to a new position. Minimizing f (ˆ gi ) is equivalent to maximizing the probability that gˆi originated from the same multivariate Gaussian distribution assumed for the training profiles. Σg thus weights the variation among corresponding points in the training profiles. Instead of using the Mahalanobis distance and the gradient operator alone, several other alternatives could be considered,12–14 including non-linear classifiers and specialized local image descriptors.15–18 In any of the aforementioned approaches to construct the appearance model for each node in the shape model, we need to introduce some changes with respect to the 2-D version. In 3-D is impractical to have a gray-level model for each point in the shape model (e.g. more than 25000 for the full-heart and 4-chamber representations). Furthermore, as we will see shortly, adapting the shape model to an image volume involves the intersection of the model instance with the image planes at every iteration of the matching procedure. These intersections will rarely occur on the shape model nodes directly but in some locations in between. The solution we have adopted is to build a mean profile and a covariance matrix of normalized GLPs for corresponding subparts of the model in all training samples. To have smaller regions, the LVM Y O subpart was further subdivided in: (1) apical, (2) septal, (3) lateral, (4) superior, (5) inferior and (6) posterior regions. Such definition corresponds to neighboring organ locations and tissues, and thus changes in texture.

4.2. Matching Procedure We have adapted the constructed shape models to the Sparse Active Shape Model (SPASM) matching procedure of van Assen et al.19 This methodology allows for using sparsely and arbitrarily oriented data. On each iteration of the algorithm, the current shape model instance intersects the image planes at a set of intersection points that do not necessarily belong to the model’s mesh. The intersection points can be assigned by closest proximity to each of the shape model subparts. Following, new candidate displacements are searched in every image

Figure 6. Multi-phase segmentation of a CTA study. In each row, the volume rendering and automated segmentation are clipped at two axial levels. Each column corresponds to a different temporal phase. Only 5 of the 15 available phases are shown (phases 1, 3, 5, 9 and 12).

Figure 7. Multi-phase segmentation of a CTA study. Volume rendering and views from different clipping planes.

slice by evaluating the Mahalanobis distance between the new gray-level profiles found in the image, and the corresponding mean profiles stored for each subpart. In order to match the shape model to the data, the resulting selected points dictated by the appearance model and lying on the intersected image planes, act as forces that are propagated to the rest of the mesh nodes. After force propagation, each node in the mesh has an associated displacement vector. A new model instance is retrieved by deriving the transformation matrix T that incorporates scaling s, rotations θx , θy , θz and translations tx , ty , tz . For this purpose, the Procrustes method9 is again employed to find the parameters that minimize the sum of squared distances (or residues) between the feature points Y and the current model shape x, by minimizing (Y − T [x])T (Y − T [x]) Finally, the shape parameters b are calculated by minimizing the sum of squared distances between the current model shape and the projection of the shape Y onto the model centered coordinate system: (T −1 [Y ] − (¯ x + P b))T (T −1 [Y ] − (¯ x + P b)). These steps are repeated either for a fixed number of iterations or until convergence is achieved, given some suitable decision rule.

4.3. Model Initialization In the matching procedure of ASMs, the first shape model instance is usually the mean shape plus a proper initial position, orientation, and scale. Initialization is necessary for the robustness (capture range) and speed of

convergence of the algorithm. Although the 3-D ASM has shown to be relatively insensitive to initialization,19 when a shape model more complex than the left ventricle is used, like a 4-chamber or 2-chamber configuration, a convenient mechanism must be employed. In our implementation, the operator defines the following anatomical landmarks: (1) centroid of the aortic valve, (2) centroid of the mitral valve, (3) endocardial apex of the LV, (4) centroid of the tricuspid valve (optional), (5) centroid of the pulmonary valve (optional). Corresponding anatomical locations in the mean shape are aligned to these points using a similarity transformation. The use of more points e.g. in the atria, ascending aorta, coronary sinus, etc., did not improve the alignment enough to justify the additional interaction required. For multi-phase segmentation, the initialization of the shape model is obtained from the segmentation of the previous temporal phase. Figure 6 and 7 illustrates two examples of whole-heart segmentation in CTA image series.

5. SUMMARY AND CONCLUSIONS In this paper we have presented a statistical shape model of the whole heart built from a large training set of CTA studies corresponding to 100 asymptomatic and pathologic patients. The iterative procedure employed for its construction has shown to converge after four iterations. The nonrigid registration algorithm based on FFDs and NMI could cope satisfactorily with the large inter-subject and inter-phase variability present in the training population. A total of 20 subparts were defined on the model and assigned an individual label. Depending on the targeted application, several shape model configurations can be assembled by their combination. Experiments were carried out to establish the ability of the different shape model configurations to generalize to samples not present in the training set. The average reconstruction error of completely unseen cases (all temporal phases of the same individual left out) was below 2 mm for a quite complete model configuration (fourchamber ). Moreover, all generated instances of the shape model were anatomically valid, without overlapping or intertwinement between neighboring anatomical regions. Although the reconstruction error would probably continue to decrease with more training images, 100 individuals seem to be representative of the human heart under healthy and pathologic conditions. In the last experiment, we found a good accuracy and precision of the point correspondence between the average heart anatomy and the training samples. An example of the application of the shape model for cardiac image segmentation in CTA was presented. Future research will go into the line of large scale validation of the segmentation robustness and accuracy, and the development of improved appearance models that can cope with texture non-linearities, like those encountered in the neighborhood of the heart in the different imaging modalities.

ACKNOWLEDGMENTS This work has been partially supported by MEC TEC2006-03617, ISCIII FIS2004/40676, and CDTI CENITCDTEAM grants. SO and EO both hold MEC-FPU grants. The work of AFF is supported by the Spanish Ministry of Education and Science under a Ram´on y Cajal Research Fellowship. The CILab is part of the ISCIII CIBER-BBN (CB06/01/0061). For all image registrations, the Image Registration Toolkit was used under Licence from Ixico Ltd.

REFERENCES 1. P. Hunter, A. Pullan, and B. Smaill, “Modelling total heart function,” Annu. Rev. Biomed. Eng. 5, pp. 147– 177, 2003. 2. E. Crampin, M. Halstead, P. Hunter, P. Nielsen, D. Noble, N. Smith, and M. Tawhai, “Computational physiology and the physiome project,” Experimental Physiology 89, pp. 1–26, 2004. 3. A. F. Frangi, D. Rueckert, J. A. Schnabel, and W. J. Niessen, “Automatic construction of multiple-object three-dimensional statistical shape models: Application to cardiac modeling,” IEEE Trans Med Imaging 21(9), pp. 1151–1166, 2002. 4. C. Studholme, D. Hill, and D. Hawkes, “An overlap invariant entropy measure of 3D medical image alignment,” 32(1), pp. 71–86, 1998.

5. D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill, M. O. Leach, and D. J. Hawkes, “Non-rigid registration using free-form deformations: Application to breast MR images,” IEEE Trans Med Imaging 18, pp. 712–721, Aug. 1999. 6. J. Schnabel, D. Rueckert, M. Quist, J. Blackall, A. Castellano Smith, T. Hartkens, G. Penney, W. Hall, H. Liu, C. Truwit, F. Gerritsen, D. Hill, and D. Hawkes, “A generic framework for non-rigid registration based on non-uniform multi-level free-form deformations,” in MICCAI, W. J. Niessen and M. A. Viergever, eds., Lecture Notes in Computer Science 2208, pp. 573–581, Springer, 2001. 7. D. Rueckert, A. F. Frangi, and J. Schnabel, “Automatic construction of 3D statistical deformation models of the brain using non-rigid registration.,” IEEE Trans Med Imaging 22(8), pp. 1014–25, 2003. 8. U. Grenander and M. I. Miller, “Computational anatomy: An emerging discipline,” in Quart. Appl. Math., 56(4), pp. 617–694, 1988. 9. J. Gower, “Generalized procrustes analysis,” Psychometrika 40, pp. 33–50, 1975. 10. M. Styner, K. Rajamani, L. Nolte, G. Zsemlye, G. Szekely, C. Taylor, and R. Davies, “Evaluation of 3D correspondence methods for model building,” pp. 63–75, 2003. 11. T. Cootes, C. Taylor, D. Cooper, and J. Graham, “Active shape models - their training and application,” Computer Vision and Image Understanding 61(1), pp. 38–59, 1995. 12. G. Behiels, F. Maes, D. Vandermeulen, and P. Suetens, “Evaluation of image features and search strategies for segmentation of bone structures in radiographs using active shape models,” Med Image Anal (6), pp. 47– 62, 2002. 13. S. Li, L. Zhu, and T. Jiang, “Active shape model segmentation using local edge structures and adaboost,” in Medical Imaging and Augmented Reality (MIAR), 3150, pp. 121–128, LNCS, 2004. 14. J. Peters, O. Ecabert, and J. Weese, “Feature optimization via simulated search for model-based heart segmentation,” in Computer Assisted Radiology and Surgery – CARS 2005, p. 3338, 2005. 15. B. van Ginneken, A. F. Frangi, J. Staal, B. M. ter Haar Romeny, and M. Viergever, “Active shape model segmentation with optimal features,” IEEE Trans Med Imaging 21(8), pp. 924–933, 2002. 16. M. de Bruijne, B. van Ginneken, W. J. Niessen, A. Maintz, and M. A. Viergever, “Active shape model based segmentation of abdominal aortic aneurysms in CTA images,” in SPIE Medical Imaging, M. Sonka and J. Fitzpatrick, eds., 5032, p. 463474, SPIE, 2003. 17. S. Ordas, L. Boisrobert, M. Huguet, and A. F. Frangi, “Cardiac MRI segmentation using 2D-ASMs with invariant optimal features,” in IEEE Computers in Cardiology, Thessaloniki, Greece, pp. 633–636, 2003. 18. F. Sukno, S. Ordas, C. Butakoff, S. Cruz, and A. F. Frangi, “Active shape models with invariant optimal features IOF-ASMs,” in Audio- and Video-Based Biometric Person Authentication, New York, USA, pp. 365–375, Springer, 2005. 19. H. C. van Assen, M. G. Danilouchkine, A. F. Frangi, S. Ordas, J. J. M. Westenberg, J. H. C. Reiber, and B. P. F. Lelieveldt, “SPASM: a 3D-ASM for segmentation of sparse and arbitrarily oriented cardiac MRI data,” Med Image Anal 10(2), pp. 286–303, 2006.