# Automatic Diagnosis of Elbow Arthritis Based on Edge Algorithm

Authors: Minwei Xia; Peng Ao; Bin Zhang; Yongjun Liao; Huixue Zhao
Journal: Journal of Sensors
 (2022)

---

## Abstract

Osteoarthritis is an age-related degenerative joint disease; it is mainly because the cartilage tissue between bones is worn and thinned, which leads to the damage of the periosteum and bone including the surrounding ligaments. Clinically, its manifestations are mainly joint pain, swelling, stiffness, and even partial loss of function, which seriously affects the quality of life of patients. The main clinical manifestations are elbow joint pain and limited movement. Elbow articular cartilage degenerates and falls off, and the more serious manifestation is subchondral hyperosteogeny and sclerosis, which leads to unsmooth articular surface and narrow joint space. Finally, elbow joint pain is severe with different degrees of mobility disorder, elbow joint extension and flexion range is getting smaller and smaller, and elbow joint pain is getting more and more serious. In this paper, the segmentation of left and right elbow images is completed based on gray projection through the analysis of image gray distribution. After obtaining the region of interest of elbow joint, the extraction algorithm of elbow joint hard bone edge is studied. Firstly, the extraction of elbow joint hard bone contour edge is completed based on active shape model algorithm combined with image characteristics. Finally, according to the extraction results of hard bone contour edge, this paper realizes the automatic diagnosis of multiple elbow arthritis indexes and compares with the results given by the image set, which proves that the whole algorithm has good adaptability and accuracy.


---

## Body

## 1. Introduction

In recent years, with the rapid development of imaging, various modern medical imaging techniques have been widely used in the diagnosis of elbow osteoarthritis, such as X-ray, ultrasound, CT, MRI and so on [1]. These imaging techniques have their own areas of expertise. Sometimes the diagnosis of disease needs to combine multiple imaging techniques to analyze different images of the same disease site. However, X-ray imaging equipment is cheap, and X-ray films are suitable for basic diagnosis, so it is still the most common to observe X-ray films [2].In today’s medical field, both scientific research and clinical application are inseparable from image processing and recognition technology [3]. In the low-level application field, the current medical image-aided diagnosis system has realized simple processing such as image enhancement, geometric transformation, and lesion location selection. These processes enable doctors to see the characteristics of lesions more clearly and improve the diagnosis rate. Furthermore, in the middle-level application, in the field of automatic segmentation of organs and tissues, many techniques in image processing have been applied and achieved good results, such as threshold-based segmentation method, region-based segmentation method, and graph theory-based segmentation method [4]. However, due to the complexity of medical images, the edges of organs and tissues are often blurred, so it is difficult to have a unified segmentation algorithm that can segment all organs and tissues with good results [5]. In the aspect of 3D reconstruction, ITK toolbox, which is developing rapidly at present, can complete the conversion from 2D image sequence to stereo image. This stereo vision conversion technology can not only visually observe the interested tissues but also rotate the three-dimensional structure, so that medical staff are no longer limited to observing the lesions from a certain angle but observe them in all directions to make more accurate judgments [6]. Although image processing technology has made great progress in the field of medical image, most of the technologies are still difficult to be called intelligent, and it needs to be further explored in the advanced application field of intelligent diagnosis.


## 2. Research Results of Osteoarthritis

The research on automatic diagnosis of osteoarthritis by using image processing and recognition technology originated in the 1980s. The initial research is mainly based on traditional image processing methods to complete elbow joint region extraction or effective edge detection. Common algorithms include Sobel edge detection, binarization based on area threshold, and so on. However, the difficulty of this kind of image segmentation lies in the blurred edges, and the boundaries between organizations are not obvious [7]. It is difficult to obtain real edges by using traditional digital image processing algorithms alone.Later, with the development of intelligent algorithms, the research direction of elbow arthritis diagnosis based on image information is divided into two directions: one direction is to study more accurate algorithms to achieve accurate edge extraction of hard bone and cartilage of elbow joint. This direction believes that as long as accurate edges of hard bone and cartilage can be extracted, a series of related parameters can be obtained, and it is not difficult to complete the diagnosis of elbow arthritis according to these parameters [8]; on the other hand, the diagnosis and analysis of elbow arthritis is obtained by doctors according to a large number of film reading experiences, and its diagnostic criteria are vague, so it is not necessary to obtain various index parameters [9]. The algorithm is successful as long as it can classify different symptoms of elbow joint images.The first direction focuses on region segmentation and edge extraction. In a broad sense, edge extraction and region segmentation are both part of image segmentation algorithm. There are many algorithms for image segmentation. From the initial threshold segmentation to the recent machine learning segmentation methods, threshold segmentation algorithm is a basic segmentation algorithm in image segmentation, which can achieve good results in extracting the region of interest of elbow joint. With the progress of image algorithm research, segmentation methods based on region growing method and active contour model have made new progress. Some scholars have realized the cartilage region segmentation in nuclear magnetic resonance images by using the image segmentation method based on region growing method. In the field of medical image segmentation, region growing method is often combined with traditional edge detection methods, which complement each other and often make the segmentation effect more ideal. Active contour model is used to segment cartilage region in 3D image [10]. Active contour model is based on the concept of energy functional to find the boundary curve of region segmentation. According to different definitions, it can be divided into geometric active contour active model and parametric active contour model [11].The research in the second direction comes into being with the rise of deep learning research, which does not advocate edge extraction and does not need to set features artificially. Instead, with the help of neural network’s powerful learning ability, it can classify images directly. Scholars put forward a new idea: abandoning the features of manual design and assuming that computers can observe weaker changes than human eyes, so as to realize data-driven classification. In this paper, several statistical features are specified to observe their sensitivity to the diagnosis of elbow arthritis. After calculating all the feature values of a given test image, the obtained feature vectors are classified by using a simple weighted nearest neighbor rule, thus realizing the diagnosis of elbow arthritis [12]. Some scholars have published a paper on quantifying the severity of elbow arthritis by using deep convolution neural network. This paper introduces the research results of deep learning into the diagnosis of elbow X-ray films and adopts the mean square error loss function, taking 70% of the pictures in the data set as training sets and 30% as test sets.According to the diagnostic requirements of elbow arthritis, the problem of the first research direction is that the accuracy of edge extraction is not enough. The research index mainly selects elbow joint space distance as an index to diagnose the severity of elbow arthritis. It is undeniable that the distance between elbow joints is the most important index to judge the severity of elbow arthritis, but besides, the amount of osteophyte, edge sclerosis, and meniscus cartilage calcification also plays a vital role in the diagnosis of elbow arthritis, which should not be ignored when designing the diagnosis system. The second research direction lies in the coupling of multiple indicators. Although the classification accuracy is improved through deep learning algorithm, it is more important for clinical diagnosis to obtain various indicators.


