1 Introduction
The World Health Organization (WHO)^{1}^{1}1http://http://www.who.int/mediacentre/factsheets/fs317/en/estimated 17.5 million deaths from cardiovascular diseases in 2012, representing of all mortalities, rendering cardiovascular conditions the main cause of death globally. Hence, the timely diagnosis and treatment followup of these pathologies is crucial. High image quality, good tissue contrast, and no ionizing radiation has established MRI as a standard clinical modality for noninvasive assessment of cardiac performance. Cardiac contractile function quantified via the systolic and diastolic volumes, ejection fraction, and myocardial mass represents a reliable diagnostic value and can be computed by segmenting the left (LV) and right (RV) ventricles from cardiac cineMRI. Although manual delineation of the ventricle is deemed as the goldstandard approach, it requires significant time and effort and is highly susceptible to inter and intraobserver variability. These limitations suggest a need for fast, robust, and accurate semi or fullyautomatic segmentation algorithms.
Various segmentation techniques for cardiac MR images have been proposed in the literature [1]. The imagebased approaches with weak or no prior information, such as thresholding, edgebased and regionbased approaches, or pixelbased classifications methods, require user interaction for proper segmentation of the illdefined regions. On the other hand, shape prior deformation models, active shape and appearance models, and atlasbased approaches are more likely to overcome this problem at the expense of manually building a training set.
Multiatlas based approaches have shown promising results in biomedical image segmentation [2]. However, they rely on a number of computationally demanding and time limiting nonrigid image registration steps followed by label fusion. Hence despite its accuracy, it has experienced minimal to no adoption in actual clinical applications primarily due to its complexity, high dependence on parameters variability, and computational demands.
On the other hand, combinatorial optimization based graphcut techniques are fast and guaranteed to produce results within a known factor of the global minimum, for some special classes of functions (termed as regular functions)
[3] and have proved to be powerful tools for image segmentation. Moreover, adding a shape constraint into the graph cut framework has been shown to improve the cardiac image segmentation results significantly [4, 5, 6]. However, these methods require a manual input to introduce a shape constraint at the right location in the image.In this work, we leverage the performance of the graph cut framework and augment it by incorporating shape constraints in the form of an average atlasbased segmentation of the anatomy whose label was generated and propagated using a single affine registration. Subsequently, we iteratively refine the segmentation using techniques similar to those described in [7, 8], to obtain an accurate and robust segmentation of the myocardium. Hence, we do not require any manual input to introduce shape constraint into the graphcut framework and simultaneously take advantage of the prior knowledge in the form of atlasbased segmentation requiring affine as opposed to nonrigid registration, which is more computationally efficient and less sensitive to parameters variability.
2 Methodology
Whole heart cineMRI images are generated by stacking 2D+T shortaxis slices acquired during a single breath hold. Since this acquisition approach introduces an intensity difference between the slices, as well as slice misalignments, we can follow one of two approaches to segment tha data: one approach is to implement a slice motion correction protocol to realign the slices into a coherent 3D volume. The other approach, also implemented here, resorts to slicewise processing and segmentation instead of a 3D segmentation.
Another challenge is the illdefined contrast of the LV myocardium in MR images, which makes the imagedriven segmentation difficult. As such, to obtain better segmentation of the apical and basal regions, we exploit the prior knowledge in the form of an average atlas. The proposed methodology formulates the segmentation problem in the context of a graph based energy minimization framework. The blood pool is first segmented using an iterative graph cut technique; then, this information is used to segment the myocardium.
2.1 Data Preprocessing
This study is conducted on 30 cardiac cineMR images taken from the DETERMINE [9] cohort available as a part of the STACOM Cardiac Atlas Segmentation Challenge Project database^{2}^{2}2http://www.cardiacatlas.org. The semiautomatically segmented images obtained by applying the method described in [10] accompany the dataset and serves as goldstandard for assessing the proposed segmentation technique.
We select a reference patient volume with good contrast, average size, and preferred LVRV orientation. All patient volumes are rotated about the zaxis (i.e., sliceencoding direction) to roughly align their orientation with that of the reference patient using the DICOM Image Orientation Patient (IPP) field. The region of interest (ROI) (in the xyplane) enclosing the left and right ventricles is extracted using the method described in [11] by correlating the 2D motion images generated from the 3D volumes across the cardiac cycle. The only manual input required by our algorithm is the start and end slices of the LV, such that the ROI is restricted in the zdirection, preventing over/under segmentation of slices that do not belong to the desired anatomy. The patient volumes are cropped to the above ROI, and, to compensate for any intensity differences (due to the slicewise acquisition), each slice is normalized (0255) prior to further processing.
2.2 Atlas Generation
The cropped 3D volumes (at the end diastole phase) for all patients are first histogram matched and then affinely registered to the reference patient image volume using the intensity based NelderMeade downhill simplex algorithm [12] available in SimpleITK. The resulting 3D affine transforms are applied to the respective ground truth segmentations. The transformed volumes and transformed ground truths are then averaged to obtain an average appearance atlas and a probabilistic atlas, respectively (Fig. 1).
The average appearance atlas is registered to a test volume using intensitybased affine registration. The resulting registration transformation is used to transform the myocardial probabilistic label to the test data, which, in turn, serves as a shape constraint for the graph cut framework.
2.3 LV Blood Pool Segmentation using Iterative Graph Cuts
To leverage the 3D LV geometry, we use the blood pool (BP) segmentation of a given slice to help refine the BP ROI in the neighboring slices. As such, we first segment the BP from the midslice, followed by its neighboring slices, and proceed accordingly, until the complete volume is segmented.
2.3.1 Intensity Distribution Model:
The myocardium probability map for each slice is normalized and inverted to produce the probability map corresponding to the blood pool (BP) and background (BG). The resulting BP/BG probability map is thresholded at 0.5 and the inner connected component is isolated to obtain the high confidence BP ROI. Otsu thresholding
[13]is applied within this ROI to obtain the initial BP region. The intensity values within this extracted BP region are then fitted to a Gaussian distribution to generate the BP intensity model.
A binary mask enclosing the myocardium is obtained by thresholding the myocardial probability map at a very small value (i.e. 0.1). Holes in the binary mask are filled to obtain a ROI enclosing the BP, myocardium, and BG. To generate the BG intensity model, we fit the intensity values within the ROI, excluding the initial BP region, to a Gaussian Mixture Model (GMM) comprising two Gaussians.
Fig. 2a shows the resulting BP loglikelihood map.Note that we propose the Gaussian distribution for modeling intensity noise in MR images instead of a more appropriate Rician distribution [14]; this simplifies our model and is a good approximation when the signaltonoise ratio is high.
2.3.2 Blood Pool/Background Probabilistic Map:
To obtain a ROI that includes the myocardium and BP, we threshold the myocardial probability map at 0.5, fill in the blood pool, and erode the resulting ROI by (selected empirically) of the radius of its smallest circumscribed circle to obtain the BP ROI. The BP/BG probability map masked by the BPROI represents the BP probability map, and its inverse represents the BG probability map. Fig. 2b shows the BP loglikelihood map.
2.3.3 GraphCut Segmentation:
We construct a graph with each node (i.e., pixel) connected to its east, west, north, and south neighbors. Two special terminal nodes representing two classes — the source (blood pool), and the sink (background) — are added to the graph and all other nodes are connected to each terminal node. The segmentation is formulated as an energy minimization problem over the space of optimal labelings :
(1) 
where the first term represents the data energy that reduces the disagreement between the labeling given the observed data at every pixel , and the second term represents the smoothness energy that forces pixels and defined by a set of interacting pair (in our case, the neighboring pixels) towards the same label.
The data energy term is represented by the terminal link (tlink) between each node and the source (or sink), which is defined as the weighted sum of the log probabilities of the intensity distribution model and the probabilistic map corresponding to the BP (or BG):
(2) 
where, is the iteration number, is the likelihood of observing the intensity given that pixel belongs to class , and
is the prior probability for class
obtained from the BP/BG probability map. The loglikelihood difference between BP and BG labels for is shown in Fig. 2c. The intensity likelihood term (first term) allows the expansion of the BP region in the first few iterations, whereas the prior probability map (second term) restricts its “spilling” (due to oversegmentation) in subsequent iterations.The smoothness energy term is computed over the links between neighboring nodes (nlinks), which are weighted based on their intensity similarity:
(3) 
where is the pixel intensity. To avoid the “spilling” of the BP into the myocardium or BG, the smoothness term changes with each iteration, such that, in order for the neighboring pixels to be assigned to the same label during the current iteration, their intensities must be closer than in the previous iteration.
Once weights are assigned to all edges in the graph, the minimum cut equivalent to the maximum flow is identified via the expansion algorithm described in [15]. This approach yields the labeling (graphcut) that minimizes the global energy of the graph that corresponds to the optimal segmentation (Fig. 2d). Lastly, the convex hull applied to the graphcut result constitutes the final BP segmentation, such that, the papillary muscles are included within the BP (Fig. 2e).
2.3.4 Myocardial Probability Map Refinement
The myocardial probability map is thresholded at 0.5, and the inner hollow circular region representing the BP is extracted. The signed distance map corresponding to the boundary of the extracted BP region is affinely registered to the signed distance map generated from the boundary of the graphcut extracted BP (2.3) segmentation. The optimum affine transformation that minimizes the sum of squared differences between the two distance maps is applied to the myocardial probability map, such that, it fits the shape of the segmented BP.
2.3.5 Iterative Refinement:
The latest BP segmentation obtained from the graph cut is used to update the intensity distribution model. The refined myocardial probability map is used to construct a new BP/BG probability map. The pixels within the latest BP segmentation are assigned very high likelihood (for belonging to the BP), and hence their labels do not change. An updated BP segmentation is obtained via another graph cut operating on the new graph energy configuration. This iterative process is repeated until the changes in the affine transform parameters for the myocardium probability map are below a predefined threshold; this iterative process usually converges within three iterations. Upon convergence, the convex hull defined by the latest segmentation result constitutes the final BP segmentation. Fig. 3 illustrates the iterative refinement process.
2.4 Myocardium Segmentation
The information from the BP segmentation along with the refined myocardial probability map is used to segment the myocardium.
2.4.1 Intensity Distribution Model:
We select a ROI in each slice based on the refined probability map, and we match the histogram of the pixel intensities within this ROI to the histogram of the midslice. We select the midslices (i.e. no apical/basal slices) to obtain a single intensity distribution model for the whole volume. The intensities of the pixels within the refined myocardial mask with probability higher than 0.5 are fitted to a single Gaussian GMM to obtain the myocardium intensity distribution model. Similarly, the intensities of the remaining pixels are fitted to a three Gaussian GMM to obtain the BG intensity distribution model. Fig. 4a shows the loglikelihood map for the myocardium.
2.4.2 Distance from the Endocardial Border:
The endocardial border is obtained from the outer edge of the final BP segmentation (Fig. 2e). The knowledge that myocardium should be closer to the endocardial border is encoded in the data term represented by the truncated distance map (empirically selected as 10 pixels). This constraint increases the likelihood of pixels near the endocardial border to be labeled as myocardium, while reducing this likelihood for the pixels located further away. Furthermore, to prevent the BP region from being labeled as myocardium, it is assigned the lowest likelihood value (Fig. 4c).
2.4.3 GraphCut Segmentation:
A graph is constructed similar to the formulation described in 2.3.3
, but this time to classify the myocardial rather than blood pool pixels. The data term is defined as the weighted sum of the intensity distribution model, refined myocardial probability map (as described in
2.3.4 and Fig. 4b), and the distance from endocardial border, with increasing relative influence, respectively. The smoothness term varies spatially according to the intensity difference between the neighboring pixes, as discussed in 2.3.3. The minimum cut in the graph yields the final myocardium segmentation (Fig. 4e).3 Results
The proposed algorithm was implemented in Python and required 45 seconds on average to segment the BP and myocardium from cine MRI volumes on an Intel Xenon 3.60 GHz 32GB RAM PC.
Adhering to the collated results reported for the LV segmentation challenge in [16]
, we evaluated our segmentation on 30 patient datasets according to the following metrics: dice index, jaccard index, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV)
[17]. To maintain approximately equal number of myocardium and nonmyocardium pixels for evaluation, such that the NPV conveys some useful information, we dilated each slice of the myocardium region, for the provided gold standard segmentation, by one fourth of the radius of the disk with equivalent area. The segmentation results for a patient dataset are overlaid onto each slice of the patient volume and shown in Fig. 5a. Fig. 5b shows a visual comparison of our segmentation results visàvis the provided semiautomated segmentation serving as a goldstandard. The metrics are summarized in Table 1 for all slices together, as well as for the midslices and apical/basal slices (first and last two slices, respectively) separately.Assessment Metric  MidSlices  Apical/BasalSlices  All Slices 

