Drift-Free 3D Orientation and Displacement Estimation for Quasi-Cyclical Movements Using One Inertial Measurement Unit: Application to Running

A Drift-Free 3D Orientation and Displacement estimation method (DFOD) based on a single inertial measurement unit (IMU) is proposed and validated. Typically, body segment orientation and displacement methods rely on a constant- or zero-velocity point to correct for drift. Therefore, they are not easily applicable to more proximal segments than the foot. DFOD uses an alternative single sensor drift reduction strategy based on the quasi-cyclical nature of many human movements. DFOD assumes that the quasi-cyclical movement occurs in a quasi-2D plane and with an approximately constant cycle average velocity. DFOD is independent of a constant- or zero-velocity point, a biomechanical model, Kalman filtering or a magnetometer. DFOD reduces orientation drift by assuming a cyclical movement, and by defining a functional coordinate system with two functional axes. These axes are based on the mean acceleration and rotation axes over multiple complete gait cycles. Using this drift-free orientation estimate, the displacement of the sensor is computed by again assuming a cyclical movement. Drift in displacement is reduced by subtracting the mean value over five gait cycle from the free acceleration, velocity, and displacement. Estimated 3D sensor orientation and displacement for an IMU on the lower leg were validated with an optical motion capture system (OMCS) in four runners during constant velocity treadmill running. Root mean square errors for sensor orientation differences between DFOD and OMCS were 3.1 ± 0.4° (sagittal plane), 5.3 ± 1.1° (frontal plane), and 5.0 ± 2.1° (transversal plane). Sensor displacement differences had a root mean square error of 1.6 ± 0.2 cm (forward axis), 1.7 ± 0.6 cm (mediolateral axis), and 1.6 ± 0.2 cm (vertical axis). Hence, DFOD is a promising 3D drift-free orientation and displacement estimation method based on a single IMU in quasi-cyclical movements with many advantages over current methods.


Introduction
Activities like walking, running, swimming, rowing, and skating are all quasi-cyclical in nature. The repetitiveness of these movements, and their associated loads on the human body, can result in overuse injuries [1]. Repetitive movements are often studied inside movement analysis laboratories for insight into overloading of the human body and performance enhancement, among other applications. With the introduction of wearable systems, motion analysis is no longer restricted to a controlled lab setting [2,3]. This opens up new possibilities of analyzing movements that are difficult to measure in a lab, due to technical constraints of optical motion capture systems (OMCS).
Inertial measurement units (IMUs) are widely used in wearable motion capture systems due to their small size and ease of use [4]. IMUs are composed of accelerometers, gyroscopes, and are often combined with magnetometers. The acceleration, orientation, and displacement of a sensor are of interest for many motion analysis applications such as impact analyses, monitoring the range of motion (ROM), or inclination of a body segment, e.g., the lower leg [5,6]. To obtain orientation and displacement from sensor accelerations and angular velocities, the orientation of the sensor in the global coordinate system (CS) (Ψ gl ) is required. Displacement can then be obtained via strapdown inertial navigation, although this process is prone to errors [7]. Drift in estimated 3D orientations should be minimized as it strongly influences the estimated displacement in Ψ gl . Drift can be compensated for by incorporating other sensors (i.e., magnetometer). However, drift reduction and 3D orientation estimation become more challenging during highly dynamic movements, prolonged measurements, or when magnetic distortions are present [8]. Drift can alternatively be reduced by applying domain-specific assumptions, such as the zero-velocity update method [9,10].
The zero-velocity update method assumes that the velocity of the foot is zero and the orientation of the foot is known during the stance phase. This information is used to reset drift in orientation, velocity, and position during each gait cycle [9,10]. Similar assumptions have been used in running in specific conditions. Bailey and Harle corrected for positional drift of an IMU on the foot by using a constant-velocity update in runners with a heel strike [11]. However, constant-or zero-velocity assumptions are not suitable for more proximal segments or runners with a forefoot strike, since a constant-or zero-velocity point is often not present [12].
To estimate orientation and displacement using a single IMU placed on body segments without a constant-or zero-velocity point, the quasi-cyclical nature of numerous movements can be used. Kalman filtering or analytical integration of acceleration in combination with assumptions about the quasi-cyclical nature of movements have been used to estimate displacements of, for example, the pelvis during walking [13][14][15]. These studies involved relatively slow movements, required multiple sensors, a calibration procedure, or prior information about the movements. As a solution to many of these drawbacks, we propose to directly use the quasi-cyclical nature of movements to estimate 3D orientation and displacement using a single IMU without the need for Kalman filtering.
Hence, the research question of this study was: How to estimate 3D orientation and displacement of a single IMU on the lower leg using the quasi-cyclical nature of running?
To answer this research question, a method is proposed in which drift-free 3D orientation and displacement of a single IMU are estimated using the quasi-cyclical nature of numerous human movements. We call this method Drift-Free Orientation and Displacement estimation (DFOD). DFOD assumes that the movement is quasi-cyclical, occurs in a quasi-2D plane, and has an approximately constant cycle average velocity. DFOD will be demonstrated in treadmill running, although it is expected to generalize to many quasi-cyclical quasi-2D movements.