## 3. Automatic Diagnosis Model of Elbow Joint Based on Edge Algorithm

### 3.1. Image Preprocessing

Due to the difference of X-ray energy or the narrow nonlinear dynamic range of imaging equipment, the exposure degree of images is different, the exposure of some images is too unbalanced, the whole image is dark, and the gray distribution range is narrow, which is difficult to identify visually and seriously affects the subsequent processing links. In order to make the visual effect of the image stronger and the detail resolution easier, this paper extends the gray level of the image to 0-255. There are many common gray level transformation algorithms, such as linear transformation, exponential transformation, logarithmic transformation, and histogram equalization. Because the gradient of each point must be calculated in the later processing, if the gray level transformation algorithm such as nonlinear or histogram equalization is adopted, it is easy to change the relative gradient of the image. In order to keep the gradient of each pixel relatively invariant, this paper chooses to use linear transformation to stretch the gray level of the image to 0-255, which is described as follows.If the gray distribution interval of the original image ism,n and the transformed gray distribution interval is s,t, the formula can be used:
(1)gx,y=t−sfx,y−mn−m+s.In image processing, the core of Gaussian smoothing is a two-dimensional convolution kernel using the characteristics of Gaussian function. The convolution kernel is used to suppress noise and blur the image. According to the characteristics of Gaussian distribution, the value of the central pixel in Gaussian distribution is the largest, so its weight ratio is also the highest. With the increase of the distance from the central pixel, the weight ratio of other pixels in the neighborhood of the central pixel is getting lower and lower. In theory, all pixels around the center point need to contribute weights, but the window size of convolution kernel is limited in actual calculation. According to the definition of small probability events in probability statistics, the contribution of pixels outside the center of the mean value can not be considered. Therefore, in image processing, the same idea can be adopted to intercept the square with finite side length centered on the point to be smoothed as the boundary area to lock those pixels that contribute weights to the pixels to be smoothed. In this paper, Gaussian convolution kernels of different sizes are repeatedly tried in the processing process to determine the optimal convolution kernel for the subsequent algorithm. Finally, it is found that when the Gaussian convolution kernel is 7∗7, the processing effect for the subsequent algorithm is the best. The convolution kernel is shown in Figure 1.Figure 1
Gaussian convolution kernel.


### 3.2. Automatic Diagnosis Model of Elbow Joint Based on Edge Algorithm

#### 3.2.1. Sample Selection and Feature Labeling

In the segmentation method based on active shape model, the construction of initial model is very important, and to build a representative initial model, it is necessary to select appropriate samples, so the selection of samples is the foundation and an important step. First, selectn sample images as training samples. The selection principle should follow diversity and representativeness; for example, the orientation of elbow joint is oblique to the left and oblique to the right and vertical. The position of elbow joint is left, right, and middle, and the size of elbow joint is also slightly different. Various diseases of elbow joint lead to the change of elbow joint space and edge, all of which should be fully considered to avoid the limitation of sample selection. After selecting samples, we need to label the feature points of these samples and manually mark K key feature points of the target object in N samples, and the selection of feature points is based on the contour edge. Because the edge detected by Canny edge detection is the finest and the integrity is relatively good, the results of Canny edge detection are used to assist the selection of feature points. Taking the sampling and selection of feature points at the lower edge of femur as an example, sampling should not only ensure that the number of sampling points can describe the shape of the contour but also avoid the redundancy of information caused by too many sampling points. Generally, the points with high curvature and T-shaped joints are selected first, and the other feature points are arranged among the above points at reasonable intervals in order. As shown in Figure 2, the selection process of feature points of elbow joint femoral edge is described.Figure 2
Sampling schematic diagram of femoral edge of elbow joint.In Figure1, the larger dots are typical feature points, and the smaller dots are feature points added equidistantly between typical feature points to depict a complete shape. After selecting feature points, the coordinates of K feature points in an image are sequentially formed into a shape vector such as formula
(2)Xi=xi1,xi1,xi2,xi2,⋯⋯,xik,xik,where xik,xik is the coordinates of the k-th feature point in the i-th image.Thus, the length of each shape vector is 2 K, andn samples can form n shape vectors, thus obtaining a sample set J=X1,X2,⋯⋯Xn.


#### 3.2.2. Sample Shape Alignment and Matching

