Improved graph cut model with features of superpixels and neighborhood patches for myocardium segmentation from ultrasound image

1.   Introduction

Myocardial infarction (MI) is a severe cardio-vascular disease that threatens human health, leading to the subsequent death of cardiomyocytes and vascular cells in the vicinity site of the infarction. Ultrasound cardiogram provides a non-invasive way to diagnose and monitor heart conditions. In the typical B-mode ultrasound (US) images, the structure of the heart and the morphology of the myocardium can be checked. For example, the left ventricle, an important pumping chamber of the heart, can be scanned in the parasternal long-axis (PLAX) and parasternal short-axis (PSAX) views showing the cavity (the hypoechoic area) and the myocardial tissue (isoechoic area surrounding the cavity). However, some inherent drawbacks of US imaging, such as low contrast, speckle noise, signal dropout, acoustic shadow, cause the myocardial tissue indistinguishable from the background. It is challenging to investigate accurate and effective segmentation algorithms of myocardium ultrasound (MUS) images. Segmentation of target tissue from other tissues or background is an essential phase of ultrasound computer-aided diagnosis ^[1]. Myocardial segmentation and assessment are the keys in the morphology and function study of myocardium ^[2].

Previous studies focused on the endocardial segmentation of MUS images. Tao et al. ^[3] proposed an evolution strategy called as "tunneling descent", which was capable of escaping from spurious local minima and used in the maximum-a-posterior (m.a.p) active contour segmentation of the endocardium. Hamou et al. ^[4] proposed an external energy for a gradient vector flow (GVF) snake. The external energy could provide additional information to the active contour model to segment the endocardium. However, there are relatively less studies on the segmentation of the epicardium or full myocardium in MUS images. Zhu et al. ^[5] proposed a coupled deformable model to extract the intact myocardial tissues. They evolved the endocardial and epicardial surfaces simultaneously while maintaining the coupling of the myocardial volume. In spite of the sophistication of the model, it is still challenging for segmentation in the low-quality MUS images with low contrast and fuzzy myocardial contours.

Many approaches have been developed for US image segmentation. Huang et al. ^[6] reviewed the literatures summarizing segmentation techniques into several classes such as thresholding-based, clustering-based, watershed-based, graph-based, active contour model, Markov random field and neural network. Neural network or deep learning has become a hot research topic. Some researchers investigated the semantic segmentation algorithms ^[7,8]. Others applied the traditional machine learning methods combined with prior information and intuitive feeling ^[9,10].

Notably, graph cut model as a traditional machine learning method and classic semi-supervised segmentation method for natural images ^{[11,12,13,14]} has developed well in US image segmentation ^{[15,16,17,18,19,20,21,22]}. The model is a flexible framework where appropriate image features and similarity measure can be selected to improve the performance of the segmentation. In addition, interactive labels provided by user determine the source terminal (S) and sink terminal (T) of the graph and ensure the robustness of the segmentation.

The graph cut model is an energy minimization problem which can be efficiently solved via optimization techniques ^{[23,24,25,26]}. In 2014, Ciurte et al. ^[19] proposed a semi-supervised segmentation method of fast Neighborhood Patches based Continuous Min-Cut (fP-CMC) to segment the US images of prostate, fetus head, and tumors of the liver and eye. In their study, neighborhood patches ^{[27,28,29,30]} were selected as image features and Pearson distance ^[31] was selected as similarity measure. Neighborhood patches as promising image features well represented and dealt with the textures to help reduce the effect of speckle in US images ^[14]. Pearson distance based on the Bayesian formulation was made under the assumption of the general speckle model. It was shown to be robust to the speckle noise presented in US images, and the higher performance in distinguishing different tissues compared with l₂ norm. Then, the continuous graph cut model was adopted and solved via a fast minimization algorithm based on augmented Lagrangian and operator splitting ^[32,33]. The good results were obtained. However, the fP-CMC method exists the following limitations:

(1) Due to the overlap of the neighborhood patches, a high probability that the two neighbor pixels were incorrectly assigned to the same class but they belong to different classes, causes the segmentation failure of the target contour with abnormal roughness. Therefore, only the neighborhood patches as image features is not sufficient, and it is desirable to incorporate with other features.

