Correlation of Subchondral Bone Density and Structure from Plain Radiographs with Micro Computed Tomography Ex Vivo

Osteoarthritis causes changes in the subchondral bone structure and composition. Plain radiography is a cheap, fast, and widely available imaging method. Bone tissue can be well seen from plain radiograph, which however is only a 2D projection of the actual 3D structure. Therefore, the aim was to investigate the relationship between bone density- and structure-related parameters from 2D plain radiograph and 3D bone parameters assessed from micro computed tomography (µCT) ex vivo. Right tibiae from eleven cadavers without any diagnosed joint disease were imaged using radiography and with µCT. Bone density- and structure-related parameters were calculated from four different locations from the radiographs of proximal tibia and compared with the volumetric bone microarchitecture from the corresponding regions. Bone density from the plain radiograph was significantly related with the bone volume fraction (r = 0.86; n = 44; p < 0.01). Mean homogeneity index for orientation of local binary patterns (HIangle,mean) and fractal dimension of vertical structures (FDVer) were related (p < 0.01) with connectivity density (HIangle,mean: r = −0.73, FDVer: r = 0.69) and trabecular separation (HIangle,mean: r = 0.73, FDVer: r = −0.70) when all ROIs were pooled together (n = 44). Bone density and structure in tibia from standard clinically available 2D radiographs are significantly correlated with true 3D microstructure of bone. Electronic supplementary material The online version of this article (doi:10.1007/s10439-015-1452-y) contains supplementary material, which is available to authorized users.


INTRODUCTION
Plain radiography is a cheap, fast, and widely available imaging method. Especially bone tissue can be well seen from plain radiographs which are significantly contributing the diagnostics of diseases that affect bone density and structure. However, as the plain 2-dimensional (2D) radiograph is a projection (summation) through the actual 3-dimensional (3D) structure, the obtained structural information is always limited compared to true 3D imaging modalities, e.g., computed tomography (CT) and magnetic resonance imaging.
Osteoarthritis (OA) causes changes in the articular cartilage and subchondral bone. Typical OA changes in the subchondral bone include bone sclerosis (thickening), osteophytes, and bone cysts. 6 Diagnosis of OA is routinely based on clinical examination and changes on plain radiographs. Typically, severity of OA is evaluated from radiographs using Kellgren-Lawrence grading scale. 14 However, Kellgren-Lawrence grading is based on subjective evaluation, it is semi-quantitative, and its inter-rater and intra-rater reliability varies from moderate to substantial. 11,35,36 To exploit all available information from the 2D radiographs, there is a need for development and use of quantitative and user-independent image analysis algorithms.
Previously, quantitative evaluation of OA changes from knee radiographs has included measurement of joint space width 5,34 as well as estimation of bone density 15,23,42 and structure. 4,17,[24][25][26][27]40 It is known that image acquisition parameters and post-processing algorithms affect the bone density evaluation from radiographs. 15 To overcome this issue, calibration of the grayscale values in an image using an aluminum step wedge have been proposed. 15,23,42 Texture analysis of bone is not as dependent on the imaging conditions as the direct evaluation of grayscale values. Bone structural analysis has been performed on plain radiographs using many different algorithms. For instance, the progression of OA has been assessed from digital knee radiography using signature dissimilarity measure method 40 and fractal signature analysis (FSA). 17 Fractal-based algorithms have also been applied to macro-radiographs 4,[24][25][26][27] and to standard film radiographs from OA knees. 16,31,41 Although macro-radiographs have a better spatial resolution, changes in bone structure can also be detected in standard radiographs. 25 Recently, local binary patterns (LBP) method for evaluation of bone structure from plain knee radiographs has also been introduced. 12 Basic LBP methods are computationally efficient and insensitive to monotonic grayscale variations. 28 Although differences in bone structure between subjects with different stages of OA and controls using LBP-based entropy have been reported, 12 the LBP-based methods have not been validated against true 3D microarchitecture of bone yet.
Previous studies with bone samples from human cadavers have shown that textural parameters from 2D high-resolution radiographs correlate significantly with 3D trabecular bone parameters. 19,29,33,37 However, these studies used small specimens harvested from human femur, not the entire bone, and cortical bone was removed from the 3D analyses. Furthermore, image analysis algorithms in these studies were mainly developed for quantification of the osteoporosis-related changes, not specifically for OA-related changes. For example, the fractal parameter (H mean ) that has been used in some osteoporosis-related studies is the average of all possible directions 19,29 while FSA provides fractal signatures in the horizontal and vertical directions.
The degree of relationship between 2D image texture parameters from standard clinical radiographs and 3D micro-CT (lCT) bone parameters from full thickness human tibia is unknown. Therefore, the aim of the current study is to investigate the correlation between bone density-and structure-related parameters from 2D plain radiograph and 3D bone parameters assessed from the lCT scans. We hypothesize that bone density and structure assessed from a 2D image are significantly correlated with the 3D microstructure of bone, possibly providing more sensitive diagnosis of bony changes in OA from conventional radiography.

