A Digital Model to Simulate Effects of Bone Architecture Variations on Texture at Spatial Resolutions of CT, HR-pQCT, and μCT Scanners

The quantification of changes in the trabecular bone structure induced by musculoskeletal diseases like osteoarthritis, osteoporosis, rheumatoid arthritis, and others by means of a texture analysis is a valuable tool which is expected to improve the diagnosis and monitoring of a disease. The reaction of texture parameters on different alterations in the architecture of the fine trabecular network and inherent imaging factors such as spatial resolution or image noise has to be understood in detail to ensure an accurate and reliable determination of the current bone state. Therefore, a digital model for the quantitative analysis of cancellous bone structures was developed. Five parameters were used for texture analysis: entropy, global and local inhomogeneity, local anisotropy, and variogram slope. Various generic structural changes of cancellous bone were simulated for different spatial resolutions. Additionally, the dependence of the texture parameters on tissue mineralization and noise was investigated. The present work explains changes in texture parameter outcomes based on structural changes originating from structure modifications and reveals that a texture analysis could provide useful information for a trabecular bone analysis even at resolutions below the dimensions of single trabeculae.


Introduction
Quantitative computed tomography (QCT) is an advanced method to measure bone mineral density (BMD) in vivo at various skeletal sites [1]. However, to date the in vivo quantitative analysis of the trabecular bone network remains challenging. For peripheral locations such as the distal radius or tibia dedicated high resolution peripheral QCT (HR-pQCT) equipment with an isotropic spatial resolution of about 130 m exists [2], but long scan times result in frequent motion artifacts and disturb the analysis of trabecular bone structure [3,4].
Analysis of the trabecular network imaged with in vivo techniques, predominantly not only with CT but also with MRI or X-ray films, has received a fair amount of attention in the past. In the majority of reported studies, binarization methods were used to separate bone from soft tissue prior to the calculation of histomorphometric parameters [5][6][7]. However, the spatial resolution of almost all in vivo imaging modalities exceeds the diameter of single trabeculae of about 100 m to 200 m [8][9][10]. Therefore, binarization techniques were avoided in the present work. Instead, texture parameters directly calculated from the gray value distribution of datasets were used.
The texture analysis of trabecular bone is not a new topic and a variety of texture parameters have been used [11][12][13]. However, a systematic validation of accuracy was rarely included in those studies. Such a validation should include the examination of the dependence of texture on different aspects like structure, BMD modifications, image resolution, and noise. In existing studies usually only the aspect of image resolution was discussed. Textural features were compared among datasets with different spatial resolutions acquired either from different imaging systems like micro-CT ( CT) and HR-pQCT [14], from MR acquisitions with different resolutions [15], or by downsampling of high resolution datasets [16].
In the present work, a digital trabecular bone model consisting of rods and plates is introduced to examine quantitatively the ability to assess the trabecular bone structure at spatial resolutions obtainable with CT, HR-pQCT, and CT scanners. The model is generic and can be used to simulate typical architectural alterations occurring in osteoporosis or arthritis including the effects of noise and image resolution. As an example in this contribution, we apply the model to quantify entropy, global and local inhomogeneity, local anisotropy, and variogram slope [17]. The aim of the present work is not to simulate architectural changes for a particular disease.
The basic question investigated in the present work is whether variations in bone architecture can be quantified by the use of texture parameters. For this purpose, various structural modifications were applied to the digital bone model. The influence of noise and spatial resolution was included in the investigation. A key topic in the use of texture parameters is their ability to distinguish differences in trabecular architecture from those in BMD. Differences in BMD can easily be measured in vivo using DXA (reproducibility error: 1%-2% [18,19]) or QCT (reproducibility error: 1%-4% [20][21][22]) techniques but these techniques cannot differentiate whether BMD changes are caused by changes in tissue mineralization or bone architecture. Thus, two further questions are as follows: Can the use of texture parameters differentiate changes in BMD from changes in BV/TV? And can their use differentiate variations in trabecular architecture that result in the same BV/TV?
The hypothesis of this study is twofold. First, the use of texture parameters is useful to characterize trabecular bone structure using clinical CT equipment, where the spatial resolution is typically inadequate to separate individual trabeculae and, second, based on the dependence of texture parameters on structure at a voxel size of 10 m, the corresponding characteristics at larger voxel sizes, that is, lower spatial resolution, can be derived.

