Novel strategies for the characterization of cancellous bone morphology: Virtual isolation and analysis

Abstract Objectives The advent of micro‐computed tomography (μCT) made cancellous bone more accessible than ever before. Nevertheless, the characterization of cancellous bone is made difficult by its inherent complexity and the difficulties in defining homology across datasets. Here we propose novel virtual methodological approaches to overcome those issues and complement existing methods. Materials and methods We present a protocol for the isolation of the whole cancellous region within a μCT scanned bone. This method overcomes the subsampling issues and allows studying cancellous bone as a single unit. We test the protocol on a set of primate bones. In addition, we describe a set of morphological indices calculated on the topological skeleton of the cancellous bone: node density, node connectivity, trabecular angle, trabecular tortuosity, and fractal dimension. The usage of the indices is shown on a small comparative sample of primate femoral heads. Results The isolation protocol proves reliable in isolating cancellous structures from several different bones, regardless of their shape. The indices seem to detect some functional differences, although further testing on comparative samples is needed to clarify their potential for the study of cancellous architecture. Conclusions The approaches presented overcome some of the difficulties of trabecular bone studies. The methods presented here represent an alternative or supporting method to the existing tools available to address the biomechanics of cancellous bone.

Cancellous bone is composed of intertwined trabecular elements forming an intricate structure with no analogy to regular solid shapes and, therefore, hard to describe through traditional morphometrics (Hildebrand et al., 1999;Odgaard, 1997). Furthermore, the developing trabecular lattice is modeled throughout life by the complex interaction between endogenous and exogenous factors (Cooper, 1990;Little et al., 2011). Given the multifaceted nature of its development, the characterization of cancellous bone should consider complexity as a feature useful to understand its variability and function.
Previous work focused on subsamples of trabecular regions to analyze the cancellous bone (e.g., Fajardo & Müller, 2001;Moon et al., 2004;Räth et al., 2008;Ryan & Ketcham, 2002). These subsamples bear only local information (Fajardo & Müller, 2001;Georgiou et al., 2019;Kivell et al., 2011) and their position and orientation are difficult to define homologously, especially across species with different external bone morphology (Georgiou et al., 2019). Studies have also observed that some measurements on cancellous bone depend on the size of the subsample  and that cancellous properties are not homogeneous across the whole epiphysis (Sylvester & Terhune, 2017). Because of these limitations, an approach focusing on the whole cancellous region could provide additional functional information. Indeed, prior research has shown that extended cancellous areas within an epiphysis can lead to more meaningful statistical comparisons (Sylvester & Terhune, 2017). In addition, characterizing the cancellous properties across the whole epiphysis allows mapping the differences across the bone and articular surfaces. This approach has been successfully used in previous research (Georgiou et al., 2019;Sylvester & Terhune, 2017;Tsegai et al., 2013;Tsegai et al., 2018b).
In summary, the study of cancellous bone can be made difficult by (I) its inherent complexity and (II) the lack of clear homology between cancellous regions in interspecific datasets. In this work we address both issues. First, we present a reproducible protocol for isolating cancellous structures from stacks of μCT images, thus allowing the analysis of cancellous regions without subsampling. Second, we present a set of indices measured on the topology of the cancellous structure that quantify new aspects of cancellous complexity.

| MATERIALS AND METHODS
Most of the operations used in this work are currently available in the R environment (R Core Team, 2020). The isolation of cancellous bone and the calculation of complexity indices are available in the package "indianaBones," currently available on Github (github.com/ AlessioVeneziano/indianaBones). A working example is provided in the Supplementary information. Following, we describe the general workflow making use of the R environment whenever possible but also highlighting the steps performed in the Amira/Avizo software (FEI Visualization). Figure 1 shows a diagram of the workflow.
2.1 | Orientate, crop the μCT stack, and define axes The μCT stack is cropped and orientated to facilitate further processing. This can be easily performed in software packages designed for handling image stacks. In Amira/Avizo, cropping can be performed in the Volume editor and orientation of the image stack is performed using the obliqueSlice function. Once cropped and orientated, the μCT stack can be exported in several file formats (e.g., TIFF F I G U R E 1 The workflow of the methodological approach introduced in this work. Actions are divided based on the environment where they are performed (R or Amira/Avizo) stack, DICOM, NIfTI). For the later calculation of some of the complexity indices, measurements are needed to provide scaling factors or reference axes for orientation. These measurements and axes can be defined in any software for image stack processing by collecting landmarks. In Amira/Avizo this can be performed using the Landmark editor: the recorded landmarks can be then exported in ASCII format and imported into R using the read.amira.set function in the "Arothron" R package (Profico et al., 2020).