Material
Tibial bones from right legs without soft tissue from eleven human cadavers (29-77 years of age) with no history of joint diseases were included in this study. 18 The cadaver knees were earlier obtained from Jy-va¨skyla¨Central Hospital, Jyva¨skyla¨, Finland, as approved by the national authority (National Authority for Medicolegal Affairs, Helsinki, Finland, Permission 1781/32/200/01).

Radiography
All bones were imaged using digital radiography (Ysio, Siemens Healthcare, Erlangen, Germany) with constant imaging parameters (63 kVp, 6 mAs, pixel size: 139 9 139 lm 2 , source-detector distance: 151 cm). Subsequently, the bones were immersed into a water bath (radius of the round plastic container: 6 cm) to simulate the effect of soft tissue and imaging was repeated using the aforementioned imaging settings.

Micro Computed Tomography (lCT)
After the radiography, the bones were cut into halves and both the medial and lateral condyles were imaged with lCT scanner (SkyScan 1176, Bruker Mi-croCT, Kontich, Belgium, 80 kV, 300 lA, 445 ms exposure, 2 frames averaged, isotropic voxel size of 17.4 lm, 0.04 mm copper + 0.5 mm aluminum filter) separately. Cutting of bones was required to fit the bones into the lCT scanner. To align the lCT slice stacks with the plain radiographs, the lCT slices were manually re-oriented by comparing the 2D coronal projection of the slice stack with the corresponding plain radiograph from the same bone.
For the lCT data, the regions of interest (ROIs) were placed in the 2D coronal projection image from the lCT slices (see: selection of regions-of-interest). After that, ROIs were extracted from every slice of the lCT stack separately. These lCT data, i.e. volumes of interest (VOIs), were then evaluated using SkyScan CTAn software (Bruker MicroCT, Kontich, Belgium). Before calculating the conventional 3D parameters, 3D median filtering (radius 2) and global thresholding (8bit grayscale value: 95) was applied to extract bone tissue from background. The calculated 3D parameters included bone volume fraction (BV/TV, the ratio of 3D total bone volume to total volume of VOI, in %), average trabecular thickness (Tb.Th, in lm), trabecular separation (Tb.Sp, mean thickness of the non-bone areas, in lm), trabecular number (Tb.N, average number of trabeculae per unit length, in 1/mm), structure model index (SMI, the relative prevalence of rods and plates), and connectivity density (Conn.Dn, the degree of connectivity of trabeculae normalized by total volume of VOI, 1/mm 3 ). The definition and calculation of these bone 3D parameters has been thoroughly described in the paper by Bouxsein et al. 3 Furthermore, to simulate plain radiography, all binarized lCT slices were summed together to construct a 2D coronal projection image from the 3D lCT data (Fig. 1). This high-resolution 2D projection image was analyzed using the same algorithms as used for plain radiographs to evaluate bone density and structure.

Selection of Regions of Interest
Four rectangle-shaped ROIs were extracted from the tibia (Fig. 1). Two ROIs (size: 344 9 803 pixels in lCT, 43 9 100 pixels in plain radiographs) were placed into the subchondral bone in the center of the medial and lateral condyles of tibia and two ROIs (803 9 803 pixels in lCT, 100 9 100 pixels in plain radiographs) immediately below the dense subchondral bone in the trabecular bone. Trabecular bone ROIs were aligned horizontally with subchondral bone ROIs. Anatomical landmarks for the ROIs were tibial spine, subchondral bone plate, and outer borders of the proximal tibia. A custom-made MATLAB software (version R2014b, The MathWorks, Inc., Natick, MA, USA) was used for the manual placement of the ROIs.