Basic Trabecular Bone Model.
The basic trabecular bone model consists of 1000 × 1000 × 1000 isotropic voxels with an edge length of 10 m, resulting in a total volume of 1 cm 3 . The edge length of 10 m represents the typical voxel size obtainable in CT datasets. A CT value of 800 HU is assigned to bone voxels and a value of −50 HU to soft tissue voxels. These values are realistically found in trabecular bone. For comparison also CT values of 200 HU were assigned to bone to investigate the texture accuracy at low bone to soft tissue contrasts. The CT value of soft tissue is between CT values of fat (−100) and water (0) [23], as soft tissue surrounding cancellous bone consists mostly of fat and tissue with water equivalent absorption characteristics.
The basic model was built as a combination of rods and plates representing an average human trabecular bone structure [10,24]. It consists of 11 × 11 cylindrical rods with a diameter of 200 m and a spacing of 700 m. The rods are equidistantly interleaved by nine parallel plates with a thickness of 200 m and a spacing of 1000 m, which are arranged perpendicular to the rods. Some of the cylinders were cut in half ( Figure 1) in order to partially break the symmetry of the structure and to make the model more realistic. This resulted in a bone to tissue volume ratio (BV/TV) of about 20%, which is a typical value found in human epiphyses [25,26]. Obviously, this model is a simplification of human trabecular bone. However, a more detailed model would not provide generality, limiting its applicability.
Different spatial resolutions (here voxel sizes) were simulated by resampling the structure using a bilinear interpolation implemented in ImageJ [27]. Subsequently, Gaussian noise with a standard deviation of 30 HU, a typical noise level found in in vivo CT acquisitions with medium reconstruction kernels [23], was added to the model. The effect of different reconstruction kernels was simulated by varying spatial resolution and image noise. (1) A more plate-like structure was simulated. Five models (PLM2-PLM6; PLM: plate-like model) were created in addition to the basic model (PLM1 = basic model) with decreasing numbers of rods ( Figure 2 PLM): 20%, 40%, 60%, 80%, and 100% of the rods, respectively, were removed from the models. The trabeculae being removed were chosen randomly (uniform distribution), separately for each layer. Plate thickness was increased with decreasing rod number to keep the overall BV/TV constant. This loss of trabeculae oriented in a specific direction resembles the situation in osteoporosis where in the vertebrae greater bone loss occurs in horizontal compared to vertical trabeculae [28].
(2) An increased bone formation with increased BV/TV was simulated. Five models (BFM2-BFM6, BFM1 = basic model, and BFM = bone formation model) were built with dilated rods and plates (Figure 2 BFM). With every new model, plate thickness and rod diameter were increased by 20 m, leading to an increase in BV/TV from 20% of the basic model to 32% of the most dilated model (Table 1) with plate thicknesses and rod diameters of 300 m each. As the global appearance pattern of the bone structure remains quite unchanged, these alterations are expected to have only a little effect on the outcome of the texture parameters, except those which predominantly depend on BV/TV. An increase in BV/TV has been found in osteoarthritis [29].
(3) A more rod-like structure was simulated. Five models (RLM2-RLM6, RLM1 = basic model, and RLM: rodlike model) were created with decreasing numbers of plates and increasing rod thickness (Figure 2 RLM), while BV/TV was kept constant. For the first new model, the central plate was removed and the diameter of all rods was increased. For the models with even fewer plates, the two next central plates were removed. In this way, five additional models with 8, 6, 4, 2, and no plates, respectively, and rod diameters between 200 m and 500 m were created. Transformations between plate-and rod-like trabecular structures are present in a variety of musculoskeletal diseases.    Table 1). The remodeling of trabecular bone is highly load driven (Wolff 's law [30]). As changes in loading are a fundamental factor in bone diseases, the detection of superficial trabecular changes may improve the understanding of the bone state during the disease progress.

Variation of Tissue Mineralization. BMD is an important
parameter, in particular, as it can be measured in vivo. The underlying causes for BMD changes are a change in trabecular architecture characterized by BV/TV or a change in the mineralization of the trabeculae (tissue mineralization). Both effects can be simulated in our model. It is easy to change the mineralization. In the basic model, CT values for bone were varied between 200 HU and 1200 HU in steps of 200 HU. Generally, tissue mineralization is measured in mg/cm 3 , but in this work we stick to HU values to be consistent with the values used in the basic model.

Variation of Noise.
Pixel noise is a fundamental feature of imaging techniques like CT. Noise increases with increasing absorption by high-attenuating objects, lower mAs settings, and smaller slice thicknesses and highly depends on the reconstruction kernel [23]. As noise strongly affects texture, different noise levels, defined here by the standard deviation of a Gaussian distribution, were applied to the basic model ranging from 5 HU to 50 HU in steps of 5 HU.

Texture Analysis.
Five different 3D texture parameters were calculated directly from the gray value distributions of the datasets: entropy, global and local inhomogeneity, local anisotropy, and variogram slope. The inhomogeneity and local anisotropy parameters were already described in [21] but are added here for completeness.

Entropy.
Here, the term entropy refers to the Shannon entropy [31], which measures the information content and is given in bits: where is the gray value range of the image which is divided into 100 equally distributed partitions ( ), called bins, between the minimum and the maximum gray value. With = 100 the number of empty bins or bins with just a few voxels is probably small. This may not be the case for considerably higher as the gray value range of imaging datasets is high (12 bits in our case).
is the probability of one partition , given as the ratio of the number of voxels in bin and the total number of voxels in the image VOI .

Global Inhomogeneity.
Global inhomogeneity (Inhom global ) measures gray value ( ) fluctuations and is equal to the standard deviation of [21]: where is the mean gray value and iterates over all voxels.

Local Inhomogeneity.
In contrast to global inhomogeneity, local inhomogeneity (Inhom local ) measures local gray value variations which are calculated in a 6-neighborhood [21]: 2.5.4. Local Anisotropy. Local anisotropy ( ) represents the variance of directedness in a local neighborhood [21]. It measures the mean angle difference between the gray value gradient ⃗ ] , , , of a single voxel and the mean gray value gradient of its 26neighborhood. Therefore, a local mean gray value gradient ⃗ ] local, , , has to be calculated first: The local anisotropy is then given by Journal of Medical Engineering 5 2.5.5. Variogram Slope. The variogram Var( ) describes the mean gray value difference between voxels in a distance to each other [32]: For every voxel V with gray value , the absolute gray value difference to all Nei, neighbor voxels V with gray value is calculated. The total sum of absolute gray value differences over all voxels V is finally normalized by twice the total number of neighbors as all neighbor pairs are considered twice. The variogram is calculated as a function of increasing voxel distance . The analyses of several CT datasets revealed a plateau of the variogram for > 3. Therefore, here, the slope of Var( ) is calculated from linear least squares fit in the interval 1 ≤ ≤ 3 voxels. Obviously, in case of smaller voxel sizes, for example, in CT data, higher distances can be used. Nevertheless, as the ultimate aim of the present work was the validation of texture parameters measured in clinical CT data d was kept constant even when analyzing data with different resolutions. Moreover, the consideration of larger d rapidly exceeds acceptable calculation times.
According to their definition and concerning only structural but not mineral changes, it can be expected that entropy and global inhomogeneity mainly depend on BV/TV rather than on the spatial distribution pattern of the trabeculae. Local inhomogeneity, local anisotropy, and variogram slope on the other hand are expected to be rather independent of BV/TV and almost exclusively describe the pattern of the structure.
These five texture parameters remained after a preselection of a higher amount of parameters (e.g., fractal dimension and lacunarity). Excluded parameters showed irregular behavior and high sensitivity on small structure variations. Figure 3 shows the effects of structural modifications on texture parameters at a voxel size of 10 m and a noise level of 30 HU. Results for the basic model are encircled. Obviously, for a given parameter, the basic model values are identical in all graphs, for example, 4.14 for entropy. For better comparison, for a given texture parameter, the scaling is kept the same for all models. The results look almost identical if bone HU values of 200 instead of 800 are used, although absolute values differ.

