Medical Engineering & Physics 36 (2014) 1246–1252

Contents lists available at ScienceDirect

Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy

Comprehensive evaluation of PCA-based finite element modelling of the human femur Lorenzo Grassi a , Enrico Schileo b,∗ , Christelle Boichon c , Marco Viceconti d , Fulvia Taddei a a

Laboratorio di Tecnologia Medica, Istituto Ortopedico Rizzoli, Bologna, Italy Laboratorio di Bioingegneria Computazionale, Istituto Ortopedico Rizzoli, Bologna, Italy c Ansys France, Lyon, France d Department of Mechanical Engineering and INSIGNEO Institute, University of Sheffield, Sheffield, UK b

a r t i c l e

i n f o

Article history: Received 22 May 2013 Received in revised form 9 June 2014 Accepted 28 June 2014 Keywords: Statistical shape modelling Femur Bone biomechanics Principal component analysis

a b s t r a c t Computed tomography (CT)-based finite element (FE) reconstructions describe shape and density distribution of bones. Both shape and density distribution, however, can vary a lot between individuals. Shape/density indexation (usually achieved by principal component analysis—PCA) can be used to synthesize realistic models, thus overcoming the shortage of CT-based models, and helping e.g. to study fracture determinants, or steer prostheses design. The aim of this study was to describe a PCA-based statistical modelling algorithm, and test it on a large CT-based population of femora, to see if it can accurately describe and reproduce bone shape, density distribution, and biomechanics. To this aim, 115 CT-datasets showing normal femoral anatomy were collected and characterized. Isotopological FE meshes were built. Shape and density indexation procedures were performed on the mesh database. The completeness of the database was evaluated through a convergence study. The accuracy in reconstructing bones not belonging to the indexation database was evaluated through (i) leave-one-out tests (ii) comparison of calculated vs. in-vitro measured strains. Fifty indexation modes for shape and 40 for density were necessary to achieve reconstruction errors below pixel size for shape, and below 10% for density. Similar errors for density, and slightly higher errors for shape were obtained when reconstructing bones not belonging to the database. The in-vitro strain prediction accuracy of the reconstructed FE models was comparable to state-of-the-art studies. In summary, the results indicate that the proposed statistical modelling tools are able to accurately describe a population of femora through finite element models. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction Biomechanical properties of bones are often investigated by means of three-dimensional (3D) reconstruction and finite element (FE) model generation, usually from computed tomography (CT) data. This is because 3D FE models of bones can reproduce both shape and material properties, which are acknowledged to be key factors in determining their biomechanical behaviour. The variability of shape and densitometry among individuals is high, and the possible changes induced by pathologies may even enhance it. In order to properly and generally answer many research questions it would be helpful to have available large collections of models that can completely describe a whole population

∗ Corresponding author. Tel.: +39 051 6366554; fax: +39 051 6366974. E-mail address: [email protected] (E. Schileo). http://dx.doi.org/10.1016/j.medengphy.2014.06.021 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

(i.e. a group that can be categorized by a unique general determinant, be it ethnicity, pathological status, etc.). Conversely, one of the most common limitations of FE bone modelling studies [1–4] is the limited size of the bone database investigated, which usually cannot be stated to be representative of a population, and may fail to achieve a good statistical power if used to discriminate between different conditions [5]. Unfortunately, the availability of large bone databases is often scarce, due to different reasons (CT scans cannot be executed on volunteer subjects, having access to clinical archives is difficult for privacy reasons, in-vitro cadaver specimens are scarcely available). One alternative is to generate large datasets of bones using statistical models. Limiting the view to the femoral bone, the availability of a database of realistic femora generated from statistical models allowed researchers to overcome limitations due to the small sample size when investigating bone fracture risk [6], analyzing the influence of anatomy on biomechanics

L. Grassi et al. / Medical Engineering & Physics 36 (2014) 1246–1252

