Estimating Lower Limb Kinematics using a Reduced Wearable Sensor Count

IEEE TRANSA CTIONS 1 Estimating Lo wer Limb Kinematics using a Reduced W earable Sensor Count Luke Sy Student Member , IEEE , Michael Raitor Student Member , IEEE , Michael Del Rosario Member , IEEE , Heba Khamis Member , IEEE , Lauren Kark, Nigel H. Lovell F ellow , IEEE , Stephen J. Redmond Senior Member , IEEE Abstract — Goal: This paper presents an algorithm for accu- rately estimating pelvis, thigh, and shank kinematics during walking using only three wearable inertial sensors. Methods: The algorithm makes novel use of a constrained Kalman filter (CKF). The algorithm iterates through the pr ediction (kinematic equation), measurement (pelvis position pseudo-measurements, zero velocity update, flat-floor assumption, and covariance lim- iter), and constraint update (formulation of hinged knee joints and ball-and-socket hip joints). Results: Evaluation of the algo- rithm using an optical motion capture-based sensor-to-segment calibration on nine participants ( 7 men and 2 women, weight 63 . 0 ± 6 . 8 kg, height 1 . 70 ± 0 . 06 m, age 24 . 6 ± 3 . 9 years old), with no known gait or lower body biomechanical abnormalities, who walked within a 4 × 4 m 2 capture area shows that it can track motion relative to the mid-pelvis origin with mean position and orientation (no bias) root-mean-squar e err or (RMSE) of 5 . 21 ± 1 . 3 cm and 16 . 1 ± 3 . 2 ◦ , r espectively . The sagittal knee and hip joint angle RMSEs (no bias) were 10 . 0 ± 2 . 9 ◦ and 9 . 9 ± 3 . 2 ◦ , respecti vely , while the corr esponding correlation coefficient (CC) values were 0 . 87 ± 0 . 08 and 0 . 74 ± 0 . 12 . Conclusion: The CKF-based algorithm was able to track the 3D pose of the pelvis, thigh, and shanks using only three inertial sensors worn on the pelvis and shanks. Significance: Due to the Kalman- filter -based algorithm’s low computation cost and the relative con venience of using only three wearable sensors, gait parameters can be computed in real-time and remotely for long-term gait monitoring. Furthermore, the system can be used to inform real- time gait assistive devices. Index T erms —Constrained Kalman filter , Gait analysis, Motion capture, Pose estimation, W earable devices, IMU I . I N T RO D U C T I O N Gait analysis (GA), part of the study of human mov ement, is a valuable clinical tool that can help diagnose mov ement disorders and assess surgical outcomes [1]. Its value has been demonstrated in studies related to: knee and hip osteoarthri- tis [2]; falls risk [3]; and children with cerebral palsy [4]. Such studies quantitativ ely analyze gait parameters which are typically categorized as either spatiotemporal (e.g., speed and stride length), kinematic (e.g., knee flexion), or kinetic (e.g., joint moments) parameters [1]. GA is also used in interactiv e and possibly unsupervised rehabilitation, where measured gait parameters can driv e real-time feedback (e.g., [5]). Gait parameters can be objectiv ely measured by human motion capture systems (HMCSs). HMCSs can be organized into two broad groups: non-wearable and wearable. Non- wearable systems, predominantly camera-based, can estimate position with up to millimeter accurac y , if well-configured and calibrated. Optical motion capture (OMC) systems are currently the industry standard for motion capture in the laboratory setting. Howe ver , the capture volume is usually restricted to tens of cubic meters, thereby limiting the duration and range over which a person’ s gait can be assessed. The unfamiliarity or unusualness of this en vironment may also lead the person to move in manners different to their normal gait [6]. W earable systems primarily use inertial measure- ment units (IMUs) which measure the linear acceleration, angular velocity , and magnetic field. Recent technological advancements have dramatically miniaturized IMUs, making it arguably the most promising technology for long-term and continuous tracking of human movement during ev eryday life [7]. Furthermore, wearable systems are typically cheaper and easier to configure, enabling the capture of large number of participants in research studies (e.g., [8]). W earable systems can measure gait parameters deriv ed from one (e.g., one foot) or multiple body segments. This paper will focus on HMCSs that monitor the kinematics of the lower body (e.g., hip and knee angles). Commercial wearable HMCSs attach one sensor per body segment (OSPS) (i.e., fiv e sensors to track the pelvis, thighs, and shanks) [9, 10]. A sensor-to-segment calibration relates the IMU frame with the body segment frame, and then each body segment’ s 3D orientation, angular rate, and acceleration are estimated based on the assumption that the sensor is rigidly attached to the body segment it is monitoring. The sensor -to-segment calibra- tion can be achiev ed through either careful manual placement, use of a relativ ely tedious set of calibration mov ements, or an automatic online calibration method which lev erages assump- tions about joint kinematics [9, 11]. In general, position is inferred through double integration of the measured accelera- tion and/or by vector addition of segment ax es where segments make a kinematic chain, while orientation is inferred through integration of measured angular rate. In practice, position estimates will experience sev ere drift if not provided with at least an intermittent external reference of position or velocity , and orientation estimates will drift if not corrected using an accelerometer (using acceleration due to gra vity to correct the pitch and roll) or a magnetometer (using magnetic north to correct yaw) [12, 13]. Acceptable measurements v ary between clinical applications. Howe ver , McGinley et.al. believ e that in most common clinical situations an RMSE of 2 ◦ or less is considered e xcellent, and RMSE between 2 ◦ and 5 ◦ is acceptable b ut may require consideration in data interpretation [14]. OSPS typically ha ve acceptable results ( < 5 ◦ RMSE no bias and adjusted for repeatability) [15]. The number of IMUs required to be attached by OSPS configurations might be considered too cumbersome or incon venient for routine daily use, and would currently be considered too expensi ve by the consumer market, costing tens of thousands of US dollars for a state-of-the-art full-body system. W e expect that a reduced- IEEE TRANSA CTIONS 2 sensor-count (RSC) configuration, where sensors are placed on a subset of the total number of body segments, will improve user comfort while also reducing setup time and system cost. Howe ver , reducing the number of sensors inherently reduces av ailable kinematic information, from which body kinematics are inferred. In the following, we briefly revie w works which infer body kinematics for RSC configurations, arranged under data-driv en and model-based approaches. A. Appr oaches to kinematic inference for RSC systems 1) Data-driven kinematic infer ence: These approaches at- tempt to reconstruct the missing data by comparing the data obtained to an existing motion database (DB) or by some learned model [16 – 18]. It assumes that kinematic parameters are correlated, such that a pattern of available measurements allows the missing kinematic parameters to be inferred from previously acquired distributions representing a broad range of common human mov ements. T autges et al . searched for the best motion from their DB through a nearest-neighbor search of poses at multiple time steps [16]. W ouda et al . utilized fiv e IMUs attached to the pelvis, arms, and legs, and inferred the orientation of the body se gments without sensors attached using a k-nearest-neighbour algorithm and artificial neural network [17]. Huang et al . dev eloped a deep neural network with a bi-directional recurrent neural network architecture that reconstructs full body pose using six IMUs located on the head, wrists, pelvis, and ankles [18]. Data- driv en approaches have been shown to efficiently reconstruct realistic motions making them useful for animation-related applications [16 – 18]. Howe ver , these approaches naturally hav e a bias toward motions already contained in the DB, inherently limiting their use in pathological gait monitoring. 2) Model-based kinematic inference: These approaches at- tempt to reconstruct body motion using kinematic and biome- chanical models. Hu et al . modeled each leg as three linked segments confined in the sagittal plane. It tracked pelvis and ankle positions from four IMUs (two at pelvis and one for each foot) while making velocity assumptions at certain gait cycle stages, and calculated hip, knee, and ankle angles using in verse kinematics [19]. Howe ver , the system could only track motion in the sagittal plane, which can cause issues when tracking activities of daily li ving (ADLs) such as side or diagonal steps. Marcard et al . modeled the body as a kinematic chain of rigid segments linked by 24 ball and socket joints [20]. A window-based optimization algorithm then calculates the model pose that best matches the orientation and acceleration measurements from six IMUs located on the head, wrists, pelvis, and ankles, while satisfying multiple joint angle con- straints. Howe ver , it has long computation time, and hence only suitable for offline analysis applications. B. Novelty This paper describes a nov el model-based algorithm based on a constrained Kalman filter (CKF) to estimate lower body kinematics using a RSC configuration of IMUs. A Kalman Filter (KF) is the optimal state estimator for linear dynamic systems. It is capable of handling uncertainty in measurements (e.g., IMU measurements) and indirect observ ations (e.g., step detection and biomechanical constraints) [21]. It also outputs model state estimates at ev ery time step, and is therefore a suitable algorithm for real-time applications. This design was motiv ated by the need to dev elop a gait assessment tool using as few sensors as possible, ergonomically-placed for comfort, to facilitate long-term monitoring of lower body movement. I I . A L G O R I T H M D E S C R I P T I O N The overall objective is to estimate the orientation of the pelvis, thighs, and shanks with respect the world frame, W (see Fig. 2); their orientation completely defines their relativ e positions, since they form a mechanical linkage through hip and knee joints. At each time step, we attempt to predict the location of the shanks and pelvis in 3D space through double integration of their linear 3D acceleration, as measured by the IMUs attached to these three body segments; note, a pre-processing step fuses the IMU accelerometer, gyroscope, and magnetometer readings to estimate the IMU (and hence associated body segment) orientation so that gravitational acceleration can be removed, gi ving 3D linear acceleration of each IMU (and associated body segment) in the world frame. T o mitigate positional drift due to sensor bias and errors in the double integration of acceleration, the 3D velocity and the height above the floor of the ankle is zeroed whenev er a footstep is detected on that leg. This zeroing of ankle velocity and correction of ankle position is implemented within the KF framework by making a pseudo-measurement; to control the otherwise ev er-gro wing error covariance for the pelvis and ankle positions, a pseudo-measurement equal to the current position state estimate with a fixed covariance is made. Next, we make two novel contributions which take us beyond state-of-the-art: (1) W e mitigate position drift in the pelvis, caused by error accumulation during double integration of acceleration, by making a noisy pseudo-measurement of 3D position of the pelvis as being the length of the unbent leg(s) abov e the floor , with an XY position of the pelvis being mid- way between the two ankle locations; in practical terms, to say we make a noisy pseudo-measurement means that we assume there is large noise cov ariance associated with this pseudo-measurement, which we set empirically to represent our uncertainty in this assumption. This effecti vely puts a soft constraint on the 3D pelvis position relativ e to the two ankles. (2) W e ensure that the positions of the shanks and pelvis are such that they obey the assumed biomechanical constraints of the lower body , being that the hip joints are ball-and-sockets and the knee joint is a hinge with one degree of freedom and a limited range of motion (R OM). These constraints result in three equations (one of which is non-linear) in the XYZ position coordinates of the shanks and pelvis which ensure that: (i) the knee-to-hip distance is equal to the thigh length; (ii) long axis of the thigh is orthogonal to the mediolateral axis of the shank (i.e., knee is a hinge), and; the knee joint angle is within a biomechanically-realistic range. Again, these constraints can be implemented in the KF framework as pseudo-measurements, although the non-linear constraint equation must be linearized first; furthermore, since we are IEEE TRANSA CTIONS 3 linearizing non-linear constraint equations, it is necessary to iterate this pseudo-measurement procedure at each time step until these non-linear constraints are satisfied to within some tolerance. A final important point is that throughout the entire algorithm the orientation estimated using the IMU during the pre-processing stage is tak en to be free of error , which is obviously untrue and will be the topic of future work. The three stages of the algorithm (i.e., prediction, measurement, and constraint updates) are shown in Fig. 1. Fig. 1. Algorithm overvie w which consists of pre-processing, KF based state estimation, and post-processing. Pre-processing calculates the inertial acceleration of each IMU (fixed on various body segments), body segment orientation, and step detection from raw acceleration, a k , angular velocity , ω k , and magnetic north heading, h k , measured by the IMU. The KF-based state estimation consists of a prediction (kinematic equation), measurement (pelvis position pseudo-measurements, covariance limiter , intermittent zero- velocity update, and flat-floor assumption), and constraint update (thigh length, hinge knee joint, and knee range of motion). Post-processing calculates the thigh orientation, q lt and q rt . In the next subsection, the chosen model and notation are defined, followed by an illustration of the selected IMU placement, and a detailed explanation of the algorithm. A. Physical model and mathematical notation The model makes use of several important joints and points of the lo wer body; namely , the mid-pelvis (taken as mid-point between hip joints), hip joints, knee joints, and ankles (i.e., P = { p mp , p lh , p rh , p lk , p rk , p la , p ra } ). D p C ∈ R 3 denotes the 3D position of joint C with respect to frame D . If frame D is not specified, assume reference to the world frame, W . Additionally , the lower body is comprised of several seg- ments, namely , the pelvis, thighs, and shanks. Each seg- ment has its own length, measured joint-to-joint (i.e., L = { d p , d lt , d rt , d ls , d rs } ), and orientation, represented using quaternions (i.e., O = { q p , q lt , q rt , q ls , q rs } ). Quaternion D q C ∈ R 4 denotes the orientation of frame C with respect frame D and will be represented as bold-faced and italicized. If frame D is not specified, assume reference to the world frame W . Note that unit quaternions can be interchangeably expressed as a rotation matrix, i.e., q S ⇐ ⇒  r S x r S y r S z  where r S x , r S y , and r S z ∈ R 3 are the basis vectors of the segment’ s local frame, S , with respect to the world frame, W . Details on how to con vert between quaternion and rotation matrix representations of 3D orientation can be found in [22]. B. IMU Setup The system will be instrumented with an IMU at each of the three locations: (1) the sacrum, approximately coincident with the mid-pelvis; (2 & 3) the left and the right shanks, just above the ankles, as shown in Fig. 2. Note that the noise from IMU measurements (e.g., acceleration) when stationary can be approximated by a normal distribution [23, Ch. 2], but more sources of error , not necessarily normally distributed as assumed by the KF , comes into play during dynamic mo- tion [24]. Nev ertheless, several publications have successfully tracked pose while approximating IMU measurements as a normal distribution [9, 13]. Fig. 2. Physical model of the lower body used by the algorithm. The circles denote the joint positions P . The solid lines denote instrumented body segments, whilst the dashed lines denote segments without IMUs attached. Estimating the orientation and position of the uninstrumented thighs is the main contribution of this paper . C. Pr e-pr ocessing Pre-processing inv olves the calculation of orientation and inertial acceleration of instrumented body segments from raw IMU measurements, and step detection for both feet. The information obtained is then used by the KF state estimator . 1) Body segment orientation: The reference frame of each IMU may not be coincident with that of the attached body segment. This misalignment is corrected using a sensor-to- segment realignment similar to [9, Sec. II-B Eq. (1)] where the orientation of the instrumented body segment at time step k , i.e., IO = { ˘ q p k , ˘ q ls k , ˘ q rs k } , is calculated from (i) measured sensor orientation of sensor frame, S , in the world frame, W , and (ii) a kno wn body segment orientation (e.g., measured by an OMC system). 2) Inertial acceleration of IMUs in the world frame: is due to body movement, as distinct from gravitational acceleration which is measured by the IMU’ s accelerometer . It is calculated by subtracting acceleration due to gravity and then expressing the acceleration of the instrumented body se gment in the world frame similar to [9, Eq. (3)]. 3) Step Detection (SD): detects if the left or right foot is in contact with the ground. It can be used to gain knowledge about the velocity and position of the ankle, treated as an intermittent measurement in the KF . This component can use any third-party SD algorithm. A revie w on different SD algorithms can be found in [25]. IEEE TRANSA CTIONS 4 D. System, measur ement, and constraint models The state estimator is based upon a CKF , which combines a system model (Section II-E1), a measurement model (Section II-E2), and a kinematic constraint model (Section II-E3). The interaction between these models is illustrated in Fig. 1. Here, we first introduce these models before presenting the CKF . x k = Fx k − 1 + Gu k − 1 + w k − 1 (1) y k = H k x k + v k , D k x k = d k (2) where k is the time step; x k is the state vector; y k is the measurement vector; w k and v k are zero-mean process and measurement noise vectors with covariance matrices Q k and R k , respectively (denoted as w k ∼ ( 0 , Q k ) and v k ∼ ( 0 , R k ) ; F , G , and H k are the state transition, input, and measurement matrices, respectiv ely; and D k x k = d k are the equality constraints that state x k must satisfy . Eqs. (1) and (2) are the standard system model used in CKF . A more detailed description can be found in [21]. The state variables in x k models the position and velocity of the instrumented body segments (i.e., the 18 × 1 vector x k =  ( p mp k ) T , ( p la k ) T , ( p ra k ) T , ( v mp k ) T , ( v la k ) T , ( v ra k ) T  T ). E. Constrained Kalman filter (CKF) The a priori (predicted), a posteriori (updated using mea- surements), and constrained state (satisfying the constraint equation) for time step k are denoted by ˆ x − k , ˆ x + k , and ˜ x + k , respectiv ely . The KF state error a priori and a posteriori co- variance matrices are denoted as P − k and P + k , respectively . An implementation will be made av ailable at: https://git.io/Je9VV. 1) Pr ediction step: estimates the kinematic states at the next time step. This prediction does not respect the kinematic constraints of the body , so joints may become dislocated after this prediction step. The input is the inertial acceleration, expressed in the world frame, as described in Section II-C2 (i.e., u k =  ( a mp k ) T ( a la k ) T ( a ra k ) T  T ). The matrices F , G , and Q are defined follo wing the kinematic equations of motion for the pelvis and ankles as shown below; that is v elocity is the integral of acceleration, and position the integral of velocity: F =  I 9 × 9 ∆ t I 9 × 9 0 9 × 9 I 9 × 9  , G =  ∆ t 2 2 I 9 × 9 ∆ t I 9 × 9  (3) Q = G diag ( σ 2 acc ) G T (4) where σ 2 acc is a 9 × 1 vector of v ariances in the ac- celerometer measurements (assumed equal variance for all accelerometry signals), and ∆ t is the sampling interv al. The a priori state estimate, ˆ x − k , and state estimate error cov ariance, P − k , are calculated following the standard KF equations [26]: ˆ x − k = F ˜ x + k − 1 + Gu k , P − k = FP + k − 1 F T + Q (5) 2) Measur ement update: estimates the next state by: (i) enforcing soft pelvis position constraints based on ankle x and y position, and initial pelvis height, and by; (ii) utilizing zero ankle velocity and flat floor assumptions whenev er a footstep is detected. For the cov ariance update, there is an additional step to limit the a posteriori cov ariance matrix elements from growing indefinitely and from becoming badly conditioned. Firstly , the soft mid-pelvis position constraints encourage: (i) the pelvis x and y position to approach the av erage of the left and right ankle x and y positions; and (ii) the pelvis z position to be close to the initial pelvis z position as time k = 0 (i.e., standing height z p ). These constraints are implemented by H mp and y mp as shown in Eqs. (6) and (7) with measurement noise v ariance σ 2 mp ( 3 × 1 vector). Secondly , if a left step is detected, the left ankle velocity is encouraged to approach zero, and the left ankle z position to be close to the floor le vel, z f . These assumptions are implemented using H ls and y ls as shown in Eq. (8) with measurement noise variance σ 2 ls ( 4 × 1 vector). Note that H rs and y rs can be constructed in a similar fashion to Eq. (8). H mp =  pelvis, left and right ankle pos. columns z }| { I 3 × 3 − 1 2 I 2 × 2 0 − 1 2 I 2 × 2 ... 0 1 × 2 0 0 1 × 2 ...  (6) y mp =  0 1 × 2 z p  T (7) H ls =  ... left ankle pos. and vel. columns z }| { 0 3 × 2 0 3 × 1 I 3 × 3 ... ... 0 1 × 2 1 0 1 × 3 ...  , y ls =  0 3 × 1 z f  (8) As floor contact (FC) varies with time, so too H k varies with time, as shown in Eq. (9). Measurements y k and state variances σ 2 k are constructed similarly to Eq. (9). H k =          [ H T mp ] T no FC [ H T mp H T ls ] T left FC [ H T mp H T rs ] T right FC [ H T mp H T ls H T rs ] T both FC (9) The a priori state ˆ x + k is calculated follo wing the traditional KF equations: R k = diag ( σ k ) , K k = P − k H T k  H k P − k H T k + R  − 1 (10) ˆ x + k = ˆ x − k + K k  y k − H k ˆ x − k  (11) Lastly , the cov ariance limiter prevents the cov ariance from growing indefinitely and from becoming badly conditioned, as will happen naturally with the KF tracking the global position of the pelvis and ankles without any global position reference. At this step, a pseudo-measurement equal to the current state ˆ x + k is used (i.e., y lim = ˆ x + pos,k , implemented by H lim as shown in Eq. (12)) with some measurement noise of variance σ 2 lim ( 9 × 1 vector). The cov ariance P + k is then calculated through Eqs. (13)-(15). H lim =  I 9 × 9 0 9 × 9  (12) H 0 k =  H T k H T lim  T , R 0 k = diag ([ σ 2 k σ 2 lim ]) (13) K 0 k = P − k H 0 T k  H 0 k P − k H 0 T k + R 0  − 1 (14) P + k = ( I − K 0 k H 0 k ) P − k (15) 3) Satisfying biomechanical constraints: After the predic- tion and measurement updates of the KF , above, the body joints may ha ve become dislocated, or joint angles extend beyond their allowed range. This part of the algorithm corrects IEEE TRANSA CTIONS 5 the kinematic state estimates to satisfy the biomechanical constraints of the human body . This amounts to projecting the current a posteriori state ˆ x + k estimate onto the constraint surface, guided by our uncertainty in each state variable en- coded by P + k . The constraint equations enforce the following biomechanical limitations: (i) the length of estimated thigh vectors ( || τ lt || and || τ rt || ) equal the thigh lengths d lt and d rt ; (ii) both knees act as hinge joints; and (iii) the knee joint angle is confined to realistic range of motion and is only allo wed to decrease (i.e., prohibit bending the knee further) to prevent the pelvis drifting towards the floor . For the sake of brevity , only the left leg formulation is shown; the right leg formulation can be deriv ed in a similar manner . Firstly , the constraint for the length of the estimated thigh vector is shown in Eq. (17) where τ lt z ( ˜ x + k ) is the thigh vector . The thigh vector (Eq. (16)) was obtained by subtracting the knee joint position from hip joint position. The hip joint position, p lh , is obtained by adding half the pelvis length along the pelvis y axis, r p y , from the mid pelvis position, p mp , while the knee joint position, p lk is obtained by adding the shank length along its z axis, r ls z , starting from the ankle joint position, p la . As this equation is non-linear, it must be linearized to fit in our KF estimator . Linearization is done via T aylor series approximation, as described in [27, Sec. 3 Eqs. (66) and (67)] leading to a system of linear equations D ltl,k ˜ x + k = d ltl,k (Eq. (18) and (19)). τ lt z ( ˜ x + k ) = hip joint pos. z }| { p mp + d p 2 r p y − knee joint pos. z }| { ( p la + d ls r ls z ) (16) c ltl ( ˜ x + k ) = q τ lt z ( ˜ x + k ) T τ lt z ( ˜ x + k ) − d lt = 0 (17) D ltl,k = ∂ c ltl ( ˜ x + k ) ∂ x = h τ lt z ( ˜ x + k ) || τ lt z || − τ lt z ( ˜ x + k ) || τ lt z || 0 1 × 12 i (18) d ltl,k = − c ltl ( ˜ x + k ) + ∂ c ltl ( ˜ x + k ) ∂ x ˜ x + k (19) Secondly , the constraint for the hinge knee joint enforces the long ( z ) axis of the thigh to be perpendicular to the mediolateral axis ( y ) of the shank, as shown in Eq. (20). This formulation is similar to [28, Sec. 2.3 Eqs. (4)]. Using Eq. (16) to expand the thigh vector τ lt z in Eq. (20) and rearranging it in terms of D lkh,k ˜ x + k = d lkh,k giv es us Eqs. (21) and (22). Since the constraint is linear , no linearization is needed. c lkh ( x k ) = τ lt z · r ls y = 0 (20) D lkh,k =  r ls y − r ls y 0 1 × 12  (21) d lkh,k = − ( d p 2 r p y − d ls r ls z ) · r ls y (22) Thirdly , the constraint for the knee range of motion (ROM) is enforced if the knee angle is outside the allowed ROM or is increasing (i.e., prohibit bending the knee further). This implementation is similar to the acti ve set method used in optimization. The knee bending prohibition was implemented to pre vent the constraint update state ˜ x + k from resulting in pose estimates where the pelvis sinks tow ards the floor . Note that the knee can still bend during the prediction and measurement update. Mathematically , this is implemented by setting the constrained knee angle α 0 lk to Eq. (23) where α lk,min = 0 ◦ to prev ent knee hyperextension and α lk,max = min (180 ◦ , ˆ α + lk ) . The knee angle α lk is calculated by taking the inv erse tangent of the thigh vector , r lt z , projected on the z and x axes of the shank orientation as shown in Eq. (24). It ranges from − π 2 to 3 π 2 as enforced by the chosen sign inside the tangent in verse function and the addition of π 2 . Rearranging Eq. (24) and setting α lk = α 0 lk in terms of D lkr,k ˜ x + k = d ltr,k as sho wn below gives us Eqs. (27) and (28). Note that rearranging Eq. (25) and substituting r lt z = τ lt z ( ˜ x + k ) || τ lt z ( ˜ x + k ) || giv es Eq. (26). α 0 lk = min ( α lk,max , max ( α lk,min , α lk )) (23) α lk = tan − 1  − r lt z · r ls z − r lt z · r ls x  + π 2 (24) − r lt z · r ls z − r lt z · r ls x = sin( α 0 lk − π 2 ) cos( α 0 lk − π 2 ) (25) τ lt z · Denoted as ψ z }| { ( r ls z cos( α 0 lk − π 2 ) − r ls x sin( α 0 lk − π 2 )) = 0 (26) D lkr,k =  ψ − ψ 0 1 × 12  (27) d lkr,k = ψ · ( − d p 2 r p y + d ls r ls z ) (28) In summary , the linearized constraint equation D k ˜ x + k = d k has D k =  D T L,k D T R,k  T and d k =  d T L,k d T R,k  T . D L,k is shown in Eq. (29) where a bounded knee angle, α lk , is defined as α lk,min < α lk < α lk,max . D R,k , d L,k and d R,k can be deriv ed similarly . D L,k = ( [ D T ltl,k D T lkh,k ] T bounded α lk , [ D T ltl,k D T lkh,k D T lkr,k ] T otherwise. (29) The non-linear constraints for both legs as defined in D k ˜ x + k = d k are then enforced through an iterativ e projection scheme based on the Smoothly Constrained KF (SCKF) [29]. SCKF applies the standard KF measurement update procedure iterativ ely , starting with a constraint cov ariance, ˜ P + k , equal to the a posteriori covariance, P + k . The iteration stops when the termination criteria for all constraints are satisfied. The termination criterion for each constraint is shown in Eq. (30) where D i,j,k and ˜ P + j,j,k denotes the scalar value at the i th row and j th column of D k and ˜ P + k , respectiv ely; and D i,k denotes the i th row of D k . Please refer to [29] for a detailed description of SCKF . Note that SCKF can be replaced with other constrained optimization algorithms ((e.g., FMINCON [30]). The goal is to satisfy the biomechanical constraints before the next time step, by whatev er means. max ( D i,j,k ˜ P + j,j,k D i,j,k ) D i,k ˜ P + k D T i,k ≥ s thresh,k (30) F . P ost-pr ocessing The orientation of the pelvis and shank are obtained from the state ˜ x + k . The orientation of the left thigh, R lt , can be cal- culated using R lt =  r ls y × r lt z r ls y r lt z  where r lt z = τ lt || τ lt || . The orientation of the right thigh, R rt , is calculated similarly . I I I . E X P E R I M E N T The goal of the e xperiment is to ev aluate the algorithm described in this paper against benchmark systems, focusing only on the kinematics of the pelvis, thighs, and shanks. IEEE TRANSA CTIONS 6 A. Setup The experiment had nine healthy subjects ( 7 men and 2 women, weight 63 . 0 ± 6 . 8 kg, height 1 . 70 ± 0 . 06 m, age 24 . 6 ± 3 . 9 years old), with no kno wn gait or lower body biome- chanical abnormalities. Our system was compared to two benchmark systems, namely an OMC (i.e., V icon) and OSPS (i.e., Xsens) systems. The V icon V antage system consisted of eight cameras covering approximately 4 × 4 m 2 capture area with millimetre accuracy . V icon data were captured at 100 Hz and processed using Nexus 2.7 software. The Xsens A winda system consisted of sev en MTx units (IMUs). Xsens data were captured at 100 Hz using MT Manager 4.8 and processed using MVN Studio 4.4 software. The V icon and Xsens recordings were synchronized by having the Xsens A winda station send a trigger pulse to the V icon system at the start and stop ev ent of each recording [31]. Each subject had reflectiv e V icon markers placed according to the Helen-Hayes 16 marker set [32], sev en MTx units attached to the pelvis, thighs, shanks, and feet according to standard Xsens sensor placement, and two MTx units attached near the ankles, as shown in Fig. 3. The MTx units were all factory calibrated. Note that there are many published works that hav e validated the Xsens system against OMC systems [4, 33]. Fig. 3. Sample front and back vie w of the OMC (i.e., V icon) reflectiv e mark ers and the OSPC (i.e., Xsens) IMU setup used in validation experiments. Each subject performed the movements listed in T able I twice (i.e., two trials). The subjects stood still before and after each trial for ten seconds. The experiment was approved by the Human Research Ethics Board of the Uni versity of New South W ales (UNSW) with approv al number HC180413. T ABLE I T Y PE S O F M OVE M E N TS D O N E I N T H E V A LI D A T I ON E X P E RI M E NT Movement Description Duration (s) Static Stand still ∼ 10 W alk W alk straight and back ∼ 30 Figure of eight W alk in figures of eight ∼ 60 Zig-zag W alk zigzag ∼ 60 5-minute walk Undirected walk, side step, and stand ∼ 300 B. Calibration 1) F rame alignment: A common reference frame is needed for objecti ve pose comparison. In this experiment, comparison was done in the world frame, W , which was set equal to the Xsens frame. The pose estimate from the OMC system which was in the vicon frame, V , was con verted to the world frame, W , using a fixed rotational offset, q V 0 , which was calculated from a calibration procedure where the direction of gra vity and the Earth’ s magnetic field were measured using a compass and pendulum with optical markers attached [34]. 2) Y aw offset: Due to v arying magnetic interference be- tween the ankle and pelvis, an additional calibration procedure was executed to ensure that the x axes (north) of all sensors were aligned. The sensors in the ankle were especially af- fected, with yaw offset (i.e., rotation offset around the z axes) errors due to this magnetic disturbance. A grid search was performed to find the yaw offset that produces the least x and y axis RMSE between the OMC inferred acceleration and the IMU acceleration measurement after adjusting for each candidate yaw offset (as shown in Eq. (31) for the left ankle where ˚ τ =  0 τ T  T , ⊗ and − 1 are the quaternion and in version operator, respectiv ely). The optimal yaw offset was then used to adjust the corresponding shank orientations. Note that the calibration was done against the OMC benchmark system to minimize pose differences due to calibration and ensure an objectiv e ev aluation of the CKF algorithm. RMSE x,y ( ˚ a la omc − ( q ya w ⊗ ˚ ˘ a la imu ⊗ q − 1 ya w )) (31) C. System parameters The proposed algorithm used measurements from the three MTx units attached to the pelvis and ankles. The algorithm and calculations were implemented using Matlab 2018b. 1) Body se gment orientation: Since our goal is to ev aluate the tracking algorithm, not the initialization procedure, the initial pose of the body segments, W q B 0 , was obtained from the OMC system where that subject was asked to stand in an N pose (standing upright, with hands by sides). The measured orientations, IO , were obtained using the orientation estimates provided by the Xsens system [9]. 2) Step detection: The algorithm used for this experiment detects a step if the variance of the free body acceleration across the 0 . 25 second windo w is belo w a threshold. The threshold was set to 1 m.s − 2 . The step detection sometimes detect steps even if the foot/ankle was moving. T o prev ent such false positi ves, the step detection output was manually revie wed and corrected. 3) State estimator: The position and velocity of ˜ x + 0 were set to the initial position and velocity obtained from the OMC benchmark system. P + 0 was set to 0 . 5 I 30 × 30 . The SCKF threshold, s thresh,k , was set to 100 . The variance parameters used to generate the process and measurement error covariance matrix Q and R are shown in T able II. T ABLE II V A R IA N C E P A RA M E T ER S F O R G E N E RAT IN G T H E P RO C E S S A N D M E AS U R E ME N T E RR OR C OV A RI A N C E M ATR I C E S , Q A N D R . Q Parameters R Parameters σ 2 acc σ 2 mp σ 2 ls and σ 2 rs σ 2 lim (m 2 .s − 4 ) (m 2 ) ([m 2 .s − 2 m 2 ]) (m 2 ) 10 2 1 9 [10 2 10 2 0 . 1] [0 . 01 1 3 10 − 4 ] 10 2 1 9 where 1 n is an 1 × n ro w vector with all elements equal to 1 D. Evaluation Metrics In order to compare the algorithm’ s accuracy against an existing motion capture system that uses fewer sensors, the following metrics were chosen. IEEE TRANSA CTIONS 7 1) Mean position RMSE: is a common metric in video based human motion capture systems (e.g., [20]), commonly denoted as d pos . In this paper , it is denoted as e pos to pre vent confusion with body segment length symbols. Howe ver , as our system only captures the lower body , the metric has been mod- ified as shown in Eq. (32) where DP = { lh, r h, l k , r k, l a, r a } , p i k is the benchmark position, and ˜ p i k is the estimated position. Six body points, N pos = 6 , were ev aluated instead of thirteen body points (full body). Furthermore, as the global position of the estimate is still prone to drift due to the absence of an external global position reference, the root position of our system was set equal to that of the benchmark system (i.e., the mid-pelvis is placed at the origin for all RMSE calculations). e pos,k = 1 N pos P i ∈ DP || p i k − ˜ p i k || (32) 2) Mean orientation RMSE: e ori are the exponential co- ordinates of the rotational offset between the estimated and benchmark orientation of the body segments which do not hav e IMUs attached to them (i.e., the thighs, DO = { q lt k q rt k } and N ori = 2 ). This metric was calculated similarly to [20], as shown in Eq. (33), where q i k is the benchmark orienta- tion, ˜ q i k is the estimated orientation, and v-operator extracts the coordinates of the ske w-symmetric matrix obtained from the matrix logarithms operation. In practical implementation, the quaternion inside the logarithm function is con verted to rotation matrix representation. e ori,k = 1 N ori P i ∈ DO || log  q i k ⊗ ( ˜ q i k ) − 1  v || (33) 3) Hip and knee joint angles RMSE and coefficient of corr elation (CC): The hip joint angles in the sagittal, frontal, and transverse plane and knee joint angle in the sagittal plane are commonly used parameters in gait analysis. Hence, the angle RMSE and CC are useful indicators of the accuracy of the system in clinical applications. CC of joint i is obtained using Eq. (34), where θ i k and ˜ θ i k are the benchmark angle and the estimated angle of joint i at time step k , respectiv ely . CC i = n P θ i k ˜ θ i k − ( P θ i k )( P ˜ θ i k ) q n ( P ( θ i k ) 2 ) − ( P θ i k ) 2 q n ( P ( ˜ θ i k ) 2 ) − ( P ˜ θ i k ) 2 (34) 4) T otal travelled distance (TTD) deviation: is the TTD error with respect to the actual TTD. TTD is calculated by summing the distance between the XY positions of the pelvis and ankles at specific event times. The ev ents for the pelvis and left ankle are when left footsteps are detected, while the ev ents for the right ankle are when right footsteps are detected. I V . R E S U L T S A. Mean position and orientation RMSE Fig. 4 shows the mean position and orientation RMSE of our algorithm for each mov ement type. The comparison in volved the output from three configurations of interest, all compared against the OMC gold standard output: i) our algorithm using acceleration (double deriv ati ve of position) and orientation obtained from OMC (denoted as CKF-OMC ); ii) our algorithm using raw 3D acceleration in sensor frame and orientation in world frame obtained from Xsens MTx IMU measurements (denoted as CKF-3IMU ), and; iii) the black box output from the MVN Studio software (denoted as OSPS ). CKF-OMC gives insight into the best possible output of the CKF algorithm when given acceleration and orientation which corresponds perfectly with the position measured by the OMC benchmark system, while the OSPS illustrates the performance of a widely accepted commercial wearable HMCS with an OSPS configuration. Both biased and unbiased (i.e., for unbi- ased, the mean difference between the angles ov er each entire trial was subtracted) e ori are presented to account for different anatomical calibration offset errors between the OMC and OSPS systems [4, 33]. Fig. 4. The mean position and orientation RMSE, e pos (top) and e ori (bottom), for CKF-OMC , CKF-3IMU , and OSPS of each motion type. The prefix b denotes biased, while nb denotes no or without bias. B. Hip and knee joint angle RMSE and CC Fig. 5 shows the knee and hip joint angle RMSE and CC for CKF-3IMU . Y , X, and Z refers to the plane defined by the normal vectors y , x , and z axes, respectively , and are also known as the sagittal, frontal, and transversal plane in the context of gait analysis. Fig. 6 sho ws a sample W alk trial. Fig. 5. The joint angle RMSE (top) and CC (bottom) of knee and hip joint angles for CKF-3IMU at each motion type. Y , X, and Z denotes the sagittal, frontal, and transv ersal plane, respecti vely . The prefix b denotes biased, while nb denotes no bias. C. T otal travelled distance (TTD) deviation T able III shows the TTD deviation for CKF-3IMU at the pelvis and ankles along with related literature [35, 36]. Refer to the supplementary material for a tvideo reconstruction of sample trials [37]. V . D I S C U S S I O N In this paper , a CKF algorithm for lower body pose esti- mation using only three IMUs, ergonomically placed on the IEEE TRANSA CTIONS 8 Fig. 6. Knee and hip joint angle output of CKF-3IMU in comparison with OSPS and the benchmark system (OMC) for a W alk trial. The subject walked straight from t = 0 to 3 s, turned 180 ◦ around from t = 3 to 5 . 5 s, and walked straight to original point from 5 . 5 s until the end of the trial. T ABLE III T OT A L TR A V E L L ED D I S T A NC E D E VI ATI O N F OR CKF-3IMU AT T H E P ELV IS A N D A N K L ES A L O NG W I T H R E L A T E D L IT E R A T U R E Pelvis Left Ankle Right Ankle [35] [36] 4.93% 3.81% 3.60% 0.3-1.5% 0.21-0.35% ankles and sacrum to facilitate continuous recording outside the laboratory , was described and ev aluated. The algorithm utilizes fewer sensors than other approaches reported in the literature, at the cost of reduced accuracy . A. Mean position and orientation RMSE comparison The mean position and orientation RMSE, e pos and e ori , of CKF-OMC , CKF-3IMU , OSPS , and related literature for walking related motions (i.e., W alk , F igure of eight , Zig-zag , and 5-minute walk ) are shown in T able IV [20]. The e pos and e ori (no bias) difference between CKF-OMC and CKF-3IMU were around 1 cm and 5 ◦ , respectiv ely , indicating that our estimator w as able to handle sensor noise well when taking ori- entation and acceleration readings from the IMUs, rather than being deri ved from the OMC position measurements. Note that the relati vely large orientation error for CKF-OMC was caused by the constant body segment length and hinge knee joint assumption; measuring body segments as the distance between adjacent joint centre positions, as estimated by the V icon software, results in segment lengths which vary slightly throughout the gait cycle. The OSPS errors were comparable to the results of CKF-OMC algorithm. There was significant improv ement in the e ori results of the OSPS when bias was remov ed from this error variable; this bias probably due to the difference in the sensor-to-body-segment rotational offset, B q S 0 , caused by OSPS’ (i.e., Xsens’) black-box pre-processing of body segment orientation (Sec. II-C1) [4, 33]. CKF-3IMU ’ s performance for motion 5-minute walk in comparison to other shorter duration mov ements shows its tracking capability even ov er long durations; note all errors are measured locally , with the mid-pelvis as the origin. CKF-3IMU was compared to two algorithms from the study of Marcard et al . [20]; namely , sparse orientation poser ( SOP ) and sparse inertial poser ( SIP ). Note that SOP used orientation measured by IMUs and biomechanical constraints, while SIP used similar information but with the addition of acceleration measured by IMUs. Both SOP and SIP were benchmarked against an OSPS system while our algorithm was benchmark ed against an OMC system. The e pos and e ori (no bias) dif fer- ence between SOP and CKF-3IMU were ~ 0 . 2 cm and 1 ◦ , respectiv ely (i.e., comparable). Howe ver , the e pos and e ori (no bias) of SIP were smaller than CKF-3IMU by ~ 2 . 2 cm and 5 ◦ , respectiv ely . This impro vement was expected, as SIP optimizes the pose over multiple frames whereas CKF-3IMU optimizes the pose for each individual frame. Comparing processing times, CKF-3IMU was dramatically faster than SIP ; CKF- 3IMU processed a 1,000-frame sequence in less than a second on a Intel Core i5-6500 3.2 GHz CPU, while SIP took 7.5 minutes on a quad-core Intel Core i7 3.5 GHz CPU. Both set-ups used single-core non-optimized MA TLAB code. The CKF-3IMU could hence be used as part of a wearable HMCS which could provide real-time gait parameter measurement to inform actuation of assistive or rehabilitation robotic devices. T ABLE IV M E AN P O S IT I O N AN D O R I EN TA T I O N R M S E O F CKF-OMC , CKF-3IMU , OSPS , S PA RS E O R IE N TA T I ON A N D I NE RT I A L P O S E R ( SOP A N D SIP ) [ 2 0 ] CKF-OMC CKF-3IMU OSPS SOP SIP e pos (cm) 3 . 9 ± 0 . 8 5 . 2 ± 1 . 3 6 . 3 ± 1 . 8 ~ 5 . 0 ~ 3 . 0 e ori ( ◦ ) no bias 11 . 2 ± 1 . 6 16 . 1 ± 3 . 2 12 . 9 ± 4 . 0 ~ 15 . 0 ~ 11 . 0 mean* ( ◦ ) 0 . 6 ± 1 . 4 1 . 2 ± 1 . 9 − 2 . 5 ± 6 . 7 - - e ori ( ◦ ) biased 13 . 4 ± 3 . 6 18 . 9 ± 5 . 3 23 . 0 ± 12 . 0 - - * aggregated on per-trial level B. Hip and knee joint angle RMSE and CC comparison The knee and hip joint angle RMSEs and CCs of CKF- 3IMU , OSPS , and related literature are shown in T able V [19, 33, 38]. The RMSE dif ference ( < 6 ◦ ) between CKF-3IMU and OSPS indicates that our estimator is robust to handle sensor noise. The CC of CKF-3IMU , specifically the knee and hip joint angles in the sagittal plane ( 0 . 87 and 0 . 74 ), indicates that these angles hav e good correlation with the OMC benchmark system . Similar qualitative observations can be made in Fig. 6. The dif ference between CKF-3IMU and the OMC benchmark on the frontal and transverse plane were probably due to the strict hinged knee joint assumption, which is not perfectly anatomically accurate; the knee has a moving center of rotation more accurately approximated as a rolling contact between the shank and thigh segments. The hinge joint does not capture the moving center of rotation and contributes to estimation error . V iolations of this assumption are apparent at t = 3 to 5 . 5 s in Fig. 6, where the subject makes a 180 ◦ turn. After the increased angular error due to the 180 ◦ turn, CKF-3IMU was able to recover during the straight walking ( t = 5 . 5 to 9 . 74 s of Fig. 6). Notice that the bias between OSPS and OMC can be observed in Fig. 6 at t = 0 of the hip angle waveforms. The results also agree with the study of Cloete et al ., namely the smaller CC and larger RMSEs (both biased and no bias) for hip joint angles in the frontal and transverse planes [33]. IEEE TRANSA CTIONS 9 Considering the study of Hu et al . and T adano et al ., the errors were 1 ◦ to 6 ◦ smaller than CKF-3IMU [19, 38]. Hu et al . used 4 IMUs (two at pelvis and one for each foot) and T adano et al . used an OSPS configuration. Ho wever , their trials were limited to a single gait cycle, in contrast to the longer - duration walking trials in this paper which had at least fiv e gait cycles inv olving turns and capture duration of 1 to 5 minutes. Lastly , it is important to note that although CKF-3IMU achiev es good joint angle CCs in the sagittal plane, the joint angle RMSE ( > 5 ◦ ) makes its utility in clinical applications uncertain. T o achieve clinical utility , one may lev erage long- term recordings by averaging out cycle-to-c ycle variation in estimation errors over many gait cycles, or consider only body mov ements from which kinematics that can be estimated more accurately . The extent to which these possible solutions can bridge the gap to clinical application for the proposed system remains to be seen in future work. T ABLE V K N EE A N D H IP A N G LE R M S E ( T O P ) A N D C C ( B OT T OM ) O F CKF-3IMU , OSPS , A N D R EL ATE D L I TE R A T U R E Joint Angle RMSE ( ◦ ) knee sagittal hip sagittal hip frontal hip transverse CKF- 3IMU biased 11 . 1 ± 3 . 0 11 . 8 ± 3 . 3 7 . 5 ± 3 . 1 17 . 5 ± 4 . 7 mean − 1 . 2 ± 4 . 2 − 4 . 3 ± 4 . 4 − 2 . 2 ± 4 . 2 − 4 . 0 ± 9 . 7 no bias 10 . 0 ± 2 . 9 9 . 9 ± 3 . 2 6 . 1 ± 1 . 8 13 . 9 ± 2 . 4 OSPS biased 7 . 9 ± 3 . 3 12 . 4 ± 6 . 0 6 . 2 ± 2 . 6 19 . 8 ± 6 . 6 mean 0 . 1 ± 6 . 1 − 10 . 9 ± 7 . 4 0 . 2 ± 2 . 5 8 . 8 ± 8 . 8 no bias 5 . 0 ± 1 . 8 3 . 6 ± 1 . 7 4 . 1 ± 2 . 2 11 . 9 ± 4 . 3 Cloete et al .[33] biased 11 . 5 ± 6 . 4 16 . 9 ± 3 . 6 9 . 6 ± 5 . 1 16 . 0 ± 8 . 8 no bias 8 . 5 ± 5 . 0 5 . 8 ± 3 . 8 7 . 3 ± 5 . 2 7 . 9 ± 4 . 9 Hu et al .[19] 4 . 9 ± 3 . 5 6 . 8 ± 3 . 0 - - T adano et al .[38] 10 . 1 ± 1 . 0 7 . 9 ± 1 . 0 - - Joint Angle CC knee sagittal hip sagittal hip frontal hip transverse CKF-3IMU 0 . 87 ± 0 . 08 0 . 74 ± 0 . 12 0 . 62 ± 0 . 12 0 . 33 ± 0 . 12 OSPS 0 . 97 ± 0 . 04 0 . 95 ± 0 . 06 0 . 72 ± 0 . 19 0 . 26 ± 0 . 20 Cloete et al .[33] 0 . 89 ± 0 . 15 0 . 94 ± 0 . 08 0 . 55 ± 0 . 40 0 . 54 ± 0 . 20 Hu et al .[19] 0 . 95 ± 0 . 04 0 . 97 ± 0 . 04 - - T adano et al .[38] 0 . 97 ± 0 . 02 0 . 98 ± 0 . 01 - - C. T otal travelled distance (TTD) deviation CKF-3IMU was able to estimate global position at 3 . 6 - 4 . 9 % TTD deviation, which is not as good as state-of-the-art dead reckoning applications [35, 36]; some of the reason for this also lies in the assumption that the velocity of the shank IMU is zero when the associated foot touches the floor , but of course this IMU continues to move with some small velocity during the stance phase. Howe ver , pelvis drift has been re- duced substantially compared to [20]. Note that the experiment was performed in a small room ( 4 × 4 m 2 ) to capture lower body pose (the main parameter the algorithm aims to solve). Note, dead reckoning algorithms are typically tested on longer courses (e.g., ~ 125 − 300 m long, sometimes with different floor levels), which must be consider when interpreting the global positioning results presented here. D. P elvis drift and numerical issues W ithout the pelvis position pseudo-measurements during the measurement update of the KF (Section II-E2) and without the non-decreasing knee angle during the constraint update (Section II-E3), the resulting pose estimate may drift into a crouching pose as shown in Fig. 7. (a) Front (b) Side Fig. 7. Snapshot of crouching pose estimate (pelvis drift) in approximate front and side perspective views. The red-green-blue and pink-yellow-c yan lines refer to CKF-3IMU estimated and ground truth poses, respectively . W ithout the cov ariance limiter during the measurement update (Section II-E2), the covariance of the CKF algorithm becomes badly conditioned for the motion 5-minute walk as the position uncertainty grows at a faster rate for the pelvis, as the ankle velocities are zeroed at the time of each footstep, which subsequently limits the rate at which the positional uncertainly of the ankles grows. Expectantly , this was not an issue for the shorter walking motions. As part of future w ork, methods to update the state er- ror covariance matrix using information obtained during the constraint update will be explored. In the current frame work, the update of the a posteriori state error covariance matrix, P + k , during the constraint step expectedly resulted in rank- deficiency (i.e., singular) due to the eigenv ectors of this cov ariance matrix being projected onto the constraint surface, which spans a subspace of the state space. This rank-deficienc y ultimately resulted in pose estimates that are physically impos- sible or numerically undefined. Furthermore, methods which can better track states and covariances of non-linear systems (e.g., body pose) will also be explored. Indeed, this weakness is a well-known limitation of the KF when applied to non- linear system models, but can be better addressed using an unscented KF or a particle filter [21, 39]. Finally , in addition to position and velocity , orientation of the body segments will also be included as state v ariables to allow the KF to make corrections to segment orientations, compensating for errors in the IMU orientation measurements. E. CKF limitations: towar ds long-term trac king of ADL The assumptions in the measurement and constraint update of the CKF algorithm impose limitations to its general usabil- ity . The pelvis position pseudo-measurements used during the measurement update helps prev ent the pelvis drifting away from the feet. Howe ver , it will also pre vent accurate pose estimation of motions such as sitting and lying down, where the pose is maintained for a duration much longer than that of a typical gait cycle. As part of future work, the constraint may be replaced by balance-based constraints that take account of pelvis contact with an external object (e.g., detect if pelvis has contact with a chair or bed). The flat-floor assumption during the measurement update limits the system to environments IEEE TRANSA CTIONS 10 were the altitude is not changing and is unable to track situations such as walking up ramps or climbing stairs. Lastly , the sensors are assumed rigidly attached to corresponding body segments, whereas in reality sensors are attached using tight velcro straps which may mov e during dynamic movements like jumping or running. The CKF algorithm is expected to work on pathological gait where the biomechanical assumptions of the constraint model are not violated. For example, subjects with knee osteoarthritis hav e demonstrated reduced knee excursion (peak knee flexion angle minus angle at initial contact in degrees) in the sagittal plane [40]. The CKF should be able to capture such motion giv en the subject has no varus or valgus deformity (hinge knee joint assumption), e ven if it includes an atypical spatiotempo- ral gait pattern; for example, intermittent stopping during the 5-minute walking trial. Howe ver , this paper provides an initial proof-of-concept ev aluation of the proposed algorithm using an RSC configuration on participants with normal gait and lower body biomechanics. Future work validating the accuracy of this system in a larger and more diverse cohort would help quantify its ultimate clinical utility . The CKF algorithm was also highly dependent on the step detection algorithm and the initial calibration step which defines the relative orientation between each IMU and its body segment. A step detection error can cause a large velocity and displacement correction on the same ankle at the next detected step which adds to the global position drift. Of course, in everyday life, the OMC reference orientation and acceleration will not be av ailable for the sensor-to-body orientation alignment and yaw offset alignment. A practical initial or e ven online calibration procedure that takes into account magnetic interference at the ankle sensors close to the floor will therefore be needed. There are different ways of sensor-to-segment calibration in existing literature ranging from static to online/adaptiv e calibration. For example, the subject may be asked to do a predefined posture at the start of pose estimation for orientation alignment [9, 41]. In addition, the subject may be asked to walk in a straight line and then back to the starting point for yaw offset alignment. The global yaw correction resulting in the highest y axis (sagittal plane) range of motion may be assumed to be aligned with the pelvis anteroposterior axis. The magnetic field is kno wn to be inhomogeneous in indoor environments, typically with stronger disturbances closer to the floor [42]. Our algorithm would benefit from orientation estimation algorithms which are rob ust to magnetic disturbances, such as those presented by Del Rosario et al. [13]. Unfortunately , recent adv ances in automatic sensor-to-se gment calibration and magnetometer correction that utilizes two IMUs and a hinge joint may not be applicable to RSC configurations (e.g., our algorithm only has one IMU on one side of the hinge joint) [43]. V I . C O N C L U S I O N This paper presented a CKF-based algorithm to estimate lower limb kinematics using a reduced sensor count config- uration, and without using any motion database. The mean position and orientation (no bias) RMSEs (relati ve to the mid- pelvis origin) were 5 . 21 ± 1 . 3 cm and 16 . 1 ± 3 . 2 ◦ , respectively . The sagittal knee and hip angle RMSEs (no bias) were 10 . 0 ± 2 . 9 ◦ and 9 . 9 ± 3 . 2 ◦ , respectiv ely , while the CCs were 0 . 87 ± 0 . 08 and 0 . 74 ± 0 . 12 . The performance is acceptable for animation applications. Howe ver , accuracy needs to improv e for clinical applications. T o improve performance, additional information is needed to prev ent pelvis drift (e.g., utilize sensors that giv e pelvis distance, position, or angle of arri val relativ e to the ankle). Lastly , the paper contrib utes a new OMC (i.e., V icon) and IMU dataset which contains longer duration walking and will be made av ailable at https://git.io/Je9VV. V I I . A C K N O W L E D G E M E N T The authors would like to thank Neuroscience Research Australia (NeuRA) for the technical support and access to their Gait laboratory; Lucy Armitage for the crash course on anatomical landmarks. Michael Raitor was supported through a Fulbright scholarship. This research was supported by an Australian Government Research T raining Program (R TP) Scholarship. The source code for the CKF algorithm and a sample video will be made av ailable at https://git.io/Je9VV. R E F E R E N C E S [1] P . B. Shull et al. , “Quantified self and human movement: A revie w on the clinical impact of wearable sensing and feedback for gait analysis and intervention, ” Gait P osture , vol. 40, no. 1, pp. 11–19, 2014. [2] P . Ornetti et al. , “Gait analysis as a quantifiable outcome measure in hip or knee osteoarthritis: a systematic review, ” Jt. Bone Spine , vol. 77, no. 5, pp. 421–425, 2010. [3] J. M. Hausdorff et al. , “Gait variability and fall risk in community- living older adults: A 1-year prospectiv e study, ” Ar ch. Phys. Med. Rehabil. , vol. 82, no. 8, pp. 1050–1056, 2001. [4] J. C. V an Den Noort et al. , “Gait analysis in children with cerebral palsy via inertial and magnetic sensors, ” Med. Biol. Eng. Comput. , v ol. 51, no. 4, pp. 377–386, 2013. [5] P . Shull et al. , “Haptic gait retraining for knee osteoarthritis treatment, ” in 2010 IEEE Haptics Symp. , IEEE, 2010, pp. 409–416. [6] J. B. Dingwell et al. , “Local dynamic stability versus kinematic vari- ability of continuous overground and treadmill walking, ” J. Biomech. Eng. , vol. 123, no. 1, p. 27, 2001. [7] S. T edesco et al. , “A review of activity trackers for senior citizens: Research perspectives, commercial landscape and the role of the insurance industry , ” Sensors , vol. 17, no. 6, p. 1277, 2017. [8] C. Tudor -Locke et al. , “Peak stepping cadence in free-living adults: 2005-2006 NHANES., ” J. Phys. Act. Health , vol. 9, no. 8, pp. 1125–9, 2012. [9] D. Roetenber g et al. , “Xsens MVN: Full 6DOF human motion tracking using miniature inertial sensors, ” Xsens Motion T echnol. BV , T ech. Rep , vol. 1, 2009. [10] BioSenics, LEGSys , 2018. [Online]. A vailable: http://www .biosensics. com/products/legsys/ (visited on 06/27/2018). [11] D. Laidig et al. , “Automatic anatomical calibration for IMU-based elbow angle measurement in disturbed magnetic fields, ” Curr . Dir . Biomed. Eng. , vol. 3, no. 2, pp. 167–170, 2017. [12] M. B. Del Rosario et al. , “Quaternion-based complementary filter for attitude determination of a smartphone, ” IEEE Sens. J. , vol. 16, no. 15, pp. 6008–6017, 2016. [13] M. B. Del Rosario et al. , “Computationally efficient adaptive error- state Kalman filter for attitude estimation, ” IEEE Sens. J . , vol. 18, no. 22, pp. 9332–9342, 2018. [14] J. L. McGinley et al. , “The reliability of three-dimensional kinematic gait measurements: A systematic revie w, ” Gait P osture , vol. 29, no. 3, pp. 360–369, 2009. [15] M. Al-Amri et al. , “Inertial measurement units for clinical mov ement analysis: Reliability and concurrent validity, ” Sensors (Switzerland) , vol. 18, no. 3, p. 719, 2018. [16] J. T autges et al. , “Motion reconstruction using sparse accelerometer data, ” A CM Tr ans. Graph. , vol. 30, no. 3, p. 18, 2011. eprint: 1006 . 4903. IEEE TRANSA CTIONS 11 [17] F . F . W ouda et al. , “Estimation of full-body poses using only five inertial sensors: an eager or lazy learning approach?” Sensors , vol. 16, no. 12, p. 2138, 2016. [18] Y . Huang et al. , “Deep inertial poser: Learning to reconstruct human pose from sparse inertial measurements in real time, ” in SIGGRAPH Asia 2018 T ech. P ap. SIGGRAPH Asia 2018 , Association for Comput- ing Machinery , Inc, 2018. eprint: 1810.04703. [19] X. Hu et al. , “Performance ev aluation of lower limb ambulatory measurement using reduced Inertial Measurement Units and 3R gait model, ” in IEEE Int. Conf. Rehabil. Robot. , 2015, pp. 549–554. [20] T . von Marcard et al. , “Sparse inertial poser: Automatic 3d human pose estimation from sparse IMUs, ” in Comput. Graph. F orum , W iley Online Library, v ol. 36, 2017, pp. 349–360. eprint: 1703.08014. [21] D. Simon, Optimal state estimation: Kalman, H infinity , and nonlinear appr oaches . John W iley & Sons, 2006. [22] J. B. Kuipers, Quaternions and r otation sequences : a primer with ap- plications to orbits, aer ospace, and virtual r eality . Princeton Uni versity Press, 1999, p. 371. [23] M. Kok et al. , “Using inertial sensors for position and orientation estimation, ” F ound. T rends Signal Pr ocess. , vol. 11, no. 1-2, pp. 1–153, 2017. eprint: arXi v:1704.06053v2. [24] O. J. W oodman, “An introduction to inertial na vigation, ” Uni versity of Cambridge, Computer Laboratory, T ech. Rep., 2007. [25] I. Skog et al. , “Zero-velocity detection - an algorithm ev aluation., ” IEEE T rans. Biomed. Eng. , vol. 57, no. 11, pp. 2657–2666, 2010. [26] R. E. Kalman, “A new approach to linear filtering and prediction problems, ” J. Basic Eng. , v ol. 82, no. 1, pp. 35–45, 1960. [27] D. J. Simon et al. , “Kalman filtering with state constraints: a survey of linear and nonlinear algorithms, ” IET Contr ol Theory Appl. , vol. 4, no. 8, pp. 1303–1318, 2010. [28] X. L. Meng et al. , “Biomechanical model-based displacement estima- tion in micro-sensor motion capture, ” Meas. Sci. T echnol. , vol. 23, no. 5, p. 055 101, 2012. [29] J. De Geeter et al. , “A smoothly constrained Kalman filter, ” IEEE T rans. P attern Anal. Mach. Intell. , vol. 19, no. 10, pp. 1171–1177, 1997. [30] Matlab, F ind minimum of constrained nonlinear multivariable function - MATLAB fmincon , 2018. [Online]. A v ailable: https :/ /au. mathworks. com/help/optim/ug/fmincon.html (visited on 03/26/2018). [31] Xsens, Syncr onising Xsens Systems with V icon Nexus , 2015. [On- line]. A v ailable: https : / / www . xsens . com / download / pdf / SynchronisingXsenswithV icon.pdf (visited on 05/07/2019). [32] M. P . Kadaba et al. , “Measurement of Lower Extremity Kinematics During Level W alking, ” J . Orthop. Res. , v ol. 8383, no. 3, pp. 383–392, 1990. [33] T . Cloete and C. Scheffer , “Benchmarking of a full-body inertial motion capture system for clinical gait analysis, ” in 2008 30th Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. , IEEE, 2008, pp. 4579–4582. [34] S. O. H. Madgwick et al. , “An efficient orientation filter for inertial measurement units (IMUs) and magnetic angular rate and gravity (MARG) sensor arrays, ” Department of Mechanical Engineering, T ech. Rep., 2010, pp. 113–118. [35] A. R. Jimenez et al. , “Indoor Pedestrian navigation using an INS/EKF framew ork for yaw drift reduction and a foot-mounted IMU, ” Proc. 2010 7th W ork. P ositioning, Navig. Commun. WPNC’10 , no. April, pp. 135–143, 2010. [36] W . Zhang et al. , “A foot-mounted PDR System Based on IMU/EKF+HMM+ZUPT+ZAR U+HDR+compass algorithm, ” 2017 Int. Conf. Indoor P osition. Indoor Navig. IPIN 2017 , vol. 2017-Janua, no. September , pp. 1–5, 2017. [37] Gait-tech/gt.papers . [Online]. A vailable: https: // github.com /gait - tech / gt.papers/tree/master/%2Bckf2019 (visited on 12/17/2019). [38] S. T adano et al. , “Three dimensional gait analysis using wearable acceleration and gyro sensors based on quaternion calculations., ” Sensors , v ol. 13, no. 7, pp. 9321–9343, 2013. eprint: NIHMS150003. [39] S. Julier et al. , “A new method for the nonlinear transformation of means and covariances in filters and estimators, ” IEEE T rans. Automat. Contr . , vol. 45, no. 3, pp. 477–482, 2000. [40] J. D. Childs et al. , “Alterations in lower extremity movement and muscle activ ation patterns in individuals with knee osteoarthritis, ” Clin. Biomech. , vol. 19, no. 1, pp. 44–49, 2004. [41] X. Meng et al. , “Hierarchical information fusion for global displace- ment estimation in microsensor motion capture, ” IEEE T rans. Biomed. Eng. , vol. 60, no. 7, pp. 2052–2063, 2013. [42] W . H. de Vries et al. , “Magnetic distortion in motion labs, implications for validating inertial magnetic sensors, ” Gait P ostur e , vol. 29, no. 4, pp. 535–541, 2009. [43] D. Laidig et al. , “Exploiting kinematic constraints to compensate magnetic disturbances when calculating joint angles of approximate hinge joints from orientation estimates of inertial sensors, ” IEEE Int. Conf. Rehabil. Robot. , pp. 971–976, 2017.