Dependence on Structure.
For each graph in Figure 3, a linear regression was calculated. An insignificant slope or a nonmonotonic dependency of a texture parameter on structural variations indicates that this parameter may not be suited to pick up structural changes.
The regression slopes of entropy and global inhomogeneity do not significantly differ from zero for PLM and RLM model modifications, in which BV/TV was constant. All other regression slopes from Figure 3 significantly differ from zero. Consequently, and in accordance with their definitions, entropy and global inhomogeneity are independent of structure as long as BV/TV is constant. However, as can be seen from the BFM results, these two texture parameters change with varying BV/TV. There is also a small but significant increase of entropy and global inhomogeneity in the SMM models, in which BV/TV also slightly increases.
In contrast, local inhomogeneity, local anisotropy, and variogram slope depend less on BV/TV as is obvious from the variations for the BFM models. They reflect structural changes for which BV/TV is constant. Slopes are significant for all models but the change is much smaller for the BFM models for which BV/TV changes are the highest. Therefore, there is evidence for the predominant dependence of these three parameters on structure patterns rather than on BV/TV. Consequently, they are able to differentiate differences in trabecular architecture that result in the same BV/TV. Of course, the regression slopes of the four model changes (PLM, BFM, RLM, and SMM) are not directly comparable with each other because their independent variables are different.
Further, a decrease (increase) in local inhomogeneity is always accompanied by an increase (decrease) in local anisotropy and a decrease (increase) in variogram slope. This redundancy suggests the use of only one of these three parameters for the texture analysis of cancellous bone.

