Computers in Biology and Medicine 55 (2014) 1–10

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Segmentation of anterior cruciate ligament in knee MR images using graph cuts with patient-specific shape constraints and label refinement Hansang Lee a, Helen Hong b,n, Junmo Kim a a b

Department of Electrical Engineering, Korea Advanced Institute of Science and Technology, Daejeon, South Korea Department of Multimedia Engineering, Seoul Women's University, Seoul, South Korea

art ic l e i nf o

a b s t r a c t

Article history: Received 9 May 2014 Accepted 12 September 2014

We propose a graph-cut-based segmentation method for the anterior cruciate ligament (ACL) in knee MRI with a novel shape prior and label refinement. As the initial seeds for graph cuts, candidates for the ACL and the background are extracted from knee MRI roughly by means of adaptive thresholding with Gaussian mixture model fitting. The extracted ACL candidate is segmented iteratively by graph cuts with patientspecific shape constraints. Two shape constraints termed fence and neighbor costs are suggested such that the graph cuts prevent any leakage into adjacent regions with similar intensity. The segmented ACL label is refined by means of superpixel classification. Superpixel classification makes the segmented label propagate into missing inhomogeneous regions inside the ACL. In the experiments, the proposed method segmented the ACL with Dice similarity coefficient of 66.4777.97%, average surface distance of 2.24770.869, and root mean squared error of 3.53871.633, which increased the accuracy by 14.8%, 40.3%, and 37.6% from the Boykov model, respectively. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Knee MRI ACL segmentation Graph cuts Shape constraints Superpixelization

1. Introduction The anterior cruciate ligament (ACL) is a ligament located in human knee joint, which plays an important role in maintaining knee joint stability [1]. For diagnosis of ACL injury and its surgical guidance, the segmentation of the ACL in knee magnetic resonance (MR) images has received an increasing amount of attention. Numerous works have described the manual segmentation of the ACL as assisted by orthopedic experts for clinical purposes [2–4]. However, because manual segmentation is time-consuming and suffered by intra- and inter-observer variability, an automatic ACL segmentation technique is necessary. In knee MRI, an automatic segmentation of ACL is made challenging by several imaging characteristics of the ACL, including 1) adjacent soft tissues, e.g., the posterior cruciate ligament (PCL) and articular cartilage, which share similar intensity distribution with the ACL; and 2) intensity inhomogeneous regions inside the ACL, especially the region attached to the tibia. As shown in Fig. 1, during the ACL segmentation process, these adjacent soft tissues with similar intensity distribution cause the leakage of the segmented ACL label into these soft tissues, while the inhomogeneity inside the ACL results

n

Corresponding author. Tel.: þ 82 2 970 5756; fax: þ82 2 970 5981. E-mail address: [email protected] (H. Hong).

http://dx.doi.org/10.1016/j.compbiomed.2014.09.004 0010-4825/& 2014 Elsevier Ltd. All rights reserved.

in loss of the inhomogeneous region from the segmented ACL label when an intensity-based segmentation technique is applied. Due to these difficulties, only few works have been published with regard to the segmentation of the ACL. Ho et al. [5] suggested an ACL segmentation method based on morphological operations and the Chan–Vese active contours model [6]. Given that the Chan–Vese model is intensity-based, both leakage and loss of the inhomogeneous region inside the ACL can occur in the segmented ACL label such that experiments on this method demonstrate low performance; the specific results showed an overall Dice similarity coefficient (DSC) of 38.1 79.1% for four knee MRI sets. Lee et al. [7] suggested an ACL segmentation method based on adaptive thresholding and graph cuts [8]. Because graph cuts, which are originally based on intensity, are also associated with leakage problems, patient-specific shape constraints for graph cuts have been proposed to avoid leakage. As a result, most of the leakage was reduced with the help of shape constraints; however, the missing inhomogeneous region inside the ACL was not fully restored. In this paper, we propose a graph-cut-based segmentation method for the anterior cruciate ligament (ACL) in knee MRI with a novel shape prior and label refinement. To avoid leakage of the segmented ACL label into regions with similar intensity, two cost terms for graph cuts, referred to as fence and neighbor shape constraints, are proposed. These shape constraints depict the shape of the ACL from the graph cut seed and not from a previously segmented dataset such that they reflect patient-specific information onto the graph cuts.

2

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

To restore inhomogeneous regions inside the ACL which are not fully segmented via graph cuts, refinement by classifying over-segmented superpixels adjacent to the ACL label with location and intensity priors is proposed. As a result, the proposed method enables the segmented ACL label to avoid leakage into regions of similar intensity and to recover missing inhomogeneous regions inside the ACL. The remainder of the paper is organized as follows. In Section 2, we present the proposed method in detail. In Section 3, experimental results and their analyses are presented, while Section 4 discusses the contribution and effect of the proposed method. In Section 5, we conclude with a brief summary of the work.

background candidates are extracted by adaptive thresholding with Gaussian mixture model fitting and outlier removal on the coronal and sagittal planes to replace the interactive seeding in the Boykov model. Second, in the graph cuts step, iterated graph cuts on the sagittal plane are designed based on the pixel-wise graph nodes and edges with the two proposed cost terms of the “fence” shape cost and the “neighbor” shape cost, which respectively reflect the planar and spatial shape information of the ACL to the graph cuts. Finally, in the label refinement step, the segmented label is refined on the sagittal plane to restore the tibia-attached part of the ACL by means of superpixel classification.

2. Proposed method

