International Journal of Automation and Computing  2018, Vol. 15 Issue (6): 673-688 PDF
Automated Segmentation of Left Ventricle Using Local and Global Intensity Based Active Contour and Dynamic Programming
G. Dharanibai1, Anupama Chandrasekharan2, Zachariah C. Alex1
1 School of Electronics Engineering, Vellore Institute of Technology (VIT) University, Vellore 632014, Tamil Nadu, India;
2 Department of Radiology, Sri Ramachandra University, Chennai 600116, Tamil Nadu, India
Abstract: The aim of this work is to develop an improved region based active contour and dynamic programming based method for accurate segmentation of left ventricle (LV) from multi-slice cine short axis cardiac magnetic resonance (MR) images. Intensity inhomogeneity and weak object boundaries present in MR images hinder the segmentation accuracy. The proposed active contour model driven by a local Gaussian distribution fitting (LGDF) energy and an auxiliary global intensity fitting energy improves the accuracy of endocardial boundary detection. The weightage of the global energy fitting term is dynamically adjusted using a spatially varying weight function. Dynamic programming scheme proposed for the segmentation of epicardium considers the myocardium probability map and a distance weighted edge map in the cost matrix. Radial distance weighted technique and conical geometry are employed for segmenting the basal slices with left ventricle outflow tract (LVOT) and most apical slices. The proposed method is validated on a public dataset comprising 45 subjects from medical image computing and computer assisted interventions (MICCAI) 2009 segmentation challenge. The average percentage of good endocardial and epicardial contours detected is about 99%, average perpendicular distance of the detected good contours from the manual reference contours is 1.95 mm, and the dice similarity coefficient between the detected contours and the reference contours is 0.91. Correlation coefficient and the coefficient of determination between the ejection fraction measurements from manual segmentation and the automated method are respectively 0.978 1 and 0.956 7, for LV mass these values are 0.924 9 and 0.855 4. Statistical analysis of the results reveals a good agreement between the clinical parameters determined manually and those estimated using the automated method.
Key words: Cardiovascular magnetic resonance     left ventricle     endocardium     epicardium     myocardium     segmentation     active contour     dynamic programming
1 Introduction

Cardiovascular disease (CVD) is the leading cause of death claiming 17.7 million lives a year globally[1]. The diagnosis and prognosis of heart condition rely upon various radiological image acquisition modalities, among which cardiac magnetic resonance (CMR) imaging using cine magnetic resonance image (MRI) sequences is well established and provides both anatomical and functional characterization of the heart[2]. The functional assessment of heart failure, congenital disease, and cardiomyopathies require the calculation of different clinical parameters such as ejection fraction (EF), myocardium mass and stroke volume. These parameters depend on the slice-by-slice and phase-by-phase delineation of endo and epicardial boundaries of the left ventricle (LV) from the CMR images. Manual tracing of these boundaries is time consuming and tedious task and leads to intra and inter observer variability resulting from the inaccuracies in tracing the complex myocardial structures. Therefore, the study of automated segmentation approaches that are accurate and requiring as little user interaction as possible is necessary.

Despite the wide variety of techniques and approaches proposed during the past two decades, automatic segmentation of left ventricle in cardiac MRI is still considered a difficult problem. A comprehensive review of left ventricle segmentation techniques can be found in [2, 3]. Among the various techniques presented in the literature, active contours represented by snakes, geometric active contours represented by level sets and graph based techniques are important, continue to be widely used, and have become focus of the contemporary research.

Segmentation techniques such as thresholding, edge detection and pixel based classification depend on pixel intensity. Intensity alone may not be sufficient to segment cardiac magnetic resonance (MR) images accurately as quality of these images are often affected by noise, low contrast, blurry edges, intensity overlapping and motion artifacts. These approaches are often combined with other segmentation methods to provide solution in two steps: segmentation of LV inner boundary in the first step followed by epicardium segmentation in the second step. Model based techniques such as active shape model, active appearance model (ASM/AAM)[48] and active contour model (ACM)[917], rely on the shape or relative distance between the endo and epicardium contours to segment both in a unified way. These approaches deform the contours under the influence of minimizing energy functional offering a great framework for the incorporation of prior information. ASM/AAM models require a manually built training dataset to constrain the segmentation, under the conditions of smeared boundaries or non-existence data. Active contour model (ACM) does not require any training set and drives the deforming contour towards the object boundary under the influence of internal and external forces. The external force deriving techniques from the image data are classified as edge based and region based. Edge based techniques use image gradient as a criterion to stop the evolving contour, which may leak into the background through the diffusive boundaries[1719]. Traditional region based ACMs use global intensity statistics and generally have better performance in the presence of weak or blurred boundaries than edge based methods. These models work on the assumption that the given image is homogeneous and hence lack the ability to deal with intensity inhomogeneity. Recently, some localized region based models have been proposed to overcome the difficulty of intensity inhomogeneity[1922]. However, these models easily get stuck in local minimum for the most contour initializations. Hybrid segmentation frameworks are now gaining popularity which integrates local edge or region information with the global region information in one generalized energy function for minimization[14, 15, 23]. While active contour model is fast but not accurate using global only information, it is accurate but not fast using local only information. The hybrid model segments the objects accurately and quickly.

Optimization methods such as dynamic programming (DP) and graph based algorithms transform the border detection process into an optimal path searching via weighted graph[2430]. Performance of DP is sensitive to the energy functional and is greatly improved by the choice of the weighting parameters of the data terms, smoothness terms and the priors incorporated in the fundamental formulations of algorithms.

In this study, due to the different natures of endo and epicardium, and the different segmentation challenges presented by them, two different techniques: 1) hybrid region based active contour that integrates local and global intensity characteristics in the energy functional and 2) DP framework incorporating probability information of the myocardium and the shape information, are proposed for segmenting LV endocardium (LVendo) and the LV epicardium (LVepi).

The main contributions of the proposed work include 1) formulation of new hybrid energy functional for handling intensity inhomogeneity and noise, 2) design of spatially varying weight function for controlling the auxiliary global intensity fitting term, 3) radial distance weighted approach for segmenting slices with left ventricle outflow tract (LVOT), 4) most apical slice segmentation using cone geometry, 5) design of new cost function for the DP incorporating myocardium probability information.

2 Materials and methods 2.1 Dataset and framework of segmentation algorithm 2.1.1 Dataset

Public data set of 45 subjects from medical image computing and computer assisted interventions (MICCAI) 2009 LV grand challenge[31] which is freely available on internet for the research purpose is used in this study for the design and testing of the proposed method. The dataset includes, 12 patients with ischemic heart failure (HF-I), 12 patients with non-ischemic heart failure (HF-NI), 12 patients with LV hypertrophy (HYP), and 9 normal (N) patients. It is also randomly divided into three groups: 15 online subjects, 15 training subjects and 15 validation subjects. Each data set was acquired as cine steady state free precession (SSFP) short-axis images during 10 – 15 seconds breath-holds. Temporal resolution of the images is 20 phases over a heart cycle with 6 – 12 short axis slices from base to apex. Field of view, image matrix size and slice thickness are respectively: 320 × 320 mm, 256 × 256 pixels and 8 – 10 mm. In plane, resolution is 1.25 – 1.52 mm.

2.1.2 Framework of segmentation algorithm

The proposed segmentation technique is composed of two stages: segmentation of LVendo in the first stage using hybrid region based active contour followed by LVepi segmentation using DP in the second stage. The whole process of slice by slice segmentation scheme is presented in the following:

