Validity, Test-Retest Reliability and Long-Term Stability of Magnetometer Free Inertial Sensor Based 3D Joint Kinematics

The present study investigates an algorithm for the calculation of 3D joint angles based on inertial measurement units (IMUs), omitting magnetometer data. Validity, test-retest reliability, and long-term stability are evaluated in reference to an optical motion capture (OMC) system. Twenty-eight healthy subjects performed a 6 min walk test. Three-dimensional joint kinematics of the lower extremity was recorded simultaneously by means of seven IMUs and an OptiTrack OMC system. To evaluate the performance, the root mean squared error (RMSE), mean range of motion error (ROME), coefficient of multiple correlations (CMC), Bland-Altman (BA) analysis, and intraclass correlation coefficient (ICC) were calculated. For all joints, the RMSE was lower than 2.40°, and the ROME was lower than 1.60°. The CMC revealed good to excellent waveform similarity. Reliability was moderate to excellent with ICC values of 0.52–0.99 for all joints. Error measures did not increase over time. When considering soft tissue artefacts, RMSE and ROME increased by an average of 2.2° ± 1.5° and 2.9° ± 1.7°. This study revealed an excellent correspondence of a magnetometer-free IMU system with an OMC system when excluding soft tissue artefacts.


Introduction
Marker-based optical motion capture (OMC) systems are commonly used in clinical movement analysis [1] and are therefore considered the gold standard. However, despite high resolutions and sub-millimeter accuracy, the application of OMC is expensive, time-consuming, and restricted to a laboratory environment. Therefore, body-worn inertial measurement units (IMUs) present a mobile alternative [1]. IMUs incorporate 3D accelerometers, 3D gyroscopes, and, typically, 3D magnetometers, measuring 3D linear acceleration, 3D angular velocity, and 3D magnetic field, respectively. Using sensor fusion algorithms, e.g., variations of the Kalman filter or optimization based methods [2], it is possible to estimate the IMUs' orientation in reference to a global coordinate system [3]. Combining more IMUs attached to linked body segments, it is possible to estimate the joint kinematics of the specified segments [1,4,5].
There are drawbacks concerning IMU systems that have to be addressed when measuring human motion. First, IMU-based orientation estimation suffers from drift due to the integration of noisy gyroscope measurements [6]. This is particularly challenging when omitting magnetometer data, which provide a global heading reference and can therefore be used to compensate for drift in the transversal plane [2].

