2D-3D reconstruction of the proximal femur from DXA scans: Evaluation of the 3D-Shaper software

Introduction: Osteoporosis is currently diagnosed based on areal bone mineral density (aBMD) computed from 2D DXA scans. However, aBMD is a limited surrogate for femoral strength since it does not account for 3D bone geometry and density distribution. QCT scans combined with finite element (FE) analysis can deliver improved femoral strength predictions. However, non-negligible radiation dose and high costs prevent a systematic usage of this technique for screening purposes. As an alternative, the 3D-Shaper software (3D-Shaper Medical, Spain) reconstructs the 3D shape and density distribution of the femur from 2D DXA scans. This approach could deliver a more accurate estimation of femoral strength than aBMD by using FE analysis on the reconstructed 3D DXA. Methods: Here we present the first independent evaluation of the software, using a dataset of 77 ex vivo femora. We extend a prior evaluation by including the density distribution differences, the spatial correlation of density values and an FE analysis. Yet, cortical thickness is left out of this evaluation, since the cortex is not resolved in our FE models. Results: We found an average surface distance of 1.16 mm between 3D DXA and QCT images, which shows a good reconstruction of the bone geometry. Although BMD values obtained from 3D DXA and QCT correlated well (r 2 = 0.92), the 3D DXA BMD were systematically lower. The average BMD difference amounted to 64 mg/cm3, more than one-third of the 3D DXA BMD. Furthermore, the low correlation (r 2 = 0.48) between density values of both images indicates a limited reconstruction of the 3D density distribution. FE results were in good agreement between QCT and 3D DXA images, with a high coefficient of determination (r 2 = 0.88). However, this correlation was not statistically different from a direct prediction by aBMD. Moreover, we found differences in the fracture patterns between the two image types. QCT-based FE analysis resulted mostly in femoral neck fractures and 3D DXA-based FE in subcapital or pertrochanteric fractures. Discussion: In conclusion, 3D-Shaper generates an altered BMD distribution compared to QCT but, after careful density calibration, shows an interesting potential for deriving a standardized femoral strength from a DXA scan.


1
Materials and Methods

Neural network for QCT segmentation
The neural network used for the QCT segmentation is a 3D version of the well-known U-Net architecture (Ronneberger et al., 2015). With respect to the original model, the number of convolutional layers was reduced from 23 to 14 to decrease the number of trainable parameters. The detailed architecture is illustrated in the Supplementary Figure S1. For the convolutional layers, the original kernel size (3x3) was extended to 3D (3x3x3). The convolution layers were padded to preserve the image size and followed by a rectified linear unit (ReLU). The max pooling layers had a kernel size 2x2x2 and stride 2 to halve the image size in all three dimensions. A dropout layer with rate 0.5 was introduced in the middle of the network to prevent overfitting. The upsampling layers had a size 2x2x2 to double the image size in all three dimensions. The output layer consisted of a final convolutional layer with kernel size 1x1x1, followed by a sigmoid activation function.
The model was trained with image-segmentation pairs for 80 epochs, using the dice coefficient as loss function and the Adam algorithm as optimizer (Kingma and Ba, 2017).

Implicit coordinate system of the proximal part of the femur
In order to transform the mask of the femur to the orientation needed for the FE simulations an implicit coordinate system similar to the one by Chandran et al. (2017) was introduced.
Initially four points were placed on the surface of the femoral head by approaching and intersecting a plane in four different directions. An initial sphere was then fitted based on the four points. The center of the so-found sphere was used to define a polar coordinate system, which was the basis to create a set of 200 rays pointing from the center to the surface of the sphere. The intersection of the 200 rays with the sphere surface were recorded and used to refine the fit of the sphere. This way the femoral head center and head radius were found. In a next step a series of progressively dilated spheres around the femoral head were intersected with the bone mask. The centers of mass of the consecutive intersections were computed and the linear regression of the so-found centers of mass resulted in a first crude estimate of the femoral neck axis. The initial estimate was refined by fitting a line to the centers of mass of the intersections between the bone mask and a series of disks perpendicular to the initial neck axis.
To standardize the proximal shaft axis detection without precise knowledge of the distal cut location, a two-step procedure was adopted. Approaching in infero-superior direction, the number of pixels was computed in each section, and the first "complete" section of the shaft was defined as the one with less than 5% variation in number of pixels with respect to the former. The centroids of the slices within 1cm above the first complete section were used to fit an initial shaft axis. A closest point between the (generally not intersecting) femoral neck and shaft axes was computed and used to correct both axes.
The neck and proximal shaft axes were used to estimate the location of the lesser trochanter (LT). The LT surface was then segmented based on a measure of local curvature of the slice-wise contoured bone and the peak of the LT was placed at the slice with the largest curvature. The section of the shaft extending distally from the base of the LT over 15 mm (scaled with the head radius to account for inter-subject variability) was used as standard location for the estimation of the proximal femoral shaft axis. Centroids of each slice in this section were computed. In addition, the slice centroid of the section with the smallest eccentricity above the LT was detected. Two candidate guesses for the shaft axis were defined based on this information. The first guess was based on a linear regression of the computed slice centroids. The second guess was based on a RANSAC fit (Chandran et al., 2017). Based on the assumption that the initial estimate of the shaft axis was reasonably good, the candidate closer to the initial guess was kept. In cases where the distal length was not sufficient to adopt the described procedure, the initial shaft axis was kept or replaced by the axis defined by the most distal slice centroid and the centroid of the section with the smallest eccentricity above the LT.
The closest point between the neck axis and the so-found shaft axis were computed again. The soobtained closest point in conjunction with the femoral head center and the slice center of the most distal slice were used to correct the femoral neck and shaft axis. The definition of the implicit coordinate system can thus be reduced to three points (head center, neck-shaft intersection point and most distal slice center).
A formal validation of this procedure was difficult to obtain and could thus not be performed. However, we tested the procedure on about 2600 3D DXA images and around 200 CT images and it worked reliably as assessed by visual control.

