1. Introduction

Cine MRI is a suitable modality in analyzing cardiac functions by segmenting the endo- and epicardial contours of the left ventricle (LV) [1]. Given the unclear boundaries between these structures and against the background, prior geometric knowledge needs to be incorporated. A typical idea is to maintain a close-to-constant distance between the endo- and epicardial contours (a “ring” structure) with either hard values [2] or by a soft constraint [3]. However, such static/pre-defined distance may fail to capture large LV shape variations across slices and phases among different subjects. More recently, it has been proposed to use a two-level contours of a single level set function to represent the endo- and epicardium [4], with their enclosing relations maintained by a smoothly varying distance. A local intensity model has been utilized to handle the overlapped distributions. However, representing the endo- and epicardium with a single level set function may hinder flexible incorporation of shape prior specific to individual structures [5,6,7]. The local intensity model could also be overly sensitive to contour initializations and kernel parameter settings [8].

In this study, we propose to represent the endo- and epicardium with two separate level set functions, but introduce a novel coupling setup embedded in the length regularization, with a symmetric “dilation/erosion” morphological relation imposed, to achieve a smoothly varying distance between the endo- and epicardial contours. The coupled level set representations further allows us to incorporate a sparse composite shape descriptor to boost the performance. Gaussian mixture models are utilized to represent overlapped intensity distributions in each subregion. Compared to local intensity models, it is more robust to contour initializations, noise, and selection of setup parameters.

2. Method

2.1. Variational segmentation formulation with a coupled length regularization and sparse composite shape prior

Figure 1 illustrates the proposed cardiac region separation scheme. Specifically, given an image $I:\Omega \to \mathbb{R}$ and two dependent level set functions ${{\phi }_{1}}, {{\phi }_{2}}:\Omega \to \mathbb{R}$, we partition $\Omega$ into three subdomains: ${\Omega _1} \buildrel \Delta \over = \left\{ {{\rm{X }}|{\rm{ }}{\phi _1}\left( {\rm{X}} \right) \le 0} \right\}$, ${\Omega _1} \buildrel \Delta \over = \left\{ {{\rm{X }}|{\rm{ }}{\phi _1}\left( {\rm{X}} \right) \le 0} \right\}$, ${\Omega _3} \buildrel \Delta \over = \left\{ {{\rm{X }}|{\rm{ }}{\phi _2}\left( {\rm{X}} \right) > 0} \right\}$, representing the left ventricular cavity, myocardium and background, respectively. A typical variational segmentation formulation consists of two components: E($\Phi$) = E_fidelity + E_reg, a data fidelity and a regularization, which we will elaborate bellow.

2.1.1. E_fidelity

We choose to model intensity distributions with Gaussian mixtures for each ${{\Omega }_{i}}$ [9]:

$$P\left( {I\left( {\rm{X}} \right)|{\Omega _i},{\theta _i}} \right) = \sum\limits_{j = 1}^{{n_i}} {\frac{{{w_{j,i}}}}{{\sqrt {2\pi \sigma _{j,i}^2} }}} \exp ( - \frac{{{{\left( {I\left( {\rm{X}} \right) - {\mu _{j,i}}} \right)}^2}}}{{2\sigma _{j,i}^2}}),$$

(1)

where n_i represents the number of Gaussian components for each ${{\Omega }_{i}}$, ${{\theta }_{i}}=\left\{ {{w}_{j, i}}, {{\mu }_{j, i}}, {{\sigma }_{j, i}} \right\}$ represents the weights and parameters for the jth Gaussian component. E_fidelity is then constructed as a weighted log-likelihood as:

$$\begin{align} & {{E}_{fidelity}}\left( \Phi , \Theta \right)=-\int_{\Omega }{\left( 1-H\left( {{\phi }_{1}} \right) \right)}\log \left( {{c}_{1}}P\left( I\left( \text{x} \right)|{{\theta }_{1}} \right) \right)d\text{x} \\ & \qquad \qquad \qquad \quad -\int_{\Omega }{H\left( {{\phi }_{1}} \right)\left( 1-H\left( {{\phi }_{2}} \right) \right)}\log \left( {{c}_{2}}P\left( I\left( \text{x} \right)|{{\theta }_{2}} \right) \right)d\text{x} \\ & \qquad \qquad \qquad \quad -\int_{\Omega }{H\left( {{\phi }_{2}} \right)}\log \left( {{c}_{3}}P\left( I\left( \text{x} \right)|{{\theta }_{3}} \right) \right)d\text{x, } \end{align}$$

(2)

where H represents the heaviside function, c₁, c₂, c₃ weights the importance/gains of correctly segmenting each corresponding regions.

