Sparse Visual-Inertial Measurement Units Placement for Gait Kinematics Assessment

—This study investigates the possibility of estimating lower-limb joint kinematics and meaningful performance indexes for physiotherapists, during gait on a treadmill based on data collected from a sparse placement of new Visual Inertial Measurement Units (VIMU) and the use of an Extended Kalman Filter (EKF). The proposed EKF takes advantage of the biomechanics of the human body and of the investigated task to reduce sensor inaccuracies. Two state-vector formulations, one based on the use of constant acceleration model and one based on Fourier series, and the tuning of their corresponding parameters were analyzed. The constant acceleration model, due to its inherent inconsistency for human motion, required a cumbersome optimisation process and needed the a-priori knowledge of reference joint trajectories for EKF parameters tuning. On the other hand, the Fourier series formulation could be used without a speciﬁc parameters tuning process. In both cases, the average root mean square difference and correlation coefﬁcient between the estimated joint angles and those reconstructed with a reference stereophotogrammetric system was 3.5deg and 0.70, respectively. Moreover, the stride lengths were estimated with a normalized root mean square difference inferior to 2% when using the forward kinematics model receiving as input the estimated joint angles. The popular gait deviation index was also estimated and showed similar results very close to 100, using both the proposed method and the reference stereophotogrammetric system. Such consistency was obtained using only three wireless and affordable VIMU located at the pelvis and both heels and tracked using


I. INTRODUCTION
G AIT analysis is the most popular method used to quantify movement disorders and to assess rehabilitation effects over time [1].It is a clinical routine used in broad variety of pathologies in rheumatology, orthopedics, endocrinology, and neurology [2].When assessing gait, kinematics parameters such as joint angles, that are the basis of most analyses, or spatiotemporal parameters (e.g., velocity and stride length) are often used.When accurately estimated, those parameters can be combined to build assessment indicators to assist the therapists in their diagnosis.For example, the Gait Deviation Index (GDI) introduced by Schwartz et al. is increasingly used by the clinical community to estimate the level of gait pathology among disabled persons [3].The GDI compares gait kinematics data of a given patient to a healthy group, by estimating the distance between healthy strides and patient's strides.This distance (with respect to the healthy group), which reflects the pathology level, requires an estimation of the kinematics parameters to be assessed properly.Stereophotogrammetric Systems (SS) are considered as the reference tools in kinematics motion analysis.However, their use require a large financial investment and a complex experimental protocol.An alternative solution proposed by Cappozzo was called the minimum measured input model approach, that is aiming at maximizing the functional information extracted from simplified and affordable experimental protocols [4].

A. Related works
In this context, mainly two technologies were proposed: Inertial Measurement Unit (IMU) and RGB/RGB-Depth camerabased systems.Small, low-cost IMUs measure three accelerations and three angular velocities.Thus, in principle, it is possible to estimate their position and orientation through double and single integration, respectively.Unfortunately, drift jeopardizes the time integration of raw signals [5].To overcome this problem, Kalman adaptive filters are often proposed in the literature to ensure real-time estimates of the IMU orientation by fusing accelerometer and gyroscope data.Since only two orientation angles can be obtained with this approach [6], several studies have proposed to use at least one IMU for each segment of interest together with a multi-body model representing the human skeleton as a mean to reduce driftinfluence.For example, El Gohary et al. proposed to use a zero velocity update approach with an Unscented Kalman Filter relying on a kinematics model for arm motions tracking [7].Remarkable results were obtained when operating with a rigid robotic arm so that human soft-tissue-artifact and sensor-tosegment calibration errors were not considered in this study.To reduce the sensors count, Sy et al. proposed recently a method to estimate lower limbs 3D joint angles during gait [8].Despite a very interesting approach based on the use of an Extended Kalman Filter (EKF) including assumptions and constraints such as rigid body constraints, zero velocity update at the foot contact, assuming flat floor or making use of the total pelvis height, the presented results displayed a large Root Mean Square Difference (RMSD) superior to 10deg when compared to the same estimates using SS.These relatively poor results might be due to the fact that a proper sensor-to-segment calibration was missing and that a constant acceleration model was used.Indeed, IMUs require a calibration phase to align each sensor with respect to its segment.Static and functional calibration techniques based on a set of pre-defined postures and/or motions are commonly used in literature.However, these methods are time-consuming, might be limited by the subject's ability to perform pre-defined postures or motions, and thus can be hardly used with some population or in clinical settings for example [9].Large errors might also be due to sub-optimal weights tuning leading to an emphasized influence of IMU drift on joint angle estimates.Joukov et al. showed the influence of parameter's tuning when using an adaptive filter such as the EKF [10].Moreover, during gait, they showed the superiority of using Fourier series over using standard constant acceleration model of joint temporal evolution.The assumption that the joint trajectories can be modelled by using low-order Fourier series has already been made in the literature [10]- [12].Fourier representation was also largely used for spatio-temporal parameters assessment using IMUs [13].Fourier series coefficients were identified using a feedback adaptive frequency phase oscillator.As Joukov et al. discussed, the learning rate and parameters tuning of the oscillator was decisive in the performance of their method [10].Nevertheless, by performing a grid search, it was possible to fine-tune the filter and obtain a low RMSD of 2.4deg for the hip and knee joint angles.However, it should be noted that the proposed method could not estimate ankle angles, required a cumbersome sensor-to-segment calibration process and was based on one IMU per investigated segment.To cope with the IMU drawbacks, some authors have proposed to use camera based systems.Recently, markerless visual motion capture based on a RGB camera and machine learning algorithms were developed for estimating human motion in real-life scenarios [14], [15].However, such approaches are not yet real-time or accurate enough to reliably estimate dependent subject's joint angles and thus are not ready to be used yet in clinical settings for example.Moreover such methods need often to be trained on large scale datasets requiring a large computational power.Devices, based on RGB-Depth sensors such as the Kinect camera (Microsoft), with embedded skeleton tracking algorithms are available.However, it was shown that the Kinect camera is not accurate enough for rehabilitation assessment mainly due to segment length variations [16].For gait analysis, the Kinect sensor can only estimate timing and spatial characteristics [17] as joint angles are not reliable enough for a consistent analysis.Numerous methods have been developed in the literature to improve kinematic estimates of RGB-D like sensors.Most of them add the use of a multi-body model to obtain a more robust joint angle estimate [15], [16].However, acceptable differences can only be obtained for motions of very large amplitude [16].
Another camera based approach is homography, which is used in Augmented Reality (AR) to accurately estimate 3D rotations and translations of a given AR marker after a proper, yet simple, camera calibration process.
Nagymate et al. proposed the use of AR markers for gait analysis [18].They have used a commercial Multi-body Kinematics Optimization (MKO) approach that was fed with segments 3D pose obtained from very large AR markers that were located at each investigated segments.They displayed a RMSD of 2.3cm for the step length, and a joint angle estimate varying from 2.5deg to 6.7deg.The authors state that they have chosen to use such large AR markers to avoid occlusion and undetected markers.In a recent survey [19], it was suggested that merging RGB-Depth and IMU data could prevent occlusions and improve joint angle estimate accuracy [20].Ahmed et al. for example were able to estimate head and feet trajectories during overground walking, with a RMSD that was less than 2% of the walked distance, based on the fusion of only three sensors: a head mounted IMU-camera and two feet mounted IMU [21].However, the joint kinematics were not assessed in this study.Moreover, it is worthy noting that MKO and EKF were alternatively used with IMU [10] and AR markers [18], [22] (as well as with SS markers [23]).Both methods seem to provide similar results however the MKO can hardly be implemented in real-time [24].
In this context, it seems that there is a consensus in the community on the fact that a multi-body model based approach should be used to perform inverse kinematics with a sparse placement of sensors.Fusing camera based information, obtained from AR markers, and IMU data could lead to reliable and robust 3D sensor pose estimate.Proper sensor-to-segment calibration and EKF parameters tuning are also essential to provide accurate joint kinematics estimation.Our group has already proposed such an approach in a recent study [22] but it was dealing with simple arm tasks.Moreover, it required one sensor module per investigated segment.Nevertheless it allowed to show that it is possible to reduce IMU drift influence and to cope with visual occlusion problems.