Step 1. Load MRI data, and standardize the intensity.

Endocardium segmentation

Step 2. Select middle slice of the image volume. Localize the LV region and compute (region of interest) region of interest (ROI) using Fourier transform of temporal intensity characteristics.

Step 3. Use difference of Gaussian (DoG) filter and multi class Otsu threshold for the initial estimate of LV. Create a new ROI bounding the estimated LV with extension. Initialize level set using estimated LV.

Step 4. Segment LV endocardium using hybrid region based active contour. Smooth the convex hull boundary of the detected region.

Step 5. Repeat Step 4 to segment all frames in the given slice sequentially. In each frame initialize level set using LV segmentation from the preceding frame.

Step 6. Generate new ROI for the next slice using LV segmentations from all temporal frames of the preceding slice. Initialize level set using LV segmentation from the previous slice.

Step 7. Repeat Steps 4–6 to segment all slices sequentially from mid to base and then from mid to apex.

Step 8. Using spatiotemporal continuity detect the most basal slice and the most apical slice. Estimate LV in the basal slice using radial distance weighted method and in the apical slice using cone geometry.

Epicardium segmentation

Step 9. Transform the given slice into polar space. Estimate probability map of myocardium. Create a mask combining LV segmentation and the binarized probability map.

Step 10. Get edge map of polar image and exclude edge pixels belonging to non-myocardial structures using the binary mask from Step 9.

Step 11. Get gradient map of polar image. Design the cost function using gradient map, edge map and the distance constraint.

Step 12. Detect epicardial boundary using dynamic programming.

Step 13. Repeat Steps 9–12 to segment all slices and temporal phases.

2.2 Preprocessing

Cardiac MRIs present significant intensity variation across patients and scanners. Hence, intensity standardization method is required to reduce the inter subject intensity variation. The intensity standardization method presented by Nyúl et al.[32] is adopted in this work for transforming the intensity histogram of each image into a standard histogram scale [0, 255].

2.3 Computation of ROI

The main steps in localizing the heart region and region of interest (ROI) computation are shown in Fig. 1. Given the 3D+t cardiac cine MR image, let I (m, n, s, t) be the gray level of (m, n)-th voxel in the s-th slice at the t-th time frame, where 1 < (m, n, s, t) ≤ M, N, S, T. First step in the ROI detection is to select an image slice from the mid position s at end diastole (ED) phase. The selected image is blurred using Gaussian kernel of size (4×sigma+1) with sigma 60 mm and binarized using Otsu threshold segmentation. These steps detect the torso by suppressing the background as shown in the first row of Fig. 1.

Heart being the only moving organ in the field of view with static structures in the background, each pixel in the heart region shows a strong variation in intensity over time. The intensity variation in the heart region can be characterized using Fourier transform (FT) and recognized using harmonic images of cine sequences. FT of sequence of intensities over time at every pixel position is calculated as follows:

 $F(m,n,k) = \sum\limits_{t = 1}^T {I(m,n,s,t)} {{\rm{e}}^{ - {\rm{j}}2\pi \frac{k}{T}t}}{\rm{,}}\;\;{\rm{ }}k = 0,1, \cdots, T - 1.$ (1)

For each cine slice, the first harmonic (H1) magnitude image F(m, n, k), where k=1 is computed. Smoothing the H1 image in the detected torso using Gaussian kernel with sigma 50 mm followed by Otsu threshold segmentation generates a new region locating the heart as shown in second row in Fig. 1. For further precise localization of the cavity, H1 image given by the pixels of the newly detected region is smoothed with sigma 25 mm and then thresholded creating a binary image whose boundary tightly fits around the heart cavity as shown in the third row in Fig. 1.

To detect the LV cavity from this region, the mid-slice image is processed using difference of Gaussian (DoG) filter and its complement is converted into binary. This enhances the saliency of the ventricle blood pools. The two largest objects lying within the foreground region identified as right ventricle (RV) and LV are extracted after performing morphological opening, closing and region filling. Region of interest (ROI) is specified as a 90 × 90 pixel window, centered on the object appearing more circular. The size of the ROI is chosen to cover the entire LV and a small portion of RV and the window is propagated to all slices and phases for cropping. Results of these steps are shown in the fourth row in Fig. 1.

 Download: larger image Fig. 1. ROI computation. First row: mid-slice image, Gaussian-smoothing (σ = 60 mm) of mid-slice and its binary. Second row: H1 image, Gaussian smoothing (σ = 50 mm) of H1 image and its binary. Third row: Gaussian smoothing (σ = 25 mm) of non-background H1 image, its binary and the H1 image with region boundary. Fourth row: binary of DoG filtered mid-slice image, the two largest connected components in the detected region and the ROI window.

2.4 LV endocardium segmentation

In this section, to cope with weak object boundaries and intensity inhomogeneity, a new hybrid region based active contour that combines local and global intensity statistics in the energy functional is proposed. The influence of these two terms on the curve evolution is complementary. When the evolving curve is close to the object boundaries, the local intensity based force becomes dominant and drives the contour towards the object boundaries. Local energy term plays a key role in accurately locating object boundaries in images with intensity inhomogeneity. When the curve is far away from object boundaries, the global intensity fitting term becomes dominant and prevents the contour from falling into local minimum.

To have a better performance in the presence of noise, more complete statistical characteristics have to be considered. The local intensity information described using Gaussian distribution with spatially varying means and variances is incorporated in the active contour for locating the boundaries in the presence of noise and intensity inhomogeneity.

2.4.1 Level set formulation

Given a gray level image I: Ω $\subset$ R2, C is a closed contour, which separates the image domain Ω into two regions: Ω1 and Ω2. The contour C is represented implicitly by the zero level set of Lipschitz function ϕ: Ω → R, with $C = \{ x \in \Omega |\phi \left( x \right) = 0\}$ , $inside\left( C \right) = \{ x \in {\Omega _1}|\phi \left( x \right) > 0\}$ and $outside\left( C \right) = \{ x \in {\Omega _2}|\phi \left( x \right) < 0\}$ . The implicitly represented contour C is evolved according to the minimization of energy functional. The global and local intensity distribution model proposed in this work consists of three terms in the energy functional: global term Eglobal, local term Elocal, and the regularization term ER.

 ${E_{GLDF}} = \underbrace {\eta {\rm{ }}P(\phi ) + v{E_L}(\phi)}_{{E_R}} + {E_{local}}(\phi ) + \omega {E_{global}}(\phi )$ (2)

where η and υ are constants and ω is an adaptive weighting coefficient. The internal term P (ϕ) ensures that ϕ is not deviating from the signed distance function and the arc length term EL(ϕ) controls the smoothness of the contour.

 $P(\phi ) = \int_\Omega {\frac{1}{2}} {(\left| {\nabla \phi } \right| - 1)^2}{\rm{d}}x,\;{E_L}(\phi ) = \int_\Omega {\left| {\nabla H(\phi )} \right|} {\rm{d}}x.$ (3)

