## Materials and Methods
### Image Processing and Deblurring Algorithm  
\(\upmu\)CT images were cropped to only include five slices from the top face of the BSE slides to facilitate comparison and reduce computational expense [34]. 3D Slicer software (3D Slicer 4.8.1) was used to transform, resample, and register the five \(\upmu\)CT slices and individual BSE images of each sample. A single \(\upmu\)CT plane was extracted from registration with the BSE image. Registration error was measured using fiducials placed on paired BSE and \(\upmu\)CT images. Five fiducials were placed on all image pairs, and the distance between fiducials was averaged for each sample.  
An established deblurring algorithm (based on Richardson-Lucy deconvolution) was applied to the enhanced \(\upmu\)CT images [18,26]. The PSF was modeled as a Gaussian and using a custom module ([[https://bitbucket.org/OrthopaedicBiomechanicsLab/deblurring](https://bitbucket.org/OrthopaedicBiomechanicsLab/deblurring)]([https://bitbucket.org/OrthopaedicBiomechanicsLab/deblurring](https://bitbucket.org/OrthopaedicBiomechanicsLab/deblurring)), accessed on 20 October 2020), the in- and out-of-plane blurring were estimated by placing profiles (implemented as rulers) across thin regions of cortical bone. Profile locations were chosen as thin cortical bone structures with non-bone regions on both sides (air or marrow) such that there were three distinct layers, or greyscale values, across the profile. In-plane blurring was estimated by profiles placed in the sagittal plane, and out-of-plane blurring was estimated by profiles placed in the coronal and axial planes. Theoretical PSF should be equal in identically acquired images (same \(\upmu\)CT scanner using the same parameters). Thus, an average was calculated based on the estimated in-plane blurring in each image (\(0.0067\pm 0.0007\) mm). Due to their thinness (\(\sim\)500 \(\upmu\)m), the slides did not have enough visible structures in the coronal and axial directions to estimate out-of-plane blurring. A whole fresh frozen vertebra (age-matched, identical HeLa cell injection protocol) was \(\upmu\)CT imaged with the enhanced scan parameters. Out-of-plane blurring was estimated from this image as 0.0068 mm. All enhanced \(\upmu\)CT images were deblurred using the average in-plane and single out-of-plane blurring parameters as inputs to the Richardson-Lucy deconvolution algorithm.  
The \(BaSO_{4}\) contrast was segmented from the \(\upmu\)CT (baseline, enhanced (Figure 1a), and deblurred) and BSE images (3D Slicer segmentation editor). First, the whole bone was segmented with automated thresholding using the Otsu method [35, 36, 37]. The segmentation was shrunk by 50 \(\upmu\)m from the outer edge (Figure 1b, light green) of the cortical shell to remove non-specific contrast not caused by microdamage [10,11]. \(BaSO_{4}\) was segmented inside the shrunken label field (Figure 1c, light blue) using thresholding at a greyscale intensity of 253 for the BSE images and \(\sim\)2500 mgHA/cm\({}^{3}\) for the \(\upmu\)CT images. The BSE images were 8-bit images with no intensity calibration for \(BaSO_{4}\). As the difference in density between barium and calcium is high, and the greyscale intensity bins are limited to the 8-bit range for the BSE images, the pixel intensities corresponding to \(BaSO_{4}\) were concentrated around the upper intensities of the image. As such, pixel intensities of 253 or above were sufficient in segmenting \(BaSO_{4}\) from the BSE images.