Geometric heart modeling is not a trivial task since the heart is a complex nonrigid organ, with multiple moving parts. The model must be anatomically accurate, accept manual editing, and provide sufficient information to guide automatic detection and segmentation. Except for a four-chamber heart model from [14] and [31], most of the previous work was focused on the Left Ventricle (LV) and/or the Right Ventricle (RV) . A closed mesh was often used to represent heart chambers [18, 24, 33]. Nevertheless, it is not clear how the atria interacted with the ventricles around the mitral and tricuspid valves in [33]. The heart model in [31] was more accurate in anatomy and it also included trunks of the major vessels connected to heart chambers. However, artificial patches were added at all valves to close the chamber meshes. These artificial patches were not processed with a special care in the segmentation algorithm; therefore, they could not be delineated accurately [14]. In our heart model, we keep the mesh open at a valve. Mesh points at the valve rim are labeled as control points , being treated differently to the normal mesh points during automatic segmentation.

The statistical shape model [9] is used in nonrigid object segmentation to enforce shape constraints and to make the system more robust. However, to build a statistical shape model, it is necessary to establish point correspondence among a group of shapes [9]. There are a few papers on building a 3D statistical shape model automatically using pair-wise or group-wise registration based approaches [16, 32], which are time consuming. Another approach is to establish correspondence among shapes through mesh resampling . Although this is difficult for a generic 3D shape, we can consistently resample the surface to establish correspondence for a few simple shapes (e.g., a tube and a parabola) [19, 24].

Given the heart model, the segmentation task is to fit the model onto an unseen volume. Since the heart is a nonrigid shape , the model fitting (or heart segmentation) procedure can be divided into two steps: object localization and boundary delineation. Most of the previous approaches were focused on boundary delineation using Active Shape Models (ASM)  [2], Active Appearance Models (AAM)  [1, 37], or deformable models  [3, 10, 14, 20, 36, 40]. The above segmentation technologies suffer from two limitations: (1) Most of them are semi-automatic and proper manual initialization is demanded. (2) Gradient descending based search in some of the approaches are likely to converge to a sub-optimal solution .

Object localization, a task largely ignored by early work in heart segmentation, is required for an automatic segmentation system. The MSL computational framework for object detection and pose estimation thus fills an important gap in a 3D heart segmentation pipeline, by efficiently providing automatic detection of the heart parts, and estimating their 3D pose and initial mean shape.

A closed mesh has been used to represent the Left Ventricle (LV) [18, 24, 33]. Due to the lack of object boundary on the image, it is hard to consistently delineate the interfaces among the LV main body, the Left Ventricular Outflow Tract (LVOT) , and the basal area around the mitral valve . The mesh often cuts the LVOT and the mitral valve at an arbitrary position [18, 24, 33]. In our heart model, the mitral valve is explicitly represented as a closed contour along its border. Since we exclude the moving valve leaflets (which are hardly visible) from the model, the basal area can be delineated more consistently. Both the endo- and epi-cardiums are delineated for the LV. The commissure contour of both surfaces corresponds to the mitral valve annulus on one side and the aortic valve annulus (lying at the bottom edge of the Valsalva sinuses ) on the other side, as shown in Fig. 7.2d. Three curves are formed by control points around the mitral valve, namely, the mitral lateral (16 points), the mitral septum (15 points), and the LVOT lateral (16 points). As shown in Fig. 7.1a, they define two closed regions, one for the interface between the LV and the Left Atrium (LA) , and the other for the LV and the LVOT. The aortic valve (annotated with 32 points) is approximated as a plane, which cuts the valve at the bottom of the Valsalva sinuses . The LA is represented as an open mesh with an open area enclosed by the mitral septum and the mitral lateral control points (as shown in Fig. 7.1b). Figure 7.2 shows the LV/LA meshes with the control points (dots that are connected appropriately to form contours in the figure). The LV endocardium , epicardium , and LA meshes are represented with 545 points and 1,056 triangles each. The LVOT mesh is represented with 64 points and 64 triangles.