Evaluation of Bone Density from Radiographs
To evaluate bone density from the plain radiographs, mean grayscale value of the ROI (=GV) and aluminum step wedge thickness that corresponds to the measured GV (=GV mmAl ) were determined. The step wedge thickness was calculated using linear interpolation between grayscale values of consecutive steps in the wedge. The step wedge was present in all images.

Texture Analysis of Bone
Before bone texture analysis, radiographs were median filtered (3 9 3 pixels) to remove high frequency noise from the images. Bone structure was evaluated using Laplacian-based methods, 12,39 LBP-based methods, 12,38 and using FSA. 20,21 Laplacian-Based Analysis Laplacian-based image was constructed as previously described. 12 The Laplacian-based method enhances the appearance of bone trabeculae and quantifies the variation in the grayscale values of the Laplacian-based image. Laplacians were calculated in the horizontal and vertical directions and summed into one matrix. Subsequently, the unprocessed ROI was multiplied with square root of the Laplacian matrix to enhance the bone and grayscale values were expanded to full dynamic range to obtain the final Laplacian-based image. To measure the randomness of the grayscale values in the Laplacian-based image, entropy of the image (E Lap ) was calculated using the following equation: where P i contains the normalized count of the grayscale value i occurring in the image. If E Lap = 0, all pixel values in the Laplacian-based image are the same, whereas higher values indicate higher variation in the pixel values of the image.

Local Binary Patterns (LBP)-Based Methods
To measure randomness of local patterns and variation in the orientation of adjacent local patterns, LBP-based methods were modified from the methods developed initially for lCT data. 38 The LBP value of a studied pixel is assessed from the grayscale level of its surrounding, while ignoring the differences in magnitudes. In our method, the image was initially divided into bone and non-bone regions by determining a local threshold for every pixel in the image using Otsu method 30 with 9 9 9 pixels (36 9 36 for 2D projection image from lCT) window size. Next, LBP operator (8neighborhood on a circle with a radius of 1) was applied in the bone regions and in the non-bone regions next to bone, i.e. in the bone edge (pixel was considered as an edge pixel if at least one of the 8 neighbors of the center pixel was a bone pixel). To reduce the amount of irrelevant patterns, grouping of patterns was performed by determining the main orientation and the number of valid neighbors (i.e. markers) for each pattern. The main orientation angle was calculated using principal component analysis. The angle (0°, 45°, 90°, and 135°) was calculated only for the patterns that consisted of 2-5 consecutive markers, otherwise the pattern was assigned as nonuniform. Furthermore, to measure the randomness of the patterns occurring in the image, entropy of the grouped patterns (E LBP ) was determined using the Eq. (1). If E LBP = 0, there is only single pattern occurring in the image. Finally, the homogeneity index for the orientation of the valid patterns (HI angle ) was derived from the co-occurrence matrix of the angles. Co-occurrence matrices were calculated in 0°, 45°, 90°, and 135°directions with one pixel distance. The nonuniform and non-bone area was excluded from the cooccurrence matrices. HI at horizontal (HI Hor angle ) and vertical (HI Ver angle ) directions and mean HI (HI angle,mean ) of the four possible directions were used in the analyses. If all adjacent patterns have similar orientation, HI angle is one, while a large variation in the orientation of local patterns results a low HI angle value.

Fractal Signature Analysis (FSA)
To estimate fractal dimension, related to complexity and roughness of an image, FSA method was used. 20,21 The method produces the fractal signatures in the horizontal and vertical directions at individual scales. In brief, the original 3 9 3 median filtered image was dilated and eroded in horizontal and vertical directions with a rod-shaped one-pixel wide structuring element. The volume, V, between dilated and eroded images was then calculated. Calculations were repeated by varying the element length r from 2 to 4 pixels. The surface area, A(r), was obtained from the Eq. (2): After that a log-log plot was constructed by plotting log of A(r) against log of r. Finally, the fractal dimension was estimated using a regression line to points between 2 and 4 (between 2 and 32 for the 2D projection image from lCT). When the structuring element is pointing in the horizontal direction, fractal dimension of vertical structures (FD Ver ) is produced and vice versa. 20 High fractal dimension values are associated with high complexity of the image, whereas low complexity results in low fractal dimension values.

