Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Med Biol. Author manuscript; available in PMC 2017 September 21.
Published in final edited form as:
PMCID: PMC5158528

Trabecular Bone Characterization on the Continuum of Plates and Rods using In Vivo MR Imaging and Volumetric Topological Analysis


Osteoporosis is associated with increased risk of fractures, which is clinically defined by low bone mineral density. Increasing evidence suggests that trabecular bone (TB) micro-architecture is an important determinant of bone strength and fracture risk. We present an improved volumetric topological analysis algorithm based on fuzzy skeletonization, results of its application on in vivo MR imaging, and compare its performance with digital topological analysis. The new VTA method eliminates data loss in the binarization step and yields accurate and robust measures of local plate-width for individual trabeculae, which allows classification of TB structures on the continuum between perfect plates and rods. The repeat-scan reproducibility of the method was evaluated on in vivo MRI of distal femur and distal radius, and high intra-class correlation coefficients between 0.93 and 0.97 were observed. The method’s ability to detect treatment effects on TB micro-architecture was examined in a two-year testosterone study on hypogonadal men. It was observed from experimental results that average plate-width and plate-to-rod ratio significantly improved after 6 months and the improvement was found to continue at 12 and 24 months. The bone density of plate-like trabeculae was found to increase by 6.5% (p = 0.06), 7.2% (p = 0.07) and 16.2% (p = 0.003) at 6, 12, 24 months, respectively. While the density of rod-like trabeculae did not change significantly, even at 24 months. A comparative study showed that VTA has enhanced ability to detect treatment effects in TB micro-architectural as compared to conventional method of digital topological analysis for plate/rod characterization in terms of both percent change and effect-size.

Keywords: Trabecular bone, volumetric topological analysis, fuzzy skeletonization, plate-width, magnetic resonance imaging

1 Introducation

Osteoporosis is a degenerative bone disease of high prevalence affecting older adults. The disorder is associated with high risk of low-trauma fractures leading to increased mortality and morbidity, and it entails significant healthcare costs. An estimated 8 million women and 2 million men in the United States meet the diagnostic criteria for osteoporosis [1]. The percentage of patients with osteoporosis increases progressively with age – among US women, 13% in their 50s, 27% in their 60s, 47% in their 70s, and 67% in their 80s have bone density levels in the diagnostic range of osteoporosis [2]. Clinically, osteoporosis is defined by low bone mineral density (BMD), which explains 65–75% of the variance in bone strength [3, 4], while the remaining variance is due to the cumulative and synergistic effects of various factors, including trabecular bone (TB) micro-architecture. There is evidence in literature [5, 6] demonstrating that effects of therapeutic agents are greater on TB as compared to overall BMD. Results from several studies [69] support the hypothesis that quality of TB micro-architecture is an important determinant of bone strength and fracture risk. Effective measures of TB micro-architecture from in vivo imaging are useful to assess bone strength and fracture risk in the realm of clinical therapy and treatment guidance, including growth and peak bone accrual, aging, post-menopausal bone loss, cancer-related bone loss and conditional bone loss, such as eating disorders, renal osteodystrophy, osteogenesis imperfacta, osteoarthritis, rheumatoid arthritis, corticosteroid medicinal intake, etc.

The conventional method for quantitative assessment of TB micro-architectural quality is a two-dimensional (2-D) cross-sectional histomorphometry from bone biopsies [10]. Over the past two decades, remarkable progress has been achieved on high-resolution in vivo bone imaging techniques [11], including magnetic resonance (MR) imaging (MRI) [4, 12, 13], high resolution peripheral quantitative computed tomography (HR-pQCT) [14, 15] and multi-row detector CT (MD-CT) [16, 17] imaging, to assess three-dimensional (3-D) TB micro-architecture. As compared to other bone imaging modalities, a major advantage of MRI is that it involves no ionizing radiation, which is beneficial for longitudinal studies.

Various topologic and geometric analytic approaches have been reported [1823] for characterizing TB microarchitecture. Parfitt et al. [18] conceived a parallel interconnected plate model of TB yielding bone area fraction, TB volume fraction, spacing, and number from 2-D histomorphometric sections. Vesterby et al. [19] conceived a new stereologic parameter, called star volume, which is the average volume of an object region that can be seen from a point inside that region unobscured in all directions. Hahn et al. [20] introduced the “trabecular bone pattern factor” which captures TB connectivity in terms of the convexity property of the TB surface defined as the ratio of the differences in perimeter and area under dilation. Hildebrand et al. [21] developed a 3-D structure model index, a function of global plate-to-rod ratio, based on the observation that the rate of volume change with thickness for a plate is different from that for a rod. Majumdar et al. [22] adopted apparent TB number, thickness, spacing and fractal dimension measures to quantify TB structural quality. Feldkamp et al. [23] showed that the makeup of TB networks can be expressed in terms of topological entities such as the 3-D Euler number.

There are several indirect MRI methods, which do not directly apply structural analysis to high-resolution intensity images. Instead, these methods exploit different contrast mechanisms to derive metrics that relate to TB microstructure. High-resolution MRI of internal gradient diffusion weighting provides the bulk measurement of the surface-to-volume ratio based on the principle that, in the highly porous TB system, the largest internal gradients exist close to TB surfaces perpendicular to the applied field [24]. MR relaxometry is another indirect method that characterizes TB orientation using magnetic field inhomogeneity at the interface between TB and marrow [25]. The method yields information on trabecular density and orientation but cannot provide detailed insight into structural make-up of the bone [26]. Bone marrow perfusion MRI was identified as a biomarker for bone quality. Dynamic contrast-enhanced MRI measures bone marrow perfusion indices representing enhancement peak and maximum slope [27, 28]. These MRI-based indirect methods provide bulk measures bone microstructure related to TB and marrow interface. In addition, morphologic interpretation of these indirect measures is not straightforward, thereby adding difficulties in relating these measures to detailed microarchitecture in clinical studies.

