Adaptive image segmentation reveals substantial cortical bone remodelling during early fracture repair

ABSTRACT

The goal of this study was to develop an image analysis algorithm for quantifying the effects of remodelling on cortical bone during early fracture healing. An adaptive thresholding technique with boundary curvature and tortuosity control was developed to automatically identify the endocortical and pericortical boundaries in the presence of high-gradient bone mineral density (BMD) near the healing zone. The algorithm successfully detected boundaries in more than 47,000 microCT images from 12 pairs of healing ovine osteotomies and intact contralateral tibiae. Resampling techniques were used to achieve data dimensionality reduction on the segmented images, allowing characterisation of radial, circumferential, and axial distributions of cortical BMD. Local (transverse slice) and total (whole bone) remodelling scores were produced. These surrogate measures of cortical remodelling derived from BMD revealed that cortical changes were detectable throughout the region covered by callus and that the localised loss of cortical BMD was highest near the osteotomy. Total remodelling score was moderately and significantly correlated with callus volume and mineral composition (r > 0.64, p < 0.05), suggesting that the cortex may be a source of mineral needed to build callus.

KEYWORDS:

• Bone fracture healing
• cortical remodelling
• adaptive thresholding
• micro-computed tomography
• ovine osteotomy

Introduction

Secondary bone healing is widely described as a four-stage process: 1) inflammation, 2) soft callus formation, 3) hard callus formation, and finally 4) remodelling (Schindeler et al. 2008; Pivonka and Dunstan 2012; Claes et al. 2012; Einhorn and Gerstenfeld 2015; Ghiasi et al. 2017). However, radiographic and mechanical analysis have revealed that the mineralised collagen matrix undergoes dynamic remodelling during the active repair stage (Tobita et al. 2012). In addition, cortical bone adjacent to the callus undergoes a remodelling process in parallel with – not after – callus formation. Nanoindentation testing of sectioned samples and non-destructive image analysis both reveal localised increases in porosity and decreases in bone mineral density (BMD) in cortical bone immediately adjacent to the fracture line (Preininger et al. 2011; Ren et al. 2022).

Previous investigators have speculated about the mechanistic origins of cortical adaptation during fracture repair. The calcium and phosphate needed for callus mineralisation could be provided from the immediately adjacent intact bone matrix (Preininger et al. 2011). Cortical adaptation could also arise in response to mechanoregulatory signals that induce resorption at the fractured bone ends to reduce gap strain in early healing (Augat et al. 2021). In ovine osteotomies, histology suggests that cortical porosification and bone end resorption may be occurring simultaneously, but these observations are qualitative and restricted to single slices with a limited field of view (Kaspar 2005; Manjubala et al. 2009; Peters et al. 2010; Einhorn and Gerstenfeld 2015; Inglis et al. 2022; Ren et al. 2022).

A major barrier to the advancement of fundamental research on coupled bone repair and remodelling is the technical challenge of measuring remodelling without destroying a sample. Previously, we used micro computed tomography (µCT) data from osteotomized sheep to measure an average 23% drop in cortical BMD just proximal to the osteotomy compared to the contralateral limb (Ren et al. 2022). In these samples, there was almost no difference between the limbs at the proximal and distal ends of the bones, far from the injury site. This finding suggests that remodelling intensity decays as the distance from the osteotomy increases. However, the global extent of cortical remodelling could not be evaluated using our previous image processing method because the remodelling activity itself presents a significant challenge for image segmentation and boundary detection between old cortical bone and callus.

The existing methods in the literature for segmenting the image of a fractured bone have drawbacks. One previous study (Lujan et al. 2010) measured callus size accurately, but was not automatic and required multiple inputs, such as manually selecting a region of interest (ROI), and relied on a fixed threshold. Further, it was applicable to only coronal/sagittal X-rays, rather than axial slices from CT scans. Automatic segmentation of cross sectional images was successful in other studies (Buie et al. 2007; Bissinger et al. 2017; Ren et al. 2022; Hopkinson et al. 2023), but these were again based on selection of fixed segmentation thresholds and therefore applied to a relatively short span of bone near fracture line or to an intact bone, leaving their applicability to more proximal and distal images of fractured bones undetermined. In our exploratory attempts to implement methods based on these prior studies, we confirmed that methods that rely on a constant set of thresholds struggle to accurately detect tissue boundaries in the presence of variable cortical resorption or callus mineralisation.

Accordingly, the technical objective of this study was to develop a new algorithm for adaptive image segmentation to measure the axial dependency of cortical remodelling in osteotomized sheep. Our goal was to automatically segment images across the entire diaphyseal region. To this end, challenges needed to be addressed regarding threshold selection, boundary detection and correction, data sampling, and dimensionality reduction for extraction of summary parameters that define remodelling. The scientific objective of this study was to assess whether the extent of cortical BMD loss during fracture repair is explained by variations in the mineral composition of the callus formed. Specifically, we hypothesised that image-derived measures of cortical remodelling are positively correlated with callus mineral density.

Materials and methods

Animal specimens and imaging

Data from 12 adult female Swiss alpine sheep with an age range of 24–30 months and an average bodyweight of 63–77 kg was obtained in a previously completed experiment (Schwarzenberg et al. 2021). Briefly, all animals had a tibial osteotomy with a 3-mm interfragmentary gap in one limb at mid-diaphysis and an intact contralateral tibia. The fractures were stabilised using stainless steel 12-hole straight veterinary locking compression plates and 3.5-mm locking screws. Animals were maintained for 9 weeks before sacrifice. Implants were removed prior to µCT scanning using an XtremeCT II Micro-CT scanner (Scanco Medical AG, Bruettisellen, Switzerland) with an X-ray voltage of 68 kVp and X-ray current of 1470 µA). The diaphyseal segments (length 150–162 mm) of the operated and intact tibiae were scanned at an isotropic resolution ($L_p$) of 60.7 µm. Each square slice image had 1654 or 1660 pixels along its row/column. Grayscale image pixels contained a 16-bit integer as pixel value ($\mathit{v}$). A phantom (Scanco KP70 phantom, QRM) was scanned in the same scanner at identical settings to quantify the constants needed for converting Hounsfield units (HU) to calibrated BMD. Animal experiments were conducted according to the Swiss laws of animal protection and welfare and authorised by the local governmental veterinary authorities (License No. ZH 183/17).

