Osteoarthritis classification using self organizing map based on gabor kernel and contrast-limited adaptive histogram equalization.

Localization is the first step in osteoarthritis (OA) classification. Manual classification, however, is time-consuming, tedious, and expensive. The proposed system is designed as decision support system for medical doctors to classify the severity of knee OA. A method has been proposed here to localize a joint space area for OA and then classify it in 4 steps to classify OA into KL-Grade 0, KL-Grade 1, KL-Grade 2, KL-Grade 3 and KL-Grade 4, which are preprocessing, segmentation, feature extraction, and classification. In this proposed system, right and left knee detection was performed by employing the Contrast-Limited Adaptive Histogram Equalization (CLAHE) and the template matching. The Gabor kernel, row sum graph and moment methods were used to localize the junction space area of knee. CLAHE is used for preprocessing step, i.e.to normalize the varied intensities. The segmentation process was conducted using the Gabor kernel, template matching, row sum graph and gray level center of mass method. Here GLCM (contrast, correlation, energy, and homogeinity) features were employed as training data. Overall, 50 data were evaluated for training and 258 data for testing. Experimental results showed the best performance by using gabor kernel with parameters α=8, θ=0, Ψ=[0 π/2], γ=0,8, N=4 and with number of iterations being 5000, momentum value 0.5 and α0=0.6 for the classification process. The run gave classification accuracy rate of 93.8% for KL-Grade 0, 70% for KL-Grade 1, 4% for KL-Grade 2, 10% for KL-Grade 3 and 88.9% for KL-Grade 4.

INTRODUCTION

Osteoarthritis (OA)is the most general joint disorder, which consists of a group of heterogeneous syndromes that affect all joint tissues and are characterized by articular cartilage and the adjacent bone which often show the most notable changes. OA reflects a prominent musculoskeletal health and has been considered as a burden for socioeconomic for society worldwide, causing progressive and irreversible articular cartilage tissue damage and consequently, the failure of the joint as a whole organ [1-3]. OA usually occurs in the finger and thumb joints, particulary at the joints near fingertips and the joint connecting the thumb the hand, near the wrist, the knee, the smaller joints of the other toes, spine, and hip. In OA, the joints terminate at the bones that was previously separated by cartilage, and then rub againt each other. It causes severe pain in the joint due to damage in the tissue and the underlying bones. The joint space is at the hip because of bone spurs formation and loss of cartilage.

OA population in Indonesia in 2004 was 60,000,000 in the year 2007 [4]. The most widely used methods to assess the OA joints status are X-ray, magnetic resonance imaging (MRI) and arthroscopy [5]. Some cases were measured by osteo-CT [6]. In general, Indonesian doctors are familiar to use X-ray imaging or MRI.

Localization is the first step in OA classification. It is well-known, however, that manual selection is time-consuming, tedious, and expensive. Even if an expert radiologist or a highly trained person is available to select regions, high inter- and intra-ob-server variabilities are still possible [7]. The ultimate goal of this research is to develop an automated system for OA classification. Furthermore, detection and segmentation of joint space area are components of the automated system which are outlined in detail.

Several methods have been applied to classify knee OA, e.g., using 20 pre-selected images as the gold standard with a 150×150 window as a center of a joint in order to find the region of interest (ROI), which is a factor 10. A scan was performed afterwards with a 15×15 shifted window. The respective position, the distances between the 15×15 of the shifted window and each 15×15 pre-defined 20 joint images, were calculated using the eucledian distance. For predicting the KL-grade a simple weight nearest-neighbor rule was employed [8]. A developed method consists of image preprocessing, and delineation of cortical bone plates (active shape model) has also been proposed to perform detection using the edge based segmentation method [7]. The main feature to be extracted is the distance between femur and tibia bone through rotated images, arrays of mean intensity, rotation parameter(s), horizontally translated images and vertically translated images [9]. The latter reseach method was implemented, however, the classification result is still unable to detect OA and was classified as KL-Grade 4.

A machine vision system for OA classification was also proposed. This can be implemented on X-ray images of healthy and unimpaired knees, using CLAHE for preprocessing, the gabor kernel, row sum graph and gray level center of mass method for segmentation process, and SOM based GLCM feature extraction for the classification.

MATERIAL AND METHOD