Saha and his colleagues developed digital topological analysis (DTA) [2932] as a means to characterize topological plates, rods, and junctions at individual TB voxels, which has been successfully applied in several studies [13, 3335]. More recently, volumetric topological analysis (VTA) was proposed to characterize individual TB voxels on the continuum between a perfect plate and a perfect rod. The method has been validated on clinical MD-CT imaging at peripheral sites [17]. Here, we improved the VTA algorithm using a recently developed fuzzy skeletonization algorithm [36] to eliminate data loss in the binarization step and evaluated it on in vivo MRI. The improved method yields accurate and robust measure of plate-widths. We propose a new characterization of plate-like structures using plate-width information. The repeat scan reproducibility of the method was examined using two different in vivo MRI protocols at distal femur and distal radius. In addition, we examined the method’s ability to detect testosterone treatment effects on hypogonnadal men in a two-year follow-up study and the results are presented. The performance of the method is compared with the conventional methods of digital topological analysis for plate/rod characterization.

2 Methods

2.1 Volumetric Topological Analysis

The overall objective of VTA was to measure the plate-width of individual trabeculae in units of micrometers and locally classify individual trabecular type on the continuum between a perfect plate and a perfect rod (Figure 1(a)). Essentially, the VTA algorithm computes local trabecular plate-width (Figure 1(b)), which is different from either trabecular thickness or skeletal width. The original VTA algorithm [17] was developed for binary objects requiring thresholding on fuzzy representations of TB images, which is a sensitive and undesired step at in vivo resolution [11]. Here, we present an improved algorithm, which is directly applicable to fuzzy TB images eliminating the thresholding step—a major source of error, especially, at in vivo image resolution.

Figure 1
(a) Trabecular bone plate/rod classification of trabecular bone using volumetric topological analysis. (b) Various measures on a schematic drawing of a plate-like structure. (c) Color-coded geodesic distance transform from the surface edge, providing ...

The VTA algorithm consists of five sequential steps: (1) fuzzy skeletonization, (2) digital topological analysis, (3) geodesic distance transform, (4) geodesic scale computation, and (5) volumetric feature propagation. The principle of VTA is described in Figure 1(c,d). Let us assume that (c) represents the surface skeleton of a plate-like object with varying plate-widths. A surface representation of an object is obtained by surface skeletonization [31, 36]. After surface skeletonization, the VTA algorithm identifies different topological entities, including surface edges (the thick outline in Figure 1(c,d)), on the skeleton using digital topological analysis. The geodesic distance between two points on a surface skeleton is the length of the shortest path between the two points such that the entire path lies on the skeletal surface. VTA computes the geodesic distance at every surface interior point from its nearest surface edge, which is referred as the geodesic distance transform (GDT) (Figure 1(c)).

The central dotted line in both (c) and (d) represents the arc-skeleton, and the points on the arc-skeleton are referred to as axial points. Note that, essentially, the computed GDT is equal to the half of the local plate-width at an axial point, e.g., the point a in (d). However, GDT fails to determine the local plate-width at non-skeletal points, e.g, the point p in (c). A feature propagation step is applied, where a non-axial point inherits the local plate-width from its most-representative axial point. For example, the point p is assigned the same plate-width as that of a after feature propagation as illustrated in (d). Finally, another level of feature propagation is applied where the plate-width values are propagated from the surface-skeleton to the entire object volume, e.g., from Figure 2(d) to (e). In this paper, we apply fuzzy skeletonization in Step 1 instead of binary skeletonization. Steps 2 to 5 are implemented using methods described in our previous work [17]. The feature propagation in Step 5 is accomplished using the principle established by Liu et al. [37] that is independent of scan or processing order. Results of intermediate steps of VTA are presented in Figure 2.

Figure 2
Intermediate steps of VTA. (a) A TB region from a μCT image of a cadaveric distal tibia. (b) Classified topological entities including plates (green), rods (red), edges (light colors), and junctions (blue) on the fuzzy skeleton of (a). (c–e) ...

The following notations are used in this paper. Let O = {(p, BVF(p)) | p [set membership] Z3} be the bone volume fraction (BVF) representation of a TB image, where Z is the set of integers, Z3 is the rectangular grid, p is a voxel location, and BVF(p) is the BVF value at p. Let O denote the set of voxels with non-zero BVF, i.e., O = {p | p [set membership] Z3BVF(p) ≠ 0}. In the following, we describe the basic principle of the fuzzy skeletonization algorithm, which eliminates data loss in the binarization step and improves the performance.

2.2 Fuzzy Skeletonization

The performance of our TB plate-width computation method is highly dependent on the quality of the surface skeleton generated from the volumetric representation of a TB network acquired by in vivo imaging. In our previous algorithm, a binary skeletonization algorithm was used, which is always associated with binarization related data loss adding skeletal inaccuracies such as disruption of trabecular rods, perforation of plates, and filling small marrow holes. Recently, Jin and Saha [36] developed a fuzzy skeletonization algorithm based on fuzzy grassfire propagation, which eliminates binarization related data loss. Blum’s grassfire transform, originally defined for binary objects, was modified for fuzzy representation of TB images where the BVF is used to define the instantaneous speed of the grassfire front at a given bone voxel. The speed of grassfire front at a given TB voxel p is inversely proportional to its BVF. Using this notion, it can be shown that the fuzzy distance transform (FDT) [38] value at p defines the time that the fire front reaches p. Therefore, the propagation time of a fire front at a TB voxel p to its neighbor q is equal to the local BVF-weighted distance between p and q, and this equality is violated only at skeletal or quench points where the propagation process is interrupted and stopped. Thus, a TB voxel p is a fuzzy quench voxel (Figure 3(c)) if the following inequality holds for every neighbor q of p

