Virtual supersampling as post-processing step preserves the trabecular bone morphometry in human peripheral quantitative computed tomography scans

In the clinical field of diagnosis and monitoring of bone diseases, high-resolution peripheral quantitative computed tomography (HR-pQCT) is an important imaging modality. It provides a resolution where quantitative bone morphometry can be extracted in vivo on patients. It is known that HR-pQCT provides slight differences in morphometric indices compared to the current standard approach micro-computed tomography (micro-CT). The most obvious reason for this is the restriction of the radiation dose and with this a lower image resolution. With advances in micro-CT evaluation techniques such as patient-specific remodeling simulations or dynamic bone morphometry, a higher image resolution would potentially also allow the application of such novel evaluation techniques to clinical HR-pQCT measurements. Virtual supersampling as post-processing step was considered to increase the image resolution of HR-pQCT scans. The hypothesis was that this technique preserves the structural bone morphometry. Supersampling from 82 μm to virtual 41 μm by trilinear interpolation of the grayscale values of 42 human cadaveric forearms resulted in strong correlations of structural parameters (R2: 0.96–1.00). BV/TV was slightly overestimated (4.3%, R2: 1.00) compared to the HR-pQCT resolution. Tb.N was overestimated (7.47%; R2: 0.99) and Tb.Th was slightly underestimated (-4.20%; R2: 0.98). The technique was reproducible with PE%CV between 1.96% (SMI) and 7.88% (Conn.D). In a clinical setting with 205 human forearms with or without fracture measured at 82 μm resolution HR-pQCT, the technique was sensitive to changes between groups in all parameters (p < 0.05) except trabecular thickness. In conclusion, we demonstrated that supersampling preserves the bone morphometry from HR-pQCT scans and is reproducible and sensitive to changes between groups. Supersampling can be used to investigate on the resolution dependency of HR-pQCT images and gain more insight into this imaging modality.