(2) The similarity measured by Pearson distance is easily affected by the brightness variation. In the case of the target contour partially shadowed by acoustic shadow, a good segmentation result would not be obtained.

(3) The graph-based method reduces computation cost with the expense of down-sampling. Consequently, the target contour is segmented with low resolution.

To solve the limitations of neighborhood patches, this study proposes a semi-supervised segmentation method of fast Superpixels and Neighborhood Patches based Continuous Min-Cut (fSP-CMC). Superpixels ^[34] are the patches consisting of neighbor pixels with similar intensity, which are similar to the neighborhood patches, but have the different attractive points of shape irregularity and no overlap with each other. Therefore, this study applies the superpixels and incorporates the location information of each superpixel as coarse features to constraint the segmented contour. Thus, the one novelty of our fSP-CMC method is the complementary combination of the location information of superpixels as coarse features with the grayscale information of neighborhood patches as fine texture features.

Pearson correlation coefficient measuring the correlation between the features is irrelevant to brightness variation. This study defines a novel similarity measure by combination with Pearson distance and Pearson correlation coefficient to capture and enhance the features correlation even for the images with low quality. This study expects to obtain more exact similarity and thus further improve the segmentation performance. That is the other novelty of fSP-CMC. Finally, to solve the low resolution problem of the segmented contour, Non-Uniform Rational B-Spline (NURBS) curve fitting ^[35,36,37] is applied to restore the original resolution.

2.   Materials and methods

2.1.   Data acquisition

The 20 Sprague-Dawley (SD) rats with normal heart and myocardial infarction were used as experimental subjects in this study. After anesthesia, the rat was laid in supine position on the physiological information monitoring platform. The four soles of the rat were fixed on the copper sheet by the conductive adhesive and the tape. Its electrocardiogram (ECG) signal was recorded and its body temperature was maintained by the platform.

An US imaging system with ultra-high resolution for small animal experiments (Vevo2100, Visual Sonics, Canada) was used to scan the heart of the rat. The hair at the imaging target area were removed by the depilatory paste. As shown in Figure 1, the PLAX and PSAX views of left ventricle were scanned, using the MS-250 probe with a transmit frequency of 21 MHz. The resolution of the imaging system was approximately 70 μm. The physical size of image acquisition was 21.00 mm for depth and 18.04 mm for width. The frame rate of acquisition was 128 Hz. The MUS images with a size of 1000 × 1000 pixels were reconstructed from IQ (in-phase Quadrature) data and enhanced by gray stretch method ^[38].

Figure 1.  The PLAX view (a) and PSAX view (c) of left ventricle of a rat scanned using MS-250 probe. (b) A typical MUS image in PLAX view. (d) A typical MUS image in PSAX view.

DownLoad: Full-Size Img PowerPoint

This study collected 200 B-mode MUS images of left ventricle of the rats (5 images of PLAX and PSAX views, respectively, for each rat). The five images of each view were sequentially recorded every two frames from End-Diastole (ED) to End-Systole (ES) in one cardiac cycle. These MUS images were processed to validate the accuracy and effectiveness of the proposed fSP-CMC method.

Before the animal experiments, ethical approval was obtained from the Medical Ethics Committee of Nanfang Hospital, Southern Medical University, China. Experiments on rats were performed in accordance with the Guidelines for the Care of Laboratory Animals of the National Institutes of Health.

2.2.   The proposed fSP-CMC segmentation method

Figure 2 shows the flow chart of this study including the proposed fSP-CMC segmentation method, which contains (1) graph construction of graph cut model, (2) feature extraction of superpixels and neighborhood patches, (3) definition of novel similarity measure, (4) setting of interactive labels, and (5) solution of continuous graph cut model.

Figure 2.  Flow chart of this study using the proposed fSP-CMC method to segment the myocardium in MUS image.

DownLoad: Full-Size Img PowerPoint

2.2.1.   Graph construction of graph cut model

Given the MUS image I, a weighted undirected graph G = (V, E, W) is constructed. The node set V is described as {x_i, i∈{1, 2, ..., N}}, where each node x_i of the graph G corresponds to a pixel in the image I and N is the number of the nodes in the graph G. The undirected edge set E = {(x_i, x_j), i, j∈{1, 2, ..., N}} connects two nodes, x_i and x_j. The similarity between x_i and x_j is measured by the edge weight W(x_i, x_j) (Equation 1).