| Import the μCT stack and binarize
In R, the μCT stack can be imported in several ways. TIFF images can be imported using the R package "tiff" (Urbanek, 2020) and the function readTiffStack available in the package "indianaBones" extends this to image stacks. DICOM and NIfTI stacks can be imported using the R package "oro.dicom" and "oro.nifti" (Whitcher et al., 2011), respectively. The stack is imported as a 3D array of dimensions K x M x N, where K and M are the image width and height in pixels, and N is the number of images in the stack. After importing, the binarization takes place, returning images with the bone and the background only (respectively as white and black pixels). This transformation is necessary because the protocol for the isolation of cancellous bone works on binary images only. Image binarization can be achieved via any suitable image segmentation method (Pham et al., 2000). One method widely used is the Otsu thresholding (Otsu, 1979;Vala & Baxi, 2013), which is available in the R package "EBImage" (Pau et al., 2010). The binarization can be performed in other software packages designed for image stack processing.

| Isolation of the cancellous bone
To isolate the cancellous from the cortical bone, the binarized stack is processed using the function splitBone in the R package "indianaBones." This function implements a protocol using a combination of the image processing operators "Dilation" and "Erosion" (Serra, 1982;Urbach & Wilkinson, 2007), which respectively enlarge and shrink features in binary images according to a pattern specified by a "structuring element." Dilation and erosion are performed using the package "EBImage," called by the "indianaBones" function splitBone. Structuring elements are matrices of odd dimensions that identify the pixel in the image being processed (Urbach & Wilkinson, 2007). Each pixel in the image is the center of the structuring element. The neighboring pixels are the cells of the element that surround its center. Each pixel is modified based on the value of its surrounding pixels, according to the pattern in the structuring element. In the case of dilation and erosion, each white pixel (region of interest) in the image will grow and shrink over its neighboring pixels in the fashion specified by the values in the structuring element. This process is illustrated in the supplementary information ( Figure S1).
Structuring elements can be defined using the makeBrush function in the "EBImage" R package.
The protocol consists of five sequential operations alternating dilation/erosion, applied iteratively, to subtractions between images and it is applied sequentially along the Z direction of the μCT stack.

| Skeletonization
The result of the isolation protocol is an image stack that is exported as a NIfTI 3D volume using the "oro.nifti" R package (although it could be exported in any format read by Amira/Avizo). The NIfTI format can be then imported inside the Amira/Avizo software, where the skeletonization takes place. Skeletonization returns the minimal geometric descriptor of an image, usually referred to as the "topological skeleton," by reducing it to a set of connected nodes and branches (Zhou & Toga, 1999) ( Figure S2). In the cancellous bone, branches represent the trabeculae while nodes are the points of connecting contiguous trabeculae. Currently, no suitable skeletonization procedure is available in the R environment, while the Amira/Avizo software provides a topological skeleton through the Auto Skeleton tool. In the isolated cancellous bone, this tool calculates a distance map of the trabeculae followed by thinning. The topological skeleton is returned as a list of node and branch coordinates, and indices of the connections between them. The skeleton is then exported from Amira/Avizo in ASCII format.

