A Geometrical Approach for Automatic Shape Restoration of the Left Ventricle

This paper describes an automatic algorithm that uses a geometry-driven optimization approach to restore the shape of three-dimensional (3D) left ventricular (LV) models created from magnetic resonance imaging (MRI) data. The basic premise is to restore the LV shape such that the LV epicardial surface is smooth after the restoration and that the general shape characteristic of the LV is not altered. The Maximum Principle Curvature () and the Minimum Principle Curvature () of the LV epicardial surface are used to construct a shape-based optimization objective function to restore the shape of a motion-affected LV via a dual-resolution semi-rigid deformation process and a free-form geometric deformation process. A limited memory quasi-Newton algorithm, L-BFGS-B, is then used to solve the optimization problem. The goal of the optimization is to achieve a smooth epicardial shape by iterative in-plane and through-plane translation of vertices in the LV model. We tested our algorithm on 30 sets of LV models with simulated motion artifact generated from a very smooth patient sample, and 20 in vivo patient-specific models which contain significant motion artifacts. In the 30 simulated samples, the Hausdorff distances with respect to the Ground Truth are significantly reduced after restoration, signifying that the algorithm can restore geometrical accuracy of motion-affected LV models. In the 20 in vivo patient-specific models, the results show that our method is able to restore the shape of LV models without altering the general shape of the model. The magnitudes of in-plane translations are also consistent with existing registration techniques and experimental findings.


Introduction
Breath-hold cine Magnetic Resonance Imaging (MRI) is an advanced imaging technique for cardiac morphological and functional assessment in clinical practice. While conventional methods of evaluation are based on MRI images, several recent methods [1], [2], [3] have been developed to utilize threedimensional (3D) models reconstructed from the MRI data. A 3D model of the LV provides a more comprehensive and accurate description of the ventricular shape and function as properties can be extracted from a combination of in-plane and out-of-plane information, as compared to analysis done on 2D methods. This has a direct impact on the evaluation of LV chamber properties such as its chamber volume, local wall curvature, myocardial wall thickness and wall stress. In our prior study, we have compared the results obtained by conventional 2D methods and 3D methods [2] and found that 3D-based quantification of regional wall stress provides more precise evaluation of cardiac mechanics.
However, factors such as respiration and patient movement contribute to misalignments in the MRI data which results in inaccuracies in the 3D models. MRI data acquired over different breath-hold positions also induce errors in the reconstructed models. [4] aims to correct for misaligned cardiac anatomy, caused by differing breath-hold positions, in multi-slice short-axis (SA) images, by rigidly registering stacks of two slices to a high-resolution 3D MR axial cardiac volume. Existing registration methods also include multi-modal dynamic cardiac image registration [5], external skin marker-based techniques [6], landmark-based techniques [7] and thorax surface-based techniques [8]. A number of cardiac image registration methods are reviewed in [9]. They are categorized into the geometric image feature approach and voxel similarity measure approach. Recently, some post-processing methods have also been proposed. The method proposed by Elen et al. [10] demonstrates the use of constrained optimization, on the assumption that the similarity of gray values at the intersection lines of different slices is higher when the relative positioning of the slices is correct than when the slices are misaligned.
The main dilemma of using an image registration approach to restore the shape of the LV is that errors induced by motion are already embedded in the images. While multi-view image registration techniques could potentially reduce such errors, these methods are essentially using data which contain error for selfcorrection. In our work, we make use of morphological knowledge of the LV to drive the shape restoration. Instead of using imagebased parameters, such as gray values, our LV shape restoration method is based on geometrical consideration. The basic premise is that the LV epicardial surface must be smooth after the restoration, which is a reasonable assumption based on observed morphology, that the myocardial surface tension or force of the LV wall forces the hearts to minimize the wall surface area, creating a smooth epicardial surface. In addition, the general shape of the LV cannot be lost in the process such as its skewness and asymmetrical configuration. The Maximum Principle Curvature (k 1 ) and Minimum Principle Curvature k 2 of the LV epicardial surface are used as the geometric measure to quantitate the local shape characteristics. Figure 1(a) shows a reconstructed 3D LV mesh model from MRI data containing motion artifacts. We aim to achieve a smooth LV mesh as shown in Figure 1(b) after shape restoration. This is achieved by a dual-resolution semi-rigid deformation, followed by a free-form geometric deformation process. The semirigid deformation process involves shifting of the myocardium contours by translating each slice in the in-plane direction (xyplane), while the free-form deformation process involves translating each individual vertex of the LV model in the x -, y -and zdirections. We formulate a smoothness objective function based on k 1 and k 2 , and solve the problem using a limited-memory quasi-Newton optimization algorithm, L-BFGS-B [11]. The L-BFGS-B algorithm is an adaptation of the BFGS algorithm with limited matrix update and it is adept at solving multivariate nonlinear bound constrained optimization problems. This paper is organized as follows: First, we provide detailed explanation of the methodology of our LV shape restoration algorithm. Next we describe the experiments done on the 30 simulated samples and the 20 in vivo patient-specific models to test the performance of the algorithm, followed by a discussion on the implications of the experimental results, and finally conclude the paper.