[7], or exploring possible correlations between biomechanics and pathology [8]. These datasets are generated starting from reduced-parameter representations of the variability of shape and density distribution [9–12], obtained from different sources (2D or 3D images, 3D models) often using methods based on principal component analysis (PCA) [13–15], on which this study will focus. The reliability of FE models generated from statistical models may be affected, depending on the generation technique, by several factors: (i) the representativeness of the starting database, from which the statistical model is inferred, with respect to the population to be studied; (ii) the adequateness of the number of statistical modes chosen to represent the population model; (iii) the ability of the synthesizing algorithm to replicate the characteristics of a specimen not belonging to the starting database. Moreover, these statistically-generated FE models should be comparable to the state-of-the-art models in terms of accuracy of the results [3,16,17]. A systematic approach towards the evaluation of statistical bone modelling has been recently proposed [18]. However, to the authors’ best knowledge a comprehensive validation in terms of: (i) representativeness of the database, (ii) accuracy in reproducing shape and material properties of specimens not belonging to the starting database, and (iii) accuracy of the produced finite element models, is lacking. The aim of the present work was to evaluate, on a large database of human femoral anatomies derived from CT scans, the ability of a PCA-based statistical modelling algorithm to accurately represent shape, bone mineral density (BMD) distribution, and strain distribution. 2. Material and methods A database (DB) of femoral CT datasets was collected and morphologically characterized. Isotopological meshes with material properties mapped from CT were built for each femur. Shape and BMD indexation procedures (developed by ANSYS) were performed on the mesh database. The completeness of the database was evaluated through a convergence study on the number of indexation modes to be used. The statistical representation of shape and BMD was evaluated (i) through leave-one-out tests to assess the accuracy in reconstructing femora not belonging to the indexation database; (ii) through comparison of simulated vs. in-vitro experimentally measured strains to assess the mechanical reliability of the reconstructed femora. 2.1. Femora DB A large database of bi-lateral whole femur CT datasets, collected for pre-surgical planning of total hip replacement on osteoarthritic patients was available at Rizzoli Orthopaedic Institute, Bologna, Italy. From that database, we identified a collection of 115 femora (44 males, 71 females) that according to an experienced surgeon had no pathological deformities. All CT datasets were obtained with a standardized protocol [19] and densitometrically calibrated [20]. CT voxel resolution ranged from 0.488 × 0.488 × 1.5 mm to 0.781 × 0.781 × 3 mm. Femoral bones were segmented from the CT images using Amira (v4.0, Visage Imaging Inc., USA), and a polygonal geometry in stereolithography file format was obtained for the external contour of each bone. The anatomical variability was characterized on the 3D reconstructed geometries, using an inhouse developed software [21], through the following anatomical descriptors: femoral neck length, femoral head diameter, caputcollum-diaphyseal (CCD) angle (all detailed in [22]), anteversion angle [23], and epicondyle length (defined as the linear distance between medial and lateral epicondyle). Basic descriptive statistics of the measurements conducted on the database are reported in

1247

Table 1 Descriptive anatomical and anthropometrical parameters of the collection of 115 femora used to perform the indexation. Mean (SD) Anatomical parameters Biomechanical length [mm] Neck length [mm] Head diameter [mm] Epicondyle length [mm] Anteversion angle [◦ ] CCD angle [◦ ] Anthropometrical data Age [years] Height [cm] Weight [Kg]

Minimum

Maximum

406 (28) 39 (4) 44 (4) 81 (13) 13 (9) 126 (8)

356 27 36 69 0 104

483 51 52 96 46 145

58 (15) 166 (9) 73 (14)

26 192 50

84 147 118

Table 1. Bone mineral density was evaluated in the femoral neck of all femora from calibrated CT images, yielding a mean volumetric BMD of 0.308 g/cm3 , with maximum and minimum values of 0.542 g/cm3 and 0.151 g/cm3 , respectively. The large range spanned by most measurements shows that our database largely encompasses the anatomical variability reported in [22,24]. 2.2. FE model generation Each femur of the database was meshed, morphing the mesh template presented in [25] (56809 nodes and 298866 10-noded tetrahedral elements, average element edge size 2 mm) to obtain a collection of 115 subject-specific isotopological finite element meshes. The morphing algorithm adopted (developed by ANSYS) is based on radial basis functions and has been recently reported and validated [25]. Material properties were mapped onto each FE model of the database using Bonemat V3 algorithm [26] with the configuration parameters identified in previous validation studies [27]. 2.3. PCA-based modelling of shape A pre-processing and an indexation step were necessary to set both the shape and the BMD models, as explained in details in Section 2.3.1. 2.3.1. Shape pre-processing step A pre-processing step [28] was applied, in which all the femora were converted to left ones (mirroring the right femur anatomies) and normalized in terms of rigid transformations (3 translations and 3 rotations) and scaling. The following operations were iterated: • Calculate a mean bone shape, averaging nodal coordinates for each node in the 115 isotopological meshes. • Adjust the mean shape to a default scaling, orientation and origin. To define the default scaling, the mean bone was scaled so that the distance between two anatomical landmarks used in the mesh morphing process (fovea of the femoral head and most superior point of the greater trochanter) was constant and equal to the mean value of the database of bones. To define the default origin and orientation, the mean shape was translated so that the centroid has (0,0,0) coordinates, and the rotations were adjusted so that the two above cited landmarks lay on the Z axis (fovea) and on the XZ plane (most superior point of the greater trochanter). • Once a mean shape was determined and its position adjusted to the default position, the scaling, origin, and orientation of each bone were optimized to minimize the distance with that mean shape, and at the end of this process a new “mean shape” was calculated.