2.1. Initial seed extraction

The proposed method consists of three main steps: initial seed extraction, graph cuts with patient-specific shape constraints, and label refinement. An outline of the proposed method is illustrated in Fig. 2. First, in the initial seed extraction step, ACL and

Replacing manual seeding for graph cuts in the Boykov model [8], the initial seed extraction step extracts object and background seeds in a semi-automatic manner. Because these extracted seeds will be used as the initial cues for the graph cuts, our objective in

Fig. 1. Anterior cruciate ligament in knee MRI. (a) ACL (arrowed) in left knee MRI viewed in coronal plane, (b) the case of ACL adjacent to PCL in coronal plane, (c) ACL (arrowed) viewed in sagittal plane, and (d) the case of inhomogeneity between femur-attached ACL (top-arrowed) and tibia-attached ACL (bottom-arrowed) in sagittal plane.

Fig. 2. The pipeline of anterior cruciate ligament segmentation algorithm.

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

this step is to make the object and background seeds represent candidates for the ACL and the background, respectively. The initial seed extraction step consists of two steps: adaptive thresholding with Gaussian mixture model fitting and outlier removal along the coronal and sagittal planes. In the first step, because the ACL usually consists of dark signal components in knee MRI, intensity thresholding is initially performed to extract preliminary candidates for the ACL and the background. To determine the threshold values adaptively for the given MRI, we estimate optimal threshold values by means of Gaussian mixture model (GMM) fitting [9]. As shown in Fig. 3(a) and (d), the signal distribution of knee MRI can be separated into two regions, a dark region which mainly represents soft tissues such as the ACL, PCL, meniscus, and cartilage, and a bright region which mainly represents knee bones and muscle. Thus, we perform GMM fitting on each knee MRI slice on the coronal plane with the number of Gaussians set to 2. As a result, two Gaussians NðμL1 ; σ L1 Þ and NðμH1 ; σ H1 Þ for dark and bright regions are extracted, where μL1 and μH1 , and σ L1 and σ H1 are the mean values for the dark and bright regions and the standard deviation values for the dark and bright regions, respectively. With these Gaussian parameters, two thresholds T O1 and T ℬ for the ACL and background candidates are defined as follows: T O1 ¼ μL1 ; T ℬ ¼ μH1 (

ð1Þ

p A O1

where

I p o T O1

pAℬ

where

Ip 4 T ℬ

ð2Þ

3

a set of connected components. With these clusters, because the ACL is most commonly located at the middle part of the VOI, the clusters for which most of their pixels are located at the top or the bottom parts of the VOI are rejected. The outlier removal is then performed again on the sagittal plane to reject the PCL and parts of the cartilage which are not fully rejected on the coronal plane from the ACL candidate, as shown in Fig. 4(d). Throughout the initial seed extraction process, we have two binary labels for the ACL candidate and the background candidates for each MRI slice. We use these ACL and background candidates as the initial object and background seeds, respectively, for the following graph-cut part. 2.2. Graph cuts with patient-specific shape constraints With the ACL and background candidates extracted during the initial seed extraction step as the object and background seeds, graph cuts are performed to segment the ACL from the knee MR images. However, the Boykov model, which is the original graphcut method proposed by Boykov et al. [8], is limited when used to cut the leakage to the adjacent soft tissues due to its intensitybased cost function. To cut the leakage clearly, we propose additive cost terms for graph cuts which reflect the patient-specific shape information of the ACL onto the graph cuts. For the entire set of labels assigned to pixels A ¼ fA1 ; :::; Ap ; :::; AjPj g, the proposed total cost function EðAÞ is stated as follows: EðAÞ ¼ ω1 U RðAÞ þ ω2 U BðAÞ þ ω3 U S2 ðAÞ þ ω4 US3 ðAÞ

ð5Þ

where I p is the intensity value of an arbitrary image pixel p, and O1 and ℬ are the sets of the ACL and background candidates, respectively. With these threshold values, we perform thresholding on the image to extract the ACL and background candidates. Fig. 3 (d) shows the first GMM fitting result, and Fig. 3(b) shows an image with the extracted ACL and background candidates. These are shown in red and blue, respectively. As shown in Fig. 3(b), the dark region is extracted so thickly that the actual ACL candidate is strongly connected to the adjacent soft tissues and bone boundaries. We refine this thresholded dark region to make the ACL candidate separable in the following steps by applying adaptive thresholding with Gaussian mixture model fitting again for the extracted dark region O1 . With the first extracted dark region O1 , for which the histogram is shown, in Fig. 3(e) we perform the second GMM fitting on this dark region O1 with the number of Gaussians set to 2 again. As a result, as shown in Fig. 3(e) two new Gaussians NðμL2 ; σ L2 Þ and NðμH2 ; σ H2 Þ for dark and relatively bright regions are separated where μL2 and μH2 , and σ L2 and σ H2 are the mean values and the standard deviation values for these regions, respectively. With these parameters, as shown in Fig. 3(c) a new threshold value T O for the ACL candidate and the new ACL candidate O thresholded by T O is defined as

here, RðAÞ; BðAÞ are the region and boundary cost terms, which are equivalent to those in the Boykov model; S2 ðAÞ; S3 ðAÞ are the proposed cost terms, referred to here as the “fence” shape cost and the “neighbor” shape cost, respectively; and ωi ; i ¼ 1; 2; 3; 4 are weighting parameters for the corresponding cost terms. In the experiments, the weights are equally set, as ω1 ¼ ω2 ¼ ω3 ¼ ω4 ¼ 1. The region cost term RðAÞ in Eq. (5) is defined as

T O ¼ μL2 þ λ U σ L2

ð3Þ

δðAp ; Aq Þ ¼

p A O where I p oT O

ð4Þ

! ðI p  I q Þ2 Bfp;qg ¼ exp  ; 2σ 2

where λ A f1; 2; 3g is the weighting coefficient for the standard deviation of the new object threshold; it is experimentally determined for each instance of 3D MRI data. With this ACL candidate, outlier removal, which rejects nonACL regions such as PCL and cartilage, is performed on both coronal and sagittal planes to create a reliable ACL seed for the graph cuts. To cut or weaken the connectivity between ACL and adjacent soft tissues and bone boundaries, the morphological operations of erosion and dilation are performed on the ACL candidate. We then cluster the pixels of the ACL candidates into

RðAÞ ¼ ∑ Rp ðAp Þ

ð6Þ

Rp ð}obj}Þ ¼ jI p  I O j; Rp ð}bkg}Þ ¼ jI p  I ℬ j;

