 Research
 Open Access
 Published:
Surface structure feature matching algorithm for cardiac motion estimation
BMC Medical Informatics and Decision Making volume 17, Article number: 172 (2017)
Abstract
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 QuasiNewton (QN). Furthermore, to resolve the disconvergence issue of QuasiNewton for highdimensional 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 intersubject 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 QuasiNewton 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 QuasiNewton is performed to solve the cost function. We experimentally demonstrate that our algorithm shows improved registration accuracy, and is convergent when used in highdimensional 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 stateoftheart 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 intensitybased, geometrical and segmentation modelbased, and point set matchingbased.
Image intensitybased 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 shortaxis and longaxis images to recover the complete threedimensional motion of the myocardium. Ebrahimi et al. [8] introduced and evaluated the performance of a nonrigid joint multilevel 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 informationtheoretic 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 fourdimensional pseudoanatomical 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 modelbased 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 nonregistration 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 computedtomography (CT) images. Papademetris et al. [12] segmented MR images and used a shapetracking approach to establish correspondence of objects. Macan et al. [14] segmented the LV and extracted characteristic boundary points by investigating curvature distribution along the threedimensional 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 matchingbased method is a type of featurebased registration method. Mathematically, point set matching can be represented as a problem to establish the pointtopoint 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 onetoone 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–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 freeform deformation model in robust point matching to analyze LV motion. The FFD model was based on arbitrarilyshaped lattices instead of parallelepipedicallyshaped lattices.
Furthermore, various metrics for determining the alignment of point sets are proposed. A correlationbased 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 groupwise rigid alignment of pointsets using a mixture of Student’s tdistribution. Their method reduces alignment errors significantly in the presence of outliers. Du et al. [23] integrated a Gaussian probability model into the boundedscaled registration problem to deal with the alignment of point sets with noise.
Here, we employ point set matchingbased 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 disconvergence problem of the optimization in highdimensional parameter space, the Stochastic Gradient Descent (SGD) is combined with QuasiNewton method in order to estimate optimal transformation parameters.
Methods
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 enddiastole represent a target image and the one obtained in endsystole 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 QuasiNewton 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 different phases differ. Two point sets in endsystole and enddiastole phases are illustrated slice by slice in Fig. 2 a 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. 2 b. All the points of this section are interpolated. At first, the original point sets are interpolated by Bspline 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 2 c 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. 3 a 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 \(\overrightarrow {BC}\) exists based on the Lagrange mean value theorem. The normal direction of A can be approximated as the perpendicular line to \(\overrightarrow {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. 3 b.
Let S=(s _{1},s _{2},...s _{ n })^{T} be the scene point set with n points, and \(M_{0}=\left (m_{1}^{0}, m_{2}^{0},...m_{m}^{0}\right)^{T}\) be the model point set with m points, where s _{ i },i=1,2,...,n and \(m_{j}^{0},j=1,2,...,m\) are two dimensional points. Supposing s _{ i } and s _{ j } are the nearest two points of s _{ k }, the normal direction vector of s _{ k } can be represented as \(\scriptsize v_{k}=H_{k}\begin {bmatrix}0 & 1\\ 1 & 0\end {bmatrix}\), where H _{ k }=s _{ j }−s _{ i }. Then, the surface structure description V _{ S }=[v _{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 follows:
where ϕ(r) = −r ^{2} l o g(r ^{2}) is the TPS basis function in 2D; \(\m_{i}^{0}q_{j}\\) is the Euclidean distance between \(m_{i}^{0}\) 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:
where [1M _{0}] is a m×3 matrix; A is a 3×2 matrix; W is a c×2 matrix. Φ is n×c matrix with \(\phi _{ij}\,=\,\phi \left (\left \left m_{i}^{0}q_{j}\right \right \right)\). Let N be the left null space of [1Q], which is a c×(c−3) matrix. A new (c−3)×2 matrix τ is defined by satisfying W= N τ, which is used to represent the elastic parameter. The transformation parameter θ contains the affine parameter A and elastic parameter W, which can be redefined as θ = [A N τ]. Then, the mapped model point set M is related to the model set M _{0} 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 \({u_{k}}\,=\,(m_{j}\,\,m_{i}) \left [\begin {array}{cc} 0 & 1\\1&0 \end {array}\right ]\). Based on above analysis, it is known m _{ i } =[1M _{0} Φ N]_{ i } θ, where [1M _{0} Φ N]_{ i } is the ith row of the matrix [1M _{0} Φ N]. We denote B _{ i } as [1M _{0} Φ N]_{ i }, and the normal direction of m _{ k } can be expressed as \({u_{k}}\,=\,T_{k}\theta \left [\begin {array}{cc}0 & \,\,1\\1&0 \end {array}\right ]\), where T _{ k } = B _{ j } − B _{ i }. Therefore, the surface structure description V _{ M } =[u _{1},⋯,u _{ m }]^{T} of the mapped model point set under parameters θ is
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 correlationbased 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, \(p\left (x \right)\,=\,\sum _{i=1}^{k}\alpha _{i}\phi (x\mu _{i},\sigma _{i})\), \(\sum _{i=1}^{k}\alpha _{i}\,=\,1\), where ϕ(xμ _{ i },σ _{ i }) 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 g m m(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 features is \(\int (gmm(v,V_{M})\,\,gmm(v,V_{S}))^{2}dv\). 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. \(\tiny \frac {\lambda }{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 enddiastole (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 closedform 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 gradientbased optimization techniques can be used to optimize J, such as the QuasiNewton method.
Optimization using the QuasiNewton method
QuasiNewton 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 \(\frac {\lambda }{2}W^{T}KW\), and J _{ v } represents the other two items.
Then,
where \(G\,=\,\frac {\partial J}{\partial M}\) is a m×2 matrix. J _{ v } can be written as \(J_{v} \,=\,\frac {\beta }{m^{2}} J_{1,v}\,\,\frac {2\beta }{mn}J_{2,v}\), where
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 QuasiNewton method can solve the nonlinear optimization problem, it might be divergent in high dimensional parameter space. In other word, QuasiNewton 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, QuasiNewton method is performed further to improve the accuracy of solution when SGD is convergent.
Optimization using Stochastic Gradient Descent
QuasiNewton 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 \( \frac {\partial f(\theta,i)}{\partial \theta } = \frac {\partial f_{1} }{\partial \theta } +\frac {\partial f_{2} }{\partial \theta }\,+\,\frac {\partial f_{3} }{\partial \theta }\,+\,\frac {\partial f_{4} }{\partial \theta }\,+\,\frac {\partial f_{5} }{\partial \theta } \). The derivative of each term is:
Based on the SGD procedure, Algorithm 1 represents details of optimization using SGD.
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, QuasiNewton 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 MICCAI 2009s 3D Cardiac Segmentation Challenge dataset [28]. In our study, the QuasiNewton method and the combination of SGD and QuasiNewton 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 QuasiNewton 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 pixelspacing of 0.931.64 mm [26].
We evaluate the registration accuracy of our algorithm by registering cardiac images at the enddiastole phase and the endsystole phase. The point sets of LV at enddiastole and endsystole phases are registered. Here, these points are provided by experts [26] in order to eliminate the effects of cardiac segmentation. The point set at enddiastole phase is the scene point set, and the one at endsystole phase is the model point set. The endsystole points are interpolated along the long axis to make the slice numbers equal to that in the enddiastole 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 endsystole and enddiastole 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 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 endsystole and enddiastole 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 QuasiNewton 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_SGD_QN further, 50^{2} control points are used to register the slices between endsystole and enddiastole 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 intersubject 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 shortaxis, enddiastolic cardiac MR images, acquired using a wholebody 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 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 MICCAI 2009s 3D Segmentation Challenge for Clinical Applications [28], provided by the Sunnybrook Health Sciences Centre. Cine steady state free precession (SSFP) MR short axis (SAX) images are obtained using 1.5T GE Signa MRI, during 1015 s breathholds with a temporal resolution of 20 cardiac phases over the heart cycle, and scanned from the enddiastolic phase. Six to 12 SAX images are obtained from the atrioventricular ring to the apex (thickness = 8 m m, gap = 8 m m, FOV = 320 m m×320 m m, matrix = 256×256) [28].
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 SavitzkyGolay 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=\frac {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 (HFI) or without (HFNI) 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 SCHYP37, 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 SCHFI05 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 QuasiNewton 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 QuasiNewton 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.
Abbreviations
 APD:

Average perpendicular distance
 CPD:

Coherent point drift
 DM:

Dice metric
 GMM:

Gaussian mixture model
 Icp:

Iterative closest point
 LV:

Left ventricle
 MR:

Magnetic resonance
 QN:

QuasiNewton
 ROI:

Region of interest
 SGD:

Stochastic gradient descent
 SSD:

Surface structure description
 TPS:

Thin plate splines
References
 1
Dharmakumar R. Assessment of regional myocardial oxygenation changes in the presence of coronary artery stenosis with balanced ssfp imaging at 3.0t: Theory and experimental evaluation in canines. J Magn Reson Imag. 2008; 27(5):1037–45.
 2
Tavakoli V, Amini AA. A survey of shapedbased registration and segmentation techniques for cardiac images. Comp Vision Image Underst. 2013; 117(9):966–89.
 3
Wang H, Amini AA. Cardiac motion and deformation recovery from mri: a review. IEEE Trans Med Imaging. 2012; 31(2):487–503.
 4
Makela T, Clarysse P, Sipila O, Pauna N, Pham QC, Katila T, Magnin IE. A review of cardiac image registration methods. IEEE Trans Med Imaging. 2002; 21(9):1011–21.
 5
LedesmaCarbayo MJ, Kybic J, Desco M, Santos A, Unser M. Cardiac motion analysis from ultrasound sequences using nonrigid registration. In: Medicae Computing and ComputerAssisted Intervention  Miccai 2001, International Conference, Utrecht, the Netherlands, October 1417, 2001, Proceedings. Utrecht: Springer;2001. p. 889–96.
 6
Chandrashekara R, Mohiaddin R, Rueckert D. Cardiac motion tracking in tagged mr images using a 4d bspline motion model and nonrigid image registration. In: IEEE International Symposium on Biomedical Imaging: From Nano To Macro. Arlington: IEEE;2004. p. 468–71.
 7
Chandrashekara R, Mohiaddin RH, Rueckert D. Analysis of 3d myocardial motion in tagged mr images using nonrigid image registration. IEEE Trans Med Imaging. 2004; 23(10):1245–50.
 8
Ebrahimi M, Kulaseharan S. Deformable image registration and intensity correction of cardiac perfusion mri. In: 5th International Workshop, STACOM 2014 Held in Conjunction with MICCAI 2014. Boston:2014. p. 13–20.
 9
Oubel E, De Craene M, Hero A, Pourmorteza A, Huguet M, Avegliano G, Bijnens B, Frangi AF. Cardiac motion estimation by joint alignment of tagged mri sequences. Med Image Anal. 2012; 16(1):339–50.
 10
Shi W, Zhuang X, Wang H, Duckett S, Luong DV, TobonGomez C, Tung K, Edwards PJ, Rhode KS, Razavi RS, et al. A comprehensive cardiac motion estimation framework using both untagged and 3d tagged mr images based on nonrigid registration. IEEE Trans Med Imaging. 2012; 31(6):1263–75.
 11
Bai W, Shi W, O’Regan DP, Tong T, Wang H, JamilCopley S, Peters NS, Rueckert D. A probabilistic patchbased label fusion model for multiatlas segmentation with registration refinement: application to cardiac mr images. IEEE Trans Med Imaging. 2013; 32(7):1302–15.
 12
Papademetris X, Sinusas AJ, Dione DP, Duncan JS. Estimation of 3d left ventricular deformation from echocardiography. Med Image Anal. 2001; 5(1):17–28.
 13
Escalanteramírez B, Moyaalbor E, Barbaj L, Cosio FA. Motion estimation and segmentation in ct cardiac images using the hermite transform and active shape models. Proc SPIE. 2013; 8856:219–24.
 14
Macan T, Loncaric S. 3d cardiac motion estimation by pointconstrained optical flow algorithm. In: International Symposium on Image and Signal Processing and Analysis. Pula:2001. p. 255–9.
 15
Shi P, Sinusas AJ, Constable RT, Ritman E, Duncan JS. Pointtracked quantitative analysis of left ventricular surface motion from 3d image sequences. IEEE Trans Med Imaging. 2000; 19(1):36–50.
 16
Besl PJ, McKay ND. A method for registration of 3d shapes. IEEE Trans Pattern Anal Mach Intell. 1992; 14(2):239–56.
 17
Chui H, Rangarajan A. A new point matching algorithm for nonrigid registration. Comput Vis Image Underst. 2003; 89(23):114–41.
 18
Myronenko A, Song X. Point set registration: Coherent point drift. IEEE Trans Pattern Anal Mach Intell. 2010; 32(12):2262–75.
 19
Jian B, Vemuri BC. Robust point set registration using gaussian mixture models. IEEE Trans Pattern Anal Mach Intell. 2011; 33(8):1633–45.
 20
Lin N, Duncan JS. Generalized robust point matching using an extended freeform deformation model: application to cardiac images. In: IEEE International Symposium on Biomedical Imaging: From Nano To Macro. Arlington: IEEE;2004. p. 320–3.
 21
Liu T, Liu W, Qiao L, Luo T. Point set registration based on implicit surface fitting with equivalent distance. In: IEEE International Conference on Image Processing. Quebec City:2015. p. 2680–4.
 22
Ravikumar N, Frangi AF. Robust groupwise rigid registration of point sets using tmixture model. Medical Imaging. 2016;9784.
 23
Du S, Liu J, Bi B, Zhu J, Xue J. New iterative closest point algorithm for isotropic scaling registration of point sets with noise. J Vis Commun Image Represent. 2016; 38(C):207–16.
 24
Tsin Y, Kanade T. A correlationbased approach to robust point set registration. In: European Conference on Computer Vision. Proceedings. Prague: Springer;2004. p. 558–69.
 25
Klein S, Pluim JP, Staring M, Viergever MA. Adaptive stochastic gradient descent optimisation for image registration. Int J Comput Vis. 2009; 81(3):227–39.
 26
Andreopoulos A, Tsotsos JK. Efficient and generalizable statistical models of shape and appearance for analysis of cardiac mri. Med Image Anal. 2008; 12(3):335–57.
 27
Stegmann MB. Active appearance models: Theory, extensions and cases. Inf Math Modell. 2000; 1(6):748–54.
 28
http://smial.sri.utoronto.ca/LV_Challenge. Accessed 22 Oct 2016.
 29
Wu Y, Ling H, Yu J, Li F. Blurred target tracking by blurdriven tracker. In: International Conference on Computer Vision. Barcelona: IEEE Computer Society;2011. p. 1100–7.
 30
Gorry PA. General leastsquares smoothing and differentiation by the convolution (savitzkygolay) method. Anal Chem. 1990; 62(6):570–3.
Acknowledgements
Not applicable.
Funding
Publication costs were funded by National Natural Science Foundation of China (U1301251) and Shenzhen Fundamental Research Program.
Availability of data and materials
Not applicable.
About this supplement
This article has been published as part of BMC Medical Informatics and Decision Making Volume 17 Supplement 3, 2017: Selected articles from the IEEE BIBM International Conference on Bioinformatics & Biomedicine (BIBM) 2016: medical informatics and decision making. The full contents of the supplement are available online at https://bmcmedinformdecismak.biomedcentral.com/articles/supplements/volume17supplement3.
Authors’ contributions
ZRZ participated in the design of the research, analysis of the data, and involved in drafting and revising of the manuscript. XY and GLC provided theoretical guidance and reviewed this paper. XY also conceived of the study, and participated in its design and data analysis. CT and WG contributed to the derivative process and data analysis. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Zhang, Z., Yang, X., Tan, C. et al. Surface structure feature matching algorithm for cardiac motion estimation. BMC Med Inform Decis Mak 17, 172 (2017). https://doi.org/10.1186/s129110170560z
Published:
Keywords
 Gaussian mixture model
 Surface structure feature
 Point set matching
 Stochastic gradient descent