1248

L. Grassi et al. / Medical Engineering & Physics 36 (2014) 1246–1252

The iteration was repeated until convergence of the mean shape, i.e. when the vector of the differences of mean shape nodal coordinates between iterations j and j−1 had a norm below a threshold that was set to 10−8 .

patient bone (Yi,mean,patient dimension = number of elements of the morphed mesh).

2.3.2. Shape indexation procedure A matrix “A” was built to include the shape of all isotopological FE meshes (115 columns = number of bones; rows = nodal coordinates of each mesh). The nodal coordinates of the mean femur were subtracted from the matrix “A” and a PCA was applied [29]: it consists in determining the singular value decomposition of the matrix (or the eigenvalue decomposition of the covariance matrix) and it leads to the determination of shape modes (the eigenvectors). Each bone of the DB could therefore be approximated through a scaling coefficient, 6 rigid body modes (translations and rotations) and the following expression:

2.5.1. Selection of shape and BMD modes First, the appropriate number of modes to be used for the statistical representation of shape and BMD modes was set. Therefore, shape and BMD indexation were performed on all the femora of the database using an increasing number of modes, until the mean reconstruction error on the whole database was comparable to the resolution of the CT data (for shape) and below 10% of relative percentage error (for BMD). A relative error metric was chosen for the BMD because of the wide range usually spanned by this scalar value within a bone. The relative BMD error was defined as the ratio between the Euclidean norm of the error vector and the norm of the real BMD values vector. The effect of the first three shape modes on the femur geometry was represented by generating femur instances with the minimum and maximum eigenvalues recorded among the projection of all database specimens.

Xmean +

n 

˛i Xi

i=1

where n is the number of modes considered, Xmean is the mean shape and the Xi vector contains the orthonormal modes. In order to project a new femur in the indexed space, a mesh isotopological with the indexed database should be obtained through a mesh morphing procedure. Then the scaling and rigid body modes should be removed, and finally a standard inner product should be applied to get the coefficients that minimize the distance between the target anatomy and the projection output. The coordinates on the mode Xi are determined by:





˛i = Xi × Xpatient − Xmean =





3∗nnodes



Xi [j] × Xpatient − Xmean [j]

j=1

where × denotes the inner product, where Xi is the vector of the ith mode, Xmean is the vector of the mean shape of the initial database of bones, Xpatient is the vector of the morphed mesh of the patient bone after removing the scaling and rigid body modes. (Xj,mean,patient dimension = 3* number of nodes). 2.4. PCA-based modelling of BMD The statistical modelling of the BMD distribution was similar to the shape statistical modelling. A matrix “B” which contains the BMD value of each element for all the femora in the database was built (115 columns = number of bones; rows = BMD of each element in a mesh). For each element ID of the isotopological mesh, a mean BMD value was calculated over the 115 bones. Then for each element of each bone in the population the mean BMD value was subtracted to the element BMD value. A PCA process was finally applied to the resulting matrix. Analogously to the shape projection, the BMD projection procedure for a target femur BMD consisted in first obtaining a morphed mesh, isotopological with the indexed database and including an element-wise BMD definition, and then projecting the vector of the patient BMD distribution in the base of modes (inner product): The coordinate on the ith BMD mode Yi is ˛i , determined by:





˛i = Yi × Ypatient − Ymean =



nelem





Yi [j] × Ypatient − Ymean [j]

j=1

where Yi is the vector of the ith mode, Ymean is the vector of the mean BMD of the initial database of bones (dimension = nelem), Ypatient is the vector of the BMD values in the morphed mesh of the

2.5. Evaluation procedure

