Quantifying Subresolution 3D Morphology of Bone with Clinical Computed Tomography

The aim of this study was to quantify sub-resolution trabecular bone morphometrics, which are also related to osteoarthritis (OA), from clinical resolution cone beam computed tomography (CBCT). Samples (n = 53) were harvested from human tibiae (N = 4) and femora (N = 7). Grey-level co-occurrence matrix (GLCM) texture and histogram-based parameters were calculated from CBCT imaged trabecular bone data, and compared with the morphometric parameters quantified from micro-computed tomography. As a reference for OA severity, histological sections were subjected to OARSI histopathological grading. GLCM and histogram parameters were correlated to bone morphometrics and OARSI individually. Furthermore, a statistical model of combined GLCM/histogram parameters was generated to estimate the bone morphometrics. Several individual histogram and GLCM parameters had strong associations with various bone morphometrics (|r| > 0.7). The most prominent correlation was observed between the histogram mean and bone volume fraction (r = 0.907). The statistical model combining GLCM and histogram-parameters resulted in even better association with bone volume fraction determined from CBCT data (adjusted R2 change = 0.047). Histopathology showed mainly moderate associations with bone morphometrics (|r| > 0.4). In conclusion, we demonstrated that GLCM- and histogram-based parameters from CBCT imaged trabecular bone (ex vivo) are associated with sub-resolution morphometrics. Our results suggest that sub-resolution morphometrics can be estimated from clinical CBCT images, associations becoming even stronger when combining histogram and GLCM-based parameters. Electronic supplementary material The online version of this article (10.1007/s10439-019-02374-2) contains supplementary material, which is available to authorized users.


INTRODUCTION
Subchondral bone sclerosis, causing increases in bone volume fraction (BV/TV, ratio of bone volume and tissue volume), trabecular thickness (Tb.Th.), trabecular number (Tb.N.) and a decrease in trabecular separation (Tb.Sp.), has been associated with osteoarthritis (OA)-driven cartilage defects. 2,5,10 These alterations in subchondral bone are often linked to the later stages of OA. However, in early OA, contrary structural alterations in subchondral bone have been reported within several animal models. 1,9,16 Furthermore, an increase in subchondral bone resorption, due to abnormally high bone turnover, has also been observed to occur in progressive OA. 23 The contradictory results in subchondral bone alterations in animal studies, and limited evidence of early OA-related subchondral bone alterations in human tissue, supports the need for further research on early OA-induced alterations of human subchondral bone at the micro-and nanostructural level.
Computed tomography (CT) modalities provide spatial resolution starting from 100 nm with synchrotron radiation nano-CT up to a 0.2-0.5 mm with clinical CT, enabling the hierarchical imaging of subchondral bone from the sub-cellular level to the organ level. 21 Desktop micro-computed tomography (lCT) imaging has become the gold standard for quantification of bone morphology and microstructure in 3D. 24 When approaching to the spatial resolution of clinical CT, trabecular bone structure has been reported to be quantifiable up to 100 lm voxel size. 17 On the other hand, with dental CT-imaged bone, strong associations between the average of the grey-level values and microstructural trabecular bone morphometrics from lCT have been reported in the literature. 15 In order to quantify bone microarchitecture from clinical X-ray based imaging modalities, several texture analysis algorithms have been applied in osteoporosis and OA research applications. Notably, fractal analysis and run-length distribution texture analysis methods from X-ray radiographs, and variogram from dual-energy absorptiometry, have been reported to correlate with morphometric parameters of bone. 4,19,26 Textural analysis of 2D radiographs has been shown to be sensitive to OA-related structural alterations and changes in bone density of subchondral bone. 13,14,19 For volumetric data, grey-level co-occurrence matrix (GLCM) based texture analysis, first introduced by Haralick, 11,12 has been previously applied to medical CT imaging of bone. 8,20,22 Moreover, GLCM parameters have shown good correlations to the BV/TV and biomechanical properties of bone from lCT data when the volumetric data is re-projected in 2D. 30 GLCM-based textural analysis is a method to extract second order statistical features from grey-level images. In comparison to histogram based parameters (i.e. grey-level mean and standard deviation, and histogram skewness and kurtosis), GLCM is not related to the absolute value of the grey-levels but to the relationship of the pixel's (or in 3D voxel's) grey-level value to the grey-level values of its neighboring pixels/ voxels. Thus, calculating GLCM and histogram-based parameters are complementary to each other. In the case of X-ray radiography and clinical CT imaging of subchondral bone, where resolution is insufficient for calculating the true morphometrics of bone, the grey-level values can still be reflective of the microarchitecture.
As CT modalities are used to study bone morphometrics at several hierarchical length scales, histogram and GLCM-based analysis of the grey-level values are likely to produce further information from the subresolution features. This will provide new tools to researchers, to help them analyze their data beyond the resolution limits of the CT system. This study aims to quantify sub-resolution bone morphometrics, which are related to OA, from clinical cone-beam CT (CBCT) ex vivo. We utilize GLCM texture and histogrambased parameters of the CBCT-imaged subchondral bone with various OA severities, and compare them with the morphometric parameters quantified from lCT. We hypothesize that trabecular bone morphometrics (BV/TV, Tb.Th., Tb.N., etc.) quantified from lCT associate with the GLCM parameters and histogram parameters calculated from the clinical CBCT data.