Statistical Analysis
To evaluate relationships between different parameters, Pearson's correlation analysis (together with 95% confidence intervals 1 ) was applied using IBM SPSS Statistic for Windows (Version 22.0, IBM Corp., Armonk, NY, USA). Table 1 shows the mean and standard deviation (SD) values for the 3D lCT parameters. The mean and SD of the 3D lCT parameters in each ROI separately are shown in the Supplementary Table 6.

Bone Density
Correlation between bone density evaluated from the plain radiograph using GV mmAl parameter and BV/TV was strong and statistically significant ( Table 2; Fig. 2). Correlation remained significant also when the bone was immersed in the water bath during radiography ( Table 2). High correlations were observed even when the ROIs were considered separately (Supplementary Table 7).
Furthermore, a strong correlation between GV from 2D projection image from lCT data and BV/TV was obtained ( Table 2 and Supplementary Table 7).

Bone Structure
Significant correlations between bone structural parameters from plain radiograph and 3D bone architectural parameters from lCT were obtained, yet, the degrees of correlations varied depending on the parameter and the direction on which the directional texture measures were calculated (Tables 3, 4, 5 and  Supplementary Tables 8, 9, 10, and 11). HI Hor angle and FD Ver were more strongly related with Conn.Dn and Tb.Sp than HI Ver angle and FD Hor in the trabecular bone ROIs (Table 5 and Supplementary Tables 10 and 11). In the Fig. 2, significant correlations between HI angle,mean and Tb.Sp and between FD Ver and Tb.Sp are shown.
Significant correlations between bone structure parameters assessed from both the original 3D lCT data and its 2D projection were also obtained (Tables 3, 4, 5 and Supplementary Tables 8, 9, 10, and 11).