Introduction
Quantitative assessment of the trabecular bone microstructure is a valuable tool in bone research because a number of bone diseases act directly on the trabecular bone surface which causes alterations in the bone microstructure.
The current standard technology to quantify human three-dimensional bone morphology is micro-computed tomography (micro-CT) ex vivo where no radiation issues have to be respected. With advances in in vivo imaging technologies, the quantitative determination of the bone morphology has entered the clinical setting. High-resolution, peripheral quantitative computed tomography (HR-pQCT) is a promising clinical tool for the monitoring of the microstructure in bone diseases and their treatments in vivo in patients. HR-pQCT provides an image resolution where single trabeculae can be resolved and the bone microstructure can be quantified using the techniques originally developed for micro-CT [1,2]. Trabecular and cortical microstructural and biomechanical parameters gained from HR-pQCT scans have been validated in comparison to scans obtained using micro-CT at distal radius, tibia and calcaneus [1,[3][4][5][6][7]. Metcalf et al. [3] found, comparing trabecular structural parameters of HR-pQCT with micro-CT, very strong correlations of BV/TV, and moderate correlations for Tb.Th and Tb.N. Christen et al. [4] developed a biomechanical analysis tool for micro-CT resolutions [5] which provides detailed insight into the reasons for a certain form of the bone microstructure. Assuming bone forms tissue at high-load locations and resorbs tissue at low-load locations, the approach finds a set of load cases that lead to the most uniform tissue loading using optimization [6][7][8]. Christen et al. investigated the tool's voxel size dependency, reproducibility, and sensitivity on HR-pQCT data and showed with this that their tool can be applied at HR-pQCT resolutions.
Further image processing techniques that were developed for micro-CT images at nominal image resolutions of about 10 μm require voxel sizes smaller than provided by HR-pQCT. This concerns the assessment of three-dimensional bone formation and resorption rates that has been validated for in vivo micro-CT [9]. With the availability of time-lapsed scans of a patient, it would be great to have access to his dynamic local bone changes per week or month. However, assuming human mineral apposition rates of about 10 μm per week [10], it would take 8 weeks of continuous bone growth until a full new voxel (82 μm) was filled with bone mass. Only dividing the voxel size by two would generate a huge benefit in reducing the error caused by partial volumes effects.
Also, having access to in silico simulation tools that use a patient's microstructure as input and predict how the bone evolves over time under a certain treatment can be of great use for both physician and patient. Such subject-specific in silico simulations to predict and possibly adapt a subject's treatment regime have been explored on preclinical in vivo micro-CT datasets at resolutions of about 10 μm [11].
A number of studies have used downscaling to investigate on the resolution dependency of the bone morphometry [12][13][14][15][16]. Little attention has been paid to the reverse effect, i.e. the upscaling of reconstructed images. Such upscaling could provide a technique so that pQCT scans get access to tools currently not applicable at clinical image resolutions. Further, the investigation of such analysis can help gain more insight into the direct effects of the image resolution on bone morphometry. Here we propose a technique to supersample HR-pQCT images. We hypothesize that the structural bone morphometry is preserved by the proposed technique and we investigate the effect on bone morphometry in terms of precision and sensitivity.

Precision
To test for the precision of bone morphometry parameters on supersampled images, a study consisting of repeated HR-pQCT scans of 14 human radii was re-used [17]. In short, the 14 formalin-fixed intact human cadaveric forearms had been imaged each 3 times with repositioning after each measurement with a HR-pQCT scanner (XtremeCT, Scanco Medical AG, Brüttisellen, Switzerland) providing an isotropic nominal resolution of 82 μm.
All raw datasets underwent a 3D Laplace-Hamming filter (Hamming cut-off frequency = 0.4, weighting factor = 0.5) [17]. The repeated scans were registered to allow evaluation of the same regions of interest. The registered images were supersampled by a factor of 2 (using trilinear interpolation), resulting in a virtual 41 μm image resolution (Fig 1). Both the original (82 μm) and the supersampled (41 μm) grayscale images were segmented at a global threshold of 400/1000 of the maximal gray scale value.
The reproducibility of bone morphometric indices at both image resolutions was expressed in precision errors PE SD , PE %CV , the intraclass correlation coefficient ICC and the 95%-confidence intervals ICC lower and ICC upper [19]. Furthermore, the bone morphometry of images at virtual 41 μm was compared to the original image resolution of 82 μm. Significant differences to the measured image resolution of 82 μm were calculated by paired Students t-test. Furthermore, linear regression analysis was performed between bone morphometry from the 82 μm and 41 μm image resolutions. All statistical comparisons were performed using R (Statistical Software package, Vienna, Austria [20]). P-values less than 0.05 were considered significant.

Sensitivity
Sensitivity of the proposed method was assessed by reanalysis of a prior study, where 100 postmenopausal women experienced a Colles' fracture at median age 63 years [21]. The fracture cases were frequency-matched to 105 age-matched control women at median age 62 years. Written informed consent was obtained from all subjects before participation in the study. The study was approved by the Mayo Clinics Institutional Review Board and the present reanalysis was based on de-identified data. The distal radius of the study subjects was measured by HR-pQCT at an isotropic voxel size of 82 μm (XtremeCT, Scanco Medical AG, Brüttisellen, Switzerland). In all study subjects, the non-fractured side (fractured subjects) or the non-dominant arm (controls) was assessed.
All raw datasets underwent a 3D Laplace-Hamming filter (Hamming cut-off frequency = 0.4, weighting factor = 0.5) [17]. Afterwards, the images were upscaled (supersampled) by a factor of 2 using trilinear interpolation. Both the original and the interpolated grayscale images were segmented at a global threshold of 400/1000 of the maximal gray scale value [21]. A trabecular mask was used to extract morphometric parameters in the trabecular compartment. The structural parameters BV/TV, BS/BV, Tb.Th, Tb.Sp, Tb.N, SMI, DA and Conn.D were calculated. Unpaired Student's t-test was used to detect differences in means between fracture cases and controls. Linear correlations between 82 μm and 41 μm resolutions were assessed. Pairwise comparisons between resolutions were calculated.

Precision
All precision values including mean and standard deviation, PE SD , PE %CV with its upper and lower confidence intervals as well as the ICC values including their upper and lower confidence intervals can be found in Table 1 Table 2 lists the coefficients of determination R 2 , the absolute %-errors between 41 μm and 82 μm resolution and the intercept functions by which the 41 μm resolution parameters have to be scaled to fit the original 82 μm voxel size parameters. Table 3 and Fig 4 show the results of the clinical data (fractured /non-fractured) measured at 82 μm image resolution and supersampled at virtual 41 μm resolution. At 82 μm resolution, bone morphometry in the fractured group performed in all morphometric parameters significantly worse than in the control group, i.e. lower BV/TV, Tb.Th, Tb.N and higher BS/BV, Tb.Sp and Conn.D (p < 0.05). Tb.Th differed only 2.71% between groups but was still significant. At 41 μm resolution, the same structural parameters were assessed. In the between-groupscomparison (fractured vs. non-fractured), the same significances as in the 81 μm comparison were detected. Only the difference of Tb.Th was with 1.48% not significant any more.

Discussion
In this study, we showed that linear supersampling of HR-pQCT scans at 82 μm voxel size to virtual 41 μm voxel size as a post-processing step is reproducible and sensitive to changes between groups. Human radius trabecular BV/TV and other trabecular parameters assessed at an isotropic 82 μm resolution were used as input structures. They were supersampled at virtual 41 μm resolution by trilinear interpolation.  Virtual supersampling preserves bone morphometry from peripheral quantitative computed tomography 82 μm HR-pQCT BV/TV and supersampled BV/TV had a strong correlation (R 2 = 0.99-1.00) and an absolute error of 4.30% -4.31%. BV/TV at 41 μm voxel size was consistently higher than at 82 μm voxel size. This phenomenon supports the findings of previous studies [3,13,22,23] which validated HR-pQCT BV/TV with 'gold standard' micro-CT at distal radius, Table 3. Sensitivity. Comparison between 100 postmenopausal women and 105 controls in terms of their trabecular bone structure parameters; values are given as mean ± standard deviation and mean %-difference between groups ( � p < 0.05, �� p < 0.01, ��� p < 0.001).

BV/TV [%]
15.15±5.31 18  Virtual supersampling preserves bone morphometry from peripheral quantitative computed tomography tibia and calcaneus. Metcalf et al. [3] suggested that the lower BV/TV at HR-pQCT resolution results either from the global threshold, beam hardening artifacts or the lower signal-to-noiseratio resulting from a higher number of partial volume voxels. In our study, we can exclude both beam hardening artifacts and threshold-related differences. Thus we can specify the reason for the lower BV/TV in HR-pQCT images to the lower signal-to-noise-ratio where those parts of the bone volume that do not fill up a whole voxel get lost. In our study, Tb.Th was consistently overestimated (4.2% -4.56%) at virtual 41 μm resolution and Tb.N was consistently underestimated (-7.47%--9.3%). Again, these findings fit with previous studies which found that HR-pQCT obtained Tb.Th is lower and Tb.N is higher than micro-CT obtained parameters [3,24]. As our study scaled from lower to higher image resolution, again we can point the reason for those findings to the lower sampling frequency resulting where the partial volumes take more effect.
Conn.D was not affected by the sensitivity study but increased at 41 μm resolution at the precision study. A reason for this can be that fine structures connected with each other can be grasped better at a higher sampling frequency which is obtained by supersampling of the image resolution.
We recognize there are several limitations to the study. First, we validated the virtual 41 μm resolution by comparison to the original 82 μm HR-pQCT input. A comparison to 'gold standard' micro-CT scans would potentially have given even more insight into possible benefits in terms of accuracy of the proposed technique. However, several studies have investigated on the differences between micro-CT and HR-pQCT and found strong correlations [3,13,23,29,30]. Nevertheless, with this constraint, our conclusions can only be drawn with respect to HR-pQCT image quality. Second, we used the first generation Xtreme-CT device which has an isotropic voxel size of 82 μm [31] while the second generation Xtreme-CT would have an isotropic resolution of 61 μm itself [32]. Third, our analysis was restricted to structural parameters and further evaluation is needed to investigate whether such supersampling can also be applied for mechanical, dynamic or other structural parameters than those presented.
In conclusion, we have investigated the precision and sensitivity of a supersampling method for human HR-pQCT scans at 82 μm voxel sizes. The investigation showed minimal Table 4. Sensitivity. Linear regression analysis between 82 μm and 41 μm resolution is shown as coefficient of determination R 2 , the intercept functions of the structural parameters and mean %-errors ( � p < 0.05, �� p < 0.01, ��� p < 0.001). differences between image resolutions of the static bone morphometry which can be assigned to the sampling frequency. The advantage of this technique is that other disruptive factors than the sampling frequency can be denied. In terms of sensitivity, all structural parameters except trabecular thickness revealed the same differences between patient groups as at original HR-pQCT resolution. We conclude that the proposed supersampling technique is reproducible and sensitive to changes between groups.

S1 File Data. All relevant data including statistical analyses and R scripts.
(ZIP)