2.1.2. E_reg

Our design of regularization consists of two components: E_length and E_shape.

E_length: For a smoothly varying distance between the endo- and epicardial contours at each location, we introduce a novel coupling in the length regularization with

$${{E}_{length}}=\int_{\Omega }{{{g}_{1}}\delta \left( {{\phi }_{1}} \right)|\nabla {{\phi }_{1}}|d\text{x}}+\int_{\Omega }{{{g}_{2}}\delta \left( {{\phi }_{2}} \right)|\nabla {{\phi }_{2}}|d\text{x}}+\frac{\beta }{2}\int_{\Omega }{|\nabla d{{|}^{2}}d\text{x}},$$

(3)

where $\delta$ represents the Dirac delta function, ${{g}_{1}}={{\left( {{\phi }_{2}}\left( \text{x} \right)+d\left( \text{x} \right) \right)}^{2}}$ and ${{g}_{2}}={{\left( {{\phi }_{1}}\left( \text{x} \right)-d\left( \text{x} \right) \right)}^{2}}$ are the dilated/eroded geometry indicator functions that traps ${{\phi }_{1}}$ to be the distance d from ${{\phi }_{2}}$ and vice versa.

E_shape: Inspired by a recently developed shape regularization [7], we construct two shape library ${D_{endo}} = \left[ {\psi _1^{endo},\psi _2^{endo}, \cdots \psi _m^{endo}} \right]$ and ${D_{epi}} = \left[ {\psi _1^{epi},\psi _2^{epi}, \cdots \psi _m^{epi}} \right]$, representing the training shapes. The shape regularization is designed as:

$${{E}_{shape}}=\int_{\Omega }{\left\{ {{\left( {{\phi }_{1}}-{{D}_{endo}}\text{w} \right)}^{2}}+{{\left( {{\phi }_{2}}-{{D}_{epi}}\text{w} \right)}^{2}} \right\}}d\text{x+}\gamma {{\left\| \text{w} \right\|}_{1}},$$

(4)

which dynamically regularizes the current estimates of ${{\phi }_{1}}$, ${{\phi }_{2}}$ towards a sparse linear combinations of the training shapes weighted by w. In this study, we utilize the single weight w to ensure the coupling of training shape selections, since the corresponding column/shape in D_endo and D_epi comes from the same patient. However, our regularization also permits a more general form that has different library size with different weights (w_endo and w_epi) under the situation when training shapes do not pair up or come from heterogeneous sources.

With the introduced E_length and E_shape, our E_reg is designed as the following:

$$\begin{align} & {{E}_{reg}}\left( \Phi , d, \text{w} \right)=\frac{{{\lambda }_{1}}}{2}\int_{\Omega }{{{\left( {{\phi }_{2}}+d \right)}^{2}}\delta \left( {{\phi }_{1}} \right)|\nabla {{\phi }_{1}}|d\text{x}+ \frac{{{\lambda }_{2}}}{2}\int_{\Omega }{{{\left( {{\phi }_{1}}-d \right)}^{2}}\delta \left( {{\phi }_{2}} \right)|\nabla {{\phi }_{2}}|d\text{x}+\frac{\beta }{2}\int_{\Omega }{|\nabla d{{|}^{2}}d\text{x}}}} \\ & \qquad \qquad \qquad +\frac{\alpha }{2}\int_{\Omega }{\left\{ {{\left( {{\phi }_{1}}-{{D}_{endo}}\text{w} \right)}^{2}}+{{\left( {{\phi }_{2}}-{{D}_{epi}}\text{w} \right)}^{2}} \right\}d\text{x}+}\gamma {{\left\| \text{w} \right\|}_{1}}+\frac{\mu }{2}\sum\limits_{i=1}^{2}{\int_{\Omega }{{{\left( \nabla {{\phi }_{i}}-1 \right)}^{2}}d\text{x, }}} \end{align}$$

(5)

where $\sum _{i = 1}^2\int_\Omega {{{\left( {\nabla {\phi _i} - 1} \right)}^2}} d{\rm{x}}$ is designed to preserve the signed distance property of each level set functions [10]. Smooth approximations of the heaviside and Delta functions are used for numerical implementation [11]. Algorithm 1 presents the block coordinate descent/minimization scheme used to minimize the variational energy.