The shapes of elbow joints are quite similar. For example, there will be an upward groove in the middle of the lower end of the femur, and the medial femur is usually larger than the lateral femur and so on. However, due to the slight inclination angle and orientation change of femur in training samples, in order to eliminate the interference of nonshape information caused by the standing posture and bone size of the subject, it is necessary to alignn shape vectors, which is called aligning shape vectors. The method adopted is Procrustes alignment. Generally speaking, general alignment is an affine transformation between shapes sought by using the idea of least square method. The basic steps are as follows:
(1)
After aligning to the origin, the shape vectors of the target objects in all training images are aligned to the shape vectors of the target objects in the first image, and the first shape vector should be normalized to the unit size(2)
Calculating the average valueX¯=1/n∑i=1nXi of the shape vector of the target object(3)
Aligning all the shape vectors generated in the first step toX(4)
Repeat steps 2 and 3 until convergence (the difference between two adjacent times is within a given threshold)E the expression to get the minimum value:
(3)E=Xi−Fs,θXj−mTWXI−Fs,θXj−m.BecauseFS,O stands for rotation and contraction transformation, we can deduce the formula
(4)Fs,θxkyk=scosθxk−ssinθykscosθxk−ssinθyk.The diagonal element of the weight diagonal matrixW is the weight Wk of each marker point; it represents the stability of the marking point. Taking the femoral edge as an example, some characteristic points are relatively stable, such as the groove vertex in the middle of the lower femur and the turning point on the left and right sides. Improving the weight of these points can effectively improve the accuracy of matching. Mark the distance of the k-th feature point and the i-th feature point in an image as Rki (still in the order of the origin mark), and the n images can get n distances, and VRki represents the variance of nRki, then
(5)Wk=∑i=1nVRki−1.Order, available formula is as follows :ax=scosθ,ay=ssinθ(6)X2−Y2W0Y2X20WZ0X2Y20Z−Y2X2axaymxmx=X1Y1C1C2.


#### 3.2.3. Dimension Reduction of PCA

PCA is principal component analysis, assuming that the number of feature points in each sample in the training set isk; then, a 2k-dimensional space is needed to process these shape vectors. Generally speaking, when K is large, it is difficult to deal with such a high-dimensional vector. In order to reduce the difficulty of processing, it is necessary to select some important features to describe the shape, that is, dimension reduction. The purpose of dimension reduction is to greatly reduce the redundancy of information while retaining most effective information. The most commonly used dimension reduction method is PCA (principal component analysis) processing, and the steps are as follows.Calculate the average shape vector:(7)X¯=1n∑i=1nXi.Calculate the covariance matrix:(8)S=1n−1∑i=1nXi−X¯TXi−X¯.The eigenvalueλi and the corresponding eigenvector pi of the covariance matrix are found, and the eigenvalues are arranged in a decreasing order, and the eigenvector corresponding to the previous eigenvalue will have the main effect on the shape vector. We select the first t eigenvectors to satisfy the formula
(9)∑i=1tλi∑i=1mλi>k.Selecting feature vectors like this can keep the original shape and reduce the computational complexity. In this way, any shape vector participating in training in the training image can be expressed as a linear combination of feature vector sets, as shown in the formula:(10)Xi=X¯+PB.


#### 3.2.4. Constructing Local Gray Model

Principal component analysis adjusts the shape as a whole. Next, the shape needs to be iterated to the real contour edge, and each iteration needs to find a new position for each feature point, while the local gray model can provide the direction of iterative development. Here, the concept of gray gradient is used to define local features for each feature point. For thei-th feature point, the creation process of local features is shown in Figure 3.Figure 3
Description diagram of local feature creation.In thej-th training image, m pixels on each side of the i-th feature point are selected on the normal line of the i-th feature point (here, the normal refers to the vertical line connecting the previous feature point of the feature point and the second feature point), thus forming a vector with 2m+1 component. As shown in the following equation, each pixel on this vector is first-order guided, and the result of the operation is called local texture gij. The i-th feature point on the other images in the training set does the same work, and all the local texture about the i-th feature point gi1,gi2,⋯⋯,gin.
(11)gij=gij1,gij2,⋯⋯,gij2m+1T.


#### 3.2.5. Model Initialization Positioning

When the model covers the target image, it needs initialization positioning. The accuracy and efficiency of ASM algorithm depend on the accuracy of initialization positioning. If the initialization positioning is far away from the real contour edge, the target may not be searched. If the initial location is close to the real contour edge, it can not only increase the accuracy of matching but also improve the speed of search, so it is necessary to transform the average shape vector properly to get a reasonable initial position.The standardized model uses initialization and positioning methods, and generally uses the maximum threshold method to binarize the segmented image shown in Figure4. The average shape model is randomly overlaid on the binary image, but to ensure that the longitudinal lowest point of the average shape model should not be lower than the elbow joint space, then, the head and tail point vectors of the average shape model are connected to form a closed area, and the average shape model is rotated and translated to make the white pixels in the closed area account for the largest proportion of the total pixels in the area. At this time, it can be considered that the initialization and positioning of the model are completed. The operation of rotation is limited by angle, while the operation of translation is unlimited in left and right directions and limited in up and down directions.Figure 4
Maximum binary diagram.The initial location of the posterior edge of the superior tibial plane is similar to that of the inferior femoral edge, while the initial location of the anterior edge of the superior tibial plane relies on the posterior edge and iterates on the basis of the posterior edge.


#### 3.2.6. Local Search Processing

After the preliminary model is established, the initial ASM model is covered on the image, and the next step is to find the new iterative position of the feature points by using the local gray model of the feature points. The schematic diagram of the search process is shown in Figure5.Figure 5
Schematic diagram of search process.The big mark points on the outer side of the graph are the initial feature points, and the small mark points on the lower side of the graph are the best matching points of the model. The normal direction of the feature points refers to the perpendicular line connecting the front and back feature points of the initial feature points, and the iteration of each feature point moves along the normal direction. The thick line segment on the left side of the normal is2m+1 pixels selected in the local gray model, and the thick line segment on the right side of the normal represents the search area when the feature points are iterated. The distance between the initial point and the best matching point is the distance dX that the initial feature points should move.A local feature is obtained by normalization, which contains2v−m+1 sublocal features in total. Because the normalized gray vector approximately obeys Gaussian distribution, Mahalanobis distance can be used to express the similarity between the new feature g of a feature point and its trained local feature. The calculation formula of Mahalanobis distance is as follows:
(12)Fs=g−gi¯Si−1g−gi¯T.The smaller the Mahalanobis distance, the greater the similarity, and the closer the boundary pointI in the model is to the boundary point I in the target, so the best position of the boundary point can be determined by this method. Then, calculate the Mahalanobis distance between these sublocal features and the local features of the corresponding feature points in the gray model.When the Mahalanobis distance is minimized, the center point of the corresponding sublocal feature is the new position of the current feature point, and all feature points undergo the same steps, so each point will have a displacement, and these displacements will be arranged into vectors, such as the formula:(13)DX=DX1,DX2,DX3,⋯⋯,DXk.After all feature points are updated iteratively, a new curve is formed.