Where d(x_i, x_j) measures the difference between x_i and x_j inversely proportional to W(x_i, x_j) and depends on the image features in this study. A low value of d(x_i, x_j) means high similarity and high weight value W(x_i, x_j) between x_i and x_j, and the two nodes may belong to the same class with a high probability. N_n(x_i) is the searching window centered at x_i with a size of n × n. It is worth noting that too narrow window (e.g., n = 1 indicating that the similarity is measured only from the nearest nodes) is not robust enough to the speckle noise for lack of statistical information. Contrarily, too wide window puts pressure on the computation. In this study, n = 5 to ensure a good compromise between nonlocality and efficiency. σ is the scaling parameter of the W(x_i, x_j) and σ² = 7e−3 for high sensitivity.

A graph cut model actually cuts two subsets A and A^c from the node set V of the graph G (s.t. A∪A^c = V, A∩A^c = Ø). A and A^c represent the target region and the background region respectively. The cost of the cut is defined as Equation 2, equal to the sum of edge weight W(x_i, x_j) (Equation 1) between two neighbor nodes x_i and x_j, where x_i∈A, x_j∈A^c, and x_j∈N_n(x_i). The minimum cut (Min-Cut) is obtained when the cut cost reaches minimum, i.e., the total similarity between A and A^c reaches weakest.

2.2.2.   Feature extraction of superpixels and neighborhood patches

To well construct the graph G, the d(x_i, x_j) in Equation 1 is essential to be solved. In this study, superpixels and neighborhood patches are chosen as image features. Figure 3 shows the sketch map of graph construction with the superpixels and neighborhood patches. A superpixel is larger than a neighborhood patch and thus contains more information. Due to the attractive points of shape irregularity of superpixels and no overlap with each other, the texture patterns of superpixels are different from those of neighborhood patches. Therefore, in this study, the coarse features from superpixels are complementarily combined with the fine features from neighborhood patches.

Figure 3.  Sketch map of graph construction with superpixels and neighborhood patches.

DownLoad: Full-Size Img PowerPoint

To generate superpixels, simple linear iterative clustering (SLIC) ^[39,40,41], adapted from k-means clustering, is used in this study. Different from k-means clustering in the whole image region, the search region of SLIC is restricted in a 2C × 2C neighborhood, where C is the distance between neighbor cluster centers sampled initially. Then, the algorithm convergence can be speeded up.

Each cluster center is carried with intensity and location information, and compared with all pixels in its 2C × 2C neighborhood by the dedicated distance measure D ^[40]. Each pixel is assigned to the cluster with the minimum D value, and each cluster center with the information follows to be updated. By iterative process, final shaped clusters are the generated superpixels. Then, two coarse features (i.e. location information) including the row and column values of the center of each superpixel are assigned to the inner pixels.

To extract fine features, neighborhood patches with a size of n_p × n_p centered on each pixel are generated. The n_p² grayscale values (i.e., grayscale information) of each neighborhood patch are assigned to the centered pixel as its fine features. n_p is linked to the integrity of texture pattern and preservation of details information. In this study, n_p = 3 to ensure that the each patch contains a complete texture pattern and preserves details information ^[42]. So, nine fine features are extracted from neighborhood patches.

2.2.3.   Definition of novel similarity measure

For each node x_i, the corresponding feature vector with a length of F_num consists of the aforementioned coarse and fine features, F(x_i) = (F¹(x_i), ..., F^F_num(x_i)). F₁_num is the number of coarse features and F₂_num is the number of fine features. Thus, F_num = F₁_num + F₂_num. Then, a novel similarity measure d(x_i, x_j) is defined as Equation 3, where Pd(x_i, x_j) represents Pearson distance (Equation 4) measuring the difference between F₂_num fine features, while Pc(x_i, x_j) represents Pearson correlation coefficient (Equation 5) measuring the correlation between F_num features. In this study, F₁_num = 2, F₂_num = 9, and F_num = 11.