Methods
The input to our algorithm is an initial 3D LV mesh model M~fC k jk~1,2,3:::,Ng reconstructed from a set of contours {C} representing the myocardial borders delineated from the SAplanes, such that N is the total number of contours. Each contour C k~f V k,i ji~1,2,3,:::m k g consists of a set of closed connected vertices {V} where m k is the total number of vertices in the k -th contour. The convention used is such that the SA-slices are parallel to the xy-plane; the contours are arranged from the apex to basal region in increasing z-values; and the vertices in the contours are cyclic (i.e., V k,m k is connected to V k,1 , since they represent a closed connected curve).
Shape Quantification Using Principle Curvatures k 1 and k 2 In this section, we briefly outline the formulation of k 1 and k 2 . In order to interrogate the geometrical properties of the LV epicardial surface mesh, we use a quadric fitting method to approximate the underlying geometry at every vertex of the mesh. A quadric surface S in 3D space can be expressed in the parametric form where u and v are the surface parameters and fa,b,c,d,eg are the quadric coefficients. To fit S at a vertex p , we select a neighborhood around p which represents the region over which S is to be fitted. The extent of this neighborhood is quantified by a n -ring value. The quadric coefficients of S are then obtained by solving a system of linear equations associated with the n -ring neighborhood using a least square method [12]. The surface S approximates the local geometry in the vicinity of a point p on the 3D mesh model. In differential geometry, the curvature of a surface S(u,v) at a point p(u,v) is evaluated with respect to a normal section. This is done by constructing a plane p such that it passes through the unit surface normaln n and unit tangent vector The intersection of p with S results in a curve called the normal section. The normal curvature k(_ v v) can be evaluated by where G~½ S u :S u S u :S v S u :S v S v :S v ~½ E F F G and D~½ S uu :n n S uv :n n S uv :n n S vv :n n ~½ L M M N are the first and second fundamental matrices of the surface, respectively. The unit surface normal can be calculated by In terms of the quadric coefficients, the equations to calculate k 1 and k 2 are where A~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d 2 ze 2 z1 p and B~azae 2 zczcd 2 zbde . The value of n -ring used in the quadric fitting affects the value of k 1 and k 2 because it determines how sensitive the method is to the effect of geometrical variation. With a bigger n -ring value, shape of the surface over a larger extent is interrogated. This takes into account the general variation of the shape, ignoring the high frequency variation in the geometry. With a smaller n -ring value, the shape of the surface over a localized region is inspected. This captures the inter-slice variations in shape.

Shape-based Optimization Objective Function
The basic premise of the shape restoration is the assumption that the LV epicardial surface is smooth. From a shape characterization perspective, this implies that it should have minimum local concavity. The objective is then to find the optimal modifications in the LV model such that its total concavity is at a global minimum. Computationally, we can calculate the maximum and minimum principle curvatures (k 1 and k 2 , respectively) for every point on the LV epicardial surface mesh to assess the amount of concavity or convexity of the surface. When k 1 or k 2 is negative, it implies that the surface at which a point lies on is concaved. Therefore, to minimize concavity, the objective function F only takes into account the summation of k 1 and k 2 values of all points with negative k 1 or negative k 2 , that is, where k 1,i is the maximum principal curvature and k 2,i is the minimum principal curvature at vertex i and m is the total number of vertices in the epicardial surface mesh. This non-linear objective function can be minimized using the L-BFGS-B algorithm [11] which is adept at solving multivariate nonlinear bound constrained optimization problems. It is based on the gradient projection method and uses a limited-memory BFGS matrix to approximate the Hessian of the objective function. The algorithm does not store the results from all iterations but only a user-specified subset. Its advantage is that it makes simple approximations of the Hessian matrices which are still good enough for a fast rate linear convergence and requires minimal  storage [11]. This will result in geometrical kinks being smoothed out but not at the expense of creating more kinks in other locations. The result is an improvement in the overall smoothness of the LV shape.
In order to maintain the overall shape characteristic of the LV model, such as its skewness and asymmetrical configuration, the optimization is performed in a dual-resolution semi-rigid geometrical deformation whereby the LV model is modified from a global and regional consideration, followed by a freeform geometrical deformation process for localized optimization.