B. Contribution
The objective of the present study is to demonstrate the validity of using sparse calibrated newly designed Visual-Inertial Measurement Units (VIMU) and adaptive filtering, during treadmill gait, to estimate relevant kinematics quantities that may be of great benefit to the clinical community.Moreover, besides new constraints implemented in the EKF, two formulations of the proposed state-vector are studied and compared: a classical constant acceleration model associated to an optimal tuning process of the EKF's parameters, and a model based on Fourier series that allows a simple tuning of the EKF's parameters based on a priori knowledge.Thus one of the sub-objective of this study is to assess if by having a more realistic state variable evolution model the tuning of the EKF parameters will be simplified.

II. METHOD
The overall principle of the proposed approach is depicted in Fig. 1   tracked with two affordable RGB cameras.The raw data collected from both IMUs and AR markers require a prior and unique calibration procedure that is based on a process detailed in our previous study [22] and briefly described in section II.B.A set of three AR markers is mounted on the top of each IMU sensor to form a so-called VIMU (see Fig. 3).Then, a practical subject-specific sensor-to-segment static calibration (section II.B.2) is done to define each VIMU pose in the segment coordinate systems (subject's segment axes and joint centres are therefore defined at this step).Finally, an EKF based on an eighteen Degrees-of-Freedom (DoF) multi-body kinematics model that is described in section II.A, fuses all measured quantities to estimate lower-limbs joint angles.The tuning of the EKF's parameters along with its state-vector expressed in two different versions with respect to the evolution of the joint trajectories are also compared and discussed (section II.C and section II.D).

A. Lower-limbs Multi-Body Model
A multi-body model of the lower-limbs is devised to relate the joint kinematics with the three VIMU measurements.It consists of = 7 segments articulated with = 18 DoF.As represented in Fig. 2.a, it is composed of two spherical joints (hip), two hinge joints (knee) and two universal joints (ankle).Moreover, the model floating base is located with respect to a global coordinate system ( 0 ) using three prismatic and three rotational virtual joints.Three VIMUs are assumed to be rigidly attached to the sacrum, and to the left and right heels, respectively (Fig 2 .b).A VIMU is composed of three AR markers and a single IMU (see Fig. 3 and Fig. 4).The pose (position and orientation) of each VIMU marker is directly estimated by the off-the-shelf ArUco library [46].The orientation provided by ArUco is expressed using the quaternion formalism as presented by Diebel [25].Thus, each VIMU allows measuring the absolute 3D position p 0 and orientation (quaternion) q 0 with respect to 0 of each AR marker as well as the 3D linear acceleration a and angular velocity with respect to the IMU coordinate system.Three coordinate systems are therefore used to represent VIMU different elements: a coordinate system corresponding to both IMU and central AR marker and two additional coordinate systems for the left and right AR markers, respectively.The forward kinematics model (FKM) and its first and second derivative (J and J) relating joint kinematics to the VIMU measurements are given as follows: where , , are the ( × 1) joints positions, velocities and accelerations vectors, respectively.R 0 (3 × 3) is the rotation matrix expressing the VIMU's central AR marker orientation in 0 .It is converted into a (4×1) quaternion vector q 0 through a custom mapping function f ().Details on the implementation of f () can be found in section 6.5 of the review paper proposed by Diebel [25].The Jacobian matrix and b is a (3 × 1) acceleration bias vector.

B. Calibration
Two easy-to-use calibration procedures are required to use the raw measurements effectively: • a calibration for each VIMU, independent of the subject, • a subject-specific sensor-to-segment calibration.On the one hand, the VIMU calibration is done without participation of the subject by positioning and moving the unit in the field of view of the RGB camera.On the other hand, thanks to the use of AR markers, an anatomical calibration of the VIMU can simply be done [26].

1) VIMU calibration:
The gyroscope and accelerometer were calibrated, without the use of any external device, as proposed by Tedaldi et al. to take into account inaccuracies in scaling factors, misalignment and non-orthogonality of axes, and non-zero biases [27].Moreover, as shown in our previous study [22], the coordinate systems of the IMU and of the central AR markers were aligned thanks to a dynamic calibration step making use of centripetal acceleration, to calculate the rotation matrix R (3 × 3).Two rigid transformations also exist between the side coordinate systems of the AR markers and of the central one: where R (3 × 3), q (4 × 1), r (3 × 1) represent the rotation matrix, its corresponding quaternion formulation, and position vector of = (left) or (right) AR marker expressed in the central one and obtained after calibration, respectively.These relative transformations were easily determined using 50 different static postures of the VIMU and solving an overdetermined system [28].These rigid transformations between the coordinate systems of the AR markers will be considered as rigid body constraints in the EKF (section II.C).Finally, a VIMU is fully calibrated given that all its measurements are expressed with respect to the coordinate system attached to the central AR marker frame .2) Sensor-to-Segment Calibration: As exemplify in Fig. 4 the information from an additional AR marker mounted on a large calibration wand, making its detection accurate and robust, can be used to pin-point anatomical landmarks similarly to what is done when using a SS [26].The calibration wand is used to pin-point eleven lower-limbs anatomical landmarks in 0 .These are the midpoint between the postero-superior illiac spines, the right and left antero-superior illiac spines, the lateral and medial femoral epicondyles, and the lateral and medial malleolus.This operation can be done in less than a minute.Joint center positions are then calculated using classical linear regression methods [29].The norm between consecutive joint center positions is used to determine the segments lengths in the multi-body model (section II.A, Fig. 2).Note that this sensor-to-segment calibration is valid as long as the VIMUs location on their corresponding segments are not modified.Otherwise, it has to be repeated.

