Body-IMU autocalibration for inertial hip and knee joint tracking

Sensor to body calibration is a key requirement for capturing accurate body movements in applications based on wearable systems. In this paper, we consider the specific problem of estimating the positions of multiple inertial measurement units (IMUs) relative to the adjacent body joints. To derive an efficient, robust and precise method based on a practical procedure is a crucial as well as challenging task when developing a wearable system with multiple embedded IMUs. In this work, first, we perform a theoretical analysis of an existing position calibration method, showing its limited applicability for the hip and knee joint. Based on this, we propose a method for simultaneously estimating the positions of three IMUs (mounted on pelvis, upper leg, lower leg) relative to these joints. The latter are here considered as an ensemble. Finally, we perform an experimental evaluation based on simulated and real data, showing the improvements of our calibration method as well as lines of future work.


INTRODUCTION
In recent years the attention towards mobile sensing for applications in sport and fitness monitoring has been drastically increasing.This is due to several factors, such as fitness awareness, the raise of motivation through gamification and competition in social networking, the proliferation of smartphones, and most importantly the advances in wearable technology and smart textiles.These wearable systems give the opportunity of self reflection to not only professional athletes but also ordinary people who are willing to pursue a healthy life and stay motivated in their training routine.
In contrast to smart watches and gadgets which contain single sensing units and mostly provide general information, such as walking speed and distance, energy expenditure, and type of activity [1,2], smart textiles with multiple embedded sensing units can provide tracking of detailed body postures, e.g. in terms of joint angles and body segment positions [4,18].The latter is required during training in order to monitor and improve exercise technique and avoid injuries [5].Moreover, body pose estimation can be used for model-based physical activity recognition, as demonstrated in [16,23].Such precise information, however, requires the knowledge of how the embedded sensors are mounted on each body segment when the user is wearing the textile.Using assumptions and information from anthropometric tables could result in inaccurate tracking, since there are different types of target users with different body shapes.In addition there are various types of clothing material which are used in the design of the wearables and result in varying sensor placements even for one specific individual.Manual measurements of these parameters are also cumbersome and error-prone.Therefore, an accurate, robust and autonomous sensor to body autocalibration method is required for a wearable training suit.
In [8,14], two functional calibration procedures are proposed in order to obtain joint axes of the lower body based on fusion of the measurements of each IMU.The estimated joint axes are then used to track the joint angles in a reliable and clinically interpretable way.However, these procedures require passive and controlled movements, which are difficult to perform without supervision.Palermo et al. in [15] propose an easy and repeatable method without the need for specific and accurate movements in order to obtain the relative orientations of the IMUs with respect to the body.However, none of the above contributions provide a method for estimating the sensor positions with respect to the adjacent joints, which is a critical information in joint angle estimation using accelerometer measurements, especially during fast rotations [7] and when modeling kinematic chains [12].In order to compensate for the resulting errors, in [14] a rotational mapping based on the integration of gyroscope measurements is used.Thus, the joint angle track-ing accuracy is limited to time due to drift.In [10,11], the IMU positions are used during tracking, however, these parameters are obtained through manual measurements, which is not accurate.In [24], a Kalman filter based method for minimizing the positional errors is proposed.However, this method requires multiple IMUs on each segment with controlled positional errors with respect to each other.Moreover, the large state vector, which is established for this method, is not efficient considering the computational power available in mobile applications.
In [21], Seel et al. propose a novel and practical method for estimating the joint axis of a hinge joint, as well as, the IMU positions with respect to a spheroidal joint, based on the gyroscope and accelerometer measurements of two IMUs mounted on the adjacent segments.In [20], the applicability of these methods is also discussed when considering joints with different degrees of freedom (DOF), and simple extensions are presented.In [21,20], the method is applied to the knee and ankle joints, while in [13] it is adopted for a home-based clinical knee rehabilitation system.The methods of Seel et al. are drift-free, since they are based directly on the measurements, and work without the need for dedicated calibration movements.However, Seel et al. didn't provide an investigation of the types of movement that lead to the most accurate results.With insufficient variations in the different DOFs, the position calibration method leads to inaccurate results.This became apparent when applying the method to the hip joint, where it is difficult to perform movements with sufficient variation.This paper builds upon the IMU position calibration method presented in [21,20] with the goal to improve its accuracy and robustness under suboptimal movement conditions and make it applicable to the IMU-body calibration for an existing smart training pants [18], where the hip and knee joints are in focus.These are important, e.g. for monitoring exercises, such as squats [19].We first provide a preliminary observability analysis of Seel's position calibration method (Section 3), and, based on this, propose a new calibration approach which, in contrast to Seel, considers three linked segments with IMUs (pelvis, upper leg, lower leg) and, respectively, two joints (hip, knee) in one estimation problem (Section 4).This makes it possible to benefit from an additional constraint, which is shown in experiments to provide more robust and accurate results under suboptimal movement conditions (Section 5).Note that this paper focuses on IMU position calibration, assuming the joint axes to be known, e.g. based on the method in [20].

