Correlation of X-Ray Vector Radiography to Bone Micro-Architecture

Besides the overall mass density, strength of trabecular bone depends significantly on its microstructure. However, due to dose constraints in medical CT imaging, it is impossible to gain sufficient information about very fine bone structures in vivo on the micrometer scale. Here we show that a recently developed method of X-ray vector radiography (XVR), an imaging method which uses X-ray scattering information to form an image, allows predictions on the bone microstructure without the explicit need to spatially resolve even individual trabeculae in the bone. We investigated thick human femoral bone samples and compared state-of-the-art μCT data with XVR imaging. A model is presented which proves that XVR imaging yields information directly correlated with the trabecular microstructure. This opens up possibilities of using XVR as a tool to help early diagnosis of bone diseases, such as osteoporosis.

O steoporosis is an under-diagnosed, serious bone disease that affects a substantial amount of the population. Estimated 200 million people worldwide are affected by osteoporosis today 1 . An early diagnosis is crucial to ensure a proper treatment before fractures might occur. The gold standard used in today's diagnostics of osteoporosis is dual-energy X-ray absorptiometry (DXA) and the T-score value calculated from the resulting bone mineral density (BMD) 2 . One often recommended location for diagnosis is the proximal femur area [3][4][5][6] , as it is subject to large mechanical loads and poses an area of high risk for a fracture when suffering from osteoporosis. Even though DXA is able to assess bone quality to some degree, it is only able to give information about the BMD, whereas bone strength also depends on other factors, such as contributions caused by its structure 7 .
Grating-based phase-contrast X-ray imaging is a recently developed method that is able to provide complementary information compared to conventional radiography [8][9][10] . Using grating-based X-ray imaging, one can also obtain information about the local scattering strength of microstructures within a sample on a scale smaller than the actual pixel size. In analogy to visible light microscopy, this has been called dark-field X-ray imaging 11 . The recorded small angle scattering (SAS) signal often varies under rotation of the sample and we have previously called this direction-dependent dark-field imaging X-ray vector radiography (XVR) 12 . As the microstructure of the trabeculae plays an essential role in determination of bone strength (e.g. in pathologies like osteoporosis), we used trabecular bone in this study 12 . In bone research the so-called degree of anisotropy (DA B ) is an important parameter to characterize the structural properties of trabcular bone [13][14][15][16] . One property of the XVR signal is that it often shows a preferred scattering orientation and thus an anisotropy of the signal can be measured 17 . To avoid confusion, the term anisotropy will always refer to the anisotropy of the XVR signal (DA XVR ) in this article.
We have previously shown that thin slices of human trabecular bone do indeed display an anisotropic scattering behaviour under rotation and that XVR is able to record the different scattering orientations. In contrast to absorption based imaging modalities, such as mCT, XVR maintains its ability to give structural information even at coarse resolution and consequently low radiation exposure 12 . Naturally, human bone does not only consist of thin slices, but rather of a complex assembly of trabeculae in various directions. While it is possible to obtain information about this trabecular structure ex-vivo using mCT, this method is not suitable for diagnostics due to the very high dose required to directly resolve the individual trabeculae in great detail. For possible use of XVR as a diagnostical tool, it is vital to have a good understanding of how the XVR signal is being formed within a whole volume of trabecular bone. In this study we compare state-of-the-art mCT of thick human femoral bone samples with the respective XVR results.  Figure 1 sketches the fundamental differences between a mCT system and an experimental XVR setup. For the remainder of this article, we use a Cartesian coordinate system with its z-axis parallel to the propagation direction, y-axis parallel to the grating lines and the x-axis orthogonal to both of them. By using an elaborate numerical analysis, we correlate the XVR signal to the present trabecular microstructure within the sample.