#### 3.2.7. Active Shape Model Parameter Update

After the initialization of the shape model and the local search strategy are established, the shape can be updated. Assume that the mapping function from the average shape vector to the initial model is shown by the following formula:(14)X=Ms,θx+m.The obtained model is further similar to the target model, and further matching is carried out by iterative method. The iterative updating process of parameters is as follows:(15)mx⟶mx+dmx,my⟶my+dmy,θ⟶θ+dθ,s⟶s+ds,B⟶B+dB.As the iteration progresses, when the difference between the two shape models is less than a certain threshold or reaches a given number of iterations, the algorithm ends, stops searching, and thinks that the model converges, and the whole search and matching process ends, and the feature points at this time are the contour edges to be searched. Figure6 shows the detection results of the lower edge of femur after different iterations, which are 3 iterations, 15 iterations, and 40 iterations from left to right.Figure 6
Detection results of iterative process.


## 3.1. Image Preprocessing

Due to the difference of X-ray energy or the narrow nonlinear dynamic range of imaging equipment, the exposure degree of images is different, the exposure of some images is too unbalanced, the whole image is dark, and the gray distribution range is narrow, which is difficult to identify visually and seriously affects the subsequent processing links. In order to make the visual effect of the image stronger and the detail resolution easier, this paper extends the gray level of the image to 0-255. There are many common gray level transformation algorithms, such as linear transformation, exponential transformation, logarithmic transformation, and histogram equalization. Because the gradient of each point must be calculated in the later processing, if the gray level transformation algorithm such as nonlinear or histogram equalization is adopted, it is easy to change the relative gradient of the image. In order to keep the gradient of each pixel relatively invariant, this paper chooses to use linear transformation to stretch the gray level of the image to 0-255, which is described as follows.If the gray distribution interval of the original image ism,n and the transformed gray distribution interval is s,t, the formula can be used:
(1)gx,y=t−sfx,y−mn−m+s.In image processing, the core of Gaussian smoothing is a two-dimensional convolution kernel using the characteristics of Gaussian function. The convolution kernel is used to suppress noise and blur the image. According to the characteristics of Gaussian distribution, the value of the central pixel in Gaussian distribution is the largest, so its weight ratio is also the highest. With the increase of the distance from the central pixel, the weight ratio of other pixels in the neighborhood of the central pixel is getting lower and lower. In theory, all pixels around the center point need to contribute weights, but the window size of convolution kernel is limited in actual calculation. According to the definition of small probability events in probability statistics, the contribution of pixels outside the center of the mean value can not be considered. Therefore, in image processing, the same idea can be adopted to intercept the square with finite side length centered on the point to be smoothed as the boundary area to lock those pixels that contribute weights to the pixels to be smoothed. In this paper, Gaussian convolution kernels of different sizes are repeatedly tried in the processing process to determine the optimal convolution kernel for the subsequent algorithm. Finally, it is found that when the Gaussian convolution kernel is 7∗7, the processing effect for the subsequent algorithm is the best. The convolution kernel is shown in Figure 1.Figure 1
Gaussian convolution kernel.


## 3.2. Automatic Diagnosis Model of Elbow Joint Based on Edge Algorithm

### 3.2.1. Sample Selection and Feature Labeling

In the segmentation method based on active shape model, the construction of initial model is very important, and to build a representative initial model, it is necessary to select appropriate samples, so the selection of samples is the foundation and an important step. First, selectn sample images as training samples. The selection principle should follow diversity and representativeness; for example, the orientation of elbow joint is oblique to the left and oblique to the right and vertical. The position of elbow joint is left, right, and middle, and the size of elbow joint is also slightly different. Various diseases of elbow joint lead to the change of elbow joint space and edge, all of which should be fully considered to avoid the limitation of sample selection. After selecting samples, we need to label the feature points of these samples and manually mark K key feature points of the target object in N samples, and the selection of feature points is based on the contour edge. Because the edge detected by Canny edge detection is the finest and the integrity is relatively good, the results of Canny edge detection are used to assist the selection of feature points. Taking the sampling and selection of feature points at the lower edge of femur as an example, sampling should not only ensure that the number of sampling points can describe the shape of the contour but also avoid the redundancy of information caused by too many sampling points. Generally, the points with high curvature and T-shaped joints are selected first, and the other feature points are arranged among the above points at reasonable intervals in order. As shown in Figure 2, the selection process of feature points of elbow joint femoral edge is described.Figure 2
Sampling schematic diagram of femoral edge of elbow joint.In Figure1, the larger dots are typical feature points, and the smaller dots are feature points added equidistantly between typical feature points to depict a complete shape. After selecting feature points, the coordinates of K feature points in an image are sequentially formed into a shape vector such as formula
(2)Xi=xi1,xi1,xi2,xi2,⋯⋯,xik,xik,where xik,xik is the coordinates of the k-th feature point in the i-th image.Thus, the length of each shape vector is 2 K, andn samples can form n shape vectors, thus obtaining a sample set J=X1,X2,⋯⋯Xn.


### 3.2.2. Sample Shape Alignment and Matching