ð7Þ

pAP

IO ¼

1 1 ∑ Ip ; Iℬ ¼ ∑ Ip jOjp A O jℬjp A ℬ

ð8Þ

where I p is the intensity value of pixel p; I O and I ℬ are the mean intensity values for the object and background seeds, respectively, and O and ℬ are the sets of pixels in the object and background seeds, respectively. The boundary cost term BðAÞ is defined in a manner equivalent to that of the Boykov model, i.e. BðAÞ ¼

∑

fp;qg A N

where



Bfp;qg U δðAp ; Aq Þ

1

if Ap a Aq

0

otherwise

;

ð9Þ

ð10Þ

ð11Þ

and N is the entire set of neighboring edges in the graph. The “fence” shape cost term S2 ðAÞ reflects the planar shape information of the ACL based on the object ACL candidate seed and the background candidate seed for a single MRI slice. Based on the assumption that the object ACL candidate seed is extracted with a proper ACL shape for a single MRI slice, one objective of this region-based cost term is to prevent the segmented label from leakage into the PCL or cartilage by maintaining the general shape

4

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

of the object ACL candidate seed. To achieve this goal, we formulate the “fence” shape cost based on a distance map, where the sources of the distance map are from the extracted object and background seeds instead of being an atlas or pre-segmented training data. For the set of entire segmented labels A, the fence shape cost S2 ðAÞ is defined as S2 ðAÞ ¼ ∑ S2 ðAp Þ

ð12Þ

S2 ð}obj}Þ ¼ α1 minO A O distðp; OÞ  α2 minB A ℬ distðp; BÞ þ C 2

ð13Þ

S2 ð}bkg}Þ ¼ 0

ð14Þ

pAP

Fig. 3. Adaptive thresholding using GMM fitting; (a) an original image in the VOI, (b) initially thresholded labels for ACL (red) and background candidates (blue) on the original image, (c) finally thresholded labels for ACL (red) and background candidates (blue) on the original image, (d) histogram of the original image and fitted Gaussian distributions, and (e) histogram of the extracted dark region (top) and its fitted Gaussian distributions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

where distðp; qÞ is the Euclidean distance between ðp; qÞ; α1 and α2 are weighting parameters for distance costs; and C 2 is shifting constant to make the cost S2 ðAp Þ nonnegative. From Eq. (13), for an arbitrary pixel p, the fence shape cost is defined as the difference between the closest distance of the object seed from p and the closest distance of the background seed from p, as shown in Fig. 5(c). If a pixel p is much farther from the object seed than it is from the background seed, then the S2 ð}obj}Þ cost is high such that the labeling of p as an object is prevented. On the other hand, if a pixel p is much closer to the object seed than to the background seed, then the cost of S2 ð}obj}Þ is low such that labeling p as an object is encouraged, as shown in Fig. 5(d). The “neighbor” shape cost S3 ðAÞ reflects the spatial shape information of the ACL to the graph cuts based on all of the initial object and background seeds for the 3D knee MRI data. Based on the assumption that the ACL evolves smoothly along the z-axis and the object ACL candidate seed represents this evolution well, one goal of this additionally region-based cost term is to enhance the segmentation by representing the shape of the ACL in the neighboring MR slices, where the ACL is not properly segmented into specific MR slices due to the insufficient information in the initially extracted seed. To achieve this goal, we formulate the “neighbor” shape cost based on a probability map by accumulating all of the slices of each ACL and background seed. For an arbitrary pixel p of an arbitrary image slice I, we first define the label

Fig. 4. Outlier removal on coronal (top) and sagittal (bottom) planes: (a) original images in the VOI, (b) initial thresholded binary masks for ACL candidate (red) on the original image, (c) ACL candidate (red) refined via outlier removal in coronal plane, and (d) final ACL candidate (red) refined via outlier removal in sagittal plane. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

5

Fig. 5. Shape constraints for graph cuts: (a) an original image in the VOI scale, (b) object (red) and background (blue) seeds on the original image, (c) the computation of “fence” shape cost corresponding to pixel p, (d) fence shape costs represented as a distance map, and. (e) “neighbor” shape costs represented as a probability map. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

frequencies of p to be in ACL seed O or background seed ℬ via Prðp A OÞ ¼

1 ∑ I p A OjI m IAM

ð15Þ

Prðp A ℬÞ ¼