Materials and Methods
Validation of DFOD was part of a larger study. For sake of clarity, only measurement systems and trials required for validation of DFOD will be described.

Participants
Four healthy recreational runners participated in this study (2M/2F, age: 30.6 ± 9.2 years, height: 181 ± 4 cm, body mass: 65.0 ± 5.4 kg). The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of METC Twente. All participants gave written informed consent before participating in the study.

Protocol
Subjects ran for 2 min on a level treadmill at 3.6 m/s. To validate DFOD with OMCS, a calibration procedure was performed in which subjects stood still in a neutral pose and flexed and extended their leg four times while keeping their upper leg horizontal. This calibration procedure was not required for DFOD but was used to convert the OMCS orientation and displacement estimates to the same CS as used in DFOD for comparison purposes.

Measurement Systems
Subjects ran on a treadmill (C-Mill, ForceLink, Culemborg, The Netherlands) while 3D angular velocities and accelerations were captured by a single IMU on the lower leg at 240 Hz (MVN Link, Xsens, Enschede, The Netherlands). The ranges of the accelerometer and angular velocity sensor were ±16 g and ±2000 • /s, respectively. Positional data from a cluster marker set on the lower leg were captured for reference measurements with an eight-camera optical motion capture system (OMCS) at 100 Hz (Vantage, Vicon, Oxford, UK). The cluster marker set consisted of four individual markers attached to a rigid plate. The IMU was placed medially to the tibial tuberosity and the cluster marker set was placed below the IMU, both on the flat surface of the tibia to ensure measurements of tibia motion, see Figure 1. Both systems were attached to the skin with double-sided tape and covered with stretched strapping tape. Overview of IMU and cluster marker set placement. "M1","M2", and "M3" refer to the individual markers of the cluster marker set. The shown coordinate system is the functional coordinate system Ψ f . The X-axis points forward (running direction), the Y-axis mediolateral, and the Z-axis upward.

Data Preparation
Optical and inertial data of the left or right lower leg were selected based on minimal OMCS marker occlusion. Optical data were upsampled to 240 Hz with linear interpolation and low-pass filtered with a recursive fourth-order 20 Hz Butterworth filter [16]. Inertial data were not filtered.
Data were segmented into gait cycles based on falling edge angular velocity zerocrossings in the sensor CS (Ψ s ) Y-axis, which was directed in the global CS (Ψ gl ) forward/mediolateral direction. These zero-crossings occur shortly before initial contact. Data were cropped to include all complete gait cycles during one minute of running, hereby excluding around 30-45 s of data in which the subject increased their running speed from standing still up to 3.6 m/s. Sensor acceleration and angular velocity in Ψ s (i.e., input signals for DFOD) as a function of the time-normalized gait cycle for a representative subject are shown in Figure 2. Data analysis was performed in MATLAB R2021a.

Figure 2.
Three-dimensional sensor acceleration (i.e., including gravity) (left figure) and sensor angular velocity (right) in Ψ s as a function of the time-normalized gait cycle for a representative subject. Solid lines represent the mean while shaded areas represent the standard deviationaround the mean over one minute of running. Note that these two signals are the input for the orientation and displacement estimation algorithm. Positive acceleration values represent an acceleration into the upward, sideward (left), and forward direction of Ψ s . Positive angular velocity values represent anti-clockwise rotations in Ψ s .

