Surface structure feature matching algorithm for cardiac motion estimation

Background Cardiac diseases represent the leading cause of sudden death worldwide. During the development of cardiac diseases, the left ventricle (LV) changes obviously in structure and function. LV motion estimation plays an important role for diagnosis and treatment of cardiac diseases. To estimate LV motion accurately for cine magnetic resonance (MR) cardiac images, we develop an algorithm by combining point set matching with surface structure features of myocardium. Methods The structure features of myocardial wall are described by estimating the normal directions of points locating on the myocardium contours using an approximation approach. The Gaussian mixture model (GMM) of structure features is used to represent LV structure feature distribution. A new cost function is defined to represent the differences between two Gaussian mixture models, which are the GMM of structure features and the GMM of positions of two point sets. To optimize the cost function, its gradient is derived to use the Quasi-Newton (QN). Furthermore, to resolve the dis-convergence issue of Quasi-Newton for high-dimensional parameter space, Stochastic Gradient Descent (SGD) is used and SGD gradient is derived. Finally, the new cost function is solved by optimization combining SGD with QN. With the closed form expression of gradient, this paper provided a computationally efficient registration algorithm. Results Three public datasets are employed to verify the performance of our algorithm, including cardiac MR image sequences acquired from 33 subjects, 14 inter-subject heart cases, and the data obtained in MICCAI 2009s 3D Segmentation Challenge for Clinical Applications. We compare our results with those of the other point set registration methods for LV motion estimation. The obtained results demonstrate that our algorithm shows inherent statistical robustness, due to the combination of SGD and Quasi-Newton optimization. Furthermore, our method is shown to outperform other point set matching methods in the registration accuracy. Conclusions We provide a novel effective algorithm for cardiac motion estimation by introducing LV surface structure feature to point set matching. A new cost function is defined to measure the discrepancy between GMMs of two point sets. The GMM of point positions and the GMM of surface structure descriptor are defined at the same time. Optimization by combining SGD and Quasi-Newton is performed to solve the cost function. We experimentally demonstrate that our algorithm shows improved registration accuracy, and is convergent when used in high-dimensional parameter space.