1 ∑ I p A ℬjI m IAM

ð16Þ

where M ¼ fIg is the entire set of 3D MRI data, represented as a set of m sagittal MRI slices; and I is an indicator function of pixel labels for a given MRI slice, defined as follows:  I p A OjI ¼

1

if p A O for given I

0

otherwise

:

ð17Þ

With these probability-map-based label frequencies, the fence shape cost S3 ðAÞ is defined as S3 ðAÞ ¼ ∑ S3 ðAp Þ

ð18Þ

S3 ð}obj}Þ ¼ β 1 Prðp A ℬÞ  β 2 Prðp A OÞ þ C 3

ð19Þ

S3 ð}bkg}Þ ¼ 0

ð20Þ

pAP

where β 1 and β2 are weighting parameters for probability costs, and C 3 is a constant to make the cost S3 ðAp Þ nonnegative. From Eq. (19), for an arbitrary pixel p, the neighbor shape cost is defined by the difference between the probability of p to be in the background seed and the probability of p to be in the object seed. For a pixel p, if the probability of p to be in the object seed is much higher than that of p to be in the background seed, then the cost of S3 ð}obj}Þ is low such that labeling p as an object is encouraged, as shown in Fig. 5(e). With the constructed graph, we perform the max-flow min-cut computation [10] to obtain the min-cut of the graph. To stabilize the segmented label, an iterated graph-cut scheme is incorporated by updating the object and background seeds conservatively with the segmented results of each step. The termination criteria for the iterative graph cuts include the volume stability of the segmented ACL label and the maximum number of iterations, i.e., 10 in our experiment. In all of the patient-specific shape constraints, the distancemap-like shape cost allows the segmentation to avoid leakage into the adjacent regions with similar intensity, and the probabilitymap-like shape cost allows the segmentation to overcome the lack of information due to the insufficient seeds in specific MRI slices.

2.3. Label refinement with superpixel classification In the ACL label segmented by the graph cuts, the true ACL is shown to be under-segmented in the bottom part of the ACL, in which the ACL is presented as particularly inhomogeneous, especially in the region attached to the tibia. To recover this missing region, a method that refines the ACL label by means of label extension based on superpixel classification is proposed. The label refinement step consists of two steps: the superpixelization of sagittal knee MR images and the classification of superpixels according to intensity and location priors. Fig. 6 shows an example image and its results corresponding to each step in the label refinement process. A superpixel [11] is defined as a cluster of pixels of an image which usually represents an over-segmented region. Because the objectives of superpixelization are to reduce the degree of image complexity and to avoid under-segmentation via the smoothing of images, the use of superpixels in the graph cuts tends to make the segmentation result robust to image inhomogeneity. Among a number of superpixelization works [12,13], we use “TurboPixels” [12] to create superpixels for the bottom half region of the sagittal knee MRI, which mainly includes the tibia-attached ACL. Formulating the oversegmentation task as a set of locally interacting region growing problems, TurboPixels perform rapid over-segmentation by means of level-set-based curve evolution. Fig. 6(b) shows a result of the superpixelization of sagittal knee MR images. In the experiments, the number of superpixels for each slice was set to 100. Classification of the superpixels is then performed considering the intensity and location conditions. Because the objective of superpixel classification is to label the superpixels which belong to the part of the ACL attached to the tibia, we construct an iterative deterministic classification scheme with two prior conditions: (1) the missing ACL label is located between the tibia and the segmented ACL label specifically on the left-bottom side, and (2) the mean intensity of the missing label is brighter than that of the segmented ACL labels but darker than that of the background labels. With these assumptions, we classify the superpixels, which are adjacent to the segmented ACL label, into the ACL and background in an iterative sense. Fig. 6(c) shows a result from the first iteration of superpixel classification, and Fig. 6(d) shows the finally refined ACL label. In the experiments, the terminal condition of iteration was determined as the stability of the label. With the superpixelization and the classification scheme according to the location and mean intensity conditions, superpixels corresponding to missing inhomogeneous regions of the ACL are correctly detected without leakage into other soft tissues or bones. As a result, the final ACL label restores the inhomogeneous region at

6

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