C. Extended Kalman Filter
An EKF is proposed to estimate lower-limbs joint kinematics, gathered in the state-vector x, that minimizes at each time step the least-square difference between the VIMU measurements, gathered in y, and their estimate from the measurement model h, gathered in h, including the estimate of the VIMUs 3D positions p 0 , orientations q 0 , linear accelerations a and angular velocities (section II.A).
(3) The vector of estimated VIMUs measurements h is based on the forward kinematics of the multi-body model but also incorporates two type of constraints terms (c1 and c2).The elements of h 1 ( × 1) are used to solve kinematics indetermination.Since the number of VIMU is reduced, i.e. a single VIMU per kinematics chain, there exist multiple kinematics solutions leading to the same VIMU measurements in the Cartesian space.However, there is not an infinite number of solutions leading to the same VIMU pose as the proposed leg model has 6DoF similarly to the number of estimated quantities from a VIMU.This kinematics indetermination is well known as the elbow-up/elbow-down problem in robotics.It has been solved in different ways.In this study, we propose to push the joint solution toward the average joint angle value, , reported for walking in the literature [30], [31]: 13deg, 8deg, 11deg, 16deg, 1deg, 6deg for 7 , 8 , 9 , 10 , 11 , 12 angles of the left leg and, symmetrically for the right leg, respectively.In the context of the EKF, this is done by adding − 6 virtual measurements aiming to minimise the square difference with at each sample of time by setting the corresponding measurements in y to zero.This is equivalent to minimizing a cost function in an optimization process [32].
The elements in h 2 (42 × 1) additionally ensure rigid body constraints between the three AR markers in the VIMU.The position and orientation of the right and left AR markers with respect to the central AR marker should match their values obtained during the VIMU calibration (section II.B).
the measurement model was established, the joint kinematics estimation was conducted, first, through the prediction of the a priori state-vector x− thanks to a state evolution model , that is approximated by a linear form F, and a process covariance matrix Q: where P − is the a priori estimation of the error covariance matrix.Two cases were investigated for the state-vector formulation: • case 1 EKF based on a Constant Acceleration (EKF-CA) model: A first classical state-vector formulation proposes to estimate the joint angles, velocities and accelerations and assumed a constant acceleration model [23]: where , 0 and I are the sampling time, and the null and identity matrices, respectively.To reduce the drift effect, a random acceleration bias b (9 × 1) was added to better track the 3D measured accelerations of each VIMU (see Eq. 1) [7].
• case 2 EKF based on a Fourier Series (EKF-FS) model: Since gait can be assumed to be a quasi-periodic task, the angle of the ℎ joint can be represented using a Fourier series: where is the number of harmonics, and is the main motion pulsation.
As described in the experimental validation (see section III), the gait activity takes place on a treadmill at constant velocity.Consequently, Fourier series coefficients are expected to converge towards near constant values and identity matrix can be used as state-vector evolution model.The state-vector and its evolution model were defined as follows: with = 1, ..., and = 1, ..., Secondly, depending on the measurement error, v, and the measurement covariance matrix, R, the Kalman gain K was calculated to update the state-vector: where H is the Jacobian matrix ℎ x and S is the measurement covariance matrix.