In Equation 3, p₁ and p₂ are the weighting factors of Pd(x_i, x_j) and Pc(x_i, x_j). p₁ > 0 and p₂ > 0 to make sure log(Pc(x_i, x_j) + p₂) > 0. Increasing p₁ value or decreasing p₂ value can increase the weight of Pd(x_i, x_j), while decreasing p₁ or p₂ value can increase the weight of Pc(x_i, x_j). The determination of the values of p₁ and p₂ depends on the contour of the target tissue. For a target tissue such as myocardium in a fuzzy background with relatively smooth contour, p₁ should decrease, otherwise high weight of Pd(x_i, x_j) is easily affected by acoustic shadow. Here it is not recommended to increase p₂, because low weight of Pc(x_i, x_j) results in the segmentation contour with abnormal roughness. Too large p₂ value causes no effect of Pd(x_i, x_j) and Pc(x_i, x_j) on d(x_i, x_j). So p₁ = 0.2 and p₂ = 10 in this study.

The solved d(x_i, x_j) is substituted in Equation 1, then the edge weight W(x_i, x_j) is obtained.

2.2.4.   Setting of interative labels

Our pilot study showed that a better segmentation result can be obtained using the enclosed inner and outer boundaries due to the larger areas of source terminal (S) and sink terminal (T) for the enclosed boundaries in comparison with the non-enclosed boundaries. Therefore, as shown in Figure 4, an enclosed inner boundary is roughly sketched inside the myocardial tissue and an enclosed outer boundary is initialized inside the background surrounding the myocardial tissue in a PLAX-viewed image. Then, y(x_i) is defined as a binary labeling function. In this study, y(x_i) is taken on a continuous domain [0, 1] (i.e., a convex set) to relax the binary constraint. The nodes inside the inner boundary assigned to the S are labeled to 1 and y(x_i) = 1. The nodes outside the outer boundary assigned to the T are labeled to 0 and y(x_i) = 0. The others between the inner and outer boundaries are unassigned (Figure 4b). These unassigned nodes are initially labeled to random values in (0, 1). Then the discrete graph cut model is transformed to the continuous one. Accordingly, the subsequent solution problem becomes computationally tractable, otherwise it would be hard to proceed due to non-convexity. For a PSAX-viewed image, the initial drawing is performed twice. First, the enclosed inner and outer boundaries are drawn inside the heart cavity and the myocardial tissue, respectively, for extraction of the endocardium. The second enclosed inner and outer boundaries are then drawn inside the myocardial tissue and the background, respectively, for extraction of the epicardium.

Figure 4.  Sketch maps of semi-supervised segmentation based on constructed graph G and the drawing of initial boundaries in a typical PLAX-viewed image. (a) Rough sketches of enclosed inner and outer boundaries manually initialized. (b) The assigned nodes with the interactive labels and the unassigned nodes. (c) Minimum cut (Min-Cut) obtained by solving continuous graph cut model. (d) The initial inner and outer boundaries in a typical PLAX-viewed image. Red contour represents the inner boundary. Green contour represents the outer boundary. Note: In a, b, c, the edges are omitted to show a clear illustration. The graph consists of nodes and edges.

DownLoad: Full-Size Img PowerPoint

2.2.5.   Solution of continuous graph cut model

After two terminals S and T of the constructed graph G are set, the segmentation is automatically performed by solving continuous graph cut model. The cost function of cut (Equation 6) consists of data item E_data(y) and smooth item E_smooth(y). β is a smoothing parameter and helps to eliminates the misclassification of small subsets of nodes. β∈[1e-4, 1e-1] is adapted to the segmentation of US images ^[19], thus β = 5e−2 in this study.

E_data(y) is defined as a measure of the similarity between A and A^c, i.e., the cut cost (Equation 2). It is also equivalently defined as Equation 7 using a non-local H¹ (NL-H¹) operator ^[14,19,43], the square of l² norm of the continuous graph gradient of labels y.

E_smooth(y) (Equation 8) measures the discontinuity of y between A and A^c, and is used to enforce spatial coherency by penalizing. However, over-penalizing should be avoided (i.e., discontinuity preserving ^[23]). Thus, total variation (TV) ^[44], the l¹ norm of the gradient of y, is chosen as smooth priors. Then, the cost function (Equation 6) is rewritten as Equation 9.