Background
Cardiovascular diseases (CVDs) are the leading cause of death in the developed world, as reported by the World Health Organization. Detecting the abnormalities in myocardial functions can assist the establishment of the early diagnosis of cardiomyopathy. The LV changes in structure and function during the development of cardiovascular diseases. Analysis of LV structure using imaging instrument is shown to be effective in reducing CVDs mortality and morbidity. Magnetic Resonance Imaging (MRI) is a state-of-the-art technique for the direct examination of changes in myocardial structure [1], with good spatial resolution and high signal to noise ratio. It allows the analysis of the structure alterations of LV and can be used to measure the functional change of LV.
Image registration technique can be employed to estimate LV motion for MR images, assisting the diagnosis of cardiac diseases [2,3]. For instance, the myocardial hypertrophy disease can be diagnosed by detecting parameters such as blood flow, ejection fraction, stroke volume, and so on [4]. Additionally, this technique allows the investigation of cardiac pathology as well.
A large number of methods for cardiac functional and motion modeling have been developed. At present, the methods of LV motion estimation can be classified into three groups: image intensity-based, geometrical and segmentation model-based, and point set matching-based.
Image intensity-based registration method optimizes a similarity function between images at different phases. A spatial transformation is used to compute the displacement of the myocardium [5,6]. Chandrashekara et al. [7] used normalized mutual information between the short-axis and long-axis images to recover the complete three-dimensional motion of the myocardium. Ebrahimi et al. [8] introduced and evaluated the performance of a non-rigid joint multi-level image registration and intensity correction algorithm, which integrated intensity change compensation and motion correction into a unified model. Furthermore, Oubel et al. [9] applied an information-theoretic metric to measure the similarity between frames of tMRI for heart motion estimation. Shi et al. [10] used a spatially weighted similarity measure between the untagged cine and the four-dimensional pseudo-anatomical MR image over time for myocardial motion estimation, while Bai et al. [11] combined the similarity metric between two images and the similarity between two probabilistic label maps of images to register atlases and segment cardiac MR datasets.
The geometrical and segmentation model-based methods segment the myocardial wall and use active contours or surfaces to model the geometrical and mechanical structure of heart. In this model, geometric features are tracked to derive the motion of the heart walls and the dense displacement field was extracted by a non-registration model [12]. Escalanteramírez et al. [13] estimated the motion of the heart based on the optical flow and image structure information that was extracted from the steered Hermite transform coefficients in sequential computed-tomography (CT) images. Papademetris et al. [12] segmented MR images and used a shape-tracking approach to establish correspondence of objects. Macan et al. [14] segmented the LV and extracted characteristic boundary points by investigating curvature distribution along the three-dimensional surface of cardiac wall. Points in two consecutive frames were matched by comparing curvature for LV motion estimation. Shi et al. [15] represented the LV surface shape using Delaunay triangulation and established correspondence between surface features to recover dense field motion from the tagged MR data.
The point set matching-based method is a type of feature-based registration method. Mathematically, point set matching can be represented as a problem to establish the point-to-point correspondence and the spatial transformation between two point sets at the same time. A cost function is defined to measure the discrepancy between two point sets based on a transformation function. The optimal transformation parameters are obtained by optimizing cost function. The iterative closest point (ICP) method [16] is the most commonly used point matching method, which estimates the global transformation between two point sets. ICP supposes a one-to-one correspondence between two point sets based on the nearest neighbor criterion, and estimate an affine transformation between two point sets. Algorithms were subsequently developed [17][18][19] to deal with elastic matching between two point sets based on the idea of ICP. Chui et al. [17] presented the point matching problem as a joint optimization problem over the parameterization of the elastic spatial mapping and the softassign for the correspondence. However, the transformation estimated by [17] was based on the correspondence relationship between a virtual point set and a real point set, instead of correspondence relationship between two real point sets. This model was extended by Lin et al. [20], who used free-form deformation model in robust point matching to analyze LV motion. The FFD model was based on arbitrarily-shaped lattices instead of parallelepipedically-shaped lattices.
Furthermore, various metrics for determining the alignment of point sets are proposed. A correlation-based point set matching method was proposed, where point sets were represented by kernel correlation, such as Gaussian Mixture Model (GMM) [19]. In [19], two GMMs of two point sets are defined and the discrepancy of these GMMs is defined to describe the alignment of two point sets. Myronenko et al. [18] modelled the point set distribution with GMM, and constrained the motion of the point set in the temporal direction for displacement estimation of 3D echo images of LV. Liu et al. [21] defined a novel, more accurate and meaningful equivalent distance to measure the position discrepancy between two point sets. Ravikumar et al. [22] proposed a probabilistic framework for group-wise rigid alignment of point-sets using a mixture of Student's t-distribution. Their method reduces alignment errors significantly in the presence of outliers. Du et al. [23] integrated a Gaussian probability model into the bounded-scaled registration problem to deal with the alignment of point sets with noise.
Here, we employ point set matching-based method to estimate motion of LV. Considering that existed point set matching methods are primarily based on the spatial relationship between two point sets, we introduce the structure information of LV to point set matching to improve LV registration accuracy at different time points. In this study, we develop a novel point set matching algorithm by considering the surface structure of LV.
The key contributions of this study are as follows: (1) The normal direction is computed as the surface structure feature to describe the structure of myocardial wall. Additionally, a cost function for the determination of the discrepancy of the GMMs of positions and GMMs of surface structures of two point sets is developed. (2) To resolve the dis-convergence problem of the optimization in high-dimensional parameter space, the Stochastic Gradient Descent (SGD) is combined with Quasi-Newton method in order to estimate optimal transformation parameters.

Overview of the methodology
Two cardiac slices taken at different times are considered a target image and a source image. For example, the slice obtained in the end-diastole represent a target image and the one obtained in end-systole represent the source image. Point sets are extracted from these images, and the point set in the target image is considered a scene set, while the other represent a model set. Since the slice numbers of two 3D cardiac images along long axis are different, we interpolate the point sets along this axis to obtain an equal number of slices of two corresponding 3D cardiac images. Afterwards, corresponding slices are considered source image and target image. The surface structure features of points located on myocardial wall are estimated approximately. A new cost function is defined as the combination of the GMM of spatial locations and the GMM of surface structures between two point sets, which is optimized by SGD and Quasi-Newton methods.