MATERIALS AND METHODS
53 osteochondral cores (Ø = 4 mm) were extracted from total knee arthroplasty (TKA) patients and from cadavers without OA diagnosis under the approval of ethical committees of Northern Ostrobothnia Hospital District (Finland, Approval No. 78/2013) and Research Ethics Committee of Northern Savo Hospital District (Finland, Approval No. 58/2013 and 134/ 2015). Dental drill with a trephine blade was used for extracting the cores. From the 53 cores, 46 cores were extracted from the tibial plateaus of 2 TKA patients (females aged 68 and 73) and 2 cadavers (males aged 68 and 68). From cadavers, cores from both tibiae were included (see Table 1). Only one tibia per TKA patient was subjected to the core extraction, because sample harvesting from these patients is conducted from the remnant tibial pieces of the TKA operation. Cores were extracted from 8 manually determined locations on the tibial plateau (Fig. 1). The remaining 7 cores (out of 53) were extracted from the weight-bearing area of TKA patient femora (N =7, 5 females and 2 males, age range: 66-86) to ensure high variability in OA severity. After core extraction, samples were preserved in a 2 80°C freezer until thawed for imaging. Prior to imaging, samples were immersed in 4% saline buffered formaldehyde for fixation and imaged in the fixation media.
The CBCT imaging was conducted with the clinical extremity CBCT scanner (Planmed Verity, Planmed Inc., Helsinki, Finland) with 80 kV, 12 mA, 200 lm voxel size, 2 min scan time, 300 projections and 20 ms pulse time. Reconstruction was conducted with the manufacturer's own reconstruction software with ''standard'' reconstruction filter.
CBCT and lCT volumes were co-registered using rigid transformations with the Dataviewer-software (v.1.5.4, Bruker microCT). Each co-registered dataset was visually evaluated to confirm the successful coregistration. For the CBCT and lCT comparison, cylindrical volume of interest (VOI) of 3000 9 3000 9 Z lm 3 (where Z is the depth from the subchondral bone plate, range: 605-4435 lm) was selected from the trabecular bone. Z was restricted based on how much trabecular bone was left in biopsies after the surgeries. The VOI was visually evaluated to be sufficient for each sample so that the minor alterations to the edges of trabecular bone core due drilling would be excluded from the analysis. The bone masks were utilized for the morphometric analyses of lCT-imaged trabecular bone. Figure 2 presents the flowchart of the methods and analyses used.
CTAn software (v.1.17.7.2, Bruker microCT) was used to calculate the relevant 3D morphometrics of the bone from the lCT data. We considered that the following trabecular bone morphometrics primarily affect to the grey-level values of CBCT imaged subchondral trabecular bone: BV/TV, Tb.N., Tb.Th., Tb.Sp., and local bone surface complexity (fractal dimension, FD). All the aforementioned morphometrics were calculated within the VOI of the lCT-imaged trabecular bone data (Fig. 1).
The histogram and GLCM texture parameters were calculated from the CBCT data within the VOI described above. The histogram parameters (mean, standard deviation, skewness, kurtosis and image entropy) and GLCM texture parameters were calculated with Matlab (v.2017b, Mathworks, Natick, MA, USA).
To briefly summarize the calculation of the GLCM, the GLCM is a N g 9 N g matrix where N g is the number of quantized grey levels. The GLCM consists of elements P(i,j), representing the number of occurrences in grey-levels i and j within a certain window defined by the displacement d and angle h for 2D images (and for 3D matrix [h, u]). The GLCM is then normalized by calculating the second order statistical probability values p(i,j) with Eq. 1: In our study, 13 GLCMs from 13 individual angles of the 3D dataset (Supplementary Table 1) were calculated and summed to provide rotationally invariant texture features. The selection of parameters d and N g affect to the calculated texture parameters. 7 We selected the voxel displacement d to be 1, because in CBCT the morphometric features, to which we correlated with the texture parameters, were below our imaging resolution. For 8-bit grey-level images, the N g can obtain values up to 256, reducing N g from 256 can decrease the noise and calculation time of the analysis in the expense of diminution of the information. Previous studies 7 have noted that the classification accuracy can be retained when N g ‡ 24. In this study, N g was selected to be 256 for the GLCM generation from the CBCT, to address the strong impact of the partial volume effect occurring in CBCT imaging. As the partial volume effect sums quantifiable micromorphometrics in relatively large voxels (as obtained from CBCT), it is necessary to keep N g high to avoid additional loss of morphometric information, which occurs when grey-levels are grouped during the GLCM calculation.
From the available textural features, inverse difference moment (IDM), angular second moment (ASM), contrast, variance, entropy, and cluster shade were calculated from the summed GLCM. While several other parameters could be calculated from the GLCM, in this study, only parameters which have often been reported in literature were considered ( Table 2).
The equation for the image entropy of the histogram based parameters is the same as for the GLCM parameters is described in Table 2. In the case of image entropy the p(i,j) is replaced by normalized histogram counts. Thus, as GLCM entropy describes randomness in texture (within the window limited by the displacement values), image entropy describes randomness of the grey-level values.
After lCT and CBCT imaging, the osteochondral cores were subjected to histological analyses. Histological sections (3 lm thick) were stained with Safranin-O and imaged with a digital pathology slide scanner (940 magnification and 0.25 lm pixel size; Aperio AT2, Leica Biosystems, Wetzlar, Germany). Subsequently, severity of osteoarthritis was evaluated using the OARSI histopathological grading system 27 from three consecutive sections by two independent graders. In the case of disagreements, a consensus grade was used for each sample.
Pearson's correlation coefficients were calculated to evaluate the association of morphometric parameters calculated from lCT with individual histogram and GLCM texture parameters calculated from CBCT. As the OARSI grade is ordinal, correlations between the morphometric, GLCM and histogram parameters were evaluated using the Spearman's rank correlation coefficient. Furthermore, we tested whether multiple histogram and GLCM parameters could estimate trabecular structure better compared to individual histogram or GLCM parameters by using stepwise linear regression with morphometric parameters as dependent variables (separate models) and all histogram and GLCM parameters as independent variables. Multicollinearity in stepwise linear regression was taken into account by calculating the variance inflation factor (VIF) and excluding parameters with VIF > 5 29 from the models. Furthermore, homoscedasticity of the data and residual distribution were qualitatively evaluated from the P-P-plot of regression standardized residuals and from the residuals vs. fitted plot. These statistical tests for comparing structural information with GLCM, histogram parameters and OARSI scores were conducted for all the 53 samples (osteochondral cores).
To evaluate whether the location of core extraction affects the results, in addition to correlation coefficients for the whole dataset, the correlation coefficients for subgroups based on the origin (TKA patients and cadavers), compartmental location (medial tibial plateau and lateral tibial plateau, femora excluded), and areal location (central tibial plateau, anterior tibial plateau, posterior tibial plateau, and distal tibial plateau) are also reported (see Table 1 for subgroup sample sizes). The dependency of the location was tested also with the models generated in stepwise linear regression described above. The effect of core location was evaluated by including the subgroup information to the model to see whether it would have an effect on the correlation (change in adjusted R 2 ). The number of cores per subgroup is shown in Table 1. All the sta-