Figure 3
Results of intermediate steps of fuzzy skeletonization. (a) 3-D display of TB region in a μ-CT image of a cadaveric distal tibia specimen. (b) A sagittal image slice displaying the fuzziness in the image. (c) All quench voxels before filtering. ...

Following the above definition, a quench voxel fails to pass the grassfire front to any of its neighbors, and thus, grassfire is extinguished at quench voxels.

Although the fuzzy quench voxel captures medial axes or symmetry-structures using the notion of fuzzy grassfire transform, it creates a large number of spurious quench voxels. Spurious quench voxels generate noisy skeletal branches, which are not necessary to describe the overall object shape and topology. Often, such noisy branches add undesired complexities in the skeleton and reduce the effectiveness of skeleton-derived measures. Therefore, it is necessary to remove spurious quench voxels, which is accomplished using their significance measures. Here, a function that resembles the “collision impact” of individual TB voxels is used to determine the significance of a quench voxel. Collision impact ξ of a quench voxel p is defined as follow:


where the function f+(x) returns the value of x, when x ≥ 0, and ‘0’, otherwise, and N*(p) is the set of 26-neighbors of p. There are two types of quench voxels – a surface quench voxel where two opposite fire fronts meet and a curve quench voxel where fire fronts meet from all directions on a plane. These unique geometric properties of a surface or a curve quench voxel, are optimally used to determine its significance based on the collision impact values of its neighboring voxels. For example, a significant surface quench voxel ensures a 3 × 3 structure of quench voxels with high collision impact in its neighborhood, while a significant curve quench voxel requires a three-voxel curve in its neighborhood with high collision impact. A quench voxel with its significance value above a threshold value thrsignificance, is referred as a significant quench voxel or a fuzzy axial voxel. As illustrated in Figure 3(d), the set of fuzzy axial voxels captures the essential geometry of the original TB object; also, it reduces a large set of non-significant quench voxels of (c).

In this work, the fuzzy skeletonization step starts with identifying the set of all fuzzy quench voxels followed by computation of fuzzy axial voxels using the above filtering method for significant quench voxels. Subsequently, TB object voxels, which are not fuzzy axial voxels, are sequentially removed in the increasing order of their FDT values while preserving the topology of the object. The notion of 3-D simple points [29] is used for topology preservation. It should be noted that, despite the filtering of quench voxels, a few spurious branches may survive in the skeleton. Therefore, a post-skeletal-pruning step is applied to prune noisy branches based on their global significance computed as the collision-impact-weighted branch length. Figure 3 shows results of intermediate steps of fuzzy skeletonization.

Fuzzy skeletonization reduces binarization related data loss, which improves the preservation of trabecular network connectivity especially at relatively thin trabeculae. Figure 4(a) shows the BVF map of TB network in an axial image slice of distal radius acquired in vivo at 1.5 T MRI. Results of binary and fuzzy skeletonization are presented in (b,c), respectively. A few examples of trabecular connectivity loss in the binary skeleton shown in (b) are highlighted in green; the matching regions in the fuzzy skeleton are also highlighted. It should be noted that the examples of in-plane connectivity loss of TB network as shown in (b) do not contradict with 3D topology preservation using simple points [29], which guaranteed by most skeletonization algorithms [39].

Figure 4
Connectivity preservation of TB network using binary and fuzzy skeletonization methods. (a) The BVF map in an axial image of distal radius acquired in vivo at 1.5 T MRI. The regions with relatively thin trabeculae are highlighted with arrows. (b,c) Results ...

2.3 Trabecular Bone Measures

The output of the VTA algorithm is a local plate-width (PW) function PW: OR+, where PW(p) yields the local TB plate-width in units of micrometers at a given bone voxel p [set membership] O, R+ are real numbers. Three different measures of average bone density, namely, bone volume by total volume (BV/TV), plate-like BV/TV, and rod-like BV/TV, over a target region of interest (ROI) V are defined as follows:




The threshold THplate was determined based on observed results of TB BMD distribution over the range of TB plate-width (see Figure 7). A clear deflation of BMD distribution at the plate-width of ~450μm was noted on all curves in Figure 7. This deflation point of 450μm was used to determine the threshold separating plate-like and rod-like trabeculae.

Figure 7
The baseline trabecular BVF distribution at various TB plate-widths among hypogonadal men from the testosterone treatment study. The shaded area represents the mean ± std of trabecular BVF values at corresponding TB plate-width.

The platelikeness and rodlikeness measures of a TB voxel p are defined using the eccentricity between the average TB thickness and the local plate-width as follows:



The parameter avth was determined as the average TB thickness reported in [40]. Note that the platelikeness measure is modelled using the eccentricity of an ellipse; here the smaller diameter captures the structure thickness, while the larger one represents TB plate-width. It should be emphasized that the above formulation of platelikeness and rodlikeness requires no threshold value. The average plate-width and plate-to-rod ratio measures are defined as follows:



Besides the above measure, two measures of TB plate/rod micro-architecture, namely SCRDTA (surface-to-curve ratio) and EIDTA (erosion index), using the conventional method of digital topological analysis were computed using the protocol prescribed in [32].

3 Experiments

The aims of our experimental plans were—(1) evaluation of reproducibility for in vivo repeat MRI scans using VTA measures at wrist and knee and (2) assessment of the method’s ability to detect treatment effects in hypogonadal men subjected to testosterone treatment. For the first experiment, ultra high-field 7 Tesla (T) MRI of the distal femur at the right knee as well as 1.5 T MRI of the distal radius at the right wrist were used. For the second experiment, 1.5 T MRI of distal tibia was used.

MR Imaging