In this work, we observed 303 unimpaired images, 50 data for training and 258 data for testing processes. All data were obtained from the Osteoarthritis Initiative (OAI) dataset for knees [[10], Fig. (). The fixed-flexion knee X-rays were acquired with the beam angle at 10 degrees. The anterior wall of the SynaFlexer positioning frame (provided by Synarc) must be in direct contact with the bucky, cassette holder or reclining table top of the radiographic unit such that there should be no angle or gap between them. To ensure these requirements, the bucky or cassette holder is lowered so that the center of the film should be at the level of the participant’s tibiofemoral joint line. The center line of the positioning frame is set to the center of the bucky or cassette holder [11].

Hardware specification used in this research is Dell Inspiron N3010 with Processor Intel(R) Core I5 CPU M430 @2,27 GHz with installed memory (RAM) 2,00 Gb. All digital radiographs were normalized into 2828×2320 images before segmentation. To classify the OA into KL-Grade 0, KL-Grade 1, KL-Grade 2, KL-Grade 3 and KL-Grade 4, there are 4 steps: preprocessing, segmentation, feature extraction, and classification. CLAHE is used for preprocessing step, i.e. for normalization of varied intensities. The segmentation process employ the Gabor kernel, template matching, row sum graph and gray level center of mass method. Here, GLCM (contrast, correlation, energy, and homogeinity) features of x-ray image have been used as training data. After completing the training process through SOM, the tests for the 258 data were run. The procedure to classify OA into the grades is presented in Fig. ().

CONTRAST-LIMITED ADAPTIVE HISTOGRAM EQUALIZATION (CLAHE)

CLAHE overcomes the limitations of standard histogram equalization [12]. There is a special case of the histogram equalization technique [13].

CLAHE technique is explained below:

Step 1: Each cell-image is divided into a number of 2+2 non-overlapping contextual region.

Step 2: Histogram of every contextual region is calculated.

Step 3: Setting the limitation of clip limit, the clip limit is a threshold parameter whereupon the contrast of the cell-image could be effectively changed.

Step 4: Each histogram is distributed in such a way that the height does not limit the clip.

Step 5: The transformation function using probability density of the input image grayscale value is used to modifyall histograms, where n is the total number of pixels in the input image and nj is the input pixel number of grayscale value j.

Step 6: The cell-image grayscale values are changed on the basis of the results of modified histograms and bilinear interpolation which are used to combine the neighboring cell-image [14].

The aim of histogram equalization mapping is to provide the input image intensity values in such a way that the histogram of the resulting images must have an approximately uniform distribution [15].

The histogram of a digital image with gray levels in the range [0, L-1] is a discrete function where ris the kth gray level, nis the number of pixels in the image with such gray level, n is the total number of pixels in the image, and k =0, 1, 2, …, L-1. Basically, p(rgives an estimate of the occurrence probability of gray level r[16].

GABOR KERNEL

The complex formula of the Gabor function in space domain is where s(x,y) is a complex sinusoid (carrier), and wr(x,y) is a two-dimensional Gaussian-shaped function (envelope).

The complex sinusoid is defined as whereupon (u0,v0) is defined as the spatial frequency and P is the phase of the sinusoid.

The real and the imaginary parts of the sinusoid are:

Parameters u0 and v0 define the spatial frequency of the sinusoid in the Cartesian coordinates. The spatial frequency expressed in polar coordinates in terms of magnitude F0 and direction ω0 is as follows: The complex sinusoid is then: The gaussian envelope is as follows: where K is the scales for the magnitude of the gaussian envelope, (x) is the peak of the function, a and b are scaling parameters of the gaussian envelope, while r denotes the rotation operator [17].

Whereas:

Fig. () shows an example of Gabor kernel with several N parameter values.

ROW SUM GRAPH

We use row sum graphs to determine the junction area. The sum graph represents the row-wise summation of the gray values [18, 19]. Considering the input image be F(i, j), the row sum graph S(i, j) can be given by

GREY LEVEL CENTER-OF-MASS

The center of mass can be determined in terms of gray level in the localized junction space area, whereas the junction area can be determined by the given formula

FEATURE EXTRACTION

GLCM (also called gray tone spatial dependency matrix) is a tabulation of the frequencies. Different combinations of pixel brightness values (gray levels) occur in an image [20]. The term of ‘GLCM texture’ has been introduced to recall the relation between two pixels, i.e. the reference and the neighbor pixels, at a time.

There are two occurrences of pixel intensity 0 with pixel intensity 1 as its neighbor (in the horizontal direction or the direction of 0°). Therefore, prior to calculate texture measures the required measures are that each GLCM cell should contain a probability. This process is called normalizing the matrix. Normalization equation is as follows: where i is the row number and j is the column number.

Haralick [21] extracted 14 features from the co-occurrence matrix, but in this research only 4 features, namely the contrast, correlation, energy and homogeneity with 4 directions and distance d=1 have been extracted. The features considered are:

SELF ORGANIZING MAP (SOM)

SOM has been progressing rapidly both in its substances and applications to cluster in various fields [22-26]. Yet, there are a large number of researches about other SOM. More in-depth analysis of SOM algorithm can be found other reference [16]. In principle, the algorithm of SOM can also be illustrated through Fig. (). The steps of learning to use SOM are initialization, determining the number of classes, setting a learning rate parameter, determination of the iteration value, and calculating the distance between random data into each input weight.

The classic SOM network learning algorithm can be formulated as follows [[27], Fig. ():

Setting learning iteration number (t), initializing all weights (wij) with small random values, setting the initial topological neighborhood (d0), determining learning rate (α0), and determining the total number of iteration (T).

If iteration number (t) is less than the total number of iterations (T), then Step 3 and Step 6 should be repeated.

Choosing a random input vectors in the training set.

Finding the winner neuron by determining the neuron j so that its weight vector wj should become closest to the input vector according to

The weight vectors of both the winner neurons j and its neighbors are updated as follows: where α(t) is a learning rate function that exponentially decreases with time. Moreover, a neighborhood order function d(t) has been defined which exponentially decreases with time. In this paper, the following equations are used for α(t) and d(t):

Setting t = t+1, if t

EVALUATION

A receiver operating characteristic (ROC) curve is most frequently used for illustrating the performance of a system because of its comprehensive and fair evaluation ability [18]. An ROC curve is a plot of true positive fraction (TPF) as a function of false positive fraction (FPF) [19, 20].

The area under the ROC curve (AUC) can be used as a criterion. Table shows the classifying level of accuracy based on AUC [21].

where TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives and FN is the number of false negatives.

RESULTS AND DISCUSSIONS

At present, techniques have been developed to classify knee OA x-ray image into KL-Grade 0, KL-Grade 1, KL-Grade 2, KL-Grade 3 and KL-Grade 4. The segmentation algorithm has been successfully implemented for knee X-ray image. Following are the results of the implementation.

THE FIRST STAGE: PREPROCESSING

The results for knee OA are shown in Fig. (). Fig. () is the original image and Fig. () is the result of CLAHE. The CLAHE was used to improve the performance of X-ray image. This experiment was designed to normalize the potential for varying intensities of data. CLAHE method seems to have the potential to be applicable in knee OA dataset, on the basis of the experiments CLAHE was found to normalize the whole image therefore it can be further processed. Without using CLAHE, the intensities of image were varying enough that it became difficult to proceed.

THE SECOND STAGE: SEGMENTATION

The position of joint space was detected by the gabor filter, template matching and row sum graph process. Gabor filter is particularly suitable for texture segmentation because it possess various properties. Experiments have shown that Gabor function can be tuned to a narrow set of frequency anywhere in the frequency domain.

The Gabor experiment is applied to the preprocessed images. More precisely, the output of Gabor channels with the following settings were applied with parameters α=8, θ=0, Ψ[0 π/2], γ=0,8 with N variation value 1 to 9. The results of the preprocessing using the gabor kernel are presented in Fig. (), where some textures are detected as part of the dot pattern texture along with the boundary but are still plausible.

The gabor filter was tuned to a specific frequency and direction in order to obtain the local frequency and orientation information so that it optimally captures both local orientation and frequency information from the image. Table shows that experiment 4 with parameter α=8, θ=0, Ψ=[0 Π/2], γ=0,8 and N=4 gives an accuracy of 95.12% for the right knee and 76., 83% for the left knee.

The row sum graphs were used to determine the junction area. The sum graph represents the row-wise summation of the gray values [18]. The result of row sum graph was based on the gabor kernel with parameters α=8, θ=0, Ψ=[0 Π/2], γ=0,8 and N=9 is shown in Fig. () for right knee and Fig. () for the left knee. The peak (maximum value) is detected in the junction space area.

There were some unsuccessful images detected to provide the junction space area. The failures were due to the presence of screws or bond plates which were detected as junction space area. Another failure in junction detection was might be due to noise and deviations of the angle of the knee joint from the horizontal.

THE THIRTH STAGE: FEATURE EXTRACTION

The 16 X-ray image features were exracted with the help of GLCM, with each direction 0°, 45°, 90°, and 135° and distance = 1. Minimum, average, and maximum values are described in Table . On the basis of the table, some extracted features for KL-Grades 0, 1, 2, 3 and 4 were overlapped. Due to this reason it becomes imposible to make a linear separation between each class, therefore, a solution was proposed to use SOM to classify KL-Grades 0, 1, 2, 3, and 4.

THE FOURTH STAGE: CLASSIFICATION

The SOM algorithm has been applied and the obtained results are briefly presented in this section. SOM parameters used in this research are:

number of iterations: 5000

momentum value: 0.5

learning rate formula: where α is the learning rate, α0=0.6, j is the epoch and is the number of total iterations.

Classification results are displayed by a confusion matrix in Table , where KL-Grade represents the grades given by human experts, and result denotes that SOM network outputs have been compared. In the matrix, each cell contains the number of samples classified for the corresponding SOM classified outputs. The value of cell (i,j) represents the number of classifications of OA determination class j into class i. If all the entries are located along the main diagonal, then it shows that the classification results of SOM are altogether according to human experts. KL-Grade 0 classified as KL-Grade 0 was 93.75% and classified as KL-Grade 2 was 6.25%. While KL-Grade 1 classified as KL-Grade 0 was 28.33%, and classified as KL-Grade 2 was 1.67%.

Accuracy values for each KL-Grade are presented in Fig. (). The accuracy values vary with the Grade, i.e. 93.8% for the KL-Grade 0, 70% for KL-Grade 1, only 4% for KL-Grade 2, 10% for KL-Grade 3 and finally, 88.9% for KL-Grade 4. Accuracy for KL-Grade 2 was only 4% because 60% of the result was classified as KL-Grade, 1, 6% was classified as KL Grade-3 and 30% was classified as KL-Grade 4.

Accuracy evaluation performance by statistical prediction model can also be done by Receiver Operating Characteristics (ROC) curve analysis. ROC curve is a graphical plotting with the y-axis expressing sensitivity (true positive rate) and the x-axis expressing false positive rate. Fig. () showed the ROC curve for the classification using SOM on the basis of GLCM.

The ROC curve contains most of the information about the accuracy of a continuous predictor; it is sometimes desirable to produce quantitative summary measures of the ROC curve. The most common implementation of such measure by far is the area under the ROC curve (AUC), where its values are shown in Fig. (). Based on Table and Fig. 10, it can be inferred that KL-Grades 0 and 4 are excellently classified, while KL-Grade 1 is well classified and KL-Grade 2 and KL-Grade 3 are failures for classification.

Fig. (10)

Graphic of AUC value for each group. The highest AUC is 0.978 from KL Grade 4, AUC value of KL Grades 0, 1, 2, 3 are 0.962, 0.8, 0.075, 0.142 and 0.978.

CONCLUSIONS

A procedure has been described to classify knee OA using SOM into KL-Grade 0, KL-Grade 1, KL-Grade 2, KL-Grade 3 and KL-Grade 4. The proposed technique employed the morphological process based on CLAHE and the gabor kernel, where it has then been successfully implemented. Experimental results have been provided to show the effectiveness of the proposed technique, and may act as a support system for medical doctors to determine the region of interest of visual characteristics found in an OA knee and hence as a decision support system to classify the severity of such OA.

Data learning comprised of 16 features and based on GLCM (contrast, correlation, energy, and homogeneity, each for direction 0o, 45°, 90°, and 135° and distance = 1) and the KL-Grade manually classified by expert was used as the standard for testing process.

Experimental result have shown the best performance when using number of iterations being 5000, momentum value as 0.5 and α0=0.6 which give classification accuracy rate of 93.8% for KL Grade 0, 70% for KL-Grade 1, 4% for KL-Grade 2, 10% for KL-Grade 3 and 88.9% for KL-Grade 4. On the basis of the AUC values, it was found that KL-Grade 0 and KL-Grade 4 were classified as exellent, KL-Grade 1 as good and KL-Grade 2 and KL-Grade 3 wre failures.

The accuracy of KL-Grade 2 and KL-Grade 3 classification cannot be considered strong, because the radiograph is often difficult to distiguish between these grades. KL-Grade 4 and KL-Grade 0 were differentiated accurately from other grades. It is considered that this system will have an impact on the decision support system when determining the OA classification.

Table 1.

Classifying Level of Accuracy Based on AUC

AUC Value	Classified as
0.90 – 1.00	Excellent
0.80 – 0.90	Good
0.70 - 0.80	Fair
0.60 – 0.70	Poor
0.50 – 0.60	Fail

Table 2.

Accuracy of The Segmentation Stage

Experiment	Parameter	Accuracy
		Right Knee	Left Knee
1	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=1	2.44	7.32
2	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=2	7.32	10.98
3	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=3	58.54	48.78
4	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=4	95.12	76.83
5	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=5	90.24	71.95
6	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=6	68.29	59.76
7	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=7	87.8	68.29
8	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=8	91.46	75.61
9	α=8, θ=0, Ψ=[0 π/2], γ=0,8 and N=9	87.8	71.95

Table 3.

Tabulated Minimum, Average and Maximum Feature Value for KL Grades 0, 1, 2, 3 and 4 Extracted from GLCM

Feature	Maximum value					Average Value					Maximum Value
	Grade					Grade					Grade
	0	1	2	3	4	0	1	2	3	4	0	1	2	3	4
Direction = 00, Distance = 1
Contrast	0,364	0,426	0,367	0,338	0,364	0,493	0,586	0,480	0,448	0,544	0,633	0,718	0,715	0,697	0,897
Correlation	0,890	0,864	0,868	0,875	0,850	0,904	0,887	0,900	0,909	0,887	0,917	0,921	0,934	0,945	0,918
Energy	0,081	0,067	0,085	0,067	0,059	0,105	0,084	0,112	0,126	0,123	0,149	0,118	0,146	0,191	0,203
Homogeneity	0,781	0,753	0,776	0,753	0,731	0,802	0,776	0,804	0,812	0,797	0,832	0,796	0,835	0,841	0,849
Direction = 450, Distance = 1
Contrast	0,803	0,778	0,689	0,604	0,342	1,292	1,566	1,113	0,967	1,032	1,809	2,168	2,114	1,912	1,709
Correlation	0,630	0,582	0,616	0,666	0,680	0,749	0,700	0,770	0,804	0,791	0,841	0,854	0,879	0,899	0,922
Energy	0,051	0,044	0,058	0,046	0,045	0,068	0,060	0,081	0,090	0,097	0,096	0,090	0,102	0,147	0,205
Homogeneity	0,654	0,639	0,646	0,649	0,644	0,689	0,676	0,710	0,723	0,719	0,727	0,740	0,757	0,779	0,852
Direction = 900, Distance = 1
Contrast	0,506	0,653	0,489	0,364	0,222	1,046	1,290	0,875	0,754	0,836	1,506	1,967	1,881	1,617	1,497
Correlation	0,693	0,650	0,695	0,737	0,744	0,816	0,773	0,838	0,863	0,849	0,917	0,892	0,923	0,944	0,958
Energy	0,056	0,043	0,058	0,050	0,050	0,073	0,064	0,084	0,095	0,100	0,104	0,090	0,109	0,156	0,198
Homogeneity	0,690	0,645	0,673	0,691	0,674	0,736	0,720	0,755	0,775	0,760	0,792	0,777	0,815	0,840	0,904
Direction = 1350, Distance = 1
Contrast	0,863	0,915	0,698	0,613	0,570	1,335	1,604	1,137	1,020	1,160	1,815	2,447	2,296	2,028	2,006
Correlation	0,630	0,543	0,604	0,644	0,629	0,743	0,695	0,767	0,793	0,767	0,831	0,817	0,878	0,899	0,873
Energy	0,051	0,038	0,053	0,043	0,044	0,067	0,057	0,077	0,089	0,092	0,104	0,085	0,102	0,142	0,172
Homogeneity	0,645	0,581	0,628	0,625	0,634	0,674	0,664	0,695	0,715	0,707	0,731	0,712	0,738	0,755	0,805

Table 4.

Classification Confusion Matrix of System

	K-L Grade					Total
	0	1	2	3	4
KL Grade 0	93,75	6,25	0,00	0,00	0,00	100,00
KL Grade 1	28,33	70,00	1,67	0,00	0,00	100,00
KL Grade 2	0,00	60,00	4,00	6,00	30,00	100,00
KL Grade 3	12,00	34,00	2,00	10,00	42,00	100,00
KL Grade 4	5,56	5,56	0,00	0,00	88,89	100,00