To further investigate the potentiality of the MCP, the present work considers three different sets of trabecular bone 3D image samples: two sets from distal radius, one of them containing 15 ex vivo μCT samples, and the other one containing 103 in vivo MRI samples; the third one containing 29 ex vivo μCT L3 vertebral samples. The final isotropic resolutions are 34 μm to the μCT and 90 μm to the MRI images, and the main analyzed direction was the axial one (craniocaudal to the vertebrae and distal-proximal to the radius).

The μCT distal radius samples, with lateral size 12 mm, were harvested with a mean distance of 9.75 mm from the distal extremity, and volumes of interest (VOI) were selected with sizes which vary according to the material’s clinical analysis. They were imaged with the scanner microCT-20 (Scanco Medical, Brüttisellen, Switzerland) and, to the noise removal, the μCT 3D images were filtered with a Gaussian 3D filter. In each case, the grayscale histogram of the filtered images has two peaks, corresponding to marrow and bone; so, they were binarized using a global threshold equal to the minimum between the two peaks. The 15 image sets have 239 slices each, with 2D ROIs 212 × 212, 237 × 237, 242 × 242, 252 × 252 e 257 × 257 pixels; the 10 other samples have 268 × 268 pixels. Additional details concerning to the sample’s preparation and acquisition protocols are described in [27].

A set of 29 μCT vertebral samples were supplied by the Department of Forensic Medicine, Jagiellonian University Medical College. The specimens were taken from female individuals without metabolic bone disease or vertebral fractures. Mean and standard deviation of the individuals age were equal to years, respectively. Immediately after dissection, all soft tissue was cleaned out and the samples were placed in containers filled with ethanol. An X-tek Benchtop CT160Xi high-resolution CT scanner (Nikon Metrology, Tring, UK) was used to scan the vertebral bodies. The images were segmented into bone and marrow cavity phases with a global thresholding method. The segmentation threshold was selected automatically based on the MaxEntropy algorithm [26], such that the information entropy consistent with a two-phase model be maximal. The final 3D binarized images have size that vary from 770 until 1088 pixels in x, from 605 until 876 pixels in y and from 413 until 713 slices (z direction), being the size average .

The elasticity study with these vertebral samples was performed by mechanical test. An MTS Mini Bionix 858.02 loading system with a combined force/torque transducer with range of 25 kN/100 N.m was used to perform the compression tests. The specimens were located between two stiff steel plates which were firmly mounted to the force/torque transduced and to an upper jaw of the loading system. Prior to mechanical testing each probed specimen was glued with a self-curing denture base acrylic resin between two polycarbonate sheets at its endplate surfaces. This procedure was chosen to create two surfaces which will be as parallel as possible above each endplate to transmit the compressive load from the loading system to each specimen in an uniform way. The polycarbonate sheets were removed from the vertebra endplates before the testing. Each vertebra was loaded in compression with a loading rate of 5 mm/min to a certain level of engineering deformation (at most of the original height of the specimen). The compressive force was monitored during the test with sample rate of 20 Hz. All data that were measured during the compression tests were transformed to plots of applied force and displacement for each specimen. Compliance of the loading system was measured as well, so, during the post-processing, it was possible to gain a true relation between an applied force and deformation of a vertebra body. The stiffness in the linear part of a loading path for each specimen was evaluated and the Young modulus, E, was defined as the ratio of the product of the stiffness and the vertebral height to the mean cross section area of the vertebral body.

A set of 103 MRI radius samples were considered from the distal metaphysis and from a group including healthy subjects and a mix of disease stages. The MRI acquisitions were performed in a 3 Tesla system and scanned in 3D using a T1-weighted gradient echo sequence (TE/TR/a=5 ms/16 ms/25 Âº). The MRI images were acquired with a nominal isotropic resolution of 180 μm. MR image processing and analysis were performed with MATLAB R2012a (The MathWorks, Inc., Natick, MA). The image preparation steps consisted of an initial segmentation using a rectangular region of interest, image intensities homogeneity correction, interpolation and binarization. All the steps were applied as in [2], with the exemption of the interpolation, which was performed by applying a 3D non-local upsampling algorithm, achieving final resolution of 90 μm [29]. It has 65 samples with 80 slices, 10 with 120 and, the other ones, vary from 30 until 200 slices, predominantly between 50 and 100. Each 2D image has laterals dimensions varying from 38 up to 206 pixels, predominantly around 70 × 100 pixels.

Finite element method simulations were conducted to estimate Young modulus in all the 103 distal radius samples as well as for the 15 μCT distal radius samples. For that, a mesh was created based on the 3D trabecular bone images using an optimized algorithm [1] implemented in Matlab R2011a, which converts each voxel to an hexahedron element (brick element). Compression stress-strain tests were numerically simulated by a finite element linear-elastic-isotropic analysis performed in Ansys v11.0 (Ansys Inc., Southpointe, PA). The bulk material properties were set to , a common value assumed to compact bone, and Poisson’s coefficient . A deformation of 1 % of the edge length was imposed in all the distal radius compression simulations. Computational cost of the simulations was approximately of 5 h per sample on a computer workstation (Quad Core at 2.83 GHz and 8 GB of RAM). After applying the homogenization theory [23], apparent Young modulus results were obtained.