Subjects and Data Acquisition
Twenty-eight healthy subjects (15 females, 13 males; age 24 ± 2.70 years; 70 ± 12.70 kg and 1.76 ± 0.09 m in height) participated in the study. Each of the subjects performed two test sessions on two days (6.75 ± 2.26 days in between). A test session consisted of one static neutral zero position (n-pose) sequence and a 6 min walk test [24]. The study was approved by the ethical committee of the Technische Universität Kaiserslautern (TUK) and meets the criteria of the declaration of Helsinki. After receiving all relevant study information, the participants signed an informed consent to the study including a permission to publish data.
IMUs were activated at least 20 min before measurement start. A static trial was performed before each subject was instrumented, with the sensors lying still for a period of at least 10 s, to estimate and subtract the gyroscope bias. These steps were conducted in accordance with the recommendations of Bergamini et al. [21].
Thirty-two retroreflective markers were attached to anatomical landmarks (AL) according to Leardini et al. [25] and the OptiTrack recommendations. Each IMU was secured in matched 3D printed boxes to which four markers were rigidly attached. These markers were used for unique identification in the optical point cloud as well as for orientation estimation. Using the OptiTrack Software, the origin of the boxes was moved to the center of the attached sensor casing. These box/sensor compounds were fixed to the body segments using straps and double-sided adhesive tape. IMUs were attached on the right and left dorsum of the foot approximately atop the base of metatarsal II-IV, on the right and left lateral aspect of the shank, due to better visibility, on the right and left lateral aspect of the lower third of the thigh and between the Spinae Iliacae Posteriores Superiores approximately atop the sacral base ( Figure 1). Inertial and optical data were simultaneously recorded at 60 Hz using XSens MVN Biomech (Version 4.3.7) and OptiTrack Motive (Version 1. 10.0) which were hardware synchronized using a standard 5 V TTL signal. The alignment orientations between the IMUs and the rigid boxes were calculated using the method described in [7]. The biomechanical model according to Cappozzo et al. [26] and the IMU-to-segment calibrations were extracted from the OMC data of the n-pose sequence. The joint centers were also calculated from the OMC data during the n-pose sequence according to the definitions of Visual3D (C-Motion, Inc, Germantown, MD, USA), a widely used software tool for 3D biomechanics research. The first OMC frame of each walking sequence was used as initialization for the IMU-based kinematics estimation. Both systems used the same biomechanical model.
The inertial data was processed with an iterated extended Kalman filter (IEKF) approach based on [14] while omitting magnetometer information. The gyroscope biases were extracted from a static sequence (see above), while the accelerometer biases were estimated in the IEKF along with the kinematics estimation using the model described in [11,27]. The same sequence was processed twice: initially to obtain a converged estimate of the acceleration bias, which was then used as initial guess in the second run. The estimated segment orientations were used to derive relative joint orientations. These were decomposed into joint angles using Euler angle decomposition [28]. The sensor fusion method is detailed in Appendix A.
To minimize STA, the OMC-based joint angles were derived from marker clusters on the rigid boxes (condition 1). For secondary analyses the joint angles were calculated based on the markers attached to the anatomical landmarks (condition 2). Initial contact (IC) was detected based on the left and right heel marker [29]. Turning phases in the 6 min walk tests were omitted. In order to investigate drift behavior, 10 left and right steps (one trial) were identified at three sections, i.e., beginning (A), middle (B), and end (C) of the test. All joint angle curves were normalized to 100 percent gait cycle (GC).

Statistical Analysis
To evaluate the IMU system, the root mean squared error (RMSE) and range of motion error (ROME), as well as 95% confidence interval (CI) were calculated for hip, knee, ankle joint, and the global pelvis orientation per section per GC. Further, Bland-Altman analysis (BA) was conducted to evaluate the limits of agreement (LoA) between the mean joint angle waveforms over all 28 subjects for both systems, considering only the normalized GC of section A. The results of the BA analysis are presented in the form 0.0°-0.0° ± 0.0°-0.0°. The first two numbers indicate the minimum and maximum of the mean differences between the systems. The last two numbers indicate the Inertial and optical data were simultaneously recorded at 60 Hz using XSens MVN Biomech (Version 4.3.7) and OptiTrack Motive (Version 1.10.0) which were hardware synchronized using a standard 5 V TTL signal. The alignment orientations between the IMUs and the rigid boxes were calculated using the method described in [7]. The biomechanical model according to Cappozzo et al. [26] and the IMU-to-segment calibrations were extracted from the OMC data of the n-pose sequence. The joint centers were also calculated from the OMC data during the n-pose sequence according to the definitions of Visual3D (C-Motion, Inc, Germantown, MD, USA), a widely used software tool for 3D biomechanics research. The first OMC frame of each walking sequence was used as initialization for the IMU-based kinematics estimation. Both systems used the same biomechanical model.
The inertial data was processed with an iterated extended Kalman filter (IEKF) approach based on [14] while omitting magnetometer information. The gyroscope biases were extracted from a static sequence (see above), while the accelerometer biases were estimated in the IEKF along with the kinematics estimation using the model described in [11,27]. The same sequence was processed twice: initially to obtain a converged estimate of the acceleration bias, which was then used as initial guess in the second run. The estimated segment orientations were used to derive relative joint orientations. These were decomposed into joint angles using Euler angle decomposition [28]. The sensor fusion method is detailed in Appendix A.
To minimize STA, the OMC-based joint angles were derived from marker clusters on the rigid boxes (condition 1). For secondary analyses the joint angles were calculated based on the markers attached to the anatomical landmarks (condition 2). Initial contact (IC) was detected based on the left and right heel marker [29]. Turning phases in the 6 min walk tests were omitted. In order to investigate drift behavior, 10 left and right steps (one trial) were identified at three sections, i.e., beginning (A), middle (B), and end (C) of the test. All joint angle curves were normalized to 100 percent gait cycle (GC).