The minimization of Equation 9 is performed by a fast minimization algorithm ^[19] (Algorithm 1) based on the splitting step (Equation 10) and the augmented Lagrangian step (Equation 11). v and q in Equations 10, 11 are the coupling terms forced to be almost equal to y. λ₁^k, λ₂^k, r₁, and r₂ in Equation 11 are the penalty parameters. In each iteration of the minimization process, each unassigned node either keeps the current label or is re-assigned a new label. A and A^c are updated according to the updated labels y^k. The minimum of the cost function (Equation 9) is obtained when the residual error between the updated labels y^k and last labels y^k−1 is less than or equal to the threshold ε.

2.3.   Post-processing of NURBS curve fitting

The NURBS curve fitting method is used to restore the original resolution of segmented contour in this study. NURBS curve C(u), a parametric approximation curve, is defined as Equation 12.

where N_{i, p}(u), i∈{0, 1, ..., n − 1} are the p-degree B-spline basis functions defined as Equation 13 and n is the number of pixels on the segmented contour.

P_i, i∈{0, 1, ..., n − 1} are control points corresponding to the pixels on the segmented contour, ω_i, i∈{0, 1, ..., n − 1} are the weights of the control points. u_i, i∈{0, 1, ..., n + p} are the knots constituting a knot vector U defined as Equation 14, the U determines how the control points affect the C(u).

In Equation 14, the values of the first (p + 1) knots are equal to a (0 by default), the values of the last (p + 1) knots are equal to b (1 by default), and the middle (n – p − 1) knots are assigned certain values. The values in the U are subject to two conditions. One is that the values remain non-decreasing from a to b, the other is that the number of repetitions (called Multiplicity) for each value is not allowed greater than (p + 1). The knot with Multiplicity equal to (p + 1) is called Full-Multiplicity Knot (e.g., a and b). The knot with Multiplicity equal to 1 is called Simple Knot.

When entering a Simple Knot, a new control point is activated while the old is deactivated, and then the B-spline basis function would interpolate between the two control points. When entering a knot with Multiplicity larger than 1, multiple control points are activated at the same time, and the B-spline basis function would not interpolate between the multiple control points. This affects the continuity of the C(u).

In this study, p = 2, ω = 1. The knot vector U begins with Full-Multiplicity Knots (a), followed by Simple Knots whose values subject to an arithmetic sequence starting from u_p+1 = 1/(n−p) with an increment of 1/(n−p), and ends with Full-Multiplicity Knots (b) to interpolate and well enhance the continuity of the C(u).

2.4.   Segmentation evaluation

In this study, the segmentation results of the fSP-CMC are evaluated on the agreement with the ground truth (GT) by medical expert and compared with the results of the fP-CMC ^[19].

As shown in Figure 5, C_Seg = {Seg₁, Seg₂, …, Seg_P} represents the set of pixels on the contour of the myocardium by the segmentation method. C_GT = {GT₁, GT₂, …, GT_Q} represents the set of pixels on the contour of the myocardium by the GT. P and Q are the number of pixels on C_Seg and C_GT, respectively. The region enclosed by C_Seg is Ω_Seg, equal to TP + FP. The region enclosed by C_GT is Ω_GT, equal to TP + FN. TP denotes the true positive pixels in the truly segmented myocardial region. TN denotes the true negative pixels in the truly segmented background region. FP denotes the false positive pixels in the falsely segmented myocardial region. FN denotes the false negative pixels in the falsely segmented background region.

Figure 5.  Sketch map of two contours C_Seg (red contour) and C_GT (blue contour) from the segmentation result and the ground truth (GT), respectively. TP, TN, FP, and FN denote the true positive pixels, the true negative pixels, the false positive pixels, and the false negative pixels, respectively.

DownLoad: Full-Size Img PowerPoint

Then, three parameters are calculated for segmentation evaluation. Dice coefficient ^[45] is defined as Equation 15 reflecting the agreement on the segmented region by the segmentation method with the GT. The value of Dice ranges from 0 to 1. 1 represents full agreement.

The mean absolute distance (MAD) ^[46] (Equation 16) and Hausdorff distance (HD) ^[47] (Equation 17) are calculated to evaluate the agreement on contours of the segmentation result with the GT.