ALGORITHM 1 Minimization scheme
1:	while |Ek - Ek-1| < tol do
2:	minimization w.r.t. ${{\theta }_{1}}$, ${{\theta }_{2}}$, ${{\theta }_{3}}$ with Expectation-Maximization (EM) algorithm [12]:    $\theta _{1}^{k+1}, \theta _{2}^{k+1}, \theta _{3}^{k+1}\leftarrow EM\left( \phi _{1}^{k}, \phi _{2}^{k}, I, {{\text{w}}^{k}}, {{d}^{k}} \right)$
3:	minimization w.r.t. w with alternating direction method of multipliers (ADMM) [13]:    ${{\text{w}}^{k+1}}\leftarrow ADMM\left( \frac{\alpha }{2}\left( \left\| \phi _{1}^{k}-{{D}_{endo}}\left. \text{w} \right\| \right._{2}^{2}+\left\| \phi _{2}^{k} \right.-{{D}_{epi}}\left. \text{w} \right\|_{2}^{2} \right)+\gamma {{\left\| \text{w} \right\|}_{1}} \right)$
4:	gradient descending w.r.t. d    $d_{1}^{k+1}\leftarrow d_{iter\_num}^{k}$ while m < iter num do    $d_{m+1}^{k+1}=d_{m}^{k+1}+\Delta t\left\{ -{{\lambda }_{1}}\left( \phi _{2}^{k}+d_{m}^{k+1} \right)\delta \left( \phi _{1}^{k} \right)|\nabla \phi _{1}^{k}|-{{\lambda }_{2}}\left( \phi _{1}^{k}-d_{m}^{k+1} \right)\delta \left( \phi _{2}^{k} \right)|\nabla \phi _{2}^{k}|+\beta \Delta d_{m}^{k+1} \right\}$    end while
5:	gradient descending w.r.t. ${{\phi }_{1}}$, ${{\phi }_{2}}$    $\phi _{1, 1}^{k+1}\leftarrow \phi _{1, iter\_num}^{k}, \phi _{2, 1}^{k+1}\leftarrow \phi _{2, iter\_num}^{k}$    while i < iter_num do $\begin{align} & \phi _{1, i+1}^{k+1}=\phi _{1, i}^{k+1}+\Delta t\{-\delta \left( \phi _{1, i}^{k+1} \right)\log \left( {{c}_{1}}P\left( I\left( \text{x} \right)|\theta _{1}^{k+1} \right) \right)+\delta \left( \phi _{1, i}^{k+1} \right)\left( 1-H\left( \phi _{2, i}^{k+1} \right) \right)\log \left( {{c}_{2}}P\left( I\left( \text{x} \right)|\theta _{2}^{k+1} \right) \right) \\ & \qquad \quad +\frac{{{\lambda }_{1}}}{2}\delta \left( \phi _{1, i}^{k+1} \right)\nabla \cdot \left( {{\left( \phi _{2, i}^{k+1}+d_{m+1}^{k+1} \right)}^{2}}\frac{\nabla \phi _{1, i}^{k+1}}{\left| \nabla \phi _{1, i}^{k+1} \right|} \right)-{{\lambda }_{2}}\left( \phi _{1, i}^{k+1}-d_{m+1}^{k+1} \right)\delta \left( \phi _{2, 1}^{k+1} \right)|\nabla \phi _{2, 1}^{k+1}| \\ & \qquad \quad -\alpha \left( \phi _{1, 1}^{k+1}-{{D}_{endo}}{{\text{w}}^{k+1}} \right)+\mu \left( \phi _{1, 1}^{k+1}-\nabla \cdot \left( \frac{\nabla \phi _{1, i}^{k+1}}{\left| \nabla \phi _{1, i}^{k+1} \right|} \right) \right)\} \\ \end{align}$ $\begin{align} & \phi _{2, i+1}^{k+1}=\phi _{2, i}^{k+1}+\Delta t\{-\delta \left( \phi _{2, i}^{k+1} \right)H\left( \phi _{1, i+1}^{k+1} \right)\log \left( {{c}_{2}}P\left( I\left( \text{x} \right)|\theta _{2}^{k+1} \right) \right)+\delta \left( \phi _{2, i}^{k+1} \right)\log \left( {{c}_{3}}P\left( I\left( \text{x} \right)|\theta _{3}^{k+1} \right) \right) \\ & \qquad \quad -{{\lambda }_{1}}\left( \phi _{2, i}^{k+1}+d_{m+1}^{k+1} \right)\delta \left( \phi _{1, i+1}^{k+1} \right)\left| \nabla \phi _{1, i}^{k+1} \right|+\frac{{{\lambda }_{2}}}{2}\delta \left( \phi _{2, i}^{k+1} \right)\nabla \cdot \left( {{\left( \phi _{1, i+1}^{k+1}-d_{m+1}^{k+1} \right)}^{2}}\frac{\nabla \phi _{2, i}^{k+1}}{\left| \nabla \phi _{2, i}^{k+1} \right|} \right) \\ & \qquad \quad -\alpha \left( \phi _{2, 1}^{k+1}-{{D}_{epi}}{{\text{w}}^{k+1}} \right)+\mu \left( \phi _{2, 1}^{k+1}-\nabla \cdot \left( \frac{\nabla \phi _{2, i}^{k+1}}{\left| \nabla \phi _{2, i}^{k+1} \right|} \right) \right)\} \\ \end{align}$    end while
6:	end while