All MR datasets used in this paper were collected in vivo and several properties have previously been published [5, 7, 33, 41]. The knee MRI was collected at New York University [41]. The right distal femur of each subject was scanned on a 7T whole body MR scanner (Siemens Medical Solutions, Erlangen, Germany) using a quadrature knee coil (18 cm diameter, transmit-receive). A high-resolution 3-D-fast low angle shot (FLASH) sequence was employed to obtain all images (TR/TE = 20/4.5; flip angle 10 degrees; bandwidth 130 Hz/pixel; one signal acquired; 130 axial images with resolution 0.195 mm × 0.195 mm × 1 mm). Scanning time was ~12 min total.

The wrist MRI data for assessment of reproducibility had been acquired previously at the University of Pennsylvania [33]. The right distal radius of each subject was acquired on a Siemens 1.5 T MAGNETOM Sonata MR scanner (Siemens Medical Solution, Erlangen, Germany), with a protocol described in detail in [33]. In brief, 3-D fast large-angle spin echo (3-DFLASE) sequence (36) yielded 0.137 mm × 0.137 mm × 0.410 mm voxel size. The imaging ROI of a 70 mm × 42 mm × 13 mm slab centered 7 mm proximal to the most proximal tip of the epiphyseal line was located during the first visit using a number of low resolution orthogonal 2-D localizer sequences. For follow-up exams, the 13 mm slab was selected by means of prospective registration [42].

For the second experiment assessing the sensitivity of the new bone measures to detect structural abnormalities, prior MRI data of the right distal tibia were re-examined to evaluate the TB micro-structural implications of hormone therapy in hypogonadal men. Full details of the prior studies are in [5, 7].


For the knee reproducibility experiment, MR bone data of four healthy male subjects (35±7.8 years) [41], without history of fracture, musculoskeletal disorder, or any history of intake of bone altering medication were recruited in the Department of Radiology at New York University,. Each subject was scanned twice within the same day with repositioning between scans. For the wrist reproducibility experiment, previously acquired MR data of twenty women (62±8.1 years; 17 postmenopausal, three premenopausal), were used [33]. Exclusion criteria included a history of fracture or treatment for osteoporosis, a body mass index (BMI) greater than 30 kg/m2, and the presence of primary bone cancer or metastases to the bone. Each subject had been scanned three times over the course of 8 weeks with average interval between scans being 20.2 days (standard deviation = 14.5 days).

For the second experiment, TB MR data of a previous longitudinal study [7] on hypogonadal men under testosterone treatment were used. Specifically, 10 previously untreated hypogonadal men (age range: 18 to 80 years) were recruited under that study, and their TB MRI data at baseline, 6, 12, and 24 month visit under testosterone treatment [7] were collected. Full details of human subjects and imaging protocols are given in [7]. In brief: all 10 men had secondary hypogonadism for a duration of more than 2 years with average duration of hypogonadism of 7.9 years (S.D. = 8.2 years). Men were excluded if they had been consuming less than 750 mg of calcium per day as determined by a food frequency questionnaire.

Image processing

Image processing steps are as follows: (1) BVF computation, (2) interpolation to 150 μm isotropic voxel size, (3) analysis using the improved VTA method, (4) ROI selection, and (5) computation of TB measures. BVF images were computed using intensity connectivity and a local marrow intensity computation approach without requiring global thresholding. In a BVF image, pixel intensity corresponds to the fractional occupancy of bone. BVF images were resampled using the windowed-sinc interpolation method producing 150 μm isotropic voxel size, and the images were subjected to VTA analysis. We first draw the ROI manually on the first MR scan and generated ROIs of follow-up scans using transformation matrices that register original image of the first scan to follow-up scans. The process involves only ROI transformation and avoids interpolation related resolution loss of original BVF images. Finally, we compute proposed measures for each BVF image.

Statistical analyses

Statistical analyses that quantitatively assess reproducibility of the VTA derived TB measures in vivo repeat MR scans are: the mean, the root-mean-square of coefficients of variation, and the intraclass correlation coefficient. For each subject, the coefficient of variation (CV) of a given TB measure was calculated across repeat scans. For a given study and a given measure, an average CV was obtained as the root mean square average of CV across all subjects (denoted RMS-CV) to measure the global variability across repeat scans. In addition, the intra-class correlation coefficient (ICC) was computed as a measure of reliability via one-way analysis of variance (ANOVA). To examine a method’s ability to detect longitudinal changes in TB micro-architecture in the second experiment, we computed both p-values and effect-size of changes. The effect-size was computed as the mean paired change of a measure between baseline and follow-up normalized by the standard deviation of the paired changes among individual subjects.

4 Results

4.1 Reproducibility of TB Plate/Rod Measures in MR Repeat Scans

Matching image slices of distal radius from post-registered wrist MR repeat scans along with parametric images are shown in Figure 5. Color-coded renditions of VTA plate-rod classification over matching ROIs in three repeat scans of a subject from the wrist study is presented in Figure 6. Quantitative results of the reproducibility study are summarized in Table 1. RMS-CV values of repeat scans for knee and wrist MRI are in the range of [1.0%, 4.7%] while the intra-class correlation coefficient (ICC) values are in the range of [0.93, 0.97]. Similar values of ICC for reproducibility of DTA parameters SCRDTA and EIDTA were reported in [33]. For the knee MRI data, we computed the ICC values for SCRDTA and EIDTA and the observed values are 0.91 and 0.93, respectively.

Figure 5
Reproducibility of VTA-based characterization of TB plate/rod micro-architecture for in vivo repeat MR scans of the distal radius on a 1.5T whole-body MR scanner: (a) axial image and ROI (red). (b) Bone volume fraction (BVF). (c) VTA-based color-coded ...
Figure 6
Reproducibility of the TB plate-width measure for in vivo MR repeat scans of distal radius. (a) Color-coded 3-D rendition of local TB plate-width in post-registered virtual image cores from three MR repeat scans of a subject. (b,c) Same as (a) for two ...
Table 1
Results of in vivo repeat scan reproducibility analysis for the knee and the wrist MRI studies.