| Indices of cancellous complexity
The topological skeleton in ASCII format is read into the R environment using the readAmiraSkeleton function in the "indianaBones" R package. The nodes, branches, and connections of the skeleton are processed in R for calculating five indices: node density, trabecular angle, trabecular connectivity, trabecular tortuosity, and fractal dimension. A 2D representation of these measurements is illustrated in Node density is the 3D spatial density of skeleton nodes which is a proxy for trabecular spatial density and the relative proximity of trabecular connections. Node density has been previously used to address bone response to osteoporosis (e.g., Chappard et al., 1988).
The link between node density and function is straightforward: higher stress is counteracted by higher density of connections between trabeculae. Here, it is measured using a kernel density approximation (Venables & Ripley, 2002) over a regular 3D grid and it is expressed as number of nodes per cm 3 . The density grid can be plotted in 3D appreciate it visually within individuals. To reduce the effect of size on the calculation of spatial node density, the 3D coordinates of the skeleton nodes can be scaled using a measurement homologous across the dataset. The calculation of node density is performed using the skelDensity function in the "indianaBones" R package.
Trabecular angle is the 3D angle in degrees between a reference axis and the unitary resultant of all trabecular directions obtained by vector sum in 3D. The idea is that the main direction of trabeculae could detect the trajectory along which the mechanical load is dispersed (Hayes & Snyder, 1981). The direction of single trabeculae is calculated as the difference between the starting and ending nodes of each branch. The reference axis has to be homologous across the sample analyzed and this can be achieved by collecting anatomical or geometric landmarks that define the starting and ending point of the axis. The trabecular angle is calculated using the skelDirection function in the "indianaBones" R package. This function also calculates a major axis if a reference axis is not supplied. Trabecular connectivity has been measured via multiple approaches (Ding et al., 2002;Kabel et al., 1999;Odgaard & Gundersen, 1993). Here we define it as the average number of branches connected to each node of the topological skeleton. Only nodes with at least two connections (nonterminal nodes) are considered to calculate the average. Higher connectivity can be expected when cancellous structures are subject to large loads because more connections and more trabeculae allow to spread the load over a wider surface, thus releasing stress on localized areas (Silva & Gibson, 1997). The calculation of trabecular connectivity is performed using the skelConnectivity function in the "indianaBones" R package. Bayarri, 2015). Therefore, tortuosity reflects flexibility when the bone is subject to load. Tortuosity is the ratio between the arc length of a branch and the linear distance between its starting and ending nodes (Roque & Alberich-Bayarri, 2015). Tortuosity runs from one to infinity, where one corresponds to a straight branch (branch length equals node distance); as the branch gets more convoluted (and its length increases), tortuosity rises and tends to infinity. Infinity represents the absence of a theoretical maximum because the branch has no geometric limit to its length or convolution. Nevertheless, trabeculae are limited to an unknown maximum by the surrounding space and biomechanical requirements, so that infinity is a purely theoretical condition. Because tortuosity is the ratio between two lengths, it is dimensionless. The trabecular tortuosity is calculated using the skelTortuosity function in the "indianaBones" R package.
Fractal dimension is an index of complexity. It measures the change in detail over different scales of observation (Falconer, 2004).
Fractal dimension measured on μCT images or radiographs has been previously applied to the study of cancellous bone (Fazzalari & Parkinson, 1997;Feltrin et al., 2004;Haire et al., 1998;Messent et al., 2005). The rationale is that more complex cancellous structures are more interconnected, which allows spreading the load over a wider surface (Silva & Gibson, 1997). Fractal dimension is here measured on the 3D coordinates of the skeleton branches using the boxcounting algorithm (Annadhason, 2012). In this approach, 3D grids of decreasing cell size (decreasing cell side length, increasing number of cells) are superimposed over the cancellous skeleton. The number of cells overlapping the structure are counted for each subsequent grid: the fractal dimension is the slope of the line fitting the number of cells that overlap the skeleton versus the inverse of the cell size.
Fractal dimension is not expressed in units because it measures a fractional dimension (e.g., a fractal dimension of 1.5 is halfway from the one-dimension of a line and the two dimensions of a square). To reduce size effects, the 3D coordinates of the skeleton branches can be scaled using a measurement homologous across the dataset. The fractal dimension is calculated using the est.boxcount function in the "Rdimtools" R package (You, 2020).