Statistical Analysis
To evaluate the IMU system, the root mean squared error (RMSE) and range of motion error (ROME), as well as 95% confidence interval (CI) were calculated for hip, knee, ankle joint, and the global pelvis orientation per section per GC. Further, Bland-Altman analysis (BA) was conducted to evaluate the limits of agreement (LoA) between the mean joint angle waveforms over all 28 subjects for both systems, considering only the normalized GC of section A. The results of the BA analysis are presented in the form 0.0 • -0.0 • ± 0.0 • -0.0 • . The first two numbers indicate the minimum and maximum of the mean differences between the systems. The last two numbers indicate the minimum and maximum of the limits of agreement (95% CI) of the two systems. The coefficient of multiple correlation (CMC) was calculated for each parameter per section per GC according to Ferrari et al. [30]. In [30] they showed that if a joint angle waveform reveals a similar ROM compared to the overall offset between the two signals, the CMC can become a complex number. If that happened in the current calculations for individual subjects, these results were not considered for further analysis. All these calculations were conducted for both, condition 1 and condition 2.
A paired t-test was performed to identify significant differences in the RMSE and ROME of the joint angles of all sections between condition 1 and condition 2. Alpha level was set a priori to 0.05. The Chi-square goodness-of-fit test was carried out to check for normal distribution in the data.
For the evaluation of the test-retest reliability, the intraclass correlation coefficient (ICC) for inter-day reliability was calculated for test day one and test day two for both systems for every joint and section according to McGraw and Wong [31]. CMC and ICC values were rated according to Koo and Li [32].
The global heading direction error of the pelvis in the transversal plane (pelvis rotation error) was examined at minute 0, 1, 2, 3, 4, 5, and 6. In addition, the RMSE and the ROME of the joint angles and pelvis flexion/obliquity were examined regarding potential linear trends over time. Therefore, lines were fitted via linear least squares regression to the RMSE and ROME values of each GC for the abovementioned angles and all test persons (Matlab function "fit" with "fittype = poly1"). The slopes of the fitted lines were computed and plotted to evaluate potential trends. Processing of the joint angles and statistics were conducted in Matlab 2015 (Mathworks Inc., Natick, MA, USA).

Condition 1-Marker Clusters
RMSE and ROME for all parameters over all sections are shown in Table 1. RMSE and ROME between the two systems revealed mean values lower than 2.40 • and lower than 1.60 • respectively in all joints. The poorest outcome concerning the RMSE was evident in knee rotation (1.75 • -2.34 • ) and knee abduction for ROME (1.11 • -1.58 • ).  Table 1. Mean root mean squared error (RMSE) and mean range of motion error (ROME) of condition 1 over all subjects ± standard deviation (SD); brackets contain 95% confidence interval (CI). A, B, C indicate beginning, middle, and end of the 6 min walk test. RMSE (deg) ± SD (95% CI) ROME (deg) ± SD (95% CI)    The CMC showed very high waveform similarity for the joint angles in the sagittal plane with mean values ranging from 0.99 to 1. Concerning the frontal and transversal plane, the CMC showed slightly lower correspondence with mean values ranging from 0.88 to 0.99. All CMC values are mapped in the box-and-whisker-plot in Figure 4a.