Dice Index  
Jaccard Index  
Sensitivity  
Specificity  
PPV  
NPV 
4 Discussion, Conclusion, and Future Work
Our validation experiments show that the overall segmentation results are comparable to those reported in [17]. Specifically, the mean values for reported assessment metrics were: dice index — 0.68 to 0.88, jaccard index — 0.53 to 0.80, sensitivity — 0.63 to 0.90, specificity — 0.73 to 0.99, PPV — 0.66 to 0.96, and NPV — 0.81 to 0.94. However, it should be noted that the metrics reported in [17] were evaluated against the consensus segmentation estimated based on the participating seven raters (manual and automatic) obtained using the STAPLE algorithm on 18 test patient datasets, whereas our reported metrics were assessed against the provided semiautomated gold standard segmentation of 30 training patient datasets. Hence, the metrics provide only an approximate estimate of our algorithm’s performance compared to the ones that participated in the challenge. Moreover, the average segmentation time of 45 secs per volume achieved using an unoptimized implementation in Python shows great potential for our algorithm to achieve, following minimal optimization, the near realtime performance demanded by some of the intended clinical applications.
Since the BP region in the midslices are better defined than in the apical/basal slices, the segmentation results are consistently better for the midslices. We also observed that the slicewise processing and iterative refinement might compromise the segmentation of the apical/basal slices due to illdefined BP regions, suggesting the need for special processing for these slices.
As part of our future work, we plan to automate the ROI detection in zdirection to eliminate the manual input required by our algorithm. In addition, instead of using a constant truncating endocardial distance constraint, we plan to use imagederived edge information to enable spatially varying truncating distances to improve the myocardium segmentation. Similarly, we will study the effect of selecting different thresholds for the probability maps, weight variability on the likelihood terms and, in turn, on the final myocardium segmentation. Lastly, we plan to extend the work and evaluate the segmentation performance on all 100 patient datasets and report performance according to the metrics outlined above.
References
 [1] Petitjean, C., Dacher, J.N.: A review of segmentation methods in short axis cardiac {MR} images. Medical Image Analysis 15(2) (2011) 169 – 184
 [2] Iglesias, J.E., Sabuncu, M.R.: Multiatlas segmentation of biomedical images: A survey. Medical Image Analysis 24(1) (2015) 205 – 219
 [3] Kolmogorov, V., Zabin, R.: What energy functions can be minimized via graph cuts? IEEE Transactions on Pattern Analysis and Machine Intelligence 26(2) (2004) 147–159
 [4] Grosgeorge, D., Petitjean, C., Dacher, J.N., Ruan, S.: Graph cut segmentation with a statistical shape model in cardiac MRI. Computer Vision and Image Understanding 117(9) (2013) 1027–1035