Results
In analogy to conventional absorption-based imaging, which is governed by the Lambert-Beer-Law, the visibility signal recorded during an XVR scan, V, can be described with an exponential function. The oscillating term is used to include anisotropic scattering within the sample 18 : In eq. (1) and (2), (v 2 Q) is the relative angle between the grating lines and the main scattering orientation in the sample. The degree of anisotropy (DA XVR ) of the XVR signal can be defined as the amplitude of the oscillating part (a 1 ) to the constant part (a 0 ), which by this definition assumes values between 0 and 1: For our measurements, human femur head bone samples were acquired from three different donors at the pathological institute of Klinikum Rechts der Isar, Munich. Eight cubes (P1-P8) of 1 cm 3 each were cut out around a defined axis through the centre of the femoral head as shown in fig. 2(a). Furthermore, a 3-D mCT scan rendering (b) as well as two cuts (c,d) are displayed. In agreement with local laws and institutional requirements, the samples were harvested to be used for educational and research purposes in consent with their respective donors. In the following, our aim is to connect the X-ray vector radiography signal to the actual microstructure. Starting from raw highresolution CT absorption data from trabecular bone samples, we have a three dimensional volume of linear attenuation coefficients m(x, y, z). We assume that scattering occurs mainly at borders between different materials in the sample, i.e. bone and soft tissue/ air, which are perpendicular to the incident X-ray beam. In the case of the investigated bones, the different materials could be easily distinguished by their linear absorption coefficient since the predominant materials were bone, soft tissue and air.
From the two gradients of m n perpendicular to the incident beam, G x,n and G y,n , we can calculate the direction of biggest change Q n (x, y) and its magnitude A n (x, y) in different planes n: Using the superposition principle for XVR imaging 18,19 to add the individual planes, we can calculate a function y(x, y), from which the orientation and anisotropy of the XVR signal can be derived: An elaborate description of the model is given in the Methods section. XVR measurements were performed on all eight cubes from three orthogonal sides each, resulting in a grand total of 24 independent scans. Likewise, the described model was applied to the mCT data along the three corresponding planes for each sample. A comparison between these two methods for specimen P4 can be seen in fig. 3. The main structural orientation is coded in the colour wheel, i.e. a horizontal structure that scatters in the vertical direction shows as light blue whereas a vertical structure appears as red. The brightness corresponds to the anisotropy, a higher anisotropy of the signal resulting in a brighter colour. As a consequence of this, completely isotropic scattering shows as black on the images. Additionally, the XVR and mCT calculation results of the remaining samples P1-P8, excluding P4, are visualized in fig. 4 and fig. 5, respectively.
We calculated the average anisotropy over a square area in the centre of the bone as indicated in fig. 3(a). The values for XVR measurements and calculations based on mCT data are listed in table 1. The correlation between the averaged anisotropy of both methods is illustrated in fig. 6(a). Each dot represents a pair of corresponding XVR/mCT data from one of the eight bone cubes and one of the three available sample orientations. The linear regression line has coefficients and respective standard errors a 5 (0.018 6 0.019) and b 5 (0.827 6 0.089). Figure 1 | Schematics of the micro-CT system (a) and experimental XVR setup (b) used for the measurements. The micro-CT system utilizes a microfocus X-ray tube and divergent illumination to magnify a sample. The XVR setup employs a grating interferometer consisting of three gratings and is run with a conventional X-ray tube. Contrary to the micro-CT system, the sample is rotated in the plane perpendicular to the optical axis.  The mean orientation w was calculated over the same area and the discrepancy D w of every pair of measurement and calculation was averaged and a mean error of the orientation D w~13 0 was found. Furthermore, to demonstrate the preservation of the scattering information at very coarse resolutions, we rebinned the raw data before doing the stepping curve analysis and XVR processing. This corresponds to using a detector with pixel size of 508 mm 2 in the XVR measurement. The resulting images for one sample orientation are shown exemplarily in figures 3(d) and 3(h). For the mCT calculation the the rebinning parameter was changed accordingly (also see Methods section). An average anisotropy was calculated for the rebinned XVR data and corresponding mCT calculation and the correlation is shown in figure 6(b). The coefficients of the linear regression line are a 5 (0.022 6 0.017) and b 5 (0.816 6 0.085).

