Automatic Generation of Personalised Skeletal Models of the Lower Limb from Three-Dimensional Bone Geometries

The generation of personalised and patient-specific musculoskeletal models is currently a cumbersome and time-consuming task that normally requires several processing hours and trained operators. We believe that this aspect discourages the use of computational models even when appropriate data are available and personalised biomechanical analysis would be beneficial. In this paper we present a computational tool that enables the fully automatic generation of skeletal models of the lower limb from three-dimensional bone geometries, normally obtained by segmentation of medical images. This tool was evaluated against four manually created lower limb models finding remarkable agreement in the computed joint parameters, well within human operator repeatability. The coordinate systems origins were identified with maximum differences between 0.5 mm (hip joint) and 5.9 mm (subtalar joint), while the joint axes presented discrepancies between 1° (knee joint) to 11° (subtalar joint). To prove the robustness of the methodology, the models were built from four datasets including both genders, anatomies ranging from juvenile to elderly and bone geometries reconstructed from high-quality computed tomography as well as lower-quality magnetic resonance imaging scans. The entire workflow, implemented in MATLAB scripting language, executed in seconds and required no operator intervention, creating lower extremity models ready to use for kinematic and kinetic analysis or as baselines for more advanced musculoskeletal modelling approaches, of which we provide some practical examples. We auspicate that this technical advancement, together with upcoming progress in medical image segmentation techniques, will promote the use of personalised models in larger-scale studies than those hitherto undertaken.


Introduction
Musculoskeletal models have proven to be powerful computational tools to study muscle function and internal forces in healthy (Hamner et al., 2010;Saxby et al., 2016) and clinical populations (Barber et al., 2017;Fox et al., 2018;Montefiori et al., 2019b). Recent technical progress in predictive simulation approaches Falisse et al., 2019) has enabled the investigation of "what if?" scenarios that could support planning and execution of physical interventions. However, applications of personalised medicine often require highly accurate representations of the anatomy of the musculoskeletal system based on medical images such as magnetic resonance imaging (MRI) or computed tomography (CT) scans. For example, personalized bone geometries are essential in orthopaedics for planning surgeries and designing personalised surgical equipment (Clarke et al., 2018;Victor and Premanathan, 2013) and for creating subject-specific musculoskeletal models.
Previous studies presented methods to generate subject-specific models of the entire lower limb (Marra et al., 2015;Modenese et al., 2018) or individual joints Brito da Luz et al., 2017;Montefiori et al., 2019a;Nardini et al., 2020) and dedicated modelling tools like NMSBuilder (Valente et al., 2017a) or specialized features in the AnyBody software (Damsgaard et al., 2006) are available to implement these workflows. Nevertheless, patient-specific musculoskeletal models are currently employed in small-sized clinical applications (Falisse et al., 2020;Montefiori et al., 2019b;Pitto et al., 2019;Taddei et al., 2012;Valente et al., 2017b), mostly because the generation of each model is a time-demanding operation requiring manual intervention by specialized operators. For example, a codified approach proposed in Modenese et al. (2018) reported around 10 hours to build a complete bilateral musculoskeletal model of the lower limbs from segmented bone geometries (around two hours to create an ipsilateral skeletal model), while Scheys et al. (2006) reported on average 65 minutes to define the lower limb musculature using an atlas-based semi-automated approach. We believe that validated and fully automatic workflows are of paramount importance to enable largescale use of these computational models.
Multiple studies with orthopaedic focus have explored the possibility of defining anatomical coordinate systems (ACSs) in the lower extremity bones based on key geometrical features. Miranda et al. (2010) and Rainbow et al. (2013) proposed automatic methods for defining ACSs for the distal femur, proximal tibia and patella, that showed minimum variability with the bone morphology. Kai et al. (2014) developed an automatic approach to identify the reference systems of the pelvis, femur and tibia based on principal axes of inertia, principal component analysis and longitudinal slicing, obtaining ACSs compatible with those created by human operators, except for the pelvis where anterior tilt was up to 18.8° higher. More recently, Renault et al. (2018) proposed multiple algorithms based on the automatic identification of the articular surfaces at the hip and knee joints, showing high repeatability of these methods when applied on 24 CT scans by three operators. However, in these previous works the ACSs were not defined consistently across publications and none of these methods has been employed for creating articulated skeletal models yet.
Statistical shape modelling workflows have recently demonstrated high potential for reconstructing bone geometries from sparse anatomical datasets (Davico et al., 2019;Nolte et al., 2016;Suwarganda et al., 2019) and landmarks digitized in the gait lab (Nolte et al., 2020;Zhang et al., 2016), but to the best of the authors' knowledge they do not yet offer methods to generate articulated skeletal models of the complete lower limb. The bone reconstructions are limited to the long bones (Nolte et al., 2020;Nolte et al., 2016) or omit the talus and foot bones (Davico et al., 2019;Suwarganda et al., 2019;Zhang et al., 2016), and in musculoskeletal modelling contexts they have been employed to perform non-linear scaling of pre-existing muscle attachments (Nolte et al., 2016) with scarce focus towards joint modelling. Hence, a comprehensive approach to generate entire lower limb models from personalised bone geometries is still missing.
The aim of this paper is to present a tool to create models of the complete lower limb from threedimensional bone geometries in a fully automatic way. The tool implements a sequence of operations normally performed in manual workflows, executes them in negligible computational time and generates models usable immediately in kinematic and kinetic analyses or employable as baselines for fully featured musculoskeletal models. The models produced by this tool were evaluated against manually created models employed in previous research and the joint parameters computed by competing algorithms were compared to assess their interchangeability. Differences in kinematics and kinetics curves calculated using manual and automated models were also preliminarily quantified on a set of gait simulations. Examples of further technical developments, such as joint articular mechanics and integration with anatomical models of musculature, are finally provided to demonstrate the tool's potential for enabling large-scale studies and broader musculoskeletal research.