RESULTS
BV/TV from lCT had strong correlation (|r| > 0.7) with the histogram parameters (mean, standard deviation), and GLCM texture parameters (correlation, cluster shade, and variance) from CBCT (Table 3 and Fig. 3). Because BV/TV is highly associated with Tb.Th., Tb.Sp., Tb.N. and FD, it is not surprising that GLCM and histogram parameters had similar associations with these other morphometric parameters too. However, only few histogram (mean and standard deviation) and GLCM parameters (correlation and cluster shade) had strong or moderate correlations (|r| > 0.5) with all the morphometric parameters. BV/ TV and Tb.N. had the highest correlations with mean (r = 0.907 and r = 0.834, respectively), Tb.Th. with standard deviation (r = 0.676), Tb.Sp. with GLCM variance (r = 2 0.803), and FD with GLCM correlation (r = 0.759). Regarding associations with OA severity, OARSI grade had moderate correlations with mean (r > 0.612), standard deviation (r > 0.663), and with GLCM cluster shade (r > 0.612).
Subgrouping revealed that the histogram mean, GLCM cluster shade, GLCM variance and GLCM correlation had strong correlation (|r| > 0.7) with BV/ TV both when the subgroups were considered separately (Supplementary Tables 2-5) or pooled together. Since BV/TV has dependency with other morphometrics, we conducted the stepwise linear regression only with BV/TV. Stepwise linear regression (Table 4) revealed that the model with mean and GLCM IDM parameters yielded highest coefficient of determination (adjusted R 2 = 0.864) with BV/TV, and when compared with the best individual parameter histogram mean the adjusted R 2 was lower (adjusted R 2 = 0.818, adjusted R 2 change = 0.046). The model was also robust to the core extraction location as the addition of location did not improve the adjusted R 2 values (TKA/Cadaver grouping: adjusted R 2 change = 2 0.001; Compartmental grouping: adjusted R 2 change = 2 0.002; Areal grouping: adjusted R 2 change = 2 0.001).
With regard to subgrouping of the samples for cadavers and TKA patients, all the trabecular morphometric parameters did have weak to moderate correlations with the OARSI grade in the cadaver and TKA groups (Supplementary Table 3). From the morphometric parameters, in both groups, BV/TV yielded the highest correlation with the OARSI grade (cadaver group: r = 0.418, TKA group: r = 0.650). Subgrouping to medial and lateral tibial plateau groups (Supplementary Table 4), revealed that all the trabecular morphometric parameters in medial tibial plateau group had significant correlations ranging from weak to moderate (0.432 < |r| < 0.684) with the OARSI grade. BV/TV and Tb.Th. (r = 0.681 and r = 0.684, respectively) had the highest correlations with the OARSI grade in the medial subgroup. However, in the lateral tibial plateau only weak correlations of Tb.Sp. and Tb.N. with the OARSI grade were observed (r = 2 0.462 and 0.408, respectively). Compartmental subgrouping (Supplementary Table 5) showed no significant correlations between the trabecular bone morphometrics and the OARSI grade in distal and anterior groups. However, in central tibial plateau group, the BV/TV had moderate and significant correlation with the OARSI grade (r = 0.673), and in posterior tibial plateau group BV/TV, Tb.Th., and Tb.Sp. had similarly significant and moderate correlations with the OARSI grade (r = 0.655, r = 0.658, and r = 2 0.598, respectively).