| Application of the isolation protocol and complexity indices
To show the results of the protocol for cancellous isolation, we use μCT scans of five skeletal regions from five species of primates: the mandibular condyle, the brow ridge, the humerus, the femur and the fibula. Additional details about the specimens are reported in the supplementary information (Table S1). Prior to the application of the protocol, the μCT stacks is binarized by Otsu thresholding using the otsu function in the "EBImage" R package. A circular structuring element of 5 × 5 pixel size is used for the dilation/erosion operators. The number of dilation and erosions varied at each step and across bones but never exceeded six.
The usage of complexity indices is shown for a small comparative sample of μCT scanned femoral heads of specimens belonging to seven species of catarrhine primates. Additional details are reported in the supplementary information (Table S2). The aim is to demonstrate the usage, feasibility, and interpretation of the indices in comparative analyses and functional frameworks. Each femoral head is processed using the protocol described above for isolation of cancellous bone and calculation of the indices. Binarization is performed using Otsu thresholding in R (package "EBImage"). The segmented cancellous regions undergo skeletonization using the Amira 5.4.5 software package (FEI Visualization). The indices are then measured in the R environment using the "indianaBones" package. The trabecular angle is calculated with reference to the mediolateral axis, defined as the major axis passing through the femoral neck. Before calculating node density and fractal dimension, the skeleton is scaled on the height of the femoral head. The mediolateral axis and the femoral height are defined in Amira 5.4.5 using the "Landmark" editor. Node density is illustrated using color maps over a 100 × 100 × 100 3D grid used to estimate the kernel density: the colors represent differences in node per cm 3 , increasing from blue to red. F I G U R E 3 Graphical intuition of the indices measured on the topological skeleton of cancellous bone. For ease of visualization, the indices are shown for a 2D topological skeleton. Node density is represented by the number of nodes per unit area and it is calculated using a kernel density approximation over a discretized space. The trabecular angle (degrees) is measured between a reference axis (not shown) and the unitary resultant (red, double-headed arrow) of all trabecular directions (blue, double-headed arrows) obtained by vector sum. Connectivity is the mean number of branches connected to nonterminal nodes. Tortuosity is the ratio between the arc length of a branch and the linear distance between its starting and ending nodes (a/b). Fractal dimension is an index of complexity measured on the coordinates of the skeleton using the box-counting algorithm. In this approach, discrete regular grids of decreasing cell size are superimposed over the cancellous skeleton and the number of cells occupied by the skeleton are counted for each grid. Fractal dimension is the slope of the line fitting the number of cells that overlap the skeleton versus the inverse of the cell size 3 | RESULTS

| Cancellous bone isolation
In all the bones tested, the isolation protocol returns the cancellous lattice with, at most, only few small areas of the compact bone left attached ( Figure S3). The 2D and 3D results for the mandibular condyle are shown in Figure 2, where this region is used to present the steps of the protocol. For the other skeletal regions, the results are shown in Figure 4.