4.2 Effects of Testosterone Treatments on T B Micro-Architecture in Hypogonadal Men

The baseline distribution of trabecular BVF at various TB plate-width among the hypogonadal men used in the testosterone study is presented in Figure 7. Changes in TB plate-rod micro-architecture at different follow-ups under the hormone therapy are presented in Figure 8. Figure 8(a–d) presents color-coded three-dimensional TB plate/rod micro-architecture over matching ROIs from post-registered distal tibia MRI of a hypogonadal man at baseline, 6-, 12-, and 24-month follow-ups. Figure 8(e) depicts the mean changes in trabecular BVF distributions across TB plate-widths at three follow-up time points; see Table 2 for quantitative results.

Figure 8
Changes in TB plate/rod micro-architecture over a 24-month MRI study monitoring effects of hormone therapy in hopogonadal men. (a–d) Color-coded display of TB plate/rod classification from distal tibia MRI at baseline, 6-, 12-, and 24-month follow-ups. ...
Table 2
Progressive changes in VTA-based TB measures derived from a two-year follow-up distal tibia MRI study during testosterone hormone therapy in hypogonadal men.

Figure 9 compares the effect-size of changes in different DTA and VTA measures of TB plate/rod micro-architecture under hormone therapy at baseline, 6-, 12-, and 24-month follow-ups. High effect-size was observed for the two VTA measures PRRVTA and PWVTA than the two DTA measures surface-to-curve ratio (SCRDTA), erosion index (EIDTA), especially, at 6- and 24-month follow-ups.

Figure 9
Effect-size of changes in different TB measures from MRI at 6-, 12-, and 24-month follow-ups under the testosterone hormone therapy in ten hopogonadal men.

5 Discussion

5.1 TB Plate/Rod Measures and Their MR Repeat Scan Reproducibility

The significance of plate/rod distribution in trabecular bone has long been recognized. A large number of histologic studies [8, 9, 43] have confirmed the relationship between erosion of trabeculae from plates to rods and higher fracture risk. Kleerekoper et al. [8] found lower mean TB plate density in individuals with osteoporotic vertebral compression deformities compared with a BMD-matched control group without fractures. Recker [43] found differences in TB plate density on two BMD-matched groups with (lower density) and without (higher density) vertebral crush fractures. Various approaches have been reported to distinguish rod-like from plate-like trabeculae. Hahn et al. [20] expressed the relation of trabecular plates to rods in terms of the ratio of concave to convex surfaces of the bone pattern in two-dimensional histologic bone sections. Specifically, the ratio of trabecular plates and rods was defined in terms of the ratio of the differences in perimeter and area under a simulated dilation. This line of thought was extended into a direct 3-D measure of structure model index [21], where the global plate-to-rod ratio was expressed in terms of the ratio of volume change with thickness under simulated 3-D dilation. Notably, however, these prior studies were all based on bone biopsies.

In previous work, Saha and his colleagues developed DTA [2932] that characterizes topological plates, rods, and junctions at individual TB voxels. This method was applied to in vivo studies assessing implications of hormonal changes [7, 44] or renal osteodystrophy [45] on TB plate-rod micro-architecture. DTA based assessment of TB micro-architecture revealed reduced trabecular plates in subjects with vertebral fracture [46] or deformity [34]. However, a major limitation of the DTA method is that resulting classifications are inherently discrete, failing to distinguish between narrow and wide plates. The balance between plates and rods during bone formation at younger ages, as well as during bone loss or anti-resorptive treatment, changes gradually and, therefore, demands classification of TB micro-architecture on the continuum between a perfect plate and a perfect rod. Also, theoretically, the discrete TB plate/rod classification by DTA is sensitive to image voxel size. For example, a trabecula classified as rod at one voxel size may be classified as a plate at a smaller voxel size. In contrast, the VTA algorithm solved the problem by computing local plate-width in the units of micrometers for individual trabeculae. Color-coded illustrations of the TB network in Figure 5 and Figure 6 confirm the method’s ability to classify individual trabeculae on the continuum between a perfect plate and a perfect rod. Regional agreements of TB plate/rod classification in repeat MR scans are visually observed in these figures. The three rows in Figure 6 represent repeat scans of TB from three different subjects with distinct trabecular plate/rod composition. The method successfully represents the large difference across the three subjects and shows much smaller variance at repeat scans. In Sec 5.2, we also compared DTA with VTA in measuring quality changes of TB micro-architecture.

It is noted in Table 1 that RMS-CV values are small for all TB measures for both the knee and the wrist MRI studies. The lowest RMS-CV value of 1.0 was observed for the core VTA measure PWVTA, which indicates that the PWVTA measure is reproducible with 1.0% variability relative to its mean value. The ICC, i.e., the variability of the measure relative to the range of values obtained across the population, was found to be 0.93 and 0.95 for the same measure for the two studies. PRRVTA shows somewhat greater RMS-CV since the variability in both platelikeness and rodlikeness measures are combined in this ratio measure (see Eqn. 9). However, the observed ICC values of PRRVTA are similar to those of PWVTA, which suggests that the across-population variation of PRRVTA is higher than that of PWVTA. High reproducibility of the VTA measures establishes their suitability for clinical or research studies including multi-center studies and longitudinal ones where multiple follow-up scans are needed.

Beyond the measures given in Table 1, the method yields unique information on the distribution of BV/TV at different values of trabecular plate-width (Figure 7). These distributions provide insight into the pathogenesis of bone degeneration (or, development) process at the level of TB micro-architecture. For example, besides early detection of trabecular bone conversion from plates to rods, the method is potentially capable of distinguishing gradual conversion of trabecular plates to rods from rapid perforation of trabecular plates. We observed that there is an apparent separation gap at around 450μm, which differentiates rod- and plate-like trabeculae.

5.2 Effects of Testosterone Treatments on TB Micro-Architecture in Hypogonadal Men