Discussion
This study provides evidence that in trabecular bone the edges of the trabeculae in beam direction are indeed responsible for the formation of the X-ray vector radiography signal. We have shown a correct prediction of the preferred scattering orientation and anisotropy using mCT data. Comparing the XVR data to those calculated from the mCT measurement, the structure orientations and anisotropy values are in good agreement. It is to note that an area in which the trabeculae display an outstanding orientation compared to the bulk can both be seen in the XVR measurement and calculated by the presented method. Fig. 4(f) and fig. 5(f) depict XVR measurement and mCT calculation results of a sample that possesses these scattering characteristics in direction 2. Comparison with the raw mCT data reveals that there is indeed a strongly irregular layer of bone present throughout most of the sample in this particular area (see fig. 2(c)). The observed orientation can be traced back to the epiphyseal plate.
The quantitative analysis of the average anisotropy values shown in fig. 6 shows a good correlation between the two methods. We have shown that the correlation does not suffer as the pixelsize of the XVR data is increased and the presented model can explain the formation of the XVR signal even for very low resolutions. Furthermore, the anisotropy was calculated correctly not only for the average values but also for the distribution over an image. The distribution is not necessarily equal over the entire sample, e.g. the XVR measurement shown in fig. 3(e) exhibits a strong decline in anisotropy towards the right, a feature that is also found in the corresponding mCT calculation (see fig. 3(b)). This study of anisotropy also confirms that our initial assumption of completely anisotropic scattering is valid, as isotropic scattering in the sample, i.e. from microspheres of size below the resolution of the mCT scan would lower the anisotropy measured by the XVR setup substantially.
Limitations of the presented methods concern the fact that the dark field signal depends on both the linear absorption coefficient m and the real part of the refractive index n 20 . In the calculations we assumed a parallel beam geometry, whereas a laboratory setup like the used XVR setup in general has a cone beam geometry (although with a small cone angle of 2u). In a more divergent geometry, a summation of the individual slices along one axis would not be correct any more and needs to be kept in mind for large samples. The presented study uses sufficiently small samples to warrant the assumption of a parallel beam addition, as the beam offset within the sample stays below the length scale of an individual pixel. Vector Radiography (XVR) images (e-h). Three orthogonal sample orientations of the same specimen are shown. The colour encodes the orientation of the strongest signal, the brightness codes the anisotropy. A structure only scattering into the horizontal plane appears red with maximum brightness, whereas a completely isotropic signal appears dark. A decline in anisotropy towards the right side can be seen in images (b) and (f) as a decline in brightness. Images (c) and (g) show closely matching orientations of mCT calculation and XVR imaging for varying structural orientation. The area over which the anisotropy was averaged for a quantitative analysis is indicated in (a). The artefacts around the sample are caused by the sample holder. The effects of a greatly increased pixel size is shown for XVR measurements (h) and the corresponding mCT calculation (d).
Figure 4 | mCT surface rendering and three XVR images from three orthogonal directions of specimen P1-P8 (a-g). Sample P4 is shown separately in figure 3. Various features can be seen in the different samples and orientations. For example, an area with structural orientation deviating by almost 90u from the surrounding area is visible in the samples shown in (f) and (g) from direction 2 and 3, respectively. The raw mCT data show a strongly irregular trabecular structure present in these particular areas (also see fig. 2(b-d)  In conclusion, we have shown that the XVR signal of trabecular bone is primarily given by the existing trabecular structure in thick samples. This validation supports the idea of using X-ray vector radiography to gather supplementary information about the trabecular structures of femoral bone. Contrary to mCT, this additional structural information can be obtained without the need of having to directly resolve the structures. The most striking advantage of XVR is that it is possible to acquire this scattering information at very coarse resolutions, comparable to DXA. Information about the bone microstructure, next to the bone density information obtained with DXA, is currently not easily accessible and, as we have shown, the scattering orientation and anisotropy of the XVR signal are two properties related to this microstructure. This shows up the possibilities of XVR radiography as a possible supplementary tool in the field of bone quality assessment, and we believe that XVR might help improve diagnosis of osteoporosis in the future.

Methods
X-ray absorption micro CT. High resolution X-ray computed tomography (mCT) data of the samples was obtained using a vjtomejx s 240 system by GE Sensing & Inspection Technologies GmbH (Germany), installed at the Technische Universität München. A sketch of the system is shown in figure 1(a).
For this purpose, the X-ray tube was operated at an acceleration voltage of 90 kV and current of 150 mA. The geometry of the setup and resulting constraints in resolution lead to an isotropic voxel size of 18.53 mm for every sample. This is sufficient to resolve the trabecular structure of human vertebral bones 21 . The duration of a single CT scan was approximately 70 minutes with 2000 projections being taken over 360u. Figs. 2(b-d) show a rendered 3D volume of sample P7 as well as two cuts through its middle. As can be seen, the trabecular structure is well resolved.
X-ray vector radiography. The used setup is sketched in figure 1(b). A detailed description of a typical XVR setup and the measurement procedures can be found in e.g. 17 .
Our experiments were performed at the Technische Universität München at a laboratory setup operating a high power X-ray tube by Comet AG, Switzerland. For the experiments, the X-ray tube was operated at 60 kV and 30 mA with a 3.0 mm aluminum filter. The two absorption gratings consisted of 160 -170 mm thick gold lines on 500 mm (G 0 ) and 150 mm (G 2 ) thick silicon substrates, filled with SU-8 photoresist. Complementing the set of gratings is a phase grating (G 1 ) made of 8 mm nickel lines on a 200 mm thick silicon substrate. This results in a phase shift of p/2 at the design energy of 45.7 keV. The absorption gratings had a periodicity of 10 mm and the phase grating had a periodicity of 5.0 mm. All gratings had a duty cycle of 0.5. The setup is built to operate at the first fractional Talbot distance with equidistant grating locations G 0 to G 1 and G 1 to G 2 of 92.7 cm. The samples were placed 28 cm behind G 1 towards G 2 . A Varian PaxScan 2520D with a CsI scintillator and pixel pitch of 127 mm was used as X-ray detector.
Micro CT data processing. There is a variety of edge detection filters in image processing, we used the Roberts Cross Operator to locate differences in density which shall be introduced here briefly 22 : Let m n be the slice with index n of the X-ray absorption data perpendicular to the beam direction in the XVR setup (z-axis). Then we can define for every slice the two gradients G x and G y as a convolution of the absorption data with the respective filter kernel: From these two gradients of m n we can calculate the direction of biggest change Q n (x, y) and its magnitude A n (x, y): We assume that every edge within the sample produces a fully anisotropic signal (a 0 5 a 1 ) in direction Q n (x, y) with magnitude A n (x, y). The XVR signal of the whole sample can be seen as a superposition of each slice's contribution 18,19 . Consequently, it is necessary to find a way to combine the calculated signals derived from the actual microstructure and originating from different slices in the propagation direction of the beam. It is also known that any XVR signal can be described as a linear superposition of various fully anisotropic signals a 1,i cos (2Q ,i ), which do not necessarily have the same phase, Q ,i , or strength, a ,i The factor 2 in the cosine's argument in eq.(8) accounts for the fact that due to the grating structure the signal does not vary with respect to a 180u rotation of the sample around the beam. With the previously mentioned assumption that scattering occurs fully anisotropic at each gradient and using the complex notation for a wave we can calculate a function that is expected to be proportional to the oscillating part of the XVR signal: From these oscillating functions y(x, y) we can extract the main structural orientation one would derive from the XVR signal Q tot (x, y), which is equal to their argument in any given pixel P(x, y). The second important information one can obtain using XVR is the anisotropy of the signal. As defined earlier, it is the ratio of the oscillating part (a 1 ) to the average signal (a 0 ) after taking the negative logarithm of the recorded XVR signal. The absolute value of y(x, y) is proportional to a 1 (x, y), whereas the average part a 0 (x, y) is proportional to the sum of the amplitudes A n (x, y) along the beam axis under the assumption of completely anisotropic scattering. Therefore, we can state that: a 0 x,y ð Þ! X n A n x,y ð Þ: ð12Þ Since the proportionality constants in eq.(11) and eq.(12) depend only on the relative values given by the mCT software and thus are the same, we can get the anisotropy by combining these equations: When comparing the anisotropy values calculated using this method with the ones directly obtained by a XVR measurement one needs to consider two effects which reduce the recorded anisotropy. Firstly, most pixelated detection mechanisms such as the used detector possess a certain lateral smearing of their response to a point signal, which can be described in terms of a point spread function (PSF). For this reason, a PSF was applied to the two-dimensional array of y(x, y) before calculating the orientation and anisotropy values by using the convolution operator (fl): The unknown PSF of the detector was approximated with a two-dimensional Gaussian function of defined full width at half maximum (FWHM) at 10 pixels, corresponding to 185 mm. Secondly, there is a difference in resolution between the mCT and the XVR scan. For the purposes of this work, the resolution of the mCT scans was much higher compared to the XVR setup. This needs to be considered in the calculations, as a single detector pixel of the XVR setup integrates the recorded intensity over an area of 127 3 127 mm 2 . A binning of 5 3 5 pixels (20 3 20 pixels for the data set presented in Figure 6(b)) into a single pixel was applied on the complex valued function y(x, y) before calculating the orientation and anisotropy values. The size of the binning was chosen to account for the magnification introduced in the used XVR setup. Both of these effects can be seen as a lateral superposition of the signals in the individual pixels. These corrections affect mainly the calculated average anisotropy, the average orientation of the calculated signal remains nearly constant with varying parameters of the corrections. Samples with lower average anisotropy generally show a bigger decrease in anisotropy than those with a higher average anisotropy. This can easily be explained by the fact that samples with an overall lower anisotropy also had a greater variance in their scattering orientation. The reduction in anisotropy is stronger when signals with greatly different scattering orientations are summed up. As the width of the PSF is increased, a summation is done over a greater number of individual scattering orientations. Consequently, samples with a rapidly varying scattering orientation show a stronger decline in average anisotropy than very homogeneously scattering samples. In figure 3, orientation (a) of the sample resembles a very inhomogeneous scattering orientation, whereas orientation (b) shows a very homogeneous distribution of the scattering orientation. The overall average decrease of the average anisotropy for a PSF with FWHM of 5 pixels compared to the uncorrected case was 0.025 with a standard deviation of 0.013. For a PSF with FWHM of 10 pixels, these values nearly double with a decrease of 0.041 and standard deviation of 0.024, respectively. The mean orientation of a sample was calculated by applying equations (9) and (10): Using this notation, every pixel is weighted equally, although different weighting factors may be used, e.g. the average scattering strength a 0 .