PROBLEM FORMULATION
The problem of IMU position estimation can be defined as follows (cf. Figure 1) [20]: Assume two IMUs, A and B, being mounted on two segments connected via a joint m.By using sequences of measurements from A and B, the goal is to determine the two IMUs' positions, lmS, S ∈ {A, B}, relative to the joint m.More precisely, lmS denotes the vector from the joint center to the IMU center, given in the local coordinate frame of the IMU.Consider, for instance, the case where m is the knee joint.Then, IMU A is placed on the upper leg segment, IMU B is placed on the lower leg segment, and the goal is to derive the two vectors from the knee joint center to the IMUs in their respective coordinate frames.
For the above defined problem, an implicit, stochastic measurement model for timestep i can be formulated as: Here, zmi refers to the measurement vector at timestamp i ∈ 1 . . .k.This measurement vector includes 3D acceleration, aSi, angular velocity, ωSi, and angular acceleration, αSi, of both IMUs, S = {A, B}.Note that αSi is assumed to be derived from ωSi, e.g.via a third order approximation [20].The quantity xm refers to the parameter vector comprising the IMU positions.Finally, for simplicity, we assume additive Gaussian measurement noise emi ∼ N (0, Σ).
Given (1), a point estimate of xm can be obtained by maximizing the maximum likelihood function: Using the logarithm formulation and then eliminating all constant terms from the optimization, this results in a weighted least squares problem, which can be solved with standard techniques: Considering a spheroidal joint, Seel et al. introduces the following deterministic measurement model for estimating the IMU positions [21]: where and emi denotes the residual to be minimized.This model can be directly used in Equation ( 3), even if the latter reduces to a least squares formulation due to its deterministic nature.
The following section presents a preliminary observability analysis of the optimization problem (3) in combination with model ( 4), in which we identify major failure cases.