Dual-resolution Semi-rigid Geometric Deformation
As the 3D LV models are reconstructed from contours of the myocardial borders, the first stage of the shape restoration process works by progressively translating these contours in the plane of their respective SA-slice. Figure 2 illustrates the shifting of the contour on a SA-slice in the in-plane direction. One can view this as a semi-rigid mesh modification since we are keeping the shape of the contours constant while shifting their position. For every slice, a centroid (X C,i ,Y C,i ) is calculated by averaging the x -and y -coordinates of all the points from that particular slice, where i is the slice index. The optimization algorithm, L-BFGS-B, will solve for the optimal (X 0 C ,Y 0 C ) for all the slices to satisfy the objective function in Equation (5).
To set up the optimization problem, we can write Equation (5) as F (x) with n variables, such that x contains the centroid coordinates (X C ,Y C ) of the contours on the SA-slices, i.e x~x 1~XC,1 To retain the skewness and asymmetry of the dataset, we fix the centroid coordinates of the top and bottom SA slices as boundary conditions, i.e.
. Hence, the number of variables in the optimization problem is n~2|(N{2) . Each of the variables x i in F(x) is subjected to the bounded-constraints lb i ƒx i ƒub i k~1,2,3,:::,n ð6Þ where lb i and ub i are the lower and upper bounds of x i , respectively. In this work, the variables are constrained to translate within a bound of + 20 mm. This value is consistent with what was observed experimentally [13] (maximum displacement recorded upon inhalation is 23.5 mm). The constraint on the translation in the x-direction is where X 0 C,k is the solution and X C,k s is the initial x-coordinate of the centroid position of the k -th contour. Similarly, the constraint on the translation in the y -direction is where Y 0 C,k is the solution and Y C,k is the initial y -coordinate of the centroid position of the k-th contour. In addition, the gradient g k associated with each variable x k must also be defined such that Since F (x) is in a non-analytical form, we need to approximate g i using finite differences. In this work, we use the forward difference method to approximate the gradient, i.e., where Dx k is a small increment in x k . In order to retain the general variation of the LV shape, we perform the optimization in dual stages -a global stage and a regional stage. As the objective function in the optimization is computed using a n -ring setting, we employ that to determine the nature of the shape characterization of the LV. In general, a large n -ring will result in shape characterization from a global perspective while a smaller n-ring will result in a regional/local shape characterization.
In the global stage, the value of n -ring selected is adaptive to the sampling resolution of the LV model. For our application, it is half the number of MRI stacks constituting the whole LV. However, the n -ring selected cannot be too huge as computation time increases with the value of n -ring. Therefore, we enforce that In the example shown in Figure 3(a), it shows an original mesh with motion artifact, where the n-ring selected is 4 due to it having 8 MRI stacks. When n -ring = 4, k 1 and k 2 are calculated by taking into account points from 4 layers above and below the current SA slice, and 4 points to the right and left of the point of interest. All the slices will shift to minimize the objection function in Equation (5). Next, to further minimize surface concavity over a smaller region of consideration, the intermediate mesh (updated with the previously obtained solution using n -ring = 4) is subjected to a second pass of optimization using a fixed n -ring = 2. This second pass is essential to further minimize the concavity over a localized region. The results from setting n-ring value = 4 and then n-ring = 2 are shown in Figure 3(b), intermediate mesh after optimization using n -ring = 4 and Figure 3(c), final mesh after optimization using n -ring = 2. As observed, the LV smoothness was restored.

Free-form Geometric Deformation
In the previous section, we discussed a dual-resolution semirigid deformation process, whereby only translational displacement of the SA-slices is considered. However, in actual fact, motion artifacts are also generated by motions in and out of the SA-planes. Therefore, in the next stage, a local free-form geometric deformation is performed. Here, the output from the semi-rigid geometric deformation process forms the input to the free-form deformation process. The shape restoration is done by progressively translating each individual vertex of the LV model in the x -, y -and z -directions. Using the same assumption that the LV epicardial surface is smooth, our objective is to find the optimal translations in the x -, y -and z-directions for each individual vertex such that the total concavity of the whole LV is at its global minimum. The same objective function F in Equation (5)     x i in F (x) is subjected to the bounded-constraints. The constraints in the x -and y -directions are determined by the distance between the vertex and its immediate neighboring vertices on the same slice/contour: where X 0 V k,i is the solution and X V k,i is the initial x-coordinate of the position of vertex V k,i . Similarly, the constraint on the translation in the y-direction is The constraint on the translation in the z -direction is determined by the inter-slice distance such that where j(Z V k,i {Z V k{1,i )j is the original distance between contour C k and its lower adjacent contour C k{1 ; j(Z V k,i {Z V kz1,i )j is the original distance between contour C k and its upper adjacent contour C kz1 ; and h is the angle between the vertex normal and the z-axis.
In Figure 4(a), we illustrate the impact of h on the volume of the LV mesh with respect to the vertical shift of the vertices. Two vertices V 1 and V 2 on the same contour C k are such that h 1 .h 2 . Given the same allowance of vertical shift such that their new positions become V 0 1 and V 2 , we observed that the deviation of V 0 2 from the original surface of the mesh is greater than V 0 1 , indicating that when h is smaller, the resulting volume change due to the vertical vertex shift is larger. Therefore, the sin h function in Equation (14) is used to constrain the vertices such that if there is greater deviation between the vertex normal from the SA-plane (i.e., h is small), the allowable vertical translation in the z-direction will be less. This will prevent unduly large change in the volume of the restored mesh. Figure 4(b) illustrates an LV epicardial surface modified by the free-form deformation process. Figure 5 shows the flow chart of the whole restoration process.

Results and Discussion
In vivo Patient-specific Models In this section, we tested our algorithm on 20 patient-specific 3D LV models reconstructed from MRI data containing motion artifacts. This study was approved by the SingHealth Centralised Institutional Review Board (CIRB No: 2009/705/C) for Human Research. All enrolled participants gave written informed consent. The MRI scan was performed using breath-held steady-state free precession technique on a 1.5T Siemens scanner (Avanto, Siemens Medical Solutions, Erlangen). TrueFISP (fast imaging with steadystate precession) MR pulse sequence with segmented k-space and retrospective electrocardiographic gating were used to acquire 2D cine images of the LV in the long-axis (LA) plane, as well as a parallel stack of 2D cine images of the LV in the SA plane, from the LV base to apex (8 mm inter-slice thickness, no inter-slice gap). Each slice was acquired in a single breath hold, with 25 temporal phases per heart cycle. The epicardial borders of contiguous SAslices were manually delineated by an experienced cardiologist using commercially-available software CMRtools (Cardiovascular Imaging Solution, UK). Both SA-and LA-views were utilized to carry out 3D LV reconstruction at the end-diastole phase. Figure 6 plots the absolute values of k 1 and k 2 before and after restoration of one of the patient sample. Quantitatively, the absolute values of k 1 and k 2 were reduced considerably after the shape restoration. Visually, we observed that the asymmetry of the LV geometry was preserved while the geometrical kinks on the surface were significantly reduced. Figure 7 shows that the shape of 2 other patient samples became smoother and retained their asymmetry after the whole dual and free form geometric deformation process. The results indicated that average contour displacements in the SA-planes were 1.36 mm and 1.30 mm, with a maximum translation magnitude of 8.68 mm and 8.65 mm in the xand y-directions respectively. The vertices are further displaced in the free-form deformation by an average of 0.10 mm, 0.11 mm and 0.07 mm, with a maximum translation magnitude of 1.73 mm, 1.95 mm and 2.93 mm in the x-, yand z-directions, respectively. As it is difficult to obtain the corresponding ground truth for each of the 20 in vivo patient-specific model, we assess the performance of our restoration algorithm by comparing the mean contour displacement values of our method with those of existing image registration techniques [14], as shown in Table 1, whereby the mean translations in the x-direction are 3 mm [15], [16] and y-direction is 3.25 mm [17] and z-direction is 1.5 mm [18]. The maximum experimental translation observed in [13] is 23.5 mm. Our results are observed to lie within the range of existing literatures.

Models with Simulated Motion
To further validate the applicability of our method, we study the performance of our method on a set of LV data with simulated motion. This is done by selecting a gold standard to act as a reference or ground truth for comparison. Ten healthy volunteers underwent breath hold cine-MRI scan. The sample with the smoothest epicardial surface is selected as the gold standard. To do this, the respective 3D models of the LV were reconstructed for each of the 10 volunteers and we subject these 3D models to a computational algorithm to assess the quantitative value of the local curvature. In addition, we apply our shape restoration method to these 3D models. The model which yields the minimum concavity and has the least modification in the restoration process is selected as the Ground Truth. This Ground Truth model is then used as a reference to investigate our algorithm's ability in performing LV shape restoration. We proceed to construct a set of 30 random samples from the ground truth model by translating three randomly selected but consecutive slices to simulate motion artifacts. They are translated within the range of 620 mm. These 30 randomized samples are then restored using our proposed algorithm and compared against the ground truth model for validation. To evaluate the accuracy of the restoration, we make use of The Metro Software [19] a tool designed to evaluate the difference between two triangular meshes. Metro adopts an approximated approach based on surface sampling and point-to-surface distance computation. The Hausdorff distance evaluated by the Metro Software measures the dissimilarity between two shapes, and it is used to quantify how different the LV mesh is before and after restoration, as compared to the ground truth. From Table 2, we can see that the Hausdorff distance after restoration w.r.t the  ground truth model is much smaller than that before restoration, except for Sample 1. This indicates that the restored meshes are more similar to the ground truth after the restoration process. In Sample 1, it is noticed that the Hausdorff distance w.r.t the ground truth model before restoration is just 1.03 mm. This implies that the randomly generated sample did not differ much from the ground truth model. The huge percentage difference of 21% is due to this small Hausdorff distance w.r.t the ground truth model before restoration. In addition to the analysis using Hausdorff distance, we evaluated the restoration using curvedness comparison. We have previously shown that curvedness is an important shape characterization measure for LV models and that it is associated with important LV functions [2]. Table 3 shows the curvedness results of the ground truth and the mean curvedness results of the 30 samples before and after restoration over 16 regions. Table 4 shows the absolute difference in curvedness with respect to the ground truth over 16 regions. The absolute difference in value and percentage of the curvedness with respect to the ground truth after restoration improved significantly for each region, which implies that curvedness value is closer to the ground truth after restoration.

Limitations and Challenges
As mentioned in the above section, the limitation of this study is in the validation of the restoration results of the 20 in vivo patientspecific models. While we do not have the ground truth for each of the patient-specific model, we compared the results of the mean contour displacements and observed that the values lie within the range reported in existing literatures. In addition, we attempted to account for this limitation by performing an experiment with simulated motion and observed significant improvement in terms of Hausdorff distance and curvedness.
We would also like to highlight that the derivation of the optimization function that drives the shape restoration is based on the observed morphology that the surface tension of the LV wall physically forces the heart to minimize the wall surface area, creating a smooth epicardial surface. However there are questions if this assumption still holds true for severe pathological conditions such as acute myocardial infarction. In a study of ventricular shape after myocardial infarction by Mitchell et al. [20], global LV geometry was assessed on cine angiography by means of a sphericity index. This index is the ratio of LV volume and a sphere with similar circumference. After myocardial infarction, the heart assumes a more spherical conformation (i.e., sphericity index increases), implying an overall diminution of surface concavity. This argument is also supported by calculating the total ventricular force or tension (T) from T~s x S, wherein s is the ventricular wall stress, and S is the surface area. For any given pressure and wall stress s, S being smallest would offer an optimal energy solution for moving blood. From clinical observation, the ventricles of patients with pathological condition remodel towards a spherical configuration (i.e., for any given ventricular volume, the surface area S is the smallest when the ventricular geometry is spherical), implying a ''smoother'' epicardial surface.

Conclusion
In this paper, we presented an automatic algorithm to restore the shape of a 3D LV mesh model using a geometry-driven optimization approach. The method used an analytical surface fitting method to approximate the geometry of the LV mesh and computed the minimum principal curvature as a quantification of the surface smoothness. Next, a limited memory quasi-Newton algorithm, L-BFGS-B, was used to correct the positions of all SAslices to achieve an optimal shape with minimal concavity. To retain the overall shape of the LV mesh, such as its asymmetrical configuration, we performed optimization using n-ring to be half the number of contours representing the myocardial borders. Next, to achieve localized smoothing, we performed a second pass of optimization using n-ring = 2. 30 sets of simulated data and 20 in vivo LV epicardial datasets at end-diastole are used as inputs to investigate the performance of our shape restoration algorithm. The results showed that there were significant improvements in the smoothness of the LV mesh both visually and quantitatively (in terms of the magnitude of maximum and minimum principal curvatures). Also, our algorithm was successful in preserving the overall shape of the LV mesh without over smoothing. Our results are also consistent with results of existing image registration techniques. Another 30 samples of data sets generated from a smooth patient set (i.e. Ground Truth), were restored and the difference in Hausdorff distance and curvedness values w.r.t the Ground Truth were smaller after restoration.