Where d(Seg_i, C_GT) (Equation 18) represents the shortest distance from the pixel Seg_i to the contour C_GT and d(GT_i, C_Seg) (Equation 19) represents the shortest distance from the pixel GT_j to the contour C_Seg. The value of MAD reflects the global agreement between two contours C_Seg and C_GT, while HD reflects the local agreement between the two contours.

One-way analysis of variance (ANOVA) ^[48] is used to test the statistical difference in Dice, MAD, and HD for two methods and two views, respectively. Statistical significance is considered at p < 0.05. In addition, box-plot is used to present the central trend and variation of these three parameters.

Additionally, this study used the receiver operating characteristic (ROC) curve ^[49] with the area under the ROC curve (AUC) to evaluate the segmentation (pixelwise classification) performance of single image. The predicted values of all the pixels obtained by a segmentation method are ranked in ascending order of value. The ROC curve starts as the discrimination threshold is equal to the smallest predicted value. FP rate (FPR) and TP rate (TPR) are calculated. FPR and TPR are defined as Equation 20. Then, the discrimination threshold changes to the second predicted value. FPR and TPR values are calculated again. This procedure repeats till all the predicted values are used as discrimination threshold. After ROC curve is obtained, AUC is calculated to quantify the segmentation ability expressed by the curve. The larger the AUC value, the better the segmentation ability. AUC∈(0.5, 1) suggests that the segmentation method is with good segmentation ability.

Finally, Bland-Altman plot ^[50] with the average of MA (myocardial area) vs. the difference of MA is plotted to overall evaluate the agreement of the segmentation method with the GT. The 95% confidence interval is (−1.96SD, 1.96SD).

The proposed fSP-CMC method compared with the fP-CMC method was implemented using MATLAB R2016b and Visual Studio 2015 hybrid programming on a personal computer with 3.30GHz Intel (R) Core i5-4590 CPU and 8GB RAM.

3.   Results

The myocardial segmentation results of the proposed fSP-CMC method were in high agreement with the GT results in the PLAX-viewed and PSAX-viewed images of left ventricle of rats (Figure 6). The deformed myocardial tissues could be extracted during systole of a cardiac cycle. Moreover, Figure 6 shows that the fSP-CMC obtained smoother segmentation contours in comparison with the fP-CMC.

Figure 6.  The MUS images and the corresponding segmentation results of left ventricle of a rat. a, d: The images in PLAX and PSAX views respectively. From left to right, the images were collected from End-Diastole (ED) to End-Systole (ES) in one cardiac cycle. b, e: The segmentation results of the fP-CMC. c, f: The segmentation results of the fSP-CMC. In b, c, e, and f, yellow contours represent the segmentation results by the fP-CMC or fSP-CMC, while blue areas indicate the segmentation results of the GT.

DownLoad: Full-Size Img PowerPoint

3.1.   ROC analysis with AUC

For each MUS image, the agreement of fSP-CMC with the GT in comparison with fP-CMC was evaluated using ROC curve with the AUC. Figure 7 shows the ROC curves for the MUS images in Figure 6 at ED moments. It was found that the ROC curves obtained from the fSP-CMC segmented images were located closer to the upper left corner indicating a higher agreement in the segmentation results with the GT. The AUC values of the fSP-CMC were higher than those of the fP-CMC, especially for the PLAX-viewed image in which part of myocardium is darkened by acoustic shadow (0.9610 > 0.9456). Figure 8 shows the variation of AUC values for all the PLAX-viewed and PSAX-viewed MUS images. The AUC results indicate that the fSP-CMC method is more stable and its segmentation results are with the smaller fluctuation.

Figure 7.  ROC curves to evaluate the segmentation ability of fP-CMC and fSP-CMC methods for the PLAX-viewed (a) and PSAX-viewed (b) MUS image at ED moment.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Line charts of AUC for 100 PLAX-viewed (a) and 100 PSAX-viewed (b) MUS images segmented by the fP-CMC (blue) and fSP-CMC (red).

DownLoad: Full-Size Img PowerPoint

3.2.   Quantitative analysis of Dice, MAD, and HD with statistical difference and box-plots

Table 1 shows the quantitative results of Dice, MAD, and HD with the statistical significant difference at level p < 0.05 and Figure 9 shows the box-plots of the three parameters.