Condition 2-Skin Markers
Most error measures showed results inferior to condition 1, when considering STA. RMSE and ROME for all joint angles over all sections are shown in Table 2. RMSE showed values lower than 6.00° over all planes and joints. ROME showed values lower than 6.10° for hip, knee and pelvis in all planes, and ankle joint in transversal and frontal plane. However, the movement in the sagittal plane in the ankle joint revealed a ROME of up to 10.66°. Figure 2b exemplary shows the LT ankle flexion of a representative subject for condition 2. Table 2. Mean RMSE and ROME of condition 2 over all subjects ± SD; brackets contain 95% CI. A, B, C indicate beginning, middle and end of the 6 min walk test.  The CMC showed very high waveform similarity for the joint angles in the sagittal plane with mean values ranging from 0.99 to 1. Concerning the frontal and transversal plane, the CMC showed slightly lower correspondence with mean values ranging from 0.88 to 0.99. All CMC values are mapped in the box-and-whisker-plot in Figure 4a.

Condition 2-Skin Markers
Most error measures showed results inferior to condition 1, when considering STA. RMSE and ROME for all joint angles over all sections are shown in Table 2. RMSE showed values lower than 6.00 • over all planes and joints. ROME showed values lower than 6.10 • for hip, knee and pelvis in all planes, and ankle joint in transversal and frontal plane. However, the movement in the sagittal plane in the ankle joint revealed a ROME of up to 10.66 • . Figure 2b exemplary shows the LT ankle flexion of a representative subject for condition 2.  Concerning condition 2 BA analysis revealed mean differences inferior to condition 1. Sagittal joint angles showed biases of −4. The paired t-test revealed significant differences between condition 1 and condition 2 of the RMSE in all joint angles and sections. ROME showed significant differences between condition 1 and 2 in all joint angles and sections excepting the right hip flexion of section A (p = 0.08). For detailed results see Table 3.

Test-Retest Reliability
The ICC revealed moderate to excellent correlations over all joints in the frontal and transversal plane (0.52-0.93) and excellent values in the sagittal plane (0.94-0.99) ( Table 4). The knee joint in the frontal and transversal plane, the pelvis in the frontal plane, and the hip joint in the transversal plane

Test-Retest Reliability
The ICC revealed moderate to excellent correlations over all joints in the frontal and transversal plane (0.52-0.93) and excellent values in the sagittal plane (0.94-0.99) ( Table 4). The knee joint in the frontal and transversal plane, the pelvis in the frontal plane, and the hip joint in the transversal plane showed the lowest values (0.52-0.76). This tendency consented with the results from ICC calculation of the optical system (0.63-0.83). However, overall ICC values of the OMC system (0.63-0.99) were higher compared to the IMU system.

Drift
The global heading direction of the pelvis in the transversal plane (pelvis rotation) drifted linearly but not consistently between subjects (45 • ± 58 • ). Global heading errors after six minutes ranged from 0.21 • up to 230 • ( Figure 5).
To analyze the drift in the joint angle data, changes in the RMSE and ROME of condition 2 over time were evaluated through linear least squares regression (line fitting). The slopes of the fitted lines of RMSE and ROME for all joint angles of the left side and global pelvis flexion and obliquity of all test persons are shown in Figures 6 and 7. Consistent positive slopes (i.e., increasing RMSE/ROME error values over time) would indicate a systematic drift over time. However, as visible in the figures, the slopes of both RMSE and ROME over time reside above as well as below zero, so that there is no clear trend over all test persons visible. Moreover, the slopes are in a range where the errors cannot be distinguished from noise given disturbing effects such as STA. visible in the figures, the slopes of both RMSE and ROME over time reside above as well as below zero, so that there is no clear trend over all test persons visible. Moreover, the slopes are in a range where the errors cannot be distinguished from noise given disturbing effects such as STA.   visible in the figures, the slopes of both RMSE and ROME over time reside above as well as below zero, so that there is no clear trend over all test persons visible. Moreover, the slopes are in a range where the errors cannot be distinguished from noise given disturbing effects such as STA.