Heaviside function H(ϕ) indicates the region enclosed by C. The auxiliary global region term Eglobal(ϕ) is formulated using the last two terms of Chan–Vese[33] model.

 $\begin{split}{E_{global}}(\phi ) = & \,{\lambda _1}\int_{{\Omega _1}} {\left| {I - {c_1}} \right|} {{\rm{ }}^{\rm{2}}}H{\rm{(}}\phi {\rm{) d}}x+ \\ & {\lambda _2}\int_{{\Omega _2}} {\left| {I - {c_2}} \right|} {{\rm{ }}^{\rm{2}}}{\rm{(1}} - H{\rm{(}}\phi {\rm{)) d}}x.\end{split}$ (4)

where λ1 and λ2 are positive constants and c1, c2 are respectively the average of image intensities inside and outside of the contour C. The local intensity properties, described using local Gaussian distribution model by Wang et al.[22] are adopted to formulate the term Elocal(ϕ).

 $\begin{split}{E_{local}}(\phi ) = & \int_\Omega {\left( {\sum\limits_{i = 1}^2 {\int_\Omega { - {K_\sigma }(x - y)} } } \right.} \\ & \Bigg. {\log ({P_{i,x}}(I(y))){M_i}(\phi (y)){\rm{d}}y} \Bigg){\rm{d}}x\end{split}$ (5)

