Validation of a model-based inverse kinematics approach based on wearable inertial sensors

Abstract Wearable inertial measurement units (IMUs) are a promising solution to human motion estimation. Using IMUs 3D orientations, a model-driven inverse kinematics methodology to estimate joint angles is presented. Estimated joint angles were validated against encoder-measured kinematics (robot) and against marker-based kinematics (passive mechanism). Results are promising, with RMS angular errors respectively lower than 3 and 6 deg over a minimum range of motion of 50 deg (robot) and 160 deg (passive mechanism). Moreover, a noise robustness analysis revealed that the model-driven approach reduces the effects of experimental noises, making the proposed technique particularly suitable for application in human motion analysis.


Introduction
Inertial measurement units (IMUs) are becoming popular components of devices that need to measure orientation in space. An IMU is composed by two triaxial sensors: a gyroscope to measure angular velocity and an accelerometer to sense linear accelerations. These quantities are then generally fused together using a sensor fusion algorithm (SFA) to estimate the IMU orientation in a global reference frame. If a triaxial magnetometer is also integrated in the IMU then the global reference frame is Earth-fixed, defined using gravitational and magnetic field directions. The IMU orientation is therefore an estimated quantity, resulting from the fusion of measurements taken in different domains (Ligorio et al. 2016) and naturally affected by undesired experimental factors such as measurement noises, external disturbances, sensor biases, etc. These factors were demonstrated to negatively impact on the accuracy and the reliability of the orientation estimates (Khaleghi et al. 2013).
Preliminary studies on motion tracking using a combination of accelerometers and gyroscopes dates back to the early 70 s (a good historical review is provided in (Picerno 2017)). However, the advent of Micro Electro-Mechanical Systems (MEMS) technology has led to a broad adoption of IMUs in a variety of fields and, in the latest 25 years, opened new frontiers in wearable human motion analysis (Picerno 2017;Cuesta-Vargas et al. 2010).
Human motion analysis was defined by Cappozzo et al. (Cappozzo et al. 2005) as the science that "aims at gathering quantitative information about the mechanics of the musculoskeletal system during the execution of a motor task". Joint kinematics, one of the key descriptors of human motion, are routinely measured in laboratory settings, where a set of stereophotogrammetric cameras track the 3D position of passive reflective markers placed on well-defined subject's bony landmarks (Wu et al. 2002;Cappozzo et al. 1995). In decades of use of stereophotogrammetric systems, experimental protocols (Ferrari et al. 2008), data processing pipelines (Kadaba et al. 1990; Davis et al. 1991), and joint kinematics estimation techniques (Grood and Suntay 1983) contributed to the success of this technology that became the de facto gold standard in biomechanics. Model-based simulations further enhanced estimated joint kinematics accuracy and, at the same time, enabled insight on musculoskeletal function (Arnold et al. 2010). The use of accurate musculoskeletal models that provide kinematics constraints helps to reduce the effects of experimental sources of errors (Duprey et al. 2010;Clement et al. 2015;Lamberto et al. 2016). Moreover, these models enable the estimation of other internal quantities such as muscle length and muscle forces (Delp et al. 2007;Saxby et al. 2016), even in real-time applications (Van den Bogert et al. 2013;Pizzolato et al. 2017). Stereophotogrammetric technology, however, is not the optimal solution for every application. Optoelectronic systems are bulky, expensive, require trained personnel to be correctly and efficiently used, and finally they can only be used in a laboratory environment. The latter issue is critical, since tasks performed inside a laboratory may not reflect real-life movements as in the case of outdoor sports or daily life tasks executed by patients.
For all these reasons, the biomechanics community looks toward IMUs as a potential alternative. IMUs are nowadays wearable, relatively cheap, easy to use and reliable enough for the requirements of most human motion analysis applications. Provided with on-board long-lasting batteries, long range communication protocols and/or on-board data logging features, IMUs could enable continuous kinematics analysis in almost every environment (i.e. clinical, outdoor, daily life, industry).
A common approach in most studies using IMUs to assess human kinematics is the so called "direct kinematics" one. This method directly prescribes the IMU orientations to the segments to which they are attached to in order to estimate the joint kinematics. Another possible approach directly uses the raw measurements provided by each IMU (i.e. angular velocity, linear acceleration and, optionally, magnetic field) together with a kinematic model of the human body inside an extended Kalman filter, to compute joints and segments kinematics (Kortier et al. 2014;Van den Noort et al. 2016;Seel et al. 2014). However, such solutions might be difficult to be used and tuned for other applications by non-experts in filtering techniques. In both these approaches, possible misalignments between sensors and joint axes can be mitigated using ad-hoc calibration procedures (Palermo et al. 2014;Van den Noort et al. 2013). However, to the best of authors' knowledge, very few studies investigated the use of a model-based inverse kinematics estimation approach (Koning et al. 2015;Borb ely and Szolgay 2017;Kortier et al. 2014;Karatsidis et al. 2018). Furthermore, none of them presented a throughout methodological assessment of accuracy and robustness to measurement noise for the proposed algorithms.
This study proposes a model-based Inverse Kinematics (IK) approach to assess the motion of multi-link systems from the orientation of IMUs placed on the constituting bodies. In this approach, the joint constraints included in the model have to be respected when calculating joint kinematics (Lu and O'Connor 1999). General applicability and ease of use motivated the choice of using IMU orientations as input for the developed methodology, even if potentially affected by the inaccuracies previously described, instead of raw sensor data. Furthermore, the proposed approach has been implemented to be model-independent, allowing users to select the most appropriate kinematic model (lower limb, upper limb, spine, etc.) according to their specific needs. Moreover, since musculoskeletal models are essentially chains of rigid bodies, the use of robots or limb-like mechanisms reduces the effects of non-methodological sources of errors when the focus is the assessment of the performances of a new IK approach. In this study, the proposed approach to calculate joint angles from IMUs was evaluated in two experimental scenarios, using respectively a robot and a passive plastic planar mechanism. This choice allows evaluating the proposed methodology without the confounding effect of errors that would have been present in human testing, e.g. soft tissue artifacts, and led to the design of ad hoc test-benches. The developed orientation-based IK is freely available as a plug-in for OpenSim (Delp et al. 2007) at the SimTK project's page 1 .

Inverse kinematics analysis based on orientation data
Orientation-based Inverse kinematics (OB-IK) is similar to the marker-based inverse kinematics (MB-IK) available in OpenSim (Delp et al. 2007), a popular model-driven global optimization procedure that allows to estimate joint kinematics starting from marker data (Lu and O'Connor 1999). The main benefit of an inverse kinematic method, assuming enough experimental kinematic measurements are available to be tracked, is the use of a multi-link mechanical model with associated joint constraints (Kainz et al. 2016). Indeed, the constraints in the model can prevent physiologically unfeasible configurations, such as joint dislocation or joint angles outside their physiologycal ranges of motion, possibly resulting from experimental errors and noises. In order for OB-IK to work, "virtual" orientation sensors should be placed on the model links matching the experimental configuration.
To programmatically retrieve the relative transformation between sensors and segments coordinate frames we developed a calibration procedure which assumes the joint angles to be known in at least one frame, usually in a predefined static configuration. The calibration procedure then locks the model in that pose and moves each virtual orientation sensor around the model segment to which it is attached to so that its orientation coincides with the experimentally measured one. Thanks to this procedure, no accurate experimental alignment between sensors and segment axes (Favre et al. 2006;Liu et al. 2009) is required. This in turns allows also to avoid experimental functional calibration trials (Seel et al. 2014).
Once the calibration phase is completed, the developed OB-IK algorithm takes as input the calibrated model with virtual orientation sensors correctly placed on the segments and the orientations provided by experimental IMUs, expressed as unitary quaternions. This representation allows to minimize data size and, at the same time, avoids the singularity implied in more compact representations (i.e. Euler Angles) (Diebel 2006). The goal of the OB-IK is to calculate the whole-model joint angles that determine the best match between the orientations of the experimental IMUs and those of the corresponding virtual orientation sensors attached to the model ( Figure 1). The algorithm always considers the model as a whole, therefore no local steps considering kinematic subchains are performed. In order to quantify the orientation mismatch between one experimental IMU and its virtual correspondent the Euler axis-angle representation was used and the angle (a) given by this representation chosen as the parameter to minimize. The developed computational tool was based on the implementation available in the Simbody (Sherman et al. 2011) source code.
The minimization problem was defined as a stateless whole-body optimization, where each time frame is solved independently from the previous ones. For each time frame, a gradient descent algorithm iteratively looks for the global minima of a weighted quadratic function of the angular tracking error a. Depending on a priori knowledge of the experimental setup (e.g., sensors' placement or hardware characteristics), a different weighting factor could be assigned to each sensor to represent the level of confidence we expect for its measurements.
Being w i the weighting factor associated to the i-th orientation sensor and a i the orientation mismatch for the i-th pair (i.e. the angular error coming from the Euler angle-axis representation of the relative orientation between the real and the virtual orientation sensor), the cost function to be minimized can be written as In Eq.1 the dependency from the set of the model generalized coordinates q has been made explicit.

Framework 1: validation against encoder measurements
A 6-DoF actuated robotic arm UR-10 (Universal Robots A/S, Denmark) was used in this experimental setup (Figure 2a). Four Cometa WaveTrack IMUs (Cometa Systems, Italy) were positioned on the four links around the three most proximal joints of the robot (i.e. shoulder-pan, shoulder-lift and elbow joints). The desired motion was defined by manually moving the robot's end-effector within the threedimensional reachable workspace, while at the same time measuring the individual joint angles by means of the embedded encoders. The recorded joint trajectories were then prescribed to the robot controller programmed to exactly repeat them four times for each task in order to ensure repeatability of the performed task at varying speed of execution. Two task at two different speeds were recorded. During the manipulation, care was taken to simultaneously involve in the motion all the three joints of interest, while also spanning a wide portion of their ranges of motion (i.e. approximately 90 deg for the shoulder_pan_joint and 45 deg for the others). Data were collected, using a common trigger signal, from both the robot encoders (125 Hz) and from IMUs (286 Hz). Two different movement speeds were selected for the assessment, respectively the 50% (TR_50) and the 100% (TR_100) of the robot maximum speed (i.e. maximum 120 deg/s for the shoulder and 180 deg/s for the elbow) in order to test robustness of joint angle estimation to various angular velocities. One trial was recorded for each speed.
A model of the UR-10 was implemented in OpenSim ( Figure 2c) porting the URDF model available as part of the ROS-Industrial package 2 .
Virtual orientation sensors were placed on model links by aligning them to known reference points, as done during the preparation of the experimental setup. After the model calibration procedure, the joint angles were computed by the developed OB-IK tool using as input the orientations provided by the real IMUs attached to the robot links.
Results obtained from the OB-IK were compared to the experimental joint angles measured from the robot encoders (Figure 3), in terms of squared Pearson correlation coefficient (r 2 ), root mean square error (RMSE) and maximum absolute error (MAE) over the full trial. Within this framework, designed as the most controlled scenario, our aim was to validate OB-IK estimates of joint angles against data measured from robot encoders, considered as gold standard for accuracy.

Framework 2: validation against markerbased kinematics
For this framework, a rigid mechanism (phantom) consisting of four links (lengths from 110 to 150 mm) connected by three co-planar hinge joints was designed and 3D printed using plastic material (Figure 2b). Four Cometa WaveTrack IMUs (Cometa Figure 2. Experimental setups. Picture of the UR-10 robot with real IMUs placed (a). Picture of the custom-designed mechanism with both IMUs and passive markers placed (b). OpenSim model of the UR-10 robot (c) and of the custom-designed mechanism (d). In (c) and (d) virtual orientation sensors are placed on the models and numbered respectively A1 to A4 and B1 to B4. Joint names are also identified in both the models. UR-10 links' length: Systems, Italy) and fourteen passive reflective markers were positioned on the mechanism links. Marker trajectories were collected using a Vicon T160 with 10 cameras (Vicon Motion Systems Ltd., UK). A common trigger signal was used to synchronize the acquisition systems. Markers data were collected at 100 Hz, IMUs data at 286 Hz.
Three different trials were recorded involving respectively one (TR_1, j-1 joint involved), two (TR_2, only j-3 joint locked) and three (TR_3) degrees of freedom of the mechanism at the same time. During all trials the phantom was manually moved on a planar surface. At the beginning of each trial the mechanism was aligned to a reference to guarantee consistency of the starting position.
A model of the phantom mechanism (Figure 2c) was developed in OpenSim matching the CAD model used to design and print it. Virtual orientation sensors were placed on the model and their orientations, with respect to the segments, were refined using the described calibration procedure. Then, the measured IMU orientations were processed using the developed OB-IK tool. Marker trajectories were low-pass filtered with a 6 Hz, 4th order, zero-lag Butterworth filter. Marker data were then processed using the standard OpenSim (v.3.3) marker-based IK tool. Simulation quality for MB-IK was evaluated using tracking metrics such as RMSE and maximum tracking errors (reported as mean ± standard deviation).
Joint angle estimates from OB-IK were compared against MB-IK results (Figure 4) in terms of r 2 , RMSE, and maximum absolute error (MAE). The latter, computed for each frame, has been classified into three classes (i.e. lower than 6 deg, between 6 and 12 deg, and higher than 12 deg). Then the percentage of frames in each class has been computed. Finally, the same classification has been performed excluding from the trial the frames corresponding to joint accelerations and joint velocities higher than the 110% of the maximum reference values reported for human gait in Appendix 1 of Winter (2009). The final aim of this additional computation was to preliminary assess the performance of the proposed methodology in conditions comparable to the ones of the final targeted applications, i.e. human motion analysis.

Robustness of joint angle estimation to noisy input data
Framework 2 was then used to assess the robustness of the OB-IK to the experimental noise. As starting point for this analysis, marker trajectories and IMUs data from TR_3 were used. Using custom Matlab (v2016-b, The MathWorks, USA) code, Gaussian  noise was added to each component (i.e. X, Y, and Z) of each marker 3 D trajectory. The characteristics of the noise distribution were chosen to approximate realistic experimental noise (mean ¼ 0 and S.D. ¼ 3 mm) during a standard data collection (Di Marco et al. 2017). This procedure was repeated 20 times with different seeds of the random noise generator, obtaining 20 noisy versions of the original trial.
A similar procedure was used to generate noisy IMUs data. Since IMUs data were stored in quaternion form, and the quaternion space is not linear, it was not possible to directly add noise to each component. Therefore, three independent noise signals were generated, one for each IMU axis, and treated as if they were Euler angles defining a "noise" rotation in space, so that they could be converted into quaternion form and finally multiplied to the experimental quaternions. This procedure has the physical meaning of applying "noise" to the original orientation expressed in quaternion form. The amplitude of the Gaussian distribution (mean: 0 deg, std: 2 deg) was chosen equal to the worst orientation error declared by IMU manufacturers in a dynamic scenario. Same as for the markers case, 20 noisy trials were generated.
The obtained noisy data were then processed according to the procedure described in Framework 2. Similarly, the obtained results were then compared against the original data using the same metrics, i.e., r 2 , RMSE, and MAE. This procedure enabled to compare the robustness of our methodology to both MB-IK and to the forward use of IMU orientations in estimating joint angles.

Validation against encoder measurements
The goal of this framework was to assess the performances of the developed OB-IK tool with respect to directly measured robot joint angles. Obtained results were similar for both recorded trials (Table 1).
For TR_50 a correlation coefficient r 2 >0:999 was obtained for all joints. The highest RMSE ¼ 0.83 deg was obtained for the shoulder_pan_joint (see Figure 2). At the same joint was recorded also the highest MAE ¼ 1.76 deg. The higher amplitude of errors at the shoulder_pan_joint with respect to the other joints could be explained by the wider range of motion that it spanned during the task.
During the second trial ( Figure 5), the robot was moving at its maximum speed. It can be noticed that the maximum error amplitude increased up to around 6 deg at the extreme position of the range of motion, when higher linear accelerations occur on adjacent links. These results can be explained by the fact that IMU orientation is not directly measured but estimated using sensor fusion techniques, which are deeply affected by filter settings. In fact, the set of parameters chosen by the IMU manufacturer, appropriate for slow movements, were not suitable for high speed ones. The slow filter behavior explains the overshoot effects in Figure 5. In fact, it was realized during the data collection that the set of parameters chosen by the IMU manufacturer, neither accessible nor modifiable through the data collection system, were appropriate for slow movement but not for high speed ones.

Validation against marker-based kinematics
Joint angle estimates from OB-IK were compared with the MB-IK results ( Figure 6). Quantitative evaluation parameters, for all the framework trials, are presented in Table 2.
For joints actively involved in trial motion a good agreement between the two estimation methods was found in terms of both RMSE (< 5.8 deg) and r 2 (> 0.98) for all the trials, with at least the 66% of the complete trial characterized by absolute errors lower than 6 deg. Only for the most distal joint we found 31.4% of frames leading to error up to 12 deg, and in less than 2.5% of frames the tracking error was larger than 12 deg with a MAE of 17 deg in the worst frame. Values related to the unmobilized joints during each trial were omitted from Table 2, which only reports metrics for mobilized joints.
In evaluating the outcomes of this framework, it is worth to remember that both methods are affected by issues that could negatively influence their outputs. IMUs are sensitive to environmental noises and their dynamic behavior strongly depends on their internal filter settings, as shown in framework 1. On the other hand, joint angle estimates from MB-IK are sensitive to experimental marker placement and segment size. In this specific case, however, the metrics from MB-IK (RMSE <1:760:5 mm, maximum tracking error <3:561:0 mm in all the trials) allow us to consider the joint angle estimation of acceptable quality. The high amplitude of MAE obtained in this framework, with spikes up to 21 deg, could be due to the effect of the mechanism's size on the two IK algorithms. For the OB-IK if two IMUs are too close to each other a cross-talk effect, generated by magnetometers, could emerge leading to inaccurate orientation estimates. For MB-IK instead the smaller the body dimensions and the distance between joints and tracked markers, the larger will be the angular offset generated by the same marker tracking error. Furthermore, consequently to the general Figure 5. Orientation-based IK (dashed black) and encoder measured joint angles (cyan) during the trial TR_100 performed at maximum robot speed. Figure 6. Left column: orientation-based IK (dashed black) and marker-based IK (cyan) joint angle estimates for TR_02 when joint j-2 and j-3 were moving contemporaneously and joint j-1 was manually kept steady. Right column: percentages of frames for each absolute error (AE) class over all the frames of the trial (blue) and over the frames corresponding to frames corresponding to joint accelerations and velocities compatible with the ones which characterize the human gait according to literature data (Winter 2009). characteristics and assumptions of the algorithms employed inside each IMU to estimate its orientation, the closer the linear acceleration due to the motion is to the gravitational acceleration, the more prone to errors is the estimated orientation. In this case, for both TR_2 and TR_3, the linear accelerations of the segments reach values close to 1.5 g. However, such values are unlikely to be reached in biomechanical applications. Excluding from the trial the frames characterized by joint accelerations and velocities higher than the maximum values reported for a typical human gait (Winter 2009), it was found that at least the 83% of each trial has errors lower than 6 deg with a maximum of 16% of frames for which the AE is in between 6 and 12 deg. Concluding, the global metrics were reasonably aligned with the results of Framework 1 and could be considered a promising starting point to investigate the performances of OB-IK in estimating human kinematics. In that case indeed body segments are larger, therefore size-effects In TR_2 joint j-2 and j-3 were moved simultaneously and joint j-1 was kept steady. During TR_3 all the three joints of the mechanism were moved at the same time. Root mean squared error (RMSE), maximum absolute error (MAE) and correlation coefficient r 2 are reported. Moreover, percentages of frames included in each absolute error class are listed (in bold only the trial frames corresponding to joint accelerations and velocities compatible with the ones which characterize the human gait, according to literature data, Winter 2009, have been considered). should be less pronounced, and accelerations should be lower, leading to more accurate IMU orientation estimates.

Noise robustness analysis
The joint angle estimates obtained with the developed OB-IK tool for the 20 noisy trials were summarized in terms of mean and standard deviation (S.D.). Comparing that mean with the results obtained for the original data, no relevant differences arose (RMSE < 0.13 deg and a MAE < 0.67 deg were obtained for all joints). The mean standard deviation of the 20 noisy estimates was lower then 0.53 ± 0.1 deg for all joints.
In the worst case scenario of forward use of IMUs orientation to estimate joint angles, for a planar movement, a maximum joint angle error of 4 deg is expected, resulting from an orientation error (of opposite sign) for both joined links of 2 deg. This value is approximately five times greater than the one obtained with OB-IK, suggesting that including model constraint within the optimization framework results in a more accurate estimation of joint kinematics. Figure 7 reports the standard deviation of the 20 noisy trials outputs for both the markers case (cyan lines) and IMUs case (dashed black lines). The amplitude of the standard deviation of OB-IK outputs over the 20 trials is constant over joints, time, and movement speed. The amplitude of the standard deviation of MB-IK outputs is instead larger for the joint connecting the smallest links (j_3) and increases at higher movement speeds (as during the second half of the trial). This suggests that OB-IK produces more consistent estimation of angles in the face of varying dynamic conditions than the MB-IK.

Conclusions
The OB-IK methodology presented in this paper combines the benefits of model-driven simulations with those of inertial sensors in the challenge of accurately estimating joint kinematics with wearable technologies. The proposed methodology was tested both against joint angles directly measured from encoders and against joint angle estimates from commonly adopted multibody optimization procedure. Robustness to noise was also evaluated.
The first testing framework demonstrated the capabilities of the developed OB-IK in estimating joint angles with a very good accuracy when compared with experimentally measured angles. The important aspect of IMU internal filter settings emerged clearly.
Even if it is an aspect related to the acquisition devices and not to the developed tool, it should be taken into account by users. Indeed, to accurately estimate joint kinematics it is necessary to select the adequate parameter set for the filter (Mazz a et al. 2012), which in this investigation was not varied.
In the second framework, the OB-IK tool was validated using synchronous measurements from a stereophotogrammetric system combined with OpenSim's MB-IK algorithm. Obtained results are promising even if characterized by larger errors than those found in the first framework. The main source of errors was the most distal segment, which was also the smallest and most mobile in the tested movements. For this segment, positioning errors of markers and IMUs, together with experimental and methodological errors like its dimensions and its motion characteristics (highest accelerations and velocities), could have affected the estimated kinematics more heavily than other joints. Nevertheless, once considering only the frames characterized by dynamic parameters comparable with the ones of a typical biomechanical application (i.e. clinical gait analysis), the global metrics could be considered fairly good. Furthermore, they are comparable with the results obtained for the first framework.
Outcomes from the noise robustness analysis suggested that the use of a model-based inverse kinematics approach could reduce the effects of experimental noises and IMU non-idealities on the final joint angle estimates.
The presented assessments constitute a necessary step before moving to the application of the developed orientation-based Inverse Kinematics tool to human motion analysis. The calibration procedure developed to adjust the placement of virtual orientation sensors in the model can be easily applied for other mechanisms similar to those investigated in this work potentially obtaining fairly good results. However, before application to the study of human motion, the calibration methodology will have to be thoroughly assessed and might require specific modifications.
Finally, the developed methodology has been implemented as a plug-in for OpenSim and made freely available via SimTK project's page and Github repository 3 . The underlying research materials for this article can be accessed at https:/doi.org/10.15131/shef.data.7097744.