Discussion
This paper evaluated the performance of a sensor fusion algorithm for estimating 3D kinematics from gyroscope and accelerometer data.
The BA analysis revealed excellent agreement between the systems. Mean difference results were similar to the results of Robert-Lachaine et al. [17]. They consent with the present findings of the sagittal joint angles showing the best agreement and the transversal joint angles the poorest agreement between the two systems. However, the limits found by Robert-Lachaine et al. [17] ranged from 2.8°-7.0°, compared to limits of 0.3°-1.3° in the present findings.
CMC values showed excellent correspondence between the systems in the sagittal plane and good to excellent correspondence in the frontal and transversal plane. These findings are in accordance with Ferrari et al. [19] and outperform the results of Zhang et al. [18]. However, the former performed an offset correction to increase waveform similarity and did not consider transversal and frontal plane of the knee joint. Zhang et al. used different coordinate systems for joint angle calculation [18]. The findings of Robert-Lachaine et al. [17] showed better CMC values concerning knee rotation and abduction (0.91 and 0.97). For ankle joint rotation and abduction, they report their lowest CMC values (0.89 and 0.77). In this case, the current system achieved higher correspondence (0.96-0.98). However, Robert-Lachaine et al. [17] do not report whether these results are mean values over the entire time period of 32 min. Furthermore, Robert-Lachaine et al. [17] analyzed a complex combination of movements rather than a standardized motion.
CMC values were lowest for knee abduction (0.88-0.93). The majority of literature reports lowest CMC outcomes considering motion in the transversal and frontal planes [17][18][19]. However, during walking, movements in these planes show smaller ranges of motion compared to the sagittal plane. It has been shown that the CMC decreases as the amplitude of movement decreases [33]. This

Discussion
This paper evaluated the performance of a sensor fusion algorithm for estimating 3D kinematics from gyroscope and accelerometer data.
The BA analysis revealed excellent agreement between the systems. Mean difference results were similar to the results of Robert-Lachaine et al. [17]. They consent with the present findings of the sagittal joint angles showing the best agreement and the transversal joint angles the poorest agreement between the two systems. However, the limits found by Robert-Lachaine et al. [17]  CMC values showed excellent correspondence between the systems in the sagittal plane and good to excellent correspondence in the frontal and transversal plane. These findings are in accordance with Ferrari et al. [19] and outperform the results of Zhang et al. [18]. However, the former performed an offset correction to increase waveform similarity and did not consider transversal and frontal plane of the knee joint. Zhang et al. used different coordinate systems for joint angle calculation [18]. The findings of Robert-Lachaine et al. [17] showed better CMC values concerning knee rotation and abduction (0.91 and 0.97). For ankle joint rotation and abduction, they report their lowest CMC values (0.89 and 0.77). In this case, the current system achieved higher correspondence (0.96-0.98). However, Robert-Lachaine et al. [17] do not report whether these results are mean values over the entire time period of 32 min. Furthermore, Robert-Lachaine et al. [17] analyzed a complex combination of movements rather than a standardized motion.
CMC values were lowest for knee abduction (0.88-0.93). The majority of literature reports lowest CMC outcomes considering motion in the transversal and frontal planes [17][18][19]. However, during walking, movements in these planes show smaller ranges of motion compared to the sagittal plane. It has been shown that the CMC decreases as the amplitude of movement decreases [33]. This might also explain the rather good results of Robert-Lachaine et al. concerning the CMC of the knee in frontal and transversal plane (0.91-0.97). In this study, subjects had to perform lifting and turning tasks which could have led to increased ROM in the mentioned planes.