Point interpolation along the long axis
Cardiac MR spatial resolution is low along the long axis and relatively high along the short axis, as shown in Fig. 1. Moreover, slice numbers of two 3D cardiac MR images in  Fig. 2a respectively, which demonstrates that the slices in these two images do not matched each other along the long axis. To resolve this problem, points are interpolated along long axis to equalize the number of slices between two images. Next, point sets of the corresponding slices are registered.
To interpolate point sets along the long axis, we construct a section traversing the center of LV image and being parallel to the long axis, as shown in Fig. 2b. All the points of this section are interpolated. At first, the original point sets are interpolated by B-spline interpolation on the xy plane slice by slice to ensure that there are dense sampled points on the section. For a given section, all points located on the section are interpolated to make the slice numbers of two 3D point sets to be equal. Secondly, the section is rotated around the long axis to obtain a new section, and the point interpolation on the new section is performed again. Repeating the described procedure, we obtain a number of points along the long axis and make sure that the number of slices in two 3D LV images are equal. Figure 2c shows an example of an interpolated point set, where red points represent the interpolated points.

Surface structure description
The heart is composed of a muscular contractile organ (myocardium) surrounded by two layers of connective tissue, endocardium and epicardium. The LV and right ventricle (RV) are separated by endocardium and epicardium. The cardiac LV has a thick muscular wall, and its structure can be approximated by a curved surface. Here, we describe the curved surface of myocardial wall using structure features, described by the normal directions of points located on the myocardial wall. Fig. 3a shows the extracted points from a slice, and B and C are two points nearest to point A. As we know, a point with tangential direction − → BC exists based on the Lagrange mean value theorem. The normal direction of A can be approximated as the perpendicular line to − → BC when B and C are close to A enough. That implies, dense sampled points can ensure the accuracy of this approximation. An example of normal direction estimation is presented in Fig 1 , · · · , v n ] T of scene point set can be presented as For the model point set M 0 , we denote its mapped point set as M = (m 1 , m 2 , ..., m m ) under a spatial transformation. Following this, we will describe the surface structure description V M of M. Here, we employ the thin plate splines (TPS) to be the spatial transformation. Let Q = q 1 , q 2 , ..., q c be the control point sets with c points, q j = q jx , q jy , j = 1, · · · , c. The mapped model point m i = m ix , m iy , i = 1, · · · , m, is expressed by TPS as where φ(r) = −r 2 log(r 2 ) is the TPS basis function in 2D; m 0 i −q j is the Euclidean distance between m 0 i and q j ; a 0x , a 1x , a 2x are affine coefficients; w jx and w jy , j = 1, · · · , c, are elastic coefficients in x and y axes respectively.
The Eq. (2) can be rewritten as: Assuming m i and m j are the nearest two points of m k , the normal direction of m k can be approximated as Based on above analysis, it is known

Point set matching using GMMs
GMM is used for the point set matching to describe point distribution [18,19,24]. Point set matching based on GMM is a correlation-based approach, which represents point sets as probability densities. Gaussian kernel is commonly used to estimate the probability density. Bingjian et al. [19] minimized the discrepancy of two Gaussian mixtures to align two point sets. In [19], they describe the gaussian mixture density function as accumulation of k gaussian functions, is a gaussian function with mean vector μ i and variance σ i .
As we know, when the number of Gaussian model is large enough, almost any probability density can be well approximated by this model. The GMM of the point set M can be represented as a gaussian mixture density function as: In this model, all variances are same to σ ; each point m i in M is the mean of the ith gaussian function.
In the method presented by Bingjian et al. [19], only spatial information is employed for point set matching, instead of introducing additional information, such as the object structure. Surface structure description is a kind of structure features to represent the anatomical structure of LV myocardial wall. Following this, we improve point set matching accuracy by introducing the surface structure feature to the gaussian mixture density. Based on a previous study [19], we expand the point set matching model by introducing the discrepancy of the GMM of surface structure features.
Similar to the GMM of point positions, the GMM of surface structure features of V M is defined as gmm(v, V M ), where v and V M are similar to x and M in (6). We select L 2 distance to determine the similarity between two Gaussian mixtures of surface structure descriptions, then the discrepancy between two GMMs of surface structure Putting the discrepancy of all GMMs together, a new cost function of point set matching is defined as, The first term in Eq. (7) is the discrepancy of GMMs based on point positions, and the second term is the discrepancy of GMMs based on surface structure descriptions of point sets. λ 2 W T KW is a penalty term to regularize the transformation to be smooth, where the ijth entry of matrix K is K ij = φ( q i −q j ), i, j = 1, . . . , c. λ and β are coefficients to balance the terms in Eq. (7). For LV motion estimation, the model set M 0 is extracted from the endsystole (ES) phase and the scene set S is extracted from the end-diastole (ED) phase. The goal of LV motion estimation is to find the parameter θ of a spatial transformation by minimizing the cost function (7).
The Eq. (7) can be rewritten as: The closed-form expression between two gaussian mixtures can be easily derived as: based on the following formula: Following this, Eq. (8) can be formulated in detail, Removing irrelevant terms in Eq. (11), the final cost function J is: It is noteworthy that J is convex and is differentiable, which implies that gradient-based optimization techniques can be used to optimize J, such as the Quasi-Newton method.