OBSERVABILITY ANALYSIS
The observability of the general optimization problem in (3) can be evaluated by computing the Fisher Information Matrix (FIM), which in our setting can be defined as [22]: Here, J h denotes the Jacobian of the measurement model ( 1) with respect to the parameter vector xm: If J h has full column rank, then F IM has full rank [22], which indicates observability.
Using the measurement model ( 4), the partial derivatives in ( 6) for timestep i are: where Here, in order to simplify, the timestep subscript is dropped and the following notations are introduced: By substituting ( 8) in (7), we obtain: where with S ∈ {A, B} and WSpq, p, q ∈ {1, 2, 3} refers to the components of WS.
Now, in order to analyze the observability of (3), one can examine the situations, where the components in (9) are dependent.Considering each IMU (the respective blocks of J h ) separately, two cases leading to rank deficiency are: 1.When there is no rotation in at least two DOFs.For instance, assume ωSx = ωSy = 0 for a sequence of measurements, and consequently αSx = αSy = 0, which results in CS3 = 0.Then, the third column in the respective block and (measurement) rows of J h equals zero, which reduces the rank.This is the case for the knee joint, since the dominant rotation is flexion/extension, i.e.only in one DOF.
2. When rotation appears with the same angular velocity in three DOFs for a sequence of measurements.In this situation ωSx = ωSy = ωSz in (10), which results in αSx = αSy = αSz.By substituting these equalities, the columns in ( 9) are related by CS1 + CS2 + CS3 = 0.
3. Additionally, when considering both IMUs, from (9), another case leading to rank deficiency is, when two IMUs rotate with the same angular velocities in all DOFs for a sequence of measurements.
In the case of the hip joint, due to the limited movements of the pelvis, it is in fact difficult to perform movements with sufficient but different angular velocities in all the DOFs of the pelvis.Additionally, as it is discussed e.g. in [6], pelvis rotations are in-phase with upper leg swings during relatively high velocities.Such movements are typical movements performed in the calibration phase, as it was visible in our experiments.Thus, it is difficult to avoid situations ( 2) and (3).Note that the movements on the two segments related to a joint should be simultaneous, since otherwise situation (1) occurs.
As indicated in Section 1, Seel et al. propose a correction of the estimated positions of two IMUs B and C in the degenerate case (1) of a hinge joint n (cf. Figure 2) [20].In this case, every point on the hinge joint axis rn is a solution of (4).The correction of the estimated vectors, lnS, S ∈ {B, C} corresponds to a shift of the joint center, n, on the known joint rotation axes, rnS, S ∈ B, C, as being represented in the coordinate frames of both IMUs.The shift computation is based on the assumption that the true joint center is the point on the joint axis that is closest to both IMUs.This can be formalized as: where S enCoR denotes the shift of the joint Center of Rotation (CoR) represented in the local coordinate frame of S. The latter is indicated by the left superscript.Moreover, rnS, S ∈ B, C are assumed known through calibration [20].The corrected IMU positions are then: The above correction method has been included in the experimental evaluation in Section 5.

PROPOSED METHOD
The previous section describes different cases of movements, which can lead to inaccurate results when applying the IMU position calibration method of [20] to the hip and knee joint, which is considered in our work.Therefore, this section proposes a novel method, which improves upon the existing one through the following extensions: The estimation of the positions x = [xm, xn] T of the three IMUs S ∈ {A, B, C} (on pelvis, upper leg, lower leg) with respect to the two joints m, n (hip, knee) (cf. Figure 2), from the IMU measurements z = [zm, zn] T is modeled jointly in one optimization problem analogously to (3), but with a new measurement model, which is introduced in the following.

Proposed measurement model
First, instead of considering only the lengths of the acceleration vectors in (4), the proposed, stochastic measurement model incorporates the relative orientations of the IMUs, leading to: where [emi, eni] T ∼ N (0, Σ) and the relative rotation RY X from X to Y can be obtained as RY X = R T GY RGX .Here, RSG, S ∈ {X, Y } is the orientation of IMU S with respect to a fixed global frame G, which is typically aligned with gravity and magnetic north.The global orientations can be estimated from the IMU measurements using a standard sensor fusion technique, e.g.[9].A respective implementation is typically provided with commercially available integrated IMU sensor packages, such as the ones integrated in the smart training pants [18].

Constraints of three connected segments
We consider the hip and knee jointly in one optimization problem.This allows us to extend the measurement models in (13) with additional constraints modeling the fact that these joints, m, n in Figure 2, are linked via the upper leg segment lmn.
We assume the flexion/extension joint axes of the hip and knee, rmy and rn, as being approximately coplanar.In Figure 3 (left), the geometry is illustrated for an ideal coplanar case with intersecting axes, where the latter span the plane P1 with normal vector nP 1.Furthermore, lmB and lnB are by definition coplanar, spanning the plane P2 with normal vector nP 2. Assuming that rmy and rn are not parallel, the normal vectors can be obtained using cross products: As is visible in the figure, in the ideal case, P1 and P2 intersect at a line lP 12 := nP 1 × nP 2, which is parallel to lmn = lmB − lnB.Moreover, nP 1, which is perpendicular to lP 12, is also perpendic-ular to lmn.These facts can be formalized leading to the following constraints: where emn1 and emn1 account for the fact that the assumed coplanarity of rmy, rn and lmn is an approximation (cf. Figure 3).
Note, when combining ( 13) with ( 15), the respective Jacobian for one timestep i is: and ∂e ni ∂l nS is computed similarly to (7).When considering the second and third block of J hi , it can be seen that the rank is not only depending on the angular velocities but also on the relative position vectors, as well as, the joint rotation axes.This gives a more complex expression, which is difficult to simplify.The rank of the first and last three columns can still be reduced as discussed in Section 3.However, by incorporating the constraint on lmB, the error emi as well as the estimation error of lmA are assumed to be reduced