Fig. 6. Label refinement via superpixel classification: (a) an original image with the ACL label (red) segmented by the graph cuts, (b) over-segmented superpixels at the bottom-half of knee MRI, (c) superpixelized ACL label (red) and ACLclassified superpixels (yellow) added to the ACL label in 1st iteration, and (d) finally refined ACL label (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the part of the ACL attached to the tibia, which are not fully detected during the graph cut step.

3. Experimental results The proposed method was tested using proton-density (PD) weighted images of knee joints on the coronal plane acquired from four normal subjects using a Philips Achieva 3.0T system with a repetition time (TR) of 1800 ms and an echo time (TE) of 35 ms. Subjects are all males with the ages of 33, 25, 42, and 16 in which the first two scanned the left knee and the other two scanned the right knee. Each MRI scan included 250 slides of MR images with a matrix size of 512  512, a pixel size of 0.31 mm  0.31 mm, and a slice thickness of 1 mm. The performance of our method was evaluated by means of a visual assessment and with respect to accuracy. Figs. 7 and 8 show the ACL segmentation results in sagittal and coronal planes, showing comparisons of the proposed method and three other methods. Figs. 7 and 8(a) show the sagittal knee MR images in the experimental data. Figs. 7 and 8(b)–(f) show the ACL labels segmented by the different segmentation methods: (b) graph cuts with only the intensity-based cost terms in the Boykov model (IB); (c) the typical shape-based graph-cut model in which the shape costs are generated in a probability map with the training-ground truth labels (TBSC); (d) the proposed method via graph cuts with patientspecific shape constraints (PSSC), (e) the proposed method via graph cuts with PSSC and label refinement, and (f) the ground truth ACL labels which were manually segmented by an orthopedist. The initial seeds for the ACL and background results in the three compared methods are identical to those of the proposed method. In the TBSC, for each MRI data, the ground truth labels for the training set are selected from other MRI data in a leave-one-out manner. In Fig. 7(b), leakage into regions of similar intensity, including the PCL, occurs in the segmented ACL label such that the resulting label includes inaccurate, dark regions. In Fig. 7(c), the leakage into the PCL shown in (b) is greatly reduced with the help of shape constraints. In Fig. 7(d), the segmentation accuracy is improved at the regions attached to the femur due to the patientspecific shape information, which reflects the local information associated with the patient such that the ACL at the femur boundaries is effectively detected compared to the TBSC results. In the part of the ACL in the region attached to the tibia, however, due to the intensity inhomogeneity represented in the ACL, as shown in

the fifth through seventh rows in Fig. 7(a) and given the lack of information in the initial seed resulting from the inhomogeneity, the ACL in the region attached to the tibia is still under-segmented, as shown in the fifth through seventh rows in Fig. 7(d). In Fig. 7(e), the resulting ACL labels show that the proposed method accurately segments the ACL in a manner similar to the ground truth as compared to the other methods, as the superpixel classification technique mostly recovers the missing parts of the ACL in the region attached to the tibia. The comparative results on sagittal plane can also be found on coronal plane. In Fig. 8(b), the heavy leakages into similar intensity region occur in IB-segmented labels. In Fig. 8(c) and (d), most leakages are avoided with the help of shape constraints in TBSC and PSSC. However, since the shape constraints in these graph cuts act like restriction from extending the label into similar intensity regions, they have side effect of sensitiveness to intensity inhomogeneity. As shown in the arrowed regions in Fig. 8(d), this side effect results in heavy losses of labels in inhomogeneous region inside ACL especially at the tibia-attached region. By applying the superpixel classification to overcome this limitation, the inhomogeneous region lost in the PSSC results are mostly recovered so that the ACL labels segmented by the proposed method in Fig. 8(e) are represented the most similar to those of ground truth in Fig. 8(f). To validate the segmentation accuracy of the proposed method, we provided a quantitative analysis of the experimental results with two evaluation criteria: the Dice similarity coefficient (DSC), the average surface distance (ASD), and the root mean square distance error (RMSE). The ACL labels segmented manually by an orthopedist are depicted as the ground truth. DSC is a similarity measure which represents the overlapping ratio of two regions, the segmented region and the ground truth in this case: DSC ¼

2jX \ Yj jXjþ jYj

Here, X is the segmented label, Y is the ground truth label, and jU j refers to the area of the label. The DSC values are distributed as [0,1], indicating that the segmentation accuracy increases as DSC increases. ASD is defined as the average Euclidean distance between two surfaces including segmented region X and the ground truth region Y, as follows: ! 1 ASD ¼ ∑ dðsX ; SðYÞÞ þ ∑ dðsY ; SðXÞÞ jSðXÞj þ jSðYÞj sX A SðXÞ sY A SðYÞ where SðXÞ and SðYÞ indicates the sets of surface points sX ,sY in X and Y, respectively, and dðv; SÞ refers the Euclidean distance between the point v and the surface S, which is computed by dðv; SÞ ¼ minjjv  sjj2 : sAS

RMSE is defined as the root mean square distance between two surfaces including segmented region X and the ground truth region Y, as follows: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ∑ ðdðsX ; SðYÞÞÞ2 þ ∑ ðdðsY ; SðXÞÞÞ2 RMSE ¼ jSðXÞj þ jSðYÞj sX A SðXÞ sY A SðYÞ Table 1 shows the DSC, ASD, and RMSE computed on the experimental dataset for the proposed method and the three compared methods. The average DSC, the average ASD, and the average RMSE of the ACL segmented by the intensity-based graph cuts (IB) are estimated to be 57.85 74.45%, 3.7657 0.707, and 5.67071.486 respectively, due to heavy leakage despite the fact that the segmented labels include most of the ACL in the ground truth. When the graph cuts with the training-based shape constraints (TBSC) are applied, the average DSC is estimated to be 57.90 7 8.00%, which is similar to that of the IB, the average ASD is

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

7

Fig. 7. The results of ACL segmentation in sagittal plane by the proposed method and three comparative methods: (a) original knee MR images with ACL (arrow ) and PCL (arrow head), (b)–(f) ACL labels (red) are segmented by (b) IB, (c) TBSC, (d) PSSC, (e) proposed method, and (f) ground truth. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

estimated to be 3.508 71.082, showing a decrease of 6.8% compared to that of the IB, and the average RMSE is estimated to be 5.760 71.928, showing an increase of 1.6% compared to that of the IB. The label segmented by the IB has high sensitivity, as the label captures most of the actual ACL; however, it also contains a large false region, similar to leakage, such that it has low specificity and low overall accuracy. On the other hand, the label segmented by the TBSC has high specificity given that since the label rejects most of the leakage. However, it fails to continue to capture the actual ACL in the parts attached to the femur and tibia. Thus, the sensitivity decreases and consequently TBSC still shows low overall accuracy. When the graph cuts with the patient-specific shape constraints (PSSC) are applied, the leakage is mostly avoided, like the TBSC, and it captures part of the ACL which is lost in the TBSC results, such that the average DSC is estimated to be 61.57710.78%, showing an increase of 3.3% as compared to the TBSC and IB results.

In addition, the average ASD is estimated to be 2.67671.194, which is decreased by 23.7% compared to that of the TBSC and by 28.9% compared to that of the IB, and the average RMSE is estimated to be 4.24172.099, which is decreased by 26.4% compared to that of the TBSC and by 25.2% compared to that of the IB. It is important to note that the PSSC reflects the patient-localized information of the ACL shape to the graph cuts and in so doing enhances ACL segmentation relative to that by the IB and the TBSC methods. In the ACL segmented by the proposed method, the average DSC is estimated to be 66.4777.97%, which is increased by 5.3% compared to that of the PSSC, while the average ASD is estimated to be 2.24770.869, showing a decrease of 16.0% compared to that of the PSSC, and the average RMSE is estimated to be 3.53871.633, showing a decrease of 16.6% compared to that of the PSSC, thus also showing the finest performance compared to the other methods. It appears that superpixel classification captures the actual ACL in the part of the ACL

8

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

Fig. 8. The results of ACL segmentation in coronal plane by the proposed method and three comparative methods: (a) original knee MR images with ACL (arrow), (b)–(f) ACL labels (red) are segmented by (b) IB, (c) TBSC, (d) PSSC, (e) proposed method, and (f) ground truth. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

attached to the tibia, which is not properly segmented by the PSSC method, causing the overall segmentation accuracy to increase. To validate the effect of the superpixel classification of the proposed method, we observe the two sub-sections vertically dividing ACL; the stem (ST) and the tibia-attached (TA) regions, and measure the performance separately via DSC, ASD, and RMSE, as shown in Table 2. Since the superpixel classification is specifically targeted to the under-segmented region of ACL near the tibia, we redefine the volumes of interest for two regions, in which the stem region includes the central part of ACL and the tibia-attached region includes the part of ACL attached to tibia. As shown in the fifth through seventh rows in Fig. 7(d), inhomogeneous regions inside the tibia-attached part of ACL are heavily lost by the shapeconstrained graph cuts such as TBSC and PSSC. Since the shapeconstrained graph cuts have weakness to intensity inhomogeneity, they failed to capture the inhomogeneous part of ACL so that the segmentation accuracy of them in the tibia-attached region is even worse than that of IB, as shown in Table 2. Overcoming this limitation by the superpixel classification, the proposed method shows improved performance in terms of DSC, ASD, and RMSE compared to the same results from both the shape-constrained

graph cuts (TBSC and PSSC). In the stem sub-section, the average DSC of the proposed method is estimated to be 75.1176.25%, which is increased by 5.1% compared to that of the PSSC and increased by 7.5% compared to that of TBSC, while the average ASD is estimated to be 1.55670.352, showing a decrease of 14.3% compared to that of the PSSC and a decrease of 32.8% compared to that of the TBSC, and the average RMSE is estimated to be 2.31970.549, showing a decrease of 18.0% compared to that of the PSSC and a decrease of 34.5% compared to that of the TBSC. In the tibia-attached subsection, the average DSC of the proposed method is estimated to be 48.14731.12%, which is increased by 10.2% compared to that of the PSSC and by 21.0% compared to that of the TBSC, while the average ASD is estimated to be 4.86375.318, showing a decrease of 22.6% compared to that of the PSSC and of 47.0% compared to that of the TBSC, and the average RMSE is estimated to be 6.13775.933, showing a decrease of 20.0% compared to that of the PSSC and of 44.5% compared to that of the TBSC. These results demonstrate that the superpixel classification of the proposed method preserves the overall performance in the stem part and restores the missing inhomogeneous regions of the bottom ACL in the part attached to the tibia.

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

Table 1 Comparison of different segmentation methods in terms of Dice similarity coefficient (DSC), average surface distance (ASD), and root mean square distance (RMSE).

9

Table 2 Comparison of different segmentation methods for two sub-regions i.e. stem (ST) and tibia-attached (TA) regions, in terms of Dice similarity coefficient, average surface distance, and root mean square distance (RMSE).

DSC (%)

Patients

DSC (%)

Methods

1

2

3

4

Total

Regions

IB TBSC PSSC Proposed

58.23 66.34 73.34 74.14

63.58 62.93 67.54 70.7

56.8 49.58 49.76 65.17

52.81 52.76 55.63 55.87

57.86 7 4.45 57.90 7 8.00 61.577 10.78 66.477 7.97

ST

ASD (pixels)

Patients

Methods

1

2

3

4

Total

IB TBSC PSSC Proposed

4.553 2.353 1.509 1.522

2.834 3.132 1.864 1.792

3.846 3.617 3.284 2.192

3.826 4.930 4.048 3.483

3.7657 0.707 3.508 7 1.082 2.6767 1.194 2.2477 0.869

RMSE (pixels)

Patients

Methods

1

2

3

4

Total

IB TBSC PSSC Proposed

7.729 3.697 2.255 2.269

4.180 5.692 2.837 2.803

5.348 5.307 5.053 3.153

5.423 8.345 6.820 5.926

5.6707 1.486 5.7607 1.928 4.2417 2.099 3.538 7 1.633

TA

ASD (pixels)

Methods

Patients 1

2

3

4

Total

IB TBSC PSSC Proposed IB TBSC PSSC Proposed

64.51 70.56 74.73 74.88 57.52 61.33 73.96 75.84

80.13 81.32 81.28 93.25 61.95 37.12 57.31 63.53

58.87 56.86 55.81 74.27 66.4 8.89 14.58 48.64

66.53 61.58 68.09 68.04 37.74 1.02 5.94 4.53

67.517 9.02 67.58 7 10.78 69.98 710.87 75.117 6.25 55.90 712.64 27.09 727.59 37.95 7 32.88 48.147 31.12

Methods

Patients

Regions

1

2

3

IB TBSC PSSC Proposed IB TBSC PSSC Proposed

4.084 2.204 1.530 1.569 4.241 2.880 1.481 1.475

1.581 1.454 1.157 1.066 2.516 6.432 2.187 2.087

4.209 2.691 2.213 1.670 2.219 8.582 7.225 3.113

RMSE (pixels) Regions

Methods

Patients 1

2

3

ST

IB TBSC PSSC Proposed IB TBSC PSSC Proposed

6.585 3.486 2.184 2.236 7.526 4.432 2.481 2.424

2.420 2.364 1.793 1.682 3.510 9.397 3.009 2.994

5.914 3.939 3.370 2.339 3.043 9.980 8.923 4.160

ST

TA

4 2.892 2.920 2.365 1.892 6.167 18.79 14.252 12.775

4

Total 3.1917 1.227 2.3177 0.648 1.816 70.570 1.556 7 0.352 3.786 71.821 9.171 76.830 6.286 75.894 4.863 75.318

Total

4. Discussion In this paper, we proposed a graph cut-based segmentation method for the ACL in knee MRI with patient-specific shape constraints and label refinement. ACL segmentation has two major challenges, in this case the presence of adjacent similar intensity regions, e.g. PCL and cartilage in which the leakage occurs when intensity-based segmentation is applied, and missing labels in intensity inhomogeneous regions inside the ACL, especially occurring in the regions attached to the tibia. Considering the leakage problem, because the intensity-based graph cuts are significantly weak to leakage into similar intensity regions, shape constraints intended to avoid leakage are required. Conventional shape constraints for graph cuts are mainly modeled by training-based shape constraints (TBSC). However, the TBSC method is limited in terms of local shape descriptions, such as failing to capture bone-attached ACL areas during ACL segmentation. To overcome this limitation associated with conventional shape constraints, we proposed patient-specific shape constraints (PSSC), which extract the shape information of the object from the patient image itself and reflect it onto the graph cuts. Two PSSCs, distance-map-based fence-shape constraints and probabilitymap-based neighbor shape constraints, are computed from the initially extracted object and background seeds to respectiv ely reflect the planar and spatial shape information onto the graph cuts. For the missing label problem in inhomogeneous regions inside the ACL, as TBSC and PSSC shape constraints restrict the label extension as a “fence,” they are limited by missing labels in inhomogeneous regions. To recover the missing labels in inhomogeneous regions inside the ACL, especially at the tibia-attached part of the ACL, we perform label refinement based on superpixel classification. Superpixels, clusters of pixels produced by edgepreserving over-segmentation, are used in the process of label refinement to handle the inhomogeneity by smoothing the intensities among the pixels. The number of superpixels was one parameter to be considered in the experiments. If the number of superpixels is too low, edge-preservation of superpixelization weakens such that the outlier region outside the boundary of

TA

4.395 4.377 3.960 3.020 9.044 20.383 16.27 14.970

4.828 71.849 3.541 70.865 2.8277 1.010 2.319 70.549 5.781 72.963 11.048 7 6.703 7.6717 6.434 6.1377 5.933

the ACL can be captured during the label refinement process. If the number of superpixels is too high, the smoothing of the pixel intensities among the superpixels weakens such that the handling of inhomogeneities is not effective during the label refinement process. Fig. 9 shows the performances of superpixel classification for the numbers of superpixels from 5 to 250, measured by DSC, ASD, and RMSE. Performances are represented mostly stable among the numbers of superpixels, while they tend to deteriorate where the number is less than 50 or larger than 200. From this observation, the number of superpixels in our experiments was experimentally set to 100. In superpixel classification, we essentially assume two things: 1) the missing labels at the tibia-attached ACL are located on the bottom-left side of the segmented ACL label, and 2) the intensity distribution of the missing labels is usually brighter than that of the segmented ACL and darker than that of the background seed region, e.g., femur and tibia bones. Based on these two assumptions, we perform the superpixel classification under certain location and intensity conditions of superpixels in an iterative manner. In the experiments, the proposed method was tested on four clinical datasets and was compared with the intensity-based Boykov model (IB) and the conventional TBSC method. In the graph cuts, PSSC not only cut the leakage compared to the IB method, but it also captured the parts of the bone-attached ACL which TBSC mainly missed by reflecting the patient-specific local shape information. Furthermore, through label refinement, the proposed method also recovered most of the tibia-attached ACL compared to PSSC, which failed to capture this part due to the inner-inhomogeneity issue. As a result, the ACL label segmented by the proposed method extracted the