Table 1.  Comparative results of Dice, MAD, and HD between the fP-CMC and fSP-CMC. The given values are expressed as Mean ± SD.

	fP-CMC				fSP-CMC
	Dice	MAD (mm)	HD (mm)		Dice	MAD (mm)	HD (mm)
PLAX	0.899 ±0.025	0.245 ±0.047	1.052 ±0.272		0.910 ±0.021 a	0.223 ±0.039 a	0.890 ±0.224 a
PSAX	0.865 ±0.030 b	0.266 ±0.092	0.892 ±0.200 b		0.900 ±0.024 a, b	0.194 ±0.040 a, b	0.665 ±0.181 a, b
a Statistically significant difference at level p < 0.05 vs. the fP-CMC method. b Statistically significant difference at level p < 0.05 vs. the PLAX view.

---

 | Show Table

DownLoad: CSV

Figure 9.  Box-plots of Dice (a), MAD (b), and HD (c) using the fP-CMC and fSP-CMC with respect to PLAX and PSAX views. Central red mark is the median, edges of the box are the 25th and 75th percentiles, whiskers extend to extreme values.

DownLoad: Full-Size Img PowerPoint

3.2.1.   Comparison of fSP-CMC with fP-CMC

Normally, a higher Dice is associated with a lower MAD and HD. In both PLAX and PSAX views, it was found that the segmentation of fSP-CMC had a superior performance. The Dice was significantly higher than that of fP-CMC (0.910 > 0.899 for PLAX, 0.900 > 0.865 for PSAX, on average, p < 0.05). The MAD and HD were significantly lower than those of fP-CMC (0.223 mm < 0.245 mm for PLAX, 0.194 mm < 0.266 mm for PSAX; 0.890 mm < 1.052 mm for PLAX, 0.665 mm < 0.892 mm for PSAX, on average, p < 0.05).

3.2.2.   Comparison of PLAX with PSAX

Parameters Dice, MAD, and HD reflect the agreement of region, global contour, and local contour, respectively. As shown in Table 1, for fSP-CMC, a higher Dice was obtained in the PLAX view, whereas lower MAD and HD values in the PSAX view. The results indicated that the fSP-CMC method had a superior performance in the region agreement for segmentation of the PLAX-viewed MUS images but a superior performance in the contour agreement for the segmentation of the PSAX-viewed MUS images. The results may be associated with the shape regularity of the target tissue. It is noted that the myocardium in the PLAX-viewed images is of an irregular horseshoe shape while the myocardium in the PSAX-viewed images is of a regular ring shape.

3.3.   Overall analysis of Bland-Altman

As shown in Figure 10b, the mean of the difference values between the fSP-CMC segmented results and the GT was 4.05 mm² and only 5% of the difference values were beyond the 95% confidence interval (−4.32 mm², 11.57 mm²). In comparison with the Bland–Altman plot by fP-CMC (Figure 10a), the more difference values between the fSP-CMC segmented results and the GT were within the 95% confidence interval and the mean value is closer to zero. A higher agreement between the segmentation result and GT suggests that the proposed fSP-CMC provided more consistent segmentation with the GT.

Figure 10.  Bland-Altman plots of myocardial area (MA) for all 200 MUS images segmented by the fP-CMC (a) and fSP-CMC (b). The middle horizontal red line represents the mean of the difference values (blue dots). The upper and lower horizontal green lines represent the upper and lower limits of 95% confidence interval (−1.96SD, 1.96SD).

DownLoad: Full-Size Img PowerPoint

4.   Discussion

4.1.   Contribution of NURBS curve fitting

Because the graph-based method reduces the resolution of the images, the target tissues would not be well segmented ^[19]. NURBS curve fitting has been extensively used in Computer Aided Design to clearly describe the complex contour parts ^[51]. This study firstly applies NURBS curve fitting as post-processing of up-sampling to remedy the segmented contours with the reduced resolution. The results show its feasibility.

4.2.   Improved performance profiting from the superpixels

