- Email: [email protected]
A novel methodology for generating 3D finite element models of the hip from 2D radiographs
Journal of Biomechanics 47 (2014) 438–444
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
A novel methodology for generating 3D finite element models of the hip from 2D radiographs Jérôme Thevenot a, Janne Koivumäki a, Volker Kuhn b,c, Felix Eckstein b,d, Timo Jämsä a,e,n a
Department of Medical Technology, University of Oulu, Oulu 90014, Finland Institute of Anatomy, Ludwig-Maximilians-University, Munich, Germany c Department of Trauma Surgery and Sports Medicine, Medical University, Innsbruck, Austria d Institute of Anatomy and Musculoskeletal Research, Paracelsus Medical University, Salzburg, Austria e Department of Diagnostic Radiology, Oulu University Hospital, Oulu, Finland b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 6 November 2013
Finite element (FE) modelling has been proposed as a tool for estimating fracture risk and patientspecific FE models are commonly based on computed tomography (CT). Here, we present a novel method to automatically create personalised 3D models from standard 2D hip radiographs. A set of geometrical parameters of the femur were determined from seven a–p hip radiographs and compared to the 3D femoral shape obtained from CT as training material; the error in reconstructing the 3D model from the 2D radiographs was assessed. Using the geometry parameters as the input, the 3D shape of another 21 femora was built and meshed, separating a cortical and trabecular compartment. The material properties were derived from the homogeneity index assessed by texture analysis of the radiographs, with focus on the principal tensile and compressive trabecular systems. The ability of these FE models to predict failure load as determined by experimental biomechanical testing was evaluated and compared to the predictive ability of DXA. The average reconstruction error of the 3D models was 1.77 mm (7 1.17 mm), with the error being smallest in the femoral head and neck, and greatest in the trochanter. The correlation of the FE predicted failure load with the experimental failure load was r2 ¼64% for the reconstruction FE model, which was significantly better (po 0.05) than that for DXA (r2 ¼ 24%). This novel method for automatically constructing a patient-specific 3D finite element model from standard 2D radiographs shows encouraging results in estimating patient-specific failure loads. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Hip fracture Biomechanics Finite element method Structure analysis Radiography
1. Introduction The clinical diagnosis of osteoporosis and eventual fracture risk is based on bone mineral density (BMD) assessed from dual energy X-ray absorptiometry (DXA). However, DXA-based BMD has been shown to be insufficient for predicting the individual risk of fracture (Kanis, 2002; Schuit et al., 2004); this can be explained by the fact that the mechanical strength of the hip is highly dependent on the shape and structure of the bone (Beck et al., 1990; Gnudi et al., 2002; Partanen et al., 2001; Pulkkinen et al., 2006), and fractures also occur in individuals with non-osteoporotic BMD (Cummings and Melton, 2002; Meunier, 1993; Schuit et al., 2004). The finite element (FE) method has shown promising results for the clinical assessment of hip strength (Cody et al., 1999). While some studies were concentrated on 2D FE models based on
n Corresponding author at: Department of Medical Technology, University of Oulu, Oulu 90014, Finland, Tel.: þ 358 29 448 6001. E-mail address: timo.jamsa@oulu.fi (T. Jämsä).
0021-9290/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jbiomech.2013.11.004
DXA images (Testi et al., 2002, 2004), most of them used 3D models derived from computed tomography (CT) scans to estimate hip strength (Bessho et al., 2004; Keyak and Rossi, 2000; Schileo et al., 2008). Unfortunately, building patient-specific 3D FE models requires the use of high-cost volumetric imaging coupled with high radiation doses and computational power. A solution to the time-consuming generation of 3D models from CT scans is the use of active shape modelling (ASM) (Gunay et al., 2007) or statistical shape modelling (Schumann et al., 2010), both of which have proven to be sufficiently accurate for different applications. Briefly, these methods will iterate an initial 3D object to fit a 2D input contour or geometrical parameters (Baker-LePain et al., 2011; Bryan et al., 2009; Gregory et al., 2004; Langton et al., 2009; Schumann et al., 2010; Zheng et al., 2009). However, these methods usually use contours extracted from 2D images and do not separate cortical and trabecular bone even though cortical thickness has shown to be an important predictor of hip fracture risk (Augat et al., 1996; Chappard et al., 2010; Gnudi et al., 2002; Michelotti and Clark, 1999; Pulkkinen et al., 2004, 2008; Szulc et al., 2006).
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
The aim of this study is to provide a new automatic method to generate 3D finite element models from standard 2D radiographic pictures. The method will create both cortical and trabecular bone compartments from a set of geometrical parameters as an input. This paper describes the methodology for generating the models and presents the results in estimating experimental failure load.
2. Study sample and imaging For the present study, we used data of twenty-eight cadaver femora obtained from a larger experimental study at the Institute of Anatomy at the Ludwig Maximilians University of Munich (Germany) (Eckstein et al., 2004; Pulkkinen et al., 2006). Seven of the bones were used for training, and twenty-one for validation. The preparation and storage of the samples has previously been described in detail (Eckstein et al., 2004). After the dissection, excision and cleaning of the surrounding soft tissue, the proximal femora were radiographed and CT-scanned. Radiographs were taken using a Faxitron X-ray system (Model 43885A; Faxitron, Hewlett Packard, McMinnville, OR) and the X-ray films were digitised together with a calibration scale using a scanner at 600 dpi. CT scans were performed using a 16-row Multi-detector CT (MD-CT) scanner (Sensation 16; Siemens Medical Solutions, Erlangen, Germany). MD-CT and DXA scan setups have been previously presented (Koivumäki et al., 2012). For CT scans, a high-resolution protocol with a slice thickness of 0.75 mm was used. The settings were 120 kVp and 100 mAs, a 512 512 pixels image matrix, and a field of view of 100 mm. The in-plane spatial resolution was 0.25 mm 0.25 mm. A calibration phantom (Osteo Phantom, Siemens, Erlangen, Germany) composed of hydroxyapatite (0 mg/cm3 for water and 200 mg/ cm3 for bone) was placed below the specimens. In vitro DXA scans of the femora were obtained with a standard narrow fan beam scanner (GE LunarProdigy, GE Lunar, Madison, WI) with the femoral specimens submerged in a water bath. Standard positioning was used across all specimens, and the proximal femoral BMD values were evaluated with the software provided by the manufacturer. The cortical and trabecular bone compartments were segmented from the CT scans as previously presented (Koivumäki et al., 2012) using Mimics (v12.1, Materialise, Leuven, Belgium), with the minimum thickness for the cortical bone adjusted to 1 mm. A maximum element size of 3 mm was applied. The femurs were mechanically tested simulating a fall on the greater trochanter, as previously presented (Eckstein et al., 2004; Pulkkinen et al., 2006). Briefly, femora were positioned at 101 from horizontal plane and the neck at 151 internal rotation; loads were then applied to the greater trochanter through a pad at a rate of 6.6 mm/s. The failure load was determined from the loaddeformation curve.
439
3.1. Generation of the 3D shape 3.1.1. Geometrical parameters as an input A limited set of geometrical parameters was used to generate the volumetric shape of the model. The measurements of the geometrical parameters were done using Image Tool software (version 3.00; University of Texas Health Science Center, San Antonio, TX), as previously described (Partanen et al., 2001; Pulkkinen et al., 2006). Femoral neck axis length (FNAL), neck-shaft angle (NSA), trochanteric width (TW), femoral head diameter (HD), femoral neck diameter (ND), femoral shaft (FSD) diameter, femoral shaft cortex (FSC), and calcar femoral cortex width (CFC) were measured, as shown in Fig. 1.
3.1.2. Building up the 3D shape In order to establish relationships between the geometrical parameters and the skeleton of the femurs, seven bones with a large range of geometrical parameters were selected as CT-based 3D training material. The impact of each geometrical parameter on the volumetric shape of the bone was evaluated to delimitate the outer shape of the bone by points. Curves were generated from the interpolation of these points to represent the skeleton of the bone and surfaces were eventually derived from them (Fig. 2). The femoral head was considered a half-sphere (Schumann et al., 2010; Thevenot et al., 2009). Specific attention was given to the shape of the femoral neck since this reflects the risk of cervical fractures (Pulkkinen et al., 2006). The non-circular cross-section of the femoral neck (Zebaze et al., 2005) was estimated by using the
3. Construction of volumetric finite element model The generation of the model can be divided into three parts: (1) generation of the 3D shape, (2) meshing, and (3) determination and implementation of the material properties. The algorithm was developed under Matlab (version 7.1, MathWorks, Natick, MA). The output of the algorithm includes three command script files corresponding to each part of the model generation. The files were then imported into the pre- and post-processing software Femap (version 9.2, UGS Corp., Plano, TX) and launched consecutively in order to generate the model.
Fig. 1. Geometrical parameters. (A) Neck-shaft angle, (B‐C) and (E‐F) femoral shaft and calcar femoral cortical thickness, respectively, (B‐D) femoral shaft diameter, (E‐ G) trochanteric width, (I‐H) femoral neck axis length, (J‐K) and (L‐M) head and neck diameters.
440
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
Fig. 2. Building up the shape of the upper femur: (A) femoral head and main axis, (B) skeleton of the neck area, (C) neck surfaces based on previous curves and (D) skeleton of the trochanteric area.
Fig. 3. Left: generated model with the principle tensile system (PTS) and principal compressive system (PCS) shown in red. Right: Typical distribution of material properties, dark values corresponding to low density and light to high density.
training material. The neck area was divided into 10 segments along the femoral neck axis, and for each of these segments, the radius of the neck was checked from CT data in every 451 rotational direction. A relationship was established between each radius and FNAL, HD, ND and FSD, to obtain accurate shape estimation for an individual femoral neck. 3.2. Meshing A smoothed solid was created from the surfaces of the model and then meshed with tetrahedral elements with a size increasing progressively from the surface to the centre of the bone to a maximum element size of 3 mm. The size of the elements near the surface of the bone was chosen to be proportional to the FSC, being less than half its size. The faces of elements on the surface were extruded with a thickness of 1 mm, corresponding to the selected minimum thickness of the cortical bone and representing the outer contour of the bone. Cortical and trabecular compartments were separated during the meshing process. The inner contour of the cortical bone was delimited from the elements within the cortex as defined by FSC and CFC. FSC was represented by two consecutive layers of elements from the outer surface in the shaft area, whereas CFC was defined by elements within the lower neck area. Finally, the cortical bone was re-meshed using 10-noded tetrahedral elements.
The remaining trabecular bone was re-meshed using 4-noded tetrahedral elements. The size of the elements along the model varies between 1 and 3 mm, with the smaller elements being close to the inner and outer surfaces of the cortical bone. The maximum size of the elements was chosen as previously suggested (Keyak and Skinner, 1992). 3.3. Material properties Distribution of the material properties was evaluated from trabecular structure analysis of the radiographic picture. The spatial distribution of the material properties in the trabecular bone was related to the principal tensile and compressive systems (PTS and PCS, respectively) (Rudman et al., 2006), representing the highest mineralised areas of the trabecular bone. Here, the Young's moduli were assessed using the homogeneity index (HI), adjusted to the material properties derived from Hounsfield Units (HU) of CT scans, with the relationship being based on the seven bones used as training material. 3.3.1. Localisation of the principal tensile and compressive trabecular systems The general localisation and thickness of the PTS was assessed in 2D from the radiographs and its 3D orientation from the
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
441
CT-scans, using the training material as reference. A main curve representing the PTS from the greater trochanter to the femoral head was generated in the model. Similarly, the PCS was approximated to go from the lower neck area to the top of the femoral head and was represented as a curved frustum (Fig. 3, left).
homogeneity index (HI). The first column and first row of the GLCM were cropped to remove the values corresponding to empty spaces in the bone, and to consequently evaluate the homogeneity of the mineralised trabecular bone. HI, as a measure of the closeness of the distribution of elements near the diagonal of the GLCM, was calculated as follows:
3.3.2. Structural analysis of the radiographs The image processing procedure presented in Fig. 4 was performed as previously reported (Thevenot et al., 2013). It was performed with an algorithm developed under Matlab with a graphical interface allowing the user to select a specific region of interest (ROI) (15 mm 15 mm). In image pre-processing, noise was removed by using a 3 3 median filter, followed by morphological top and bottom hat operations. A gradient-based grey-level image was then constructed, which enabled the visualisation of the trabecular structure even under varying image contrasts and brightness conditions of the radiograph (Pulkkinen et al., 2005). The ROI was chosen within the fibres of the PCS (inferomedial neck, at the neck/head transition area) and the trabecular main orientation (TMO) was then evaluated using Fourier transform. An empirically-derived threshold of 0.03 was used to remove nonsignificant frequencies in the Fourier-based image, and linear regression was applied to fit a line g(x)¼ axþb to the threshold values in order to extract the TMO in the frequency domain. Eventually, the angle α between the TMO and the shaft axis in the space domain was calculated and the ROI rotated along it. A greylevel co-occurrence matrix (GLCM) with a distance of one pixel was calculated along the TMO for determination of the
GLCMði; jÞ HI ¼ ∑ i;j 1 þ i j
Selection of the region of interest (ROI)
Filtering of the picture
Gradient-based transform
Calculation of trabecular main orientation (TMO) using Fourier transform (angle α)
Rotation α of the ROI: alignment along the TMO
Gray-level co-occurrence matrix (GLCM) along the TMO
The evaluated radiograph-based HI of the trabeculae-rich lower femoral neck area was used to derive the Young's modulus along the FE model. 3.3.3. Implementation of the material properties Young's modulus of 14.2 GPa was applied for the cortical bone (Viceconti, 1998). An intermediate Young's modulus of 7 GPa was used to the first layer of trabecular elements to simulate the transition between cortical and trabecular bones. A Poisson's ratio of 0.33 was selected for both cortical and trabecular bone (Lengsfeld et al., 1998). CT scan data of the training material were used to establish the relationship between the calculated HI and the material properties of the model. An average grey value of all of the voxels inside an element was calculated in Mimics (Koivumäki et al., 2012). The bone equivalent density (ash density, ρash) was then defined by assuming a linear relationship in which the density is proportional to the attenuation (Koivumäki et al., 2012). The Young's modulus (E, MPa) was estimated using the equation E ¼10.095ρ, where ρ is the CT density value (mg/cm3) (Duchemin et al., 2008). Yield stress was calculated using the previously proposed equations, 114ρash1.72 and 137ρash1.88, for cortical and trabecular bone, respectively (Keller, 1994; Keyak et al., 1994). The curve representing the trabecular PTS orientation was divided into 20 sub-curves and elements surrounding them within a consistent estimated PTS thickness were affected with a specific Young's modulus. An extra transition layer of elements with reduced Young's modulus was generated. According to the CT data, the Young's modulus increases from the greater trochanter to the femoral head area along the tensile system. Here, this increase was estimated by linear increase of the modulus through the PTS, with each sub-area representing different material properties. A similar method was applied for the trabecular PCS. Here, the higher Young's modulus values were located in the lower neck and the lower values in the upper part of the femoral head, differences between which were estimated using a linear distribution. The material properties of trabecular elements outside PTS and PCS were grouped depending on their location. Homogeneous material properties were affected separately for shaft, trochanteric, femoral neck, and femoral head areas, derived from the HI. As a result, a total of 66 material properties were created: - Cortical bone (1) and the transition layer between cortical and trabecular bone (1). - PTS (20) and the transition layer (20). - PCS (10) and the transition layer (10). - Trabecular areas outside PTS and PCS (shaft, trochanteric, femoral neck, and femoral head) (4).
GLCM reduction 3.4. Finite element analysis of failure load
Calculation of the homogeneity index based on the GLCM Fig. 4. Flowchart of the image processing to obtain the homogeneity index of the trabecular structure.
After applying the boundary conditions to simulate the experiment (Eckstein et al., 2004), bilinear elastoplastic FE analysis was performed, as previously presented (Koivumäki et al., 2012), using the Newton–Raphson method and Drucker–Prager yield criterion (Drucker and Prager, 1952), applying a post-yield modulus 5% of
442
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
the Young's modulus (Bayraktar et al., 2004), and an ultimate tensile stress presumed to be 0.8 times the compressive yield stress (Bessho et al., 2007). The fracture was determined to occur when a surface cortical element failed (Bessho et al., 2007), failure in tension was characterised for a maximum principal stress higher than the ultimate tensile stress and failure in compression for a minimum principal strain lower than 7300 microstrain (Koivumäki et al., 2012). The failure load predicted from the finite element analysis was compared with the experimentally measured failure load. 3.5. Statistical analysis
5. Discussion
Linear regression analysis was used to evaluate the relationship between the simulated and the experimentally measured fracture loads, and between BMD and BMC and the experimental fracture load. One-tailed Fisher z-transformation was used to test the hypothesis that the simulated fracture load had a higher correlation with the experimental fracture load than DXA. Statistical analyses were performed using SPSS software (SPSS 16.0 for Windows, SPSS Inc. Chicago, IL). A p-value o0.05 was considered to be statistically significant.
4. Results The present methodology generated 3D FE models of upper femur from 2D radiographs with visually similar appearances to the CT-based models. The average reconstruction error for the whole upper femur was 1.77 mm ( 71.17 mm), with the error being lowest in the femoral head (1.03 mm 70.45 mm) and femoral neck (1.27 mm 70.60 mm) and greatest in the trochanteric area (2.24 mm 7 1.33 mm) (Table 1, Fig. 5). The relationship between the estimated and the experimental failure load is presented in Fig. 6. The estimated failure load values Table 1 Average reconstruction error for different upper femur areas of randomly selected validation models (N ¼ 7). Data presented as mean 7 SD. Femoral head Femoral neck Trochanter Shaft area Error (mm) 1.03 7 0.45
were well correlated with the experimental ones (r2 ¼64%), SEE¼ 543 N, po 0.001). The slope was 1.20 (not significantly different from 1, p ¼0.353), with an intercept of 472 N (not significantly different from 0, p¼ 0.477). The coefficients of determination between BMD and BMC, and the experimental fracture load were r2 ¼24% (N ¼ 20, SEE¼ 793 N, p¼ 0.028) and r2 ¼35% (N ¼20, SEE ¼734 N, p ¼0.006), respectively. The coefficient of determination was significantly higher for the FE model when compared to BMD (p o0.05) but not when compared to BMC (p ¼0.11).
1.277 0.60
Total
2.247 1.33 1.577 0.79 1.777 1.17
In the present study, we developed a novel method to obtain a 3D finite element model based on a standard 2D radiographic picture. The 3D model is built up automatically using the geometric shape of the femur derived from simple measurements of geometrical dimensions. The heterogeneous distribution of material properties in the model is based on the trabecular structure analysis of the radiograph with a focus on the orientation and location of the principal tensile and compressive systems. This method appeared able to construct the patient-specific 3D shape with feasible accuracy, and the FE analysis showed promising results in estimating the individual experimental failure load. While most studies building 3D models from 2D data use the ASM method from the outer contour of the bone, we decided to create a method that also defines the inner contour of the cortical bone. It has been shown that the failure load can be predicted with reasonable accuracy from radiographic texture analysis and geometric measurements (Chappard et al., 2010; Pulkkinen et al., 2008), such as cortical thickness and NSA. Thus, creating a model based on these measurements and separating trabecular and cortical compartments was justified. Our method specifically concerns the femoral neck area, since the morphology and structure of the femoral neck play a critical role in the risk of cervical hip fracture (Duboeuf et al., 1997). Furthermore, it has been shown that low failure load fractures predominantly occur at the femoral neck (Pulkkinen et al., 2006) and that these fractures are under-diagnosed using the current clinical procedure (Pulkkinen et al., 2010). Since the directional organisation of trabecular bone may offer further improvement in fracture risk evaluation (Chappard et al., 2005), we decided to distribute the
Fig. 5. Typical distribution of reconstruction error compared to manually segmented model from CT-scans.
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
443
In conclusion, we present here a novel method to obtain a 3D finite element model based on a standard 2D radiographic picture. The heterogeneous distribution of material properties for the model is extracted from the trabecular structure analysis of the radiograph. The method appeared to construct the patient-specific 3D shape with feasible accuracy when compared to the true 3D model constructed from CT scans. In addition, the proximal femoral failure load in a sideways fall configuration can be estimated with reasonable accuracy by using the 3D FE model based on a standard 2D radiographic picture. This model appeared to be more predictive for fracture load than DXA. The accuracy of the model and the method used provide promise for effective clinical use. Conflict of interest statement The authors have no conflicts of interest related to this study. Acknowledgements
Fig. 6. Experimental failure load vs. failure load predicted from a radiograph-based 3D finite element model. N ¼21.
material properties with an emphasis on the localisation and orientation of the trabecular principal tensile and principal compressive systems. Based on the studies of Pulkkinen et al. (2008) and Chappard et al. (2010), texture parameters derived from picture analysis of plain radiographs were found to correlate well with BMD and can be used to evaluate fracture risk. In the present study, texture analysis was similar to that previously reported (Thevenot et al., 2013), with a GLCM orientated along the fibres instead of perpendicular to them. Furthermore, the homogeneity index was not normalised by the thickness of the bone, as it was suggested that it gives better correlation with the Hounsfield Unit assessed volumetrically (Thevenot et al., 2013) and, eventually, a better estimation of the material properties of the model. The function of the principal tensile trabecular system is to transfer the forces from the femoral head to the shaft area during activities by using an arch-shape, whereas the principal compressive system transmits forces to the medial aspect of the shaft (Rudman et al., 2006). Based on Wolff's law (Wolff, 1892), the bone structure changes according to the stress acting on it (Miller et al., 2005), and it has been suggested that the principal orientation used in an orthotropic model should be defined from the trabecular structure (Yang et al., 2010). Some limitations remain in our methodology. First, this study is based on cadaver material. For in vivo applications, the radiographic structure analysis would require extra filtering to remove the artefacts caused by the soft tissues and pelvis. In vivo radiography may also include inaccuracies due to the internal rotation of the hip during imaging. However, the present method offers the possibility to adjust the rotation angle, by which the bias can be corrected. Furthermore, the number of bones used for the training was quite restricted, and even if they represented different geometries, the accuracy of the generated shape was limited in some areas, especially in the top of the trochanteric area and near the minor trochanter. Nevertheless, the lack of accuracy in these areas is not a critical problem, since clinical hip fractures typically occur in the femoral neck and intertrochanteric regions. Finally, the algorithm has shown difficulties in generating bones with peculiar geometry; therefore, this problem will be solved in the future by increasing the robustness of the script.
Mr Ilja Everilä and Mr Mikko Paakkolanvaara are acknowledged for their help in image processing, and Dr Thomas Link and Dr Pasi Pulkkinen for help in the collection of original data. This study was financially supported by the Finnish Funding Agency for Technology and Innovation, the Academy of Finland, the International Graduate School in Biomedical Engineering and Medical Physics, the National Doctoral Programme of Musculoskeletal Disorders and Biomaterials, the Finnish Cultural Foundation, the North Ostrobothnia Regional Fund, and Deutsche Forschungsgemeinschaft Grant DFG LO 730/2-2. Study sponsors had no involvement in the study design, in the collection, analysis and interpretation of data; in the writing of the manuscript; and in the decision to submit the manuscript for publication. References Augat, P., Reeb, H., Claes, L.E., 1996. Prediction of fracture load at different skeletal sites by geometric properties of the cortical shell. J. Bone Miner. Res. 11, 1356–1363. Baker-LePain, J., Luker, K., Lynch, J., Parimi, N., Nevitt, M., Lane, N., 2011. Active shape modeling of the hip in the prediction of incident hip fracture. J. Bone Miner. Res. 26, 468–474. Bayraktar, H.H., Morgan, E.F., Niebur, G.L., Morris, G.E., Wong, E.K., Keaveny, T.M., 2004. Comparison of the elastic and yield properties of human femoral trabecular and cortical bone tissue. J. Biomech. 37, 27–35. Beck, T.J., Ruff, C.B., Warden, K.E., Scott Jr., W.W., Rao, G.U., 1990. Predicting femoral neck strength from bone mineral data. A structural approach. Invest. Radiol. 25 (1), 6–18. Bessho, M., Ohnishi, I., Matsuyama, J., Matsumoto, T., Imai, K., Nakamura, K., 2007. Prediction of strength and strain of the proximal femur by a CT-based finite element method. J. Biomech. 40, 1745–1753. Bessho, M., Ohnishi, I., Okazaki, H., Sato, W., Kominami, H., Matsunaga, S., Nakamura, K., 2004. Prediction of the strength and fracture location of the femoral neck by CT-based finite-element method: a preliminary study on patients with hip fracture. J. Orthop. Sci. 9, 545–550. Bryan, R., Nair, P.B., Taylor, M., 2009. Use of a statistical model of the whole femur in a large scale, multi-model study of femoral neck fracture risk. J. Biomech. 42 (13), 2171–2176. Chappard, C., Bousson, V., Bergot, C., Mitton, D., Marchadier, A., Moser, T., Benhamou, C.L., Laredo, J.D., 2010. Prediction of femoral fracture load: crosssectional study of texture analysis and geometric measurements on plain radiographs versus bone mineral density. Radiology 255 (2), 536–543. Chappard, C., Brunet-Imbault, B., Lemineur, G., Giraudeau, B., Basillais, A., Harba, R., Benhamou, C.L., 2005. Anisotropy changes in post-menopausal osteoporosis: characterisation by a new index applied to trabecular bone radiographic images. Osteoporos. Int. 16, 1193–1202. Cody, D.D., Gross, G.J., Hou, F.J., Spencer, H.J., Goldstein, S.A., Fyhrie, D.P., 1999. Femoral strength is better predicted by finite element models than QCT and DXA. J. Biomech. 32, 1013–1020. Cummings, S.R., Melton, L.J., 2002. Epidemiology and outcomes of osteoporotic fractures. Lancet 359, 1761–1767. Drucker, D.C., Prager, W., 1952. Soil mechanics and plastic analysis for limit design. Q. Appl. Math. 10, 157–165.
444
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
Duboeuf, F., Hans, D., Schott, A.M., Kotzki, P.O., Favier, F., Marcelli, C., Meunier, P.J., Delmas, P.D., 1997. Different morphometric and densitometric parameters predict cervical and trochanteric hip fracture: the EPIDOS Study. J. Bone Miner. Res. 12 (11), 1895–1902. Duchemin, L., Bousson, V., Raossanaly, C., Bergot, C., Laredo, J.D., Skalli, W., Mitton, D., 2008. Prediction of mechanical properties of cortical bone by quantitative computed tomography. Med. Eng. Phys. 30, 321–328. Eckstein, F., Wunderer, C., Boehm, H., et al., 2004. Reproducibility and side differences of mechanical tests for determining the structural strength of the proximal femur. J. Bone Miner. Res. 19, 379–385. Gnudi, S., Ripamonti, C., Lisi, L., Fini, M., Giardino, R., Giavaresi, G., 2002. Proximal femur geometry to detect and distinguish femoral neck. fractures from trochanteric fractures in postmenopausal women. Osteoporos. Int. 13, 69–73. Gregory, J.S., Testi, D., Stewart, A., Undrill, P.E., Reid, D.M., Aspden, R.M., 2004. A method for assessment of the shape of the proximal femur and its relationship to osteoporotic hip fracture. Osteoporos. Int. 15, 5–11. Gunay, M., Shim, M-B., Shimada, K., 2007. Cost- and time-effective three-dimensional bone-shape reconstruction from X-ray images. Int. J. Med. Rob. Comput. Assisted Surg. 3, 323–335. Kanis, J.A., 2002. Diagnosis of osteoporosis and assessment of fracture risk. Lancet 359, 1929–1936. Keller, T.S., 1994. Predicting the compressive mechanical behaviour of bone. J. Biomech. 27, 1159–1168. Keyak, J.H., Lee, I.Y., Skinner, H.B., 1994. Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometric measures. J. Biomed. Mater. Res. 28, 1329–1336. Keyak, J.H., Rossi, S.A., 2000. Prediction of femoral fracture load using finite element models: an examination of stress- and strain-based failure theories. J. Biomech. 33, 209–214. Keyak, J.H., Skinner, H.B., 1992. Three-dimensional finite element modelling of bone: effects of element size. J. Biomed. Eng. 14, 483–489. Koivumäki, J.E.M., Thevenot, J., Pulkkinen, P., Kuhn, V., Link, T.M., Eckstein, F., Jämsä, T., 2012. CT-based finite element models can be used to estimate experimentally measured failure loads in the proximal femur. Bone 50, 824–829. Langton, C.M., Pisharody, S., Keyak, J.H., 2009. Generation of a 3D proximal femur shape from a single projection 2D radiographic image. Osteoporos. Int. 20 (3), 455–461. Lengsfeld, M., Schmitt, J., Alter, P., Kaminsky, J., Leppek, R., 1998. Comparison of geometry-based and CT voxel-based finite element modelling and experimental validation. Med. Eng. Phys. 20, 515–522. Meunier, P.J., 1993. Prevention of hip fractures. Am. J. Med. 95 (5A), 75S–78S. Michelotti, J., Clark, J., 1999. Femoral neck length and hip fracture risk. J. Bone Miner. Res. 14, 1714–1720. Miller, Z., Fuchs, M.B., Arcan, M., 2005. Trabecular bone adaptation with an orthotropic material model. J. Biomech. 35, 247–256. Partanen, J., Jämsä, T., Jalovaara, P., 2001. Influence of the upper femur and pelvic geometry on the risk and type of hip fractures. J. Bone Miner. Res. 16, 1540–1546. Pulkkinen, P., Eckstein, F., Lochmüller, E-M., Kuhn, V., Jämsä, T., 2006. Association of geometric factors and failure load level with the distribution of cervical vs. trochanteric hip fractures. J. Bone Miner. Res. 21 (6), 895–901. Pulkkinen, P., Jämsä, T., Lochmüller, E.M., Kuhn, V., Nieminen, M.T., Eckstein, F., 2008. Experimental hip fracture load can be predicted from plain radiography
by combined analysis of trabecular bone structure and bone geometry. Osteoporos. Int. 19, 547–558. Pulkkinen, P., Luomala, M., Obregon, V., Jämsä, T., 2005. Gradient-based method for the assessment of trabecular structure from radiographs. IFMBE Proc., 12. Pulkkinen, P., Partanen, J., Jalovaara, P., Jämsä, T., 2004. Combination of bone mineral density and upper femur geometry improve the prediction of hip fracture. Osteoporos. Int. 15, 274–280. Pulkkinen, P., Partanen, J., Jalovaara, P., Jämsä, T., 2010. BMD T-score discriminates trochanteric fractures from unfractured controls, whereas geometry discriminates cervical fracture cases from unfractured controls of similar BMD. Osteoporos. Int. 21 (7), 1269–1276. Rudman, K.E., Aspden, R.M., Meakin, J.R., 2006. Compression or tension? The stress distribution in the proximal femur. Biomed. Eng. Online 5, 12. Schileo, E., Taddei, F., Cristofolini, L., Viceconti, M., 2008. Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and fracture location on human femurs tested in vitro. J. Biomech. 41, 356–367. Schuit, S., van der Klift, M., Weel, A., Laet, Cde, Burger, H., Seeman, E., Hofman, A., Uitterlinden, A., van Leeuwen, J., Pols, H., 2004. Fracture incidence and association with bone mineral density in elderly men and women: the Rotterdam Study. Bone 34, 195–202. Schumann, S., Tannast, M., Nolte, L.P., Zheng, G., 2010. Validation of statistical shape model based reconstruction of the proximal Afemur—a morphology study. Med. Eng. Phys. 32 (6), 638–644. Szulc, P., Duboeuf, F., Schott, A.M., Dargent-Molina, P., Meunier, P.J., Delmas, P.D., 2006. Structural determinants of hip fracture in elderly women: re-analysis of the data from the EPIDOS study. Osteoporos. Int. 17, 231–236. Testi, D., Cappello, A., Sgallari, F., Rumpf, M., Viceconti, M., 2004. A new software for prediction of femoral neck fractures. Comput. Methods Programs Biomed. 75 (2), 141–145. Testi, D., Viceconti, M., Cappello, A., Gnudi, S., 2002. Prediction of hip fracture can be significantly improved by a single biomedical indicator. Ann. Biomed. Eng. 30, 801–807. Thevenot, J., Hirvasniemi, J., Finnilä, M., Pulkkinen, P., Kuhn, V., Link, T., Eckstein, F., Jämsä, T., Saarakkala, S., 2013. Trabecular homogeneity index derived from plain radiograph to evaluate bone quality. J. Bone Miner. Res. , http://dx.doi.org/ 10.1002/jbmr.1987. ([Epub ahead of print]). Thevenot, J., Pulkkinen, P., Koivumäki, J.E.M., Kuhn, V., Eckstein, F., Jämsä, T., 2009. Discrimination of cervical and trochanteric hip fractures using radiographybased two-dimensional finite element models. Open Bone J. 1, 16–22. Viceconti, M., 1998. A comparative study on different methods of automatic mesh generation of human femurs. Med. Eng. Phys. 20, 1–10. Wolff, J., 1892. Das Gesetz der transformation der Knochen, Translation by Maquet P & Furlong R: The Law of Bone Remodelling. Springer-Verlag, Berlin p. 1986. Yang, H., Ma, X., Guo, T., 2010. Some factors that affect the comparison between isotropic and orthotropic inhomogeneous finite element material models of femur. Med. Eng. Phys. 32, 553–560. Zebaze, R., Jones, A., Welsh, F., Knackstedt, M., Seeman, E., 2005. Femoral neck shape and the spatial distribution of its mineral mass varies with its size: Clinical and biomechanical implications. Bone 37, 243–252. Zheng, G., Gollmer, S., Schumann, S., Dong, X., Feilkas, T., González Ballester, M.A., 2009. A 2D/3D correspondence building method for reconstruction of a patient-specific 3D bone surface model using point distribution models and calibrated X-ray images. Med. Image Anal. 13 (6), 883–899.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
A novel methodology for generating 3D finite element models of the hip from 2D radiographs Jérôme Thevenot a, Janne Koivumäki a, Volker Kuhn b,c, Felix Eckstein b,d, Timo Jämsä a,e,n a
Department of Medical Technology, University of Oulu, Oulu 90014, Finland Institute of Anatomy, Ludwig-Maximilians-University, Munich, Germany c Department of Trauma Surgery and Sports Medicine, Medical University, Innsbruck, Austria d Institute of Anatomy and Musculoskeletal Research, Paracelsus Medical University, Salzburg, Austria e Department of Diagnostic Radiology, Oulu University Hospital, Oulu, Finland b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 6 November 2013
Finite element (FE) modelling has been proposed as a tool for estimating fracture risk and patientspecific FE models are commonly based on computed tomography (CT). Here, we present a novel method to automatically create personalised 3D models from standard 2D hip radiographs. A set of geometrical parameters of the femur were determined from seven a–p hip radiographs and compared to the 3D femoral shape obtained from CT as training material; the error in reconstructing the 3D model from the 2D radiographs was assessed. Using the geometry parameters as the input, the 3D shape of another 21 femora was built and meshed, separating a cortical and trabecular compartment. The material properties were derived from the homogeneity index assessed by texture analysis of the radiographs, with focus on the principal tensile and compressive trabecular systems. The ability of these FE models to predict failure load as determined by experimental biomechanical testing was evaluated and compared to the predictive ability of DXA. The average reconstruction error of the 3D models was 1.77 mm (7 1.17 mm), with the error being smallest in the femoral head and neck, and greatest in the trochanter. The correlation of the FE predicted failure load with the experimental failure load was r2 ¼64% for the reconstruction FE model, which was significantly better (po 0.05) than that for DXA (r2 ¼ 24%). This novel method for automatically constructing a patient-specific 3D finite element model from standard 2D radiographs shows encouraging results in estimating patient-specific failure loads. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Hip fracture Biomechanics Finite element method Structure analysis Radiography
1. Introduction The clinical diagnosis of osteoporosis and eventual fracture risk is based on bone mineral density (BMD) assessed from dual energy X-ray absorptiometry (DXA). However, DXA-based BMD has been shown to be insufficient for predicting the individual risk of fracture (Kanis, 2002; Schuit et al., 2004); this can be explained by the fact that the mechanical strength of the hip is highly dependent on the shape and structure of the bone (Beck et al., 1990; Gnudi et al., 2002; Partanen et al., 2001; Pulkkinen et al., 2006), and fractures also occur in individuals with non-osteoporotic BMD (Cummings and Melton, 2002; Meunier, 1993; Schuit et al., 2004). The finite element (FE) method has shown promising results for the clinical assessment of hip strength (Cody et al., 1999). While some studies were concentrated on 2D FE models based on
n Corresponding author at: Department of Medical Technology, University of Oulu, Oulu 90014, Finland, Tel.: þ 358 29 448 6001. E-mail address: timo.jamsa@oulu.fi (T. Jämsä).
0021-9290/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jbiomech.2013.11.004
DXA images (Testi et al., 2002, 2004), most of them used 3D models derived from computed tomography (CT) scans to estimate hip strength (Bessho et al., 2004; Keyak and Rossi, 2000; Schileo et al., 2008). Unfortunately, building patient-specific 3D FE models requires the use of high-cost volumetric imaging coupled with high radiation doses and computational power. A solution to the time-consuming generation of 3D models from CT scans is the use of active shape modelling (ASM) (Gunay et al., 2007) or statistical shape modelling (Schumann et al., 2010), both of which have proven to be sufficiently accurate for different applications. Briefly, these methods will iterate an initial 3D object to fit a 2D input contour or geometrical parameters (Baker-LePain et al., 2011; Bryan et al., 2009; Gregory et al., 2004; Langton et al., 2009; Schumann et al., 2010; Zheng et al., 2009). However, these methods usually use contours extracted from 2D images and do not separate cortical and trabecular bone even though cortical thickness has shown to be an important predictor of hip fracture risk (Augat et al., 1996; Chappard et al., 2010; Gnudi et al., 2002; Michelotti and Clark, 1999; Pulkkinen et al., 2004, 2008; Szulc et al., 2006).
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
The aim of this study is to provide a new automatic method to generate 3D finite element models from standard 2D radiographic pictures. The method will create both cortical and trabecular bone compartments from a set of geometrical parameters as an input. This paper describes the methodology for generating the models and presents the results in estimating experimental failure load.
2. Study sample and imaging For the present study, we used data of twenty-eight cadaver femora obtained from a larger experimental study at the Institute of Anatomy at the Ludwig Maximilians University of Munich (Germany) (Eckstein et al., 2004; Pulkkinen et al., 2006). Seven of the bones were used for training, and twenty-one for validation. The preparation and storage of the samples has previously been described in detail (Eckstein et al., 2004). After the dissection, excision and cleaning of the surrounding soft tissue, the proximal femora were radiographed and CT-scanned. Radiographs were taken using a Faxitron X-ray system (Model 43885A; Faxitron, Hewlett Packard, McMinnville, OR) and the X-ray films were digitised together with a calibration scale using a scanner at 600 dpi. CT scans were performed using a 16-row Multi-detector CT (MD-CT) scanner (Sensation 16; Siemens Medical Solutions, Erlangen, Germany). MD-CT and DXA scan setups have been previously presented (Koivumäki et al., 2012). For CT scans, a high-resolution protocol with a slice thickness of 0.75 mm was used. The settings were 120 kVp and 100 mAs, a 512 512 pixels image matrix, and a field of view of 100 mm. The in-plane spatial resolution was 0.25 mm 0.25 mm. A calibration phantom (Osteo Phantom, Siemens, Erlangen, Germany) composed of hydroxyapatite (0 mg/cm3 for water and 200 mg/ cm3 for bone) was placed below the specimens. In vitro DXA scans of the femora were obtained with a standard narrow fan beam scanner (GE LunarProdigy, GE Lunar, Madison, WI) with the femoral specimens submerged in a water bath. Standard positioning was used across all specimens, and the proximal femoral BMD values were evaluated with the software provided by the manufacturer. The cortical and trabecular bone compartments were segmented from the CT scans as previously presented (Koivumäki et al., 2012) using Mimics (v12.1, Materialise, Leuven, Belgium), with the minimum thickness for the cortical bone adjusted to 1 mm. A maximum element size of 3 mm was applied. The femurs were mechanically tested simulating a fall on the greater trochanter, as previously presented (Eckstein et al., 2004; Pulkkinen et al., 2006). Briefly, femora were positioned at 101 from horizontal plane and the neck at 151 internal rotation; loads were then applied to the greater trochanter through a pad at a rate of 6.6 mm/s. The failure load was determined from the loaddeformation curve.
439
3.1. Generation of the 3D shape 3.1.1. Geometrical parameters as an input A limited set of geometrical parameters was used to generate the volumetric shape of the model. The measurements of the geometrical parameters were done using Image Tool software (version 3.00; University of Texas Health Science Center, San Antonio, TX), as previously described (Partanen et al., 2001; Pulkkinen et al., 2006). Femoral neck axis length (FNAL), neck-shaft angle (NSA), trochanteric width (TW), femoral head diameter (HD), femoral neck diameter (ND), femoral shaft (FSD) diameter, femoral shaft cortex (FSC), and calcar femoral cortex width (CFC) were measured, as shown in Fig. 1.
3.1.2. Building up the 3D shape In order to establish relationships between the geometrical parameters and the skeleton of the femurs, seven bones with a large range of geometrical parameters were selected as CT-based 3D training material. The impact of each geometrical parameter on the volumetric shape of the bone was evaluated to delimitate the outer shape of the bone by points. Curves were generated from the interpolation of these points to represent the skeleton of the bone and surfaces were eventually derived from them (Fig. 2). The femoral head was considered a half-sphere (Schumann et al., 2010; Thevenot et al., 2009). Specific attention was given to the shape of the femoral neck since this reflects the risk of cervical fractures (Pulkkinen et al., 2006). The non-circular cross-section of the femoral neck (Zebaze et al., 2005) was estimated by using the
3. Construction of volumetric finite element model The generation of the model can be divided into three parts: (1) generation of the 3D shape, (2) meshing, and (3) determination and implementation of the material properties. The algorithm was developed under Matlab (version 7.1, MathWorks, Natick, MA). The output of the algorithm includes three command script files corresponding to each part of the model generation. The files were then imported into the pre- and post-processing software Femap (version 9.2, UGS Corp., Plano, TX) and launched consecutively in order to generate the model.
Fig. 1. Geometrical parameters. (A) Neck-shaft angle, (B‐C) and (E‐F) femoral shaft and calcar femoral cortical thickness, respectively, (B‐D) femoral shaft diameter, (E‐ G) trochanteric width, (I‐H) femoral neck axis length, (J‐K) and (L‐M) head and neck diameters.
440
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
Fig. 2. Building up the shape of the upper femur: (A) femoral head and main axis, (B) skeleton of the neck area, (C) neck surfaces based on previous curves and (D) skeleton of the trochanteric area.
Fig. 3. Left: generated model with the principle tensile system (PTS) and principal compressive system (PCS) shown in red. Right: Typical distribution of material properties, dark values corresponding to low density and light to high density.
training material. The neck area was divided into 10 segments along the femoral neck axis, and for each of these segments, the radius of the neck was checked from CT data in every 451 rotational direction. A relationship was established between each radius and FNAL, HD, ND and FSD, to obtain accurate shape estimation for an individual femoral neck. 3.2. Meshing A smoothed solid was created from the surfaces of the model and then meshed with tetrahedral elements with a size increasing progressively from the surface to the centre of the bone to a maximum element size of 3 mm. The size of the elements near the surface of the bone was chosen to be proportional to the FSC, being less than half its size. The faces of elements on the surface were extruded with a thickness of 1 mm, corresponding to the selected minimum thickness of the cortical bone and representing the outer contour of the bone. Cortical and trabecular compartments were separated during the meshing process. The inner contour of the cortical bone was delimited from the elements within the cortex as defined by FSC and CFC. FSC was represented by two consecutive layers of elements from the outer surface in the shaft area, whereas CFC was defined by elements within the lower neck area. Finally, the cortical bone was re-meshed using 10-noded tetrahedral elements.
The remaining trabecular bone was re-meshed using 4-noded tetrahedral elements. The size of the elements along the model varies between 1 and 3 mm, with the smaller elements being close to the inner and outer surfaces of the cortical bone. The maximum size of the elements was chosen as previously suggested (Keyak and Skinner, 1992). 3.3. Material properties Distribution of the material properties was evaluated from trabecular structure analysis of the radiographic picture. The spatial distribution of the material properties in the trabecular bone was related to the principal tensile and compressive systems (PTS and PCS, respectively) (Rudman et al., 2006), representing the highest mineralised areas of the trabecular bone. Here, the Young's moduli were assessed using the homogeneity index (HI), adjusted to the material properties derived from Hounsfield Units (HU) of CT scans, with the relationship being based on the seven bones used as training material. 3.3.1. Localisation of the principal tensile and compressive trabecular systems The general localisation and thickness of the PTS was assessed in 2D from the radiographs and its 3D orientation from the
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
441
CT-scans, using the training material as reference. A main curve representing the PTS from the greater trochanter to the femoral head was generated in the model. Similarly, the PCS was approximated to go from the lower neck area to the top of the femoral head and was represented as a curved frustum (Fig. 3, left).
homogeneity index (HI). The first column and first row of the GLCM were cropped to remove the values corresponding to empty spaces in the bone, and to consequently evaluate the homogeneity of the mineralised trabecular bone. HI, as a measure of the closeness of the distribution of elements near the diagonal of the GLCM, was calculated as follows:
3.3.2. Structural analysis of the radiographs The image processing procedure presented in Fig. 4 was performed as previously reported (Thevenot et al., 2013). It was performed with an algorithm developed under Matlab with a graphical interface allowing the user to select a specific region of interest (ROI) (15 mm 15 mm). In image pre-processing, noise was removed by using a 3 3 median filter, followed by morphological top and bottom hat operations. A gradient-based grey-level image was then constructed, which enabled the visualisation of the trabecular structure even under varying image contrasts and brightness conditions of the radiograph (Pulkkinen et al., 2005). The ROI was chosen within the fibres of the PCS (inferomedial neck, at the neck/head transition area) and the trabecular main orientation (TMO) was then evaluated using Fourier transform. An empirically-derived threshold of 0.03 was used to remove nonsignificant frequencies in the Fourier-based image, and linear regression was applied to fit a line g(x)¼ axþb to the threshold values in order to extract the TMO in the frequency domain. Eventually, the angle α between the TMO and the shaft axis in the space domain was calculated and the ROI rotated along it. A greylevel co-occurrence matrix (GLCM) with a distance of one pixel was calculated along the TMO for determination of the
GLCMði; jÞ HI ¼ ∑ i;j 1 þ i j
Selection of the region of interest (ROI)
Filtering of the picture
Gradient-based transform
Calculation of trabecular main orientation (TMO) using Fourier transform (angle α)
Rotation α of the ROI: alignment along the TMO
Gray-level co-occurrence matrix (GLCM) along the TMO
The evaluated radiograph-based HI of the trabeculae-rich lower femoral neck area was used to derive the Young's modulus along the FE model. 3.3.3. Implementation of the material properties Young's modulus of 14.2 GPa was applied for the cortical bone (Viceconti, 1998). An intermediate Young's modulus of 7 GPa was used to the first layer of trabecular elements to simulate the transition between cortical and trabecular bones. A Poisson's ratio of 0.33 was selected for both cortical and trabecular bone (Lengsfeld et al., 1998). CT scan data of the training material were used to establish the relationship between the calculated HI and the material properties of the model. An average grey value of all of the voxels inside an element was calculated in Mimics (Koivumäki et al., 2012). The bone equivalent density (ash density, ρash) was then defined by assuming a linear relationship in which the density is proportional to the attenuation (Koivumäki et al., 2012). The Young's modulus (E, MPa) was estimated using the equation E ¼10.095ρ, where ρ is the CT density value (mg/cm3) (Duchemin et al., 2008). Yield stress was calculated using the previously proposed equations, 114ρash1.72 and 137ρash1.88, for cortical and trabecular bone, respectively (Keller, 1994; Keyak et al., 1994). The curve representing the trabecular PTS orientation was divided into 20 sub-curves and elements surrounding them within a consistent estimated PTS thickness were affected with a specific Young's modulus. An extra transition layer of elements with reduced Young's modulus was generated. According to the CT data, the Young's modulus increases from the greater trochanter to the femoral head area along the tensile system. Here, this increase was estimated by linear increase of the modulus through the PTS, with each sub-area representing different material properties. A similar method was applied for the trabecular PCS. Here, the higher Young's modulus values were located in the lower neck and the lower values in the upper part of the femoral head, differences between which were estimated using a linear distribution. The material properties of trabecular elements outside PTS and PCS were grouped depending on their location. Homogeneous material properties were affected separately for shaft, trochanteric, femoral neck, and femoral head areas, derived from the HI. As a result, a total of 66 material properties were created: - Cortical bone (1) and the transition layer between cortical and trabecular bone (1). - PTS (20) and the transition layer (20). - PCS (10) and the transition layer (10). - Trabecular areas outside PTS and PCS (shaft, trochanteric, femoral neck, and femoral head) (4).
GLCM reduction 3.4. Finite element analysis of failure load
Calculation of the homogeneity index based on the GLCM Fig. 4. Flowchart of the image processing to obtain the homogeneity index of the trabecular structure.
After applying the boundary conditions to simulate the experiment (Eckstein et al., 2004), bilinear elastoplastic FE analysis was performed, as previously presented (Koivumäki et al., 2012), using the Newton–Raphson method and Drucker–Prager yield criterion (Drucker and Prager, 1952), applying a post-yield modulus 5% of
442
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
the Young's modulus (Bayraktar et al., 2004), and an ultimate tensile stress presumed to be 0.8 times the compressive yield stress (Bessho et al., 2007). The fracture was determined to occur when a surface cortical element failed (Bessho et al., 2007), failure in tension was characterised for a maximum principal stress higher than the ultimate tensile stress and failure in compression for a minimum principal strain lower than 7300 microstrain (Koivumäki et al., 2012). The failure load predicted from the finite element analysis was compared with the experimentally measured failure load. 3.5. Statistical analysis
5. Discussion
Linear regression analysis was used to evaluate the relationship between the simulated and the experimentally measured fracture loads, and between BMD and BMC and the experimental fracture load. One-tailed Fisher z-transformation was used to test the hypothesis that the simulated fracture load had a higher correlation with the experimental fracture load than DXA. Statistical analyses were performed using SPSS software (SPSS 16.0 for Windows, SPSS Inc. Chicago, IL). A p-value o0.05 was considered to be statistically significant.
4. Results The present methodology generated 3D FE models of upper femur from 2D radiographs with visually similar appearances to the CT-based models. The average reconstruction error for the whole upper femur was 1.77 mm ( 71.17 mm), with the error being lowest in the femoral head (1.03 mm 70.45 mm) and femoral neck (1.27 mm 70.60 mm) and greatest in the trochanteric area (2.24 mm 7 1.33 mm) (Table 1, Fig. 5). The relationship between the estimated and the experimental failure load is presented in Fig. 6. The estimated failure load values Table 1 Average reconstruction error for different upper femur areas of randomly selected validation models (N ¼ 7). Data presented as mean 7 SD. Femoral head Femoral neck Trochanter Shaft area Error (mm) 1.03 7 0.45
were well correlated with the experimental ones (r2 ¼64%), SEE¼ 543 N, po 0.001). The slope was 1.20 (not significantly different from 1, p ¼0.353), with an intercept of 472 N (not significantly different from 0, p¼ 0.477). The coefficients of determination between BMD and BMC, and the experimental fracture load were r2 ¼24% (N ¼ 20, SEE¼ 793 N, p¼ 0.028) and r2 ¼35% (N ¼20, SEE ¼734 N, p ¼0.006), respectively. The coefficient of determination was significantly higher for the FE model when compared to BMD (p o0.05) but not when compared to BMC (p ¼0.11).
1.277 0.60
Total
2.247 1.33 1.577 0.79 1.777 1.17
In the present study, we developed a novel method to obtain a 3D finite element model based on a standard 2D radiographic picture. The 3D model is built up automatically using the geometric shape of the femur derived from simple measurements of geometrical dimensions. The heterogeneous distribution of material properties in the model is based on the trabecular structure analysis of the radiograph with a focus on the orientation and location of the principal tensile and compressive systems. This method appeared able to construct the patient-specific 3D shape with feasible accuracy, and the FE analysis showed promising results in estimating the individual experimental failure load. While most studies building 3D models from 2D data use the ASM method from the outer contour of the bone, we decided to create a method that also defines the inner contour of the cortical bone. It has been shown that the failure load can be predicted with reasonable accuracy from radiographic texture analysis and geometric measurements (Chappard et al., 2010; Pulkkinen et al., 2008), such as cortical thickness and NSA. Thus, creating a model based on these measurements and separating trabecular and cortical compartments was justified. Our method specifically concerns the femoral neck area, since the morphology and structure of the femoral neck play a critical role in the risk of cervical hip fracture (Duboeuf et al., 1997). Furthermore, it has been shown that low failure load fractures predominantly occur at the femoral neck (Pulkkinen et al., 2006) and that these fractures are under-diagnosed using the current clinical procedure (Pulkkinen et al., 2010). Since the directional organisation of trabecular bone may offer further improvement in fracture risk evaluation (Chappard et al., 2005), we decided to distribute the
Fig. 5. Typical distribution of reconstruction error compared to manually segmented model from CT-scans.
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
443
In conclusion, we present here a novel method to obtain a 3D finite element model based on a standard 2D radiographic picture. The heterogeneous distribution of material properties for the model is extracted from the trabecular structure analysis of the radiograph. The method appeared to construct the patient-specific 3D shape with feasible accuracy when compared to the true 3D model constructed from CT scans. In addition, the proximal femoral failure load in a sideways fall configuration can be estimated with reasonable accuracy by using the 3D FE model based on a standard 2D radiographic picture. This model appeared to be more predictive for fracture load than DXA. The accuracy of the model and the method used provide promise for effective clinical use. Conflict of interest statement The authors have no conflicts of interest related to this study. Acknowledgements
Fig. 6. Experimental failure load vs. failure load predicted from a radiograph-based 3D finite element model. N ¼21.
material properties with an emphasis on the localisation and orientation of the trabecular principal tensile and principal compressive systems. Based on the studies of Pulkkinen et al. (2008) and Chappard et al. (2010), texture parameters derived from picture analysis of plain radiographs were found to correlate well with BMD and can be used to evaluate fracture risk. In the present study, texture analysis was similar to that previously reported (Thevenot et al., 2013), with a GLCM orientated along the fibres instead of perpendicular to them. Furthermore, the homogeneity index was not normalised by the thickness of the bone, as it was suggested that it gives better correlation with the Hounsfield Unit assessed volumetrically (Thevenot et al., 2013) and, eventually, a better estimation of the material properties of the model. The function of the principal tensile trabecular system is to transfer the forces from the femoral head to the shaft area during activities by using an arch-shape, whereas the principal compressive system transmits forces to the medial aspect of the shaft (Rudman et al., 2006). Based on Wolff's law (Wolff, 1892), the bone structure changes according to the stress acting on it (Miller et al., 2005), and it has been suggested that the principal orientation used in an orthotropic model should be defined from the trabecular structure (Yang et al., 2010). Some limitations remain in our methodology. First, this study is based on cadaver material. For in vivo applications, the radiographic structure analysis would require extra filtering to remove the artefacts caused by the soft tissues and pelvis. In vivo radiography may also include inaccuracies due to the internal rotation of the hip during imaging. However, the present method offers the possibility to adjust the rotation angle, by which the bias can be corrected. Furthermore, the number of bones used for the training was quite restricted, and even if they represented different geometries, the accuracy of the generated shape was limited in some areas, especially in the top of the trochanteric area and near the minor trochanter. Nevertheless, the lack of accuracy in these areas is not a critical problem, since clinical hip fractures typically occur in the femoral neck and intertrochanteric regions. Finally, the algorithm has shown difficulties in generating bones with peculiar geometry; therefore, this problem will be solved in the future by increasing the robustness of the script.
Mr Ilja Everilä and Mr Mikko Paakkolanvaara are acknowledged for their help in image processing, and Dr Thomas Link and Dr Pasi Pulkkinen for help in the collection of original data. This study was financially supported by the Finnish Funding Agency for Technology and Innovation, the Academy of Finland, the International Graduate School in Biomedical Engineering and Medical Physics, the National Doctoral Programme of Musculoskeletal Disorders and Biomaterials, the Finnish Cultural Foundation, the North Ostrobothnia Regional Fund, and Deutsche Forschungsgemeinschaft Grant DFG LO 730/2-2. Study sponsors had no involvement in the study design, in the collection, analysis and interpretation of data; in the writing of the manuscript; and in the decision to submit the manuscript for publication. References Augat, P., Reeb, H., Claes, L.E., 1996. Prediction of fracture load at different skeletal sites by geometric properties of the cortical shell. J. Bone Miner. Res. 11, 1356–1363. Baker-LePain, J., Luker, K., Lynch, J., Parimi, N., Nevitt, M., Lane, N., 2011. Active shape modeling of the hip in the prediction of incident hip fracture. J. Bone Miner. Res. 26, 468–474. Bayraktar, H.H., Morgan, E.F., Niebur, G.L., Morris, G.E., Wong, E.K., Keaveny, T.M., 2004. Comparison of the elastic and yield properties of human femoral trabecular and cortical bone tissue. J. Biomech. 37, 27–35. Beck, T.J., Ruff, C.B., Warden, K.E., Scott Jr., W.W., Rao, G.U., 1990. Predicting femoral neck strength from bone mineral data. A structural approach. Invest. Radiol. 25 (1), 6–18. Bessho, M., Ohnishi, I., Matsuyama, J., Matsumoto, T., Imai, K., Nakamura, K., 2007. Prediction of strength and strain of the proximal femur by a CT-based finite element method. J. Biomech. 40, 1745–1753. Bessho, M., Ohnishi, I., Okazaki, H., Sato, W., Kominami, H., Matsunaga, S., Nakamura, K., 2004. Prediction of the strength and fracture location of the femoral neck by CT-based finite-element method: a preliminary study on patients with hip fracture. J. Orthop. Sci. 9, 545–550. Bryan, R., Nair, P.B., Taylor, M., 2009. Use of a statistical model of the whole femur in a large scale, multi-model study of femoral neck fracture risk. J. Biomech. 42 (13), 2171–2176. Chappard, C., Bousson, V., Bergot, C., Mitton, D., Marchadier, A., Moser, T., Benhamou, C.L., Laredo, J.D., 2010. Prediction of femoral fracture load: crosssectional study of texture analysis and geometric measurements on plain radiographs versus bone mineral density. Radiology 255 (2), 536–543. Chappard, C., Brunet-Imbault, B., Lemineur, G., Giraudeau, B., Basillais, A., Harba, R., Benhamou, C.L., 2005. Anisotropy changes in post-menopausal osteoporosis: characterisation by a new index applied to trabecular bone radiographic images. Osteoporos. Int. 16, 1193–1202. Cody, D.D., Gross, G.J., Hou, F.J., Spencer, H.J., Goldstein, S.A., Fyhrie, D.P., 1999. Femoral strength is better predicted by finite element models than QCT and DXA. J. Biomech. 32, 1013–1020. Cummings, S.R., Melton, L.J., 2002. Epidemiology and outcomes of osteoporotic fractures. Lancet 359, 1761–1767. Drucker, D.C., Prager, W., 1952. Soil mechanics and plastic analysis for limit design. Q. Appl. Math. 10, 157–165.
444
J. Thevenot et al. / Journal of Biomechanics 47 (2014) 438–444
Duboeuf, F., Hans, D., Schott, A.M., Kotzki, P.O., Favier, F., Marcelli, C., Meunier, P.J., Delmas, P.D., 1997. Different morphometric and densitometric parameters predict cervical and trochanteric hip fracture: the EPIDOS Study. J. Bone Miner. Res. 12 (11), 1895–1902. Duchemin, L., Bousson, V., Raossanaly, C., Bergot, C., Laredo, J.D., Skalli, W., Mitton, D., 2008. Prediction of mechanical properties of cortical bone by quantitative computed tomography. Med. Eng. Phys. 30, 321–328. Eckstein, F., Wunderer, C., Boehm, H., et al., 2004. Reproducibility and side differences of mechanical tests for determining the structural strength of the proximal femur. J. Bone Miner. Res. 19, 379–385. Gnudi, S., Ripamonti, C., Lisi, L., Fini, M., Giardino, R., Giavaresi, G., 2002. Proximal femur geometry to detect and distinguish femoral neck. fractures from trochanteric fractures in postmenopausal women. Osteoporos. Int. 13, 69–73. Gregory, J.S., Testi, D., Stewart, A., Undrill, P.E., Reid, D.M., Aspden, R.M., 2004. A method for assessment of the shape of the proximal femur and its relationship to osteoporotic hip fracture. Osteoporos. Int. 15, 5–11. Gunay, M., Shim, M-B., Shimada, K., 2007. Cost- and time-effective three-dimensional bone-shape reconstruction from X-ray images. Int. J. Med. Rob. Comput. Assisted Surg. 3, 323–335. Kanis, J.A., 2002. Diagnosis of osteoporosis and assessment of fracture risk. Lancet 359, 1929–1936. Keller, T.S., 1994. Predicting the compressive mechanical behaviour of bone. J. Biomech. 27, 1159–1168. Keyak, J.H., Lee, I.Y., Skinner, H.B., 1994. Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometric measures. J. Biomed. Mater. Res. 28, 1329–1336. Keyak, J.H., Rossi, S.A., 2000. Prediction of femoral fracture load using finite element models: an examination of stress- and strain-based failure theories. J. Biomech. 33, 209–214. Keyak, J.H., Skinner, H.B., 1992. Three-dimensional finite element modelling of bone: effects of element size. J. Biomed. Eng. 14, 483–489. Koivumäki, J.E.M., Thevenot, J., Pulkkinen, P., Kuhn, V., Link, T.M., Eckstein, F., Jämsä, T., 2012. CT-based finite element models can be used to estimate experimentally measured failure loads in the proximal femur. Bone 50, 824–829. Langton, C.M., Pisharody, S., Keyak, J.H., 2009. Generation of a 3D proximal femur shape from a single projection 2D radiographic image. Osteoporos. Int. 20 (3), 455–461. Lengsfeld, M., Schmitt, J., Alter, P., Kaminsky, J., Leppek, R., 1998. Comparison of geometry-based and CT voxel-based finite element modelling and experimental validation. Med. Eng. Phys. 20, 515–522. Meunier, P.J., 1993. Prevention of hip fractures. Am. J. Med. 95 (5A), 75S–78S. Michelotti, J., Clark, J., 1999. Femoral neck length and hip fracture risk. J. Bone Miner. Res. 14, 1714–1720. Miller, Z., Fuchs, M.B., Arcan, M., 2005. Trabecular bone adaptation with an orthotropic material model. J. Biomech. 35, 247–256. Partanen, J., Jämsä, T., Jalovaara, P., 2001. Influence of the upper femur and pelvic geometry on the risk and type of hip fractures. J. Bone Miner. Res. 16, 1540–1546. Pulkkinen, P., Eckstein, F., Lochmüller, E-M., Kuhn, V., Jämsä, T., 2006. Association of geometric factors and failure load level with the distribution of cervical vs. trochanteric hip fractures. J. Bone Miner. Res. 21 (6), 895–901. Pulkkinen, P., Jämsä, T., Lochmüller, E.M., Kuhn, V., Nieminen, M.T., Eckstein, F., 2008. Experimental hip fracture load can be predicted from plain radiography
by combined analysis of trabecular bone structure and bone geometry. Osteoporos. Int. 19, 547–558. Pulkkinen, P., Luomala, M., Obregon, V., Jämsä, T., 2005. Gradient-based method for the assessment of trabecular structure from radiographs. IFMBE Proc., 12. Pulkkinen, P., Partanen, J., Jalovaara, P., Jämsä, T., 2004. Combination of bone mineral density and upper femur geometry improve the prediction of hip fracture. Osteoporos. Int. 15, 274–280. Pulkkinen, P., Partanen, J., Jalovaara, P., Jämsä, T., 2010. BMD T-score discriminates trochanteric fractures from unfractured controls, whereas geometry discriminates cervical fracture cases from unfractured controls of similar BMD. Osteoporos. Int. 21 (7), 1269–1276. Rudman, K.E., Aspden, R.M., Meakin, J.R., 2006. Compression or tension? The stress distribution in the proximal femur. Biomed. Eng. Online 5, 12. Schileo, E., Taddei, F., Cristofolini, L., Viceconti, M., 2008. Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and fracture location on human femurs tested in vitro. J. Biomech. 41, 356–367. Schuit, S., van der Klift, M., Weel, A., Laet, Cde, Burger, H., Seeman, E., Hofman, A., Uitterlinden, A., van Leeuwen, J., Pols, H., 2004. Fracture incidence and association with bone mineral density in elderly men and women: the Rotterdam Study. Bone 34, 195–202. Schumann, S., Tannast, M., Nolte, L.P., Zheng, G., 2010. Validation of statistical shape model based reconstruction of the proximal Afemur—a morphology study. Med. Eng. Phys. 32 (6), 638–644. Szulc, P., Duboeuf, F., Schott, A.M., Dargent-Molina, P., Meunier, P.J., Delmas, P.D., 2006. Structural determinants of hip fracture in elderly women: re-analysis of the data from the EPIDOS study. Osteoporos. Int. 17, 231–236. Testi, D., Cappello, A., Sgallari, F., Rumpf, M., Viceconti, M., 2004. A new software for prediction of femoral neck fractures. Comput. Methods Programs Biomed. 75 (2), 141–145. Testi, D., Viceconti, M., Cappello, A., Gnudi, S., 2002. Prediction of hip fracture can be significantly improved by a single biomedical indicator. Ann. Biomed. Eng. 30, 801–807. Thevenot, J., Hirvasniemi, J., Finnilä, M., Pulkkinen, P., Kuhn, V., Link, T., Eckstein, F., Jämsä, T., Saarakkala, S., 2013. Trabecular homogeneity index derived from plain radiograph to evaluate bone quality. J. Bone Miner. Res. , http://dx.doi.org/ 10.1002/jbmr.1987. ([Epub ahead of print]). Thevenot, J., Pulkkinen, P., Koivumäki, J.E.M., Kuhn, V., Eckstein, F., Jämsä, T., 2009. Discrimination of cervical and trochanteric hip fractures using radiographybased two-dimensional finite element models. Open Bone J. 1, 16–22. Viceconti, M., 1998. A comparative study on different methods of automatic mesh generation of human femurs. Med. Eng. Phys. 20, 1–10. Wolff, J., 1892. Das Gesetz der transformation der Knochen, Translation by Maquet P & Furlong R: The Law of Bone Remodelling. Springer-Verlag, Berlin p. 1986. Yang, H., Ma, X., Guo, T., 2010. Some factors that affect the comparison between isotropic and orthotropic inhomogeneous finite element material models of femur. Med. Eng. Phys. 32, 553–560. Zebaze, R., Jones, A., Welsh, F., Knackstedt, M., Seeman, E., 2005. Femoral neck shape and the spatial distribution of its mineral mass varies with its size: Clinical and biomechanical implications. Bone 37, 243–252. Zheng, G., Gollmer, S., Schumann, S., Dong, X., Feilkas, T., González Ballester, M.A., 2009. A 2D/3D correspondence building method for reconstruction of a patient-specific 3D bone surface model using point distribution models and calibrated X-ray images. Med. Image Anal. 13 (6), 883–899.