D. EKF Parameters Tuning
The proposed EKF requires the tuning of the following elements: the initial state-vector value x0 , the initial estimation of the error covariance matrix P 0 , the process noise covariance matrix Q and the measurement noise covariance matrix R.
For the initial state-vector and error covariance matrices, it has been shown that their values mainly affect the initial part of the estimation [33].Thus, x0 was set with the result of the inverse geometric model calculated from the AR markers data collected at the first sample of time.P 0 was set equal to the identity matrix giving the same influence to all joints.
The tuning of the elements of Q and R is more sensitive as they impact the filter convergence rate and stability.The elements of R corresponding to the IMU data can be experimentally determined.However, the accuracy of the visual measurements is variable as a function of different aspects such as marker occlusion, motion blur, or distance and orientation with respect to the camera and thus is difficult to be assessed.The elements of R corresponding to the constraints c1 and c2 (section II.C) must also be tuned leading to a large number of parameters.Moreover, the ratio between the tuning of the elements of R and Q, that is more sensitive and task and joint dependent [23], affects the filter performance [33].The tuning of Q reflects the trust level that could be given to the state-vector evolution model .For example, the constant acceleration model is obviously incorrect for human motion leading to inconsistencies in the state-vector [35] and thus this should be reflected in Q.Similarly, despite the bracelet form, AR markers have shown some sensitivity to occlusions.Fortunately, these are easily detectable.Thus, if the loss of all three markers data was detected, the corresponding elements of the EKF noise covariance matrix R will be automatically updated to reflect the new state of input data.Consequently, two R matrices were considered: the one corresponding to the whole markers trajectories when no occlusion of at least one marker has occurred and another one when all AR markers were occluded.In the latter, the missing markers poses were replaced by their previous values and their associated R elements were given a large value (1 3 ) indicating a lack of trust in these measurements.Considering the assumption that the parameters tuning is task-dependent [23], which is commonly used in EKF for motion analysis, the data of one randomly chosen subject was used to identify the elements of Q.The same identified values were then used for the remaining subjects.Then, depending on the state-vector formulation, two methods to tune the parameters of R and Q were proposed in this study.
In the case 1, corresponding to the EKF-CA, an optimisation process aiming at tuning the parameters of the EKF-CA was used.The optimization process consists of finding V = [ Ω R ] (13 × 1) that minimizes the square difference between the estimated joint angles and the reference ones, obtained from the SS and MKO, respectively (see section IV).The optimization process is computationally costly.Thus, it is proposed to group some parameters and to set the elements corresponding to the constraints c2 to 1 −3 .The choice of this value reflects the fact that the reference values are almost exact.The thirteen parameters that need to be identified are: three parameters , , and corresponding to the joint angles, velocities and accelerations elements of the process covariance matrix Q, and ten parameters corresponding to the measured positions , quaternions , linear accelerations , angular velocities Ω , and to the six joints R = in constraint c1 with = 1...6 (assuming left and right symmetry).The problem of the optimal tuning of V then boils down to: where is the number of samples of the considered trial.V − and V + represent the lower and upper boundaries, respectively.
(18 × 1) is the vector of the reference joint angles calculated using the SS (see Section IV).It should be noted that the gradient of this formulation is difficult to calculate since the equation in the EKF needs to be integrated.Thus, Eq. 9 was solved using a local derivative-free solver proposed in the popular IPOPT library [34], leading to V = [165.233.18 0.17 0.13 0.17 1.09 0.46 [0.09 0.04 0.03 2.32 0. 1 1]] .
In the case 2 corresponding to EKF-FS, it is crucial to determine firstly the number of harmonics required to accurately estimate the joint angles, velocities and accelerations.To do so, it is proposed to a priori solve the following fitting problem for different values of : s.t.0 ≤ (10) The Fourier representation is also used to obtain an analytical estimate of the joint velocities and accelerations.Thereafter, the results of this fitting process were used to determine a good trade-off between the number of parameters and the accuracy.Then, the parameter tuning is only based on a priori knowledge of the biomechanics of the task.The diagonal elements in Q were set to a small value (0.01) to enforce their convergence to pseudo-constant values [35].This guarantees the convergence of the corresponding error covariance P to a value close to zero.The elements of R related to the measurements were all empirically set to 0.1.As it regards the constraints (c1) related to the difference between estimated and average joint angles, the elements were set accordingly to the joint angle amplitude such that R = [1 1 1 10 1 10] while the corresponding elements of the rigid body constraints (c2) were set to 1 −3 .

A. VIMU Prototype
The proposed motion capture system was developed to be affordable, user-friendly, compact and light-weight.VIMU prototypes are shown in Fig. 3.a.Each VIMU consists of one IMU (MPU-6050, Invensense) embedding a 3D gyroscope and a 3D accelerometer, and communicating at 100Hz using a Bluetooth wireless module.On the top of the IMU enclosure, measuring 4x4x2cm, three AR markers of 3.6x3.6cmare mounted in a 3D printed bracelet form.The bracelet weights 40g and has different forms depending on the targeted segment.Its form is convenient to improve the AR markers detection rate during motions involving complex rotations.VIMUs were attached to human body segments using a velcro strap.Two standard 60Hz high-definition RGB USB cameras (USBFHD01, ELP) were used to track the AR markers poses.A custom program is developed in C++ to read all VIMUs data synchronously at a frequency rate of 60 Hz.The total cost of the proposed system was below 100 e.