2.5.2. Leave one out tests Once the number of indexation modes to be used for projections was set, leave-one-out tests were performed on all the specimens to assess the accuracy of shape and BMD projection in reconstructing femora not belonging to the indexation database. Since the database is made of iso-topological FE meshes, a node-to-node distance optimization criterion was adopted for the shape projection. Mean and peak errors with respect to the model generated from CT with the standard procedure were calculated for each leave-oneout loop. 2.5.3. Strain prediction accuracy. Strain prediction accuracy was checked in eight femoral anatomies, not belonging to the training database, for which invitro experimental strain measurements obtained on the proximal half of the femur are available. Details of the experiments are available in [16], and here briefly summarised. The eight strain-evaluated femurs were four paired couples of fresh-frozen femora coming from four male donors. The donors were older, taller and heavier than the average of the population, but well within the minimum and maximum values observed in the whole population. In terms of anatomical parameters, they covered quite a large span of the observed population variability, but for most parameters they were not centred near the mean of the population (detailed data in Supplementary on-line table). All femora were tested under six reproducible quasi-axial loading conditions, mimicking the range of admissible directions for the hip joint reaction force during daily activities. The distal third was rigidly constrained. The proximal half of each femur was instrumented with 15 triaxial strain gauges, distributed around each aspect of the femoral head, neck, proximal metaphysis, and diaphysis, yielding a total of more than 1400 measurements for the whole sample. Digitisation was available for all femora, and made it possible to accurately identify in the models the loading points and the measurement locations. On the same eight femora, a high strain prediction accuracy using a standard subject-specific FE modelling chain was reported [2]. In the following, we will refer to this FE modelling technique as “standard procedure”. In the present study, the eight femora were modelled both in terms of shape and BMD using the statistical modelling tool. The element-by-element difference in BMD with respect to the original model was computed for each reconstructed bone to provide a map of BMD changes. In order to evaluate the strain prediction accuracy of the projected FE models, the FE analyses described in

L. Grassi et al. / Medical Engineering & Physics 36 (2014) 1246–1252

1249

Table 2 Convergence results (shape): mean and maximum reconstruction errors (node-tonode distance) for the shape indexation, as a function of the number of modes used. No. of modes 5 10 15 20 25 30 35 40 45 50

Mean projection error

Max projection error

Mean [mm]

Max [mm]

Mean [mm]

Max [mm]

1.42 1.04 0.88 0.77 0.68 0.60 0.53 0.48 0.43 0.39

2.63 1.73 1.34 1.27 0.99 0.88 0.75 0.67 0.63 0.56

6.90 5.68 5.02 4.47 4.03 3.63 3.18 2.94 2.70 2.43

11.72 11.08 10.17 8.69 7.27 6.84 6.31 6.11 5.52 4.84

Table 3 Leave one out results (shape): shape reconstruction errors among all “leave one out” tests. Mean, standard deviation and absolute maximum values are listed for both mean and maximum errors recorded on each “left out” femur, as reconstructed through indexation and projection.

Mean value SD Max value

Mean distance [mm]

Max distance [mm]

1.22 0.41 2.75

5.51 2.30 16.02

[27] were replicated (same boundary conditions) in the eight FE projected models. The accuracy of projected FE models with respect to in-vitro measurements was assessed with linear regression between experimental and predicted strains to quantify the global prediction accuracy, and by error metrics (root mean square (RMS) error and peak error of FE predictions). Fig. 1. Femoral shape variations induced by the first three modes. The extreme values assumed by the eigenvalues in all 115 bones of the database were considered (minimum eigenvalue is rendered with wireframe).