The Right Ventricle (RV) has a complicated shape with separate inflow and outflow portions . Using a plane (indicated by a horizontal black line in Fig. 7.3a) passing the divergence point of the inflow and outflow tracts, we can naturally split the RV into three parts with the RV main body lying below the cutting plane. We call this plane “the RV divergence plane.” During the manual labeling of ground truth, the RV divergence plane is determined in the following way. We first determine the tricuspid valve , which is approximated as a plane. We then move the plane toward the RV apex to a position where the RV inflow and outflow tracts diverge. As shown in Fig. 7.3b, two curves, tricuspid lateral (23 points) and Right Ventricular Outflow Tract (RVOT) lateral (15 points), are annotated on the RV divergence plane to define the inflow and outflow connections. On a short axis view, the RV main body is a crescent (as shown in Fig. 7.3b and c). Two cusp points on the intersection are important landmarks, and an automatic detection algorithm is able to detect them reliably. They are explicitly represented in our model (Fig. 7.3d). The tricuspid (28 points) and pulmonary (18 points) valves are approximated as a plane. Similar to the LA, the Right Atrium (RA) is represented as an open mesh with the open area defined by the tricuspid valve. Figure 7.4 shows the RV/RA meshes with the control points (dots that are connected appropriately to form contours in the figure). In our model, the RV is represented with 761 points and 1,476 triangles and the RA is represented with 545 points and 1,056 triangles.

After object pose estimation, we use the aligned mean shape as an initial estimate of the true shape. To calculate the mean shape, we demand the annotated training meshes to have intrinsic anatomical correspondence across volumes. Furthermore, we use the Active Shape Model (ASM) to enforce the shape constraint during boundary evolution. To build the statistical shape model of the ASM, we also need the point correspondence . Establishing mesh point correspondence is a difficult task for a generic 3D shape. Fortunately, for a few simple shapes, such as a tube or a parabola, we can consistently resample the surface to establish this correspondence. Since it is much easier to consistently resample a 2D curve, we use a few planes to cut the 3D mesh to get a set of 2D intersection curves. The resulting 2D curves are uniformly sampled to get a point set with built-in correspondence. Using different methods to select cutting planes, we develop two resampling schemes, the rotation-axis based method for simple parabola-like shapes such as the LV , LA , and RA , and parallel-slice based method for the more complicated RV. In both methods, the long axis of a chamber needs to be determined from the annotated mesh. Generally, we define the long axis as the line connecting the center of a valve and the mesh point farthest from the valve (which often corresponds to the apex). For example, the LV long axis is the line connecting the mitral valve center and the LV apex . The RV long axis is determined as a line passing the RV apex and perpendicular to the divergence plane .

The rotation-axis based method is appropriate for a roughly rotation symmetric shape. Cutting the mesh with a plane passing the rotation axis, we get a 2D intersection. As shown in Fig. 7.5b, the rotation axis separates the intersection into two parts, which can be uniformly resampled independently. We then rotate the plane around the axis to get another intersection and resample the intersection in the same way. Repeating the above process, we achieve a set of points with built-in correspondence. We use this approach to resample the LV, LA, and RA, and also the LVOT , and RV inflow and outflow tracts , which can be approximated as tubes. To establish the correspondence between the cutting planes, we need a consistent way to determine the first cutting plane, which is often based on an anatomical landmark. Different chambers may use different landmarks. For the LV, we pick the cutting plane passing the aortic valve center as the first cutting plane.

The rotation-axis based method is not applicable to the RV main body since it is not rotation symmetric. Instead, a parallel-slice based method is developed for the RV, where we use a plane perpendicular to the RV long axis to cut the 3D mesh, as shown in Fig. 7.6a. The shape of an RV short-axis intersection is a crescent containing two cusp points (which have a high curvature and can be reliably determined from the intersection contour). They split the contour into two parts and each can be uniformly resampled, as shown in Fig. 7.6c.

Compared to other approaches [11, 16, 32] for establishing point correspondence, our solution is simple and each shape is processed independently. No time-consuming and error-prone 3D registration is necessary. Although our approach is not as generic as [11, 16, 32], it provides superior results for the heart chambers.

The statistical shape model [9] (a.k.a., the point distribution model ) is widely used in computer vision and medical image analysis to enforce prior shape constraint during the segmentation of a nonrigid object. The statistical shape model is learned from a set of training shapes, . Each shape M _i is composed of the same number of points J and can be represented as . Point is assumed to have correspondence across volumes. That means it represents the same anatomical position in different volumes. The point correspondence can be established by various methods, such as the mesh resampling techniques presented in the previous section.