The shapes of elbow joints are quite similar. For example, there will be an upward groove in the middle of the lower end of the femur, and the medial femur is usually larger than the lateral femur and so on. However, due to the slight inclination angle and orientation change of femur in training samples, in order to eliminate the interference of nonshape information caused by the standing posture and bone size of the subject, it is necessary to alignn shape vectors, which is called aligning shape vectors. The method adopted is Procrustes alignment. Generally speaking, general alignment is an affine transformation between shapes sought by using the idea of least square method. The basic steps are as follows:
(1)
After aligning to the origin, the shape vectors of the target objects in all training images are aligned to the shape vectors of the target objects in the first image, and the first shape vector should be normalized to the unit size(2)
Calculating the average valueX¯=1/n∑i=1nXi of the shape vector of the target object(3)
Aligning all the shape vectors generated in the first step toX(4)
Repeat steps 2 and 3 until convergence (the difference between two adjacent times is within a given threshold)E the expression to get the minimum value:
(3)E=Xi−Fs,θXj−mTWXI−Fs,θXj−m.BecauseFS,O stands for rotation and contraction transformation, we can deduce the formula
(4)Fs,θxkyk=scosθxk−ssinθykscosθxk−ssinθyk.The diagonal element of the weight diagonal matrixW is the weight Wk of each marker point; it represents the stability of the marking point. Taking the femoral edge as an example, some characteristic points are relatively stable, such as the groove vertex in the middle of the lower femur and the turning point on the left and right sides. Improving the weight of these points can effectively improve the accuracy of matching. Mark the distance of the k-th feature point and the i-th feature point in an image as Rki (still in the order of the origin mark), and the n images can get n distances, and VRki represents the variance of nRki, then
(5)Wk=∑i=1nVRki−1.Order, available formula is as follows :ax=scosθ,ay=ssinθ(6)X2−Y2W0Y2X20WZ0X2Y20Z−Y2X2axaymxmx=X1Y1C1C2.


### 3.2.3. Dimension Reduction of PCA

PCA is principal component analysis, assuming that the number of feature points in each sample in the training set isk; then, a 2k-dimensional space is needed to process these shape vectors. Generally speaking, when K is large, it is difficult to deal with such a high-dimensional vector. In order to reduce the difficulty of processing, it is necessary to select some important features to describe the shape, that is, dimension reduction. The purpose of dimension reduction is to greatly reduce the redundancy of information while retaining most effective information. The most commonly used dimension reduction method is PCA (principal component analysis) processing, and the steps are as follows.Calculate the average shape vector:(7)X¯=1n∑i=1nXi.Calculate the covariance matrix:(8)S=1n−1∑i=1nXi−X¯TXi−X¯.The eigenvalueλi and the corresponding eigenvector pi of the covariance matrix are found, and the eigenvalues are arranged in a decreasing order, and the eigenvector corresponding to the previous eigenvalue will have the main effect on the shape vector. We select the first t eigenvectors to satisfy the formula
(9)∑i=1tλi∑i=1mλi>k.Selecting feature vectors like this can keep the original shape and reduce the computational complexity. In this way, any shape vector participating in training in the training image can be expressed as a linear combination of feature vector sets, as shown in the formula:(10)Xi=X¯+PB.


### 3.2.4. Constructing Local Gray Model

Principal component analysis adjusts the shape as a whole. Next, the shape needs to be iterated to the real contour edge, and each iteration needs to find a new position for each feature point, while the local gray model can provide the direction of iterative development. Here, the concept of gray gradient is used to define local features for each feature point. For thei-th feature point, the creation process of local features is shown in Figure 3.Figure 3
Description diagram of local feature creation.In thej-th training image, m pixels on each side of the i-th feature point are selected on the normal line of the i-th feature point (here, the normal refers to the vertical line connecting the previous feature point of the feature point and the second feature point), thus forming a vector with 2m+1 component. As shown in the following equation, each pixel on this vector is first-order guided, and the result of the operation is called local texture gij. The i-th feature point on the other images in the training set does the same work, and all the local texture about the i-th feature point gi1,gi2,⋯⋯,gin.
(11)gij=gij1,gij2,⋯⋯,gij2m+1T.


### 3.2.5. Model Initialization Positioning

When the model covers the target image, it needs initialization positioning. The accuracy and efficiency of ASM algorithm depend on the accuracy of initialization positioning. If the initialization positioning is far away from the real contour edge, the target may not be searched. If the initial location is close to the real contour edge, it can not only increase the accuracy of matching but also improve the speed of search, so it is necessary to transform the average shape vector properly to get a reasonable initial position.The standardized model uses initialization and positioning methods, and generally uses the maximum threshold method to binarize the segmented image shown in Figure4. The average shape model is randomly overlaid on the binary image, but to ensure that the longitudinal lowest point of the average shape model should not be lower than the elbow joint space, then, the head and tail point vectors of the average shape model are connected to form a closed area, and the average shape model is rotated and translated to make the white pixels in the closed area account for the largest proportion of the total pixels in the area. At this time, it can be considered that the initialization and positioning of the model are completed. The operation of rotation is limited by angle, while the operation of translation is unlimited in left and right directions and limited in up and down directions.Figure 4
Maximum binary diagram.The initial location of the posterior edge of the superior tibial plane is similar to that of the inferior femoral edge, while the initial location of the anterior edge of the superior tibial plane relies on the posterior edge and iterates on the basis of the posterior edge.


### 3.2.6. Local Search Processing

After the preliminary model is established, the initial ASM model is covered on the image, and the next step is to find the new iterative position of the feature points by using the local gray model of the feature points. The schematic diagram of the search process is shown in Figure5.Figure 5
Schematic diagram of search process.The big mark points on the outer side of the graph are the initial feature points, and the small mark points on the lower side of the graph are the best matching points of the model. The normal direction of the feature points refers to the perpendicular line connecting the front and back feature points of the initial feature points, and the iteration of each feature point moves along the normal direction. The thick line segment on the left side of the normal is2m+1 pixels selected in the local gray model, and the thick line segment on the right side of the normal represents the search area when the feature points are iterated. The distance between the initial point and the best matching point is the distance dX that the initial feature points should move.A local feature is obtained by normalization, which contains2v−m+1 sublocal features in total. Because the normalized gray vector approximately obeys Gaussian distribution, Mahalanobis distance can be used to express the similarity between the new feature g of a feature point and its trained local feature. The calculation formula of Mahalanobis distance is as follows:
(12)Fs=g−gi¯Si−1g−gi¯T.The smaller the Mahalanobis distance, the greater the similarity, and the closer the boundary pointI in the model is to the boundary point I in the target, so the best position of the boundary point can be determined by this method. Then, calculate the Mahalanobis distance between these sublocal features and the local features of the corresponding feature points in the gray model.When the Mahalanobis distance is minimized, the center point of the corresponding sublocal feature is the new position of the current feature point, and all feature points undergo the same steps, so each point will have a displacement, and these displacements will be arranged into vectors, such as the formula:(13)DX=DX1,DX2,DX3,⋯⋯,DXk.After all feature points are updated iteratively, a new curve is formed.