Stance configuration
Femoral strength is mostly computed in a configuration mimicking a side fall. We modeled also a more physiologic stance configuration like in Dall' Ara et al. (2013).
The mask of the bone was transformed rigidly so that the proximal shaft axis formed an angle of 20° with the vertical axis. An embedding cap was then applied on the femoral head with a depth of three to four voxels. The femur was cut at a distance corresponding to one head radius distal to the lesser trochanter. All nodes on the distal cut surface were fully constrained. A vertical displacement was applied to a node at the center of the head, which was kinematically coupled to the top surface of the embedding cap (Panyasantisuk et al., 2018). Supplementary Figure S2 shows the stance configuration with the boundary conditions.
Vertical displacements and reaction forces acting on the driving node in the femoral head center were recorded. In opposition to the fall load case a clear maximum in the force-displacement characteristic was reached in a large part of the cases in stance configuration. In these cases, strength was defined as the reaction force at the peak of the force-displacement curve. Otherwise, strength was defined as the reaction force when a displacement corresponding to 6% of the femoral head radius was reached.
After a first run where meshes based on QCT and 3D DXA images were treated independently, two different corrections for BMC differences between QCT and 3D DXA were introduced as for the side fall condition. A collective correction based on the regression equation between the BMC values of 3D DXA-based and the QCT-based meshes, and an individual correction where the BMC of each 3D DXA-based mesh was adjusted to the corresponding QCT-based mesh As for the side fall configuration the spatial distribution of damage (3D damage maps) was extracted for the QCT-based simulations and the individually corrected 3D DXA-based simulations, and an average damage distribution was computed.
The results of the simulations in stance configuration are shown in Supplementary Figure S4, Supplementary Figure S5 and Supplementary Table S2.

Constitutive model
The material model used for the FE simulations performed in this paper is an isotropic formulation of the model developed by  with yielding defined by the criterion described in .
The elastic stiffness tensor is given by Where 0 and 0 are the Lamé constants with the positive-definiteness condition 3 0 + 2 0 > 0 0 > 0 Accordingly, the compliance tensor is given by Where 0 is the reference Young's modulus and 0 the Poisson ratio. The relationship between the Lamé constants and the engineering constants is given by: The tissue function modifies the power law for trabecular bone typically represented by low density elements, so that realistic stiffness values are obtained for high density elements that typically contain a large fraction of cortical bone. The piecewise relationship is illustrated in Supplementary Figure S3.
The quadric yield criterion  coupled with a perfectly plastic post-yield behavior is given by: The origin, orientation and shape of the criterion are defined by the fourth-order tensor and second order tensor . Where 0 + and 0 − are the tensile and compressive yield stresses. In the isotropic case, yielding through shearing is given by 0 : Where 0 is the so-called interaction parameter .
The total Green-Lagrange strain is decomposed in its elastic and plastic parts.

= +
Stiffness is isotropically reduced by a damage variable ( ) resulting in the following equation for the stress: The isotropic damage variable ( ) is driven by the cumulated plastic strain : Where is given by:

= ∫ ‖̇‖
The numerical values of the constants used are given in Supplementary Table S1.  Figure S1: Network architecture of the QCT segmentation model.