10

H. Lee et al. / Computers in Biology and Medicine 55 (2014) 1–10

most accurate region of the ACL from both qualitative and quantitative perspectives. The processing times of comparative methods were 0.308, 0.308, 0.309, and 0.367 s on average per each slice of MRI for IB, TBSC, PSSC, and the proposed method, respectively. While IB, TBSC, and PSSC took nearly the same processing time each other, the proposed method took an extra processing time of superpixelization, measured as 0.058 s on average. However, the extra processing time was not significant. All of our evaluations were performed on MATLAB with 2.4GHz Intel Core i7-4700MQ CPU and 12.0GB of RAM. The segmentation of ACL can be applied for various clinical applications. For instance, during the surgical tissue graft replacement in an ACL reconstruction surgery, since clinicians have very limited knowledge of the original bone-ACL attachment sites, they usually located the graft experimentally. With the proposed method, more natural sites for graft replacement can be found by detecting the bone-ACL attachment sites of other side of the knee which is not damaged, with the ACL and knee bone segmentation, and symmetrically registering these sites to the damaged knee bones. Furthermore, it could also be applied for ACL injury detection by observing signal inhomogeneity inside the segmented ACL, since the torn ACL region is represented roughly inhomogeneous in knee MR images. Besides the extensive data acquisition with multiple resolutions, and their experiments and analyses, future works will focus on applying our method to various applications including graft site detection and ACL injury detection.