Image preparation and registration

Each µCT scan contained about 2500 slice images in a structured DICOM format. The images were analysed in Python 3.11.1 using pydicom 2.3.1 package, numpy 1.24.2 package, and opencv-python 4.7.0.72 library.

The initial processing step was a registration procedure to match slices in each operated (op) limb to slices in the corresponding intact contralateral (con) limb. Axial registration was challenging because although sheep have bilaterally symmetric skeletal anatomy, subtle left-right differences in cortical area and shape are known to emerge during development (Becker et al. 2020). The excised and scanned samples were also irregular in length. Therefore, the tibial nutrient artery canal (NAC) was identified as a common anatomical feature that can be easily paired to repeatably achieve axial registration between left-right limb pairs, before any segmentation (boundary detection) has been performed. The NAC appears as an intraosseous hole over several consecutive axial slices and its location was used to define the axial offset ($h_{\mathit{off}})$ between the operated and intact contralateral image sets:

(1) $h_{\mathit{off}}=h_{NAC}^{con}-h_{NAC}^{op}$(1)

where $h_{\mathit{NAC}}$ is the distance between the first proximal slice and the first slice showing a full NAC hole. The absolute value of the axial offset ranged from 70 slices (4 mm) to 530 slices (35 mm). Figure 1 illustrates the registration procedure for a representative operated-contralateral pair.

Figure 1. Quantitative analysis of remodelling in an operated ovine tibia relative to the intact contralateral tibia from the same animal required axial registration of µCT slices from two separate scans. Scans were aligned by matching the level of the nutrient artery canal (annotated NAC) in the two scans. In this example animal, the axial offset is $h_{\mathit{off}}$ = 135 slices. Numbers on each image refer to the slice position in image stack, with white text for intact and red text for operated.

Segmentation and cortical boundary detection