Dependence on Tissue Mineralization.
Texture results for varying tissue mineralization for the basic model at a voxel size of 10 m are shown in Figure 4. For better comparison, % changes between results for 800 HU and 1200 HU are added to each graph in Figure 4. All parameters show an almost linear dependence on tissue mineralization. Entropy decreases with increasing mineralization and global and local inhomogeneity as well as variogram slope increase. Local anisotropy remains almost constant, whereas global inhomogeneity and variogram slope change very strongly. Figures 5 and 6 show the same graphs as Figure 3 but at voxel sizes of 90 m corresponding to in vivo HR-pQCT and 250 m corresponding to in vivo QCT datasets. At larger voxel sizes, that is, at lower spatial resolutions, the absolute values of texture parameters change. Values for entropy, local inhomogeneity, and variogram slope increase in comparison to 10 m for all models whereas values for global inhomogeneity and local anisotropy decrease.

Dependence on Spatial Resolution.
One important question that can be addressed with the digital model is, Is there a limit in spatial resolution for the applicability of texture parameters and, if yes, where?
The answer to this question depends on whether (1) the interpretation of texture results shall be independent of spatial resolution or whether (2) an interpretation of texture at a given resolution suffices.
In order to better understand scenario (1) Figure 7 shows one example for anisotropy using the PLM model. If the relationship is not monotonic, the corresponding textural parameter is not suited to detect differences in texture at lower spatial resolutions. Thus, for scenario (1), high R 2 values obtained from a linear regression applied to the graphs exemplified in Figure 7 are a prerequisite though not being sufficient for an appropriate interpretation of the texture results calculated from datasets with lower spatial resolutions. R 2 values are shown in Table 2. The table also lists the ratio of slopes obtained from the linear regressions of a given graph in Figure 5 or Figure 6 and the corresponding graph in Figure 3. This is a figure of merit to investigate whether a decrease in spatial resolution causes a stronger or weaker dependence of texture on structure differences.
A stringent requirement to apply a resolution independent interpretation of texture results would entail that not only the relative change within a given structure modification, that is, for one specific graph in Figure 3, was resolution independent but also the relation among different types of changes, that is, for multiple graphs, was resolution independent. In other words all fields in Table 2 should indicate significance and all slope ratios should be positive. Obviously, this requirement is not fulfilled for any parameter. As a consequence only the resolution dependent interpretation of texture parameters can be used. At different spatial resolutions, different texture parameters or different combinations of texture parameters must be used to characterize or differentiate changes in BV/TV, mineralization, and structure.