| Application of the complexity indices
Summary statistics of the complexity indices and their SD for the specimens analyzed are detailed in Table 1  The trabecular angle measured on the femoral head was referenced onto the mediolateral axis. Figure 6 shows the resultant direction of the trabeculae for each specimen and the angles are reported in Table 1. All angles are oriented mediolaterally with only minor departures from the reference axis. Average trabecular connectivity is larger in humans and Pan (3.86 ± 1.29 SD and 3.81 ± 1.19 SD, respectively) than in the other taxa (Table 1). The femoral head of Gorilla shows an average of 3.53 ± 0.95 SD branches per node, followed by Papio (3.49 ± 0.86 SD), and Macaca (3.38 ± 0.78 SD). Symphalangus (3.28 ± 0.68 SD) and Hylobates (3.22 ± 0.58 SD) exhibit the lowest average connectivity in the sample. For tortuosity (Table 1), Macaca exhibits the lowest average values (1.11 ± 0.19 SD). In Pan (1.18 ± 0.17 SD), humans (1.21 ± 0.21 SD), and Gorilla (1.23 ± 0.24 SD), the average tortuosity is larger than in Macaca but comparable to F I G U R E 4 Semi-automatic isolation of cancellous bone in the femoral head of Symphalangus syndactylus (a), the proximal humerus of Alouatta caraya (b), the distal fibula of Cercopithecus albogularis (c) and the brow ridge of Mandrillus sphynx (d). The 3D μCT scan is cut (red line) to limit the cancellous isolation to a region of interest. The results are here shown on a single 2D slice (indicated by the blue line on the 3D scan) and on the full 3D μCT stack (the cutting planes used to isolate the 3D regions of interest is shown in red) Papio (1.23 ± 0.25 SD). The highest tortuosity in the sample is displayed by Hylobates (1.29 ± 0.36 SD) and Symphalangus (1.26 ± 0.31 SD). Fractal dimension was calculated on the topological skeletons scaled on the height of the femoral head. For Pan and humans (

| DISCUSSION
This work presents novel strategies for the processing and analysis of cancellous bone. These approaches do not aim to substitute existing methods but rather to complement them. In addition, these strategies aim to overcome certain limitations of previous studies: the constraints associated to defining homologous cancellous subsamples and the characterization of the inherent complexity of trabecular structures.
The isolation protocol provides a flexible way of separating the cancellous from the compact bone thanks to the sequential application of image processing operators. Figure 4 demonstrates that the protocol can be fruitfully applied to skeletal elements of very different morphology to isolate cancellous bone with precision and avoiding manual segmentation. Also, this tool allows focusing the analysis of cancellous bone on the whole epiphysis rather than on subsamples, thus going past the issues raised in recent literature (Georgiou et al., 2019;Sylvester & Terhune, 2017;Tsegai et al., 2013;Tsegai et al., 2018b) and highlighted in the introduction of this work.
The indices presented in this work rely on the reduction of the cancellous shape to its topological skeleton, which is known to enhance certain geometrical and topological aspects of a shape, such as T A B L E 1 Complexity indices calculated on the topological skeleton of the cancellous bone in the femoral head Note: SD is shown only for the indices for which its calculation was possible. All indices are unitless, except for node density and the trabecular angle.
Fractal dimension is here presented as scaled on the height of the femoral head. For the definition and calculation of the indices, see main text.
F I G U R E 5 Node density of the femoral head, measured using a kernel density approximation over a regular 3D grid. It is expressed as the number of nodes of the skeletonized cancellous bone per cm 3 . The node density is here shown for a small sample of primates over the coronal (L-M-S-I) and parasagittal (A-P-S-I) planes. The density increases from blue to red. (A, anterior; P, posterior; S, superior; I, inferior; L, lateral; M, medial) connectivity, length and direction (Davies, 2004). Therefore, measuring indices of connectivity, tortuosity, density and, overall, complexity on the topological skeleton can be advantageous, because it reduces the effect of features not associated with the complexity of the trabecular lattice (e.g., trabecular thickness and shape). Furthermore, the indices are borrowed from, or inspired by, measurements used in previous literature and we place them within a reproducible framework.
With regard to the functional meaning of the indices, we want to stress that the example presented in this work is meant to illustrate the usage and interpretation of the indices. Clear functional implications cannot be drawn from the results and based on such a small comparative sample. Nevertheless, the results seem to suggest the potential of some of the indices in detecting functional differences: the human, Pan and Gorilla specimens, whose locomotion produces high mechanical load on the hind limbs, exhibit average node density and connectivity higher than the quadrupedal, and brachiating taxa; the lowest values of fractal dimension are exhibited by Hylobates and Symphalangus, whose arboreal lifestyle relies consistently on the forelimbs.
Summarizing, the strategy presented in this work innovates the following aspects of cancellous bone studies: (I) it introduces a tool for the isolation of cancellous bone avoiding manual segmentation; (II) it measures indices of complexity taking into account only the topology of the cancellous bone (without confounding, localized cancellous features); (III) provides descriptors of cancellous morphology that can be used to appreciate cancellous complexity numerically (Table 1) and visually (e.g., Figure 5). Further analysis is needed to clarify the extent of functional signal detected by the indices.

ACKNOWLEDGMENTS
The authors wish to thank the curators at the various institutions where data were collected: the Natural History Museum "La Specola,"

DATA AVAILABILITY STATEMENT
Part of the μCT data used in this study is available through the online database Morphosource (https://www.morphosource.org/). Some data is property of the institutions where the specimens were scanned and cannot be made publicly available. Access to this material can be granted by those institutions upon reasonable request. Information about the scans and their identifiers in Morphosource or at their institution are provided in the supplementary information (Table S1 and S2).