DISCUSSION
Current results for full thickness tibial bones indicate that estimates for both the subchondral bone density and structure, evaluated from 2D plain radiographs, are significantly correlated with the corresponding 3D parameters from lCT.
The bone density estimated from the plain radiograph after grayscale calibration to the aluminum step wedge (GV mmAl ) was strongly related with the bone volume fraction from lCT. The correlation remained high even when the effect of soft tissue was simulated with the water bath. The water bath increased scattering and decreased the quality of the image. The mean grayscale value from the 2D coronal projection from binarized lCT slices correlated also strongly with the bone volume fraction. This finding is in line with a previous study that showed a strong correlation (r = 0.90) between 2D density estimate and bone volume fraction from 3D lCT data. 37 Therefore, based on the current results, estimation of bone volumetric density from 2D radiographs is feasible at least when the grayscale values corresponding to the bone falls into the range of the step wedge grayscale values. In this case, aluminum thickness corresponding to the mean grayscale value of bone can be linearly interpolated. However, if the grayscale values of the bone are outside the range of the step wedge (i.e. corresponding aluminum thickness needs to be extrapolated), the method may be less reliable since the detector response of the X-ray device is usually not perfectly linear.
From the texture parameters evaluated from the 2D radiograph, especially fractal dimension of vertical structures (FD Ver ) and mean homogeneity index for orientation of local patterns (HI angle,mean ) were significantly related with the connectivity density and trabecular separation in 3D. This finding for the fractal dimension is consistent with an earlier study. 22 Since the HI angle,mean parameter is the mean HI value from four different directions, it is less affected by the orientation of the image. It can be hypothesized that if the image would be oriented along the trabeculae, fractal dimensions or directional homogeneity indices would correlate even more strongly with thickness and separation of trabeculae. Our results support this hypothesis since the degree of correlation varied depending on which direction the fractal dimension or homogeneity indices of local patterns were calculated. For example, HI angle in the horizontal direction and FD Ver in trabecular bone area were significantly related with the TABLE 2. Pearson correlation coefficients (95% confidence interval) between bone densities evaluated from both plain radiographs and 2D lCT projection image and BV/TV.  trabecular separation whereas HI angle in the vertical direction and FD Hor were less related. This is because the trabeculae were aligned more vertically than horizontally in this data set and, thus, when calculating HI angle , there were less variation in the orientation of adjacent local patterns in vertical direction. However, in the previous texture analysis studies of a knee joint, the images have not been oriented along the main direction of trabeculae and therefore we decided to not orient the images either. It should be noted that when calculating FD Ver , the structuring element was actually pointing in the horizontal direction. 20 Therefore, FD Ver and HI angle in the horizontal direction are related although not directly measuring the same phenomenon.
The degrees of correlation between lCT parameters and entropies (Laplacian-based or entropy of grouped patterns) were variable. One explanation for the positive correlation between entropy of grouped patterns and connectivity density is that when the bone is highly connected, more different patterns are detected in the texture analysis and eventually the entropy of patterns is therefore higher. Laplacian-based entropy might be better suited for the analysis of femoral neck, where the orientation of the trabeculae is usually clear and the ROI can be easily aligned along the trabeculae. 39 In the current study, the Laplacians were calculated in the vertical and horizontal directions and summed together, which may have decreased the sensitivity of method for bone changes.
In our analyses, both the subchondral bone and trabecular bone VOIs contained both cortical and trabecular bone. These bone types are also superimposed in the plain radiograph and, therefore, we did not extract cortical bone layer from the lCT analyses. Furthermore, the cortical bone layer is very thin at the area from which the trabecular bone VOI was extracted. When cortical bone was removed from the lCT data, the results remained virtually the same (data not shown, tested with medial trabecular bone VOIs). This has also been confirmed in an earlier study in which bone volume fraction was calculated with and without cortical bone and a strong correlation was obtained (r = 0.73). 10 Our results show that radiographic texture analysis may serve as a complementary method in OA diagnostics since good correlations with 3D microarchitecture of bone were obtained. This conclusion is supported by the finding that the bone density-and structure-related parameters from radiographs correlated significantly with 3D lCT parameters that have been shown to change during OA. In general, osteoarthritic subchondral bone has higher bone volume fraction and trabeculae are thicker than in healthy bone. 2,4,[7][8][9]13 Structure model index has been reported to be lower in OA bone than in the control bone. 7,8 Some studies have reported higher connectivity and increased trabecular separation but fewer trabeculae in OA bone compared to controls. 8, 13 However, there are also studies suggesting that in OA bone, number of trabeculae is higher and they are closer together than in controls. 2,4,[7][8][9] This discordance is likely due to the difference in anatomical sites studied and to the different stages of OA in the samples.
The degrees of correlations between radiographic parameters and 3D microarchitecture of bone were not similar in each ROI in the current study. Several factors including differences between subchondral and trabecular bone regions as well as between medial and lateral regions have influenced these variations. For      example, organization of trabecular network is not similar over these regions, which might have affected the bone structural parameters as discussed earlier.
Furthermore, the range in bone volume fraction was smallest in the medial trabecular bone and, therefore, the degree of correlation between the radiograph-based density and bone volume fraction might have been lower in this area. As the pooling of the ROIs might have artificially increased the significance of the correlation coefficients, the correlation coefficients for each ROI are separately shown in the Supplementary Material. The most significant limitation of the current study is that the donors did not have any diagnosed joint disease at the time of death and, thus, only limited variation in bone density and structure was presumed. However, as the age range of the cadavers was 29-77 years, it is highly probable that some of the bones from older cadavers have actually had some osteoarthritic tissue-level changes. This assumption is further supported by the histological Mankin scores available for the cartilage-bone samples drilled from the contralateral knee 32 : the Mankin scores varied from 1 to 9 (healthy = 0, severe OA = 14) for the samples from the medial and lateral tibial plateaus. On the other hand, the variability of the 3D parameter values obtained from the lCT data also suggests that the study sample was relatively heterogeneous. Consequently, actual variation in bone density and structure in this sample set has probably been higher than could be presumed for completely intact samples. Nevertheless, as the visual signs of OA were not specifically evaluated from these tibiae, this limitation remains and justifies future studies with a larger sample set including both non-OA and OA subjects to further clarify the sensitivity of the methods reported here. As second limitation, our samples did not contain soft tissue that reduces quality of the radiograph and, therefore, generalization of the methods in vivo is partially restricted. However, the effect of soft tissue was still simulated by immersing the bones into a water bath during radiography.
In conclusion, estimates for the subchondral bone density and structure, evaluated from 2D plain radiographs, were significantly correlated with the corresponding 3D parameters from lCT. Therefore, evaluation of bone density and bone structure is feasible from the standard clinically available radiographs.

ELECTRONIC SUPPLEMENTARY MATERIAL
The online version of this article (doi:10.1007/s 10439-015-1452-y) contains supplementary material, which is available to authorized users. 8