Dependence on Noise.
The effect of image noise on texture parameters is shown in Figure 8  Percentage changes between values at 800 HU and 1200 HU are given as Δ. Table 2: 2 values of the linear regression analyses for the voxel sizes 90 m and 250 m relative to the reference voxel size (10 m). Significant linear correlations ( < 0.05) are highlighted in bold. In parenthesis: slope ratios obtained from the linear regressions of a given graph in Figure 5 or Figure 6 and the corresponding graph in Figure 3. noise. Local inhomogeneity and variogram slope also linearly depend on noise. Regarding percentage changes, the impact of noise is the lowest on local anisotropy and global inhomogeneity.   mineralization ( Figure 4). All corresponding R 2 values were high (>0.99). Then entropy and global inhomogeneity were calculated for the basic model, BFM3 and BFM6, and local inhomogeneity, local anisotropy, and variogram slope for basic model, PLM3 and PLM6. For all parameters, absolute differences with respect to the basic model were derived. Finally, those deviations in noise (Δ noise) relative to a noise level of 30 HU and in tissue mineralization (Δ mineralization) relative to 800 HU that would result in the same effect as the structural changes above were calculated. Results are given in Table 3. If numbers are high, then the corresponding texture parameter can pick up structural changes pretty independent of changes in noise and/or mineralization.