where M1(ϕ) = H(ϕ) and M2(ϕ) = 1–H(ϕ). I(y) is the intensity in the circular neighborhood Ox of radius ρ centered at point x, whose size can be controlled using Gaussian kernel function:

 \begin{aligned}{K_\sigma }(x - y) = \left\{ \begin{aligned}& \displaystyle\frac{1}{a}\exp \left( { - \displaystyle\frac{{|x - y{|^2}}}{{2{\sigma ^2}}}} \right),\;\;{\rm{if}}\;\;|x - y|\; \le \rho \\& \;\;\;\;\;0,\;\;\;\;\;\;{\rm{else}}\;\;\end{aligned} \right.\end{aligned} (6)

where a is a constant and $\int_{Ox} {{K_s}(x-y){\rm{d}}y} = 1$ . Probability density Pi, x(I(y)) of intensity in the partition i of the neighborhood Ox is modeled using Gaussian distribution.

 ${P_{i,x}}(I(y)) = \frac{1}{{\sqrt {2\pi } {\sigma _i}(x)}}\exp \left( { - \frac{{{{(I(y) - {u_i}(x))}^2}}}{{2\sigma _i^2(x)}}} \right)$ (7)

where ui and σi are the local intensity means and variances of each partition i. The function ϕ that minimizes EGLDF(ϕ) in (2), is solved using the gradient descent scheme on the regularized Euler-Lagrange equation, resulting in the following curve evolution:

 $\begin{array}{l}\displaystyle\frac{{\partial \phi }}{{\partial t}} = - \displaystyle\frac{{\partial E\left( \phi \right)}}{{\partial \phi }} = \delta (\phi )\left( {{F_1} + \omega {F_2}} \right) + \upsilon \delta (\phi ){\rm{div}}\left( {\displaystyle\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right) + \\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\eta \left[ {{\nabla ^2}\phi - {\rm{div}}\left( {\displaystyle\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right)} \right]\end{array}$ (8)

where δ(ϕ) is the dirac $\delta$ function. F1 and F2 are respectively the spatially varying local and global intensity fitting forces and are defined as

 $\begin{split}& {F_1}\! = -\!{\lambda _1}\!\!\int_\Omega \!\!{{K_\sigma }(x - y)} \left(\! {\log {\sigma _1}(y) + \frac{{{{(I(x) - {u_1}(y))}^2}}}{{\sigma _1^2}}} \!\right){\rm{d}}y +\!\!\! \\ & \quad\quad{\lambda _2}\int_\Omega {{K_\sigma }(x - y)} \left( {\log {\sigma _2}(y) + \frac{{{{(I(x) - {u_2}(y))}^2}}}{{\sigma _2^2}}} \right){\rm{d}}y\\ &{F_2} = - {\left| {I(x) - {c_1}} \right|^2} + {\left| {I(x) - {c_2}} \right|^2}.\end{split}$ (9)

The local intensity means $u_1$ , $u_2$ , variances σ1, σ2 and the global average intensities $c_1$ , $c_1$ are the minimizers of the functional EGLDF(ϕ) and are updated according to (10), throughout the iterations as the contour evolves.

 $\begin{split}& {u_i} = \frac{{\displaystyle\int {{K_\sigma }(x - y)I(y){M_i}\left( {\phi (y)} \right){\rm{d}}y} }}{{\displaystyle\int {{K_\sigma }(x - y){M_i}\left( {\phi (y)} \right){\rm{d}}y} }}\\ & {\sigma _i} = \frac{{\displaystyle\int {{K_\sigma }(x - y){{\left( {{u_i}(x) - I(y)} \right)}^2}{M_i}\left( {\displaystyle\phi (y)} \right){\rm{d}}y} }}{{\displaystyle\int {{K_\sigma }(x - y){M_i}\left( {\phi (y)} \right){\rm{d}}y} }}\\& {c_i} = \frac{{\displaystyle\int {I(y){M_i}\left( {\phi (y)} \right){\rm{d}}y} }}{{\displaystyle\int {{M_i}\left( {\phi (y)} \right){\rm{d}}y} }}.\end{split}$ (10)

For numerical implementation, H(ϕ) and δ(ϕ) are regularized as

 $\begin{split}& {H_\varepsilon }\left( \phi \right) = \frac{1}{2}\left[ {1 + \frac{2}{\pi }\arctan \left( {\frac{x}{\varepsilon }} \right)} \right] \\ & {\delta _\varepsilon }\left( \phi \right) = \frac{1}{\pi }\left( {\frac{\varepsilon }{{{x^2} + {\varepsilon ^2}}}} \right)\end{split}$ (11)

where ε is a constant regularizing the functions.

Images severely affected by the intensity inhomogeneity, should be segmented using the local fitting force, reducing the influence of global force. Otherwise, the global force dominates, preventing the evolving contour being stopped at the object boundaries. In order to control the influence of global term, it is necessary to select the value of ω in (2), according to the degree of intensity inhomogeneity. Rather selecting ω as a constant for the entire image domain, an adaptive function, which provides different values at different pixels is proposed. The proposed dynamically varying weight function is calculated based on the local intensity contrast:

 $\omega = \gamma {\rm{average}}\left( {{C_{ox}}} \right)\left( {1 - {C_{ox}}} \right)$ (12)

where γ is a constant and Cox is the image contrast in the local neighbourhood ox of the given image:

 $\begin{split} & {C_{ox}}(x) = \sqrt {\frac{{\displaystyle\int_{{O_x}} {{K_\sigma }\left( {x - y} \right){{\left( {\frac{{I(y) - m(x)}}{{Ig}}} \right)}^2}{\rm{d}}y} }}{{\displaystyle\int_{{O_x}} {{K_\sigma }\left( {x - y} \right){\rm{d}}y} }}}\\ & m(x) = \sqrt {\frac{{\displaystyle\int_{{O_x}} {{K_\sigma }\left( {x - y} \right)I(y){\rm{d}}y} }}{{\displaystyle\int_{{O_x}} {{K_\sigma }\left( {x - y} \right){\rm{d}}y} }}} .\end{split}$ (13)

Gaussian kernel Kσ controls the size of the neighborhood region Ox. Ig is the maximum gray-level in the image and m(x) is the local mean. Cox is smaller in regions with smooth intensity changes and larger in regions with gradients. Averaging Cox over the entire image represents the overall contrast. The factor (1 – Cox) dynamically adapts the value ω. For the regions with high local contrast, ω becomes smaller, and for low contrast regions, this becomes larger.

2.4.3 Endocardial contour extraction

LV endocardium in each slice and phase of the 4D cardiac MRI is segmented using the hybrid region based active contour proposed in the previous section. The algorithm segments the mid slice first, proceeds to segment other slices in the basal direction, and then the apical slices. Temporal phases are segmented in sequence from the first frame.

To obtain faster and accurate segmentation, the initial contour is placed close to the LV boundary[34]. If the LV region can be roughly and automatically obtained using techniques such as thresholding, then the roughly obtained LV region can be used to construct the initial level set ϕo. For this, ROI of the mid slice (Fig. 2(a)) is divided into four classes: blood, soft tissues, background, and the bright artefacts, using multilevel Otsu thresholding method (Fig. 2(b) (top)). The binary of DoG filtering of mid slice image provides a mask as shown in Fig. 2(b) (bottom). When the multi class image is masked using the binary image, it detaches the cavity from the erroneous structures, enhancing the saliency of the LV cavity. The masked image in Fig. 2(c) is decomposed into three independent images: an image with the brighter class alone, an image with two brightest classes combined, an image with third brighter class alone. By analyzing these three images, the object which is the most circular and close to the center is chosen as LV (Fig. 2(d)). The irregularities in the detected LV region are smoothed by morphological closing using disk shaped structuring element followed by region filling. The resulting smooth LV region is labelled as B and is used for the level set initialization:

 ${\phi _o} = {c_o}\left( {0.5 - B} \right){\rm{ }}$ (14)

where co is a constant. To reduce the size of the computational domain, a new ROI is defined by fitting a rectangular window around the labelled region B with 8 pixels safety margin. The new ROI is thus obtained overlaid with the initial contour ϕ0 and shown in Fig. 2(e). According to (8), the function ϕt evolves as given by Fig. 2(f) from its initial position until convergence. In addition to the cavity, the final contour is found enclosing other brighter regions. In the post processing, to get the final segmentation as given by Fig. 2(g), LV blood region appearing circular and very close to the center of the image alone is retained with the holes filled.

 Download: larger image Fig. 2. Endocardium segmentation (a) ROI; (b) Multilevel threshold (top) and binary of DoG filtered image (bottom); (c) Masked multiclass image; (d) Initial estimate of LV; (e) Initial contour overlaid on new ROI; (f) Final contour; (g) Final segmentation. (Color versions of the figures in this paper are available online)

Papillary muscles are excluded from the myocardium by computing the convex hull of the segmented object. To obtain a smooth boundary, the edge points of the convex hull are transformed into the Fourier domain. Retaining the first five harmonics and computing the inverse Fourier transform of the resulting Fourier coefficients give a smooth boundary. Fig. 3 shows the detected contours before and after smoothing superimposed on the original ROI image.

 Download: larger image Fig. 3. Final contour before (red) and after (cyan) Fourier smoothing

For the given slice position, the remaining images in the temporal sequence that covers the complete cardiac cycle are sequentially processed using the same procedure described above. The new ROI computed for the first frame in the given slice position is used for the rest of the frames. In each cardiac phase, the LV cavity segmented from the preceding phase is used for the initial estimation of LV. After converting the image of the current cardiac phase into binary, the object having maximum overlap with the LV cavity segmented in the previous phase is used as initial estimate of LV. The endocardium segmentations from all temporal phases of the present slice position are maximum intensity projected (MIP) to create the new ROI for the images in the next slice position. The rectangular region fit around the MIP image with eight pixels safety margin, provides the new ROI for the next slice. The whole process explained above is repeated on the remaining slices covering the whole heart from mid-level to basal and from mid-level to apical.

In basal slices with left ventricle outflow tract (LVOT), myocardium appears as an incomplete annular shape. In the most apical slices, blood lumen is too small to be resolved. In these slices, the level set based contour evolution leaks into the myocardium. Hence, separate techniques are proposed to segment the slices with LVOT and the most apical slices showing very small blood region.

2.4.4 Segmentation of basal slices with LVOT using radial distance weighted method

While processing the basal slices, if major axis length of blood region segmented in the present slice exceeds that of the preceding slice by 1.2 times, then the present slice is identified as basal slice with LVOT (Fig. 4(a)). For estimating the endocardial boundary, this basal slice is converted into binary using Otsu′s bi-level threshold and an edge map is computed using canny edge detector as shown in Figs. 4(b) and 4(c). The dilated LV mask from the preceding slice weighted by the radial distance from its centroid (Fig. 4(d)) is used to mask edge map as shown in Fig. 4(e). This suppresses the edge pixels outside the mask to zero while the inside pixels are weighted by its radial distance from the centre point. Edge pixels whose radial distance is greater than one standard deviation from the mean distance are considered non endocardium and are excluded (Fig. 4(f)). The convex hull of the remaining pixels forms a region as shown in Fig. 4(g) whose boundary defines the endocardium as shown in Fig. 4(h).

 Download: larger image Fig. 4. Segmentation of basal slice with LVOT. (a) Basal slice; (b) Binary image of basal slice; (c) Edge map; (d) LV mask weighted by distance; (e) Distance weighted edge map; (f) Potential edge map of endocardium; (g) Convex hull; (h) Final segmentation of LV.

2.4.5 Segmentation of most apical slices

While processing the slices from mid to apex, the segmented regions in the consecutive slices exhibit a smooth change in the area. In most apical slices, the blood lumen is too small to be detected causing the evolving contour to leak out of apex. This gives a sudden rise in the area segmented by the active contour. Apical slices showing such discontinuity in the area can be detected as follows: if the centre of mass of the current slice segmentation is displaced by 7 pixels from that of the previous slice or if the ratio between the current slice and previous slice area is greater than 1.4, then apical slice with area discontinuity is detected. For these apical slices, segmentation is corrected using the cone geometry: For the given slice near to apex, at position j and at phase k, the radius R of the LV cavity is related to the radius of its two preceding slices as follows:

 $R(I(m,n,j,k)) = 2 \times \frac{{R(I(m,n,j - 1,k))}}{{R(I(m,n,j - 2,k))}}.$ (15)

A circular binary mask created using the radius estimated from (15), gives the corrected segmentation for the most apical slices detected with area discontinuity. Figs. 5(a) and 5(b) show the segmentation of LV in the two preceding slices and Fig. 5(c) shows the contour leakage (green) in the most apical slice. The segmentation corrected using the cone geometry is shown in red.

 Download: larger image Fig. 5. Segmentation in apical slices. (a) Slice j-2; (b) Slice j-1; (c) Contour leakage in the slice j (green) corrected using conical geometry (red).

2.5 LV epicardium segmentation

Epicardial boundaries are detected using dynamic programming (DP) in the polar representation of images. In polar space, the roughly circular shape of LV can be visualized as a linear structure and the boundary tracing problem can be solved by linear searching.

Pixels along the radial scan lines drawn from the barycentre of the LV at uniformly spaced angles form the columns in the polar image. Thus, x-axis of the polar image represents the angle θ in [–π, +π], and y-axis represents the radius r in [0, R]. Given the ROI image of short axis slice (Fig. 6(a)), the radial line pointing to the left forms the first column of the polar image as shown in Fig. 6(b) (top). The boundary of endocardium which is circular in the original ROI appears approximately linear in the polar representation.

2.5.1 Myocardium probability map computation

In the dynamic programming based boundary tracing problems, design of cost matrix is crucial for the accurate image segmentation. In this work, a new cost matrix based on myocardial probability map is proposed for segmenting epicardium. The probability map of myocardium is estimated based on the assumption that intensity distribution is Gaussian in nature.

 ${P_M}(x) = \frac{1}{{\sqrt {2\pi } \sigma }}{\rm{ }}\exp \left( { - \frac{{{{(x - m)}^2}}}{{2{\sigma ^2}}}} \right)$ (16)

where m and σ are the mean and standard deviation of the gray level distribution in the myocardium. To estimate the mean and standard deviation, pixels belonging to myocardium are sampled using dilated masks of LV blood pool. The LV mask dilated using structuring element of size k is subtracted from that of the mask dilated using the structuring element of size k + 3. The resulting mask is overlaid on the gray level image to sample the intensities belonging to myocardium. The mean and standard deviation of the samples in the masked region are used in (16), to estimate myocardium probability map. The polar transform of the estimated probability map PM(r, θ) is shown in Fig. 6(b).

 Download: larger image Fig. 6. Polar form of ROI. (a) ROI with radial lines; (b) Polar map I(r, θ) of ROI (top) and probability map PM(r, θ) of myocardium (bottom).

2.5.2 New cost function

The possible boundary lines in the polar map are considered as polylines, ${V_N} = \left\{ {{v_1},{\rm{ }}{v_2}, \cdots ,{v_N}} \right\}$ having N vertices with one vertex on each column of the image grid. The coordinate of each vertex is (r, θ). The cost function Csum(VN) of the polyline VN is defined as the sum of local costs along a candidate boundary VN:

 ${C_{\rm sum}}({V_N}) = C({v_1}) + \sum\limits_{i = 2}^N {C({v_i})} + D({v_i}_{ - 1},{v_i}).$ (17)

For a node vi, local cost C(vi) is built using the following image information, edge map, probability map of myocardium (PM), and gradient magnitude. All these image features are normalized to the range of [0, 1]. The term D(vi-1, vi) is related to the geometrical distance between the vertical positions of vi-1 and vi. Larger distance yields higher cost value. This term prevents the random fluctuations along the optimal path and keeps the line smoother. The desired boundary VN is the optimal path minimizing Csum(VN).

 ${{\widehat V}_N} = \left\{ {{V_N}\left| {\min ({C_{{\rm{sum}}}}({V_N}))} \right.} \right\}.$ (18)

For building the cost function C(vi), an edge map Edge(r, θ) shown in Fig. 7(a) is computed by combining the edges detected from polar image I(r, θ) and the probability map PM(r, θ). The binary image of probability map shown in Fig. 7(b) together with the LV mask segmented in the first stage, produce a mask which is used to filter the edges belonging to non-myocardial structures as follows:

 E(r,\theta ) = \left\{ \begin{aligned}& 0,\;\;\;\;\;\;\;\;\;\;{\rm{if }}\;\;\;\;\;\;\;\;mask(r,\theta ) = {\rm{1}}\\& Edge{(r,}\;\;\theta ),\;\;\;\;\;{\rm{otherwise}}{\rm{.}}\end{aligned} \right. (19)

The resulting filtered edge map E(r, θ) is shown in Fig. 7(c). In each column of E(r, θ), except the first edge pixel, all other pixels are suppressed to zero. Given the row coordinates ${\{ {\rho _\theta }\} _{\theta = 1, \cdots ,N}}$ of the surviving edge pixels, the mean radial distance of the potential edges is represented by ρm and the standard deviation by ρσ. Every pixel in the edge map E(r, θ) is weighted by its distance from the mean value ρm according to (20), and the result Wedge(r, θ) is shown in Fig. 7(d).

 ${W_{edge}}(r,\theta ) = E\left( {r,\theta } \right)\times \exp \left( {\frac{{ - 0.5\times|{\rho _\theta } - {\rho _m}|}}{{{\rho _\sigma }}}} \right).$ (20)

The gradient map for building the cost function is derived by combining the vertical gradients of PM(r, θ) and I(r, θ). The vertical gradient of PM(r, θ) has strong negative values at the epicardial boundary while that of I(r, θ) may have positive or negative values due to different nature of tissues surrounding epicardium. The combined gradient cost image is derived as follows:

 {C_G}(r,\theta ) = \left\{ \begin{aligned} & - \left| {{\nabla _r}I(r,\theta ){\rm{ }}} \right| + {\nabla _r}{P_M}(r,\theta ), \\& {\rm{if}}\;\;\;\; ({\nabla _r}{P_M}(r,\theta ){\rm{ )}} < {\rm{0}}\\& {\rm{0}},\;\;\;\;{\rm{otherwise}}\end{aligned} \right. (21)

where $\nabla r$ represents gradient operation in the radial direction. The influence of papillary muscles and other small intensity fluctuations in the endocardium are suppressed by masking the cost image using the segmented LV region. The masked gradient cost image CG(r, θ) is shown in Fig. 7(e). Finally, the gradient based local cost is weighted by the edge map according to (22) and the result is shown in Fig. 7(f).

 $C(r,\theta ) = {W_{edge}}\times{C_G}(r,\theta ).$ (22)

If the vertical distance between the nodes vi–1 and vi is large, the smoothness constraint D(vi–1, vi) assigns large cost value, thus penalizes the big jumps along the optimal boundary. The distance cost function D(vi–1, vi) is defined as follows:

 $D({v_{i - 1}},{v_i}) = dist(l) = 1 - \exp ( - 0.5\times|l|)$ (23)

where l is the difference between the vertical positions of the nodes vi–1 and vi.

2.5.3 Dynamic programming

The optimal boundary minimizing the cost Csum(VN) can be searched recursively using DP. In recursive form, the candidate minima of the cost function in (22) can be expressed as a multistage cost accumulation process:

 \left\{ \begin{aligned}& {\rm{Acost }}(r,\theta + 1) = \mathop {\min }\limits_{ - m \le l \le m} {\rm{Acost }}(r + l,\theta ) +\\ &\qquad C{\rm{ }}(r,\theta + 1) + dist(l){\rm{,}}\;\;{\rm{for}}\;\; \theta = {\rm{2}}, \cdots, N - {\rm{1}}\\& {\rm{Acost }}(r,{\rm{ }}1) = C(r,{\rm{ }}1),\;\;\;\;{\rm{for}}\;\;\theta = {\rm{1}}\end{aligned} \right. (24)

where m is the range of nodes in each column within which a node from the previous column is allowed to move. Thus, in the forward recursion, accumulated cost of the visiting node is constructed as the sum of the minimum costs of all predecessors of the given node, the own cost, and the distance cost of the visiting node. The small values in the cumulative cost indicate highly likely edges. The optimal path traced back starting from the minimal cost of the last column to the minimal cost of the first column, provides the epicardial boundary. The recursively searched minimal cost path superimposed on accumulated cost is shown in Fig. 7(g). The detected optimal epicardial boundary superimposed on the edge map and the polar image are respectively shown in Figs. 7(h) and 7(i). The epicardial boundary is then transformed back to the Cartesian space. Epicardial boundaries are smoothed in a similar way as that of endocardium using Fourier transform. The detected epicardial boundary before and after Fourier smoothing, superimposed on the original ROI is shown in Fig. 8.

 Download: larger image Fig. 7. Cost function and epicardium segmentation. Endocardial contour superimposed in red (a) Edge map; (b) Binary of myocardial probability map; (c) Filtered edge map; (d) Distance weighted edge map; (e) Masked gradient cost; (f) Gradient and edge combined local cost; (g) Accumulated cost with the traced minimum cost path (green); (h) Epicardial contour (minimum cost path) superimposed on edge map; (i) Endo and edpicardial contours superimposed on polar form of ROI.

 Download: larger image Fig. 8. Endo and epicardial contours of LV transformed back into Cartesian space. Epicardial contours before (white) and after (green) FFT smoothing.

3 Results and discussions

The developed algorithm has been implemented in Matlab 7.2 and was run on a HP workstation with Intel Xeon quad core processor at 3.2 GHz having 12 GB RAM. In the level set formulation, the time step τ is chosen to be 5.0 satisfying the condition $\displaystyle\frac{\tau }{\eta } \le 0.2$ . The weighting parameters λ1, λ2 and the regularizing parameter ε are set to 1.0 and ν is set to 6.5. The neighborhood function Kσ is selected as a truncated Gaussian kernel of size K×K, where K = 4σ + 1 with the value of σ set to 3. The parameter γ in the adaptive weighting function is chosen as 0.001.

The two stage algorithm segments both endo and epicardium from all temporal phases at each slice position enabling the computation of LV volume and mass over the complete cardiac cycle.

Fig. 9 shows the sample segmentation of the subject SC-HF-I-05 using the proposed framework. The automatically detected contours are shown in solid white lines and the ground truth reference contours are shown in colour. The expert′s drawn reference contours for endocardium (cyan) are available at ED and end systole (ES) phases of all slices and the epicardium reference contours (red) are available only at ED phase. The volumes enclosed by endocardium (LVendo) and epicardium (LVepi) in every phase of the cardiac cycle are calculated as the sum of areas enclosed by the respective contours from base to apex. The LVM is computed as the product of myocardial volume (LVepi – LVendo) and the myocardial tissue density 1.05 g.

 Download: larger image Fig. 9. Sample volume segmentation of LV, white contours denote the automatically detected endo and epicardial contours, cyan and red contours denote the ground truth

Fig. 10 shows some segmentation outputs of the proposed method and the ground truth. For each subject, the left image is of ED phase and the right image is of ES phase. The image pairs in each row are from online, validation and training datasets. Each column presents the results from four pathology groups: ischemic (HF-I), non-ischemic (HF-NI), hypertrophy (HYP) and normal (N).

 Download: larger image Fig. 10. Segmentation outputs by the proposed method. Results in the first pair of columns are from training dataset, the second pair of columns is from validation dataset and the third pair of columns is from online dataset. In each subject, left image is from ED phase and the right image is from ES phase. The red and green contours indicate endo and epicardium results of the proposed method, cyan and yellow contours represent the ground truth.

Performance of the segmentation algorithm is evaluated using two metrics: average perpendicular distance (APD), and the dice similarity coefficient (DSC). Contours detected with an average distance less than 5 mm from the reference contours are considered good contours. Only the good contours are considered for evaluating the performance of the segmentation framework. In addition to this, the clinical parameters namely end-diastolic volume (EDV), end-systolic volume (ESV), EF and LV mass (LVM) at ED are computed to validate the performance of the algorithm. These parameters are selected based on the fact that the clinical experts have delineated the endocardium at diastole and systole and the epicardium in diastole only.

Quantitative assessment of the segmented contours performed on the 45 cases of the public data is reported in Table 1. The results categorized by the subject group and the subject pathology are summarized in Table 2. The algorithm has detected 99% of contours with APD<5 mm. APD of the detected good endo and epicardium are respectively 1.941 ± 0.452 mm and 1.958 ± 0.412 mm. Similarity measure DSC is found to be respectively 0.883 ± 0.054 and 0.931 ± 0.02 for the endo and epicardium. This shows a good agreement between manual and automatic delineation. The bad contours are mostly detected in the subjects with hypertrophy. Due to enlarged and thickened heart walls, LV cavities in these subjects appear very small causing contour leakage through the myocardial wall. The conical geometry interpolation method described in Section 2.4.5, overcomes the difficulty of segmenting apical slices to some extent.

Table 1
Segmentation results on 45 cases of MICCAI public data set

Table 2
Quantitative performance evaluation of the proposed segmentation method

Linear regression and Bland-Altman analysis are performed to compare the clinical parameter estimations between manual method and automated method as shown in Fig. 11. Statistical measures such as correlation, bias, limits of agreement, absolute error, signed error and root mean square error (RMSE) between manual and automated methods are presented in Table 3.

Table 3
Statistical analyses of clinical parameters measured using the proposed segmentation method

 Download: larger image Fig. 11. Bland-Altman (first column) and linear regression (second column) plots for EDV, ESV, EF and LVM

Performance of the proposed algorithm is compared with the state-of-the-art methods in the literature that used the same database from MICCAI 2009 segmentation challenge. Table 4 shows comparison of the segmentation results in terms of DSC, APD and percentage of good contours. Linear regression and Bland-Altman analysis of the clinical parameters measured by the proposed method and state-of-the-art methods in the literature are presented in Table 5. The proposed method exhibits good correlation with manual measurements. The regression slope in EF measurement is approximately one with a bias of –1.72%. However, the correlation coefficients are good in the EF measurement and are better than the other methods in the literature. The correlation coefficient and regression slope are good in measuring LVM, but the bias from the Bland-Altman analysis is –11.7 g. The coefficient of determination of EF and LVM are 0.956 7 and 0.855 4. Therefore, algorithm is accurate in computing EF and LVM.

Table 4
Comparison of segmentation accuracy among the state-of-the-art methods and the proposed method

Table 5
Linear regression and Bland-Altman analysis of the proposed method compared to the state-of-the-art methods

LV volumes EDV and ESV are underestimated by 11.1 ml and 7.58 ml leading to a bias in LVM and EF measurements. The bias in the estimation of clinical parameters is mainly due to the fact that, contours which are considered bad have not been included in the volume calculations. The proposed method is found good in computing EF than LVM. The bad contours in endocardium segmentation are mostly found in the subjects SC-HYP-07, SC-N-05 and SC-N-10. In epicardium segmentation, bad contours are mostly found in the subjects SC-HF-NI-14, SC-HF-NI-15, SC-HF-NI-33 and SC-N-10. Among the detected good contours, epicardium from the hypertrophy group has largest APD. This is due to wall enlargement and wall thickening causing the densely packed papillary muscles and trabeculae hardly distinguishable from myocardium. Compared to other groups, subjects from HF-I group are segmented with least APD. Overlap measure indicates that epicardium is segmented with better accuracy than the endocardium.

4 Conclusions

The algorithm proposed in this paper is a two stage process, which segments endo and epicardium in all slices and phases of the 4D cardiac MRI. A comprehensive approach based on the first harmonic of temporal intensity variation, Gaussian smoothing and DoG filtering techniques are used for the precise localization of LV region and ROI computation. The local and global energy based hybrid active contour segments LV endocardium more accurately in the first stage followed by epicardium in the second stage using dynamic programming framework. The main advantages of the proposed segmentation methods are that they are fully automatic, do not require manual intervention and are data driven.

In the dynamic programming framework, a probability map of the myocardium is constructed to suppress the edge features of non-myocardium tissues in the edge map. DP is applied to the region close to the myocardium to exclude the irrelevant tissues and therefore undesirable structures are not included in the search range.

Segmentation results of the proposed method demonstrated the conformity between the automatically estimated clinical parameters and those computed from the manual segmentations. In conclusion, the proposed method can be used to estimate cardiac functional parameters automatically during computer aided diagnosis.

Acknowledgements

This work was supported by Department of Science and Technology, Ministry of Science and Technology, India (No. DST/TSG/ICT/2010/08). The authors are grateful to the management of Vellore Institute of Technology (VIT) University and Sri Ramachandra University for their help and guidance. The authors are also thankful to the reviewers for their comments and suggestions which helped in improving the paper.

References
 [1] World Health Organization. Cardiovascular diseases (CVDs). WHO Fact Sheet, No. 317. Geneva, Switzerland: World Health Organization, 2017. [2] P. Peng, K. Lekadir, A. Gooya, L. Shao, S. E. Petersen, A. F. Frangi. A review of heart chamber segmentation for structural and functional analysis using cardiac magnetic resonance imaging. Magnetic Resonance Materials in Physics, Biology and Medicine, vol.29, no.2, pp.155-195, 2016. DOI:10.1007/s10334-015-0521-4 [3] C. Petitjean, J. N. Dacher. A review of segmentation methods in short axis cardiac MR images. Medical Image Analysis, vol.15, no.2, pp.169-184, 2011. DOI:10.1016/j.media.2010.12.004 [4] S. C. Mitchell, J. G. Bosch, B. P. F. Lelieveldt, R. J. van der Geest, J. H. C. Reiber, M. Sonka. 3D active appearance models: Segmentation of cardiac MR and ultrasound images. IEEE Transactions on Medical Imaging, vol.21, no.9, pp.1167-1178, 2002. DOI:10.1109/TMI.2002.804425 [5] S. C. Mitchell, B. P. F. Lelieveldt, R. J. van der Geest, H. G. Bosch, J. H. C. Reiver, M. Sonka. Multistage hybrid active appearance model matching: Segmentation of left and right ventricles in cardiac MR images. IEEE Transactions on Medical Imaging, vol.20, no.5, pp.415-423, 2001. DOI:10.1109/42.925294 [6] H. C. van Assen, M. G. Danilouchkine, A. F. Frangi, S. Ordas, J. J. M. Westenberg, J. H. C. Reiber, B. P. F. Lelieveldt. SPASM: A 3D-asm for segmentation of sparse and arbitrarily oriented cardiac MRI data. Medical Image Analysis, vol.10, no.2, pp.286-303, 2006. DOI:10.1016/j.media.2005.12.001 [7] H. H. Zhang, A. Wahle, R. K. Johnson, T. D. Scholz, M. Sonka. 4D cardiac MR image analysis: Left and right ventricular morphology and function. IEEE Transactions on Medical Imaging, vol.29, no.2, pp.350-364, 2010. DOI:10.1109/TMI.2009.2030799 [8] S. Zambal, J. Hladuvka, K. Bühler. Improving segmentation of the left ventricle using a two-component statistical model. In Proceedings of the 9th International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, Copenhagen, Denmark, pp. 151–158, 2006. [9] I. Ben Ayed, S. Li, I. Ross. Embedding overlap priors in variational left ventricle tracking. IEEE Transactions on Medical Imaging, vol.28, no.12, pp.1902-1913, 2009. DOI:10.1109/TMI.2009.2022087 [10] R. El Berbari, I. Bloch, A. Redheuil, E. Angelini, E. Mousseaux, F. Frouin, A. Herment. An automated myocardial segmentation in cardiac MRI. In Proceedings of the 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, IEEE, Lyon, France, pp. 4508–4511, 2007. [11] G. Hautvast, S. Lobregt, M. Breeuwer, F. Gerritsen. Automatic contour propagation in cine cardiac magnetic resonance images. IEEE Transactions on Medical Imaging, vol.25, no.11, pp.1472-1482, 2006. DOI:10.1109/TMI.2006.882124 [12] M. P. Jolly. Automatic segmentation of the left ventricle in cardiac MR and CT images. International Journal of Computer Vision, vol.70, no.2, pp.151-163, 2006. DOI:10.1007/s11263-006-7936-3 [13] M. Lynch, O. Ghita, P. F. Whelan. Segmentation of the left ventricle of the heart in 3D+T MRI data using an optimized nonrigid temporal model. IEEE Transactions on Medical Imaging, vol.27, no.2, pp.195-203, 2008. DOI:10.1109/TMI.2007.904681 [14] M. Lynch, O. Ghita, P. F. Whelan. Left-ventricle myocardium segmentation using a coupled level-set with a priori knowledge. Computerized Medical Imaging and Graphics, vol.30, no.4, pp.255-262, 2006. DOI:10.1016/j.compmedimag.2006.03.009 [15] N. Paragios. A variational approach for the segmentation of the left ventricle in cardiac image analysis. International Journal of Computer Vision, vol.50, no.3, pp.345-362, 2002. DOI:10.1023/A:1020882509893 [16] C. Pluempitiwiriyawej, J. M. F. Moura, Y. J. L. Wu, C. Ho. STACS: New active contour scheme for cardiac MR image segmentation. IEEE Transactions on Medical Imaging, vol.24, no.5, pp.593-603, 2005. DOI:10.1109/TMI.2005.843740 [17] A. Gupta, L. von Kurowski, A. Singh, D. Geiger, C. C. Liang, M. Y. Chiu, L. Adler, M. Haacke, D. L. Wilson. Cardiac MR image segmentation using deformable models. In Proceedings of Computers in Cardiology Conference, IEEE, London, UK, pp. 747–757, 1993. [18] S. Ranganath. Contour extraction from cardiac MRI studies using snakes. IEEE Transactions on Medical Imaging, vol.14, no.2, pp.328-338, 1995. DOI:10.1109/42.387714 [19] D. Geiger, A. Gupta, L. A. Costa, J. Vlontzos. Dynamic programming for detecting, tracking, and matching deformable contours. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol.17, no.3, pp.294-302, 1995. DOI:10.1109/34.368194 [20] C. M. Li, C. Y. Kao, J. C. Gore, Z. H. Ding. Implicit active contours driven by local binary fitting energy. In Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, IEEE, Minneapolis, USA, 2007. [21] S. Lankton, A. Tannenbaum. Localizing region-based active contours. IEEE Transactions on Image Processing, vol.17, no.11, pp.2029-2039, 2008. DOI:10.1109/TIP.2008.2004611 [22] L. Wang, L. He, A. Mishra, C. M. Li. Active contours driven by local Gaussian distribution fitting energy. Signal Processing, vol.89, no.12, pp.2435-2447, 2009. DOI:10.1016/j.sigpro.2009.03.014 [23] T. Chen, J. Babb, P. Kellman, L. Axel, D. Kim. Semiautomated segmentation of myocardial contours for fast strain analysis in cine displacement-encoded MRI. IEEE Transactions on Medical Imaging, vol.27, no.8, pp.1084-1094, 2008. DOI:10.1109/TMI.2008.918327 [24] H. Y. Lee, N. Codella, M. Cham, M. Prince, J. Weinsaft, Y. Wang. Left ventricle segmentation using graph searching on intensity and gradient and a priori knowledge (lvGIGA) for short-axis cardiac magnetic resonance imaging. Journal of Magnetic Resonance Imaging, vol.28, no.6, pp.1393-1401, 2008. DOI:10.1002/jmri.v28:6 [25] D. Mahapatra. Cardiac image segmentation from cine cardiac MRI using graph cuts and shape priors. Journal of Digital Imaging, vol.26, no.4, pp.721-730, 2013. DOI:10.1007/s10278-012-9548-5 [26] J. Cousty, L. Najman, M. Couprie, S. Clément-Guinaudeau, T. Goissen, J. Garot. Segmentation of 4D cardiac MRI: Automated method based on spatio-temporal watershed cuts. Image and Vision Computing, vol.28, no.8, pp.1229-1243, 2010. DOI:10.1016/j.imavis.2010.01.001 [27] S. P. Dakua, J. Abi-Nahed. Patient oriented graph-based image segmentation. Biomedical Signal Processing and Control, vol.8, no.3, pp.325-332, 2013. DOI:10.1016/j.bspc.2012.11.009 [28] A. Pednekar, U. Kurkure, R. Muthupillai, S. Flamm, I. A. Kakadiaris. Automated left ventricular segmentation in cardiac MRI. IEEE Transactions on Biomedical Engineering, vol.53, no.7, pp.1425-1428, 2006. DOI:10.1109/TBME.2006.873684 [29] H. F. Hu, H. H. Liu, Z. Y. Gao, L. Huang. Hybrid segmentation of left ventricle in cardiac MRI using Gaussian-mixture model and region restricted dynamic programming. Magnetic Resonance Imaging, vol.31, no.4, pp.575-584, 2013. DOI:10.1016/j.mri.2012.10.004 [30] H. F. Hu, Z. Y. Gao, L. M. Liu, H. H. Liu, J. F. Gao, S. Z. Xu, W. Li, L. Huang. Automatic segmentation of the left ventricle in cardiac MRI using local binary fitting model and dynamic programming techniques. PLoS One, vol. 9, no. 12, Article number e114760, 2014. DOI: 10.1371/journal.pone.0114760. [31] P. Radau Y. Lu, K. Connelly, G. Paul, A. J. Dick, G. A. Wright. Evaluation framework for algorithms segmenting short axis cardiac MRI. MIDAS Journal-Cardiac MR Left Ventricle Segmentation Challenge, London, UK, 2009. [Online], Available: http://hdl.handle.net/10380/3070, October 24, 2012. [32] L. G. Nyúl, J. K. Udupa, X. Zhang. New variants of a method of MRI scale standardization. IEEE Transactions on Medical Imaging, vol.19, no.2, pp.143-150, 2000. DOI:10.1109/42.836373 [33] T. F. Chan, L. A. Vese. Active contours without edges. IEEE Transactions on Image Processing, vol.10, no.2, pp.266-277, 2001. DOI:10.1109/83.902291 [34] H. L. Huang, X. Zuo, C. Huang. Arbitrary initialization for Chan-Vese model. Optik-International Journal for Light and Electron Optics, vol.125, no.18, pp.5257-5263, 2014. DOI:10.1016/j.ijleo.2014.06.051 [35] C. Constantinides, E. Roullot, M. Lefort, F. Frouin. Fully automated segmentation of the left ventricle applied to cine MR images: Description and results on a database of 45 subjects. In Proceedings of Annual International Conference of IEEE Engineering in Medicine and Biology Society, IEEE, San Diego, USA, pp. 3207–3210, 2012. [36] C. Constantinides, Y. Chenoune, N. Kachenoura, E. Roullot, E. Mousseaux, A. Herment, F. Frouin. Semi-automated cardiac segmentation on cine magnetic resonance images using GVF-snake deformable models. MIDAS Journal- Cardiac MR Left Ventricle Segmentation Challenge, London, UK, 2009. [Online], Available: http://hdl.handle.net/10380/3108, October 25, 2012. [37] H. F. Hu, H. H. Liu, Z. Y. Gao, L. Huang. Hybrid segmentation of left ventricle in cardiac MRI using Gaussian-mixture model and region restricted dynamic programming. Magnetic Resonance Imaging, vol.31, no.4, pp.575-584, 2013. DOI:10.1016/j.mri.2012.10.004 [38] S. Huang, J. M. Liu, L. C. Lee, S. K. Venkatesh, L. L. S. Teo, C. Au, W. L. Nowinski. An image-based comprehensive approach for automatic segmentation of left ventricle from cardiac short axis cine MR images. Journal of Digital Imaging, vol.24, no.4, pp.598-608, 2011. DOI:10.1007/s10278-010-9315-4 [39] M. P. Jolly, H. Xue, L. Grady, J. Guehring. Combining registration and minimum surfaces for the segmentation of the left ventricle in cardiac cine MR images. In Proceedings of the 12th International Conference on Medical Image Computing and Computer-Assisted Intervention, London, UK, pp. 910–918, 2009. [40] H. Liu, H. F. Hu, X. Y. Xu, E. M. Song. Automatic left ventricle segmentation in cardiac MRI using topological stable-state thresholding and region restricted dynamic programming. Academic Radiology, vol.19, no.6, pp.723-731, 2012. DOI:10.1016/j.acra.2012.02.011 [41] S. Queirós, D. Barbosa, B. Heyde, P. Morais, J. L. Vilaça, D. Friboulet, O. Bernard J. D′hooge. Fast automatic myocardial segmentation in 4d cine CMR datasets. Medical Image Analysis, vol.18, no.7, pp.1115-1131, 2014. DOI:10.1016/j.media.2014.06.001 [42] J. Schaerer, C. Casta, J. Pousin, P. Clarysse. A dynamic elastic model for segmentation and tracking of the heart in MR image sequences. Medical Image Analysis, vol.14, no.6, pp.738-749, 2010. DOI:10.1016/j.media.2010.05.009 [43] M. G. Uzunbas, S. T. Zhang, K. M. Pohl, D. Metaxas, L. Axel. Segmentation of myocardium using deformable regions and graph cuts. In Proceedings of the 9th IEEE International Symposium on Biomedical Imaging, IEEE, Barcelona, Spain, pp. 254–257, 2012. [44] J. Wijnhout, D. Hendriksen, H. Van Assen, R. Van der Geest. LV challenge LKEB contribution: fully automated myocardial contour detection. MIDAS Journal-Cardiac MR Left Ventricle Segmentation Challenge, London, UK, 2009. [Online]. Available: http://hdl.handle.net/10380/3115, October 25, 2012. [45] L. Marak, J. Cousty, L. Najman, H. Talbot. 4D morphological segmentation and the MICCAI LV-segmentation grand challenge. MIDAS Journal-Cardiac MR Left Ventricle Segmentation Challenge, London, UK, 2009. [Online], Available: http://hdl.handle.net/10380/3085, October 25, 2012. [46] C. L. Feng, C. M. Li, D. Z. Zhao, C. Davatzikos, H. Litt. Segmentation of the left ventricle using distance regularized two-layer level set approach. In Proceedings of the 16th International Conference on Medical Image Computing and Computer-Assisted Intervention, Nagoya, Japan. pp. 477–484, 2013. [47] X. H. Qian, Y. Lin, Y. Zhao, J. Wang, J. Liu, X. H. Zhuang. Segmentation of myocardium from cardiac MR images using a novel dynamic programming based segmentation method. Medical Physics, vol.42, no.3, pp.1424-1435, 2015. DOI:10.1118/1.4907993