The paired operated-contralateral slice images had differences that introduced challenges for segmentation and cortical boundary detection. The contralateral scans contained only cortical bone, with high-contrast boundaries in every slice. In contrast, the operated scans had callus of varying density inside and/or outside the cortex. Remodeling at the endocortical and/or pericortical boundaries was an obvious barrier to boundary detection (Figure 1, operated slice #1320). Slices near the osteotomy were sometimes discontinuous (operated slice #1340). Screw holes also created discontinuities in the cortical boundary (operated slices #820, #1260, and #2340). For these reasons, two types of slices were excluded in the operated scans: 1) slices within the osteotomy region and 2) slices with a screw-hole discontinuity.

Boundary detection in the operated and contralateral limbs shared some principal steps, which are depicted by Figure 2 in an example slice image from an operated limb taken just proximal to the osteotomy. First, each image (I) was convolved with a square Gaussian kernel (G) of size $\mathit{n}$ = 31 with standard deviation of $\mathit{\sigma }$ = 10 to remove noise. Next, image binarization was carried out, with pixel values higher than a specified threshold ($\mathit{\tau }$) assigned a binary 1 and the remaining pixels assigned a binary 0. The pericortical, endocortical, and callus boundaries were detected from the binarized images and a medullary centroid was defined as the mean of the endocortical boundary. The resulting boundaries were in a Cartesian coordinate system with uneven point angular spacing. To create a uniform boundary discretisation and to define radial paths from the medullary centroid outward, we resampled the as-detected boundary in polar coordinates using a linear interpolation at every integer angle from 0 to 359. A Savitzky-Golay filter was then applied on the resampled points to produce a more physiologically realistic smooth boundary.

Figure 2. Representative example of the boundary detection procedure for a single slice of an operated limb, segmented using non-adaptive thresholding. The raw slice image (a) was de-noised by applying a Gaussian filter, with a structural similarity index measure (SSIM) of 0.98 with respect to the raw image (b). In this example, preliminary pericortical (c), endocortical (d), and callus (e) boundaries were detected after segmenting with thresholds ${\tau }^{\mathit{peri}}$ = 7,000, ${\tau }^{\mathit{endo}}$ = 5950, and ${\tau }^{\mathit{call}}$ = 4,250 that differentiated higher density bone from the combined region of all mineralised tissue including callus. The preliminary pericortical boundary (c) was not correct due to cortical remodelling, so a convexity constraint was enforced to produce a more physiological shape (f). Boundaries were then uniformly resampled (g) and smoothed with a savitsky-golay filter (h), which were subsequently cleaned up to remove redundant points (i). Wherever the corrected pericortical boundary was larger than the detected callus boundary, the pericortical radius was collapsed to the detected callus radius (j). In this example, no further change was needed because callus was thicker than 15 pixels along most radial lines (k).

Images with callus required additional steps for boundary detection. Simple threshold-based segmentation produced an as-detected pericortical boundary that could be locally concave, where the true mid-diaphyseal cortex cross-section should be convex. Therefore, a convexity constraint was imposed prior to interpolation to resample with regularly spaced points. The convexity constraint was not applied to images proximal to the most proximal screw or distal to the most distal screw because the native bone outer boundary can be naturally concave near its ends. The focus of this algorithm was on detecting remodelling changes in cortical bone, so in instances of complex outer boundary of callus, only the inner point of the callus boundary was preserved along each radial path line before converting boundaries back to Cartesian coordinates.

Boundaries needed final touch-ups in two circumstances. First, the imposition of the convexity constraint on the pericortical boundary followed by linear interpolation occasionally resulted in a pericortical boundary with a larger local radius than the calculated callus boundary. To correct this error, the pericortical boundary was collapsed onto the callus boundary wherever the former extended beyond the latter. Second, far from the osteotomy, a thin pericortical false callus was sometimes detected. To correct for this artefact, the pericortical boundary was dilated onto the callus boundary when callus thickness at all the 360 profile lines was smaller than 15 pixels (~0.9 mm).

Adaptive thresholding

A critical step of the boundary detection process described above is selection of appropriate thresholds ($\mathit{\tau }$) for differentiation between higher- and lower-density regions of tissue. In preliminary testing of our algorithm, we observed that a fixed threshold ($\mathit{\tau }$ = 4,250) could be used for identifying boundaries in slice images without callus, but this fixed threshold did not correctly detect the old pericortical boundary in many slices of the operated limbs. Callus density is heterogeneous at this stage of healing, with generally higher BMD at the proximal and distal ends and lower BMD near the osteotomy line. For this reason, it was impossible to identify one threshold that could reliably distinguish between callus and bone in all images. Hence, we introduced an algorithm to adaptively determine appropriate endocortical (${\tau }^{\mathit{endo}}$) and pericortical (${\tau }^{\mathit{peri}}$) thresholds on a slice-by-slice basis.

The presence of callus inside the cortex near the osteotomy was irregular, occasionally filling the entire medullary space with a porous structure. When the endocortical threshold (${\tau }^{\mathit{endo}}$) was not high enough, this internal callus was interpreted as cortical bone. The endocortical boundary detected in this manner would be non-physiologically tortuous and a poor estimate of the true endocortical boundary. Conversely, when the endocortical threshold (${\tau }^{\mathit{endo}}$) was too high, the inner wall of the cortex was identified as callus, leading to an overly tight but smooth boundary. This observation suggested that boundary tortuosity could be used to help identify the optimal endocortical segmentation threshold (${\tau }^{\mathit{endo}}$). The adaptive algorithm needed to increase the candidate endocortical threshold (${\tau }^{\mathit{endo}}$) from an initial value until an acceptable detection of endocortical boundary was achieved. To evaluate the acceptability of a detected endocortical boundary, we defined endocortical curve tortuosity ($\mathit{\zeta }\left[L_p^{-5}\right]$) as follows:

(2) $\mathit{\zeta }=\frac{{\int }_0^{2\pi }\mathit{\kappa }{\prime }^2\mathit{d\theta }}{{\int }_0^{2\pi }\mathit{rd\theta }}$(2)

where $\mathit{\kappa }\left[L_p^{-1}\right]$ is the local curvature and $\mathit{r}[L_p]$ is the interpolated radius, just before being smoothed by a Savitzky-Golay filter, both expressed in pixel units. Adaptive endocortical thresholding started from the lower bound of endocortical threshold (${\tau }_{min}^{\mathit{endo}}$ = 4,250) and increased in increments of 50 until either tortuosity was less than or equal to an acceptable value (${\zeta }_{\mathit{acp}}$, to be determined) or it reached an upper bound (${\tau }_{max}^{\mathit{endo}}$ = 6,750). Increasing the endocortical threshold beyond its upper bound (${\tau }_{max}^{\mathit{endo}}$) consistently produced overly tight boundaries, even though possibly smoother. In the presence of dense callus in the medullary canal, the best-detected endocortical boundary was not necessarily perfectly smooth, unlike other images in which the endocortical boundaries were perfectly smooth. To address this variability, the definition of maximum acceptable tortuosity (${\zeta }_{\mathit{acp}}$) needed to be context-sensitive and dependent on the endosteal callus density. Therefore, it was taken as a linear function of the local endocortical threshold (${\tau }^{\mathit{endo}}$) divided by a constant $\mathit{c}$, empirically calibrated as $2.5\times {10}^{-2}$. Use of a constant value for the maximum acceptable tortuosity (${\zeta }_{\mathit{acp}}$) did not stop the algorithm at an ideal endocortical threshold (${\tau }^{\mathit{endo}}$) in images with substantial callus and generated an overly tight endocortical boundary. The chosen linear function definition leniently allowed for a higher boundary tortuosity as endocortical threshold rose.

Tortuosity proved to be useful for adaptive detection of endocortical boundaries, but our initial investigations indicated that there was no consistent relationship between tortuosity and the quality of pericortical boundary detection in an operated limb. This suggested the need for a different criterion for adaptive thresholding of the pericortical boundary. Slice images before the first and after the last screws typically did not have any callus, so their pericortical boundaries were detected by the same threshold used in contralateral images (${\tau }_{\mathrm{\infty }}^{\mathit{peri}}$ = 4,250). Due to the absence of callus in these images, the detected boundaries with this threshold were physiological, so the convexity constraint was not imposed on them. In fact, imposing convexity would hurt the correctly detected boundaries because slice images in the proximal and distal regions were not necessarily convex in shape (e.g. operated slices #100, #2340 and contralateral slice #235, #2475 in Figure 1). For images in the central diaphysis, the proper pericortical threshold (${\tau }^{\mathit{peri}}$) was iteratively identified by calculating a cross-sectional area ratio ($k_A$) within each single slice image:

(3) $k_A=\frac{{\mathit{A}}_{\mathit{endo}}}{A_{\mathit{peri}}}$(3)

where the enclosed area of endocortical boundary ($A_{\mathit{endo}}\left[L_p^2\right]$) was divided by the enclosed area of the pericortical boundary ($A_{\mathit{peri}}\left[L_p^2\right]$). We assumed that the ideal pericortical threshold (${\tau }^{\mathit{peri}}$) would minimize the difference in calculated area ratios between axially registered slices in the intact ($k_A^{con}$) and operated ($k_A^{op}$) limbs. The algorithm increased the threshold (${\tau }^{\mathit{peri}}$) for pericortical boundary detection from a lower bound (${\tau }_{min}^{\mathit{peri}}$ = 7,000) iteratively in increments of 50 until either the area ratio in the operated image ($k_A^{op}$) equaled/exceeded that in the contralateral image ($k_A^{con}$) or until the pericortical threshold (${\tau }^{\mathit{peri}}$) reached the upper bound (${\tau }_{max}^{\mathit{peri}}$ = 7,600). The upper and lower bounds of the endocortical and pericortical threshold search ranges (${\tau }_{min}^{\mathit{endo}},{\tau }_{max}^{\mathit{endo}},{\tau }_{min}^{\mathit{peri}},{\tau }_{max}^{\mathit{peri}}$) were empirically determined by visually examining the results during algorithm development to determine the limits of segmentation quality. Table 1 contains the values selected for the constants.

Table 1. Thresholding parameters and their roles in the algorithm for adaptive boundary detection.

Boundary	Parameter	Value	Explanation
Endocortical	${\tau }_{min}^{\mathit{endo}}$	4250	Lower bound of threshold
	${\tau }_{max}^{\mathit{endo}}$	6750	Upper bound of threshold
	$\mathrm{\Delta \tau }$	50	Increment of threshold
	$\mathit{c}$	0.025	Calibration constant for tortuosity limit
Pericortical	${\tau }_{\mathrm{\infty }}$	4250	Threshold before/after first/last screw
	${\tau }_{min}^{\mathit{peri}}$	7000	Lower bound of threshold
	${\tau }_{max}^{\mathit{peri}}$	7600	Upper bound of threshold
	$\mathrm{\Delta \tau }$	50	Increment of threshold

BMD calibration and radial sampling for dimensionality reduction

After boundary detection, the final endocortical and cortical boundaries were defined by points spaced at every integer angle. Using these points, radial profile lines were used to sample bone mineral density (BMD). Hounsfield unit (HU) values in each pixel were calculated using a linear conversion from pixel value ($\mathit{v}$) according to the information provided in the DICOM file metadata:

(4) $\mathit{HU}=0.515\mathit{v}-1000$(4)

BMD ($\mathit{\rho }$) was obtained from a linear conversion of pixel HU values based on a calibration curve defined from a radiological phantom, which in this case was:

(5) $\mathit{\rho }=0.38010\mathit{HU}-7.37444$(5)

To measure circumferential changes, the radial line average BMD (${\bar{\rho }}_{\mathit{\theta }}$) was calculated along each profile line and represents BMD at the corresponding angle ($\mathit{\theta }$). Similarly, circumferential line average BMD (${\bar{\rho }}_{\mathit{r}}$) was computed at incremental normalized radii from the endocortical to pericortical boundary and represents BMD at the corresponding radius. Figure 3 illustrates the paths used for calculating these quantities. To enable global analysis, we extracted a representative density ($\mathit{P}$) from each slice image, defined as the maximum of circumferential line average BMD (${\bar{\rho }}_{\mathit{r}}$). The representative density ($\mathit{P}$) of each operated or contralateral slice image was divided by its maximum in the contralateral limb to obtain normalized representative density ($\bar{P}$). This normalized representative density ($\bar{P}$) was a function of axial position (z) and was interpreted to assess spatial remodeling effects in 1D (axial position relative to osteotomy), thereby achieving dimensionality reduction from the 3D imaging dataset. For subsequent analysis, the normalized representative density for contralateral scans (${\bar{P}}^{\mathit{con}}$) was smoothed by a seventh-degree polynomial (${\bar{P}}_{\mathit{s}}^{\mathit{con}}$), but representative density for operated images (${\bar{P}}^{\mathit{op}}$) was not smoothed to avoid obscuring local remodeling features. A slicewise remodeling index ($\mathit{R}$) as a function of axial position ($\mathit{z}$) was calculated by:

(6) $\mathit{R}\left(\mathit{z}\right)=\frac{{\bar{P}}_{\mathit{s}}^{\mathit{con}}\left(\mathit{z}\right)-{\bar{P}}^{\mathit{op}}\left(\mathit{z}\right)}{{\bar{P}}_{\mathit{s}}^{\mathit{con}}\left(\mathit{z}\right)}$(6)

Figure 3. Examples of paths over which circumferential line average (${\bar{\rho }}_{\mathit{r}}$), maximum circumferential line average ($\mathit{P}$), and radial line average (${\bar{\rho }}_{\mathit{\theta }}$) BMD were calculated.

Summary parameters defining the spatial variation of cortical remodelling activity were then extracted from the remodelling index data. The slices corresponding to the upper 10% of calculated $\mathit{R}\left(\mathit{z}\right)$ values defined a region of high remodelling activity. Maximum remodelling index ($\mathit{M}$) and its axial position ($z_M$) were defined as the average of remodelling index ($\mathit{R}$) and axial position ($\mathit{z}$) within the region of high remodelling activity. The height of this region ($\mathit{H}$) became another global characteristic.

Finally, the summed total remodelling activity for the entire scan ($\mathit{S}$) was calculated through numerical integration as follows:

(7) $\mathit{S}=\frac{{\int }_{{\mathit{z}}_{\mathit{min}}}^{z_{\mathit{max}}}{\bar{P}}_{\mathit{s}}^{\mathit{con}}\left(\mathit{z}\right)-{\bar{P}}^{\mathit{op}}\left(\mathit{z}\right)\mathit{dz}}{{\int }_{{\mathit{z}}_{\mathit{min}}}^{z_{\mathit{max}}}{\bar{P}}_{\mathit{s}}^{\mathit{con}}\left(\mathit{z}\right)\mathit{dz}}$(7)

where $z_{\mathit{max}}$ and $z_{\mathit{min}}$ are the axial positions of the most proximal and distal paired slices for the operated and contralateral scans in each animal.

Callus morphometric measurements

To enable a global analysis of the association between remodelling activity and callus production, the following morphometric measures were also obtained for each animal: total callus volume ($V_c$) [cm³], callus median density (${\rho }_c)$ [mgHA/cm³], and callus mineral composition ($V_c\times {\rho }_c$) [mgHA]. These morphometric measurements were calculated from callus segmentation described previously (Schwarzenberg et al. 2019, 2021).

Statistical analysis

All statistical analyses were performed in IBM SPSS Statistics (v.29; IBM Corp., Armonk, NY). To assess whether cortical remodelling (loss of BMD) was related to the quantity of new mineralised tissue in the callus, we performed Pearson’s correlations between all calculated remodelling parameters and all callus morphometric measurements, a total of 21 correlations.

Results

The algorithm successfully detected boundaries in a total of 16,496 operated and 31,168 contralateral slice images from 12 sheep. Table 2 reports the number of included and excluded slice images for each scan. All available slice images from contralateral limbs were included. A few slice images in operated limbs showed completely incorrect detection of boundaries (negative cortical thickness) and were excluded automatically. One sheep (L) had a region of slice images with distinctly high and spatially irregular remodelling, which were left out manually (this animal is discussed ahead).

Table 2. The number of slice images classified according to sheep and limb.

Sheep	Contralateral		Operated
	Included	Paired	Excluded	Included
A	2476	1309	0	1309
B	2477	1221	0	1221
C	2477	1289	0	1289
D	2477	1064	3	1061
E	2671	1490	0	1490
F	2671	1447	8	1439
G	2671	1518	0	1518
H	2639	1606	0	1606
I	2631	1546	16	1530
J	2670	1324	0	1324
K	2670	1453	0	1453
L	2638	1256	0	1256
Total	31168	16523	27	16496

Representative boundary detection and data sampling results are presented in Figure 4 (example with substantial callus present) and Figure 5 (example with little callus present). Circumferential variations of radial line average BMD (${\bar{\rho }}_{\mathit{\theta }}$) were generally minimal. Notably, in the operated limb of Figure 5, slices close to screw holes showed loss of BMD close to the screws. The presence of the NAC hole also lowered the radial line average BMD (${\bar{\rho }}_{\mathit{\theta }}$) over a small angular span in a limited number of slices. These effects aside, circumferential variations in BMD were largely negligible and consistent between operated and contralateral limbs, so the subsequent analysis was focused on radial and axial variations.

Figure 4. Comparison of paired slices from operated and contralateral limbs located just proximal to the osteotomy where substantial callus is present. a/b) images after boundary detection and resampling in 1-degree circumferential increments. c/d) radial line average BMD (${\bar{\rho }}_{\mathit{\theta }}$) superimposed over the cortical cross-section shape in each image for reference. e) point cloud of BMD (ρ) distribution in each cross section. Representative densities for operated ($P^{\mathit{op}}$) and contralateral ($P^{\mathit{co}}$) slice images are indicated by the large dots. Considerable cortical remodeling was measured for this slice ($\mathit{R}$ = 10%).