B. Experimental protocol
Ten healthy male volunteers (age 25.4 ± 4.1 years old, weight 71.2 ± 16.0 Kg, height 1.7 ± 0.08m) participated in the experiments.They gave their informed consent following ethical regulation for non-invasive experiments of the University of Tokyo Agriculture and Technology (Tokyo, Japan) and of the University of Paris-Est Créteil (Créteil, France) where experiments took place.Once the VIMUs and retroreflective markers were worn by the subject, at the different locations shown in Fig2.b, every subject was told to stand statically, for less than one minute, in any comfortable posture that allows detection of VIMUs AR markers using the RGB cameras.Two RGB cameras were setup on tripods to be able to visualize the whole scene and detect all AR markers.The two RGB cameras parameters were obtained thanks to a standard calibration process based on a chArUcoBoard [36].The data collected by the two cameras were thereafter expressed in the same coordinate system using a large AR marker located on the floor in the field of view of the two cameras (Fig. 4).Meanwhile, anatomical landmarks of interest were pin-pointed using the affordable calibration wand and also captured using the SS (see Fig. 2.b).As a result, joint centres' positions, segments' lengths, as well as the local pose of each VIMU sensor relative to its corresponding segment anatomical frame, could be accurately calculated (Section.II.B).Hereafter, all subjects performed normal gait on a treadmill for one minute while wearing their running shoes.Each participant was free to choose his natural and comfortable gait speed.

Accuracy Assessment
The ability of the two proposed EKF state-vector formulations to accurately estimate the lower-limbs joint kinematics was assessed by calculating the RMSD and the Pearson Correlation Coefficient (CC) between the estimated and reference joint angles obtained by means of the SS.This choice is due to the fact that we were interested in estimating the temporal joints evolution, independently of the alignment and calibration offsets between the global reference frames of the different motion capture systems.Note that the sensor data are presented in international system of units while the outputs of the proposed system are given in degrees and centimeters since these units are more easily interpreted in a clinical context.The reference joint angles were obtained using a standard reference SS (6 Prime cameras, Optitrack) and a state of the art MKO tracking 22 reflective markers [24] with the multi-body model shown in Fig. 2. The joint angles computed with the SS follow the Conventional Gait Model [37].Moreover, the knee flexion axis was defined orthogonal to the long axis of the femur and, at the hip, the joint coordinate system corresponded to 2-by-2 mutually orthogonal axes, similarly to the multibody model as shown in Fig. 2.a.The data of the three VIMUs, located at the pelvis and both heels of the subjects as shown in Fig. 2.b, were synchronously collected.The gait deviation index (GDI) [3] was calculated using the estimated joint angles obtained from EKF-CA and from EKF-FS.Both results were compared with the mean GDI obtained from reference SS data, with a leave-one-subject-out cross validation.Estimates of nine subjects were used to build the basis required to estimate GDI of the 10 ℎ subject.A total of approximately 260 healthy strides were used to build the GDI basis matrix on 9 subjects.The strides of the 10 ℎ subject were labelled as test strides, and used to calculate the GDI relatively to healthy strides.This operation was performed for each subject.
Moreover, based on the joint angles estimates and the forward kinematic model (see Eq. 1), the sacrum and heels positions were estimated and used to compute the gait stride lengths.The occurrences of heel strikes on the treadmill were first determined using the method developed in [38] and then used to estimate the corresponding stride lengths as in [39].
Finally, we have studied the differences between EKF-CA and EKF-FS results by verifying that they were statistically nonsignificant.This was done using the corresponding RMSD and CC of all estimated gait kinematics parameters and of the GDI through a paired t-test with an alpha level set a priori to 0.05.The Pearson Chi-square goodness-of-fit test was conducted to check for normal distribution in the data.Fig. 5. Mean RMSD between the joint angles, velocities and accelerations estimated using Fourier series expansion and those obtained using the SS as a function of the number of harmonics .