Discussion
The quantitative analysis of the trabecular bone structure increasingly attracts attention in osteoarthritis [33,34], osteoporosis [12,35], rheumatoid arthritis [36], and other musculoskeletal diseases. However, it remains unclear what texture parameters really quantify cancellous bone, in particular, if they are applied to lower resolution datasets. Therefore, the aim of the present work was to investigate the ability of texture parameters to quantify trabecular bone changes at different levels of spatial resolution. In contrast to the hypothesis, the results show that the interpretation of texture parameters depends on spatial resolution, because their characteristic   response to a change in trabecular structure at 10 m differs from that at 250 m. Without a simulation of this response at the spatial resolution under consideration, that is, without a priori knowledge of the expected response, measured texture results cannot be interpreted, as an increase or decrease of the value of a texture parameter can have multiple causes. As a consequence, the use of a realistic digital trabecular bone model is vital for this differentiation task.
The main aim of a texture analysis is to provide information on trabecular structure in addition to BMD, which alone can easily be quantified by DXA and QCT. That implies the question of how strong texture parameters depend on BV/TV, which is the main determinant of BMD as newly formed bone mineralizes up to 70% of its final value within a few days [37,38]. A strong dependence would limit the additional value of a texture analysis.
At a voxel size of 10 m, a change of entropy or global inhomogeneity indicates a change in BV/TV but not in structure. As entropy is more noise sensitive global inhomogeneity may be the preferred parameter but it cannot differentiate between changes in BV/TV and mineralization. For this task both parameters must be used in combination   as entropy decreases but global inhomogeneity increases with mineralization (see Figure 4). Local anisotropy and variogram slope and to a lesser degree local inhomogeneity are all capable of differentiating two trabecular networks with identical BV/TV but different structure. The variations simulated by the PLM models, more plates and fewer rods, can be better detected than variations simulated by the RLM models, less plates and thinner rods. These are difficult to differentiate from those caused by artificial surface modifications in the SMM models. While a more quantitative analysis is required, probably the measurement of one of these three texture parameters will suffice to allow for the statement that trabecular structure is changed but BV/TV is not. As local inhomogeneity is very sensitive to noise (see Table 3) and its dependence on BFM and RLM variations is somewhat similar, the other two parameters are preferable measurements.
If the dependence on mineralization (see Figure 4) is added to the variety of possible modifications, local anisotropy remains the parameter of choice because its variation with mineralization is small compared to those shown in Figure 3. If on the other hand changes in mineralization need not be differentiated from those in BV/TV variogram slope may be a better choice as its percentage variation is much higher than that of local anisotropy. With increasing voxel size, structure details disappear and the response of texture parameters changes (Figures 5 and 6). The dependence of texture parameters on spatial resolution is affected by the coarseness of the structure. A rather fine structure leads to a more inhomogeneous appearing structure at small voxel sizes, whereas, at large voxel sizes, it leads to a more homogeneous structure due to the resampling process. These homogeneity changes in global appearance strongly affect all texture parameters. As a consequence, with a change in spatial resolution several of the patterns shown in Figure 3 such as variogram slope for the PLM variations or entropy for the RLM variations fundamentally change (compare to corresponding graphs, e.g., in Figure 6) and the question remains, which texture parameters carry structural information at low spatial resolution? In the following only the voxel size of 250 m will be discussed.
At 250 m, entropy and global homogeneity still strongly depend on structure independent BV/TV changes; however, in contrast to a voxel size of 10 m, now entropy also changes with RLM modifications and global inhomogeneity increases with PLM modifications. An increase in entropy or global inhomogeneity is no longer uniquely caused by a BV/TV increase. However, for example, for a concurrent increase of entropy, global and local inhomogeneity, and variogram slope, PLM variations are excluded because this would cause entropy (regarding PLM variations) to decrease. RLM variations are also excluded because local inhomogeneity would decrease. Thus, even at 250 m, a structure independent change of BV/TV should be identifiable, but a combination of texture parameters has to be measured.
In order to further differentiate changes in tissue mineralization from changes in BV/TV, the information in Figure 6 has to be combined with dependence on noise and mineralization. As can be seen from Figure 4, which looks similar at 250 m, increases in mineralization increase global and local inhomogeneity and variogram slope; only for entropy the regression slope is negative. Thus, the use of the former three texture parameters would not really allow separating BV/TV and mineralization changes and this may only be possible by carefully analyzing the entropy results in combination. In principle, the data presented in Table 3 are required for this purpose but these are examples only for PLM (local inhomogeneity, local anisotropy, and variogram slope) and BFM (entropy and global inhomogeneity).
Apart from changes in mineralization and BV/TV, changes in texture parameters caused by surface modifications (SMM) are small in relation to other structural changes. At a spatial resolution of 250 m, these superficial modifications cannot be detected and will not further be considered. Therefore, only two different scenarios remain: a trabecular structure becomes more plate-like (PLM) or more rod-like (RLM). Local anisotropy is well suited to detect transitions to a more rod-like structure (RLM), especially if the transition starts from a hybrid structure (e.g., from RLM 1 in Figure 6, point at the right, to RLM 4) and not from a structure that already contains very few plates (e.g., RLM 4 to RLM 6). Also, local anisotropy is rather stable with respect to changes in noise and tissue mineralization. In order to quantify changes in BV/TV and transitions to a more platelike structure originating from a hybrid structure (PLM), local inhomogeneity and variogram slope are well suited. Changes in both parameters indicate changes in BV/TV, whereas an exclusive change in variogram slope indicates a transition to a more plate-like structure.
The present study has several limitations. First, a given disease may cause multiple concurrent structural variations. In this case a more complex multivariate analysis will be required. Moreover, real bone structures were not incorporated in the study. Their investigation would complement the simulated data. Also, only five texture parameters were examined. A variety of other parameters exist.
In conclusion, a digital bone model was presented, with which texture parameters for the analysis of cancellous bone in human joints affected by musculoskeletal diseases can be simulated and validated. The presented texture parameters are capable of quantifying changes in the trabecular bone structure even at large voxel sizes of 250 m achievable in in vivo CT acquisitions. The interpretation of texture parameters strongly depends on spatial resolution. The trends of local inhomogeneity, local anisotropy, and variogram slope vary with different resolution levels. Furthermore, changes in noise and tissue mineralization have to be considered when comparing the texture analysis results among different datasets. Nevertheless, the present work demonstrates that in QCT datasets a texture analysis can complement the BMD analysis to improve the diagnosis and monitoring of patients with musculoskeletal diseases by quantifying changes in the fine trabecular network.

Disclosure
The present work was performed in (partial) fulfillment of the requirements for obtaining the degree "Dr. rer. biol. hum. " at the University of Erlangen-Nuremberg.