Optimization using the Quasi-Newton method
Quasi-Newton method represents one of the most effective ways to solve the nonlinear optimization problem, with fast convergence speed and high accuracy. Here, we derive the gradient of Eq. (12) in detail. We seperate J as J = J d + J v , where J d is the sum of the first two items and λ 2 W T KW , and J v represents the other two items. Then, where G = ∂J ∂M is a m × 2 matrix. J v can be written as Obtaining the derivatives of J 1,v and J 2,v is simple: Therefore, the derivatives of J v can be obtained as while the derivatives of Eq. (12) can be obtained as Although Quasi-Newton method can solve the nonlinear optimization problem, it might be divergent in high dimensional parameter space. In other word, Quasi-Newton method cannot possibly find solutions when there are too many control points. In order to resolve this issue, the Stochastic Gradient Descent (SGD) algorithm is employed to optimize Eq. (12) robustly. Furthermore, considering that using SGD is not good at obtaining the precise solution, Quasi-Newton method is performed further to improve the accuracy of solution when SGD is convergent.

Optimization using Stochastic Gradient Descent
Quasi-Newton method uses the full training set to update parameters at each iteration, which tends to converge to local optima easily in high dimensional parameter space. SGD addresses this issues by following the negative gradient of the cost function using only a single or a few training examples [25]. By measuring gradient changes, it is easy to construct a model of the cost function to produce a superlinear convergence. SGD can follow the negative gradient of the cost function after being exposed to only a single or a few training examples, which can overcome computation cost and lead to fast convergence. Observing that our cost function defined above is differentiable, the SGD algorithm can simply compute the gradient of the cost function using only a single moving model point m i . Let f (θ, i) be the part of J, which is influenced by m i : The five terms in Eq. (21) are denoted as f 1 , f 2 , f 3 , f 4 , f 5 , respectively. Then, the derivatives of f (θ, i) is ∂θ . The derivative of each term is: Based on the SGD procedure, Algorithm 1 represents details of optimization using SGD. Update the learning rate η . 8: until an approximate minimum is obtained SGD algorithm is robust to obtain an approximate solution, but it is not good at finding an accurate solution.
To improve the solution accuracy further, Quasi-Newton method is employed for the optimization of Eq. (12), beginning at the optimal solution obtained by SGD.

Results and discussion
In order to confirm the improved performance of our method, three cardiac datasets are used. The first dataset included cardiac MR image sequences acquired from 33 subjects [26], the second set is composed of 14 intersubject heart slices [27], while the third set is the MIC-CAI 2009s 3D Cardiac Segmentation Challenge dataset [28]. In our study, the Quasi-Newton method and the combination of SGD and Quasi-Newton method are employed to demonstrate the improved performance of optimization of the cost function. The name following "SSD" in this article represents the optimization used in our method, and for example, "SSD_QN" denotes that the Quasi-Newton method is used. To compare the performance of our method with those of the previously developed methods, GMM [19] and CPD [18] methods are evaluated as well.