DISCUSSION
In this study, we aimed to quantify sub-resolution subchondral bone morphometrics related to OA from clinical CBCT ex vivo. Our results demonstrate that sub-resolution features can be quantified, from clinical CBCT using 3D texture analysis by calculating histogram and GLCM parameters.
The primary finding of our study is that select parameters from histogram and GLCM analyses are In all equations l x , l y and r x , r y denote the mean and standard deviation of the row and column sums of the GLCM, respectively.

BIOMEDICAL ENGINEERING SOCIETY
associated with specific, micro-level morphometric structures of subchondral bone better than the others. As seen in Table 3 and Fig. 3, histogram mean, histogram standard deviation, GLCM cluster shade, GLCM correlation and GLCM variance seem to associate well with the various morphometric structures. However, our results also indicate that using only individual histogram or GLCM parameters to predict individual morphometric structures might not feasible. For example, BV/TV had strong correlation with all of the aforementioned histogram and GLCM parameters (r > 0.782). Even though histogram mean had the highest correlation with the BV/TV when all the cores were included, in some subgroups the strongest association was with other histogram/GLCM parameters (e.g. standard deviation in cadavers, medial tibial, anterior tibial, and posterior tibial plateau subgroups, and GLCM cluster shade in lateral tibial and distal tibial plateau subgroups). Since the relation of the histogram/GLCM parameters with morphometric parameters seem to vary depending on the subgrouping, the use of multiple histogram/GLCM parameters combined to predict the trabecular bone morphometrics is the most preferable option in the future.
The stepwise linear regression model combining histogram and GLCM parameters (histogram mean and GLCM IDM) was robust for the location (with all subgrouping criterias: adjusted R 2 change < 0) and predicted BV/TV better when compared to the individual parameters (histogram mean: adjusted R 2 = 0.818; histogram mean and IDM: adjusted R 2 = 0.864). As the correlation between individual his-togram mean and BV/TV is already strong, only a slight increase in coefficient of determination (adjusted R 2 change = 0.047) of combined histogram mean and GLCM IDM with BV/TV might seem redundant. However, it should be noted that this study was conducted with small ex vivo osteochondral cores. This set limitations as the imaging and reconstruction protocols are optimized for in vivo limbs. First of all, we are lacking the noise generated from the soft tissue and surrounding bone matrix which is present in vivo. Secondly, this study uses only 8-bit images as the imaging and reconstruction parameters did limit the full utilization of the 16-bit dynamic range for small osteochondral cores. The combination of the first and second order grey-level statistics could be advantageous in reducing the effect of noise in the analysis of the in vivo clinical CT data. This is true especially if 16bit dynamic range can be utilized.
In this study, the used voxel size was 200 lm. Due to the technical advancements in current CBCT systems, the high resolution scans of extremities with minimal radiation dose has been surfaced to clinical practice. 25 For example, the CBCT scanner used in this study has been reported to cause only 6 lSv effective dose (with 200 lm resolution) in a CBCT scan of the ankle, which compares to the effective dose of 4 conventional 2D radiographs from the same anatomical region. 18 The histogram/GLCM analysis might be best suited for clinical CBCT images of extremities, especially because the multidetector CT (MDCT) systems provide poorer resolution (~0.6mm voxel size) with higher radiation dose. Other approaches determining trabecular bone morphometrics from clinical resolution images have been proposed. Methods utilizing peripheral quantitative CT and lMRI can reach resolutions around 80-140 lm voxel size, and thus, direct estimates of bone strength could be calculated. 3,28 However, with those methods, sufficient resolution is required as they often require binarization steps. In our study, the resolution of CBCT 200 lm was not sufficient for binarization as the average trabecular structure thickness varied from 0.5-1 voxel side lengths. In this type of scenario, the partial volume effect prevents the differentiation of trabecular bone from the trabecular cavities and the conventional selection of the threshold values based on the known attenuation values of bone is not possible. Approaches using fuzzy skeletonization can partially overcome the problems of binarization and these methods have proven good results from MDCT scans with clinical resolution images (~200 lm voxel size). 6 Unfortunately, Chen et al. reported that the accuracy of this method decreases with VOI diameter smaller than 5.25 mm, 6 and in our study the core diameter restricted the VOI diameter to 3 mm. Therefore, comparison of those methods with GLCM/histogram analysis remains to be elucidated in future studies with larger samples.
The combined histogram/GLCM parameters yielded better association with the morphometrics than the individual histogram/GLCM parameters. As this shows promise for future studies, one must be aware of collinearity and autocorrelation when creating these types of multiparameter models from GLCM/histogram parameters. There are mathematical based collinearities and autocorrelations between the parameters, e.g. between mean and standard deviation where both are dependent on the grey-level average. Another example includes inverse relationships between some of the GLCM parameters, such as GLCM IDM and GLCM contrast, i.e., another value increases the other value tends to decrease. In this study we considered both positive and negative relations as an advantage in interpretation of our results. Often with multiple correlations, the simultaneous statistical interference might result in Type I errors, and consequently, p values require correction. Here the interpretation of the individual significant correlations of the histogram and GLCM parameters were always compared with the correlations of the variable with dependency (i.e. previously mentioned mean and b FIGURE 3. Scatter plots of the highest correlations with the histogram/GLCM parameters with the different trabecular bone morphometrics. Correlation coefficients for all parameters, including the ones presented in this figure, are presented in Table 3  standard deviation) before making conclusions about the significance of the parameter correlations. Thus, understanding of the parameters' descriptive properties in addition to the mathematical base is crucial when generating more complex multiparameter prediction models. The association of trabecular bone morphometric parameters with the histopathological OARSI grade were in line with our previous study. 10 With the increasing OARSI grade, which evaluates the progressive degeneration of articular cartilage, we found a similar decrease in trabecular bone separation and increases in trabecular bone volume fraction, trabecular thickness, trabecular number and fractal dimension. The main difference compared to our previous study is that the current associations (correlation coefficients) between the trabecular bone morphometrics and the OARSI grade are only weak to moderate, whilst in the previous study the correlations were strong (|r| > 0.7). This difference might be related to a different and smaller patient population here. Furthermore, the associations between the OARSI grade and structural parameters varied in different subgroups, which supports strong dependence on the sample location. Thus, more controlled core extraction location-wise would be preferable in future studies when investigating these kinds of associations. Interestingly, OARSI had stronger correlations with histogram mean, histogram standard deviation and GLCM cluster shade (r = 0.612, r = 0.664, and r = 0.612, respectively) while the highest morphometric parameter association with the OARSI grade yielded r > 0.584(BV/TV). This gives indication that the underlying image texture/grey level information could be used in classification of OA severity from clinical CT bone data. These associations between the histogram/GLCM parameters related to OA severity need to be investigated with a larger patient population. Still, a smaller patient population in this study does not restrict us from making general conclusions on the associations between local trabecular bone structure and histogram/GLCM parameters calculated from clinical CT images.
In conclusion, we demonstrated that GLCM texture and histogram-based parameters from trabecular bone are associated with sub-resolution morphometrics in CBCT ex vivo. Our results also suggest that sub-resolution morphometrics can be predicted from clinical CBCT images even more accurately when combining histogram and GLCM texture based parameters. These methods show great potential for deriving the microstructural subchondral bone alterations from clinical CBCT. However, further studies of the feasi-bility of the combined use of histogram and GLCM based parameters for in vivo are required.