Estimate Sensor Orientation in a Functional CS (Ψ f )
The aim of DFOD is to estimate 3D orientation and displacement of a single sensor in a functional CS (Ψ f ) of which the vertical and mediolateral axes are fixed and the origin moves with the body at the cycle average velocity. Ψ f is defined in Figure 1; Figure 3. DFOD assumes that the body segment on which the sensor is placed: • moves quasi-cyclical (i.e., cycles are similar) • moves in a quasi-2D plane (i.e., most movement occurs in a 2D plane) • has an approximately constant cycle average velocity The time-and gait cycle-dependent rotation matrix R f s,i (t) from Ψ s to Ψ f , representing the sensor orientation in a functional drift-free CS of which the vertical and mediolateral axes are fixed and the origin moves with the body at the cycle average velocity, can be written as three subsequent rotations as in Equation (1): where the sensor CS (Ψ s ) is first rotated to a sensor-fixed partly functional CS (Ψ p f ) with the time-independent rotation matrix from Ψ s to Ψ p f (R p f s ). Then, Ψ p f is rotated to a drifting partly functional CS (Ψ d f ) with the time-dependent rotation matrix from Ψ p f to Ψ d f (R  Columns represent different coordinate systems (CS). For each CS some basic information is stated: symbol, name, measurement system on which the CS is based, origin, fixed functional axes (i.e., which axes have functional meaning), and the presence of drift. Available quantities in each CS are shown in white squares (all time-dependent), rotation matrices from one CS to another are shown in blue arrows, curved arrows at the top represent the different rotations which are referred to in the text, integrations over time are shown as green arrows. At the bottom of the figure, a schematic representation of the CS with respect to the lower leg of a runner is shown, orange boxes represent the IMU, and blue dots represent the cluster marker set. DFOD is validated against an OMCS based on the quantities in the red squares. Note that the CSs in grey (two right columns) are only used for validation of DFOD and are not a part of DFOD. IMU = inertial measurement unit; OMCS = optical motion capture system; GCZM = gait cycle zero mean (mean value over each gait cycle is subtracted from the gait cycle); CS = coordinate system; Integration error accumulation can be reduced by aligning the rotation axis of a quasi-2D movement with one axis in 3D space to create a partly functional CS (Ψ p f ) [17]. One axis has functional and anatomical meaning in Ψ p f . The functional axis ( → y s p f , the y-axis of the sensor-fixed partly functional CS (Ψ p f ) expressed in the sensor CS (Ψ s )) is perpendicular to the plane of movement. Therefore, this axis is described by the first principal component of the angular velocity in Ψ s ( → ω s ), measured by the 3D angular velocity sensor of the IMU, over one minute of running [18]: To create a rotation matrix from Ψ p f to Ψ s , a temporary X-axis is defined by arbitrarily setting the X-axis of Ψ p f ( → x , s p f ) to the X-axis of Ψ s : The Z-axis of Ψ p f ( → z s p f ) was computed and → x , s p f updated to ensure an orthogonal CS according to the TRIAD algorithm [19]: The time-invariant orthonormal rotation matrix from Ψ p f to Ψ s was: The time-invariant rotation matrix from Ψ s to Ψ p f (R To go from a sensor-fixed CS to a drifting CS in which axes do not depend on the sensor orientation, the angular velocity in Ψ p f ( → ω p f ) was integrated according to Bortz [20,21].
→ ω p f was expressed as a skew-symmetric matrix ( ω p f ), and the differential equation was solved and used to obtain the rotation matrix from Note that Ψ d f drifts, predominantly around the y-axis (Ψ d f y ), due to accumulated integration errors from Equation (3a,b). This drift needs to be corrected to get a useful orientation estimate of the sensor (rotation 3 in Figure 3).

Rotation 3: From Drifting Partly Functional CS (Ψ d f ) to Drift-Free Functional CS (Ψ f )
Following an assumption of quasi-cyclical running, the lower leg keeps rotating around the same mediolateral axis. By continuously calculating this rotation axis we can correct for integration drift from Equation (3a,b). The rotation axis was again based on the first principal component of the 3D angular velocity, now in Ψ d f ( → ω d f ), over five complete gait cycles Equation (4a). Multiple gait cycles were included to obtain a more robust estimate of → y f d f ,i (see also Section 2.7 for algorithm characteristics): Following an assumption of approximately constant cycle average velocity running, the free acceleration in Ψ f will be approximately zero-mean over a complete gait cycle. Hence, the mean total acceleration (i.e., including gravity) over a complete number of gait cycles represents the gravitational acceleration and is directed vertically. Therefore, the temporary Z-axis of the functional CS Ψ f ( , f d f ,i updated to ensure an orthogonal CS according to the TRIAD algorithm [19]: The orthonormal drift-correcting rotation matrix from Ψ d f to Ψ f was: where R f d f ,i has a constant value within each cycle but varies over cycles to correct for drift. The drift-free 3D rotation matrix of the sensor in a functional CS (R f s,i ) of which the vertical and mediolateral axes are fixed, and the origin moves with the body at the cycle average velocity was then computed with Equation (1).

From Sensor Orientation to Sensor Displacement
Three-dimensional angular velocity and total (i.e., including gravity) acceleration obtained by subtracting the modulus of the gravitational acceleration ( Following an assumption of approximately constant cycle average velocity, the free acceleration in Ψ f will be approximately zero-mean over a complete number of gait cycles. Hence, the mean free acceleration value over a window of five gait cycle was subtracted to correct for drift: where → a f , f a,GZCM is the free acceleration with a gait cycle zero-mean (GCZM).
→ a f , f a,GZCM was numerically integrated (Figure 3, Column Ψ f , upper green arrow) with the trapezoidal rule to obtain the velocity ( Following an assumption of approximately constant cycle average velocity, the mean velocity in Ψ f over a complete number of gait cycles is approximately zero in all axes since the origin of Ψ f moves with the body at the cycle average velocity. Hence, the driftcorrected GCZM velocity was computed by subtracting the mean velocity over a window of five gait cycles: : Following an assumption of approximately constant cycle average velocity, the mean displacement in Ψ f approximates zero is all directions over a complete number of gait cycles. Hence, the mean displacement over five gait cycles was subtracted and GCZM displacement ( → s f ,GCZM ) was computed and used as outcome measure:

Validation of Orientation and Displacement Estimates
The steps described below are used to validate DFOD and are not part of DFOD. To compare the results of DFOD against OMCS, the orientation of the cluster marker set was computed and both the orientation and displacement were transformed to Ψ f .

Rotation 4: From Optical Motion
The orientation of the OMCS cluster marker set in Ψ gl was based on the relative positions of three of its individual markers according to the TRIAD algorithm [19]: where → p gl,m refers to the position of the individual markers of the cluster marker set in Ψ gl , see Figure 1.
The orthonormal rotation matrix of Ψ cl to Ψ gl (R gl cl ) was: To be comparable, the 3D orientation and position of the OMCS cluster marker set need to be expressed in the same functional CS as used for the sensor orientation and displacement estimate of DFOD. Therefore, the OMCS functional Y-axis (Ψ , T i being the duration of cycle i. See Section 2.7 for algorithm characteristics. By using a larger time interval, the change in rotation during this interval is relatively large compared to the errors. The rotation matrix from time point t j = t i + j × T i 7 to the next was then computed as follows: R Subsequently, R t j+1 t j of cycle i was transformed to a rotation axis ( → v rot,i,j ) which corresponds to the vector part of a quaternion that can be derived from a rotation matrix [22].
→ v rot,i,j was multiplied by a factor −1 for the extension part of each calibration movement cycle to ensure that the rotation axes were approximately equally directed for all intervals.
The functional coordinate axis → y gl f (i.e., the rotation axis of the lower leg during the flexionextension movements) was subsequently determined by averaging all resulting rectified rotation axes → v , rot,i,j for all intervals j and all cycles i: The temporary Z-axis of Ψ f ( The orthonormal time-invariant rotation matrix from Ψ f to Ψ gl (R gl f ) was: The time-dependent rotation matrix of Ψ cl in Ψ f (R f cl ) was then computed and represented the orientation of the cluster marker set in Ψ f (rotation 4 in Figure 3):

Orientation and Displacement Validation
The sensor and cluster orientation estimates in Ψ f of DFOD and OMCS were expressed in Euler angles (rotation order: YZX) for visualization purposes. To show the added driftreducing benefit of DFOD in estimating sensor orientation, sensor orientation was also computed by integrating the sensor angular velocity in Ψ s , similar to Equation (3a,b) without any drift reducing methods. This resulted in the sensor orientation with respect to the initial sensor orientation (Ψ s,init ) at the start of the first gait cycle. The position of the marker closest to the IMU was selected and displacement during each gait cycle was computed (similar procedure to Equation (7b)). OMCS and IMU data were timesynchronized based on the GCZM displacement in the forward direction of the sensor and cluster marker set ( → s f ,GCZM,x ). Three-dimensional differences in Euler angles and displacement between DFOD and OMCS over time-normalized gait cycles were quantified as root mean square errors (RMSE) and absolute mean differences. A 1D orientation error was computed by transforming the difference in orientation between DFOD and OMCS to an axis-angle representation and using the rotation angle as an outcome [23]. This 1D angle represents the rotation that is necessary to align R f s and R f cl . A 1D displacement error was defined as the root mean square of the 3D displacement errors. Additionally, differences at the first and last sample of each gait cycle, differences in minimum and maximum values, and the ROM between DFOD and OMCS for each gait cycle were computed and correlations between extrema and ROM were quantified with Pearson correlation coefficients. Correlations are interpreted as very strong for r = (0.90, 1.00), strong for r = (0.70, 0.89), moderate for r = (0.40, 0.69), weak for r = (0.20, 0.39), and very weak for r = (0.00, 0.19) [24]. The mediolateral and vertical axis of DFOD Equation (4a,b) were based on five gait cycles unless stated otherwise.

Algorithm Characteristics
DFOD assumes that the sensor on the lower leg moves quasi-cyclically and in a quasi-2D plane. To quantify how valid these assumptions are for the lower leg motion during treadmill running, respectively the mean cycle time and standard deviation and the explained variance of the first principal component of the angular velocity in Ψ s over one minute of running were computed. The effect of using data of different numbers of gait cycles to determine these axes of Ψ f and its error with respect to an OMCS was tested. Data of 1 up to 15 gait cycles were used to define the mediolateral and vertical axes of Ψ f , resulting in a total of 225 combinations which were tested. The outcome measure of this analysis was the 1D orientation and displacement estimate.
Full trust in the TRIAD algorithm [19] was given to the mediolateral functional axis Equation (4a) since this axis is not influenced by the violation of the approximately constant cycle average velocity assumption. The number of points used to estimate the rotation axis of Equation (9a,b) was based on a trial and error procedure to obtain a small variation in the obtained axes while using as few intervals as possible. Note that the results of this trial and error process were only used to validate DFOD and were not part of DFOD.
To investigate the effect of sampling frequency on the performance of DFOD, IMU data were resampled from 240 Hz to 120 Hz and 60 Hz before DFOD was used to estimate orientation and displacement. For this analysis, the vertical and mediolateral axis of DFOD were both based on five gait cycles and the 1D orientation and displacement estimates were used as outcome measures.

Results
An average of 79 gait cycles (range: 66-94) per subject were analyzed. When not stated otherwise, the mediolateral and vertical axes of DFOD Equation (4a,b) were based on data of five gait cycles.

Estimation of Orientation
Estimated lower leg sensor orientations without drift reduction, with drift reduction according to DFOD and from an OMCS are shown in Figure 4. Estimated lower leg sensor orientations of DFOD were compared to an OMCS in treadmill running. Mean RMSE for orientations in the sagittal plane were 3.1 ± 0.4 • while they were larger in the frontal (5.3 ± 1.1 • ) and transversal plane (5.0 ± 2.1 • ). The mean 1D rotation error (i.e., angle over which R f s needs to be rotated to coincide with R f cl ) was 7.5 ± 1.7 • . The 3D mean difference at the start and end of the gait cycle, absolute difference, and maximum and minimum difference in orientation together with the difference in ROM of DFOD and OMCS are shown in Table 1. Correlations between the 3D maximal angle, minimal angle, and ROM from DFOD and OMCS ranged from strong (0.768) to very strong (0.99). Mean 3D orientations of DFOD and OMCS for a representative subject are shown in Figure 5. Table 1. Mean orientation differences between DFOD and OMCS for all subjects combined. RMSE refers to the root mean square difference in 3D orientation. "Difference start cycle" and "Difference end cycle" refer to the difference between DFOD and OCMS (OMCS-DFOD) at the first and last sample of the gait cycle. "∆ Maximal angle" and "∆ Minimal angle" refer to the differences in the estimated maximal and minimal orientation during each gait cycle between DFOD and OMCS (OMCS-DFOD). "∆ ROM" refers to the mean differences in the estimated range of motion during each gait cycle for DFOD and OMCS. Pearson correlation coefficients (r) are provided between brackets.

Estimation of Displacement
Estimated lower leg sensor displacements of DFOD were compared to an OMCS in treadmill running. Mean RMSE for displacements in the forward direction were 1.6 ± 0.2 cm and similar for the mediolateral (1.7 ± 0.6 cm) and vertical direction (1.6 ± 0.2 cm). The mean 1D displacement error (i.e., length of the vector between the estimated sensor position of DFOD and OMCS) was 2.7 ± 0.4 cm. The 3D mean difference at the start and end of the gait cycle, absolute difference, maximum difference, and minimum difference in displacement, together with the difference in ROM of DFOD and OMCS, are shown in Table 2. Correlations between the 3D maximal displacement, minimal displacement, and ROM were moderate (r = 0.50) to strong (r = 0.82). Mean 3D displacements of DFOD and OMCS for a representative subject are shown in Figure 6. Table 2. Mean displacement differences between DFOD and OMCS for all subjects combined. RMSE refers to the root mean square difference in 3D sensor and cluster displacement. "Difference start cycle" and "Difference end cycle" refer to the difference between DFOD and OCMS (OMCS-DFOD) at the first and last sample of the gait cycle. "∆ Maximal displacement" and "∆ Minimal displacement" refer to the differences in the estimated maximal and minimal displacement during each gait cycle between DFOD and OMCS (OMCS-DFOD). "∆ ROM" refers to the mean differences in the estimated range of motion during each gait cycle for DFOD and OMCS. Pearson correlation coefficients (r) are provided between brackets.  Positive displacements in the X, Y, and Z-axis represent movement into the forward, sideward (left), and upward direction, respectively. The origin moves forward with the cycle average velocity.

Algorithm Characteristics
Two metrics were computed to show how valid the assumptions of a quasi-cyclical and quasi-2D movements were for treadmill running. The average cycle time was 0.68 ± 0.03 s/stride and the standard deviation ranged from 0.8-1.9% of the average cycle time. The first principal component of the angular velocity explained on average 90.2 ± 5.7% (range: 84.6-95.8%) of the variance.
The mediolateral and vertical axes of Ψ f Equation (4a,b) are based on data of five complete gait cycles. The effect of using data of more or less gait cycles to define these axes on the mean 1D orientation error is investigated and shown in Figure 7. The lowest mean 1D orientation error was found when the vertical axis was based on data of 11 gait cycles and the mediolateral axis on data of 8 gait cycles (mean error: 7.5 • ). The highest mean orientation error was found when the vertical and mediolateral axes were both based on data of 1 gait cycle (mean error: 7.7 • ). The effect of the number of gait cycles on the mean 1D displacement error is shown in Figure 8. The lowest mean displacement error was found when the vertical axis was based on data of 10 gait cycles and the mediolateral axis on data of 15 gait cycles (mean error: 2.6 cm). The highest mean displacement error was found when the vertical and mediolateral axes were both based on data of 1 gait cycle (mean error: 4.5 cm). To investigate the effect of sampling frequencies on DFOD, inertial data were resampled from 240 Hz to 120 Hz and 60 Hz before applying DFOD. Compared to a sampling frequency of 240 Hz, the 1D orientation error increased with 0.3 • for 120 Hz and 2.2 • for 60 Hz. The 1D displacement error increased with 1.2 cm for 120 Hz and 12.7 cm for 60 Hz.

Discussion
A new method, called Drift-Free Orientation and Displacement estimation (DFOD), is proposed to estimate drift-free 3D sensor orientation and displacement based on a single IMU. DFOD uses the quasi-cyclical behavior of human movements and assumes a quasi-2D movement with an approximately constant cycle average velocity. The performance of DFOD for a sensor on the lower leg was validated with an optical motion capture system (OMCS) in treadmill running. Errors in estimated sensor orientation and displacement between DFOD and OMCS were comparable to errors of other orientation and displacement algorithms. However, DFOD is independent of a constant-or zero-velocity point, a biomechanical model, a magnetometer, Kalman filtering, or a calibration procedure. Hence, DFOD is a promising method for quasi-cyclical motion analysis with a single IMU and has many advantages over current methods.

Estimation of Orientation
Estimated lower leg sensor orientations of DFOD were compared to an OMCS in treadmill running. DFOD performs best for orientation estimation in the sagittal plane, possibly because the largest ROM occurs around the axis perpendicular to this plane Equation (4a) in running.
To reduce drift in orientation estimation, a drift reducing rotation which was constant within each cycle, but varied over cycles, was applied (rotation 3). Orientation drift is relatively slow compared to the duration of a gait cycle (i.e., two min before Ψ d f drifts 90 • , or ±0.5 • /stride, around Ψ f y ). Hence, a constant drift reducing rotation for each gait cycle seemed sufficient, although this did result in small discontinuities between gait cycles. In future work, a continuous drift reducing rotation could improve the performance of DFOD.
Since we are not aware of studies that estimated lower leg orientations during running, the results of DFOD can only be compared with studies estimating foot and thigh orientations during running and walking. Foot orientations during running have mostly been based on constant-or zero-velocity updates with an additional drift reducing component (e.g., based on joint center accelerations, filtering, or an orientation reset). At speeds similar to our study, sagittal plane foot orientations could be estimated with errors varying between 2.0 • and 20.8 • [11,12,25]. Frontal plane foot orientation errors differed from 2.6 • to 4.4 • [12,25]. Upper leg orientations during walking have been estimated with an RMSE of 1.9 ± 0.5 • , although the zero acceleration and angular velocity update used in that study does not apply to continuous quasi-cyclical movements like running [26]. Orientation errors in our study are similar or slightly larger than found in literature for other body segments, although these studies used drift reducing methods unsuitable for a sensor on the lower leg in running (i.e., based on a constant-or zero-velocity point).
Tibial orientations in the sagittal and transversal plane are commonly studied with regard to running injuries [27][28][29]. The sagittal plane orientation of the tibia at initial contact has been shown to be 4.9 • larger in injured than in uninjured runners and the tibia ROM in the transverse plane is around 15 • in running [30]. With a mean difference of -0.3 ± 3.7 • at the start of the gait cycle (just before initial contact) and -5.6 ± 2.1 • in the transversal plane ROM, DFOD is capable to detect meaningful changes in tibia orientations during running.

Estimation of Displacement
Estimated lower leg sensor displacements of DFOD were compared to an OMCS in treadmill running. OMCS cluster marker placement can explain some of the errors in the forward and vertical directions. The OMCS cluster marker set is placed below the IMU (see Figure 1). Lower placement of the cluster marker set results in a larger ROM for OMCS compared to DFOD in the forward and vertical direction. Hence, actual displacement errors in the forward and vertical direction are expected to be smaller than those reported in this study.
Since we are not aware of studies that estimated lower leg displacements during running, the results of DFOD can only be compared with studies estimating foot displacements and stride length based on IMU data during running. In literature, estimates of sagittal plane foot displacement during running at a speed similar to the speed in this study had an absolute 1D positional error of 5 ± 2 cm at maximal foot height and initial contact [11]. The absolute 1D positional error in our study was 2.7 ± 0.4 cm. Previously, stride length based on an IMU in a shoe could be estimated with a mean absolute error of 7.6 cm [31]. DFOD has a mean absolute displacement error of 1.4 ± 0.2 cm in the forward direction. Hence, displacement errors of DFOD for the tibia sensor are smaller than those reported by literature for the foot segment in running.
DFOD estimates the displacement of a sensor on the lower leg. However, the displacement of each point on the tibia can be estimated based on the orientation of the sensor and the distance from the sensor to the point of interest. When the distance from the sensor to the ankle joint is known, the forward (step length) and upward (step height) displacement of the ankle can be estimated. Running velocity can then be obtained with the step length and cycle time. Hence, DFOD provides insight into the 3D trajectory of the lower leg during running and can be used to estimate step length, step height, and running velocity based on a single IMU on the lower leg.

Algorithm Characteristics
The assumptions that treadmill running is a quasi-cyclical and quasi-2D movement seem to hold based on the standard deviation of the cycle times (0.8-1.9% of the cycle time) and the explained variance of the first principal component for the angular velocity in Ψ s (84.6-95.8%). The explained variance shows that DFOD is capable of accurately estimating orientation and displacement even when 15% of the angular velocity in Ψ s occurs outside the 2D plane of a movement.
The effect of computing the functional mediolateral Equation (4a) and vertical Equation (4b) axes based on different numbers of gait cycles was found to be very small. The 1D orientation and displacement errors differed only 0.2 • and 1.9 cm between the best-and worst-performing combination of the number of included gait cycles. Hence, during indoor treadmill running at a constant velocity, the number of gait cycles for the vertical and mediolateral axes has a limited influence on the results of DFOD.
However, the goal is to apply DFOD in less controlled environments such as outdoor running. Outdoor running is likely to result in a less cyclical running pattern [32]. It is hypothesized that for outdoor running, a smaller number of gait cycles to compute the functional mediolateral Equation (4a) and vertical Equation (4b) axes is favored over a larger number since assumptions are less likely to be violated over shorter periods. Five gait cycles to define the vertical and mediolateral axes is expected to be a reasonable tradeoff between including more data to compensate for the increased variability in outdoor running while still being able to adapt to sudden changes in the gait pattern and reduce violations of assumptions. Hence, five gait cycles for both the mediolateral and vertical axes Equation (4a,b) were used in this study as the default setting for DFOD.
To investigate the effect of sampling frequencies on DFOD, inertial data were resampled from 240 Hz to 120 Hz and 60 Hz before applying DFOD. Orientation and displacement errors drastically increased when IMU data resampled to 60 Hz were used as input for DFOD. These results indicate that DFOD provides satisfactory results for a sampling frequency of 240 Hz and 120 Hz, but not for 60 Hz.

Limitations
Multiple assumptions were made to create DFOD, which can be violated by running outdoors. When runners run outside, they have a less constrained gait pattern than on a treadmill [32], and can freely change their running velocity and run up or downhill, thereby violating some assumptions of DFOD. Violation of the assumption of an approximately constant cycle average velocity does not influence the mediolateral axis of Ψ f Equation (4a) since this axis is based on the first principal component of the angular velocity of the lower leg sensor. Additionally, this axis is not influenced by taking a turn or running in circles, since it moves with the body. However, the vertical axis of Ψ f Equation (4b) is influenced by a violation of the approximately constant cycle average velocity assumption. This axis is equal to the direction of the total acceleration (i.e., including gravity) over a complete number of gait cycles when the cycle average velocity is constant. When a runner accelerates or decelerates, the free acceleration will not have a zero-mean over a complete number of gait cycles and will result in an offset in the estimated vertical axis proportional to the magnitude of the acceleration or deceleration. Since five gait cycles are included to estimate both functional axes, DFOD minimizes the effect of violated assumptions and is expected to recover from a short violation of assumptions within five gait cycles.
Similarly, the assumption of a quasi-cyclical 2D movement might be violated more often in running outdoors since impact accelerations are higher when running overground compared to a treadmill [33], due to uneven terrain, stumbling, or taking a turn. DFOD will recover from short violations of the quasi-cyclical 2D movement assumption within five gait cycles. Running-induced fatigue has been shown to increase variability in the gait pattern [34]. This increased variability and less cyclicity might cause the assumptions of DFOD to be less valid in fatigued running, resulting in larger orientation and displacement errors. Since DFOD has an origin that moves with the body at the cycle average velocity, a change in elevation caused by running on a sloped surface will cause the origin of DFOD to move up or down with the body. An elevation change will be visible over time; however, the average displacement will still be zero.
This study aimed to propose and validate a new algorithm that makes use of the quasi-cyclical nature of many movements. The algorithm was tested on treadmill running data of four runners and provided satisfactory results for all runners. Hence, to test the idea of using the quasi-cyclical nature of many human movements to estimate orientation and displacement, a limited number of subjects is appropriate. However, before DFOD can be used to study running kinematics it should be validated in more runners and different conditions. This study estimated sensor orientation and displacement during running while segment orientations might provide more insight for motion analysis. For a sensor to segment calibration, two axes that relate to both CSs are required. One of these axes is already defined in DFOD Equation (2a). The other axis could be based on the direction of the gravitational acceleration during neutral standing, in which the tibia is assumed to be vertical. However, this sensor to segment calibration does require an additional calibration procedure.

Future Research
In future work, DFOD should be validated in a less controlled setting, such as outdoor running, in multiple body segments, and different quasi-2D movements like cycling and skating. The influence of short violations of the assumptions of DFOD, increased variability in the gait pattern (i.e., caused by fatigue), less cyclical movements, and different speeds on estimated orientations and displacement should be assessed in (outdoor) running. Additionally, the effect of continuous drift reduction instead of a drift reduction during each gait cycle Equation (4e) could be evaluated to improve the performance of DFOD. As long as two functional axes can be defined, DFOD should be able to estimate sensor orientation and displacement. Hence, the generalized idea of DFOD could also be applied to quasi-cyclical 3D movements like swimming. For 3D movements, the validity of the functional mediolateral axis Equation (4a) based on the first principal component of the angular velocity should be assessed. This component is expected to be less pronounced in 3D versus 2D movements. Finally, a sensor to segment calibration procedure could be added to enable DFOD to calculate segment orientations instead of sensor orientations.

Conclusions
The Drift-Free Orientation and Displacement estimation method (DFOD) is proposed and validated. DFOD estimates drift-free 3D sensor orientation and displacement with a single IMU in quasi-cyclical quasi-2D plane movements with an approximately constant cycle average velocity. DFOD does not require a calibration procedure, biomechanical model, constant-or zero-velocity point, Kalman filtering, or magnetometer. Small errors in lower leg sensor orientation and displacement were found when DFOD was validated against an optical reference system in treadmill running. Hence, DFOD is a promising method for quasi-cyclical motion analysis, especially when using a minimal sensor setup.