Cardiac MR image registration
Initially, we use a cardiac dataset [26] to evaluate the registration accuracy of our method. In this dataset, 33 cardiac MR image sequences are provided. Each image sequence consisted of 20 frames and the number of slices acquired along the long axis ranged from 8 to 15. The distance between slices ranged from 6 to 13 mm. The size of all image slices is 256 × 256 pixels with a pixel-spacing of 0.93-1.64 mm [26].
We evaluate the registration accuracy of our algorithm by registering cardiac images at the end-diastole phase and the end-systole phase. The point sets of LV at end-diastole and end-systole phases are registered. Here, these points are provided by experts [26] in order to eliminate the effects of cardiac segmentation. The point set at end-diastole phase is the scene point set, and the one at end-systole phase is the model point set. The end-systole points are interpolated along the long axis to make the slice numbers equal to that in the end-diastole points. Afterwards, these two point sets are matched using our algorithm slice by slice. In Fig. 4, the interpolated model points and scene points are presented, and they are extracted from the slices in end-systole and end-diastole phases respectively.
Average Perpendicular Distance (APD) between the mapped model point set and the scene point set is computed to evaluate our algorithm quantitatively. The APD Fig. 4 An example of a scene point set and b model point set extracted from images at end-diastole and end-systole phases respectively represent the average distance between two point sets, as shown in Fig. 5.
We use 10 2 control points to estimate the transformation function between end-systole and end-diastole slices. Table 1 lists the APD results of 33 subjects using three different algorithms. For each subject, the APD is averaged registered results of all corresponding slices along the long axis. Noted that the APD results of SSD_QN are smaller than that of GMM and CPD for most subjects. It implies that the registration accuracy can be improved by introducing the surface structure features to the point set matching.
Analyses conducted using the Quasi-Newton method demonstrate the issue of divergence that appears with large number of control points. We increase the number of control points from 20 2 to 40 2 . The APD results are listed in Table 2, and we show that the results obtained by applying the SSD_QN and GMM methods diverge when too many control points used, while CPD is shown to be robust in convergence. For the convergent SSD_QN, it outperforms GMM and CPD in respect of APD for most cases.
To demonstrate the improved performance of SSD_S -GD_QN further, 50 2 control points are used to register the slices between end-systole and end-diastole in all samples. In Fig. 6, APD results obtained for 33 subjects using CPD and SSD_SGD_QN with 50 2 control points are presented. Since GMM is divergent, the APD results obtained by this method are not provided in Fig. 6. For the majority of cases, SSD_SGD_QN outperforms CPD, which demonstrates the stability of SSD_SGD_QN in high dimensional parameter space.

Registration of the inter-subject LV
Furthermore, we use a cardiac dataset provided by the Danish Research Centre for Magnetic Resonance (DRCMR) [27], comprising 14 gray scale 256 × 256 images. All images used are short-axis, end-diastolic cardiac MR images, acquired using a whole-body MR unit  (Siemens Impact) operating at 1.0 Tesla. The endocardial and epicardial contours of the LV are manually annotated by experts [27]. Case 1 represent a sence image, while all other cases are model images. The contour marked in other cases are mapped to the image of case 1 using SSD_QN, GMM, and CPD. Afterward, the mapped contour is compared with a contour marked by an expert, to evaluate the registration accuracy. We use the APD to evaluate the performance of three methods, since it can be used to evaluate the registration accuracy of three methods. We locate the region of interest (ROI) in heart in the original image, and employ the template matching [29] algorithm (TMA) to extract epicardium and endocardium profiles by constructing many typical LV templates.
The candidate template is generated by a particle, and the optimal particle is obtained by matching the target and the candidate templates. Following this, the point sets of epicardium and endocardium are extracted around the candidate template contour, as shown in Fig. 7.
In this experiment, 10 2 control points are employed to ensure the convergence of SSD_QN. The APD results are listed in Table 3. SSD_QN is shown to outperform GMM Fig. 6 The APD results using SSD_SGD_QN and CPD respectively and CPD in all cases except case 13. In this case, many noisy points are observed in the extracted myocardial contour, which leads to the registration errors.
To visualize the registration results further, we present the registration results obtained by three methods for case 2 in Fig. 8. The ground truth is marked by red points, and the mapped contours estimated by three methods are marked by green points. The outlines of endocardium and epicardium, obtained using SSD_QN, are shown to be close to the ground truth. Moreover, these outlines are demonstrated to be smoother than those obtained using GMM and CPD, which demonstrates the advantage of surface structure feature for the description of the circular structure of myocardium.
Furthermore, we use SSD_SGD_QN to demonstrate the optimization performance in point registration. In Fig. 9, the APD results obtained using SSD_SGD_QN, GMM, and CPD with 20 2 control points are presented. Although GMM is shown to converge in this experiment, its registration accuracy is not satisfactory. CPD is shown to be converged, but the APD results obtained using this method are shown to be larger than those determined using SSD_SGD_QN and GMM.