IV. RESULTS
The mean RMSD of the joint angles, velocities and accelerations calculated using Fourier series representation and those estimated using the SS as a function of the number of harmonics for a randomly selected subject are shown in Fig. 5.To minimise the RMSD for the joint trajectories the number of harmonics was set to = 4 accordingly to Fig. 5. Fig. 6 shows the tracking performance using the proposed EKF-CA of the data collected by the VIMU sensor attached to the left heel while walking on the treadmill.Overall the input data were correctly tracked.The RMSD on the IMU data tracking were on average of 0.017 m.s −2 and of 0.029rad.s−1 for the linear accelerations and for the angular velocities, respectively.For the AR marker positions and quaternions, the RMSD were on average of 3.8cm and 2.2 −2 , respectively.These relatively large tracking errors were due to the fact that the AR markers were not detected just after toe-off as exemplify by the zoomed windows of Fig. 6.In fact at this moment the AR markers were positioned almost perpendicularly to the cameras.Despite these tracking errors, the EKF was able to ] Samples Fig. 6.Comparison between the input data of the VIMU attached to the left heel (red) and their tracking with the EKF-CA (dashed black).These are the 3D position and quaternion, 3D linear acceleration and 3D angular velocity.Note that the vertical grey lines are used to separate data that are presented for all N samples on the same subplot.
accurately estimate the eighteen joint angles of the lower-limbs model described in Fig. 2.This is valid for both the CA and FS state-vector formulations using the optimal parameters tuning.Fig. 7 exemplifies the accurate estimates of the 6 DoF that were used to calculate the pose of the pelvis with respect to 0 for both approaches.In this figure, the corresponding average RMSD of the 3D pelvis position was equal to 0.4±0.0cm(CC 0.86 ± 0.20) and 1.6 ± 0.9deg (CC 0.75 ± 0.16) in the case of EKF-CA and 0.4 ± 0.0cm (CC 0.80 ± 0.30) and 2.1 ± 0.7deg (CC 0.76 ± 0.21) in the case of EKF-FS.Fig. 8 shows a representative comparison of the left leg joint angles estimate obtained using both EKF-CA and EKF-FS approaches.The corresponding joint angles were satisfactorily estimated with an average RMSDs and CCs of 3.3 ± 1.6deg and 0.76 ± 0.24 using the EKF-CA and of 3.3 ± 1.2deg and 0.74 ± 0.20 using the EKF-FS, respectively.Table I summarizes the results of the comparison between the lower-limb joint trajectories obtained with the proposed system and the SS, for each considered measure, based on the data from all subjects.Note that the results for the leg joints (i.e., hip, knee and ankle) were obtained over the average of left and right leg values.In general, results were similar between both EKF-CA and EKF-FS approaches.For the pelvis the average RMSDs were of 0.6±0.1cm(CC 0.84±0.21)and 2.0±0.25deg(CC 0.60 ± 0.09) when using the EKF-CA and 0.7 ± 0.0cm (CC 0.82±0.22)and 2.2±0.4deg(CC 0.63±0.18)when using the EKF-FS.For the leg joint angles, i.e., the hip, knee, and Sterephotogrammetric system EKF-CA EKF-FS Fig. 7. Representative comparison of the pelvis joint trajectories obtained for a random selected subject while walking on a treadmill using the stereophotogrammetric system (black) and the proposed approach based either on a constant acceleration model (red) or on Fourier series expansion (dashed blue).
ankle joint angles, the average RMSDs for all subjects were equal to 3.9 ± 1.0deg (CC 0.74 ± 0.19) and 3.9 ± 0.8deg (CC 0.71 ± 0.19) when using EKF-CA and EKF-FS, respectively.Table .II summarizes the mean and standard deviation RMSDs for right and left stride lengths estimates.These were calculated first using the VIMUs raw data (without processing) obtained from the AR markers pose estimation, then using both EKF-CA and EKF-FS estimated data, and were compared with respect to the SS.In both cases, the use of the EKF reduced the RMSD in stride length estimates.Specifically, the EKF-FS showed an improvement by 50% approximately passing from 2.6cm to 1.3cm.Moreover, using the joint angles estimates for the detected strides, the GDI was also calculated.Results presented in Table II show that all the mean of GDIs (the reference and the estimated ones), are very close to 100 indicating a normal walking pattern.
The paired t-test results for RMSD and CC of each joint trajectory estimated using both the EKF-CA and EKF-FS implementations are reported in Table I.For all joints, results showed nonsignificant statistical differences (p > 0.05) between EKF-CA and EKF-FS with the exception of the pelvic tilt angle (p = 0.001, p = 0.039).Nevertheless, considering the above results, the gait kinematics estimates do not seem to be deteriorated by the poor estimate of the pelvic tilt angle.Errors in pelvis tilt angle also mechanically affects hip flexion-extension.However, while pelvis tilt is of Sterephotogrammetric system EKF-CA EKF-FS Fig. 8. Representative comparison of the left leg joint angles obtained for a random selected subject while walking on a treadmill using the stereophotogrammetric system (black) and the proposed approach based either on a constant acceleration model (red) or on Fourier series expansion (dashed blue) limited amplitude during gait, hip flexion-extension is of large amplitude.Therefore, the segment angle errors do not translate into significant joint angle errors between EKF-CA and EKF-FS (p = 0.200 for RMSD and p = 0.071 for CC).A similar statistical test was performed in Table III for the RMSD of the stride lengths and for the GDI calculation.Both stride lengths' RMSD and GDI showed nonsignificant differences between EKF-CA and EKF-FS.The difference in GDI values resulting from each EKF implementation and SS data was also nonsignificant.On the other hand, the raw stride length values compared to their estimates based on the two EKFs (see Table III) revealed, as expected, significant differences.This emphasizes the improvement provided by the EKFs in either implementation.Fig. 9 shows the evolution of the nine Fourier series coefficients corresponding to the ankle flexion angle estimated using EKF-FS.The coefficients converge toward pseudo-constant values showing, as expected, the gait periodicity.