Figure 5. Comparison of paired slices from operated and contralateral limbs located ~68 mm proximal to the centre of osteotomy where little callus is present. a/b) images after boundary detection and resampling in 1-degree circumferential increments. c/d) radial line average BMD (${\bar{\rho }}_{\mathit{\theta }}$) superimposed over the cortical cross-section shape in each image for reference. e) point cloud of BMD (ρ) distribution in each cross section. Representative densities for operated ($P^{\mathit{op}}$) and contralateral ($P^{\mathit{co}}$) slice images are indicated by the large dots. Moderate cortical remodeling was measured for this slice ($\mathit{R}$ = 5%).

Radial-direction cortical remodelling depended on the axial position ($\mathit{z}$) relative to the osteotomy, with the most substantial cortical changes occurring close to the osteotomy. To compare remodeling patterns between individual slices at different axial positions, point clouds containing BMD values for all 360 radial paths were plotted using a normalized radial coordinate ($\bar{r}$) that takes a value of 0 at the pericortical boundary and −1 and the endocortical boundary. Collapsing the radial data in this way revealed the higher remodeling activity close to the osteotomy (Figure 4) compared to far from the osteotomy (Figure 5). The large red dots in Figures 4e and 5e illustrate representative density in each operated image ($P^{\mathit{op}}$), while the large black dots illustrate representative density in the matched contralateral slice ($P^{\mathit{con}}$). These quantities represent BMD ($\mathit{\rho }$) in the paired images concisely and combine to produce the remodeling index ($\mathit{R}$) at this axial location ($\mathit{z}$). Closer to the osteotomy, the representative density marker (large red dot) is farther from the pericortical surface and lower in magnitude relative to the intact representative density (large black dot) (Figure 4e versus Figure 5e), indicating greater radial erosion of cortical BMD relative to the contralateral bone.