Condition 2-Skin Markers
As mentioned, most values of errors increased. This consented with the findings of Seel et al. [4] who compared knee and ankle sagittal joint kinematics of the human leg and the transfemoral prosthesis of an above-knee amputee. They found the errors on the human leg four times higher than on the prosthesis. On average, RMSE and ROME increased by 2.20 • ± 1.50 • and 2.90 • ± 1.70 • , respectively. Nüesch et al. [23] reported RMSE of the sagittal plane for hip, knee, and ankle joint (9.60 • , 7.60 • , 4.50 • ). The present system revealed better results concerning the sagittal plane for hip and knee (3.83 • , 2.66 • ) and similar results for the ankle (5.48 • ). However, Nüesch et al. [23] conducted their examination on a treadmill, which could have affected IMU-derived data [14]. Fasel et al. [12] found in their evaluation of skiing similar hip abduction ROME (−3.3 • ± 4.1 • ), slightly higher knee abduction ROME (4.2 • ± 5.5 • ) and distinctively higher hip flexion ROME (−10.7 • ± 4.3 • ). Their findings concerning knee flexion, hip and knee rotation showed better results for ROME around −0.1 • and 0.5 • . However, note that Fasel et al. [12] investigated a different task, which limits comparability. Further, they conducted their examination on an indoor skiing carpet, comparable to a treadmill, therefore inviting the same considerations associated with Nüesch et al. [23]. Interestingly, ROME of the ankle joint flexion increased from 1.61 • in condition 1 to 10.66 • in condition 2. In the IMU system and when calculating joint angles based on marker clusters, the foot is assumed to be one rigid segment. However, the foot is a complex organization of bony segments [34]. Skin markers attached to anatomical landmarks on the foot are in fact placed on different segments rather than on one segment only. This might explain the differences in the ankle flexion between IMU and OMC system. BA analysis also showed increased biases compared to condition 1. The ankle flexion was the most affected joint angle in the sagittal plane, the knee abduction the most affected joint angle in the frontal plane and the hip rotation the most affected joint angle in the transversal plane ( Figure 3). CMC values decreased mostly for the transversal and frontal plane. This might be due to the increased uncertainty concerning CMC and smaller ranges of motion [33]. However, the ankle joint flexion still presents good to excellent values. Qualitative examination of the ankle flexion waveforms showed that there is an excellent match between waveforms of the IMU system and the OMC system at about 10 to 50% GC. However, at 60 to 100% GC an offset appeared (Figure 2b). This might explain the difference of ROM, but still similar shape of waveform and good CMC values. Al-Amri et al. [20] showed a similar shape of waveform for ankle flexion in their examination. No such differences were found between the IMU system and the OMC system based on the marker clusters (Figure 2a). Al-Amri et al. [20] examined in their recent study the validity of a commercial IMU system during one 8 m walk, taking into account errors due to STA. They found excellent CMC values concerning the sagittal plane of all three joints and the frontal plane of the hip joint but stated a poor outcome for the remaining joints in the frontal and transversal plane. However, they did not report exact values of the CMC for the frontal and transversal plane because results were complex numbers at times. In this study, CMC values for some subjects also resulted in complex numbers. On average, 3 out of 28 subjects per joint angle resulted in complex numbers. However, these subjects were ignored for the calculation of the mean CMC values shown in Figure 4. Nevertheless, concerning the findings of Al-Amri et al. [20], R 2 values indicated rather poor correlations for the joint angles in the transversal plane. That consents with the present results, showing moderate CMC values in the hip joint angle of the transversal plane (Figure 4b). Al-Amri et al. [20] used different biomechanical models for their analysis and stated possible uncertainties in the optical data due to the marker protocol. Further, Fiorentino et al. [35] showed that the 3D hip joint angles and ROM based on OMC systems are significantly influenced by STA.

Test-Retest Reliability
ICC values revealed moderate to excellent results for all joint angles. Knee abduction (0.56-0.58) and pelvis obliquity (0.52) showed the lowest test-retest correlation. These results are in accordance with Mills et al. [36]. However, the ICC calculation for the OMC system based on skin markers showed overall better outcome. This fact might be explained by expert marker placement. The IMUs were attached only approximately to the same spot on the segments. Furthermore, the IMUs were more prone to STA than the skin markers due to the rigid boxes and positioning focused on better visibility. These circumstances may have caused different amounts of errors on the two measurement days. For both systems, values of measures of reproducibility were not high. A possible explanation is the fact that inter-day variability of gait is considered higher than, for example, intra-day variability [36].