2.2. Parameter settings

For all experiments in this study, we manually chose and fixed the parameters as: c₁ = c₂ = 1, c₃ = 0.1, ${{\lambda }_{1}}$ = ${{\lambda }_{2}}$ = 2, β = 2, α = γ = 0.02, μ = 1 and time step $\Delta t$ = 0.1. We set the Gaussian component number to be n₁ = n₂ = 1 for the cavity and myocardium, and n₃ = 2 for the background. Our algorithm was implemented with matlab 2012b. It took about one minute to segment one image volume on a macbook pro laptop with quad-core i7 2.3GHz and 8GB RAM. No specific code optimization and/or paralellization schemes were used.

3. Experimental results

We evaluated the proposed method on datasets from MICCAI LV segmentation challenge (http://smial.sri.utoronto.ca/LV_Challenge/Home.html), which consists of 15 training datasets and 15 validation datasets from patients of different pathologies. The Cine steady state free precession MR short axis images were obtained with a 1.5T GE Signa MRI, with FOV = 320 × 320 mm², image size = 256 × 256 and slice-thickness = 8 mm. The manual segmented endo- and epicardial contours were provided in all slices at end-diastolic (ED) and end-systolic (ES) phases. Our training endo- and epicardium libraries were built upon 15 training datasets, and the evaluation was performed on the 15 validation datasets.

The segmentation accuracy was evaluated based on both Dice similarity coefficient (DSC) metric and average perpendicular distance (APD). DSC is defined as $DSC \buildrel \Delta \over = \frac{{2\left| {{C_{seg}} \cap {C_{truth}}} \right|}}{{\left| {{C_{seg}}} \right| + \left| {{C_{truth}}} \right|}}$, where C_seg and C_truth are the segmented foreground from the obtained and manual segmentation respectively. APD measures the perpendicular distance from the segmented contour to the corresponding manually contour, averaged over all contour points.

Figure 2 compares segmented contours from our method with manual segmented contours on one typical LV example. Figure 3 demonstrates the flexibility and robustness of the proposed method on segmenting one LV example with various thickness of its myocardium within and across different slices. Table 1 reports and compares the DSC and APD statistics of the proposed method with results from other state-of-the-art methods. The proposed method had slightly lower APD compared to [4]. It achieved competitive/better DSC and APD accuracy compared to the other state-of-the-art methods.

4. Discussion and conclusion

Our proposed coupled length/distance regularization dynamically updates and encourages a smoothly varying distance between the endo- and epicardium through the optimization process. It is more flexible and accurate to capture the large endocardium-epicardium variations across different slices and among different patients, as compared to conventional approaches where such distance is usually static and pre-determined by heuristics [2,3,16].

The proposed method is similar to the dual-layer single level set approach [4] in regularizing a smoothly varying distance between the endo- and epicardium, but with added flexibility of incorporating individual priors for each structure separately [5,6,7], such as via coupled or decoupled sparse composite shape models depends on the training source. In this study, we utilize a single weight w to ensure the coupling of epi- and endocardium selections, since those shapes come from a single patient in each shape library. However, under situations that training shapes do not pair up, our formulation would be flexible enough to provide an asymmetric and decoupled regularization. For example, when endo- and epicardial training shapes come from different groups of patients, our model is capable of providing an asymmetric regularization by enforcing different weights (w_endo and w_epi) from such different training sources to leverage all training instances. Such flexibility also permits further extension of the current model to a joint left and right ventricle (RV) segmentation framework, where training LV and RV shapes often come from heterogeneous sources. In the specific context where differentiating intensity gradient is often absent at the structure interface, the more global Gaussian mixture model we use has the advantage of being more stable, more robust to contour initialization and noise, and does not require the tuning of (heterogeneous) local kernel parameters.

In summary, we have proposed a novel variational segmentation method by introducing geometrical coupling into the length regularization and a flexible coupled/decoupled sparse composite shape prior. Our method has been evaluated on the validation datasets from MICCAI challenge and achieved comparable/better DSC and APD accuracy compared to other state-of-the-art methods. A joint left and right ventricle segmentation method is under development.

Conflict of Interest

All authors declare no conflicts of interest in this paper.