Workflow to generate automatic skeletal models
A set of computational methods proposed in previous literature to define ACSs automatically (Kai et al., 2014;Miranda et al., 2010;Renault et al., 2018) were acquired and included in a more extensive modelling workflow. The geometrical methods from Renault et al. (2018) were obtained from the public "GIBOC-KNEE" MATLAB toolbox (https://github.com/renaultJB/GIBOC-Knee-Coordinate-System) and extensively modified and expanded, while the methods described by Kai et al. (2014) were independently reimplemented and those of Miranda et al. (2010) were obtained through contacting the authors of the publication. Additional algorithms were developed ad hoc for the purposes of this investigation ( Table 1).
The implemented workflow ( Figure 1) consisted of the following steps: a) obtaining segmented threedimensional bone geometries of the pelvis, femur, patella, tibia, fibula, talus, calcaneus and the other foot bones from medical images; b) automatically processing these bone models to extract the geometrical parameters required to define ACSs and appropriate joint coordinate systems (JCSs) for the parent and child bodies of each joint of the lower limb; c) creating an articulated skeletal model of the lower limb in OpenSim format (Delp et al., 2007) using the identified JCSs and ACSs.
In step b), bone geometries were analysed starting with a transformation to the ACSs defined by their principal axes of inertia (Gonzalez-Ochoa et al., 1998;Mirtich, 1996), followed by bone-specific features extraction. The complete list of algorithms available to define each JCS is reported in Table 1 and the details of the methodologies are described in their reference publications, and for the newly developed algorithms for the pelvis, talus and foot, in the supplementary materials. The articulated skeletal models were generated leveraging the MATLAB (The MathWorks, Natick, MA, USA) application programming interface (API) of OpenSim 4.1 (Seth et al., 2018): a rigid body with appropriate inertial properties (McConville et al., 1980;Winter, 2009), was created for each leg segment and joints defined based on the JCSs identified at step c). The individual JCSs, consistent with Modenese et al. (2018), are described in Table 1. For convenience in the validation step, all rigid bodies shared a local coordinate system coincident with that of the medical images, as they would have in a model generated using NMSBuilder (Valente et al., 2017a). The lower limb models included five bodies: pelvis, femur, tibia (including fibula and a rigidly attached patella), talus and foot (including calcaneus and foot bones), and five joints: a free joint between pelvis and ground (6 degrees of freedom, or "DoF"), a ball and socket joint for the hip joint (3 rotational DoF) and hinge joints for the tibiofemoral, talocrural and subtalar joints (1 rotational DoF each). No explicit patellofemoral joint was included in the models. The models also included fourteen landmarks (Table 2) automatically identified on the bone surfaces and intended for registration with the skin markers used in standard gait analysis.
The entire set of scripts implementing this workflow was organized in a MATLAB toolbox named STAPLE (Shared Tools for Automatic Personalised Lower Extremity modelling).