The results of the validation show that the proposed fSP-CMC has significant improvement compared with fP-CMC and has potentials to be another reliable segmentation method to segment the full myocardium in the PLAX and PSAX-viewed US images of left ventricle of the small animals. In the previous study of Reference ^[19], the fP-CMC method only applied the grayscale information of neighborhood patches as features and TV (i.e., the l¹ norm of the gradient of labels y) as the smooth priors. However, it was insufficient to solve the problem of unsmoothness caused by the graph cut model. If let the l² or higher norm as the smooth priors, the abnormal smoothness with loss of useful information would be induced ^[44]. In comparison with fP-CMC, the fSP-CMC proposed in this study incorporates the location information of superpixels (coarse features) to constrain the contour of the full myocardium extracted. Not only is the desired smoothness achieved, but also the image information is not lost (Figure 6).

4.3.   Improved performance profiting from the novel similarity measure

This study proposes the fSP-CMC method with the novel similarity measure, which combines Pearson correlation coefficient and Pearson distance. The correlation between the all features is measured by Pearson correlation coefficient, while the difference between the fine features is measured by Pearson distance. By balancing the weights of Pearson correlation coefficient and Pearson distance, the more exact similarity is obtained and the performance of the fSP-CMC segmentation method is improved. As Figure 7a shown, fSP-CMC performs better than fP-CMC for the US image in which part of myocardium is darkened by acoustic shadow. The fSP-CMC method not only extracts the myocardial tissues, but could also flexibly and satisfactorily segment different targets, for example, atherothrombotic plaques and heterotopic ossification tissues ^[52].

4.4.   Subsidiarity of the interactive labels

In this study, the interactive labels were set to determine the terminals S and T of the graph G (Figure 4). The process of adding or subtracting these labels corrects segmentation imperfections by reset. During the processing, there is no need to reconstruct the graph G (Figure 2), and the performance of the method is not affected. If the segmented region does not completely include the target region, the source terminal (S) could be compensated by adding the corresponding labels. If the segmented region includes the undesired surrounding tissues, the sink terminal (T) could be enlarged in the same way. Additionally, the S or T could be updated by subtracting the labels of the wrongly assigned nodes. Therefore, the graph-based segmentation method is reproducible and particularly suited for the tissue segmentation of the US images with a poor quality. Furthermore, with the easy and quick reset, the method could be used in two-region segmentation, as well as in multi-region segmentation.

4.5.   Limitation of fSP-CMC

Using the proposed fSP-CMC method, the segmentation results might be affected by three factors. (1) One factor is the quality of MUS image. In this study, fuzzy myocardial contour and presence of acoustic shadow in most of PSAX-viewed images affected the segmentation performance in region agreement (Dice), which was inferior to the PLAX-viewed images. (2) The second factor is the shape regularity of the target tissues. The fSP-CMC method possesses superior performance with lower values MAD and HD on the PSAX-viewed images in global and local contour agreement to the PLAX-viewed images. (3) The segmentation results might be affected by manual interaction. Fortunately, the results could be improved by the reset of interactive labels.

Additionally, the proposed fSP-CMC method is semi-supervised and requires the initial boundaries initialized in proper regions. The initial inner and outer boundaries could be properly drawn in the myocardial tissue and the background in the first frame. However, these two boundaries might invade into inappropriate regions due to the motion of the myocardium in the following frames, causing incorrect segmentation of the myocardium. To solve this problem, the more insight exploration will focus on automatic contour tracking in combination with fSP-CMC. With the tracked contour, it is able to approximately estimate the two boundaries around the myocardial contour in the following frames.

Although this study scanned the heart of the small animal, which is smaller and beating faster than human heart, it is more challenging to segment the tiny myocardial tissues from ultrasound images. The good results of segmentation obtained in this study suggest the application potentials of the proposed method in human myocardium extraction.

5.   Conclusion

This study presents a method of fSP-CMC based on a continuous graph cut model for segmentation of the full myocardium from the MUS images. The complementary combination of superpixels and neighborhood patches and the definition of the novel similarity measure improve the segmentation performance of the fSP-CMC. The validation results of the MUS images indicate that the proposed fSP-CMC method with superpixels has superior performance to the fP-CMC method with no superpixels in aspects of contour smoothness, segmentation accuracy and effectiveness.

Acknowledgments

This work is supported by National Natural Science Foundation of China [No. 81371560], Guangdong Science and Technology Program [No. 2016A020216017 and No. 2013B021800039], and Guangdong Natural Science Foundation [No. 2014A030313329].

Conflict of interest

All authors declare no conflicts of interest in this paper.