To enable a global analysis of the differences between limbs and the relationship between axial position and remodelling activity, the slice-by-slice normalised representative densities ($\bar{P}$) and remodeling indices ($\mathit{R}$) were plotted with respect to the axial position of the slice images. A representative example for sheep (A) is provided in Figure 6. The intact limb cortical density was highest and relatively constant across the midsection of the diaphysis, a distance of approximately 70 mm, and dropped slightly to a normalized value of ~ 0.94 in both proximal and distal ends. In contrast, the density variations in the operated limb were complex. Furthest from the osteotomy, density was nearly equivalent between the operated and contralateral limbs, resulting in a remodeling index ($\mathit{R}$) that approached zero. Closer to the osteotomy, remodeling activity increased substantially, peaking close to the osteotomy. For this example animal (A), the remodeling zone height was $\mathit{H}$ = 19.0 mm, peaking at $\mathit{M}$ = 9.39% at $z_m$ = −8.3 mm, on the distal side. The summed total remodeling activity for this animal was $\mathit{S}$ = 3.24%. BMD also exhibited variability between screw holes, generally following a global trend of increasing remodeling activity closest to the osteotomy, but reflecting a localised loss of mineral near screws. Global comparison of limbs in the other 11 animals in this study are shown in the Supplemental Digital Content: Appendix.