Registration of cine MR cardiac images
In order to confirm the superior performance of our method in the registration of cine MR cardiac images, we analyze 15 cardiac cine MR validation datasets from the    This MR cardiac image dataset is used for MICCAI 2009s 3D Segmentation Challenge for Clinical Applications. Cardiac image segmentation is related to the image registration, and the segmentation results of one slice obtained in a single time point can be propagated to other time points using deformable registration, which takes advantage of the strong temporal correlation between phases. Here, we analyze the transformation between slices at the enddiastolic and endsystolic phases. Registrations are performed from enddiastole to endsystole and vice versa. The LV contours extracted from the endsystolic phase can be mapped to the enddiastolic phase using the estimated transformation, representing the segmented results of LV at enddiastole, and vice versa.
For this analysis, we eliminate the surrounding organs and locate the area containing LV. Considering that only the contour of endocardium is provided in MICCAI 2009s dataset, we segment the endocardium and extract the boundary points of endocardium by edge detection. We ignore the effects of the papillary muscles of endocardial wall, as they are minor. Since the outline of endocardium is irregular, we smooth out the border by using Savitzky-Golay polynomial filter [30]. The original cardiac slice and an automatically extracted endocardial outline are presented in Fig. 10.
Two principal evaluation standards in the MICCAI 2009s 3D Segmentation Challenge for Clinical Applications are the APD, which measures the average distance  error between two point sets, and the Dice Metric (DM), that shows the contour overlap proportion between the mapped contour and the target contour. DM = 2S 3 S 1 +S 2 , where S 1 and S 2 are endocardial surface areas obtained  manually by experts and by automatic methods, while S 3 represents the overlap area between S 1 and S 2 . A higher DM values indicates better registration results. The forward registration (endsystole to enddiastole) and the reverse registration (enddiastole to endsystole) are performed respectively. All these DM values represent the average of registration results in two directions. We used 10 2 control points to estimate the transformation function between two heart slices, and 15 clinical cine MR cases with four patient categories (heart failure, with (HF-I) or without (HF-NI) ischemia, hypertrophy (HYP), and normal (N)) are used to investigate the performance of SSD_QN, GMM, and CPD.
In Tables 4 and 5, the APD and DM results obtained using three described methods are presented. The APD results obtained by SSD_QN are shown to be smaller than those obtained by GMM and CPD in most cases, except for one case. In the case of SC-HYP-37, the extracted points are not able to outline the endocardium wall, which led to the deviation of the normal directions of points. Consequently, registration accuracy using SSD_QN decreased, which indicates that the initial segmentation results affect the precision of registration.
For a detailed illustration of the registration results, the mapped contour and the contour marked by experts using the data from subject SC-HF-I-05 are presented in Fig. 11. The first line shows the slices in the enddiastolic phase and the mapped endsystolic contours using SSD_QN, GMM, and CPD. The second line indicates the slices in the endsystolic phase and the mapped enddiastolic contours using three methods. The ground truths are marked by green points, and the mapped contours are denoted by blue points. Short red lines illustrate the minimum distance between the points marked by the expert and the mapped points. The shortening of the red lines indicate a higher registration accuracy. The APD obtained by SSD_QN is shown to be superior to that obtained by GMM and CPD for forward registration. However, the APD of the reverse registration obtained using SSD_QN is larger than that obtained using GMM and CPD, because the parameters used in our experiments are suitable for forward, and not reverse registration.
The APD obtained by SSD_SGD_QN, GMM, and CPD with 20 2 control points are compared as well in Fig. 12. Similar to the previously presented results, GMM and CPD are converged with these parameters as well, and the APD obtained using SSD_SGD_QN is shown to be close to those obtained using GMM and CPD.  Finally, the displacement vector graphs estimated by our algorithm SSD_QN with 10 2 control points are illustrated in Fig. 13. Displacement vectors of LV from endsystole to enddiastole in three slices are shown, where red points and blue points represent the model and the mapped model points, and green arrows illustrate the displacement vectors of these points. The estimated LV motion is shown to coincide generally with the real motion.

Conclusion
LV motion estimation is important in quantitative assessment of myocardial function and dynamic behavior of human heart, which is invaluable in the diagnosis of cardiac diseases. In this paper, a novel point set matching algorithm is proposed to estimate LV motion. The main contribution of the proposed algorithm is introducing the structure feature of LV to point set matching. The surface structure features of LV is described using normal directions, and the GMM of surface structure features is defined. By measuring the discrepancy of all GMMs of two point sets, a new cost function of point set matching is constructed. SGD and Quasi-Newton method are combined to optimize the cost function. Performance of our algorithm is verified using three cardiac image datasets. The obtained results demonstrate that when small amount of control points used, our algorithm with Quasi-Newton optimization outperforms GMM and CPD in LV motion estimation. When too many control points used, our algorithm with the combination of SGD and QN optimization is more robust than GMM and CPD. The evaluation performed using MICCAI 2009s 3D Segmentation Challenge for Clinical Applications dataset demonstrate that the applicability of our motion estimation method remains the same when analyzing different cardiac diseases. The method we develop and present here could be applied to clinicallyobtained data, to demonstrate its applicability in a clinical environment.