Before performing shape analysis , the external transformations of the training shapes need to be compensated. That means we need to transform the shapes to the same coordinate system. This process, called Generalized Procrustes Analysis (GPA)  [13], searches for a mean shape and an optimal transformation for each individual shape M _i to minimize the sum of residual errors after alignment,
The similarity transformation is used as T in GPA and is often solved using an iterative approach. Please refer to Sect. 6.​3.​2 of Chap. 6 for an in-depth discussion on GPA. Figure 7.7 shows the RV shapes before and after the GPA.

Suppose after alignment a training shape M _i is transferred to m _i . The mean shape over the whole training set after alignment is
The residual, d m _i , of each shape to the mean shape is
The covariance matrix , which is a 3J × 3J symmetric matrix, of the aligned shapes is calculated as
Using eigen-decomposition , the covariance matrix can be decomposed as
where λ _i is an eigenvalue and P _i is the corresponding eigenvector. Since is a real symmetric matrix that is semi-positive definite, all of its eigenvalues are non-negative real numbers. Without loss of generality, suppose the eigenvalues are ordered as . In the statistical shape model , the linear space represented by P is often called the shape space and each axis of the space, P _i , is called a deformation mode . Its corresponding eigenvalue λ _i is the amount of the variation of this deformation mode. The total variation λ _T is

A shape in the training set can be presented as the mean shape plus a weighted sum of the deformation modes,
where is a vector of weights. Since eigenvectors P _i are orthonormal to each other, the weight w _i can be easily calculated by the dot product of vectors m and P _i ,
In general, a shape out of the training set cannot be represented exactly with Eq. (7.7). However, if the training set is large enough and can well represent the diversity of the whole shape population of the target object, the representation residual error is small.

In practice, we do not want to use the full shape space learned from a training set. The first t deformation modes explain most amount of the shape variation and the remaining modes mainly correspond to noise in the data. Therefore, a shape m is often approximated as
This process is often called a projection into a shape subspace. The shape space is often truncated such that a certain amount α (e.g., α = 98 %) of variation is preserved. Therefore, t is the minimum integer such that
By varying the weights w _i in Eq. (7.9), a new shape can be generated. To make the generated shape similar to those in the training set, we often enforce a limit to the weights,

The covariance matrix in Eq. (7.4) has a dimension of 3J × 3J. In most of our applications, a shape is represented as a dense mesh to achieve a detailed segmentation of the object. For example, our whole heart model is represented as 3,005 points. The corresponding covariance matrix has a dimension of 9, 015 × 9, 015, which consumes a lot of memory (310 MB if the values are represented as 4-byte ‘double’ variables) and eigen-decomposition of such a huge matrix can be quite slow. Fortunately, the rank of the covariance matrix is bounded,
In practice, the number of training shapes, N, is much smaller than the number of points in a shape, therefore N < < 3J. For example, our heart chamber detectors are trained on no more than 500 volumes. That means most eigenvalues of are zero. To enforce statistical shape constraint, the deformation modes with small (or zero) eigenvalues are discarded. It is not necessary to calculate those eigenvectors at all.

In the following, we present two efficient methods to calculate the eigenvectors P _i of based on the fact that has a low rank. We can write the residual vectors of aligned training shapes d m _i into a matrix,
The dimension of matrix D is 3J × N. We then perform Singular Value Decomposition (SVD) of matrix D directly, without calculating the covariance matrix at all. Suppose matrix D can be decomposed as
where U is a 3J × 3J matrix; is a 3J × N rectangular matrix with nonzero eigenvalues at the diagonal; and V is an N × N matrix. Here, U and V are orthonormal, and . Then, can be represented as follows,
The above equation is an eigen-decomposition of . The largest N eigenvalues λ _i of is related to the eigenvalues σ _i of D as
The remaining eigenvalues are zero. Eigenvector P _i of is a left-eigenvector of D, P _i  = U _i . The SVD of D is much faster than eigen-decomposition of since matrix D has a lower dimension.

Alternatively, we can calculate another covariance matrix as
The dimension of matrix is N × N, which is much smaller than that of . Suppose the eigen-decomposition of is
That means
Left multiply both sides of the equation with D, we get
That means has the same eigenvalues as , (s _i  = λ _i ); and for an eigenvector Q _i of , if we left multiple it with D, it is an eigenvector of , (P _i  = DQ _i ).

Both alternative methods presented here are much more efficient to calculate the eigenvalues and eigenvectors of than working directly on . Among all three methods, the one based on the eigen-decomposition of the covariance matrix is the most efficient in our applications since has the smallest dimension of N × N when N < 3J.