Automatic Articular Cartilage Segmentation Based on Pattern Recognition from Knee MRI Images

Edge Detection Based on Iterative Threshold

The Canny edge detector [8] is an ideal detection operator because of its sensitivity and high signal-to-noise ratio (SNR). An improved self-adaptive threshold Canny detector was used to set the initial gradient threshold (T_l, T_h; T_l < T_h) introduced by Huo [9], which combines the advantages of the Otsu method and probability model method to obtain the low threshold T_l and high threshold T_h adaptively. T_l and T_h are the key parameters for detecting and connecting edges using the Canny detector. The quality of the results obtained using the Canny detector is closely related to the threshold value; a low threshold will generate many weak edges, and a high threshold may not be able to detect the BCI. Thus, the threshold value must be set adaptively based on the premise that BCI can be detected completely, and it is important to set a suitable gradient threshold adaptively, which can greatly reduce the number of iterations required. The edges of skin, bone, and cartilage were strong in the MRI images; therefore, other weak edges could be removed by setting a higher threshold. The Canny detector was used to extract the edges of the two MRI images, as shown in Fig. 1, and the results are shown in Fig. 3.

The main edges in Fig. 3 are the skin, femur, tibia, patella, and muscle tissue, which have a high gradient. The BCI was detected completely, which contained many pixels and was larger than the majority of other edges. However, the BCI was not always detected completely in all slices. In this case, the cartilage could not be segmented successfully when the edge detection was poor. Thus, further improvement of the threshold setting was required. The BCI was identified as a “strong edge” (many pixels, large area, and high gradient) by edge detection compared to other edges, which could be determined by further detection experiments. Therefore, the threshold could be corrected with the feedback obtained from calculating the number of “strong edges”. This process is illustrated in Fig. 4. With the initial threshold, the number of strong edges was calculated first. When the number was excessively high (e.g., more than 10), the threshold had to be increased and lowered in small increments so that a better edge detection result could be obtained using this adaptive threshold setting. The number of iterations was typically less than 50.

Edge Optimization

Other edges were also detected using the adaptive threshold (Fig. 3), e.g., example, muscle tissue, compact bone, and cancellous bone, which resulted in some edges that are broken and coarse, with many “glitches” and short branches. Therefore, we used morphological processing (connected component labels, spur removal, edge connection, and thinning) to improve the results. After the morphological processing of the results shown on the right side of Fig. 3, the optimized results shown in Fig. 5 were obtained.

Cartilage covers the bearing surface of the knee, and the curvature changes slowly; however, a large part of the BCI extracted by edge detection is not covered by cartilage shown in Fig. 5. Thus, corner detection was used to break the edges [10], which allowed for more accurate BCI detection. The corner detection results are shown in Fig. 6, except for the edge endpoints. There still are two corners that divide the patella edge into three parts. The cartilage only covers the right one, and a similar configuration can be gathered from the edges of the femur and patella. Corners are always detected in areas where the curvature changes fast, where there is no cartilage coverage; thus, the search space of cartilage is narrowed, and the identification error caused by unrepresentative interconnecting edges is reduced (shown in Fig. 6).

SVM-Based Edge Identification

The optimized edges include skin, BCI, and other strong edges, and the BCI must be identified first before calculating the distance from the BCI, which separates different cartilage from other edges. If optimized edges are not identified accurately, the cartilage would not be segmented.

The edges are classified into six classes (left and right skin, tibial edge, femoral edge, patella edge cartilage, and the other edges). The skin edge has many pixels and a high length-to-width ratio, and each of its sides (left or right) does not have foreground pixels. Similarly, five main features (shown in Table 1) are used to classify the BCI and other edges.

As shown in Table 1, the BCI can be distinguished from other edges by certain features (opening direction, location, number of pixels in the horizontal direction, etc.), which are discussed in detail in a later section. Although some of the short branches are removed by morphological processing, there are still few edges with more than two endpoints; therefore, a contour tracing method involving the eight neighboring pixels is used to obtain the edge sequence and break the multiple-branched edge into a single sequence based on which all features are calculated.

As shown in Fig. 7, the opening direction of an edge sequence is calculated by determining where the intermediate points are located along the line connecting endpoints, which determines whether the opening direction is downward, upward, or straight, and three values are sufficient to recognize the BCI. When the intermediate points are a located at a short distance (e.g., 10) from the connecting line of the endpoints, the edge is determined to be straight. The intermediate points are equally selected along the edge. The final value of the opening direction is the ratio of the number of upward, downward, and straight points to the total number of points (20, representing a balance between accuracy and computational efficiency). The location of an edge is the centroid coordinate of a finite set of points in the vertical and horizontal directions [11].

The eight-direction Freeman chain code is calculated first when the number of pixels in the horizontal direction and along a clockwise rotation angle are calculated, as shown in Fig. 8. The horizontal direction is 0 or 4 for the chain code. Thus, the number of pixels in the horizontal direction is the ratio of the chain code number, 0 or 4, to the total number of pixels in a sequence. The clockwise rotation angle is calculated as follows: In the first code direction, calculate the clockwise rotation angle of the chain code, then determine and compare the rotation angle value with that of the next chain code, which is calculated individually. When the angle is greater than 180°, it is considered negative, and the cumulative value of the angle is determined to be the final clockwise rotation angle.