### 3.2.7. Active Shape Model Parameter Update

After the initialization of the shape model and the local search strategy are established, the shape can be updated. Assume that the mapping function from the average shape vector to the initial model is shown by the following formula:(14)X=Ms,θx+m.The obtained model is further similar to the target model, and further matching is carried out by iterative method. The iterative updating process of parameters is as follows:(15)mx⟶mx+dmx,my⟶my+dmy,θ⟶θ+dθ,s⟶s+ds,B⟶B+dB.As the iteration progresses, when the difference between the two shape models is less than a certain threshold or reaches a given number of iterations, the algorithm ends, stops searching, and thinks that the model converges, and the whole search and matching process ends, and the feature points at this time are the contour edges to be searched. Figure6 shows the detection results of the lower edge of femur after different iterations, which are 3 iterations, 15 iterations, and 40 iterations from left to right.Figure 6
Detection results of iterative process.


## 3.2.1. Sample Selection and Feature Labeling

In the segmentation method based on active shape model, the construction of initial model is very important, and to build a representative initial model, it is necessary to select appropriate samples, so the selection of samples is the foundation and an important step. First, selectn sample images as training samples. The selection principle should follow diversity and representativeness; for example, the orientation of elbow joint is oblique to the left and oblique to the right and vertical. The position of elbow joint is left, right, and middle, and the size of elbow joint is also slightly different. Various diseases of elbow joint lead to the change of elbow joint space and edge, all of which should be fully considered to avoid the limitation of sample selection. After selecting samples, we need to label the feature points of these samples and manually mark K key feature points of the target object in N samples, and the selection of feature points is based on the contour edge. Because the edge detected by Canny edge detection is the finest and the integrity is relatively good, the results of Canny edge detection are used to assist the selection of feature points. Taking the sampling and selection of feature points at the lower edge of femur as an example, sampling should not only ensure that the number of sampling points can describe the shape of the contour but also avoid the redundancy of information caused by too many sampling points. Generally, the points with high curvature and T-shaped joints are selected first, and the other feature points are arranged among the above points at reasonable intervals in order. As shown in Figure 2, the selection process of feature points of elbow joint femoral edge is described.Figure 2
Sampling schematic diagram of femoral edge of elbow joint.In Figure1, the larger dots are typical feature points, and the smaller dots are feature points added equidistantly between typical feature points to depict a complete shape. After selecting feature points, the coordinates of K feature points in an image are sequentially formed into a shape vector such as formula
(2)Xi=xi1,xi1,xi2,xi2,⋯⋯,xik,xik,where xik,xik is the coordinates of the k-th feature point in the i-th image.Thus, the length of each shape vector is 2 K, andn samples can form n shape vectors, thus obtaining a sample set J=X1,X2,⋯⋯Xn.


## 3.2.2. Sample Shape Alignment and Matching

The shapes of elbow joints are quite similar. For example, there will be an upward groove in the middle of the lower end of the femur, and the medial femur is usually larger than the lateral femur and so on. However, due to the slight inclination angle and orientation change of femur in training samples, in order to eliminate the interference of nonshape information caused by the standing posture and bone size of the subject, it is necessary to alignn shape vectors, which is called aligning shape vectors. The method adopted is Procrustes alignment. Generally speaking, general alignment is an affine transformation between shapes sought by using the idea of least square method. The basic steps are as follows:
(1)
After aligning to the origin, the shape vectors of the target objects in all training images are aligned to the shape vectors of the target objects in the first image, and the first shape vector should be normalized to the unit size(2)
Calculating the average valueX¯=1/n∑i=1nXi of the shape vector of the target object(3)
Aligning all the shape vectors generated in the first step toX(4)
Repeat steps 2 and 3 until convergence (the difference between two adjacent times is within a given threshold)E the expression to get the minimum value:
(3)E=Xi−Fs,θXj−mTWXI−Fs,θXj−m.BecauseFS,O stands for rotation and contraction transformation, we can deduce the formula
(4)Fs,θxkyk=scosθxk−ssinθykscosθxk−ssinθyk.The diagonal element of the weight diagonal matrixW is the weight Wk of each marker point; it represents the stability of the marking point. Taking the femoral edge as an example, some characteristic points are relatively stable, such as the groove vertex in the middle of the lower femur and the turning point on the left and right sides. Improving the weight of these points can effectively improve the accuracy of matching. Mark the distance of the k-th feature point and the i-th feature point in an image as Rki (still in the order of the origin mark), and the n images can get n distances, and VRki represents the variance of nRki, then
(5)Wk=∑i=1nVRki−1.Order, available formula is as follows :ax=scosθ,ay=ssinθ(6)X2−Y2W0Y2X20WZ0X2Y20Z−Y2X2axaymxmx=X1Y1C1C2.


## 3.2.3. Dimension Reduction of PCA

PCA is principal component analysis, assuming that the number of feature points in each sample in the training set isk; then, a 2k-dimensional space is needed to process these shape vectors. Generally speaking, when K is large, it is difficult to deal with such a high-dimensional vector. In order to reduce the difficulty of processing, it is necessary to select some important features to describe the shape, that is, dimension reduction. The purpose of dimension reduction is to greatly reduce the redundancy of information while retaining most effective information. The most commonly used dimension reduction method is PCA (principal component analysis) processing, and the steps are as follows.Calculate the average shape vector:(7)X¯=1n∑i=1nXi.Calculate the covariance matrix:(8)S=1n−1∑i=1nXi−X¯TXi−X¯.The eigenvalueλi and the corresponding eigenvector pi of the covariance matrix are found, and the eigenvalues are arranged in a decreasing order, and the eigenvector corresponding to the previous eigenvalue will have the main effect on the shape vector. We select the first t eigenvectors to satisfy the formula
(9)∑i=1tλi∑i=1mλi>k.Selecting feature vectors like this can keep the original shape and reduce the computational complexity. In this way, any shape vector participating in training in the training image can be expressed as a linear combination of feature vector sets, as shown in the formula:(10)Xi=X¯+PB.