Figure 6. Global differences between the operated and contralateral limb of one animal were evaluated based on the changes of a) normalised representative densities ($\bar{P}$) in each limb and b) remodelling index ($\mathit{R}$) with characteristic measures obtained from its scatter.

Calculated global remodelling parameters for all animals as well as their averages are summarised in Table 3. The average remodelling zone height ($\mathit{H}$) was 31.5 mm (SD = 8.9) and accounted for ~ 21% of the diaphyseal segment (150 mm). The average maximum remodeling index ($\mathit{M}$) was 11.2% (SD = 1.6) occurring at a 4.4 mm distal to the center of osteotomy (SD = 4.9). Summed total remodeling activity ($\mathit{S}$) was 4.6% on average (SD = 1.1). Morphometric parameters for each callus are summarized in Table 4. The average total callus volume, callus median density, and callus mineral composition were 11.2 cm³ (SD = 6.2), 697 mgHA/cm³ (SD = 84.0), and 7,725 mgHA (SD = 4,024), respectively.

Table 3. Summary of the characteristic measures obtained from the global assessment of remodelling for the sheep in this study.

Animal ID	Remodeling Zone Height $\mathit{H}$ [mm]	Maximum Remodeling Index $\mathit{M}$ [%]	Axial Position of Maximum Remodeling ${\mathit{z}}_{\mathit{M}}$ [mm]	Summed Total Remodeling Activity $\mathit{S}$ [%]
A	19.0	9.4	−8.3	3.2
B	36.0	13.5	−4.5	5.9
C	25.0	12.3	−9.3	4.6
D	36.0	13.0	2.8	6.2
E	48.0	11.3	−2.4	4.6
F	41.0	13.4	−14.6	5.5
G	23.0	9.1	−0.9	3.3
H	31.0	11.1	−2.8	3.2
I	36.0	10.8	−4.9	5.0
J	29.0	10.2	0.9	5.2
K	22.0	9.3	−4.9	3.6
Average	31.5	11.2	−4.4	4.6
Std. Dev.	8.9	1.6	4.9	1.1

Table 4. Summary of the callus morphometric measures for the sheep in this study.

Animal ID	Total Callus Volume ${\mathit{V}}_{\mathit{c}}$ [cm^3]	Callus Median Density ${\mathit{\rho }}_{\mathit{c}}$ [mgHA/cm^3]	Callus Mineral Composition ${\mathit{V}}_{\mathit{c}}\times {\mathit{\rho }}_{\mathit{c}}$ [mgHA]
A	7.7	643	4,930
B	16.0	724	11,576
C	5.6	646	3,631
D	23.1	570	13,170
E	1.5	636	939
F	13.0	695	9,047
G	6.6	682	4,475
H	6.3	842	5,341
I	16.3	656	10,698
J	14.3	838	11,936
K	12.5	736	9,235
Average	11.2	697	7,725
Std. Dev.	6.2	84.0	4,024

Pearson’s correlation coefficient (r) and the corresponding p-value is provided in Table 5 for each combination of remodelling (Table 3) and callus morphometry (Table 4) parameters. All the remodelling parameters except axial position of maximum remodelling ($z_M$) were highly and significantly correlated with each other ($\mathit{r}$ >0.61, $\mathit{p}$ <0.05). Among the morphometric parameters, total callus volume ($V_c$) and callus mineral composition ($V_c\times {\rho }_c$) were highly and significantly correlated ($\mathit{r}$ = 0.961, $\mathit{p}$ >0.001), which is expected because total mineral composition is the product of callus volume and mineral density. The strongest correlation between the two groups of independent parameters (remodeling and morphometric measures) was between summed total remodeling activity ($\mathit{S}$) and total callus volume ($V_c$) ($\mathit{r}$ = 0.687, $\mathit{p}$ = 0.020), followed by the correlation between summed total remodeling activity ($\mathit{S}$) and callus mineral composition ($V_c\times {\rho }_c$) ($\mathit{r}$ = 0.644, $\mathit{p}$ = 0.033).

Table 5. Pearson’s correlation coefficients and the corresponding p-values (italicised) for all possible pairs of remodelling and callus morphometry parameters.

	$\mathit{M}$	${\mathit{z}}_{\mathit{M}}$	$\mathit{H}$	${\mathit{V}}_{\mathit{c}}$	${\mathit{\rho }}_{\mathit{c}}$	${\mathit{V}}_{\mathit{c}}\times {\mathit{\rho }}_{\mathit{c}}$
Summed Total Remodeling Activity ($\mathit{S}$ [%])	0.793**	0.080	0.618*	0.687*	−0.268	0.644*
	0.004	0.815	0.043	0.020	0.425	0.033
Maximum Remodeling Index $\mathit{M}$ [%]		−0.241	0.636*	0.372	−0.230	0.306
		0.475	0.035	0.259	0.496	0.361
Axial Position of Maximum Remodeling ${\mathit{z}}_{\mathit{M}}$ [mm]			0.017	0.254	0.089	0.240
			0.961	0.452	0.794	0.477