V. DISCUSSION
The proposed system proposes to fuse visual and inertial data of three VIMUs attached to the pelvis, left and right heels into an EKF and a lower-limb kinematic model.This allowed handling the inaccuracies that result from using affordable sensors.Visual occlusion was greatly limited thanks to the    bracelet form of the VIMUs.Rigid body constraints were integrated in the EKF formulation allowing to reconstruct marker loss based on its non-occluded partners within the same VIMU.This was possible as the rigid transformation between each two VIMU markers was priory calibrated (section II.B.1).Nevertheless, depending on the motion phases, the complete occlusion of the three AR markers could not be entirely avoided.Yet, the EKF was able to cope with such occlusion by relying more on IMU data thanks to the online update of the measurement covariance matrix R.This is shown in the zoomed capture of Fig. 6 where one can note that the EKF is not tracking inconsistent visual data.However, it is important to note that the occlusion was never longer than few seconds.
We have shown in a recent study that EKF approach was able to cope efficiently with such short occlusion [22].Thus we can conclude that the proposed method is not very sensitive to occlusion with the current setup.Moreover, constraints having the effects of joint springs pushing the joint angles toward mean values were added to the EKF as soft constraints in the measurement model h.This was motivated by the fact that a sparse placement of the sensors was used.Thus, for the same VIMUs tracking, different joint configurations, possibly within the joint limits, are possible.Basically one could imagine to have springs of different stiffness attached to each joint.This concept is similar to the popular potential field approach [40] commonly used within the robotics community.
The proposed approach has been validated with a direct comparison to the joint kinematics estimated using a reference SS.For consistency in the joint trajectories comparison, the SS data were also processed with MKO using the same multibody model.When tracking the same data, EKF and MKO are expected to provide very similar results [24].Over the fifteen lower-limbs joint angles the mean RMSD was of 3.5deg.These results show the ability of the proposed approach to handle the 3D IMU drift without magnetometers involvement while compensating a sparse sensors placement.
For both EKF implementations, the knee and hip flexionextension angles estimates presented a relatively high CC of 0.927 on average and a low RMSD ( 5deg).In contrast, the CCs of the ankle dorsi/plantar flexion (0.435) and the pelvic tilt angle (0.485) estimates are the poorest but with a satisfactory RMSD of 3.1deg and 2.0deg, respectively.Moreover the p-value of the pelvic tilt angle is the only angle displaying significant different results between the two EKF implementations.This might be explained by the small angle amplitude of this joint while walking on the treadmill.Moreover, the tracking error is not equally distributed among all joint trajectories due to constraints c1 (similar to joint springs) and to different parameters values in the covariance matrix R.
When compared to reference kinematics (bone pins or fluoroscopy), a systematic review reported errors between 1 and 22deg at the knee joint when tracking retro-reflective markers placed on both thigh and shank segments [24].The proposed method outperforms related studies dealing with a sparse sensor placement [8] (RSMD of 10deg) or when using a simple visual observation (RMSD of 9deg [10]) that is often still used by physiotherapists (see section I.A for further details on these studies).When using a single IMU per investigated segment and a quasi-periodic assumption it is possible to reach a better accuracy without monitoring the ankle joints (RMSD of 2.4deg) [10].Of course this increases the complexity and the price of the overall system.While similar results were observed by comparing the use of a constant acceleration model and a Fourier series expansion in the EKF (see Table .I), it is important to note that preprocessing to obtain these results is very different with the two methods.In the former, the model assumption is inherently inaccurate in human motion.The literature related to the use of EKF for motion analysis rarely proposes methods to model these inaccuracies and to reflect them in the parameters tuning.Cerveri et al. [41] were the sole authors to propose a frequency based approach but this requires to have a prior knowledge of the joint angle evolution, applies to cyclic motions and was not tested with low-cost sensors data.We have shown in a previous study [16], using a RGB-D camera and a markerless joint center estimate algorithm, that an optimization process outperforms frequency based and is absolutely required to tune the EKF parameters.Nevertheless, when using markerless data as input and despite an optimal tuning of the EKF parameters and a simplified experimental setup [16], it was not possible to achieve a joint estimate RMSD lower than 9.7deg on average when compared to a reference SS.Specifically, for tasks closely related to walking, a RMSD up to 20deg was observed for the hip joints.This emphasizes that markerless approaches are not yet accurate enough to estimate human motion, and the superiority of our proposed system being a good and affordable compromise between reference SS and markerless systems.
Unfortunately, an optimal parameters tuning approach is time consuming, since it requires to run the EKF several thousands of times leading to at least 10h of calculation time.Moreover, it can not be performed independently from a SS.Although this optimisation is task dependant and can performed once for different subjects, it requires that one subject performs the task wearing simultaneously VIMUs and retro-reflective markers at least once.In the latter, when using a Fourier series model, the EKF parameters can be directly set based on a priori knowledge.As exemplified in Fig. 9, the Fourier series converge towards pseudo-constant values, confirming, as expected, the quasi-periodicity of the joint trajectories.The quasi-periodicity of gait was already used for motion analysis in the literature [10]- [12] and seems to be a valid assumption even for patients.
Regarding values obtained by the optimization process for the tuning of the EKF-CA, they are difficult to be interpreted as the elements of Q and R are interdependent.However we can observe that for the noise covariance matrix they are of the same order of magnitude than the ones set empirically for the EKF-FS and that for the joint soft constraints the ratio of one for ten between some angles is also respected.Despite the differences in the two EKF implementations, the paired t-test results obtained for GDI and most gait kinematics (see Tables I and III) were not statistically significant.This suggests that both EKF-CA and EKF-FS can be used interchangeably for gait kinematics estimation while emphasizing the relevance of using an optimal tuning of the EKF parameters and/or a realistic process model that describes the joints evolution during movement.The position of the joint centers were estimated thanks to an anatomical calibration even if functional calibration movements could have been used [9].However, the proposed system targets patients that might not be able to perform exciting postures or motions [9].Nevertheless, for patients having a certain level of recovery, such as patients at home, the fact that with the proposed system position and orientation data could be directly provided by the AR markers would make the estimate of joint centers straightforward [9].Note that in both EKF formulations, the difference in joint angles estimates might be explained by the constant body segment lengths, computed statically as the distance between subsequent joint centre positions, as well as the VIMU-to-body-segment translational and rotational offsets.The latter may occur due to inaccuracies in AR markers measurements at initial static calibration, soft tissue artifacts, and the assumption of sensors-tobody-segments rigid connection.In such case, both segments lengths and VIMUs local poses could be modelled in the EKF but at the cost of increasing the state-vector dimensions.
The NRMSD of the strides lengths estimated with the SS and estimated when using the corrected EKF data was on average 1.5cm or 2% for both implementation cases (see Table II).This result is better than the 3cm reported in a recent study when using IMU data collected at both heel levels [42].Many studies have intended to estimate similar spatial gait parameters.For instance, Kose et al. reported a maximal error up to 3% for step length estimate [43].However, their approach, despite its ease of use, suffers from signal features dependency limiting its scope to normal walking.It relies directly on the IMU signal for heel strike detection and assumes a support leg in full extension.The approach proposed by Ahmed et al. uses a camera but only at the head level and thus it suffers from the same problems [21].These authors do not directly compute stride or step length, however their foot trajectory estimates, even if acquired over long travelled distance, showed a relatively large RMSD varying from 5 to 13cm.A recent study compared different spatio-temporal gait parameters between healthy and amputee patients [44].The patient's average stride length was about 1.05 m, while the normal walking was on average of 1.24 m.This difference is much larger than the accuracy provided by the proposed system.Thus, it would be easy to distinguish between the different walking modes.
Having correctly estimated joint trajectories and stride lengths allow to estimate other indexes such as the GDI.The GDI difference between each EKF and SS data was not statistically significant.Moreover, GDI results were all very close to 100.By definition, a GDI greater or equal to 100 means that the test strides are close enough to the healthy strides to be classified as normal ones.The obtained GDI illustrates that the joint trajectories, as well as stride lengths, estimated using the proposed approach have the potential to provide clinically useful indicators.In addition, from the EKF results, the evolution of the Fourier coefficients could be further used to evaluate the gait variability of patients [11].
The proposed system is composed of two RGB cameras.This number could be easily increased to match a regular motion capture system and target more complex tasks where the subject will not be facing the same set of cameras.Note that unlike classical stereophotogrammetric approaches, the proposed one is less sensitive to visual occlusion thanks to the simultaneous use of IMU data.However, as previously demonstrated [22], the subject should not be too long away from the camera(s) field of view as the IMUs drift will impact the accuracy of the joint angle estimate.This is inherent to any purely IMU based system [10].Another limitation of purely IMU based system is that the absolute position of the body with respect to the ground, and for instance to a forceplate, remains undetermined which impedes joint kinetics to be computed.The addition of externally fixed RGB camera(s) solves this issue.This way, the proposed affordable system shall be suitable for both kinematics and kinetics.The proposed system could be used outside of the laboratory but it will still require trained staff to perform the anatomical calibration.It could be used to monitor gait in small clinical centers or at doctor or physiotherapist practice for example.This is already the case of much more expensive devices that can monitor solely the knee joint angle [45].