Evaluation of the automatic models
The models produced using the automatic workflow were compared against musculoskeletal models of the lower limb generated for other purposes using NMSBuilder (manual models) and previously employed in published research (Montefiori et al., 2019b) or contributions at international conferences (Modenese et al., 2020;Modenese et al., 2019). These subject-specific models were built following the codified approach of Modenese et al. (2018) and using bone geometries available in public datasets plus an in vivo MRI dataset collected with the approval of the Imperial College Research Ethics Committee. A complete description of the datasets, characterized by quality of the bone geometry meshes ranging from very good (LHDL-CT) to low (JIA-MRI), is provided in Table 3 and Articulated skeletal models for all datasets were generated using the STAPLE toolbox (automatic models) and their JCSs compared against those of the manual models, reporting differences in their origin location and axes orientation. These quantities were evaluated in the common global coordinate system of the medical images. The automatic models were created using algorithms that matched those of the manual approach (pelvis: STAPLE-pelvis, femur: GIBOC-Cylinder, tibia: Kai-Tibia, talus: STAPLE-Talus and foot: STAPLE-Foot).

Comparison of joint parameters estimated by different algorithms
The JCSs of those joints for which more than one algorithm was available (ground-pelvis and knee joint) were then calculated using all the available options and the resulting JCSs compared to those employed in the evaluation part of the study, used as reference. Linear distances between origins and angular differences between axes were quantified and expressed in the JCS of the reference algorithm.