The computations of , CEP and τ values were done using OsteoImage, a computer program developed by one of the authors especially to TB image analyses. The statistical analyses were performed with the free software RGui [34] and the 3D image reconstructions were done with ImageJ (http://​rsbweb.​nih.​gov/​ij/​).

The TB volume fraction, , represents the quantity of TB content present in the sample volume and is obtained by the ratio:
where V _trab is the trabecular volume and V _total is the total sample volume. From a 3D binary image sample, the TB volume fraction may be computed using the number of voxels representing the trabecular bone and the total volume is the number of voxels of the whole sample.

The trabecular network connectivity can be inferred by the Euler-Poincaré characteristic, EPC, which can be estimated by automatic counting of isolated parts, I, redundant connections, C, and closed cavities, H [47]:
As the trabeculae have no closed cavities [18] and the number of isolated parts is approximately 1 in a well structured sample, the EPC value should be negative and the lower the value the higher the connectivity [8]; in this case, the connectivity is estimated by its modulus. A positive EPC value indicates that the sample has more isolated parts than connections, and, therefore, the EPC indicates that its structure has lost much of its connectedness.

As EPC is a zero-dimensional measure, it needs to be estimated by a three-dimensional test; for practical purposes, a couple of parallel 2D images can be used, forming a disector [21, 35, 43, 47], and the EPC can be estimated for each one of them inside the volume of interest. In general, the EPC is given normalized by its volume size, EPC _V . The algorithm to compute the EPC can be seen in [36].

The tortuosity, τ, characterizes how much an object departures from being straight and this concept has been extended to the trabecular bone network. Geometrically, it is defined as
where L _G is the geodesic distance between two connected points, say a and b, of the trabecular network without passing across other phases (marrow cavity); and L _E is the Euclidean distance between these points, which will be considered here as the distance between two parallel reference planes (see Fig 2) [50]. This approach allows to classify as tortuous, , any filamentous structure that is not perpendicular to the reference planes.

Gommes et al. [19] proposed a geodesic reconstruction (GR) algorithm that can be applied on binary images to estimate the geodesic length. This algorithm was implemented in a previous work [38] and was used to the solid phase of the bone samples, sweeping the image along the reference plane direction, reconstructing the trabecular bone network voxel by voxel. The number of GR necessary to recover all the trabeculae of an image depends on their sinuosities, exceeding the number of analyzed slices considered as the Euclidean distance; the equality occurs only in the case of a structure completely perpendicular to the sweeping direction.

During the GR process, the algorithm computes and stores the Euclidean, L _E , and the geodesic, L _G , lengths. A distribution of Euclidean and geodesic lengths is generated. Taking the geodesic distance average, at each Euclidean distance, the tortuosity can be estimated as the slope of the best fit line of points . This algorithm can be applied directly to 3D binarized μCT or MRI images. More details of the algorithm implementation can be found in [38, 40].

The elasticity is an important property of a material because it reflects its stiffness and flexibility when subject to load. Imposing an uniaxial strain ε to the sample, it is related with the stress σ as follows
where E is the Young modulus of elasticity. Usually, σ is obtained from the sample reaction force, divided by the area where it is being applied on. Rigorously, the trabecular structure is not isotropic [22, 44, 45], hence E is not a scalar, but a symmetric tensor; nevertheless, considering the complexity of modeling a porous structure, an isotropic model can be reasonably assumed [1, 14].

The 3D trabecular bone images were meshed to the elastic simulation using an optimized algorithm [1] implemented in Matlab R2011a (The MathWorks Inc., Natick, MA) which converts each voxel to an hexahedron element (brick element). Compression stress-strain test in each space direction was numerically simulated by a finite element linear-elastic-isotropic analysis performed in Ansys v11.0 (Ansys Inc., Southpointe, PA). The bulk material properties were set to , a common value assumed to compact bone, and Poisson’s coefficient . A deformation of 1 % of the edge length was imposed in all the compression simulations. Computational cost of the simulations was approximately of 5 h per sample on a computer workstation (Quad Core at 2.83 GHz and 8 GB of RAM). After applying the homogenization theory [23], apparent Young modulus results were obtained in each spatial direction (E _x , E _y , E _z ).

Table 1 presents the mean and standard deviation (SD) that were obtained for the distal radius μCT and MRI trabecular bone cohorts. Firstly, by a simple inspection of the data in Table 1, it is observed that in the z direction τ has the lowest mean and SD, and the E has the highest value ones, in both groups. This corresponds to the distal-proximal direction, which is normally the direction that is more frequently submitted to tensile and compressive forces, when compared to the x and y ones, corresponding to the horizontal sweeping directions. This evidence is an indication that the trabeculae get aligned to turn the structure stronger, which is in agreement with the very well known fact that the trabecular bone aligns in the direction which it is more frequently mechanically demanded [20, 45, 49].