5. Conclusions In this paper, we proposed the graph cut-based segmentation method for the ACL in knee MRI with patient-specific shape constraints and superpixel-classification-based label refinement. In the graph cuts, to avoid the label leakage into regions with similar intensity, two additional shape constraints were formulated to reflect the planar and spatial shape information of the ACL from the extracted seeds to the graph cuts. To restore the missing inhomogeneous regions of the ACL, refinement of the segmented labels was performed by classifying superpixels adjacent to ACL labels with location and intensity priors. In the experiments, the proposed method segmented the ACL with an average DSC of 66.47 77.97%, which increased the accuracy by 15%, 15%, and 9% from the Boykov model, the TBSC, and the PSSC graph cuts, respectively; the average ASD was 2.247 70.869, which decreased the error by 40%, 36%, and 16% from the Boykov model, TBSC, and PSSC graph cuts, respectively. Conflict of interest statement We declare that we have no conflict of interest. Acknowledgments This research was supported in part by Seoul R&BD Program (ST110036) and ETRI R&D Program (14ZC1400, The Development of a Realistic Surgery Rehearsal System based on Patient-Specific Surgical Planning) funded by the Government of Korea. Authors are grateful to Dr. Joon Ho Wang at the Samsung Medical Center in Seoul, Korea for providing the knee MR data sets and their manual segmentation results. References [1] C.L. Arden, K.E. Webster, N.F. Taylor, J.A. Feller, Return to sport following anterior cruciate ligament reconstruction surgery: a systematic review and meta-analysis of the state of play, Br. J. Sports Med. 45 (2011) 596–606. [2] J.C. Gardiner, J.A. Weiss, Subject-specific finite element analysis of the human medial collateral ligament during Valgus knee loading, J. Orthop. Res. 21 (2003) 1098–1106. [3] G. Limbert, J. Middleton, M. Taylor, Finite element analysis of the human ACL subjected to passive anterior tibial loads, Comput. Methods Biomech. Biomed. Eng. 7 (1) (2004) 1–8. [4] A.M.W. Chaudhari, E.A. Zelman, D.C. Flanigan, C.C. Kaeding, H.N. Nagaraja, Anterior cruciate ligament-injured subjects have smaller anterior cruciate ligaments than matched controls: a magnetic resonance imaging study, Am. J. Sports Med. 37 (7) (2009) 1282–1287. [5] J.H. Ho, W.Z. Lung, C.L. Seah, C.L. Poh, K. Sheah, D.T.T. Lie, K.S.A. Yew, Anterior cruciate ligament segmentation: using morphological operations with active contour, in: Proceedings of the 4th International Conference on Bioinformatics and Biomedical Engineering (iCBBE), 2010. [6] T.F. Chan, L.A. Vese, Active contours without edges, IEEE Trans. Image Process. 10 (2) (2001) 266–277. [7] H. Lee, H. Hong, J. Kim, Anterior cruciate ligament segmentation from knee MR images using graph cuts with geometric and probabilistic shape constraints, in: Proceedings of the Asian Conference on Computer Vision (ACCV), 2012. [8] Y. Boykov, G. Funka-Lea, Graph cuts and efficient N-D image segmentation, Int. J. Comput. Vis. 70 (2) (2006) 109–131. [9] Z.K. Huang, K.W. Chau, A new image thresholding method based on Gaussian mixture model, Appl. Math. Comput. 205 (2008) 899–907. [10] V. Kolmogorov, R. Zabin, What energy functions can be minimized via graph cuts? IEEE Trans. Pattern Anal. Mach. Intell. 26 (2) (2004) 147–159. [11] X. Ren, J. Malik, Learning a classification model for segmentation, in: Proceedings of the IEEE International Conference on Computer Vision (ICCV), 2003, pp. 10–17. [12] A. Levinshtein, A. Stere, K.N. Kutulako, D.J. Fleet, S.J. Dickinson, K. Siddiqi, TurboPixels: fast superpixels using geometric flows, IEEE Trans. Pattern Anal. Mach. Intell. 31 (12) (2009) 2290–2297. [13] S. Yu, J. Shi, Multiclass spectral clustering, in: Proceedings of the IEEE International Conference on Computer Vision, vol. 1, 2003, pp. 313–319.

Fig. 9. Performance of superpixel classification for numbers of superpixels measured by (a) DSC, (b) ASD, and (c) RMSE.