[5]
Freedman, D., Zhang, T.:
Interactive graph cut based segmentation with shape priors.
Proceedings  2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, CVPR 2005
I (2005) 755–762  [6] Mahapatra, D.: Cardiac Image Segmentation from Cine Cardiac MRI Using Graph Cuts and Shape Priors. Journal of Digital Imaging 26(4) (2013) 721–730
 [7] Slabaugh, G., Unal, G.: Graph cuts segmentation using an elliptical shape prior. Proceedings  International Conference on Image Processing, ICIP 2 (2005) 1222–1225
 [8] Vu, N., Manjunath, B.S.: Shape prior segmentation of multiple objects with graph cuts. 26th IEEE Conference on Computer Vision and Pattern Recognition, CVPR (2008)
 [9] Kadish, A.H., Bello, D., Finn, J.P., Bonow, R.O., Schaechter, A., Subacius, H., Albert, C., Daubert, J.P., Fonseca, C.G., Goldberger, J.J.: Rationale and design for the defibrillators to reduce risk by magnetic resonance imaging evaluation (determine) trial. Journal of Cardiovascular Electrophysiology 20(9) (2009) 982–987
 [10] Li, B., Liu, Y., Occleshaw, C.J., Cowan, B.R., Young, A.A.: Inline automated tracking for ventricular function with magnetic resonance imaging. JACC: Cardiovascular Imaging 3(8) (2010) 860 – 866
 [11] BenZikri, Y.K., Linte, C.A.: A robust automated left ventricle region of interest localization technique using a cardiac cine mri atlas (2016)
 [12] Nelder, J.A., Mead, R.: A simplex method for function minimization. The Computer Journal 7(4) (1965) 308–313
 [13] Otsu, N.: A threshold selection method from graylevel histograms. IEEE Transactions on Systems, Man, and Cybernetics 9(1) (1979) 62–66
 [14] Gudbjartsson, H., Patz, S.: The rician distribution of noisy mri data. Magnetic Resonance in Medicine 34(6) (1995) 910–914
 [15] Boykov, Y., Veksler, O., Zabih, R.: Fast approximate energy minimization via graph cuts. IEEE Trans PAMI 23 (2001) 1222–39
 [16] Suinesiaputra, A., Cowan, B.R., AlAgamy, A.O., Elattar, M.A., Ayache, N., Fahmy, A.S., Khalifa, A.M., MedranoGracia, P., Jolly, M.P., Kadish, A.H., Lee, D.C., Margeta, J., Warfield, S.K., Young, A.A.: A collaborative resource to build consensus for automated left ventricular segmentation of cardiac {MR} images. Medical Image Analysis 18(1) (2014) 50 – 62
 [17] Suinesiaputra, A., Cowan, B.R., Finn, J.P., Fonseca, C.G., Kadish, A.H., Lee, D.C., MedranoGracia, P., Warfield, S.K., Tao, W., Young, A.A. In: Left Ventricular Segmentation Challenge from Cardiac MRI: A Collation Study. Springer Berlin Heidelberg, Berlin, Heidelberg (2012) 88–97
Comments
There are no comments yet.