Drift
One difference of the current IMU system compared to the systems used in the referenced studies is the omission of the magnetometer information. Favre et al. [10] measured the 3D knee angle omitting magnetometer information and presented results with a mean error of 4.00 • to 8.10 • . They also introduced a functional calibration method to align the joint coordinate system. Note that the present evaluation relates to the sensor fusion algorithm, while the calibration, as mentioned in Section 2.1, was obtained from the OMC system. Thus, the results presented in this study can be considered free of calibration errors.
However, omitting global heading direction information (obtained through undisturbed magnetometer measurements) typically leads to drift [2,12]. In the present study, gait was measured for six minutes. Robert-Lachaine et al. [17] analyzed ergonomic tasks over a period of 32 min. However, they do not state results for different sections of their test procedure and thus give no hint as to whether drift appeared in the kinematic data or not.
Fasel et al. [12] measured skiing for 120 s. They reported errors in their joint angle results (based on individual per-IMU orientation estimates) due to drift and introduced methods for drift reduction [13] based on adjacent segments as a second processing stage (cf. Section 1). On the contrary, the present study revealed no systematic drift over all test persons, neither in the 3D joint angle data nor in the global pelvis flexion/obliquity. This is due to including biomechanical constraints in terms of connected segments and environmental constraints in terms of ground contacts directly into the estimation. However, Fasel et al. [13] reported that mainly highly flexed joint angles were affected from the drift. In their examination of skiing subjects reached distinctively higher peaks in hip and knee flexion compared to the present study.
Bergamini et al. [21] examined the drift in the orientation obtained from two inertial sensors mounted on wrist and pelvis during 180 s of gait. In the transversal plane, a drift of up to 40 • was measured. In the sagittal plane and frontal plane, drift was smaller with values up to approximately 5 • . The findings of Bergamini et al. [21] are comparable with this study's results for drift in the global heading direction estimate (investigated at the root segment, i.e., the pelvic rotation) measured at values of up to 230 • . Note, inter-segment constraints do not provide corrective information concerning the global heading direction, which explains the linear drift observed. Consequently, the global pelvic rotation was neglected in the evaluation. However, the evaluation and interpretation of the ROM of the pelvic rotation and its reproducibility were independent of the drift.

Conclusions
The present algorithm for the calculation of 3D joint angles based on gyroscope and accelerometer data from seven IMUs mounted on the lower body shows good to excellent agreement when compared to a common OMC system and excluding STA. However, in this study, the influence of STA was shown. Especially the ankle joint was highly affected by these and further artefacts, e.g., a limited biomechanical model. Further research has to be conducted to better compensate these effects.
In terms of reliability, the results indicate that an OMC system combined with an experienced examiner delivers a better outcome, particularly for knee abduction and rotation and the ankle joint. Better placement of the shank sensor and smaller IMUs might improve overall reliability and sensitivity to STA.
Systematic drift was observed only in the global transversal plane angle (investigated at the root segment, i.e., the pelvis rotation). There was no systematic drift observed over all test persons in the other kinematic parameters. However, in clinical gait analysis the ROM per GC is the more essential criterion, which was also measured with satisfying accuracy for the pelvis rotation.
The current sensor fusion algorithm was not only shown to be comparable to other algorithms, but also tends to outperform most algorithms examined so far in terms of its accuracy while being magnetometer-independent. However, it has to be considered that this study focused on the evaluation of the sensor fusion algorithm, while the IMU-to-segment calibration, the biomechanical model and the initialization were obtained from the OMC system. Therefore, the next step consists of evaluating the validated sensor fusion method in a setup, where all information was obtained purely from the IMU system.
Nevertheless, this examination reveals promising results of a magnetometer-independent sensor fusion algorithm that showed no systematic drift in the joint angle data. Therefore, a stand-alone system incorporating this algorithm provides potential for applications in clinical gait analysis and further implementations.
Author Contributions: W.T. contributed to the development of the study conception, data acquisition, analysis and interpretation, manuscript drafting and revision. M.M. contributed to data acquisition, data processing, manuscript drafting. B.T. contributed to data analysis, manuscript drafting and revision. M.F. contributed to study conception, data interpretation, manuscript drafting and revision. G.B. contributed to study conception, data processing, data interpretation, manuscript drafting and revision.