Remodeling Zone Height $\mathit{H}$ [mm]				0.102	−0.177	0.060
				0.764	0.602	0.862
Total Callus Volume ${\mathit{V}}_{\mathit{c}}$ [cm^3]					−0.130	0.961**
					0.702	<0.001
Callus Median Density ${\mathit{\rho }}_{\mathit{c}}$ [mgHA/cm^3]						0.124
						0.717
Callus Mineral Composition ${\mathit{V}}_{\mathit{c}}\times {\mathit{\rho }}_{\mathit{c}}$ [mgHA]

*Correlation is significant at the 0.05 level (2-tailed).

**Correlation is significant at the 0.01 level (2-tailed).

Discussion

The newly developed adaptive boundary detection algorithm successfully detected the endocortical and pericortical boundaries of intact and operated ovine tibiae and enabled calculation of summary metrics to characterise the spatial variation of cortical remodelling during fracture repair. A key observation from the development of this procedure was the need for adaptive thresholding to detect regions corresponding to old cortical bone. While fixed thresholding successfully detected boundaries in contralateral images, operated images had wide variations of bone and callus objects within and between sheep, necessitating a flexible method of boundary detection.

Figure 7 illustrates the importance of choosing the endocortical threshold (${\tau }^{\mathit{endo}}$) adaptively by showing two example images with different cross sections in an operated limb. In both cases, the adaptively-determined endocortical threshold (${\tau }^{\mathit{endo}}$) for each image detected the endocortical boundary accurately, while the swapped endocortical threshold (${\tau }^{\mathit{endo}}$) from the other image failed. Similarly, Figure 8 highlights the importance of adaptive determination of pericortical threshold (${\tau }^{\mathit{peri}}$) using the same images as in Figure 7 and their correct endocortical boundaries. Again, adaptive thresholding successfully detected the pericortical boundary, but swapping the pericortical thresholds (${\tau }^{\mathit{peri}}$) between images within this animal failed.

Figure 7. Successful detection of endocortical boundaries in different slice images within one operated limb required adaptive endocortical thresholding (${\tau }^{\mathit{endo}}$). a/b) for a slice located ~11 mm distal to the osteotomy centre with no endosteal callus, the endocortical boundary was properly detected with ${\tau }^{\mathit{endo}}$ = 4,250, but improperly detected with ${\tau }^{\mathit{endo}}$ = 6,650 (correct threshold in d). c/d) for a slice located 4 mm proximal to the osteotomy centre with endosteal callus, the endocortical boundary was incorrectly detected with ${\tau }^{\mathit{endo}}$ = 4,250 (correct threshold in a), but properly detected with ${\tau }^{\mathit{endo}}$ = 6,650.

Figure 8. Detection of pericortical boundaries in the same slice images of figure 7. a/b) pericortical boundary was properly detected with ${\tau }^{\mathit{peri}}$ = 7,600, but was improperly dilated with ${\tau }^{\mathit{peri}}$ = 7,000 (correct threshold in d). c/d) pericortical boundary was severely eroded with ${\tau }^{\mathit{peri}}$ = 7,600 (correct threshold in a), but correctly follows pericortical wall with ${\tau }^{\mathit{peri}}$ = 7,000.

In all intact tibiae, our results did not show any substantial circumferential variations of radial line average BMD (${\rho }_{\theta }$) (examples: Figures 4 and 5). This finding is in contrast with recent reports of spatial variations in femoral cortical bone of 4-month-old sheep (Manandhar et al. 2023) and of foals younger than 1 year old (Moshage et al. 2020). The difference might be rooted in the use of skeletally mature rather than juvenile animals in this study.

In the fractured tibiae, we found widespread cortical remodelling consistent with findings from prior image analysis (Ren et al. 2022) and nanoindentation measurements (Preininger et al. 2011). Cortical remodelling was substantial close to the osteotomy and almost negligible near the ends of the operated bone (Figure 6). Since all the animals had mid-diaphyseal osteotomies, our results cannot reveal how the location of the osteotomy may affect the remodelling pattern. Notably, the measure of remodelling here is different from that in (Ren et al. 2022) due our inclusion of the full-length bone image set, so the numbers reported cannot be directly compared. However, our large region of interest was useful because it revealed both global effects of remodelling related to the injury repair, and confirmed the presence of localised effects near screws that may be attributable to microdamage (G. Wang et al. 2014a; L. Wang et al. 2014b; Yu et al. 2014). BMD changes we observed near screws could be related to the pre-strain existing around bone screws due to the mechanical interaction of screw threads and bone, which is intensified upon loading the bone (MacLeod et al. 2012). This mechanical environment might have caused microdamage in the vicinity of the screw, triggering the remodelling patterns observed in our data between screw holes.

Formation of callus requires a source of minerals. We identified a significant correlation between the global remodelling score and the total mineral composition of the callus. This finding supports the hypothesis that rapid anabolic activity to form the callus is enabled in part by catabolic depletion over a large region of the cortical bone. In reviewing the analyzed images to look for this effect, we discovered an interesting pattern of remodelling in some animals. While cortical bone was expected to lose its mineral content near its walls to supply callus formation, some animals exhibited a noteworthy drop in BMD within the cortical core, rather than near the walls. This pattern of remodelling was evident in multiple animals but was extremely pronounced in one animal (Sheep L) and is shown in Figure 9. The presence of this ring-like cortical remodelling effect produced spurious pericortical segmentations in this animal close to the osteotomy because the outer cortex was too extensively remodelled to be detectable. Consequently, Table 3 lacks the remodelling measures of this animal and it was left out of the statistical analysis.