EXPERIMENTAL RESULTS
The different IMU position calibration methods (Seel et al. and proposed) were first evaluated using noisy synthetic measurements simulated at 100 Hz for a setup as illustrated in Figure 2. In order to provide sufficient variations of the input data for the optimization problem, the measurements were downsampled to 10 Hz.The simulator output comprises 3D accelerometer, gyroscope and magnetometer measurements, as well as, the global orientations of each segment with respect to a fixed reference frame G.These values are calculated based on the inputs to the simulator, which include joint angles (generated from sine functions with different scales in each DOF), IMU to body poses, and also measurement noise specifications comparable to those of the IMUs used in the real data tests below.
In order to evaluate the accuracy and robustness of the different methods for different types of movements, two sets of joint angle sequences were simulated, each with 100 randomly generated IMU to body poses: Type 1 provides low variations among the DOFs of each segment, whereas the knee joint is modeled with one DOF.This resembles a natural calibration movement.Type 2 provides high variations among the DOFs and, in addition, models both joints with three DOFs.Hence, the type 1 set provides more realistic but suboptimal movement conditions, whereas the type 2 set is expected to work well with both calibration methods.
Table 1: The angular velocity variations for the simulated and the real IMU measurements: 12 represents the variation between the first and the second DOF for each IMU A, B, C. Each number expresses the mean of squared differences between the simulated/measured angular velocities for all the timesteps.
In order to extract the joint rotation axes required in (14), an extra phase of movements has been added in which each segment rotates in all possible DOFs separately.More details about the simulated joint trajectories in terms of resulting angular velocity variations are provided in Table 1.
Our proposed method, as well as, the methods of Seel et al. were implemented in Matlab using lsqnonlin with the Levenberg-Marquardt algorithm.In the following, the behavior of these methods are evaluated through several experiments.
The effect of only introducing the relative IMU orientations as formalized in Equation ( 13) without coupling the two estimation problems is illustrated in Figure 4, where the Root Mean Squared Error (RMSE) for the estimated IMU positions is provided for 10 different mounting orientations using type 2 trajectories.This shows already an improvement of the proposed method with respect to accuracy and robustness under optimal movement conditions.
Figure 5 illustrates the RMSE when applying the different calibration methods to the two different types of movements.Here, our method includes also the constraint introduced in Section 4.2.The results show that higher movement variations during calibration improve the performance of both algorithms.However, with less variations our proposed method outperforms the ones of Seel et al..This counts for both the original method in [21] and the correction described in Equation ( 12) [20].
Moreover, the convergence of both calibration methods was tested     using 100 random initial values.Figure 6 presents the results of this test for measurements simulated from one IMU to body configuration.This shows that the two methods are not generally sensitive to the initial values.Nevertheless, for type 2 trajectories, both methods converge to similar results, with less than 1 cm error.However, our proposed method converges faster.With suboptimal type 1 trajectories, both methods require a higher number of iterations; however, the method of Seel et al. converges with a large error of 61 cm, while our proposed method converges with only 3 cm error.
In order to evaluate the proposed method on real measurements, we captured a dataset of 10 trials performed by one subject using a prototype wearable system [18] and an optical reference system [3].The test setup is shown in Figure 7.In order to obtain the IMU poses, each IMU was rigidly connected with a rigid body marker.The IMUs were interconnected via textile cables and fixed to the garment using snap buttons.In order to reduce artifacts due to movement of the garment, the IMU-marker-sets were strapped firmly on the pelvis and the legs.We also used straps in order to attach marker clusters on anatomical landmarks around the hip and knee joints, from which we determined the joint CoRs.In each trial, for deducing the joint rotation axes, the subject first performed movements in each DOF of each joint separately.Then, in order to test the proposed IMU position calibration algorithm, arbitrary movements were performed simultaneously for the different segments.The resulting angular velocity variations are included in Table 1, showing that the real data corresponds to simulated type 1 movements.The global IMU orientations required in Equation ( 13) were calculated based on the method described in [17].The estimated IMU position vectors were then compared with the references calculated from the optical markers.Note that errors due to marker positioning are present, however, similarly for all the tested methods.Figure 8 illustrates the RMSE when applying the different calibration methods to the real measurements.The proposed method provides more accurate or comparable results in 80 percent of the trials, while a slightly worse performance can be observed in two of the trials.The difference between the performance of the proposed method on the simulated data and the real data can be described by two factors: First, due to magnetic disturbances, the orientation estimation of the real data is erroneous.Second, the position errors for the real data are obtained by comparing the lengths of the estimated vectors and not comparing component-bycomponent as for the simulated data.This was mainly due to the inaccuracies in finding the joint coordinates from the optical markers.A deeper analysis of the failure cases and further investigations with a higher number of subjects are required in order to draw final conclusions.

