Segmentation of myocardium from cardiac MR images using a novel dynamic programming based segmentation method Xiaohua Qian, Yuan Lin, Yue Zhao, Jing Wang, Jing Liu, and Xiahai Zhuang Citation: Medical Physics 42, 1424 (2015); doi: 10.1118/1.4907993 View online: http://dx.doi.org/10.1118/1.4907993 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/42/3?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in Fully automated segmentation of cartilage from the MR images of knee using a multi-atlas and local structural analysis method Med. Phys. 41, 092303 (2014); 10.1118/1.4893533 Cavity contour segmentation in chest radiographs using supervised learning and dynamic programming Med. Phys. 41, 071912 (2014); 10.1118/1.4881096 Myocardial borders segmentation from cine MR images using bidirectional coupled parametric deformable models Med. Phys. 40, 092302 (2013); 10.1118/1.4817478 Improved dynamic-programming-based algorithms for segmentation of masses in mammograms Med. Phys. 34, 4256 (2007); 10.1118/1.2791034 A new 2D segmentation method based on dynamic programming applied to computer aided detection in mammography Med. Phys. 31, 958 (2004); 10.1118/1.1688039
Segmentation of myocardium from cardiac MR images using a novel dynamic programming based segmentation method Xiaohua Qian SJTUCU International Cooperative Research Center, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China and Shanghai Advanced Research Institute, Chinese Academy of Sciences, Shanghai 201210, China
Yuan Lin Division of Research and Innovations, Carestream Health, Inc., Rochester, New York 14615
Yue Zhao College of Electronic Science and Engineering, Jilin University, 2699 Qianjing Street, Changchun, Jilin 130012, China
Jing Wang Suzhou Institute of Biomedical Engineering and Technology Chinese Academy of Sciences, Suzhou, Jiangsu 215163, China
Jing Liu Shanghai Advanced Research Institute, Chinese Academy of Sciences, Shanghai 201210, China
Xiahai Zhuanga) SJTUCU International Cooperative Research Center, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
(Received 2 September 2014; revised 2 December 2014; accepted for publication 17 January 2015; published 27 February 2015) Purpose: Myocardium segmentation in cardiac magnetic resonance (MR) images plays a vital role in clinical diagnosis of the cardiovascular diseases. Because of the low contrast and large variation in intensity and shapes, myocardium segmentation has been a challenging task. A dynamic programming (DP) based segmentation method, incorporating the likelihood and shape information of the myocardium, is developed for segmenting myocardium in cardiac MR images. Methods: The endocardium, i.e., the left ventricle blood cavity, is segmented for initialization, and then the optimal epicardium contour is determined using the polar-transformed image and DP scheme. In the DP segmentation scheme, three techniques are proposed to improve the segmentation performance: (1) the likelihood image of the myocardium is constructed to define the external cost in the DP, thus the cost function incorporates prior probability estimation; (2) the adaptive search range is introduced to determine the polar-transformed image, thereby excluding irrelevant tissues; (3) the connectivity constrained DP algorithm is developed to obtain an optimal closed contour. Four metrics, including the Dice metric (Dice), root mean squared error (RMSE), reliability, and correlation coefficient, are used to assess the segmentation accuracy. The authors evaluated the performance of the proposed method on a private dataset and the MICCAI 2009 challenge dataset. The authors also explored the effects of the three new techniques of the DP scheme in the proposed method. Results: For the qualitative evaluation, the segmentation results of the proposed method were clinically acceptable. For the quantitative evaluation, the mean (Dice) for the endocardium and epicardium was 0.892 and 0.927, respectively; the mean RMSE was 2.30 mm for the endocardium and 2.39 mm for the epicardium. In addition, the three new techniques in the proposed DP scheme, i.e., the likelihood image of the myocardium, the adaptive search range, and the connectivity constrained DP algorithm, improved the segmentation performance for the epicardium with 0.029, 0.047, and 0.007 in terms of the Dice and 0.98, 1.31, and 0.21 mm in terms of the RMSE, respectively. Conclusions: The three techniques (the likelihood image of the myocardium, the adaptive search range, and the connectivity constrained DP algorithm) can improve the segmentation ability of the DP method, and the proposed method with these techniques has the ability to achieve the acceptable segmentation result of the myocardium in cardiac MR images. Therefore, the proposed method would be useful in clinical diagnosis of the cardiovascular diseases. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4907993] Key words: cardiac MRI, myocardium segmentation, likelihood image, adaptive search range, dynamic programming
1424
Med. Phys. 42 (3), March 2015
0094-2405/2015/42(3)/1424/12/$30.00
© 2015 Am. Assoc. Phys. Med.
1424
1425
Qian et al.: Segmentation of myocardium from cardiac MR
1. INTRODUCTION Cardiovascular diseases are the leading causes of deaths in the modern countries.1,2 Medical imaging technologies, such as the magnetic resonance imaging (MRI), computed tomography, echography, and nuclear medicine, are effective in reducing the mortality and morbidity of cardiovascular diseases. MRI provides information of the morphology and pathology of the heart and has emerged as a widely used noninvasive modality in clinical applications of cardiology.3 Extracting this information, for computing functional indices such as ejection fraction and wall thickness, relies on accurate delineation of the endocardial and epicardial boundaries of the left ventricle (LV). Manual delineation is labor-intensive and suffers from both inter- and intraexpert variability, therefore, automatic or semiautomatic algorithms become increasingly desirable. However, there are a number of challenges inherently from the cardiac MRI.4,5 First, due to the blood flow and partial volume effects, intensity patterns of the cardiac MRI vary across different scans. Second, structural shapes of the heart alter because of the different patients and different conditions. Moreover, the muscle tissues (with similar intensity) can be classified into anatomically different substructures, such as the myocardium, trabeculae carneae, and papillary muscles. Additionally, the myocardium and its neighboring tissues such as the liver and lung generally have the comparable intensity, due to the similar T2/T1 ratios. In recent years, a large number of segmentation approaches have been developed for the delineation of the LV endoand epicardium contours in cardiac MRI,5 including the approaches without prior6–9 and with prior.10–24 Methods without prior includes the threshold,6 clustering,7 snake/active contour model,8,9 etc. However, in order to deal with the complex cardiac images, the prior knowledge about the shape and intensity of the LV is essential for these segmentation methods. Segmentation methods that take advantage of the prior knowledge, including the active shape model (ASM),10–12 active appearance model (AAM),15–17 sparse shape composition,13,14 atlas or registration-based techniques,18–21 active contour model,22–25 and optimal graph search,26,27 have demonstrated considerable competence in the LV segmentation. The ASM method first builds a statistical shape model by applying the principal component analysis on a set of aligned shapes, and then fits the model to a target image by adjusting the model parameters to match the image content. AAM originates from the extension of the ASM to include the gray level modeling and can capture the texture variations. AAM and ASM are competent in the applications of segmenting the cardiac images with noisy structures and weak edges, but they are also sensitive to the initializations. The sparse shape representation is recently proposed to approximate the shape of a target image using sparse linear combination of other shape instances, which have been previously annotated. Unlike ASM or AAM, this method employs the sparsity constraint to derive a smooth and realistic shape. For the atlas-based methods, the labeled atlas is registered to the target image and then the resultant transformation is applied to the Medical Physics, Vol. 42, No. 3, March 2015
1425
atlas, thereby obtaining the final segmentation. Optimal graph searching is developed to simultaneously segmenting multiple interacting surfaces in which the optimality is controlled by the cost functions designed for individual surfaces and by several geometric constraints defining the surface smoothness and interrelations. Though these methods with priors can achieve good performance in many cases, the robustness can be limited because of the dependence on the training data. More details of the segmentation methods based on priors for LV segmentation can be found in the recent papers of Refs. 3–5 To reduce the dependence on the choice of a training set, several approaches28–31 for LV segmentation have been proposed. Ayed et al.28 proposed to segment the end- and epicardial boundaries simultaneously based on the manually delineated first frame. They introduced a max-flow based segmentation method by recovering subject-specific model distributions learned from the first frame. Nambakhsh et al.29 investigated fast detection of the LV endo- and epicardium surface in cardiac MRI with the convex-relaxation optimization and distribution matching. This algorithm requires a manual segmentation of the first frame for training. Wu et al.30 presented a parametric active contour model to delineate the myocardial boundary using the anatomical shape of the LV, rather than using the constraints derived from a finite training set. In their work, a shape-similarity based energy function was proposed to prevent the snake contour from being strapped into faulty edges and to preserve weak boundaries. Eslami et al.31 proposed a new segmentation framework with prior knowledge. The prior knowledge is integrated into a new formulation of the random walks method. This method does not extract any principal model of the shape or intensity of the LV. In this work, we proposed a data-driven, hybrid segmentation method based on the dynamic programming (DP) framework. The DP method can incorporate the shape constraint when searching the optimal contour. However, the DP technique has not yet been fully investigated and only a small number of works22–35 have been reported for the LV segmentation. Geiger et al.32 used the DP technique to segment the LV from the cardiac MRI. The fuzzy sets were incorporated into the DP segmentation framework by Lalande et al.33 Fu et al.,34 applied a wavelet-based image technique to enhance the left ventricular endocardial and epicardial profiles, and then used a DP-based method to detect the borders. Recently, Kurkure et al.35 developed a hybrid segmentation method that the initial segmentation was obtained by a region growing segmentation method based on the fuzzy connectedness framework, and then the final contour was delineated by the DP technique. In this work, we improved the DP framework for the segmentation of the LV, particularly for the epicardium. First, the endocardium, i.e., the LV blood cavity, is segmented using a simple but robust algorithm for initialization. This reliable result provides a shape prior for the segmentation of the epicardium, which is more challenging. Second, a likelihood image of the myocardium and an adaptive search range (ASR) are constructed for the DP, based on the estimation of the myocardium thickness. Finally, a novel connectivity constrained DP (CCDP) approach is used to search the epicardial
1426
Qian et al.: Segmentation of myocardium from cardiac MR
1426
contour that minimizes a cost function, which incorporates both the probability information of the myocardium and the shape information of the endocardium. The main methodological novelties of the proposed DP method include (1) the likelihood map of the myocardium, (2) the adaptive search range, and (3) the connectivity constrained DP algorithm.
and adaptive search range (Sec. 2.B.2), and the connectivity constraint (Sec. 2.B.3). Figure 1 provides an outline of the proposed segmentation scheme.
1. The likelihood map of the myocardium is constructed to define the external cost in the DP. By this way, the prior probability estimation can be incorporated into the cost function. Besides, the likelihood image is able to suppress the gradients of nonmyocardium tissues in the gradient map, thereby restraining the negative effect when segmenting the epicardium from the nonmyocardium tissues. 2. The adaptive search range is introduced to determine the polar-transformed image, to exclude the irrelevant tissues. It imposes a soft shape constraint on the epicardium segmentation, while the DP can exert the hard shape constraint with the segmented endocardium to overcome the low contrast regions such as the boundaries between the myocardium and neighboring tissues. 3. The connectivity constrained DP algorithm is developed to overcome the Knapsack problem,36 i.e., to maintain the continuity between the beginning and the end of the polar coordinate image.
Since the entire segmentation process relies on the robust localization of the LV, our method starts with a user interaction to set the central location of LV on the middle slice of a selected phase, normally the first phase or the end diastolic (ED) phase. Then, this center point is propagated to the other slices and phases. This is based on the assumption that the central locations of the adjacent slices and phases are very close to each other. The ROI on a slice is defined as an 80 × 90 pixel window centered at the propagated LV central location. The size of ROI is chosen to cover the entire LV and myocardium, and a small part of the right ventricle (RV). Note that in order to minimize the user input, only one mouse click is required for the whole cine magnetic resonance (MR) sequence.
2. METHOD The proposed DP-based myocardium segmentation technique is composed of two phases, i.e., an initialization phase (Sec. 2.A) and a segmentation phase (Sec. 2.B). In the initialization phase, a series of preparations are performed for the subsequent epicardium segmentation, including region of interest (ROI) localization (Sec. 2.A.1), initial segmentation (Sec. 2.A.2), likelihood image estimation (Sec. 2.A.3), and polar transformation with adaptive search range (Sec. 2.A.4). In the segmentation phase, a modified DP framework is employed to segment the epicardium. This DP method adopts a new cost function which incorporates the likelihood image
2.A. Initialization phase 2.A.1. ROI localization
2.A.2. Initial segmentation
The ROI determined in Subsection 2.A.1 is used for the initial segmentation in this subsection. The main purpose of this process is to segment the LV and RV. Based on the septum between the LV and RV, we can estimate the lower bound of the myocardium thickness, with which we can further obtain the statistical property of the myocardium and convert the image into the likelihood image. By using the bias-corrected FCM method,37 the tissues included in the ROI can be initially classified into high (e.g., blood pool), medium (e.g., myocardium), and low (e.g., lung) intensity regions. The LV blood cavity is then identified as the connected region covering the LV center identified in the LV localization step, and the other largest connected region with similar intensity is extracted as the RV. Finally, the endocardium is delineated by a convex hull operation on the extracted blood cavity, followed by a Fourier transform-based contour
F. 1. Outline of the proposed segmentation scheme. Medical Physics, Vol. 42, No. 3, March 2015
1427
Qian et al.: Segmentation of myocardium from cardiac MR
1427
F. 2. The segmentation for the endocardium. (a) The ROI determined using the ROI localization (Sec. 2.A.1); (b) the segmentation result of the left ventricle; (c) the segmentation result of the right ventricle; (d) the final segmentation result of the endocardium.
smoothing.38 Figure 2 shows the segmentation process and results of the endocardium. 2.A.3. Likelihood image estimation
In order to estimate the thickness of the myocardium, we first discretize contours of the epicardium based on the segmented LV and RV. This estimation is illustrated in Fig. 3. Without loss of generality, for each image slice, we set the origin to the barycenter of LV such that the intensity image can be expressed as I(x i ,y i )(i = 1, ..., S), where S is the total pixel number of the selected slice, and (x i ,y i ) indicates the coordinate of the ith pixel. The discretized contour of LV, R1(θ k ), and RV, R2(θ k ), can be computed, respectively, as R1(θ k ) = max x 2i + y 2i | i ∈ ALV, y i = x i tan θ k (1) and R2(θ k ) = min
x 2j + y 2j | j ∈ ARV,y j = x j tan θ k ,
(2)
where we empirically set θ k = k∆θ for k = 1,...,180 and ∆θ = 2◦ in this work. ALV and ARV, respectively, indicate the pixel sets of the LV and RV. Here, we have empirically discretized the azimuthal angle from the LV center by 180 radial lines with increment ∆θ = 2◦. Then, the myocardium thickness at different angles can be computed as d k = R2 (θ k ) − R1(θ k ), and the thinnest myocardium thickness d m can be estimated by averaging the ten smallest myocardium thickness positions. With d m and R1(θ k ), a set of myocardial pixels can be located within d m/2 distance outward from the endocardium
boundary, or formally, M (θ k ) = {i| y i = x i tan θ k and 0 2 2 ≤ x i + y i − R1(θ k ) ≤ d m/2}. The located pixels are then used to estimate the statistical properties of the myocardium in terms of the mean value (µm ) and standard deviation (σ m ). Such that, the intensity image, I(x i ,y i ), can be converted into the likelihood image, I ′(x i ,y i ), with the following formula by assuming a statistical model for pixel intensities: 2 I ′(x i ,y i ) ∝ exp − I(x i ,y i ) − µm /(σ m )2 . (3) The converted myocardium likelihood images can minimize the influence from nonmyocardium tissue structures as shown in Fig. 4(a). 2.A.4. Polar transformation with adaptive search range
To facilitate the segmentation process with the DP method, we first defined the adaptive ring-shape search region of the myocardium to be starting from the endocardium and the width being RMAX = min{1.2 d m + 4.20}. Note that the maximal search range is limited to 20 mm, since the thickness of myocardium generally ranges from 6 to 16 mm. In Fig. 4(b), the adaptive search range is marked on an original MR image. Then, the ring-shape search region of the myocardium in Cartesian coordinates is transformed to an image in polar coordinates. We defined a (quasi-) polar transformation for the likelihood image by P(k, y) = I ′ ([ y + R1(θ k )]cos θ k ,[ y + R1(θ k )]sin θ k ),
(4)
where 0 ≤ y − R1 (θ k ) ≤ RMAX. Figures 4(c) and 4(d), respectively, depict the polar-transformed images from the original MR image and the likelihood image within the adaptive search range. 2.B. Segmentation phase with dynamic programming 2.B.1. Workflow of the dynamic programming based segmentation
F. 3. Illustration of the estimation of myocardium thickness. Medical Physics, Vol. 42, No. 3, March 2015
The DP scheme searches for a globally optimal path by minimizing the path cost consisting of both image information (i.e., external energy) and smoothness constraints (i.e., internal energy).
1428
Qian et al.: Segmentation of myocardium from cardiac MR
1428
F. 4. Likelihood image estimation and polar transformation: (a) the likelihood image of the myocardium; (b) ring-shape search region for epicardium; (c) polar transformation result of the image within adaptive search region: arrow A is the right ventricle, B is myocardium, and C is lung; (d) polar transformation result of the likelihood image.
For the polar-transformed image with N angular sample steps (N = 180 in this work), the minimal cost path out of all possible paths from step 1 to step N, ( y1, ..., y N ), can be found recursively in the DP framework. In particular, the minimal energy of all possible paths from step 1 to vertical coordinates yk at step k, Ekmin ( yk ), must satisfy the following recursive relation: min ( yk−1) + Dk ( yk ) Ekmin( yk ) = min Ek−1 yk −1
(5)
with E1min = D1( y1). Dk ( yk ) is the local cost introduced in the k-th step at the coordinates yk . The minimizer mk−1 ( yk ) radial min ( yk−1) + Dk ( yk ) records the optimal ver= argmin yk −1 Ek−1 tical location at (k − 1)-th step on the path leading to yk in k steps. At the N-th step, the optimal path coordinate is min given by yˆ N ≡ m N = argmin y N Ek−1 ( y N ). The optimal path, ( yˆ1, yˆ2,..., yˆ N ), is then reconstructed by backtracking minimizers at all steps starting from yˆ N , i.e., yˆ k−1 = mk−1 ( yˆ k ). The optimal path is transformed back onto the original image, inverse transformation of Eq. (5), to generate a contour line which gives the segmentation result of the epicardium.
− ∇ y P(k, y) Dkext( yk ) = 0
yk ≥ 0.25 d m , o.w.
(7)
where ∇ y represents the gradient operator in y-direction. The internal cost Dkext ( yk ) is introduced to prevent random jumps and fluctuations along the optimal contour path. In the DP, the path smoothness at yk can be enforced by assigning large costs for high deviation from the vertical position in the previous step, yk−1. Additionally, in this study, we proposed a second constraint from the optimal vertical location at (k −2)th step (i.e., second order smoothness constraint), mk−2 ( yk−1), such that the local internal cost is a weighted combination, int Dkint( yk ; yk−1) = D int k ( yk − yk−1) + wint D k ( yk − m k−2( yk−1)), (8)
= ( yk − yk−1)2 and D int = ( yk − mk−2)2. D int , called where D int k k k the penalty term of the interval points, restricts the contour from the big jumps between the interval vertical columns, thus it can smooth the transitions of the contour path. In fact, D int and D int represent the distance constraint between the k k epicardium and endocardium, since the two boundaries are nearly parallel and have a distance within a specific range.
2.B.2. New cost functions
We consider the local cost function in the k-th step at the radial location yk that can be decomposed into an external cost Dkext ( yk ) and an internal cost Dkint ( yk ), Dk ( yk ) = wt Dkext( yk ) + (1 − wt )Dkint( yk ),
(6)
where wt represents the weight between the internal and external costs. For the epicardium segmentation in this study, the external cost Dkext ( yk ) is defined as the negative magnitude of the local gradients in y-direction on the polar-transformed likelihood image. Therefore, the minimal path is attracted to the points with strong vertical gradients on the polar-transformed images. In order to minimize the influence due to intensity fluctuations near the endocardium (i.e., yk = 0, k = 1, ..., N), the local external costs are set to zero for yk < 0.25 d m , k = 1, ..., N. So, Medical Physics, Vol. 42, No. 3, March 2015
2.B.3. Connectivity constrained dynamic programming algorithm
The optimal path finding of the traditional DP, ( yˆ1, yˆ2, . . . , yˆ N ), does not guarantee a closed contour, since the path ends yˆ1 and yˆ N are not necessarily of the same values. This is called the Knapsack problem.36 Figure 5(a) demonstrates the minimal path on the polar image with unequal vertical locations at the path ends. This Knapsack problem results in an open contour for the epicardium, which is considered unrealistic [Fig. 5(c)]. Here, we propose an iterative correction method when contour discontinuities are detected. First, an unconnected contour is determined if | yˆ1 − yˆ N | > t, where t is the tolerance for distance discontinuity. The selection of t depends on the performance requirement, in this work, we set t = 3 pixels. Second, for the discontinuous contour, the first NS steps of the polar-transformed image are concatenated to the end of
1429
Qian et al.: Segmentation of myocardium from cardiac MR
1429
F. 5. Demonstration of the connectivity constrained dynamic programming: (a) the minimum path (red line) in the polar gradient map; (b) the minimum path in the newly formed image; (c) discontinuity of the outlined epicardium contour using the traditional method; and (d) the result using the proposed method.
the original polar-transformed image [Fig. 5(b)], such that the newly formed image has Ns + N steps and the external cost on the new image is given by for k ≤ N Dkext( yk ) D˜ kext( yk ) = D ext ( yk ) for k > N . k−N
(9)
Third, a new minimal path, y˜1, y˜2,..., y˜ N +N s is determined on the newly formed polar-transformed image using DP. Finally, the starting end point of the contour satisfying continuity constraint is located on the new minimal path as min{i| 2j=0 y˜ i+ j − y˜ i+N −1+ j < 3t,i = 1,...,Ns − 1}. The value of Ns can be increased iteratively until such point is found. In Fig. 5(d), the epicardium contour is closed after the discontinuity correction.
segmented and reference boundaries, defined by 1 RMSE = * N ,
(1) Dice metric: Let Vs represents the volume segmented by the algorithm and Vr represents the reference standard. The Dice is defined as Dice =
2Vs ∩Vr . Vs ∪Vr
(10)
The value of the Dice ranges from 0, no overlap between the two volumes, to 1, a perfect overlap. (2) Root mean squared error: The RMSE computes the distance between the corresponding points on the Medical Physics, Vol. 42, No. 3, March 2015
1/2
(x s,i − x r,i ) + ( ys,i − yr,i ) 2
2+
i=1
,
(11)
-
where (x s,i , ys,i ) is a point on the segmented boundary, and the (x r,i , yr,i ) is the closest point to (x s,i , ys,i ) on the reference boundary. The lower the RMSE is the better the accuracy of the method is in determining the boundaries. In order to compare with the other works in the literature, the average perpendicular distance (APD)39 was also used in the validation study. (3) Reliability: The reliability of segmentation method is assessed by the reliability function, defined as R (d) =
2.C. Evaluation metrics
Four measures were used to evaluate the performance of the proposed segmentation method: the Dice metric (Dice), the root mean squared error (RMSE), the reliability (R),28 and the correlation coefficient (R). Here are explanations of the four measures.
N
Number of volumes segmented with Dice >d , Total number of volumes (12)
where d ∈ [0,1]. R (d) indicates the reliability in yielding Dice d. The higher R, the better performance. (4) Correlation coefficient: The R between Vs and Vr is used to evaluate the quality of a least squares fitting, given by n R=
n
Vs,iVr,i −
n
i=1
i=1 1/2
Vs,i
n
Vr,i
i=1 1/2
2 2 n n n n *.n V 2 − * V + +/ *.n V 2 − * V + +/ s,i r,i s,i r,i ,i=1 - - , i=1 , i=1 - , i=1
.
(13)
The value of R is between 0 and 1. Higher R indicates better match between segmented volumes and the reference standard.
1430
Qian et al.: Segmentation of myocardium from cardiac MR
T I. Details of the two datasets.
Number of subjects Image resolution Temporal resolution Pixel spacing (mm) Slice thickness (mm)
Private dataset
MICCAI 2009
19 256 × 256 20 1.25∼1.52 13.8∼20
45 256 × 256 20 1.25∼1.56 8∼10
3. EXPERIMENTS AND RESULTS 3.A. Data sets
We tested the proposed method on two datasets: a private clinical dataset and a public dataset from the MICCAI 2009 LV Segmentation Grand Challenge. Table I gives the details of the two datasets. The private dataset has 2036 slices, acquired from 19 subjects. This data set was mainly used to qualitatively and visually evaluate the performance of our segmentation method, since the gold standard from our manual delineation can be biased. The quantitative assessment was mainly performed on the 45 subjects of the MICCAI challenge dataset, which consists of 15 online subjects, 15 training subjects, and 15 validation subjects.
1430
results were visually acceptable. Figure 6 displays the 20 phases of the middle slice from a randomly selected subject, which is a typical segmentation result. One can observe that the delineated endocardium (green curves) and epicardium (red curves) contours accord well with the truth boundaries. Our proposed method can minimize the influence of the trabeculae carneae and papillary muscles, and tackle the epicardium segmentation difficulties caused by the indistinct boundaries between the myocardium and its adjacent tissues. Figure 7 provides the visual results of three randomly selected slices, from the apical, middle-cavity, and basal slices, from the MICCAI challenge dataset. One can see that the delineated contours using the proposed method align well with the reference boundaries (yellow discontinuous curves) provided by the Challenge organizer. The Dice and RMSE corresponding to each column of Fig. 7 are listed in Table II. The relatively poor case is the one in the left column of Fig. 7(b), which has Dice of 0.82 and RMSE of 3.31 mm. The poor segmentation performance in this case was mainly due to the small size of the structure of the cavity at the apex and the poor image quality caused by motion artifacts. Nevertheless, the two examples demonstrate that our method can reliably segment the myocardium. 3.C. Quantitative evaluation results
3.B. Qualitative evaluation
3.C.1. The performance of the proposed method
The proposed method was applied to segment the endocardium and epicardium on the private dataset, and the
We quantitatively evaluated the segmentation results on the 45 cases from the MICCAI Grand Challenge dataset using
F. 6. Results of 20 phases for a slice from the private dataset. The inner and outer curves depict the endo- and epicardium boundaries, respectively. Medical Physics, Vol. 42, No. 3, March 2015
1431
Qian et al.: Segmentation of myocardium from cardiac MR
1431
F. 7. Results of three randomly selected cases from the MICCAI 2009: (a) online dataset, (b) training dataset, (c) validation dataset. The images from top to bottom are the apex slices, the middle slices, and the basal slices. For each case, the images in the left column belong to the ES phase, and the images in the right column belong to the ED phase. The yellow discontinuous curves depict the reference contours provided by the Challenge organizer.
Dice, RMSE, the Reliability (R), and correlation coefficient (R), as provided in Table III. The mean Dice for endocardium and epicardium is 0.892 and 0.927, respectively; the mean RMSE is 2.30 mm for endocardium and 2.39 mm for epicardium, respectively. The analysis results of the two metrics confirm the good performance of the proposed method, as the LV segmentation results reported in Ref. 28 show that the mean Dice higher than 0.80 indicates an excellent agreement with the manual segmentations, and the mean RMSE less than 3.0 mm is generally considered as good conformity of the results to the reference boundaries. As shown in Table III, the correlation coefficients are 0.98 for endocardium and 0.99 for epicardium. The linear regression plotted in Fig. 8 indicates strong correlation between the results of the proposed method and the reference standard. The proposed method yields a reliability of R (0.85) = 1 for epicardium detection, namely, all these cases have a good agreement (Dice > 0.8), and R (0.85) = 0.97 for endocardium detection. Figure 9(a) plots R as a function of d (d ≥ 0.7) for the endocardium and epicardium. It further shows the acceptable performance of the proposed method. 3.C.2. The effects of the three new techniques
The three new techniques, i.e., the ASR, the likelihood image of the myocardium (PMM), and the CCDP algorithm, are introduced into the DP framework to improve the performance T II. The quantitative performance (Dice and RMSE) for each column of Fig. 7. (a) SC-HF-I-10 ES
ED
Endo Endo
(b) SC-HF-I-40 ES
Epi
ED
Endo Endo
(c) SC-HF-I-5 ES
Epi
T III. Quantitative performance evaluations on 45 cases for proposed method. The last column: reliability of the Dice of 0.85 [R (0.85)].
ED
Endo Endo
Dice Epi
Dice 0.946 0.949 0.936 0.819 0.947 0.950 0.937 0.936 0.938 RMSE (mm) 1.85 1.84 2.80 3.31 1.25 1.51 1.53 1.93 2.22
Medical Physics, Vol. 42, No. 3, March 2015
of the DP-based segmentation scheme. Three experiments were performed on the 45 Grand Challenge cases, to explore the effects of these techniques in the proposed method. To examine the advantage of the likelihood image, we used the image intensity, instead of the likelihood image, for gradients computation in Eq. (6), and the other implementation of the segmentation framework remained the same. This is referred to as the method without PMM. The results in Table IV show that the mean Dice is 0.898, while the proposed method achieved 0.927; the mean RMSE is 3.37 mm, 0.91 mm worse than the proposed method; the R (0.85) = 0.98, while R (0.85) = 1 from the proposed method. These demonstrate the proposed method based on the likelihood image of the myocardium can achieve better performance than the method without PMM. In the second experiment, we performed the segmentation using a fixed search range, 20 mm, referred to as the method without ASR. The epicardium segmentation results in Table IV show that the mean Dice, mean RMSE, and R (0.85) are 0.880, 3.70, and 0.80 mm, respectively. It is apparent that the fixed search range technique with 20 mm is inferior to the proposed method using the adaptive search range. Finally, the method using the classic DP scheme without the connectivity constraint was examined. This is referred to as the method without CCDP. Its performance, shown in Table IV, is as follows: The mean Dice is 0.920, 0.007 less than that of the proposed method, and the mean RMSE is 2.63 mm, 0.23 mm worse than the proposed method.
RMSE (mm)
45 cases
Mean
Std
Mean
Std
R
R (0.85)
Endo Epi
0.892 0.927
0.043 0.017
2.30 2.39
0.67 0.66
0.98 0.99
0.97 1
1432
Qian et al.: Segmentation of myocardium from cardiac MR
F. 8. Segmentation volumes of our method versus manual volumes: (a) Endocardium detection: the obtained correlation coefficient is 0.98; (b) epicardium detection: the obtained correlation coefficient is 0.99.
In Fig. 9(b), one can find that the reliability curve of the method without CCDP is below that of the proposed method, which confirms the effectiveness of the connectivity constraint for the DP scheme. 3.C.3. The results from the literature
We compared the results of the proposed method with the results reported in the literature9,39–41 in terms of the Dice and APD for the endocardium detection and epicardium segmentation on the validation dataset from the MICCAI 2009, in Table V. The statistics of the Dice for the proposed method are 0.886 ± 0.039 for the endocardium and 0.928 ± 0.016 for the endocardium; and the statistics of the APD of our method are also similar to the other methods. This demonstrates that our method obtains the acceptable performance in terms of the endo- and epicardium segmentation. For reference, Table VI also provides the segmentation results from the literature12,17,28,29,42 in terms of the RMSE for the myocardium segmentation. It should be noted that an objective interwork comparison would be difficult as the test data and implementation of the different works could be very different. 4. DISCUSSION The focus of this work is to explore the segmentation ability of the new DP framework for the epicardium in Medical Physics, Vol. 42, No. 3, March 2015
1432
F. 9. The reliability of our method: (a) evaluated for the endocardium and epicardium segmentation, respectively; (b) the methods without the proposed new techniques evaluated for the epicardium segmentation.
cardiac cine MRI. We proposed three new techniques to advance the DP framework, i.e., the likelihood image of the myocardium, the adaptive search range, and the connectivity constrained DP algorithm. Here, we further discuss these techniques. The main purpose of constructing the myocardial likelihood image is to suppress the gradients of nonmyocardium tissues in the gradient map for the cost function. The gradients of nonmyocardium tissues may increase the external cost for the DP and deteriorate the optimal searching results of the epicardial contour. The average Dice of the segmentation method without the likelihood image is 0.3 less than that of our proposed method, meaning the likelihood image can improve the performance of the DP framework.
T IV. Comparisons of the epicardium segmentation to show the advantages of the likelihood image of the myocardium (PMM), ASR, and CCDP algorithm. RMSE (mm)
Dice 45 cases
Mean
Std
Mean
Std
R
R (0.85)
Proposed method Method without PMM Method without ASR Method without CCDP
0.927 0.898 0.880 0.920
0.017 0.021 0.031 0.019
2.39 3.37 3.70 2.62
0.66 0.76 0.94 0.68
0.99 0.97 0.97 0.98
1 0.98 0.80 1
1433
Qian et al.: Segmentation of myocardium from cardiac MR
T V. Comparisons of the segmentation accuracy of other approaches in terms of the Dice and APD on the validation dataset of the MICCAI 2009 Grand Challenge. Dice Endo
APD (mm) Epi
Validation references
Mean
Std
Mean
Schaerer et al. (Ref. 9) Jolly et al. (Ref. 37) Marak et al. (Ref. 38) Feng et al. (Ref. 36) Ours
0.869 0.880 0.860 0.890 0.886
0.043 0.030 0.040 0.040 0.039
0.921 0.930 0.930 0.940 0.928
Endo
Epi
Std Mean Std Mean Std 0.020 0.020 0.010 0.020 0.016
2.97 2.44 3.00 1.93 2.01
0.38 0.62 0.59 0.37 0.50
3.14 2.05 2.60 1.64 1.98
0.33 0.59 0.38 0.42 0.51
The adaptive search range is introduced to exclude the irrelevant tissues, thereby reducing the risk of undesirable gradients in the searching range. Whereas, the fixed search range, in order to ensure the complete coverage of the myocardium, has to be set a large size, such as 20 mm. Table IV shows the worse performance of the method using a fixed search range. This is because in the likelihood image, some tissues with similar intensity to the myocardium, such as the lung and liver, are rendered with high probability values, and their boundaries to the myocardium are also indistinct. Hence, these tissues have negative impacts on the optimal contour search of the DP when using the fixed search range. The connectivity constraint is incorporated into the DP algorithm to overcome a so-called Knapsack problem, i.e., the starting point does not connect with the ending point as shown in Fig. 5: Demonstration of the connectivity constrained dynamic programming: (a) the minimum path (red line) in the polar gradient map; (b) the minimum path in the newly formed image; (c) discontinuity of the outlined epicardium contour using the traditional method; and (d) the result using the proposed method (a). With the connectivity constrained DP algorithm, the epicardium contour is assured to be closed, as shown in Fig. 5. Demonstration of the connectivity
1433
constrained dynamic programming: (a) the minimum path (red line) in the polar gradient map; (b) the minimum path in the newly formed image; (c) discontinuity of the outlined epicardium contour using the traditional method; and (d) the result using the proposed method (d). Note that the performance of the proposed method in Table IV is just slightly better than that of the method without the connectivity constraint, because the Knapsack problem only happened in a small number of cases in the experiment. The DP scheme can be adopted and easily integrated into the myocardium segmentation framework, thanks to its two advantages. First, DP can flexibly impose the shape constraints of the endocardium to guide the epicardium segmentation; thereby we can overcome the difficulties caused by indistinguishable boundaries. Second, DP can make a local search for the desired boundary, such that we can reduce the adverse effects from the other tissues and obtain a desirable segmentation result. The main advantage of our segmentation method is that it only needs a minimized user-input (i.e., a seed point from a mouse click) and it is data-driven. Our method does not require any prior knowledge of the epicardium (such as shape information or local intensity statistics), instead this knowledge can be estimated from the MR images themselves for each phase. Therefore, the bias of prior knowledge derived from a limited training dataset can be avoided. A potential limitation of our segmentation scheme is that the epicardial segmentation is performed based on the endocardium segmentation, and an inaccurate endocardium segmentation result can affect the performance of the epicardium segmentation. Thus, we proposed an interaction for locating the center of the LV on the middle slice in the ED phase to pursue a reliable segmentation of the endocardium. As the result, our endocardium detection rate can reach 99.58%, based on which the detection rate of epicardium is 98.23%. Apparently, reliable endocardium segmentation can further improve the overall performance of our segmentation method. Finally, it is worth noting that the endocardium segmentation is not the emphasis of this work.
T VI. Comparisons with other methods in terms of the RMSE (mean ± std). References
Endo (mm)
Epi (mm)
Zhang et al. (Ref. 15) Zhu et al. (Ref. 12)
1.69 ± 0.38
1.89 ± 0.50
Active appearance model
50 subjects
0.69 ± 0.13
1.27 ± 0.18
22 subjects
Cordero et al. (Ref. 39) Cyrus et al. (Ref. 26)
1.37 ± 0.20
1.22 ± 0.17
Manual segmentation of ED phase and propagation using active shape model Deformable model
1.35 ± 0.13
2.05 ± 0.10
Ben Ayed et al. (Ref. 25)
1.60 (pixel)
1.99 (pixel)
Ours
2.30 ± 0.67
2.39 ± 0.66
Medical Physics, Vol. 42, No. 3, March 2015
Methodology
Manual segmentation first frame and propagation using convex relaxation Manual segmentation first frame and propagation using max-flow New dynamic programming
Data
43 subjects affected by acute myocardial infraction 20 subjects
20 subjects
45 subjects
1434
Qian et al.: Segmentation of myocardium from cardiac MR
The proposed method is implemented in using unoptimized code and takes about 25.9 s to segment one cardiac phase volume running on a personal computer equipped with a 2.4 GHz CPU.
5. CONCLUSION In this study, we presented a novel dynamic programming based method for myocardium segmentation from cardiac cine MRI. Our method only needs a minimized user-input (i.e., a seed point) and is data-driven. The knowledge of the epicardium (such as the shape information or local intensity statistics) is estimated for each phase individually; therefore no prior information from training data is required. The results demonstrated that the proposed method can obtain reliable and accurate myocardium segmentation, potentially providing a useful tool for academic research and developing new clinical applications of MRI for cardiology.
ACKNOWLEDGMENTS This work was supported in part by the National Natural Science Foundation of China (81301283), and in part by the SRF for ROCS, SEM of China (47th). a)Author
to whom correspondence should be addressed. Electronic mail: [email protected] 1D. L. Kasper, E. Braunwald, and A. Fauci, Harrison’s Principles of internal Medicine, 17th ed. (McGraw-Hill, New York, NY, 2008). 2V. Fuster, R. O’rourke, R. Walsh, and P. Pool-Wilson, Hurst’s the Heart, 12th ed. (McGraw-Hill, New York, NY, 2007). 3V. Tavakoli and A. Amini, “A survey of shaped-based registration and segmentation techniques for cardiac images,” Comput. Vis. Image Understanding 117, 966–989 (2013). 4X. Zhuang, “Challenges and methodologies of fully automatic whole heart segmentation: A review,” J. Healthcare Eng. 4, 371–408 (2013). 5C. Petitjean and J. Dacher, “A review of segmentation methods in short axis cardiac MR images,” Med. Image Anal. 15, 169–184 (2011). 6X. Lin, B. Cowan, and A. Young, “Automated detection of left ventricle in 4D MR images: Experience from a large study,” in Proceedings of the 9th International Conference on Medical Image Computing and ComputerAssisted Intervention(MICCAI), Nagoya, Japan (Springer, Berlin, Heidelberg, 2006), pp. 728–735. 7C. Cocosco, W. Niessen, T. Netsch, E. Vonken, G. Lund, A. Stork, and M. A. Viergever, “Automatic image-driven segmentation of the ventricles in cardiac cine MRI,” J. Magn. Reson. Imaging 28, 366–374 (2008). 8H. Lee, N. Codella, M. Cham, J. Weinsaft, and Y. Wang, “Automatic left ventricle segmentation using iterative thresholding and an active contour model with adaptation on short-axis cardiac MRI,” IEEE Trans. Biomed. Eng. 57, 905–913 (2010). 9J. Schaerer, C. Casta, J. Pousin, and P. Clarysse, “A dynamic elastic model for segmentation and tracking of the heart in MR image sequences,” Med. Image Anal. 14, 738–749 (2010). 10S. C. Mitchell, J. G. Bosch, B. P. F. Lelieveldt, R. J. van der Geest, J. H. Reiber, and M. Sonka, “3-D active appearance models: Segmentation of cardiac MR and ultrasound images,” IEEE Trans. Med. Imaging 21, 1167–1178 (2002). 11H. C. V. Assen, M. G. Danilouchkine, M. S. Dirksen, J. H. C. Reiber, and B. P. F. Lelieveldt, “A 3-D active shape model driven by fuzzy inference: Application to cardiac CT and MR,” IEEE Trans. Inf. Technol. Biomed. 12, 595–605 (2008). 12Y. Zhu, X. Papademetris, A. J. Sinusas, and J. S. Duncan, “Segmentation of the left ventricle from cardiac MR images using a subject-specific dynamical model,” IEEE Trans. Med. Imaging 29, 669–687 (2010). Medical Physics, Vol. 42, No. 3, March 2015
1434 13S.
Zhang, Y. Zhan, M. Dewan, J. Huang, D. N. Metaxas, and X. S. Zhou, “Towards robust and effective shape modeling: Sparse shape composition,” Med. Image Anal. 16, 265–277 (2012). 14S. Zhang, Y. Zhan, and D. N. Metaxas, “Deformable segmentation via sparse representation and dictionary learning,” Med. Image Anal. 16, 1385–1396 (2012). 15A. Andreopoulos and J. K. Tsotsos, “Efficient and generalizable statistical models of shape and appearance for analysis of cardiac MRI,” Med. Image Anal. 12, 335–357 (2008). 16S. Zambal, J. Hladuvka, and 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 ComputerAssisted Intervention(MICCAI), Copenhagen, Denmark (Springer, Berlin, Heidelberg, 2006), pp. 151–158. 17H. Zhang, A. Wahle, R. K. Johnson, T. D. Scholz, and M. Sonka, “4-D cardiac MR image analysis: Left and right ventricular morphology and function,” IEEE Trans. Med. Imaging 29, 350–364 (2010). 18M. Lorenzo-Valdés, G. I. Sanchez-Ortiz, A. G. Elkington, R. H. Mohiaddin, and D. Rueckert, “Segmentation of 4d cardiac MR images using a probabilistic atlas and the EM algorithm,” Med. Image Anal. 8, 255–265 (2004). 19X. Zhuang, K. S. Rhode, S. R. Arridge, R. Razavi, D. L. G. Hill, D. J. Hawkes, and S. Ourselin, “An atlas-based segmentation propagation framework using locally affine registration–application to automatic whole heart segmentation,” in Proceedings of the 11th International Conference on Medical Image Computing and Computer-Assisted Intervention(MICCAI), New York, NY (Springer, Berlin, Heidelberg, 2008), pp. 425–433. 20X. Zhuang, K. S. Rhode, R. S. Razavi, D. J. Hawkes, and S. Ourselin, “A registration-based propagation framework for automatic whole heart segmentation of cardiac MRI,” IEEE Trans. Med. Imaging 29, 1612–1625 (2010). 21W. Bai, W. Shi, D. P. O’regan, T. Tong, H. Wang, S. Jamil-Copley, and D. Ruechert, “A probabilistic patch-based label fusion model for multiatlas segmentation with registration refinement: Application to cardiac MR images,” IEEE Trans. Med. Imaging 32, 1302–1315 (2013). 22N. Paragios, “A variational approach for the segmentation of the left ventricle in cardiac image analysis,” Int. J. Comput. Vis. 50, 345–362 (2002). 23N. Paragios, “A level set approach for shape-driven segmentation and tracking of the left ventricle,” IEEE Trans. Med. Imaging 22, 773–776 (2003). 24J. Folkesson, E. Samset, R. Kwong, and C. Westin, “Unifying statistical classication and geodesic active regions for segmentation of cardiac MRI,” IEEE Trans. Inf. Technol. Biomed. 12, 328–334 (2008). 25I. B. Ayed, S. Li, and I. Ross, “Embedding overlap priors in variational left ventricle tracking,” IEEE Trans. Med. Imaging 28, 1902–1913 (2009). 26K. Li, X. Wu, D. Z. Chen, and M. Sonka, “Optimal surface segmentation in volumetric images-a graph-theoretic approach,” IEEE Trans. Pattern Anal. Mach. Intell. 28, 119–134 (2006). 27K. Li, S. Millington, X. Wu, D. Z. Chen, and M. Sonka, “Simultaneous segmentation of multiple closed surfaces using optimal graph searching,” in Proceedings of the 19th International Conference on Information Processing in Medical Imaging (IPMI), Glenwood Springs, CO (Springer, Berlin, Heidelberg, 2005), pp. 406–417. 28I. B. Ayed, H. M. Chen, K. Punithakumar, I. Ross, and S. Li, “Max-flow segmentation of the left ventricle by recovering subject-specific distributions via a bound of the Bhattacharyya measure,” Med. Image Anal. 16, 87–100 (2012). 29C. Nambakhsh, J. Yuan, K. Punithakumar, A. Goela, M. Rajchl, T. M. Peters, and I. B. Ayed, “Left ventricle segmentation in MRI via convex relaxed distribution matching,” Med. Image Anal. 17, 1010–1024 (2013). 30Y. Wu, Y. Wang, and Y. Jia, “Segmentation of the left ventricle in cardiac cine MRI using a shape-constrained snake model,” Comput. Vis. Image Understanding 117, 990–1003 (2013). 31A. Eslami, A. Karamalis, A. Katouzian, and N. Navab, “Segmentation by retrieval with guided random walks: Application to left ventricle segmentation in MRI,” Med. Image Anal. 17, 236–253 (2013). 32D. Geiger, A. Gupta, L. A. Costa, and J. Vlontzos, “Dynamic programming for detecting, tracking, and matching deformable contours,” IEEE Trans. Pattern Anal. Mach. Intell. 17, 294–302 (1995). 33A. Lalande, L. Legrand, P. M. Walker, F. Guy, Y. Cottin, S. Roy, and F. Brunotte, “Automatic detection of left ventricular contours from cardiac cine
1435
Qian et al.: Segmentation of myocardium from cardiac MR
magnetic resonance imaging using fuzzy logic,” Invest. Radiol. 34, 211–217 (1999). 34J. Fu, J. Chai, and S. Wong, “Wavelet-based enhancement for detection of left ventricular myocardial boundaries in magnetic resonance images,” Magn. Reson. Imaging 18, 1135–1141 (2000). 35U. Kurkure, A. Pednekar, R. Muthupillai, S. D. Flamm, and I. A. Kakadiaris, “Localization and segmentation of left ventricle in cardiac cine-MR images,” IEEE Trans. Biomed. Eng. 56, 1360–1370 (2009). 36T. H. Cormen, C. E. Leiseron, and R. L. Rivest, Introduction to Algorithms, 23rd ed. (McGraw-Hill, New York, NY, 1999). 37M. N. Ahmed, S. M. Yamany, N. Mohamed, A. A. Farag, and T. Moriarty, “A modified fuzzy c-means algorithm for bias field estimation and segmentation of MRI data,” IEEE Trans. Med. Imaging 21, 193–199 (2002). 38W. Shi, X. Zhuang, H. Wang, S. Duckett, D. Oregan, P. Edwards, and D. Rueckert, “Automatic segmentation of different pathologies from cardiac cine MRI using registration and multiple component EM estimation,” in Proceedings of Functional Imaging and Modeling of the Heart (Springer, Berlin, Heidelberg, 2011), pp. 163–170.
Medical Physics, Vol. 42, No. 3, March 2015
1435 39C.
Feng, C. Li, D. Zhao, C. Davatzikos, and 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(MICCAI) Nagoya, Japan (Springer, Berlin, Heidelberg, 2013), pp. 477–484. 40M. P. Jolly, “Fully automatic left ventricle segmentation incardiac cine MR images using registration and minimum surfaces,” in MICCAI 2009 Workshop on Cardiac MR Left Ventricle Segmentation Challenge, MIDAS J. (2009). 41L. Marak, J. Cousty, L. Najman, and H. Talbot, “4D morphological segmentation and that MICCAI LV-Segmentation grand challenge,” in MICCAI 2009 Workshop on Cardiac MR Left Ventricle Segmentation Challenge, MIDAS J. (2009). 42L. Cordero-Grande, G. Vegas-Sánchez-Ferrero, P. Casaseca-de-la-Higuera, J. Alberto San-Román-Calvar, A. Revilla-Orode, M. Martín-Fernández, and C. Alberola-López, “Unsupervised 4D myocardium segmentation with a markov random field based deformable model,” Med. Image Anal. 15, 283–301 (2011).
Segmentation of myocardium from cardiac MR images using a novel dynamic programming based segmentation method Xiaohua Qian SJTUCU International Cooperative Research Center, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China and Shanghai Advanced Research Institute, Chinese Academy of Sciences, Shanghai 201210, China
Yuan Lin Division of Research and Innovations, Carestream Health, Inc., Rochester, New York 14615
Yue Zhao College of Electronic Science and Engineering, Jilin University, 2699 Qianjing Street, Changchun, Jilin 130012, China
Jing Wang Suzhou Institute of Biomedical Engineering and Technology Chinese Academy of Sciences, Suzhou, Jiangsu 215163, China
Jing Liu Shanghai Advanced Research Institute, Chinese Academy of Sciences, Shanghai 201210, China
Xiahai Zhuanga) SJTUCU International Cooperative Research Center, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
(Received 2 September 2014; revised 2 December 2014; accepted for publication 17 January 2015; published 27 February 2015) Purpose: Myocardium segmentation in cardiac magnetic resonance (MR) images plays a vital role in clinical diagnosis of the cardiovascular diseases. Because of the low contrast and large variation in intensity and shapes, myocardium segmentation has been a challenging task. A dynamic programming (DP) based segmentation method, incorporating the likelihood and shape information of the myocardium, is developed for segmenting myocardium in cardiac MR images. Methods: The endocardium, i.e., the left ventricle blood cavity, is segmented for initialization, and then the optimal epicardium contour is determined using the polar-transformed image and DP scheme. In the DP segmentation scheme, three techniques are proposed to improve the segmentation performance: (1) the likelihood image of the myocardium is constructed to define the external cost in the DP, thus the cost function incorporates prior probability estimation; (2) the adaptive search range is introduced to determine the polar-transformed image, thereby excluding irrelevant tissues; (3) the connectivity constrained DP algorithm is developed to obtain an optimal closed contour. Four metrics, including the Dice metric (Dice), root mean squared error (RMSE), reliability, and correlation coefficient, are used to assess the segmentation accuracy. The authors evaluated the performance of the proposed method on a private dataset and the MICCAI 2009 challenge dataset. The authors also explored the effects of the three new techniques of the DP scheme in the proposed method. Results: For the qualitative evaluation, the segmentation results of the proposed method were clinically acceptable. For the quantitative evaluation, the mean (Dice) for the endocardium and epicardium was 0.892 and 0.927, respectively; the mean RMSE was 2.30 mm for the endocardium and 2.39 mm for the epicardium. In addition, the three new techniques in the proposed DP scheme, i.e., the likelihood image of the myocardium, the adaptive search range, and the connectivity constrained DP algorithm, improved the segmentation performance for the epicardium with 0.029, 0.047, and 0.007 in terms of the Dice and 0.98, 1.31, and 0.21 mm in terms of the RMSE, respectively. Conclusions: The three techniques (the likelihood image of the myocardium, the adaptive search range, and the connectivity constrained DP algorithm) can improve the segmentation ability of the DP method, and the proposed method with these techniques has the ability to achieve the acceptable segmentation result of the myocardium in cardiac MR images. Therefore, the proposed method would be useful in clinical diagnosis of the cardiovascular diseases. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4907993] Key words: cardiac MRI, myocardium segmentation, likelihood image, adaptive search range, dynamic programming
1424
Med. Phys. 42 (3), March 2015
0094-2405/2015/42(3)/1424/12/$30.00
© 2015 Am. Assoc. Phys. Med.
1424
1425
Qian et al.: Segmentation of myocardium from cardiac MR
1. INTRODUCTION Cardiovascular diseases are the leading causes of deaths in the modern countries.1,2 Medical imaging technologies, such as the magnetic resonance imaging (MRI), computed tomography, echography, and nuclear medicine, are effective in reducing the mortality and morbidity of cardiovascular diseases. MRI provides information of the morphology and pathology of the heart and has emerged as a widely used noninvasive modality in clinical applications of cardiology.3 Extracting this information, for computing functional indices such as ejection fraction and wall thickness, relies on accurate delineation of the endocardial and epicardial boundaries of the left ventricle (LV). Manual delineation is labor-intensive and suffers from both inter- and intraexpert variability, therefore, automatic or semiautomatic algorithms become increasingly desirable. However, there are a number of challenges inherently from the cardiac MRI.4,5 First, due to the blood flow and partial volume effects, intensity patterns of the cardiac MRI vary across different scans. Second, structural shapes of the heart alter because of the different patients and different conditions. Moreover, the muscle tissues (with similar intensity) can be classified into anatomically different substructures, such as the myocardium, trabeculae carneae, and papillary muscles. Additionally, the myocardium and its neighboring tissues such as the liver and lung generally have the comparable intensity, due to the similar T2/T1 ratios. In recent years, a large number of segmentation approaches have been developed for the delineation of the LV endoand epicardium contours in cardiac MRI,5 including the approaches without prior6–9 and with prior.10–24 Methods without prior includes the threshold,6 clustering,7 snake/active contour model,8,9 etc. However, in order to deal with the complex cardiac images, the prior knowledge about the shape and intensity of the LV is essential for these segmentation methods. Segmentation methods that take advantage of the prior knowledge, including the active shape model (ASM),10–12 active appearance model (AAM),15–17 sparse shape composition,13,14 atlas or registration-based techniques,18–21 active contour model,22–25 and optimal graph search,26,27 have demonstrated considerable competence in the LV segmentation. The ASM method first builds a statistical shape model by applying the principal component analysis on a set of aligned shapes, and then fits the model to a target image by adjusting the model parameters to match the image content. AAM originates from the extension of the ASM to include the gray level modeling and can capture the texture variations. AAM and ASM are competent in the applications of segmenting the cardiac images with noisy structures and weak edges, but they are also sensitive to the initializations. The sparse shape representation is recently proposed to approximate the shape of a target image using sparse linear combination of other shape instances, which have been previously annotated. Unlike ASM or AAM, this method employs the sparsity constraint to derive a smooth and realistic shape. For the atlas-based methods, the labeled atlas is registered to the target image and then the resultant transformation is applied to the Medical Physics, Vol. 42, No. 3, March 2015
1425
atlas, thereby obtaining the final segmentation. Optimal graph searching is developed to simultaneously segmenting multiple interacting surfaces in which the optimality is controlled by the cost functions designed for individual surfaces and by several geometric constraints defining the surface smoothness and interrelations. Though these methods with priors can achieve good performance in many cases, the robustness can be limited because of the dependence on the training data. More details of the segmentation methods based on priors for LV segmentation can be found in the recent papers of Refs. 3–5 To reduce the dependence on the choice of a training set, several approaches28–31 for LV segmentation have been proposed. Ayed et al.28 proposed to segment the end- and epicardial boundaries simultaneously based on the manually delineated first frame. They introduced a max-flow based segmentation method by recovering subject-specific model distributions learned from the first frame. Nambakhsh et al.29 investigated fast detection of the LV endo- and epicardium surface in cardiac MRI with the convex-relaxation optimization and distribution matching. This algorithm requires a manual segmentation of the first frame for training. Wu et al.30 presented a parametric active contour model to delineate the myocardial boundary using the anatomical shape of the LV, rather than using the constraints derived from a finite training set. In their work, a shape-similarity based energy function was proposed to prevent the snake contour from being strapped into faulty edges and to preserve weak boundaries. Eslami et al.31 proposed a new segmentation framework with prior knowledge. The prior knowledge is integrated into a new formulation of the random walks method. This method does not extract any principal model of the shape or intensity of the LV. In this work, we proposed a data-driven, hybrid segmentation method based on the dynamic programming (DP) framework. The DP method can incorporate the shape constraint when searching the optimal contour. However, the DP technique has not yet been fully investigated and only a small number of works22–35 have been reported for the LV segmentation. Geiger et al.32 used the DP technique to segment the LV from the cardiac MRI. The fuzzy sets were incorporated into the DP segmentation framework by Lalande et al.33 Fu et al.,34 applied a wavelet-based image technique to enhance the left ventricular endocardial and epicardial profiles, and then used a DP-based method to detect the borders. Recently, Kurkure et al.35 developed a hybrid segmentation method that the initial segmentation was obtained by a region growing segmentation method based on the fuzzy connectedness framework, and then the final contour was delineated by the DP technique. In this work, we improved the DP framework for the segmentation of the LV, particularly for the epicardium. First, the endocardium, i.e., the LV blood cavity, is segmented using a simple but robust algorithm for initialization. This reliable result provides a shape prior for the segmentation of the epicardium, which is more challenging. Second, a likelihood image of the myocardium and an adaptive search range (ASR) are constructed for the DP, based on the estimation of the myocardium thickness. Finally, a novel connectivity constrained DP (CCDP) approach is used to search the epicardial
1426
Qian et al.: Segmentation of myocardium from cardiac MR
1426
contour that minimizes a cost function, which incorporates both the probability information of the myocardium and the shape information of the endocardium. The main methodological novelties of the proposed DP method include (1) the likelihood map of the myocardium, (2) the adaptive search range, and (3) the connectivity constrained DP algorithm.
and adaptive search range (Sec. 2.B.2), and the connectivity constraint (Sec. 2.B.3). Figure 1 provides an outline of the proposed segmentation scheme.
1. The likelihood map of the myocardium is constructed to define the external cost in the DP. By this way, the prior probability estimation can be incorporated into the cost function. Besides, the likelihood image is able to suppress the gradients of nonmyocardium tissues in the gradient map, thereby restraining the negative effect when segmenting the epicardium from the nonmyocardium tissues. 2. The adaptive search range is introduced to determine the polar-transformed image, to exclude the irrelevant tissues. It imposes a soft shape constraint on the epicardium segmentation, while the DP can exert the hard shape constraint with the segmented endocardium to overcome the low contrast regions such as the boundaries between the myocardium and neighboring tissues. 3. The connectivity constrained DP algorithm is developed to overcome the Knapsack problem,36 i.e., to maintain the continuity between the beginning and the end of the polar coordinate image.
Since the entire segmentation process relies on the robust localization of the LV, our method starts with a user interaction to set the central location of LV on the middle slice of a selected phase, normally the first phase or the end diastolic (ED) phase. Then, this center point is propagated to the other slices and phases. This is based on the assumption that the central locations of the adjacent slices and phases are very close to each other. The ROI on a slice is defined as an 80 × 90 pixel window centered at the propagated LV central location. The size of ROI is chosen to cover the entire LV and myocardium, and a small part of the right ventricle (RV). Note that in order to minimize the user input, only one mouse click is required for the whole cine magnetic resonance (MR) sequence.
2. METHOD The proposed DP-based myocardium segmentation technique is composed of two phases, i.e., an initialization phase (Sec. 2.A) and a segmentation phase (Sec. 2.B). In the initialization phase, a series of preparations are performed for the subsequent epicardium segmentation, including region of interest (ROI) localization (Sec. 2.A.1), initial segmentation (Sec. 2.A.2), likelihood image estimation (Sec. 2.A.3), and polar transformation with adaptive search range (Sec. 2.A.4). In the segmentation phase, a modified DP framework is employed to segment the epicardium. This DP method adopts a new cost function which incorporates the likelihood image
2.A. Initialization phase 2.A.1. ROI localization
2.A.2. Initial segmentation
The ROI determined in Subsection 2.A.1 is used for the initial segmentation in this subsection. The main purpose of this process is to segment the LV and RV. Based on the septum between the LV and RV, we can estimate the lower bound of the myocardium thickness, with which we can further obtain the statistical property of the myocardium and convert the image into the likelihood image. By using the bias-corrected FCM method,37 the tissues included in the ROI can be initially classified into high (e.g., blood pool), medium (e.g., myocardium), and low (e.g., lung) intensity regions. The LV blood cavity is then identified as the connected region covering the LV center identified in the LV localization step, and the other largest connected region with similar intensity is extracted as the RV. Finally, the endocardium is delineated by a convex hull operation on the extracted blood cavity, followed by a Fourier transform-based contour
F. 1. Outline of the proposed segmentation scheme. Medical Physics, Vol. 42, No. 3, March 2015
1427
Qian et al.: Segmentation of myocardium from cardiac MR
1427
F. 2. The segmentation for the endocardium. (a) The ROI determined using the ROI localization (Sec. 2.A.1); (b) the segmentation result of the left ventricle; (c) the segmentation result of the right ventricle; (d) the final segmentation result of the endocardium.
smoothing.38 Figure 2 shows the segmentation process and results of the endocardium. 2.A.3. Likelihood image estimation
In order to estimate the thickness of the myocardium, we first discretize contours of the epicardium based on the segmented LV and RV. This estimation is illustrated in Fig. 3. Without loss of generality, for each image slice, we set the origin to the barycenter of LV such that the intensity image can be expressed as I(x i ,y i )(i = 1, ..., S), where S is the total pixel number of the selected slice, and (x i ,y i ) indicates the coordinate of the ith pixel. The discretized contour of LV, R1(θ k ), and RV, R2(θ k ), can be computed, respectively, as R1(θ k ) = max x 2i + y 2i | i ∈ ALV, y i = x i tan θ k (1) and R2(θ k ) = min
x 2j + y 2j | j ∈ ARV,y j = x j tan θ k ,
(2)
where we empirically set θ k = k∆θ for k = 1,...,180 and ∆θ = 2◦ in this work. ALV and ARV, respectively, indicate the pixel sets of the LV and RV. Here, we have empirically discretized the azimuthal angle from the LV center by 180 radial lines with increment ∆θ = 2◦. Then, the myocardium thickness at different angles can be computed as d k = R2 (θ k ) − R1(θ k ), and the thinnest myocardium thickness d m can be estimated by averaging the ten smallest myocardium thickness positions. With d m and R1(θ k ), a set of myocardial pixels can be located within d m/2 distance outward from the endocardium
boundary, or formally, M (θ k ) = {i| y i = x i tan θ k and 0 2 2 ≤ x i + y i − R1(θ k ) ≤ d m/2}. The located pixels are then used to estimate the statistical properties of the myocardium in terms of the mean value (µm ) and standard deviation (σ m ). Such that, the intensity image, I(x i ,y i ), can be converted into the likelihood image, I ′(x i ,y i ), with the following formula by assuming a statistical model for pixel intensities: 2 I ′(x i ,y i ) ∝ exp − I(x i ,y i ) − µm /(σ m )2 . (3) The converted myocardium likelihood images can minimize the influence from nonmyocardium tissue structures as shown in Fig. 4(a). 2.A.4. Polar transformation with adaptive search range
To facilitate the segmentation process with the DP method, we first defined the adaptive ring-shape search region of the myocardium to be starting from the endocardium and the width being RMAX = min{1.2 d m + 4.20}. Note that the maximal search range is limited to 20 mm, since the thickness of myocardium generally ranges from 6 to 16 mm. In Fig. 4(b), the adaptive search range is marked on an original MR image. Then, the ring-shape search region of the myocardium in Cartesian coordinates is transformed to an image in polar coordinates. We defined a (quasi-) polar transformation for the likelihood image by P(k, y) = I ′ ([ y + R1(θ k )]cos θ k ,[ y + R1(θ k )]sin θ k ),
(4)
where 0 ≤ y − R1 (θ k ) ≤ RMAX. Figures 4(c) and 4(d), respectively, depict the polar-transformed images from the original MR image and the likelihood image within the adaptive search range. 2.B. Segmentation phase with dynamic programming 2.B.1. Workflow of the dynamic programming based segmentation
F. 3. Illustration of the estimation of myocardium thickness. Medical Physics, Vol. 42, No. 3, March 2015
The DP scheme searches for a globally optimal path by minimizing the path cost consisting of both image information (i.e., external energy) and smoothness constraints (i.e., internal energy).
1428
Qian et al.: Segmentation of myocardium from cardiac MR
1428
F. 4. Likelihood image estimation and polar transformation: (a) the likelihood image of the myocardium; (b) ring-shape search region for epicardium; (c) polar transformation result of the image within adaptive search region: arrow A is the right ventricle, B is myocardium, and C is lung; (d) polar transformation result of the likelihood image.
For the polar-transformed image with N angular sample steps (N = 180 in this work), the minimal cost path out of all possible paths from step 1 to step N, ( y1, ..., y N ), can be found recursively in the DP framework. In particular, the minimal energy of all possible paths from step 1 to vertical coordinates yk at step k, Ekmin ( yk ), must satisfy the following recursive relation: min ( yk−1) + Dk ( yk ) Ekmin( yk ) = min Ek−1 yk −1
(5)
with E1min = D1( y1). Dk ( yk ) is the local cost introduced in the k-th step at the coordinates yk . The minimizer mk−1 ( yk ) radial min ( yk−1) + Dk ( yk ) records the optimal ver= argmin yk −1 Ek−1 tical location at (k − 1)-th step on the path leading to yk in k steps. At the N-th step, the optimal path coordinate is min given by yˆ N ≡ m N = argmin y N Ek−1 ( y N ). The optimal path, ( yˆ1, yˆ2,..., yˆ N ), is then reconstructed by backtracking minimizers at all steps starting from yˆ N , i.e., yˆ k−1 = mk−1 ( yˆ k ). The optimal path is transformed back onto the original image, inverse transformation of Eq. (5), to generate a contour line which gives the segmentation result of the epicardium.
− ∇ y P(k, y) Dkext( yk ) = 0
yk ≥ 0.25 d m , o.w.
(7)
where ∇ y represents the gradient operator in y-direction. The internal cost Dkext ( yk ) is introduced to prevent random jumps and fluctuations along the optimal contour path. In the DP, the path smoothness at yk can be enforced by assigning large costs for high deviation from the vertical position in the previous step, yk−1. Additionally, in this study, we proposed a second constraint from the optimal vertical location at (k −2)th step (i.e., second order smoothness constraint), mk−2 ( yk−1), such that the local internal cost is a weighted combination, int Dkint( yk ; yk−1) = D int k ( yk − yk−1) + wint D k ( yk − m k−2( yk−1)), (8)
= ( yk − yk−1)2 and D int = ( yk − mk−2)2. D int , called where D int k k k the penalty term of the interval points, restricts the contour from the big jumps between the interval vertical columns, thus it can smooth the transitions of the contour path. In fact, D int and D int represent the distance constraint between the k k epicardium and endocardium, since the two boundaries are nearly parallel and have a distance within a specific range.
2.B.2. New cost functions
We consider the local cost function in the k-th step at the radial location yk that can be decomposed into an external cost Dkext ( yk ) and an internal cost Dkint ( yk ), Dk ( yk ) = wt Dkext( yk ) + (1 − wt )Dkint( yk ),
(6)
where wt represents the weight between the internal and external costs. For the epicardium segmentation in this study, the external cost Dkext ( yk ) is defined as the negative magnitude of the local gradients in y-direction on the polar-transformed likelihood image. Therefore, the minimal path is attracted to the points with strong vertical gradients on the polar-transformed images. In order to minimize the influence due to intensity fluctuations near the endocardium (i.e., yk = 0, k = 1, ..., N), the local external costs are set to zero for yk < 0.25 d m , k = 1, ..., N. So, Medical Physics, Vol. 42, No. 3, March 2015
2.B.3. Connectivity constrained dynamic programming algorithm
The optimal path finding of the traditional DP, ( yˆ1, yˆ2, . . . , yˆ N ), does not guarantee a closed contour, since the path ends yˆ1 and yˆ N are not necessarily of the same values. This is called the Knapsack problem.36 Figure 5(a) demonstrates the minimal path on the polar image with unequal vertical locations at the path ends. This Knapsack problem results in an open contour for the epicardium, which is considered unrealistic [Fig. 5(c)]. Here, we propose an iterative correction method when contour discontinuities are detected. First, an unconnected contour is determined if | yˆ1 − yˆ N | > t, where t is the tolerance for distance discontinuity. The selection of t depends on the performance requirement, in this work, we set t = 3 pixels. Second, for the discontinuous contour, the first NS steps of the polar-transformed image are concatenated to the end of
1429
Qian et al.: Segmentation of myocardium from cardiac MR
1429
F. 5. Demonstration of the connectivity constrained dynamic programming: (a) the minimum path (red line) in the polar gradient map; (b) the minimum path in the newly formed image; (c) discontinuity of the outlined epicardium contour using the traditional method; and (d) the result using the proposed method.
the original polar-transformed image [Fig. 5(b)], such that the newly formed image has Ns + N steps and the external cost on the new image is given by for k ≤ N Dkext( yk ) D˜ kext( yk ) = D ext ( yk ) for k > N . k−N
(9)
Third, a new minimal path, y˜1, y˜2,..., y˜ N +N s is determined on the newly formed polar-transformed image using DP. Finally, the starting end point of the contour satisfying continuity constraint is located on the new minimal path as min{i| 2j=0 y˜ i+ j − y˜ i+N −1+ j < 3t,i = 1,...,Ns − 1}. The value of Ns can be increased iteratively until such point is found. In Fig. 5(d), the epicardium contour is closed after the discontinuity correction.
segmented and reference boundaries, defined by 1 RMSE = * N ,
(1) Dice metric: Let Vs represents the volume segmented by the algorithm and Vr represents the reference standard. The Dice is defined as Dice =
2Vs ∩Vr . Vs ∪Vr
(10)
The value of the Dice ranges from 0, no overlap between the two volumes, to 1, a perfect overlap. (2) Root mean squared error: The RMSE computes the distance between the corresponding points on the Medical Physics, Vol. 42, No. 3, March 2015
1/2
(x s,i − x r,i ) + ( ys,i − yr,i ) 2
2+
i=1
,
(11)
-
where (x s,i , ys,i ) is a point on the segmented boundary, and the (x r,i , yr,i ) is the closest point to (x s,i , ys,i ) on the reference boundary. The lower the RMSE is the better the accuracy of the method is in determining the boundaries. In order to compare with the other works in the literature, the average perpendicular distance (APD)39 was also used in the validation study. (3) Reliability: The reliability of segmentation method is assessed by the reliability function, defined as R (d) =
2.C. Evaluation metrics
Four measures were used to evaluate the performance of the proposed segmentation method: the Dice metric (Dice), the root mean squared error (RMSE), the reliability (R),28 and the correlation coefficient (R). Here are explanations of the four measures.
N
Number of volumes segmented with Dice >d , Total number of volumes (12)
where d ∈ [0,1]. R (d) indicates the reliability in yielding Dice d. The higher R, the better performance. (4) Correlation coefficient: The R between Vs and Vr is used to evaluate the quality of a least squares fitting, given by n R=
n
Vs,iVr,i −
n
i=1
i=1 1/2
Vs,i
n
Vr,i
i=1 1/2
2 2 n n n n *.n V 2 − * V + +/ *.n V 2 − * V + +/ s,i r,i s,i r,i ,i=1 - - , i=1 , i=1 - , i=1
.
(13)
The value of R is between 0 and 1. Higher R indicates better match between segmented volumes and the reference standard.
1430
Qian et al.: Segmentation of myocardium from cardiac MR
T I. Details of the two datasets.
Number of subjects Image resolution Temporal resolution Pixel spacing (mm) Slice thickness (mm)
Private dataset
MICCAI 2009
19 256 × 256 20 1.25∼1.52 13.8∼20
45 256 × 256 20 1.25∼1.56 8∼10
3. EXPERIMENTS AND RESULTS 3.A. Data sets
We tested the proposed method on two datasets: a private clinical dataset and a public dataset from the MICCAI 2009 LV Segmentation Grand Challenge. Table I gives the details of the two datasets. The private dataset has 2036 slices, acquired from 19 subjects. This data set was mainly used to qualitatively and visually evaluate the performance of our segmentation method, since the gold standard from our manual delineation can be biased. The quantitative assessment was mainly performed on the 45 subjects of the MICCAI challenge dataset, which consists of 15 online subjects, 15 training subjects, and 15 validation subjects.
1430
results were visually acceptable. Figure 6 displays the 20 phases of the middle slice from a randomly selected subject, which is a typical segmentation result. One can observe that the delineated endocardium (green curves) and epicardium (red curves) contours accord well with the truth boundaries. Our proposed method can minimize the influence of the trabeculae carneae and papillary muscles, and tackle the epicardium segmentation difficulties caused by the indistinct boundaries between the myocardium and its adjacent tissues. Figure 7 provides the visual results of three randomly selected slices, from the apical, middle-cavity, and basal slices, from the MICCAI challenge dataset. One can see that the delineated contours using the proposed method align well with the reference boundaries (yellow discontinuous curves) provided by the Challenge organizer. The Dice and RMSE corresponding to each column of Fig. 7 are listed in Table II. The relatively poor case is the one in the left column of Fig. 7(b), which has Dice of 0.82 and RMSE of 3.31 mm. The poor segmentation performance in this case was mainly due to the small size of the structure of the cavity at the apex and the poor image quality caused by motion artifacts. Nevertheless, the two examples demonstrate that our method can reliably segment the myocardium. 3.C. Quantitative evaluation results
3.B. Qualitative evaluation
3.C.1. The performance of the proposed method
The proposed method was applied to segment the endocardium and epicardium on the private dataset, and the
We quantitatively evaluated the segmentation results on the 45 cases from the MICCAI Grand Challenge dataset using
F. 6. Results of 20 phases for a slice from the private dataset. The inner and outer curves depict the endo- and epicardium boundaries, respectively. Medical Physics, Vol. 42, No. 3, March 2015
1431
Qian et al.: Segmentation of myocardium from cardiac MR
1431
F. 7. Results of three randomly selected cases from the MICCAI 2009: (a) online dataset, (b) training dataset, (c) validation dataset. The images from top to bottom are the apex slices, the middle slices, and the basal slices. For each case, the images in the left column belong to the ES phase, and the images in the right column belong to the ED phase. The yellow discontinuous curves depict the reference contours provided by the Challenge organizer.
Dice, RMSE, the Reliability (R), and correlation coefficient (R), as provided in Table III. The mean Dice for endocardium and epicardium is 0.892 and 0.927, respectively; the mean RMSE is 2.30 mm for endocardium and 2.39 mm for epicardium, respectively. The analysis results of the two metrics confirm the good performance of the proposed method, as the LV segmentation results reported in Ref. 28 show that the mean Dice higher than 0.80 indicates an excellent agreement with the manual segmentations, and the mean RMSE less than 3.0 mm is generally considered as good conformity of the results to the reference boundaries. As shown in Table III, the correlation coefficients are 0.98 for endocardium and 0.99 for epicardium. The linear regression plotted in Fig. 8 indicates strong correlation between the results of the proposed method and the reference standard. The proposed method yields a reliability of R (0.85) = 1 for epicardium detection, namely, all these cases have a good agreement (Dice > 0.8), and R (0.85) = 0.97 for endocardium detection. Figure 9(a) plots R as a function of d (d ≥ 0.7) for the endocardium and epicardium. It further shows the acceptable performance of the proposed method. 3.C.2. The effects of the three new techniques
The three new techniques, i.e., the ASR, the likelihood image of the myocardium (PMM), and the CCDP algorithm, are introduced into the DP framework to improve the performance T II. The quantitative performance (Dice and RMSE) for each column of Fig. 7. (a) SC-HF-I-10 ES
ED
Endo Endo
(b) SC-HF-I-40 ES
Epi
ED
Endo Endo
(c) SC-HF-I-5 ES
Epi
T III. Quantitative performance evaluations on 45 cases for proposed method. The last column: reliability of the Dice of 0.85 [R (0.85)].
ED
Endo Endo
Dice Epi
Dice 0.946 0.949 0.936 0.819 0.947 0.950 0.937 0.936 0.938 RMSE (mm) 1.85 1.84 2.80 3.31 1.25 1.51 1.53 1.93 2.22
Medical Physics, Vol. 42, No. 3, March 2015
of the DP-based segmentation scheme. Three experiments were performed on the 45 Grand Challenge cases, to explore the effects of these techniques in the proposed method. To examine the advantage of the likelihood image, we used the image intensity, instead of the likelihood image, for gradients computation in Eq. (6), and the other implementation of the segmentation framework remained the same. This is referred to as the method without PMM. The results in Table IV show that the mean Dice is 0.898, while the proposed method achieved 0.927; the mean RMSE is 3.37 mm, 0.91 mm worse than the proposed method; the R (0.85) = 0.98, while R (0.85) = 1 from the proposed method. These demonstrate the proposed method based on the likelihood image of the myocardium can achieve better performance than the method without PMM. In the second experiment, we performed the segmentation using a fixed search range, 20 mm, referred to as the method without ASR. The epicardium segmentation results in Table IV show that the mean Dice, mean RMSE, and R (0.85) are 0.880, 3.70, and 0.80 mm, respectively. It is apparent that the fixed search range technique with 20 mm is inferior to the proposed method using the adaptive search range. Finally, the method using the classic DP scheme without the connectivity constraint was examined. This is referred to as the method without CCDP. Its performance, shown in Table IV, is as follows: The mean Dice is 0.920, 0.007 less than that of the proposed method, and the mean RMSE is 2.63 mm, 0.23 mm worse than the proposed method.
RMSE (mm)
45 cases
Mean
Std
Mean
Std
R
R (0.85)
Endo Epi
0.892 0.927
0.043 0.017
2.30 2.39
0.67 0.66
0.98 0.99
0.97 1
1432
Qian et al.: Segmentation of myocardium from cardiac MR
F. 8. Segmentation volumes of our method versus manual volumes: (a) Endocardium detection: the obtained correlation coefficient is 0.98; (b) epicardium detection: the obtained correlation coefficient is 0.99.
In Fig. 9(b), one can find that the reliability curve of the method without CCDP is below that of the proposed method, which confirms the effectiveness of the connectivity constraint for the DP scheme. 3.C.3. The results from the literature
We compared the results of the proposed method with the results reported in the literature9,39–41 in terms of the Dice and APD for the endocardium detection and epicardium segmentation on the validation dataset from the MICCAI 2009, in Table V. The statistics of the Dice for the proposed method are 0.886 ± 0.039 for the endocardium and 0.928 ± 0.016 for the endocardium; and the statistics of the APD of our method are also similar to the other methods. This demonstrates that our method obtains the acceptable performance in terms of the endo- and epicardium segmentation. For reference, Table VI also provides the segmentation results from the literature12,17,28,29,42 in terms of the RMSE for the myocardium segmentation. It should be noted that an objective interwork comparison would be difficult as the test data and implementation of the different works could be very different. 4. DISCUSSION The focus of this work is to explore the segmentation ability of the new DP framework for the epicardium in Medical Physics, Vol. 42, No. 3, March 2015
1432
F. 9. The reliability of our method: (a) evaluated for the endocardium and epicardium segmentation, respectively; (b) the methods without the proposed new techniques evaluated for the epicardium segmentation.
cardiac cine MRI. We proposed three new techniques to advance the DP framework, i.e., the likelihood image of the myocardium, the adaptive search range, and the connectivity constrained DP algorithm. Here, we further discuss these techniques. The main purpose of constructing the myocardial likelihood image is to suppress the gradients of nonmyocardium tissues in the gradient map for the cost function. The gradients of nonmyocardium tissues may increase the external cost for the DP and deteriorate the optimal searching results of the epicardial contour. The average Dice of the segmentation method without the likelihood image is 0.3 less than that of our proposed method, meaning the likelihood image can improve the performance of the DP framework.
T IV. Comparisons of the epicardium segmentation to show the advantages of the likelihood image of the myocardium (PMM), ASR, and CCDP algorithm. RMSE (mm)
Dice 45 cases
Mean
Std
Mean
Std
R
R (0.85)
Proposed method Method without PMM Method without ASR Method without CCDP
0.927 0.898 0.880 0.920
0.017 0.021 0.031 0.019
2.39 3.37 3.70 2.62
0.66 0.76 0.94 0.68
0.99 0.97 0.97 0.98
1 0.98 0.80 1
1433
Qian et al.: Segmentation of myocardium from cardiac MR
T V. Comparisons of the segmentation accuracy of other approaches in terms of the Dice and APD on the validation dataset of the MICCAI 2009 Grand Challenge. Dice Endo
APD (mm) Epi
Validation references
Mean
Std
Mean
Schaerer et al. (Ref. 9) Jolly et al. (Ref. 37) Marak et al. (Ref. 38) Feng et al. (Ref. 36) Ours
0.869 0.880 0.860 0.890 0.886
0.043 0.030 0.040 0.040 0.039
0.921 0.930 0.930 0.940 0.928
Endo
Epi
Std Mean Std Mean Std 0.020 0.020 0.010 0.020 0.016
2.97 2.44 3.00 1.93 2.01
0.38 0.62 0.59 0.37 0.50
3.14 2.05 2.60 1.64 1.98
0.33 0.59 0.38 0.42 0.51
The adaptive search range is introduced to exclude the irrelevant tissues, thereby reducing the risk of undesirable gradients in the searching range. Whereas, the fixed search range, in order to ensure the complete coverage of the myocardium, has to be set a large size, such as 20 mm. Table IV shows the worse performance of the method using a fixed search range. This is because in the likelihood image, some tissues with similar intensity to the myocardium, such as the lung and liver, are rendered with high probability values, and their boundaries to the myocardium are also indistinct. Hence, these tissues have negative impacts on the optimal contour search of the DP when using the fixed search range. The connectivity constraint is incorporated into the DP algorithm to overcome a so-called Knapsack problem, i.e., the starting point does not connect with the ending point as shown in Fig. 5: Demonstration of the connectivity constrained dynamic programming: (a) the minimum path (red line) in the polar gradient map; (b) the minimum path in the newly formed image; (c) discontinuity of the outlined epicardium contour using the traditional method; and (d) the result using the proposed method (a). With the connectivity constrained DP algorithm, the epicardium contour is assured to be closed, as shown in Fig. 5. Demonstration of the connectivity
1433
constrained dynamic programming: (a) the minimum path (red line) in the polar gradient map; (b) the minimum path in the newly formed image; (c) discontinuity of the outlined epicardium contour using the traditional method; and (d) the result using the proposed method (d). Note that the performance of the proposed method in Table IV is just slightly better than that of the method without the connectivity constraint, because the Knapsack problem only happened in a small number of cases in the experiment. The DP scheme can be adopted and easily integrated into the myocardium segmentation framework, thanks to its two advantages. First, DP can flexibly impose the shape constraints of the endocardium to guide the epicardium segmentation; thereby we can overcome the difficulties caused by indistinguishable boundaries. Second, DP can make a local search for the desired boundary, such that we can reduce the adverse effects from the other tissues and obtain a desirable segmentation result. The main advantage of our segmentation method is that it only needs a minimized user-input (i.e., a seed point from a mouse click) and it is data-driven. Our method does not require any prior knowledge of the epicardium (such as shape information or local intensity statistics), instead this knowledge can be estimated from the MR images themselves for each phase. Therefore, the bias of prior knowledge derived from a limited training dataset can be avoided. A potential limitation of our segmentation scheme is that the epicardial segmentation is performed based on the endocardium segmentation, and an inaccurate endocardium segmentation result can affect the performance of the epicardium segmentation. Thus, we proposed an interaction for locating the center of the LV on the middle slice in the ED phase to pursue a reliable segmentation of the endocardium. As the result, our endocardium detection rate can reach 99.58%, based on which the detection rate of epicardium is 98.23%. Apparently, reliable endocardium segmentation can further improve the overall performance of our segmentation method. Finally, it is worth noting that the endocardium segmentation is not the emphasis of this work.
T VI. Comparisons with other methods in terms of the RMSE (mean ± std). References
Endo (mm)
Epi (mm)
Zhang et al. (Ref. 15) Zhu et al. (Ref. 12)
1.69 ± 0.38
1.89 ± 0.50
Active appearance model
50 subjects
0.69 ± 0.13
1.27 ± 0.18
22 subjects
Cordero et al. (Ref. 39) Cyrus et al. (Ref. 26)
1.37 ± 0.20
1.22 ± 0.17
Manual segmentation of ED phase and propagation using active shape model Deformable model
1.35 ± 0.13
2.05 ± 0.10
Ben Ayed et al. (Ref. 25)
1.60 (pixel)
1.99 (pixel)
Ours
2.30 ± 0.67
2.39 ± 0.66
Medical Physics, Vol. 42, No. 3, March 2015
Methodology
Manual segmentation first frame and propagation using convex relaxation Manual segmentation first frame and propagation using max-flow New dynamic programming
Data
43 subjects affected by acute myocardial infraction 20 subjects
20 subjects
45 subjects
1434
Qian et al.: Segmentation of myocardium from cardiac MR
The proposed method is implemented in using unoptimized code and takes about 25.9 s to segment one cardiac phase volume running on a personal computer equipped with a 2.4 GHz CPU.
5. CONCLUSION In this study, we presented a novel dynamic programming based method for myocardium segmentation from cardiac cine MRI. Our method only needs a minimized user-input (i.e., a seed point) and is data-driven. The knowledge of the epicardium (such as the shape information or local intensity statistics) is estimated for each phase individually; therefore no prior information from training data is required. The results demonstrated that the proposed method can obtain reliable and accurate myocardium segmentation, potentially providing a useful tool for academic research and developing new clinical applications of MRI for cardiology.
ACKNOWLEDGMENTS This work was supported in part by the National Natural Science Foundation of China (81301283), and in part by the SRF for ROCS, SEM of China (47th). a)Author
to whom correspondence should be addressed. Electronic mail: [email protected] 1D. L. Kasper, E. Braunwald, and A. Fauci, Harrison’s Principles of internal Medicine, 17th ed. (McGraw-Hill, New York, NY, 2008). 2V. Fuster, R. O’rourke, R. Walsh, and P. Pool-Wilson, Hurst’s the Heart, 12th ed. (McGraw-Hill, New York, NY, 2007). 3V. Tavakoli and A. Amini, “A survey of shaped-based registration and segmentation techniques for cardiac images,” Comput. Vis. Image Understanding 117, 966–989 (2013). 4X. Zhuang, “Challenges and methodologies of fully automatic whole heart segmentation: A review,” J. Healthcare Eng. 4, 371–408 (2013). 5C. Petitjean and J. Dacher, “A review of segmentation methods in short axis cardiac MR images,” Med. Image Anal. 15, 169–184 (2011). 6X. Lin, B. Cowan, and A. Young, “Automated detection of left ventricle in 4D MR images: Experience from a large study,” in Proceedings of the 9th International Conference on Medical Image Computing and ComputerAssisted Intervention(MICCAI), Nagoya, Japan (Springer, Berlin, Heidelberg, 2006), pp. 728–735. 7C. Cocosco, W. Niessen, T. Netsch, E. Vonken, G. Lund, A. Stork, and M. A. Viergever, “Automatic image-driven segmentation of the ventricles in cardiac cine MRI,” J. Magn. Reson. Imaging 28, 366–374 (2008). 8H. Lee, N. Codella, M. Cham, J. Weinsaft, and Y. Wang, “Automatic left ventricle segmentation using iterative thresholding and an active contour model with adaptation on short-axis cardiac MRI,” IEEE Trans. Biomed. Eng. 57, 905–913 (2010). 9J. Schaerer, C. Casta, J. Pousin, and P. Clarysse, “A dynamic elastic model for segmentation and tracking of the heart in MR image sequences,” Med. Image Anal. 14, 738–749 (2010). 10S. C. Mitchell, J. G. Bosch, B. P. F. Lelieveldt, R. J. van der Geest, J. H. Reiber, and M. Sonka, “3-D active appearance models: Segmentation of cardiac MR and ultrasound images,” IEEE Trans. Med. Imaging 21, 1167–1178 (2002). 11H. C. V. Assen, M. G. Danilouchkine, M. S. Dirksen, J. H. C. Reiber, and B. P. F. Lelieveldt, “A 3-D active shape model driven by fuzzy inference: Application to cardiac CT and MR,” IEEE Trans. Inf. Technol. Biomed. 12, 595–605 (2008). 12Y. Zhu, X. Papademetris, A. J. Sinusas, and J. S. Duncan, “Segmentation of the left ventricle from cardiac MR images using a subject-specific dynamical model,” IEEE Trans. Med. Imaging 29, 669–687 (2010). Medical Physics, Vol. 42, No. 3, March 2015
1434 13S.
Zhang, Y. Zhan, M. Dewan, J. Huang, D. N. Metaxas, and X. S. Zhou, “Towards robust and effective shape modeling: Sparse shape composition,” Med. Image Anal. 16, 265–277 (2012). 14S. Zhang, Y. Zhan, and D. N. Metaxas, “Deformable segmentation via sparse representation and dictionary learning,” Med. Image Anal. 16, 1385–1396 (2012). 15A. Andreopoulos and J. K. Tsotsos, “Efficient and generalizable statistical models of shape and appearance for analysis of cardiac MRI,” Med. Image Anal. 12, 335–357 (2008). 16S. Zambal, J. Hladuvka, and 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 ComputerAssisted Intervention(MICCAI), Copenhagen, Denmark (Springer, Berlin, Heidelberg, 2006), pp. 151–158. 17H. Zhang, A. Wahle, R. K. Johnson, T. D. Scholz, and M. Sonka, “4-D cardiac MR image analysis: Left and right ventricular morphology and function,” IEEE Trans. Med. Imaging 29, 350–364 (2010). 18M. Lorenzo-Valdés, G. I. Sanchez-Ortiz, A. G. Elkington, R. H. Mohiaddin, and D. Rueckert, “Segmentation of 4d cardiac MR images using a probabilistic atlas and the EM algorithm,” Med. Image Anal. 8, 255–265 (2004). 19X. Zhuang, K. S. Rhode, S. R. Arridge, R. Razavi, D. L. G. Hill, D. J. Hawkes, and S. Ourselin, “An atlas-based segmentation propagation framework using locally affine registration–application to automatic whole heart segmentation,” in Proceedings of the 11th International Conference on Medical Image Computing and Computer-Assisted Intervention(MICCAI), New York, NY (Springer, Berlin, Heidelberg, 2008), pp. 425–433. 20X. Zhuang, K. S. Rhode, R. S. Razavi, D. J. Hawkes, and S. Ourselin, “A registration-based propagation framework for automatic whole heart segmentation of cardiac MRI,” IEEE Trans. Med. Imaging 29, 1612–1625 (2010). 21W. Bai, W. Shi, D. P. O’regan, T. Tong, H. Wang, S. Jamil-Copley, and D. Ruechert, “A probabilistic patch-based label fusion model for multiatlas segmentation with registration refinement: Application to cardiac MR images,” IEEE Trans. Med. Imaging 32, 1302–1315 (2013). 22N. Paragios, “A variational approach for the segmentation of the left ventricle in cardiac image analysis,” Int. J. Comput. Vis. 50, 345–362 (2002). 23N. Paragios, “A level set approach for shape-driven segmentation and tracking of the left ventricle,” IEEE Trans. Med. Imaging 22, 773–776 (2003). 24J. Folkesson, E. Samset, R. Kwong, and C. Westin, “Unifying statistical classication and geodesic active regions for segmentation of cardiac MRI,” IEEE Trans. Inf. Technol. Biomed. 12, 328–334 (2008). 25I. B. Ayed, S. Li, and I. Ross, “Embedding overlap priors in variational left ventricle tracking,” IEEE Trans. Med. Imaging 28, 1902–1913 (2009). 26K. Li, X. Wu, D. Z. Chen, and M. Sonka, “Optimal surface segmentation in volumetric images-a graph-theoretic approach,” IEEE Trans. Pattern Anal. Mach. Intell. 28, 119–134 (2006). 27K. Li, S. Millington, X. Wu, D. Z. Chen, and M. Sonka, “Simultaneous segmentation of multiple closed surfaces using optimal graph searching,” in Proceedings of the 19th International Conference on Information Processing in Medical Imaging (IPMI), Glenwood Springs, CO (Springer, Berlin, Heidelberg, 2005), pp. 406–417. 28I. B. Ayed, H. M. Chen, K. Punithakumar, I. Ross, and S. Li, “Max-flow segmentation of the left ventricle by recovering subject-specific distributions via a bound of the Bhattacharyya measure,” Med. Image Anal. 16, 87–100 (2012). 29C. Nambakhsh, J. Yuan, K. Punithakumar, A. Goela, M. Rajchl, T. M. Peters, and I. B. Ayed, “Left ventricle segmentation in MRI via convex relaxed distribution matching,” Med. Image Anal. 17, 1010–1024 (2013). 30Y. Wu, Y. Wang, and Y. Jia, “Segmentation of the left ventricle in cardiac cine MRI using a shape-constrained snake model,” Comput. Vis. Image Understanding 117, 990–1003 (2013). 31A. Eslami, A. Karamalis, A. Katouzian, and N. Navab, “Segmentation by retrieval with guided random walks: Application to left ventricle segmentation in MRI,” Med. Image Anal. 17, 236–253 (2013). 32D. Geiger, A. Gupta, L. A. Costa, and J. Vlontzos, “Dynamic programming for detecting, tracking, and matching deformable contours,” IEEE Trans. Pattern Anal. Mach. Intell. 17, 294–302 (1995). 33A. Lalande, L. Legrand, P. M. Walker, F. Guy, Y. Cottin, S. Roy, and F. Brunotte, “Automatic detection of left ventricular contours from cardiac cine
1435
Qian et al.: Segmentation of myocardium from cardiac MR
magnetic resonance imaging using fuzzy logic,” Invest. Radiol. 34, 211–217 (1999). 34J. Fu, J. Chai, and S. Wong, “Wavelet-based enhancement for detection of left ventricular myocardial boundaries in magnetic resonance images,” Magn. Reson. Imaging 18, 1135–1141 (2000). 35U. Kurkure, A. Pednekar, R. Muthupillai, S. D. Flamm, and I. A. Kakadiaris, “Localization and segmentation of left ventricle in cardiac cine-MR images,” IEEE Trans. Biomed. Eng. 56, 1360–1370 (2009). 36T. H. Cormen, C. E. Leiseron, and R. L. Rivest, Introduction to Algorithms, 23rd ed. (McGraw-Hill, New York, NY, 1999). 37M. N. Ahmed, S. M. Yamany, N. Mohamed, A. A. Farag, and T. Moriarty, “A modified fuzzy c-means algorithm for bias field estimation and segmentation of MRI data,” IEEE Trans. Med. Imaging 21, 193–199 (2002). 38W. Shi, X. Zhuang, H. Wang, S. Duckett, D. Oregan, P. Edwards, and D. Rueckert, “Automatic segmentation of different pathologies from cardiac cine MRI using registration and multiple component EM estimation,” in Proceedings of Functional Imaging and Modeling of the Heart (Springer, Berlin, Heidelberg, 2011), pp. 163–170.
Medical Physics, Vol. 42, No. 3, March 2015
1435 39C.
Feng, C. Li, D. Zhao, C. Davatzikos, and 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(MICCAI) Nagoya, Japan (Springer, Berlin, Heidelberg, 2013), pp. 477–484. 40M. P. Jolly, “Fully automatic left ventricle segmentation incardiac cine MR images using registration and minimum surfaces,” in MICCAI 2009 Workshop on Cardiac MR Left Ventricle Segmentation Challenge, MIDAS J. (2009). 41L. Marak, J. Cousty, L. Najman, and H. Talbot, “4D morphological segmentation and that MICCAI LV-Segmentation grand challenge,” in MICCAI 2009 Workshop on Cardiac MR Left Ventricle Segmentation Challenge, MIDAS J. (2009). 42L. Cordero-Grande, G. Vegas-Sánchez-Ferrero, P. Casaseca-de-la-Higuera, J. Alberto San-Román-Calvar, A. Revilla-Orode, M. Martín-Fernández, and C. Alberola-López, “Unsupervised 4D myocardium segmentation with a markov random field based deformable model,” Med. Image Anal. 15, 283–301 (2011).