Benito et al. [5] reported the results of a two-year follow-up study [5] examining the effects of testosterone replacement treatment on TB micro-architecture using the same imaging tools. The MR image from these two studies were re-analyzed using the improved VTA algorithm. Here, we present the result and compare it with DTA in assessing effects of testosterone treatment on TB micro-architecture among hypogonadal men.

Figure 7 presents the distribution of TB volume at different plate-widths for hypogonadal men at different time-points in the follow-up study. Here, the VTA algorithm provides unique information not obtainable by DTA or any other existing tools. Figure 8(a–d) depicts the change in TB micro-architecture over a matching ROI for a hypogonadal subject during the two-year treatment period. A gradual and progressive shift in trabecular composition toward more plate-like structures is noted in the figure. The quantitative results of Table 2 confirm this effect on TB micro-architecture of the hormone treatment. The mean plate-width PWVTA and plate-to-rod ratio PRRVTA showed statistically significant changes as early as 6 months of the treatment. After 24 months of treatment, all TB measures showed statistically significant improvements, except BV/TVrod over rod-like trabeculae. BV/TVplate was found to increase 6.5% at the 6 months but did not quite reach significance (p = 0.07). Distributions of changes in BVF at different trabecular plate-width during the two-year treatment period are presented in Figure 8(e). At 12- and 24-month follow-ups, there is apparent increase in BVF around ~625 μm plate-width, which is close to the mean TB plate-width observed in Table 1. In Table 2, quantitative results of changes for two-year follow-up study on hypogonadal men during testosterone hormone therapy are presented. Figure 9 presents the effect-size of changes of VTA and DTA measures at different follow-ups.

The above results shed light on the mechanism of hormone therapy in terms of TB micro-architecture. A net increase in BV/TVplate occurs due to gradual conversion of TB rods into plates under hormone therapy, and the shift is detectable as early as at six months. On the other hand, the relative constancy of the BV/TVrod measure is associated with the uncertainty of the balance between the rates of rod to plate conversion and generation of new TB rods. While the conversion of TB rods to plates decrease rod density, generation of new TB rods increases it. These two opposite effects nearly cancel each other.

There is histologic evidence demonstrating improvements in TB microarchitecture in hypogonadal men after testosterone treatment. Baran et al. [47] reported increases in relative osteoid volume, total osteoid surface, linear extent of bone formation, and bone mineralization in an osteoporotic hypogonadal male after 6 months of testosterone replacement. Francis et al. [48] observed increased total and free plasma 1,25-dihydroxyvitamin D in hypogonadal men on testosterone treatment, and bone biopsy at 6 weeks after treatment confirmed increase in skeletal retention of calcium and bone formation by testosterone treatment. It is in accordance with our observations of positive changes in TB micro-architectural changes after testosterone treatment. The early detection of positive changes in TB microarchitecture using VTA in hypogonadal men after 6 months of testosterone treatment (Table 2) suggests that VTA may be a more responsive or sensitive way to monitor osteoporosis or treatment effect as compared to DTA or BV/TV, which required more time to detect treatment effects. Specifically, as compared to DTA measures, two VTA measures plate-width PWVTA and plate-to-rod ratio PRRVTA show statistically significant changes under treatment at six months, while none of the DTA measures produced significant difference at six months. In addition, VTA offers resolution-independent measures in units of micrometers, and insight into the mechanism of a treatment at the level of individual trabecular plate/rod micro-architecture.

6 Conclusion

We improved the VTA algorithm to eliminate data loss in the binarization step using fuzzy skeletonization. The method offers reproducible measures of TB plate/rod micro-architecture at in vivo MRI. It computes local plate-width and characterizes individual trabeculae on the continuum of a perfect plate and perfect rod. The method yields the width information of individual trabeculae and offers an insight speculation into the bone structure at micrometer level, which existing tools are not able to provide. By studying a group of human subjects, we observe a pattern of width distribution of human TB. The method is suitable for cross-sectional and follow-up studies toward answering clinical and biological questions. The method should be equally suited for studies with other modalities providing trabecular bone images with limited spatial resolution, including high-resolution CT. However, evaluation of the method’s full potential requires drug intervention studies in larger cohorts of patients.

The method presented in this paper is based on direct microstructural analysis of TB network in high-resolution MR imaging, and provides measures of plate/rod distribution of individual trabeculae. Unlike the indirect methods, the morphologic interpretations of our measures are well defined, and histologic evidence confirms the relationship between osteoporosis and the gradual conversion of trabecular plates to rods, a process well known to increase fracture risk. Further, our method offers plate/rod characterization at the level of individual trabeculae, which may be useful as a regional biomarker, especially, in recognizing heterogeneous bone loss. Lastly, it would be worthwhile to compare the relationship between high-resolution image-based direct morphologic measures of TB micro-architecture and indirect MRI measures and their strengths and weaknesses in assessing fracture risk and their ability to monitor disease progression or response to treatment.


This work was supported by the NIH grants R01AR054439, R01AR056260, K23AR059748, R01AR41443 and R01AR053156.