The edges can then be classified via the SVM using the feature vector of each edge. However, the degree of difficulty with which different edges are recognized varies, and the order from high to low degree of difficulty is skin edge, femoral edge, tibial edge, and patella edge. According to this order, the one-versus-one SVM method [12] is used to identify edges. This method trains a classifier after every two classes. Thus, for six classes, the method will train 6 × 5 / 2 = 15 classifiers.

The kernel of the SVM model is a linear function. The edge features previously introduced for classification are determined first. Approximately, 1000 edges obtained from the training set (details in section Classifier Design) are marked manually as training data, and features that do not range from 0 to 1 are normalized first using a min-max normalization method. The classifiers are trained with the training data. The recognition accuracy is greater than 95 % when the SVM model is used to recognize unknown edges. Most of the recognized incorrectly edges are not the BCI or do not have cartilage coverage, which does not affect the cartilage segmentation procedure.

Feature Extraction

Based on the experiments and the literature [1, 3], the following features were finally determined for use in pixel classification: the intensity f(x, y); the distance from the BCI d(x, y) [1]; the gradient norm determined by the Prewitt operator |∇f(x, y)| and position o(x, y) [3]; and the intensity variance of the eight neighboring pixels σ₈(x, y). All features were calculated after denoising by a Gaussian low-pass filter in the section of image preprocessing. The distance from the BCI was the Euclidean distance between each pixel and the corresponding nearest point on the BCI. When a pixel was located inside the bone tissue, the distance was considered to be negative; otherwise, it was considered to be positive. The sign (positive or negative) was calculated by determining whether the point was in the polygonal region of bone formed by the BCI. When a pixel was located near the edge of the patella cartilage, the position was determined to be the horizontal distance between the pixel and the left-most point of the femoral cartilage edge. When a pixel was located near the edge of the femoral cartilage or tibial cartilage, the position was determined to be the vertical distance between the pixel and the lowest end of the femur. The position reflected the relationship between the different types of cartilage in the knee. The distance from the BCI and position could limit the range of pixel classification and was also useful in distinguishing the contacted cartilage from the muscle tissue close to the cartilage.

Classifier Design

Pixels were classified into two classes: cartilage and non-cartilage. Femoral cartilage, tibial cartilage, and patellar cartilage classifiers were built separately using the minimum error probability classifier of Bayesian decision theory. The Bayesian classifier only requires a small amount of training data to estimate the parameters necessary for classification and is faster speed and more effective and accurate compared to the SVM. More importantly, the accuracy of the classifier in segmenting our MRI images is higher than that of the SVM, which has been previously tested. Let $\mathit{x}\mathit{=}\left(x_{\mathsf{1}},x_{\mathsf{2}},\cdot \cdot \cdot ,x_n\right)$ be the feature vector described in section Feature Extraction, and ω_i is a class. If the probability of class membership in ω₁ is higher than the probability of class membership in ${\omega }_{\mathsf{2}}$, then x belongs to ω₁; otherwise, x belongs to ${\omega }_{\mathsf{2}}$ (see [13] for details).

Described mathematically,

where $P\left(x_{\mathsf{1}},x_{\mathsf{2}},\cdot \cdot \cdot ,x_n|{\omega }_i\right)$ is called the class-conditional probability distribution and P(ω_i) is the prior probabilities. In practice, $P\left(x_{\mathsf{1}},x_{\mathsf{2}},\cdot \cdot \cdot ,x_n\right)$ is identical to all classes and therefore can be ignored. In this case, the formula (1) can be rewritten as follows:

The features are considered to contribute independently to the probability, and a typical assumption is that the features associated with each class are continuous and distributed according to a Gaussian distribution. Therefore, only the feature parameters (mean and variance in each class) and the prior probabilities P(ω_i) need to be obtained to design a classifier. The feature vector is calculated first, but it is time-consuming to calculate all pixels. For example, there are more than 400,000 pixels in one case, which requires approximately 1 h to calculate the feature vector for each pixel. Thus, according to the anatomical features of cartilage, two strategies are used to reduce the number of training pixels, which reduces the computational efficiency by approximately 80 %:

1. Whether cartilage is included

   When an image does not have a BCI, there is no need to classify the corresponding cartilage.

2. Narrow the search space

Cartilage lies on the surface of bone and has a width of 1–6 mm. Only pixels less than 8 mm from the BCI (includes some non-cartilage pixels) need to be classified.

Pixel Classification and Post-processing

The same training strategies are used to reduce the number of pixels for classification when the cartilage is segmented with the three Bayesian classifiers. After pixel classification, the segmentation results still contain holes, repeated classifications, non-smooth areas, isolated points, and more than one connected regions, which does not match the anatomical features cartilage. Thus, post-processing is required.

Repeated classification occurs due to contact between some slices of patellar cartilage and femoral cartilage or between femoral cartilage and tibial cartilage. Pixels in these areas would be classified by two classifiers, and the pixels may be classified as two types of cartilage. In this situation, the pixels are classified into the class with the higher probability. Regarding the holes, non-smooth areas, and isolated points, a morphological operator (when a background pixel has at least five foreground pixels, then the pixel is the foreground pixel) is used. When a segmentation result has many connected regions, the small regions (30 pixels) are removed.