CONCLUSION
In this work we presented a practical IMU to body position autocalibration method, specifically developed for the lower body.The method can be used for smart training pants with multiple embedded IMUs, e.g. for application in exercise and fitness monitoring.On simulated data, the presented method shows a clear improvement in terms of accuracy and robustness compared to a previous method, in particular when it comes to suboptimal (low variation) movements during calibration.This has also been confirmed through preliminary tests with real data, though further experiments are required as part of our future work.Another route of future work relates to the fact that the proposed method uses the global IMU orientations, which are currently estimated separately and are considered as inputs to the position calibration.Here, the results might be improved by adding, as a final step, a joint optimization of the time-varying global IMU orien-  tations and the quasi-constant IMU positions with respect to the adjacent segments.Moreover, using a weighted least squares estimate, the presented method could be extended to provide uncertainty measures for the estimated calibration parameters of the different IMUs, which could then be used during tracking in order to improve the results.Here, varying uncertainties could account not only for different types of movements during calibration, but also e.g. for varying amounts of garment movements or other soft tissue artifacts due to different suit material and style, or different body shapes.

Figure 1 :
Figure 1: Illustration of the IMU position estimation problem.IMUs A and B are mounted on the adjacent segments of the joint m.

Figure 2 :
Figure 2: Illustration of the IMU position estimation problem specifically addressing the lower body.The spheroidal joint n refers to the hip (three rotation axes: rmx, rmy, and rmz), and the hinge joint n refers to the knee (one rotation axis: rn).The two joints are linked via the upper leg segment, with a fixed length of lmn.The IMUs A, B, C are mounted on the segments connected through the hip and knee joints.They are placed on the pelvis, upper leg, and lower leg, respectively.The quantities lmA, lmB, lnB, lnC are the IMU position vectors to be estimated.RGA, RGB, RGC are the orientations of the IMUs with respect to the global coordinate frame G.
(a) Idealized geometry with intersecting joint axes at hip and knee.(b) Realistic geometry with skew joint axes.

Figure 3 :
Figure 3: Illustration of the constraints formulated in (15): The left figure shows the idealized geometry, which is the basis for the constraints.It approximates the more realistic setup in the right figure with skew joint axes, leading to the plane P 1.The latter can be assumed to have a small angle α with respect to the ideal plane P 1.
Seel et al. measurement model

Figure 4 :
Figure 4: RMSE for 10 different IMU mounting orientations when applying the measurement model defined in (13) to type 2 movements.
method with type 1 trajectory our method with type 2 trajectory Seel et al. method with type 1 trajectory Seel et al. method with type 2 trajectory Seel et al. corrected method with type 1 trajectory Seel et al. corrected method with type 2 trajectory

Figure 5 :
Figure5: RMSE for 100 different IMU to body configurations when using the measurement model defined in(13) and the additional constraint in(15).

Figure 6 :
Figure 6: Convergence test with 100 random initial values for the IMU position.In this test the true value of lmA is [0.05, 0.09, 0.03] and the true value of lmB is [−0.17,0.004, 0.03].

Figure 7 :
Figure 7: Real test setup.The red arrows show three IMUs which are mounted on pelvis (1, not visible), upper leg (2) and lower leg (3).

Figure 8 :
Figure 8: Results on real measurements: RMSE for one IMU to body configuration and 10 trials performed by one subject.