Gait simulations using automatic and manual models
The sensitivity of joint angles and net joint moments to JCS differences was quantified simulating six gait trials of JIA-MRI (Montefiori et al., 2019b) with the correspondent automatic and manual models, using the OpenSim inverse kinematic and inverse dynamic analyses. The resulting sets of curves were then compared using a Statistical Parametric Mapping (SPM) two-tailed t-test (significance level:

Results
All the automatic models employed in the study were successfully generated in less than 30s each using a standard Z640 Dell Workstation (RAM: 64 GB, CPU: 2 Intel Xeon E5-2630 2.40 GHz).
The comparison of automatic and manual models (Table 4) resulted in an overall strong similarity of the joint parameters across all considered datasets. The hip and talocrural joint centres were in excellent agreement with the manual estimations (maximum differences hip: 0.5 mm, talocrural: 1.2mm), whereas the maximum difference was 2.5 mm at the knee joint and 5.9 mm at the subtalar joint in the JIA-MRI model due to the low quality of bone reconstruction (difference range in the other datasets: 0.4-3.2 mm). The cranial-caudal position of the pelvis-ground joint origin differed by up to 4.9 mm (range: 0.8-4.9 mm) causing minor differences in pelvic tilt (range: 1.8°-3.6°, Figure 3-A) that propagated to the hip-parent JCS. The axes of the knee and talocrural hinge joints (medio-lateral Z axes) were estimated with maximum differences from manual models of 1.0° degree, while the subtalar joint axes presented maximum differences up to 2.9° in the datasets with good quality bone geometries but reached 11.3° in the JIA-MRI model (Figure 3-B).
When comparing competing algorithms, we observed that differences among their JCSs were not always negligible (Table 5). The Kai-Pelvis algorithm presented larger pelvic tilt offsets than STAPLE-Pelvis, while GIBOC-Spheres provided the knee-parent JCSs closest to the reference algorithm (differences <1° for all cases except JIA-MRI). Kai-Femur and GIBOC-Ellipsoid resulted in more posteriorly and anteriorly located JCSs respectively, with maximum angular differences for the knee joint axis up to 4.7° (LHDL-CT) for the former and 5.2° (TLEM2-CT) for the latter. At the proximal tibia, all GIBOC algorithms computed more proximal origins than Kai-Tibia's (range: 5.8-12.5 mm), with angular differences in the range 0.5°-11.5° for the medio-lateral axis but smaller for the proximal-distal axis (range: 0.9°-2.9°). Overall, GIBOC-Plateau and Miranda-Tibia, which failed processing TLEM2-CT, identified similar JCSs, as expected by very similar algorithms. The JCSs from GIBOC-Ellipse and GIBOC-Centroids were also found more similar to each other than to the Kai-Tibia reference.
The joint angles from the walking simulations performed with the manual and automatic JIA-MRI models ( Figure 4) presented correlation coefficients greater than 0.99 (p<0.0001) and RMSE ≈ 1° degree for all coordinates except pelvis tilt and hip flex/extension, for which these errors (3.4°) were consistent with the JCS offsets in the sagittal plane (Table 4). The SPM t-test returned significant differences between the kinematics curves only for small portions of the gait cycle (from 0% to 13%), except for pelvis tilt and hip flex/extension (≈100%) and pelvic rotation (32%). The joint moments ( Figure 5) had correlation coefficients larger than 0.94 (p<0.0001) and mean RMSE ranging from 0.078 (hip flex/extension) to 0.007 Nm/Kg (ankle dorsi/plantarflexion). Significant differences between the curve sets were observed mostly in the swing phase of gait, where the inertial effects were more pronounced, in clusters ranging from 11% to 37% of the gait cycle. Detailed SPM plots are available in the supplementary materials.

Discussion
The aim of this paper was to present a tool to automatically create personalised models of the lower limb from three-dimensional bone geometries segmented in medical images. To evaluate the proposed methodology, we generated four automatic models and compared them against models manually created from the same data in other research projects. We found that the automatic and manual models were remarkably similar (Table 4), with largest differences observed for the pelvis-ground and subtalar joints. In the pelvis, the JCS origin was misplaced by up to 4.9 mm cranially due to the identified bony landmarks (Figure 3-A). This difference causes a systematic anterior tilt offset in the range 1.8°-3.6° that propagates to the pelvis and hip joint kinematics, as confirmed by the gait simulations. Although not negligible, this offset represents a substantial improvement (almost 15°) compared to the results of Kai et al. (2014). More recent algorithms, e.g. Fischer et al. (2019), could be considered for future additional comparisons. The automatic subtalar joint axis in the JIA-MRI dataset was also noticeably different from the equivalent manual model (11.3°, Figure 3-B). This discrepancy was attributed to the low-quality talus bone reconstruction since in all the other models the same axis was estimated within 2.9°.
JCS differences at the hip, knee and talocrural joint were well within the ranges of human interoperator repeatability reported as standard deviations in previous investigations (Hannah et al., 2017;Martelli et al., 2014;Montefiori et al., 2019b) and even the largest difference observed at the subtalar axis was close to the maximum inter-operator variability (9.6°) reported by Montefiori et al. (2019a).
These JCSs differences caused joint angles offsets in the gait simulations, as expected based on previous studies (Kainz et al., 2016;Martelli et al., 2014;Montefiori et al., 2019b;Valente et al., 2014). The largest difference in net joint moments ( Figure 5) occurred for the hip flex/extension moment in the swing phase of gait and was attributed to dissimilar shank mass between models (25% difference), based on the results of Wesseling et al. (2014). This discrepancy was due to the inertial properties being estimated using segmented soft tissues geometries (not publicly available) in the manual JIA-MRI model (Montefiori et al., 2019b) and using regression equations (McConville et al., 1980;Winter, 2009) in the automatic model. It is worth highlighting that simulation results from the STAPLE models will have zero inter-and intra-operator variability due to model construction.
In all considered models, the JCSs were created as in Modenese et al. (2018), but other algorithms were also assessed (  (2018). At the distal femur the differences among JCSs were small (<5° in the 93% of estimations) but not negligible, therefore the choice of the algorithm, e.g. fitting ellipsoids (Sholukha et al., 2011) or spheres (Yin et al., 2015) to the femoral condyles articular surface, must be justified with careful functional anatomy considerations relevant to the research question. At the proximal tibia, the mechanical axis (Y axis) was similar between Kai-Tibia and GIBOC algorithms (range: 0.9°-2.9°) but less for . Larger differences were observed in the medio-lateral axis (range: 0.5°-11.5°) and anterior-posterior axis (range: 1.9°-11.7°). The workflow of Modenese et al. (2018) connects the tibia to the tibiofemoral axis identified in the femur in the segmented pose, hence when the knee joint assumes the neutral position with zero flexion/extension angle, offsets due to standard supine imaging (Hirschmann et al., 2015) could still be present. The automatic algorithms could inform procedures to adjust the tibiofemoral alignment, but further work is required in this respect.
A current limitation of the study is that our validation datasets were only four, and despite their heterogeneity, did not include bone geometries presenting abnormalities, deformities and partial geometries due to imaging (Henckel et al., 2006) or amputations that could be encountered clinically.
This limitation is mitigated by the availability of multiple algorithms for the long bones, some of which specifically developed for incomplete geometries (Miranda et al., 2010;Renault et al., 2018), but additional validation is undoubtedly required. Also, the quality of the bone geometries can vary depending on the medical image modality and their specifications and resolution, potentially affecting the automatic algorithms. This possibility was not systematically investigated in this study, although we can report that, when processing low quality bone models, algorithms based on bone global geometrical features (Kai-algorithms) were generally more robust than those relying on the identification of the articular surfaces (GIBOC-algorithms), while the Miranda-algorithms seemed the most affected and failed to process some of the datasets.
The importance of a fully automatic approach as a step towards clinical application of modelling can be fully appreciated considering the reduction in total processing time for generating a skeletal model, from segmentation of the medical images to the final OpenSim model. For the ipsilateral models considered in this study, the bone segmentation and manual modelling steps required a comparable time (between two and four hours each, depending on the experience of the operator). For example, an expert operator segmented the ICL-MRI dataset, the only one that we processed directly, in around two hours using the semi-automatic functionalities of ITK-Snap (Yushkevich et al., 2006) and created an ipsilateral skeletal model in roughly the same amount of time. Based on the TLEM2 segmentation time (Table 3), our previous experience (Renault et al., 2018) and literature reports (Matsiushevich et al., 2019), we estimated similar processing times also for CT datasets. STAPLE executes modelling workflows in seconds, therefore practically halving the total processing time to generate a skeletal model and reducing it essentially to an image segmentation task. This means, for example, that when radiological scans are collected and segmented for planning musculoskeletal surgical interventions, obtaining further model-based biomechanical analyses becomes a quick and inexpensive option.
Considering that deep learning techniques are reducing also the segmentation time by order of magnitudes, e.g. Noguchi et al. (2020) reported ~12s for segmenting a full-body CT scan of ~600 slices, the generation of personalised lower limb models in a number of clinical applications appears technically feasible.
The STAPLE toolbox is intentionally modular and can automate the generation of entire or partial lower limb models (see supplementary materials for examples). Moreover, the GIBOC and STAPLE algorithms provide a large amount of anatomical information such as articular surfaces (Figure 6-A) and bone profiles that can be used for implementing more advanced joint models than those proposed here, e.g. contact models (Brandon et al., 2017;Conconi et al., 2015) or parallel mechanisms (Sancisi and Parenti-Castelli, 2011).
Extending skeletal models with models of muscle anatomy can also be completely automated (Modenese and Kohout, 2020). In a previous contribution (Modenese et al., 2020) we have used a non-rigid iterative closest point registration (Audenaert et al., 2019) to map the muscle attachment areas from a cadaveric dataset to the ICL-MRI participant's bones, generating highly-discretized, personalised muscle representations from segmented muscle geometries and simulating their kinematics during gait (Figure 6-B). Future efforts will focus on streamlining these methodologies towards a comprehensive, fully automated, modelling tool for generating subject-specific musculoskeletal models.
In summary, this work presents a computational tool enabling researchers to generate articulated, subject-specific skeletal models of the lower limb in negligible time through a completely automatic workflow that takes three-dimensional bone geometries as inputs. These models can be used immediately for kinematic and kinetic analyses or can serve as extendable baselines for complete musculoskeletal models including musculotendon actuators. This work is framed in a long-term plan aiming to advance the state of the art of anatomical modelling and promote large-scale clinical adoption of personalised computational models of the musculoskeletal system through complete automation of the most challenging modelling tasks.
All data, models and scripts used in this paper are available for download, as detailed in the Appendix. To facilitate the reproducibility and replication of our results, we have released our research code and data with this publication. STAPLE is openly developed and shared under the Apache License 2.0 via its repository at https://github.com/modenaxe/msk-STAPLE. All of the data and scripts needed to run the calculations reported in this work, as well as the post-processing scripts to reproduce the figures in this paper and the conference abstracts of the authors mentioned in the main text are available at https://github.com/modenaxe/auto-lowerlimb-models-paper. The version of the scripts and data used in this paper are available through Zenodo at [link will be added upon acceptance] and at https://simtk.org/projects/auto-sk-models.

CRediT authorship contribution statement
NOTE: Public links will be activated upon paper acceptance. 1 3 Figure 2 Musculoskeletal models of the lower limb from previous research used for validating the automatically generated skeletal models (first row). The models were built using bone reconstructions of variable quality (second row). Details about these models are available in Table 3. 1 8 Table 1 Joint coordinate systems and available algorithms to calculate their parameters. Kai-algorithms are described in Kai et al. (2014), GIBOC-algorithms are described in Renault et al. (2018), Miranda-algorithms are described in Miranda et al. (2010) and STAPLE-algorithms are described in the Supplementary Materials of this publication. Please note that the ground coordinate system coincides with that of the medical images in all models and that the algorithms used at the pelvis were always applied to bone geometries excluding the sacrum bone (considered challenging to segment in MRI scans). The algorithms applied to the "tibia" rigid body processed the geometry of the proximal tibia, full tibia or full tibia and fibula, depending on the algorithm. Note that the "foot" rigid body, including the geometries of calcaneus and foot bones, is called "calcn" in the OpenSim models for consistency with the other models included in the software distribution. Origin: coincident with knee parent origin. Axes: Z axis aligned with medio-lateral axis of knee parent Y axis is perpendicular to Z lying in the same plane as talocrural-child origin X axis is perpendicular to Y and Z.

Rigid Body
talocrural parent N/A Origin: coincident with the talocrural child. Axes: Z is aligned with the Z axis of the talocrural child. Y axis is perpendicular to Z, lying in the same plane of Z and the knee joint center. X axis is perpendicular to Y and Z. talus talocrural child STAPLE-Talus Origin: point at midpoint of the length on the axis of the cylinder fitter to the talar trochlea articular surface. Axes: Z is the axis of the cylinder fitted to the talar trochlea articular surface. X is perpendicular to Z, lying on a plane parallel to the foot sole XZ plane. Y is perpendicular to X and Z. subtalar parent STAPLE-Talus Origin: center of the sphere fitted to the articular surface of the talocalcaneal joint. Axes: Z axis on the line from the center of sphere fitted to the talocalcaneal articular surface to that fitted to the talonavicular articular surface. Y is perpendicular to Z, lying in the same plane of Z and the knee joint centre. X axis is perpendicular to Y and Z. foot subtalar child N/A Origin and axes defined by subtalar parent.
foot sole (auxiliary) STAPLE-Foot Origin: most distal point of the calcaneus.
Axes: X axis pointing from the most caudal point on the talus to the midpoint of the most caudal points on the 1st and 5th metatarsal heads. Y axis perpendicular to the plane identified by the points defining X. Z axis is perpendicular to X and Y.   Table 4 Differences between the joint coordinate systems in the manual and automatic skeletal models. The automatic models were created using the algorithms STAPLE-Pelvis, GIBOC-Femur, Kai-Tibia, STAPLE-Talus and STAPLE-Foot. Linear distances between the origins of the joint coordinate systems are expressed in the reference system of the medical images, for which axes directions are: Z pointing cranially, Y posteriorly, X to the left for LHDL-CT, TLEM2-CT and JIA-MRI; Z pointing cranially, Y anteriorly and X to the right for ICL-MRI.  Table 5 Comparison among the joint coordinate systems estimated by all the available algorithms and the specified reference algorithm. Components of the displacement vectors are expressed in the joint coordinate system of the reference algorithm. N/S means "not solved" and indicates that the algorithm crashed or computed a manifestly incorrect solution. 3.2 JIA-MRI -3.3 6.0 0.0 6.9 11.7 1.8 11.5 * Miranda-Femur crashed when processing JIA-MRI and fitted the shaft of the femur of ICL-MRI.

Dataset
**Miranda-Tibia identified inverted axes direction for TLEM2-CT and was not able to process the JIA-MRI tibial geometry unless the fibula was included as well.