Conflicts of Interest:
The authors declare no conflict of interest. The authors alone are responsible for the content and writing of this paper. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Appendix A
The here presented sensor fusion method for segment kinematics estimation is a slightly modified version of [14], which itself was based on [2]. Note that the contribution of [14] was the extension of the best performing extended Kalman filter (EKF)-based segment orientation estimation method of [2] with global translation estimation via ground contact estimation. For this study, the method of [14] was used while omitting magnetometer measurements (i.e., skipping measurement updates with the magnetometer data). Moreover, it was extended with online accelerometer bias estimation as proposed in [11,27] and by using an IEKF instead of an EKF [37]. For the reader's convenience, the resulting method is detailed in the following. For further information, please refer to the abovementioned publications.

Appendix A.1. Biomechanical Model
The biomechanical model is a stick-figure-like model with rigid segments. Each segment S i ∈ S has an associated segment coordinate frame localized in the global position-less coordinate system G. The latter has z pointing up (opposing gravity). The global six DOF pose of the i-th segment is represented by the unit quaternion q GS (denoting the orientation between G and S) and the position vector S G (denoting the position of frame S in G). Furthermore, each segment S i defines a set of points P S i which are represented in the segment coordinate system and are used for joint definitions: a joint J k ∈ J defines which points from which segments are connected. Hence, a joint is a quadruple (i, k, j, l), meaning that p k ∈ P S i is connected to p l ∈ P S j in the global frame. Additionally, each segment S i has an IMU I i attached, with the rigid transformation q SI , I S i . This transformation is referred to as IMU-to-segment (I2S) calibration. The full model of the lower body is depicted in Figure A1. Additionally, each segment has an IMU attached, with the rigid transformation { , } . This transformation is referred to as IMU-to-segment (I2S) calibration. The full model of the lower body is depicted in Figure A1. Figure A1. Biomechanical model of the lower body with segments (magenta lines), joint centers (red spheres), four ground contact points on each foot (green spheres), IMU placement, and involved coordinate frames. A technical coordinate system is associated to each IMU ( ). The segment coordinate systems ( ) are drawn at the proximal ends of the segments. The six DoF transformations, each in terms of an orientation (quaternion) and a translation , between the IMU coordinate frames and the associated segment coordinate frames are called IMU-to-segment calibrations. The symbol denotes the global coordinate system. The figure has been taken from [14].
is the estimated global point height, and is the estimated global point velocity. Here, S(·) denotes the skew-symmetric matrix of a vector. For information fusion, a binary discrete Bayes filter [38] is used for each potential contact point. Assuming the prediction being independent of the previous state, given a random variable X t = {x 1 := contact | x 2 := no contact} the Bayes filter simplifies to The probabilities are obtained as follows: Here, the parameters for the sigmoid function were empirically determined. For each ground contact point, these probabilities are calculated in every time step. If the largest ground contact probability exceeds a threshold p th , then two soft constraints (in terms of measurement updates) are applied to the corresponding point p k inside the kinematics estimation filter: a zero-velocity update p and δ z ∼ N (0, Σ z ) being Gaussian white measurement noises. If the threshold was not exceeded, a zero-plane update is applied to the lowest ground contact point only if this has a negative height. Otherwise no additional measurement update is conducted. Obviously, in the current state, the method assumes one plane. Table A1 contains all previously introduced noises and parameters. negative height. Otherwise no additional measurement update is conducted. Obviously, in the current state, the method assumes one plane. Table A1 contains all previously introduced noises and parameters.               Figure A6. BA diagrams of 3D joint rotations of hip, knee, and ankle of the right side for condition 2. The solid line indicates the mean difference. The dashed lines indicate the LoA (95% CI) of the mean difference. Figure A7. BA diagrams of 3D joint rotations of the pelvis for condition 2. The solid line indicates the mean difference. The dashed lines indicate the LoA (95% CI) of the mean difference.