VI. CONCLUSION
A new approach for motion analysis based on the use of VIMU sensors and an EKF has been proposed and validated in this paper.VIMU sensors can be considered a good compromise between, on the one hand, SS systems which are costly and cumbersome and, on the other hand, purely IMU based methods which requires one IMU per segment and advanced anatomical calibrations to reach an accuracy below 5deg.On the computational side, the propose EKF-FS represent one of the most advance implementations of this approach (including constraints and a clear rationale for the filter parameter tuning).

Fig. 1 .
Fig. 1.Overview of the proposed approach based on Visual-Inertial Measurement Units (VIMU) and an Extended Kalman Filter (EKF) taking advantage of the human body and task biomechanics.

Fig. 2 .
Fig. 2. (a) Lower-limbs multi-body model composed of segments, of lengths , = 1, .., , and Degrees-of-Freedoms (DoF) denoted , = 1, .., .(b) Experimental setup, three VIMU attached to the pelvis, left, and right heels in a bracelet form are collecting data during gait thanks to two RGB cameras located behind the treadmill.

Fig. 3 .
Fig. 3. (a) The VIMU prototype with three different bracelet arrangements to fit different subjects and/or segment sections.(b) Description of the coordinate systems associated to a VIMU.The coordinate system of the IMU sensor, left, and right AR markers are all linked to the central marker coordinate system through rotation matrices R and a translation vectors r.

Fig. 4 .
Fig. 4. Example of lower-limbs anatomical landmarks pin-pointed using a calibration wand base on an AR marker.

Fig. 9 .
Fig. 9. Representative evolution of the Fourier series coefficients corresponding to the ankle rotation joint ( 12 ) estimated by the EKF-FS.

-specific calibration Subject-specific calibration
, and is based on the simultaneous use of affordable IMUs, incorporating only 3D accelerometers and gyroscopes sensors, i.e., no magnetometers, and a set of AR markers

TABLE I RESULTS
OF THE COMPARISON BETWEEN EIGHTEEN LOWER-LIMBS JOINT TRAJECTORIES OBTAINED USING THE PROPOSED AFFORDABLE SYSTEM AND THE SS. RESULTS HAVE BEEN REPORTED AS MEAN ± SD OVER ALL SUBJECTS AND RESULTS OF HIP, KNEE, AND ANKLE JOINTS HAVE BEEN OBTAINED OVER THE AVERAGE OF LEFT AND RIGHT LEG VALUES.RESULTS OF THE PAIRED T-TEST BETWEEN EKF-CA AND EKF-FS FOR RMSD AND CC OF EVERY JOINT TRAJECTORY HAVE ALSO BEEN REPORTED.

TABLE II COMPARISON
BETWEEN USING THE AFFORDABLE RAW/SS DATA AND THE DATA ESTIMATED BY EKF, FOR OBTAINING THE STRIDE LENGTHS AND GAIT DEVIATION INDEX (GDI), OVER ALL SUBJECTS AND STRIDES.