3. Results 3.1. Selection of shape and BMD modes The reconstruction error (node-to-node distance) for the shape indexation of the femur database is reported in Table 2 against the number of modes used. When using 50 modes, the mean error was considerably smaller than the resolution of the CT scan images, and the maximum error among all the nodes of all the FE meshes was below 5 mm. As a consequence, 50 modes were used in the leaveone-out tests. The shape variations induced by the first three shape modes are described in Fig. 1. As to the reconstruction error for BMD indexation, the mean relative percentage error (reported in Table 4 against the number of modes used) went below 10% when using 40 modes (mean projection error of 0.04 g/cm3 , corresponding to a maximum estimation error in the elastic modulus of about 1 GPa for 1.4 g/cm3 of BMD, which is typical of very compact bone, when using the

density–elasticity relationship [30]. Consequently, 40 modes were used in the leave-one-out tests. 3.1.1. Leave one out tests The results of leave-one-out tests for all 115 femora are resumed in Table 3 for shape and Table 5 for BMD. The mean error for the node-to-node distance was 1.22 mm, with the largest errors located in the distal condyles region. The mean relative percentage error for BMD was 13% (mean projection error of 0.05 g/cm3 ). 3.1.2. Strain prediction accuracy The element-wise differences in BMD between reconstructed and original femora were in general moderate, though some peaks were observed internally (Fig. 2). Despite the existing BMD differences, experimental strains were predicted accurately (Fig. 3). The

Table 4 Convergence results (density):absolute and relative reconstruction errors for the BMD indexation, as a function of the number of modes used. No. of modes 5 10 20 30 40 50 60 70

Relative error Mean [%] 18 15 13 11 10 9 7 6

Mean projection error

Max projection error

Max [%]

Mean [g/cm3 ]

Max [g/cm3 ]

Mean [g/cm3 ]

Max [g/cm3 ]

27.76 22.58 17.16 14.96 13.72 12.81 12.22 10.61

0.068 0.060 0.050 0.044 0.039 0.034 0.030 0.026

0.086 0.075 0.061 0.052 0.044 0.040 0.036 0.034

0.688 0.636 0.548 0.477 0.414 0.354 0.303 0.257

1.217 0.915 0.866 0.799 0.758 0.730 0.687 0.550

1250

L. Grassi et al. / Medical Engineering & Physics 36 (2014) 1246–1252

Fig. 2. Exemplar case of the element-wise difference (g/cm3 ) of BMD distribution in the proximal femur between an original bone and its database-projected reconstruction, shown for one of the eight femora not belonging to the training database that were used for the assessment of strain prediction accuracy. The same plot for all eight femora used to assess the strain prediction accuracy is provided for completeness as supplementary on-line material. Table 5 Leave one out results (density): BMD reconstruction error among all “leave one out” tests. Mean, standard deviation and absolute maximum values are listed for both relative and absolute errors recorded on each “left out” femur, as reconstructed through indexation and projection.

Mean value SD Max value

Relative error [%]

Mean projection error [g/cm3 ]

Max projection error [g/cm3 ]

13 2 21

0.05 0.01 0.09

0.65 0.14 1.23

determination coefficient (R2 ) was 0.90, with slope and intercept close to 1 and 0, respectively. The percentage root mean square error was about 10%. 4. Discussion A PCA-based statistical modelling algorithm was proposed and evaluated on a large database of human femoral anatomies derived from CT scans, to see if it can accurately describe and reproduce shape, BMD distribution, and the consequent biomechanical behaviour of a population of bones. A set of 50 modes for shape and 40 modes for BMD was adopted for the indexation in order to keep reconstruction errors well below the pixel size for shape and 10% for BMD. This result is comparable to those reported for similar approaches in literature [14,31,32]. While it is true that the indexation of 50 modes implies a notable

computational weight, this has to be performed only once, while projection, and even more model synthesis, is computationally much less demanding. It has also to be pointed out that, to limit computation time, a node-to-node distance was computed, which is likely to overestimate distances with respect to a node-to-surface normal distance, which could have led to a lower number of shape modes. It is interesting to see that even limiting to the first three modes, although some characteristic deformations emerged (for example a generalized torsion leading to changes in anteversion angle in the second mode, and a generalized bending in the frontal plane leading to changes in caput-collum-diaphyseal angle in the third mode, see Fig. 1), the interpretation of shape changes is not straightforward in terms of parameters derived from traditional anatomy. The interpretation gets even more difficult for subsequent modes, so that a doubt may arise: are higher modes just representing noise? This is not likely, since projection error steadily decreases. However, a non-negligible amount of noise is surely brought in the analysis at higher modes, due to the non-avoidable reconstruction error in the source information, related to CT image resolution, and to the image segmentation process. Perhaps the correct question would be: can a good balance between accurate femur representation and noise exclusion be achieved? In the authors’ opinion the answer depends on the model purpose. The present work sets a higher boundary (50 shape modes and 40 BMD) for a reconstruction of shape and BMD that is good enough to accurately predict superficial strains. An interesting study would be to search for the simplest statistically–generated model (i.e. the lowest number of modes) needed to obtain accurate strain estimates, but this was out of the scope of our work. Using 50 modes, shape was reconstructed in out-of-database bones with an average accuracy of 1.22 mm, which is slightly higher than the average CT pixel value that was assumed as a reference. This may indicate that the actual collection of 115 femur anatomies is not yet sufficient to fully reproduce the anatomical and densitometric variability of the reference population. However, our database is already one of the largest reported [14,31,32]. Future works will try to further increase the number of specimens in the database in order to improve the description of population variability, and to investigate if gender- or age-specific models can reduce noise or variability. Maximum errors for shape reconstruction (up to around 5 mm) were located in the condyles region, but a good accuracy was obtained for biomechanically relevant regions such as femoral

Fig. 3. Comparison of the strains experimentally measured on eight femora (six different quasi-axial loading conditions) with the strains predicted on the corresponding synthesized bones.

L. Grassi et al. / Medical Engineering & Physics 36 (2014) 1246–1252

shaft and femoral neck (from preliminary elaborations on a subset of proximal femora, the maximum error ranged from 1.3 to 2.8 mm, and the mean error was approximately 0.5 mm). Using 40 modes BMD was reconstructed in out-of-database bones with an average error of 0.05 g/cm3 , corresponding to less than 10% for compact bone (BMD > 0.6, corresponding to a wet apparent density of 1 g/cm3 according to [2]), an amount not likely to significantly modify FE models results when looking at superficial bone strains. However, looking more in detail to the BMD maps of reconstructed femora at the proximal femur, we found that while low errors were diffusely reported on the surface, small peak error areas were present internally (Fig. 2), likely due to element distortions in the morphing process. Despite the above reported increase in the shape-reconstruction error and the reported peak BMD errors, good performances were obtained when computing the in-vitro strain prediction accuracy in a set of bones not belonging to the training database (Fig. 3). The obtained results were encouraging, showing a correlation of 0.9 with respect to the experimental measurements, with only a very mild degradation of all indicators with respect to the accuracy results of the reference procedure [2] and of other state-of-the-art models [33,34]. All the reported findings suggest that in the proximal femur the proposed statistical modelling technique is robust and well corroborated. Therefore, it could be prospectively used to study the variables affecting the proximal femoral fracture, the design of hip prostheses or fixation devices, or to characterize the degeneration of femoral head in hip osteoarthritis. Further studies may be needed to better characterize the distal epiphysis, and therefore address knee-related designs and pathologies. A preliminary exploration revealed that caution and further studies are needed also to assess the ability of the indexation tool to produce databases of realistic femur instances from random combinations of eigenvalues. In fact, when eigenvalues were set for all shape modes to extreme (10th and 90th percentile) but also to 1st and 3rd quartile values among those extracted from the 115 femora database, unrealistic modifications of the mean anatomy were produced. Similar results were obtained when using the extreme eigenvalues recorded for the first ten modes, so the observed unrealistic deformations should not be due to the use of noise-related modes. The fact that not all linear combinations of modes can bring to realistic bone appearance highlights the complexity of the highly multidimensional space of bone indexation, and points to the need of finding appropriate metrics or sub-space partitions. Another aspect that would need to be studied, and whose lack represents a limitation of the present work, is the coupling of shape and BMD variations. A relationship between the two is likely to exist, but since shape and BMD have different units (distance vs. density) normalization coefficients or relative weights would be needed in the indexation. The choice to address separately shape and BMD was thus taken for simplicity, and to exclude any bias. Looking at the problem from a broader perspective, the extent to which shape and BMD are coupled is a complex issue. In a population, small and big femora irrespective of high or osteoporotic density can coexist, as well as curved and straight ones, etc. In an individual, shape and BMD evolution due to aging and osteoporosis seem to be largely decoupled. Shape and BMD coupling/decoupling warrants further investigation, which was however out of the scope of the present work. Another limitation is the 3D-to-3D nature of this work, as opposed to approaches that start from images to build models [35,36], or to 2D-to-3D approaches [37,38]. The rationale of the present study was to robustly define the pre-requisites needed for (i) characterizing bones (or populations of bones) through their shape/BMD indexation/projection properties, and (ii) producing new bone instances from a shape/BMD indexation space,

1251

rather than to find a new method to derive models from imaging data, or to enrich the information available from 2D imaging sources. Nevertheless, the proposed technique could be in the future adapted to attempt a 3D reconstruction based on projection of two-dimensional data (e.g. 2D representations of the proximal femur from radiographs). To the authors’ knowledge this is the first work reporting a combined validation of shape and BMD reconstruction accuracy on a large database of clinical in-vivo CT scans, and an in-vitro validation with respect to experimental measurements for a sub-sample of the population. Moreover, the database comprised a high variability of volumetric bone mineral densities at the femoral neck, encompassing normal and osteoporotic subjects, and thus corroborating the generality of the proposed technique. Theoretically, the proposed technique is not even limited to the femoral bone, and could be applied to other bones with no changes but the anatomical landmarks to be identified for the mesh morphing. In summary, the results of the present study indicate that the proposed shape and BMD indexation algorithms are comprehensively able to reliably describe a population of femora through accurate finite element models. Declarations None. Funding This work was partially funded through the EU grant VPHOP (FP7-ICT #223865), and from the Italian Program of Donation for Research “5 per mille”, years 2009 and 2010. Ethical approval Not required. Acknowledgements The authors would like to thank Mauro Ansaloni for the valuable help in the data processing, and especially Ilaria Palmadori for her help in the revision of the manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.medengphy. 2014.06.021. Conflict of interest statement None declared. References [1] Keyak JH, Skinner HB, Fleming JA. Effect of force direction on femoral fracture load for two types of loading conditions. J Orthop Res 2001;19:539–44. [2] Schileo E, Dall’ara E, Taddei F, Malandrino A, Schotkamp T, Baleani M, Viceconti M. An accurate estimation of bone density improves the accuracy of subjectspecific finite element models. J Biomech 2008;41:2483–91. [3] Bessho M, Ohnishi I, Matsumoto T, Ohashi S, Matsuyama J, Tobita K, Kaneko M, Nakamura K. Prediction of proximal femur strength using a CT-based nonlinear finite element method: differences in predicted fracture load and site with changing load and boundary conditions. Bone 2009;45:226–31. [4] Speirs AD, Heller MO, Duda GN, Taylor WR. Physiologically based boundary conditions in finite element modelling. J Biomech 2007;40:2318–23. [5] Radcliffe IAJ, Prescott P, Man HS, Taylor M. Determination of suitable sample sizes for multi-patient based finite element studies. Med Eng Phys 2007;29:1065–72.

1252

L. Grassi et al. / Medical Engineering & Physics 36 (2014) 1246–1252

[6] Bryan R, Nair PB, Taylor M. Use of a statistical model of the whole femur in a large scale, multi-model study of femoral neck fracture risk. J Biomech 2009;42:2171–6. [7] Baldwin MA, Langenderfer JE, Rullkoetter PJ, Laz PJ. Development of subject-specific and statistical shape models of the knee using an efficient segmentation and mesh-morphing approach. Comput Methods Programs Biomed 2010;97:232–40. [8] Bredbenner TL, Eliason TD, Potter RS, Mason RL, Havill LM, Nicolella DP. Statistical shape modeling describes variation in tibia and femur surface geometry between Control and Incidence groups from the osteoarthritis initiative database. J Biomech 2010;43:1780–6. [9] Cootes TF, Taylor CJ, Cooper DH, Graham J. Active shape models-their training and application. Comput Vision Image Understanding 1995;61:38–59. [10] Davies R, Twining C, Taylor C. Statistical models of shape: optimization and evaluation. London: Springer; 2008. p. 302. [11] Lebailly F, Lima LVPC, Clairemidi A, Aubert B, Guerard S, Chaibi Y, Guise Jde, Fontaine C, Skalli W. Semi-automated stereoradiographic upper limb 3D reconstructions using a combined parametric and statistical model: a preliminary study. Surg Radiol Anat 2012;34:757–65. [12] Väänänen SP, Isaksson H, Julkunen P, Sirola J, Kröger H, Jurvelin JS. Assessment of the 3-D shape and mechanics of the proximal femur using a shape template and a bone mineral density image. Biomech Model Mechanobiol 2011;10:529–38. [13] Whitmarsh T, Fritscher KD, Humbert L, Rio Barquero LM Del, Roth T, Kammerlander C, Blauth M, Schubert R, Frangi AF. A statistical model of shape and bone mineral density distribution of the proximal femur for fracture risk assessment. In: Proceedings of MICCAI. 2011. p. 393–400, 14. [14] Bryan R, Mohan PS, Hopkins A, Galloway F, Taylor M, Nair PB. Statistical modelling of the whole human femur incorporating geometric and material properties. Med Eng Phys 2010;32:57–65. [15] Kozic N, Weber S, Büchler P, Lutz C, Reimers N, González Ballester MA, Reyes M. Optimisation of orthopaedic implant design using statistical shape space analysis based on level sets. Med Image Anal 2010;14:265–75. [16] Schileo E, Taddei F, Malandrino A, Cristofolini L, Viceconti M. Subject-specific finite element models can accurately predict strain levels in long bones. J Biomech 2007;40:2982–9. [17] Trabelsi N, Yosibash Z, Wutte C, Augat P, Eberle S. Patient-specific finite element analysis of the human femur-A double-blinded biomechanical validation. J Biomech 2011;44:1666–72. [18] Bonaretti S. Statistical models of shape and density for population-based analysis of bone mechanics with applications to fracture risk assessment and implant design. In: PhD thesis. University of Bern; 2011. [19] Lattanzi R, Baruffaldi F, Zannoni C, Viceconti M. Specialised CT scan protocols for 3-D pre-operative planning of total hip replacement. Med Eng Phys 2004;26:237–45. [20] Kalender WA, Felsenberg D, Genant HK, Fischer M, Dequeker J, Reeve J. The European Spine Phantom—a tool for standardization and quality control in spinal bone mineral measurements by DXA and QCT. Eur J Radiol 1995;20:83–92. [21] Viceconti M, Taddei F, Montanari L, Testi D, Leardini A, Clapworthy G, Sint Jan SVan. Multimod data manager: a tool for data fusion. Comput Methods Programs Biomed 2007;87:148–59.

[22] Noble PC, Box GG, Kamaric E, Fink MJ, Alexander JW, Tullos HS. The effect of aging on the shape of the proximal femur. Clin Orthop Relat Res 1995:31–44. [23] Murphy SB, Simon SR, Kijewski PK, Wilkinson RH, Griscom NT. Femoral anteversion. J Bone Joint Surg Am 1987;69:1169–76. [24] Sugano N, Noble PC, Kamaric E. A comparison of alternative methods of measuring femoral anteversion. J Comput Assist Tomogr 1998;22:610–4. [25] Grassi L, Hraiech N, Schileo E, Ansaloni M, Rochette M, Viceconti M. Evaluation of the generality and accuracy of a new mesh morphing procedure for the human femur. Med Eng Phys 2011;33:112–20. [26] Taddei F, Schileo E, Helgason B, Cristofolini L, Viceconti M. The material mapping strategy influences the accuracy of CT-based finite element models of bones: an evaluation against experimental measurements. Med Eng Phys 2007;29:973–9. [27] Schileo E, Taddei F, Cristofolini L, Viceconti M. 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 2008;41:356–67. [28] Cootes TF, Taylor CJ, Cooper DH, Graham J. Training models of shape from sets of examples. Training 1992;557:9–18. [29] Boichon C, Rochette M, Schileo E, Grassi L, Taddei F, Viceconti M. Shape indexing of human femora using morphing and principal component analysis. In: Proceedings of the first VPH congress. 2010. p. 3–5. [30] Morgan E. Trabecular bone modulus–density relationships depend on anatomic site. J Biomech 2003;36:897–904. [31] Schuler B, Fritscher KD, Kuhn V, Eckstein F, Link TM, Schubert R. Assessment of the individual fracture risk of the proximal femur by using statistical appearance models. Med Phys 2010;37:2560–71. [32] Belenguer Querol L, Büchler P, Rueckert D, Nolte LP, González Ballester MA. Statistical finite element model for bone shape and biomechanical properties. Med Image Comput Comput Assist Interv 2006;9:405–11. [33] Bessho M, Ohnishi I, Matsuyama J, Matsumoto T, Imai K, Nakamura K. Prediction of strength and strain of the proximal femur by a CT-based finite element method. J Biomech 2007;40:1745–53. [34] Trabelsi N, Yosibash Z, Milgrom C. Validation of subject-specific automated p-FE analysis of the proximal femur. J Biomech 2009;42: 234–41. [35] Fritscher KD, Grünerbl A, Schubert R. 3D image segmentation using combined shape-intensity prior models. Int J Comput Assisted Radiol Surg 2007;1:341–50. [36] Schmid J, Magnenat-Thalmann N. MRI bone segmentation using deformable models and shape priors. In: International conference on medical image computing and computer-assisted intervention. 2008. p. 119–26, 11. [37] Langton CM, Pisharody S, Keyak JH. Generation of a 3D proximal femur shape from a single projection 2D radiographic image. Osteoporos Int 2009;20:455–61 (A journal established as result of cooperation between the European Foundation for Osteoporosis and the National Osteoporosis Foundation of the USA). [38] Väänänen SP, Jurvelin JS, Isaksson H. Estimation of 3D shape, internal density and mechanics of proximal femur by combining bone mineral density images with shape and density templates. Biomech Model Mechanobiol 2012;11:791–800.