## 3.2.4. Constructing Local Gray Model

Principal component analysis adjusts the shape as a whole. Next, the shape needs to be iterated to the real contour edge, and each iteration needs to find a new position for each feature point, while the local gray model can provide the direction of iterative development. Here, the concept of gray gradient is used to define local features for each feature point. For thei-th feature point, the creation process of local features is shown in Figure 3.Figure 3
Description diagram of local feature creation.In thej-th training image, m pixels on each side of the i-th feature point are selected on the normal line of the i-th feature point (here, the normal refers to the vertical line connecting the previous feature point of the feature point and the second feature point), thus forming a vector with 2m+1 component. As shown in the following equation, each pixel on this vector is first-order guided, and the result of the operation is called local texture gij. The i-th feature point on the other images in the training set does the same work, and all the local texture about the i-th feature point gi1,gi2,⋯⋯,gin.
(11)gij=gij1,gij2,⋯⋯,gij2m+1T.


## 3.2.5. Model Initialization Positioning

When the model covers the target image, it needs initialization positioning. The accuracy and efficiency of ASM algorithm depend on the accuracy of initialization positioning. If the initialization positioning is far away from the real contour edge, the target may not be searched. If the initial location is close to the real contour edge, it can not only increase the accuracy of matching but also improve the speed of search, so it is necessary to transform the average shape vector properly to get a reasonable initial position.The standardized model uses initialization and positioning methods, and generally uses the maximum threshold method to binarize the segmented image shown in Figure4. The average shape model is randomly overlaid on the binary image, but to ensure that the longitudinal lowest point of the average shape model should not be lower than the elbow joint space, then, the head and tail point vectors of the average shape model are connected to form a closed area, and the average shape model is rotated and translated to make the white pixels in the closed area account for the largest proportion of the total pixels in the area. At this time, it can be considered that the initialization and positioning of the model are completed. The operation of rotation is limited by angle, while the operation of translation is unlimited in left and right directions and limited in up and down directions.Figure 4
Maximum binary diagram.The initial location of the posterior edge of the superior tibial plane is similar to that of the inferior femoral edge, while the initial location of the anterior edge of the superior tibial plane relies on the posterior edge and iterates on the basis of the posterior edge.


## 3.2.6. Local Search Processing

After the preliminary model is established, the initial ASM model is covered on the image, and the next step is to find the new iterative position of the feature points by using the local gray model of the feature points. The schematic diagram of the search process is shown in Figure5.Figure 5
Schematic diagram of search process.The big mark points on the outer side of the graph are the initial feature points, and the small mark points on the lower side of the graph are the best matching points of the model. The normal direction of the feature points refers to the perpendicular line connecting the front and back feature points of the initial feature points, and the iteration of each feature point moves along the normal direction. The thick line segment on the left side of the normal is2m+1 pixels selected in the local gray model, and the thick line segment on the right side of the normal represents the search area when the feature points are iterated. The distance between the initial point and the best matching point is the distance dX that the initial feature points should move.A local feature is obtained by normalization, which contains2v−m+1 sublocal features in total. Because the normalized gray vector approximately obeys Gaussian distribution, Mahalanobis distance can be used to express the similarity between the new feature g of a feature point and its trained local feature. The calculation formula of Mahalanobis distance is as follows:
(12)Fs=g−gi¯Si−1g−gi¯T.The smaller the Mahalanobis distance, the greater the similarity, and the closer the boundary pointI in the model is to the boundary point I in the target, so the best position of the boundary point can be determined by this method. Then, calculate the Mahalanobis distance between these sublocal features and the local features of the corresponding feature points in the gray model.When the Mahalanobis distance is minimized, the center point of the corresponding sublocal feature is the new position of the current feature point, and all feature points undergo the same steps, so each point will have a displacement, and these displacements will be arranged into vectors, such as the formula:(13)DX=DX1,DX2,DX3,⋯⋯,DXk.After all feature points are updated iteratively, a new curve is formed.


## 3.2.7. Active Shape Model Parameter Update

After the initialization of the shape model and the local search strategy are established, the shape can be updated. Assume that the mapping function from the average shape vector to the initial model is shown by the following formula:(14)X=Ms,θx+m.The obtained model is further similar to the target model, and further matching is carried out by iterative method. The iterative updating process of parameters is as follows:(15)mx⟶mx+dmx,my⟶my+dmy,θ⟶θ+dθ,s⟶s+ds,B⟶B+dB.As the iteration progresses, when the difference between the two shape models is less than a certain threshold or reaches a given number of iterations, the algorithm ends, stops searching, and thinks that the model converges, and the whole search and matching process ends, and the feature points at this time are the contour edges to be searched. Figure6 shows the detection results of the lower edge of femur after different iterations, which are 3 iterations, 15 iterations, and 40 iterations from left to right.Figure 6
Detection results of iterative process.


## 4. Automatic Diagnosis Experiment of Elbow Arthritis Based on Edge Algorithm

### 4.1. Determination of Diagnostic Index

Clinically, some of these indicators are quantitative, such as elbow joint space distance and related angle. Some are nonquantitative and vague, such as osteophyte amount and sclerosis degree. They are all directly related to the symptoms of elbow arthritis. For example, narrowing the distance between elbow joints will aggravate the pain of patients when sitting, standing, and walking. Meniscus calcification can cause gout symptoms in patients. The changes of the upper femoral angle, the lower tibial angle, and the tibial-femoral angle will make the joints turn inside (outside) and so on. If we can automatically measure the quantitative indexes and transform the nonquantitative indexes into statistical quantitative analysis with the help of computers, the diagnosis efficiency will be greatly improved, and the workload of doctors will be reduced. In this chapter, the elbow joint space distance is measured automatically, and the automatic diagnosis results of these three indexes are compared with the results recorded in the reference documents in the image set.