Figure 9. Noteworthy pattern of extreme remodelling in the sheep L. Remodelling at the outer wall of cortical bone (slice #1500) formed a ring of low BMD at the outer wall (#1638) extending inward to the core (#1783) and reaching the inner wall partially (#1883) before occupying an irregular region of cortex (#2032 and #2132).

This study has limitations that suggest promising future directions for continuing research. First, slices with screw hole defects were excluded, which may forfeit some useful insights about remodelling near screws. The number of excluded slices also varied slightly due to minor variations in the alignment of the plates and screws for each animal at the time of surgery. While the slicewise remodelling index ($\mathit{R}$) was not affected by exclusion of slices, the summed total remodelling activity ($\mathit{S}$) is a global measure and it could be more sensitive to slice exclusion if comparing results between animals with different implant types and numbers of screws. Second, we registered the operated and contralateral image sets using the NAC and used the intact cortical area ratio to adapt the segmentation thresholds in the operated limbs. This assumption was reasonable, but the natural asymmetry existing between the limbs potentially introduces errors in the detection of periosteum. This method is also limited in that periocortical boundary detection in the operated limb is tied to that in the contralateral limb, so both must be excised and scanned.

Future alternatives to address these limitations include active contour models that start with an initial guess of the boundary and deform until fitting an edge (Kass et al. 1988). Such a model can outline the entire cortex bifurcated by a screw (Williams and Shah 1992), unlike what happens in thresholding. Further, it can be constrained to produce physiological shapes similar to a database of boundaries (Cootes and Taylor 1992). A larger database of intact bone scans could also enable development of statistical shape models (Patil et al. 2023) for cortical boundary mapping in the presence of callus, even in cases of extreme irregular remodelling where threshold-based methods are unsuitable (e.g. Sheep L of this study).

Finally, it is important to note that the BMD-derived measures of cortical bone change reported here are surrogate measures of the net effects of remodelling activity (coupled bone resorption and formation) at 9 weeks after surgery. The voxelized BMD measurement itself is also a surrogate measure of cortical microstructure (e.g. porosity) and mineral composition. In the future, higher resolution and time-series imaging may be helpful for direct assessment of coupled remodelling at the microstructural level.

Conclusion

In this study, cortical bone density changes during ovine fracture healing were quantified using an adaptive boundary detection procedure over the entire diaphyseal segment, a much larger region of interest than in other prior studies. Global effects of cortical resorption related to the injury repair occurred concurrently with the localised effects near screws related to microdamage. The findings suggested that cortical adaptation occurs over an extensive region of bone and that cortical bone may be serving as a local source of the mineral needed to build callus. Future applications of this work include evaluation of end-stage fracture repair to assess the temporal pattern of fracture line consolidation and cortical BMD recovery, processes that are important for long-term recovery of bone strength after early stiffness has been achieved through callus bridging.

Acknowledgements

Research reported in this publication was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the National Institutes of Health under Award Number R21AR081435 and by the National Science Foundation under Grant Number CMMI-1943287. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or National Science Foundation. The preclinical studies from which data was obtained were funded by the Johnson & Johnson Family of Companies. The authors wish to thank Beat Lechmann (Johnson & Johnson Family of Companies) for agreeing to grant access to the ovine study data for our analyses. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the Johnson & Johnson Family of Companies.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/21681163.2024.2345165

Additional information

Funding

The work was supported by the DePuy Synthes; National Science Foundation [CMMI-1943287]; National Institutes of Health – U.S. [R21AR081435]; National Science Foundation - U.S. [CMMI-1943287].

Notes on contributors

Alireza Ariyanfar

Alireza Ariyanfar Currently a PhD student in Mechanical Engineering at Lehigh University, he studies bone fracture healing by processing medical images. He completed an MSc and a BSc in the same major in Iran at Sharif University of Technology and the University of Science and Technology, respectively.

Karina Klein

Karina Klein PhD is an experienced clinician and veterinary orthopedic surgeon, co-director of the Musculoskeletal Research Unit (MSRU, University of Zürich, Switzerland). She has directed and performed many studies and surgeries for musculoskeletal regeneration, in compliance with Good Laboratory Practice. Dr. Klein has extensive experience in developing animal models to answer specific research questions in a wide range of fields (spine, bone augmentation, tracheal stents, various prostheses).

Brigitte von Rechenberg

Brigitte von Rechenberg The founder and former head of the Musculoskeletal Research Unit at the University of Zurich in Switzerland, Prof. em. Dr. med. vet. Brigitte von Rechenberg is a board-certified veterinary surgeon (Dipl. ECVS) specialized in small animal surgery and in trauma and orthopedics. She has extensive experience in clinical and experimental animal surgery. Her expertise spans various orthopedic applications, especially bone and cartilage. She has, since 2020, officially retired from the University of Zürich, and has formed her own consulting company (Brigitte von Rechenberg Consulting GmbH), where she still shares her expertise as a consultant for the Musculoskeletal Research Unit.

Salim Darwiche

Salim Darwiche A bioengineer by training, Dr. Darwiche is the co-director of the Musculoskeletal Research Unit (MSRU, University of Zürich, Switzerland), a Good Laboratory Practice certified, pre-clinical in vivo studies academic research group. He is an experienced study director, with ten years of expertise in designing, implementing and managing pre-clinical in-vivo studies in musculoskeletal tissue regeneration (bone, cartilage, tendon and skin).

Hannah L. Dailey

Hannah L. Dailey Associate Professor in the department of Mechanical Engineering & Mechanics at Lehigh University. Her research interests include computational biomechanics of fracture fixation and bone healing, fracture nonunion, and mechanics of musculoskeletal tissues. Her group emphasizes imaging-driven engineering approaches to clinical problems in orthopaedics and currently collaborates with veterinarians and surgeons in hospitals worldwide. Prof. Dailey also serves as co-founder and Chief Scientific Officer of OrthoXel, DAC, an Irish-based orthopaedic device company.