1. Melton LJ. How Many Women Have Osteoporosis Now. J Bone Min Res. 1995;10:175–177. [PubMed]
2. Looker AC, Orwoll ES, Johnston CC, Lindsay RL, Wahner HW, Dunn WL, Calvo MS, Harris TB, Heyse SP. Prevalence of low femoral bone density in older US adults from NHANES III. J Bone Min Res. 1997;12:1761–1768. [PubMed]
3. Ammann P, Rizzoli R. Bone strength and its determinants. Osteoporos Int. 2003;14(Suppl 3):S13–18. [PubMed]
4. Wehrli FW, Saha PK, Gomberg BR, Song HK, Snyder PJ, Benito M, Wright A, Weening R. Role of magnetic resonance for assessing structure and function of trabecular bone. Top Magn Reson Imag. 2002;13:335–356. [PubMed]
5. Benito M, Vasilic B, Wehrli FW, Bunker B, Wald M, Gomberg B, Wright AC, Zemel B, Cucchiara A, Snyder PJ. Effect of testosterone replacement on trabecular architecture in hypogonadal men. J Bone Miner Res. 2005;20:1785–91. [PubMed]
6. Chesnut CH, III, Majumdar S, Newitt DC, Shields A, Van Pelt J, Laschansky E, Azria M, Kriegman A, Olson M, Eriksen EF, Mindeholm L. Effects of salmon calcitonin on trabecular microarchitecture as determined by magnetic resonance imaging: results from the QUEST study. J Bone Miner Res. 2005;20:1548–1561. [PMC free article] [PubMed]
7. Benito M, Gomberg B, Wehrli FW, Weening RH, Zemel B, Wright AC, Song HK, Cucchiara A, Snyder PJ. Deterioration of trabecular architecture in hypogonadal men. J Clin Endocrinol Metab. 2003;88:1497–502. [PubMed]
8. Kleerekoper M, Villanueva AR, Stanciu J, Rao DS, Parfitt AM. The role of three-dimensional trabecular microstructure in the pathogenesis of vertebral compression fractures. Calc Tiss Int. 1985;37:594–597. [PubMed]
9. Legrand E, Chappard D, Pascaretti C, Duquenne M, Krebs S, Rohmer V, Basle MF, Audran M. Trabecular bone microarchitecture, bone mineral density and vertebral fractures in male osteoporosis. J Bone Min Res. 2000;15:13–19. [PubMed]
10. Parfitt AM, Drezner MK, Glorieux FH, Kanis JA, Malluche H, Meunier PJ, Ott SM, Recker RR. Bone histomorphometry: standardization of nomenclature, symbols, and units. Report of the ASBMR Histomorphometry Nomenclature Committee. J Bone Miner Res. 1987;2:595–610. [PubMed]
11. Krug R, Burghardt AJ, Majumdar S, Link TM. High-resolution imaging techniques for the assessment of osteoporosis. Radiol Clin North Am. 2010;48:601–21. [PMC free article] [PubMed]
12. Majumdar S. Magnetic resonance imaging of trabecular bone structure. Top Magn Reson Imaging. 2002;13:323–34. [PubMed]
13. Chang G, Pakin SK, Schweitzer ME, Saha PK, Regatte RR. Adaptations in trabecular bone microarchitecture in Olympic athletes determined by 7T MRI. J Magn Reson Imaging. 2008;27:1089–95. [PMC free article] [PubMed]
14. Boutroy S, Bouxsein ML, Munoz F, Delmas PD. In Vivo Assessment of Trabecular Bone Microarchitecture by High-Resolution Peripheral Quantitative Computed Tomography. J Clin Endocrinol Metab. 2005;90:6508–6515. [PubMed]
15. Burghardt AJ, Pialat JB, Kazakia GJ, Boutroy S, Engelke K, Patsch JM, Valentinitsch A, Liu D, Szabo E, Bogado CE, Zanchetta MB, McKay HA, Shane E, Boyd SK, Bouxsein ML, Chapurlat R, Khosla S, Majumdar S. Multicenter precision of cortical and trabecular bone quality measures assessed by high-resolution peripheral quantitative computed tomography. J Bone Miner Res. 2013;28:524–36. [PMC free article] [PubMed]
16. Link TM, Vieth V, Stehling C, Lotter A, Beer A, Newitt D, Majumdar S. High-resolution MRI vs multislice spiral CT: Which technique depicts the trabecular bone structure best? Eur Radiol. 2003;13:663–671. [PubMed]
17. Saha PK, Xu Y, Duan H, Heiner A, Liang G. Volumetric topological analysis: a novel approach for trabecular bone classification on the continuum between plates and rods. IEEE Trans Med Imag. 2010;29:1821–1838. [PMC free article] [PubMed]
18. Parfitt AM, Mathews CHE, Villanueva AR, Kleerekoper M, Frame B, Rao DS. Relationships between surface, volume, and thickness of iliac trabecular bone in aging and in osteoporosis. Implications for the microanatomic and cellular mechanisms of bone loss. J Clin Invest. 1983;72:1396–1409. [PMC free article] [PubMed]
19. Vesterby A, Gundersen HJ, Melsen F. Star volume of marrow space and trabeculae of the first lumbar vertebra: sampling efficiency and biological variation. Bone. 1989;10:7–13. [PubMed]
20. Hahn M, Vogel M, Pompesius-Kempa M, Delling G. Trabecular bone pattern factor--a new parameter for simple quantification of bone microarchitecture. Bone. 1992;13:327–30. [PubMed]
21. Hildebrand T, Rüegsegger P. Quantification of bone microarchitecture with the structure model index. Comput Meth Biomech Biomed Eng. 1997;1:15–23. [PubMed]
22. Majumdar S, Newitt D, Mathur A, Osman D, Gies A, Chiu E, Lotz J, Kinney J, Genant H. Magnetic resonance imaging of trabecular bone structure in the distal radius: relationship with X-ray tomographic microscopy and biomechanics. Osteoporos Int. 1996;6:376–385. [PubMed]
23. Feldkamp LA, Goldstein SA, Parfitt AM, Jesion G, Kleerekoper M. The direct examination of three-dimensional bone architecture in vitro by computed tomography. J Bone Min Research. 1989;4:3–11. [PubMed]
24. Sigmund E, Cho H, Song YQ. High - resolution MRI of internal field diffusion - weighting in trabecular bone. NMR in Biomedicine. 2009;22:436–448. [PubMed]
25. Chung H, Wehrli F, Williams J, Kugelmass S. Relationship between NMR transverse relaxation, trabecular bone architecture, and strength. Proceedings of the National Academy of Sciences. 1993;90:10250–10254. [PubMed]
26. Wehrli FW, Song HK, Saha PK, Wright AC. Quantitative MRI for the assessment of bone structure and function. NMR in Biomedicine. 2006;19:731–764. [PubMed]
27. Griffith JF, Yeung DK, Antonio GE, Wong SY, Kwok TC, Woo J, Leung PC. Vertebral Marrow Fat Content and Diffusion and Perfusion Indexes in Women with Varying Bone Density: MR Evaluation 1. Radiology. 2006;241:831–838. [PubMed]
28. Griffith JF, Yeung DK, Tsang PH, Choi KC, Kwok TC, Ahuja AT, Leung KS, Leung PC. Compromised bone marrow perfusion in osteoporosis. Journal of Bone and Mineral Research. 2008;23:1068–1075. [PubMed]
29. Saha PK, Chaudhuri BB. Detection of 3-D simple points for topology preserving transformations with application to thinning. IEEE Trans Pat Anal Mach Intel. 1994;16:1028–1032.
30. Saha PK, Chaudhuri BB. 3D digital topology under binary transformation with applications. Computer vision and image understanding. 1996;63:418–429.
31. Saha PK, Chaudhuri BB, Majumder DD. A new shape preserving parallel thinning algorithm for 3D digital images. Pat Recog. 1997;30:1939–1955.
32. Saha PK, Gomberg BR, Wehrli FW. Three-dimensional digital topological characterization of cancellous bone architecture. Int J Imag Sys Tech. 2000;11:81–90.
33. Lam SC, Wald MJ, Rajapakse CS, Liu Y, Saha PK, Wehrli FW. Performance of the MRI-based virtual bone biopsy in the distal radius: serial reproducibility and reliability of structural and mechanical parameters in women representative of osteoporosis study populations. Bone. 2011;49:895–903. [PMC free article] [PubMed]
34. Wehrli FW, Gomberg BR, Saha PK, Song HK, Hwang SN, Snyder PJ. Digital topological analysis of in vivo magnetic resonance microimages of trabecular bone reveals structural implications of osteoporosis. J Bone Min Res. 2001;16:1520–31. [PubMed]
35. Liu XS, Sajda P, Saha PK, Wehrli FW, Bevill G, Keaveny TM, Guo XE. Complete volumetric decomposition of individual trabecular plates and rods and its morphological correlations with anisotropic elastic moduli in human trabecular bone. J Bone Miner Res. 2008;23:223–35. [PubMed]
36. Jin D, Saha PK. A new fuzzy skeletonization algorithm and its applications to medical imaging. Proc of the Int Conf Image Anal Proc (ICIAP); Naples, Italy. 2013. pp. 662–671.
37. Liu Y, Jin D, Li C, Janz KF, Burns TL, Torner JC, Levy SM, Saha PK. A robust algorithm for thickness computation at low resolution and its application to in vivo trabecular bone CT imaging. IEEE Trans Biomed Eng. 2014;61:2057–2069. [PMC free article] [PubMed]
38. Saha PK, Wehrli FW, Gomberg BR. Fuzzy distance transform: theory, algorithms, and applications. Computer Vision and Image Understanding. 2002;86:171–190.
39. Saha PK, Strand R, Borgefors G. Digital topology and geometry in medical imaging: a survey. Medical Imaging, IEEE Transactions on. 2015;34:1940–1964. [PubMed]
40. Saha PK, Wehrli FW. Measurement of trabecular bone thickness in the limited resolution regime of in vivo MRI by fuzzy distance transform. Medical Imaging, IEEE Transactions on. 2004;23:53–62. [PubMed]
41. Chang G, Wang L, Liang G, Babb JS, Saha PK, Regatte RR. Reproducibility of subregional trabecular bone micro-architectural measures derived from 7-Tesla magnetic resonance images. MAGMA. 2011;24:121–5. [PMC free article] [PubMed]
42. Rajapakse CS, Magland JF, Wehrli FW. Fast prospective registration of in vivo MR images of trabecular bone microstructure in longitudinal studies. Magn Reson Med. 2008;59:1120–6. [PMC free article] [PubMed]
43. Recker RR. Architecture and vertebral fracture. Calc Tiss Int. 1993;53(Suppl 1):S139–142. [PubMed]
44. Wehrli FW, Ladinsky GA, Jones C, Benito M, Magland J, Vasilic B, Popescu AM, Zemel B, Cucchiara AJ, Wright AC, Song HK, Saha PK, Peachey H, Snyder PJ. In vivo magnetic resonance detects rapid remodeling changes in the topology of the trabecular bone network after menopause and the protective effect of estradiol. J Bone Miner Res. 2008;23:730–740. [PubMed]
45. Wehrli FW, Leonard MB, Saha PK, Gomberg BG. Quantitative high-resolution MRI reveals structural implications of renal osteodystrophy on trabecular and cortical bone. J Magn Reson Imag. 2004;20:83–89. [PubMed]
46. Liu XS, Stein EM, Zhou B, Zhang CA, Nickolas TL, Cohen A, Thomas V, McMahon DJ, Cosman F, Nieves J. Individual trabecula segmentation (ITS)-based morphological analyses and microfinite element analysis of HR-pQCT images discriminate postmenopausal fragility fractures independent of DXA measurements. J Bone Min Res. 2012;27:263–272. [PMC free article] [PubMed]
47. Baran DT, Bergfeld MA, Teitelbaum SL, Avioli LV. Effect of testosterone therapy on bone formation in an osteoporotic hypogonadal male. Calcified tissue research. 1978;26:103–106. [PubMed]
48. Francis RM, Peacock M, Aaron J, Selby P, Taylor G, Thompson J, Marshall D, Horsman A. Osteoporosis in hypogonadal men: role of decreased plasma 1, 25-dihydroxyvitamin D, calcium malabsorption, and low bone formation. Bone. 1986;7:261–268. [PubMed]