### 4.2. Automatic Measurement of Elbow Joint Correlation Angle

In the X-ray film of elbow joint, the calculation of elbow joint-related angles including inferior femoral angle and superior tibial angle all depends on femoral physiological axis and tibial physiological axis. The femoral physiological axis is the line between the apex of the interfemoral fossa and the midpoint of the femoral axis, and the tibial physiological axis is the line between the intertibial fossa and the midpoint of the tibial axis. The inferior femoral angle is the lateral angle between the tangent line of the inferior femoral edge of elbow joint and the physiological axis of femur, which is generally between 75 and 85. If it is less than 75, it is defined as cubitus valgus. The superior angle of tibia is the lateral angle between tibial plateau and tibial physiological axis. It is generally between 85 and 95, and if it is greater than 95, it is defined as cubitus varus. As shown in Figures7(a) and 7(b), there are schematic diagrams of inferior femoral angle and superior tibial angle.Figure 7
Schematic diagram of inferior femoral angle and superior tibial angle.
(a)(b)Figures8 and 9 show the comparison between the automatic measurement results of the lower femoral angle, the upper tibial angle, and the femoral tibial angle of the right hand of 20 images in the image set and the manual measurement results in the reference document. “Solid blue six-pointed star” stands for automatic measurements, and “green circle” stands for manual measurements in reference documentation, in deg. It can be seen from the results that the automatic measurement results are in good agreement with the manual measurement results.Figure 8
Comparison of automatic and manual measurement of inferior femoral angle.Figure 9
Comparison of automatic and manual measurement of upper femoral angle.


## 4.1. Determination of Diagnostic Index

Clinically, some of these indicators are quantitative, such as elbow joint space distance and related angle. Some are nonquantitative and vague, such as osteophyte amount and sclerosis degree. They are all directly related to the symptoms of elbow arthritis. For example, narrowing the distance between elbow joints will aggravate the pain of patients when sitting, standing, and walking. Meniscus calcification can cause gout symptoms in patients. The changes of the upper femoral angle, the lower tibial angle, and the tibial-femoral angle will make the joints turn inside (outside) and so on. If we can automatically measure the quantitative indexes and transform the nonquantitative indexes into statistical quantitative analysis with the help of computers, the diagnosis efficiency will be greatly improved, and the workload of doctors will be reduced. In this chapter, the elbow joint space distance is measured automatically, and the automatic diagnosis results of these three indexes are compared with the results recorded in the reference documents in the image set.


## 4.2. Automatic Measurement of Elbow Joint Correlation Angle

In the X-ray film of elbow joint, the calculation of elbow joint-related angles including inferior femoral angle and superior tibial angle all depends on femoral physiological axis and tibial physiological axis. The femoral physiological axis is the line between the apex of the interfemoral fossa and the midpoint of the femoral axis, and the tibial physiological axis is the line between the intertibial fossa and the midpoint of the tibial axis. The inferior femoral angle is the lateral angle between the tangent line of the inferior femoral edge of elbow joint and the physiological axis of femur, which is generally between 75 and 85. If it is less than 75, it is defined as cubitus valgus. The superior angle of tibia is the lateral angle between tibial plateau and tibial physiological axis. It is generally between 85 and 95, and if it is greater than 95, it is defined as cubitus varus. As shown in Figures7(a) and 7(b), there are schematic diagrams of inferior femoral angle and superior tibial angle.Figure 7
Schematic diagram of inferior femoral angle and superior tibial angle.
(a)(b)Figures8 and 9 show the comparison between the automatic measurement results of the lower femoral angle, the upper tibial angle, and the femoral tibial angle of the right hand of 20 images in the image set and the manual measurement results in the reference document. “Solid blue six-pointed star” stands for automatic measurements, and “green circle” stands for manual measurements in reference documentation, in deg. It can be seen from the results that the automatic measurement results are in good agreement with the manual measurement results.Figure 8
Comparison of automatic and manual measurement of inferior femoral angle.Figure 9
Comparison of automatic and manual measurement of upper femoral angle.


## 5. Conclusion

In recent years, with the rapid development of imaging, various modern medical imaging technologies are widely used in the diagnosis of elbow osteoarthritis, such as X-ray, ultrasound, CT, and MRI. These imaging technologies have their own fields of expertise. Sometimes, it is necessary to combine various imaging technologies to analyze different images of the same diseased site for the diagnosis of diseases. In this paper, the automatic diagnosis of elbow osteoarthritis is completed based on X-ray image information. Firstly, the region of interest of elbow joint is extracted, a small rectangular region where the elbow joint is located is extracted from the whole X-ray image, and then, the contour extraction of hard bone edge in two-dimensional image is realized based on active shape model. Complete the automatic calculation of elbow joint related angle; the edge sclerosis of hard bone of elbow joint was quantitatively analyzed. In this chapter, two indexes, elbow joint space distance and elbow joint related angle, are selected for automatic measurement and quantitative analysis. In the aspect of elbow joint gap distance measurement, firstly, the elbow joint is corrected by rotation to make the automatic measurement of gap distance conform to the clinical medical measurement criteria, and then, the automatic measurement results are compared with the manual measurement results in the reference documents, showing a good agreement.The accuracy of elbow bone contour edge extraction needs to be improved, and the search rules in active shape model are still being tried to improve the accuracy of the algorithm from two aspects: the original gray model improvement and the subsequent search criteria improvement. The statistical description of various nonquantitative indicators of elbow arthritis needs to be improved. In the future, the quantitative analysis of all nonquantitative indicators should be completed to promote the comprehensiveness of automatic diagnosis of elbow arthritis.



---

*Source: 2199262-2022-10-11.xml*