Method for Estimating Three-Dimensional Knee Rotations Using Two Inertial Measurement Units: Validation with a Coordinate Measurement Machine

Three-dimensional rotations across the human knee serve as important markers of knee health and performance in multiple contexts including human mobility, worker safety and health, athletic performance, and warfighter performance. While knee rotations can be estimated using optical motion capture, that method is largely limited to the laboratory and small capture volumes. These limitations may be overcome by deploying wearable inertial measurement units (IMUs). The objective of this study is to present a new IMU-based method for estimating 3D knee rotations and to benchmark the accuracy of the results using an instrumented mechanical linkage. The method employs data from shank- and thigh-mounted IMUs and a vector constraint for the medial-lateral axis of the knee during periods when the knee joint functions predominantly as a hinge. The method is carefully validated using data from high precision optical encoders in a mechanism that replicates 3D knee rotations spanning (1) pure flexion/extension, (2) pure internal/external rotation, (3) pure abduction/adduction, and (4) combinations of all three rotations. Regardless of the movement type, the IMU-derived estimates of 3D knee rotations replicate the truth data with high confidence (RMS error < 4° and correlation coefficient r≥0.94).


Introduction
The human knee is susceptible to injury from multiple mechanisms including, for example, hyperextension (including varus and valgus components), over-use, and direct impact (see, for example, [1,2]). Knee injuries are also known to induce osteoarthritis and may severely limit activities of daily living [3,4]. Accordingly, the three-dimensional rotations across the knee (flexion/extension, internal/external rotation, and abduction/adduction) serve as important markers of knee health and performance and in multiple contexts such as human mobility, worker safety and health, athletic performance, and warfighter performance. For example, clinicians observe knee rotations to diagnose knee injuries and to assess the need for surgical interventions (such as knee arthroplasty), recovery, and physical therapy [5,6]. Many biomechanical analyses incorporate measurements of knee rotations including studies of the long-term effects from knee injuries [7,8], magnetic north from two IMUs may differ due to discrepancies in magnetometer data [28,29] despite these corrections, ultimately limiting this advantage. Reducing these discrepancies may follow from improving IMU hardware, updating filter parameters, or including additional (and complementary) sensors for fusion [30].
This study serves as an extension to those above by presenting, and carefully validating, a new method to estimate 3D rotations across the knee using a pair of shank-and thigh-mounted IMUs. Although this study employs IMUs with magnetometers, the method herein applies equally well to IMUs without magnetometers. Like [20][21][22], we use the medial/lateral axis of the knee as a kinematic constraint. However, in addition to estimating knee flexion/extension, our method also successfully estimates internal/external rotation and abduction/adduction. In so doing, the distinct world frames of the two IMUs are aligned by exploiting a vector form of the constraint equation in [21] that the medial/lateral axis have identical orientation in the two world frames. This constraint yields the (dynamic) correction needed to align the two world frames during, and also between, time intervals when the knee is functioning predominantly as a hinge joint, but without any requirement that the knee return to its initial orientation as in [24]. We open this paper with an overview of this new method together with a summary of companion experiments conducted on a coordinate measuring machine that yields high-precision truth data for validation.

Materials and Methods
As stated above, the objective of this paper is to present a new method for estimating 3D rotations across the knee that employs a kinematic constraint imposed by the knee's structure to remove drift and to carefully benchmark the accuracy of these estimates. To this end, we employ a coordinate measurement machine (CMM) (MicroScribe G2X, Solution Technologies, Oella, MD, USA) that functions as a knee analog and that provides truth data for validation. The CMM shown in Figure 1a embeds high precision optical encoders (0.0003 • resolution [31]) that measure rotations about three axes representing knee flexion/extension (FE), internal/external rotation (IE), and abduction/adduction (AA). Two IMUs (Opal sensors, APDM, Portland, OR, USA; Sensor characteristics and orientation estimate information available at http://www.apdm.com/wearable-sensors/), rigidly mounted to the illustrated two links, replicate the functions of thigh-and shank-mounted IMUs (T and S, respectively). The resulting apparatus enables direct comparison of IMU-estimated knee flexion/extension, internal/external rotation, and abduction/adduction to measured values from the three high precision optical encoders and over a wide range of simulated 3D knee movements.
Sensors 2017, 17,1970 3 of 15 may follow from improving IMU hardware, updating filter parameters, or including additional (and complementary) sensors for fusion [30]. This study serves as an extension to those above by presenting, and carefully validating, a new method to estimate 3D rotations across the knee using a pair of shank-and thigh-mounted IMUs. Although this study employs IMUs with magnetometers, the method herein applies equally well to IMUs without magnetometers. Like [20][21][22], we use the medial/lateral axis of the knee as a kinematic constraint. However, in addition to estimating knee flexion/extension, our method also successfully estimates internal/external rotation and abduction/adduction. In so doing, the distinct world frames of the two IMUs are aligned by exploiting a vector form of the constraint equation in [21] that the medial/lateral axis have identical orientation in the two world frames. This constraint yields the (dynamic) correction needed to align the two world frames during, and also between, time intervals when the knee is functioning predominantly as a hinge joint, but without any requirement that the knee return to its initial orientation as in [24]. We open this paper with an overview of this new method together with a summary of companion experiments conducted on a coordinate measuring machine that yields high-precision truth data for validation.

Materials and Methods
As stated above, the objective of this paper is to present a new method for estimating 3D rotations across the knee that employs a kinematic constraint imposed by the knee's structure to remove drift and to carefully benchmark the accuracy of these estimates. To this end, we employ a coordinate measurement machine (CMM) (MicroScribe G2X, Solution Technologies, Oella, MD, USA) that functions as a knee analog and that provides truth data for validation. The CMM shown in Figure 1a embeds high precision optical encoders (0.0003° resolution [31]) that measure rotations about three axes representing knee flexion/extension (FE), internal/external rotation (IE), and abduction/adduction (AA). Two IMUs (Opal sensors, APDM, Portland, OR, USA; Sensor characteristics and orientation estimate information available at http://www.apdm.com/wearablesensors/), rigidly mounted to the illustrated two links, replicate the functions of thigh-and shankmounted IMUs (T and S, respectively). The resulting apparatus enables direct comparison of IMUestimated knee flexion/extension, internal/external rotation, and abduction/adduction to measured values from the three high precision optical encoders and over a wide range of simulated 3D knee movements. Definitions of three frames of reference for a human knee associated with a shank-mounted IMU (blue) including the shank IMU frame, , the shank anatomical frame, , and the shank IMU's world frame, . Analogous frames of reference are illustrated for the thigh-mounted IMU (green). (b) Definitions of three frames of reference for a human knee associated with a shank-mounted IMU (blue) including the shank IMU frame, F S , the shank anatomical frame, F AS , and the shank IMU's world frame, F WS . Analogous frames of reference are illustrated for the thigh-mounted IMU (green).

Experimental Procedure
Data from the two IMUs are first time-synchronized to the encoder data from the CMM. The assembly is rotated by hand about the CMM's base (white axis in Figure 1a) with the three knee axes (FE, IE, AA) locked. The angle measured by the optical encoder about the base (dashed white) axis is differentiated with respect to time yielding an angular velocity signal to compare with those measured by the thigh (green) and shank (blue) IMUs. The data from the two IMUs are already time-synchronized, and their synchronization with the data from the CMM follows from measuring (and subsequently subtracting) the time delay between their respective angular rates.
Next, we complete two functional alignment movements that pre-identify the anatomical axes of the shank and thigh for the knee analog. First, the superior/inferior axes of the shank and thigh are estimated by holding each segment still (for approximately 10 s) while vertical. The measured acceleration for each segment defines the direction of gravity, which is also aligned with the superior/inferior axis of each segment. Next, the medial-lateral axis is established using essentially the same procedure outlined in [21]. In particular, the CMM is exercised purely about the flexion/extension axis with the two remaining knee axes locked. In so doing, the knee acts as a pure hinge, and one can readily compute the medial/lateral (hinge) axis with respect to the sense axes of each IMU. The resulting medial/lateral axes, so measured by the thigh-and shank-mounted IMUs, play a key role in the estimation process described below.
Finally, we consider four characteristically distinct knee movements for generating the truth data for validation. These movements, each repeated for N = 50 trials, include (1) pure flexion/extension, (2) pure internal/external rotation, (3) pure abduction/adduction, and (4) combinations of all three rotations. Each type of movement is made by hand, and one trial lasts approximately 10 s (with the appropriate CMM axes either free or locked).

Defining Segment Frames of Reference
Figure 1b illustrates three distinct frames of reference associated with each of the thigh-and shank-mounted IMUs. In particular, we call attention to the shank IMU frame, F S (defined by the IMU sense axes), the shank anatomical frame, F AS , and the shank IMU world frame, F WS . The three analogous frames of reference (F T , F AT , F WT ) associated with the thigh-mounted IMU are also illustrated. The ultimate goal is to estimate the 3D rotations (FE, IE, AA) across the knee and doing so requires estimating the orientation of the shank anatomical frame, F AS , relative to the thigh anatomical frame, F AT . That critical step, however, requires introducing a common world frame of reference for the two IMUs as described below. Prior to that, we first establish the separate world frames and anatomical frames for both segments described below.
The quaternion output (which is provided by proprietary software from APDM; for examples on how to compute quaternion output from IMU data, refer to [18,20,24,25]) from an IMU say S is used to construct a direction cosine matrix (DCM) between an IMU frame, F S , and a world frame, F WS . The world frame F WS is defined by three mutually orthogonal axes (X WS ,Ŷ WS ,Ẑ WS ) with theẐ WS axis chosen to align with gravity (using accelerometer data from S), theX WS -axis chosen to align with magnetic north (using magnetometer data from S), and theŶ WS axis computed from Y WS =Ẑ WS ×X WS and thus chosen to point west. Let R WS/S represent the resulting DCM from F S to F WS , a result that necessarily utilizes the magnetic north estimate from S. An analogous procedure holds for the thigh-mounted IMU leading to the construction of R WT/T representing the DCM from F T to F WT , a result that necessarily utilizes the magnetic north estimate from T. Note also that the location and orientation of either IMU on its respective segment are arbitrary and the orientation of each IMU relative to its respective anatomical axes is established using the aforementioned functional alignment movements as detailed next.
The two functional alignment movements establish the shank anatomical axes (X AS ,Ŷ AS ,Ẑ AS ) that define F AS and the thigh anatomical axes (X AT ,Ŷ AT ,Ẑ AT ) that define F AT . The procedure for both body segments is identical, and so we detail only that for the shank. First, the average acceleration measured during the still period yields a first estimate of the shank-fixedẐ AS axis (superior-inferior axis) of the shank anatomical frame that is approximately aligned with gravity. (Note that during a trial with a human subject standing still, the direction of gravity remains approximately aligned with the superior-inferior axis.) This still period is followed by rotations purely about the CMM flexion/extension axis (while locking the other two rotational degrees of freedom). Following Seel et al. [21], the resulting angular velocity ω S and ω T measured by S and T, respectively, is used to define a unit vectorn S aligned with the medial-lateral (hinge) axis measured in F S . The medially pointing direction of this unit vector defines the anatomicalX AS axis. The anterior-posterior axis,Ŷ AS , follows immediately fromŶ Finally, we adjust (if needed) the superior-inferior axis so that the medial-lateral axis remains orthogonal to the other two anatomical axes per, The resulting orthonormal triad (X AS ,Ŷ AS ,Ẑ AS ), which are measured with respect to the shank IMU frame F S , define the shank anatomical frame with the hinge axisn S =X AS The (constant) DCM from the shank IMU frame, F S , to the shank anatomical frame, F AS , follows from, where each row contains the components of the anatomical axes measured with respect to F S . Per ISB convention [32], the medial-lateral axis corresponds to the FE axis, the anterior-posterior axis corresponds to the AA axis, and the superior-inferior axis corresponds to the IE axis. An analogous procedure establishes the knee hinge axisn T =X AT and the thigh anatomical frame (X AT ,Ŷ AT ,Ẑ AT ), which are measured with respect to the thigh IMU frame F T . The (constant) DCM from the thigh IMU frame, F T , to the thigh anatomical frame, F AT , follows from,

Estimating 3D Knee Rotations Following Construction of a Common World Frame
The dynamic 3D knee rotations are ultimately estimated from the (time-varying) DCM R(t) AT/AS from the shank anatomical frame, F AS , to the thigh anatomical frame, F AT . One may believe that this DCM follows from the component rotations defined above per, However, this result is correct only in the rare instances when the two IMU world frames, F WS and F WT , are aligned. Frequently, the magnetometers in the two IMUs, S and T, provide distinct estimates of magnetic north (especially indoors where ferromagnetic interferences are often prevalent), and thus their respective world frames F WS and F WT will be misaligned in general. This is also true for methodologies that do not employ magnetometers, where changes in orientation are estimated from world frames constructed by other means/assumptions. Furthermore, the estimated world frames often vary with time due to sensor drift (bias) errors. In short, the two IMUs are independent sensors yielding independent and time-varying (drifting) world frames as also illustrated in the example results that follow. Therefore, the key challenge lies in constructing a common world frame of reference for the two IMUs, or equivalently, estimating the "correction" DCM C(t) WT/WS from F WS to F WT . There is no clear way of determining the angular differences in the estimates of magnetic north from the two IMUs. However, one can exploit the constraint that the hinge (medial-lateral) axes,n T andn S , should be identical in a common world frame during the time intervals when the knee acts predominantly as a hinge joint. We treat this "hinge constraint" as a full vector equation and therefore generalize the scalar treatment of this constraint employed in [21]. Whenever this hinge criterion is satisfied, the correction DCM C(t) WT/WS can be constructed fromn T andn S by computing the axis of rotation,k, and the associated angle of rotation, θ, needed to alignn T andn S in their respective world frames. To this end, the cross product, defines the axis of rotationk(t), and the dot product, defines the required angle of rotation θ(t). The requisite correction DCM follows from Rodrigues' rotation formula [33] per, where k(t) is the skew symmetric form ofk(t). The correction DCM corrects the small misalignment between the two world frames and between successive times when the hinge criterion (described below) is satisfied. The resulting corrected form of Equation (5) becomes, where R(t) AT/AS is again the needed time-varying direction cosine describing the orientation of the anatomical shank frame relative to the anatomical thigh frame. Note that the ISB recommends [31] first calculating FE, then IE, and finally AA. However, due the mechanical design of the knee analog, the reverse order is required. For the knee analog, the joints are in series order such that the rotation sequence needed to rotate from the shank link to the thigh link require rotations about the AA axis first, the IE axis second, and the FE axis third. Decomposing R AT/AS with that order in mind (analogous to the procedure in [24]) yields the 3D rotation angles across the knee analog. The above strategy for aligning the two world frames holds for the time intervals when the knee analog in this experiment is predominantly functioning as a hinge joint. To this end, one must develop criteria for determining: (1) when the knee is predominantly functioning as a hinge joint and; (2) a strategy for aligning the two world frames in between those time intervals (i.e., when the knee is no longer predominantly functioning as a hinge joint).
The time intervals when the knee is predominantly functioning as a hinge joint are identified under two conditions: Case (1a) when the knee joint is approximately stationary; and Case (1b) when the knee joint is rotating.
(1a). Functioning as a Hinge-Stationary Case: The knee analog may be considered a hinge joint in the limiting case when it is essentially locked (i.e., straight) and stationary. These time intervals are identified when the segments are stationary and also nominally aligned with gravity (such as the still period in the functional alignment movement sequence). In such instances, the knee functions as a hinge with zero hinge rotation rate. These conditions are considered satisfied in our experiments when, in which → a (0) T and → a (0) S denote the (constant) acceleration measured during the functional alignment step with the segments vertical. (1b). Functioning as a Hinge-Rotating Case: When the conditions for the stationary case above are not met, the knee joint is considered rotating. During such instances, the knee may still function predominantly as a hinge joint (with non-zero hinge rotation rate) whenever the angular velocities of the shank ω S and the thigh ω T are predominantly aligned with the hinge axis defined byn T andn S on each segment respectively. For our experiment, the knee functions as a hinge with non-zero hinge rotation rate whenever The numerical thresholds for the criteria above are stringent, but appropriate for the knee analog we employ in our study. (These thresholds would need to be adjusted for human subject testing as discussed further in the Results and Discussion section.) There are frequent periods of times when neither criteria are met and the knee analog no longer functions purely as a hinge as described next.
(2) Not Functioning as a Hinge: When neither case above holds, the knee is no longer functioning as a hinge and appreciable internal/external rotation and/or abduction/adduction exists. During such time intervals, we assume that C(t) WT/WS varies slowly and continuously and that it can be estimated by linear interpolation between two consecutive "update" times when either of the two cases above hold.

Results and Discussion
This section contains both qualitative and quantitative comparisons and discussions of the angles estimated using IMU data and the truth data obtained from the CMM. A brief discussion of extending the method for measurements on a human knee concludes the section.
We begin by noting that the DCMs between the sensor frames and their respective world frames are defined by both horizontal (or yaw) and vertical (or elevation) angular components. Studying these angular components is important for understanding the challenge (and the solution) to constructing the correction DCM C WT/WS introduced above that corrects the differing and drifting IMU world frames. Figure 2 illustrates the orientation of the shank IMU frame F S relative to its world frame F WS in terms of both yaw ψ and vertical ϕ angular components. In particular, the unit vectorsX WS andŶ WS define the "horizontal" plane for F WS , which differ from that of F S by the vertical angle ϕ. Moreover, the frames also differ by the yaw angle ψ, which is defined as the angle betweenX WS and the projection ofX S onto the horizontal plane for F WS . One way to visualize the drift error is to examine how either angle varies with time over a trial during which the knee is periodically returned to its initial orientation.
For example, Figure 3a illustrates how the yaw angles computed for both IMUs (relative to their respective world frames) vary with time for the simplest testing session when the knee undergoes pure flexion/extension. The portion of the testing session shown encompasses the two functional alignment movements (first shaded interval corresponds to the still period and the first unshaded interval corresponds to FE rotations), followed by another still period (second shaded interval), then a longer (second unshaded) interval with four trials of five repetitive knee flexion/extension movements (with 46 more trials thereafter that are not shown). The knee is approximately returned to its original position at the end of each knee flexion/extension trial. Note that this approximately 3-min sample is taken from a testing session lasting approximately 15 min, and substantial drift arises over this long period of time. Figure 3b illustrates the Boolean values for the criteria for the two cases (Case 1a and Case 1b) when the knee acts as a hinge; a value of 1 corresponds to the criteria being satisfied, and a value of 0 corresponds to the criteria not being satisfied. As expected, the knee analog functions as a hinge for the stationary case (Case 1a) during the rest periods. Similarly, the knee analog functions as a hinge during the subsequent rotations that induce pure flexion/extension for this trial (Case 1b). (There are also short intervals when Case 1b is not satisfied as the angular velocity magnitudes fall below the selected threshold values.) angles estimated using IMU data and the truth data obtained from the CMM. A brief discussion of extending the method for measurements on a human knee concludes the section.
We begin by noting that the DCMs between the sensor frames and their respective world frames are defined by both horizontal (or yaw) and vertical (or elevation) angular components. Studying these angular components is important for understanding the challenge (and the solution) to constructing the correction DCM ⁄ introduced above that corrects the differing and drifting IMU world frames. Figure 2 illustrates the orientation of the shank IMU frame relative to its world frame in terms of both yaw and vertical angular components. In particular, the unit vectors and define the "horizontal" plane for , which differ from that of by the vertical angle . Moreover, the frames also differ by the yaw angle , which is defined as the angle between and the projection of onto the horizontal plane for . One way to visualize the drift error is to examine how either angle varies with time over a trial during which the knee is periodically returned to its initial orientation. For example, Figure 3a illustrates how the yaw angles computed for both IMUs (relative to their respective world frames) vary with time for the simplest testing session when the knee undergoes pure flexion/extension. The portion of the testing session shown encompasses the two functional alignment movements (first shaded interval corresponds to the still period and the first unshaded interval corresponds to FE rotations), followed by another still period (second shaded interval), then a longer (second unshaded) interval with four trials of five repetitive knee flexion/extension movements (with 46 more trials thereafter that are not shown). The knee is approximately returned to its original position at the end of each knee flexion/extension trial. Note that this approximately 3min sample is taken from a testing session lasting approximately 15 min, and substantial drift arises  The orientation drift error manifests in the slowly changing (low frequency) components of and despite the fact that the knee returns to its nominal orientation between successive flexion/extension movements. In particular, the net change in and over the entire 15min period are 65° and 101°, respectively. This corresponds to drift rates for the thigh and shank sensors of −0.07°/s and −0.11°/s, respectively, which are consistent with previously measured drift rates for MEMS inertial sensors; see, for example, [34]. The two world frames are misaligned at the start of the trial and continue to drift apart throughout the trial and it should also be noted the drift is not exclusively about the vertical ( ) axis. The misalignment of the world frames, and its associated drift, is compounded by sources of ferromagnetic interference in the laboratory environment. Despite the observably large misalignment of the world frames, the method above correctly identifies the correction DCM needed to accurately estimate 3D rotations across the knee analog. Figure 4 illustrates the differences between the true (CMM) and estimated (IMU) flexion/extension, internal/external rotation and abduction/adduction angles for one trial of five The orientation drift error manifests in the slowly changing (low frequency) components of ψ thigh and ψ shank despite the fact that the knee returns to its nominal orientation between successive flexion/extension movements. In particular, the net change in ψ thigh and ψ shank over the entire 15-min period are 65 • and 101 • , respectively. This corresponds to drift rates for the thigh and shank sensors of −0.07 • /s and −0.11 • /s, respectively, which are consistent with previously measured drift rates for MEMS inertial sensors; see, for example, [34]. The two world frames are misaligned at the start of the trial and continue to drift apart throughout the trial and it should also be noted the drift is not exclusively about the vertical (ẑ) axis. The misalignment of the world frames, and its associated drift, is compounded by sources of ferromagnetic interference in the laboratory environment. Despite the observably large misalignment of the world frames, the method above correctly identifies the correction DCM needed to accurately estimate 3D rotations across the knee analog. Figure 4 illustrates the differences between the true (CMM) and estimated (IMU) flexion/ extension, internal/external rotation and abduction/adduction angles for one trial of five repetitive flexion/extension movements previously considered in Figure 3. These estimates closely track the truth values reported by the optical encoders; refer also to quantitative comparison that follows. Also shown are the differences between the true and (uncorrected) estimates that result from using Equation (5) instead of (9), or equivalently, assuming that C WT/WS = I. Clearly, ignoring this correction leads to large errors and for all three rotation angles. Note that throughout this trial, the knee analog functions largely as a hinge, meaning that either Case 1a or Case 1b is almost always satisfied (and the linear interpolation associated with Case 2 is not required)-refer to results illustrated in Figure 3b. However, this is not the case in the following trials that induce substantial internal/external rotation and/or abduction/adduction.   Figure 5 illustrates example results for two trials that consider (a) pure internal/external rotation and (b) pure abduction/adduction. As with the prior case, the differences between the true and estimated flexion/extension, internal/external rotation, and abduction/adduction angles remain small, meaning the estimates closely track the truth values reported by the optical encoders. For instance, observe the very slight "off-axis" errors along rotation axes that are otherwise physically constrained on the knee analog. These small off-axis rotations likely derive from very minor errors in the estimated orientations of one or both of the anatomical frames, and . Nonetheless, the agreement between the IMU estimates and the truth data remains excellent. By contrast, large errors again arise in general when one ignores the correction. The agreement with the correction for these limiting cases of 1D rotation remains undiminished when the knee analog is unconstrained and undergoes fully three-dimensional rotations as shown next. Example results from pure flexion/extension trial. The difference between CMM truth data and IMU-derived estimates with the correction (black) and without the correction (red) of flexion/extension, internal/external rotation, and abduction/adduction are plotted for a representative time period. Figure 5 illustrates example results for two trials that consider (a) pure internal/external rotation and (b) pure abduction/adduction. As with the prior case, the differences between the true and estimated flexion/extension, internal/external rotation, and abduction/adduction angles remain small, meaning the estimates closely track the truth values reported by the optical encoders. For instance, observe the very slight "off-axis" errors along rotation axes that are otherwise physically constrained on the knee analog. These small off-axis rotations likely derive from very minor errors in the estimated orientations of one or both of the anatomical frames, F AT and F AS . Nonetheless, the agreement between the IMU estimates and the truth data remains excellent. By contrast, large errors again arise in general when one ignores the correction. The agreement with the correction for these limiting cases of 1D rotation remains undiminished when the knee analog is unconstrained and undergoes fully three-dimensional rotations as shown next.
constrained on the knee analog. These small off-axis rotations likely derive from very minor errors in the estimated orientations of one or both of the anatomical frames, and . Nonetheless, the agreement between the IMU estimates and the truth data remains excellent. By contrast, large errors again arise in general when one ignores the correction. The agreement with the correction for these limiting cases of 1D rotation remains undiminished when the knee analog is unconstrained and undergoes fully three-dimensional rotations as shown next. We first present a sample time record of the yaw angles for each IMU in Figure 6a when the knee undergoes combined 3D rotations. Following Figure 3, Figure 6 includes a portion of a testing session that encompasses the two functional alignment movements (first shaded and unshaded intervals, respectively), followed by a nominally still period during which the constraints are removed (second shaded interval), then a longer interval with six repetitive 3D movements (with four more thereafter that are not shown). The knee analog is approximately returned to its original position at the end of each knee movement. Note that this is taken from a testing session lasting approximately 6 min, and substantial drift arises over this period of time. Compared to the results of Figure 3, the results of this longer time sample in Figure 6a exhibit even greater drift error as again manifested in the slowly changing (low frequency) components of and . The net changes in and over the entire 6-minute period is 46° and 28°, respectively, which correspond to drift rates of 0.14°/s and −0.08°/s. Despite the very apparent drift, the method produces excellent estimates of the 3D rotations. Figure 6b documents the Boolean values for the criteria for the two cases (Case 1a and Case 1b) when the knee acts as a hinge; a value of 1 corresponding to the criteria being satisfied and a value of 0 corresponding to the criteria not being satisfied. As before, the knee analog functions as a hinge for the stationary case (Case 1a) during rest periods and functions as a hinge during the flexion/extension functional alignment movement (Case 1b). Note that during the second shaded We first present a sample time record of the yaw angles for each IMU in Figure 6a when the knee undergoes combined 3D rotations. Following Figure 3, Figure 6 includes a portion of a testing session that encompasses the two functional alignment movements (first shaded and unshaded intervals, respectively), followed by a nominally still period during which the constraints are removed (second shaded interval), then a longer interval with six repetitive 3D movements (with four more thereafter that are not shown). The knee analog is approximately returned to its original position at the end of each knee movement. Note that this is taken from a testing session lasting approximately 6 min, and substantial drift arises over this period of time. Compared to the results of Figure 3, the results of this longer time sample in Figure 6a exhibit even greater drift error as again manifested in the slowly changing (low frequency) components of ψ thigh and ψ shank . The net changes in ψ thigh and ψ shank over the entire 6-minute period is 46 • and 28 • , respectively, which correspond to drift rates of 0.14 • /s and −0.08 • /s. Despite the very apparent drift, the method produces excellent estimates of the 3D rotations. Figure 6b documents the Boolean values for the criteria for the two cases (Case 1a and Case 1b) when the knee acts as a hinge; a value of 1 corresponding to the criteria being satisfied and a value of 0 corresponding to the criteria not being satisfied. As before, the knee analog functions as a hinge for the stationary case (Case 1a) during rest periods and functions as a hinge during the flexion/extension functional alignment movement (Case 1b). Note that during the second shaded area, the physical constraints are being removed from the CMM in preparation for the combination rotation trials to follow.
shaded interval), then a longer interval with six repetitive 3D movements (with four more thereafter that are not shown). The knee analog is approximately returned to its original position at the end of each knee movement. Note that this is taken from a testing session lasting approximately 6 min, and substantial drift arises over this period of time. Compared to the results of Figure 3, the results of this longer time sample in Figure 6a exhibit even greater drift error as again manifested in the slowly changing (low frequency) components of and . The net changes in and over the entire 6-minute period is 46° and 28°, respectively, which correspond to drift rates of 0.14°/s and −0.08°/s. Despite the very apparent drift, the method produces excellent estimates of the 3D rotations. Figure 6b documents the Boolean values for the criteria for the two cases (Case 1a and Case 1b) when the knee acts as a hinge; a value of 1 corresponding to the criteria being satisfied and a value of 0 corresponding to the criteria not being satisfied. As before, the knee analog functions as a hinge for the stationary case (Case 1a) during rest periods and functions as a hinge during the flexion/extension functional alignment movement (Case 1b). Note that during the second shaded area, the physical constraints are being removed from the CMM in preparation for the combination rotation trials to follow.  Figure 7 illustrates the estimated flexion/extension, internal/external rotation and abduction/adduction angles for a representative time period from Figure 6 during which all three angles were being simultaneously exercised. Consistent with the above results, the estimates of all three angles closely track the truth values reported by the optical encoders when the correction is employed. When the correction is not employed, the estimates can be rather poor and particularly so for abduction/adduction in this example.  Figure 7 illustrates the estimated flexion/extension, internal/external rotation and abduction/ adduction angles for a representative time period from Figure 6 during which all three angles were being simultaneously exercised. Consistent with the above results, the estimates of all three angles closely track the truth values reported by the optical encoders when the correction is employed. When the correction is not employed, the estimates can be rather poor and particularly so for abduction/adduction in this example. Figure 7 illustrates the estimated flexion/extension, internal/external rotation and abduction/adduction angles for a representative time period from Figure 6 during which all three angles were being simultaneously exercised. Consistent with the above results, the estimates of all three angles closely track the truth values reported by the optical encoders when the correction is employed. When the correction is not employed, the estimates can be rather poor and particularly so for abduction/adduction in this example. Example results from combined 3D rotation trial. The differences between the CMM truth data and the IMU-derived estimates with the correction (black) and without the correction (red) for flexion/extension, internal/external rotation, and abduction/adduction are plotted for a representative time period.

Quantitative Comparisons
The above results illustrate very close qualitative agreement between the IMU-derived estimates of the joint angles and those measured directly using the embedded optical encoders. Next, we provide quantitative comparisons for the entire data set. To start, consider Figure 8 which shows the Example results from combined 3D rotation trial. The differences between the CMM truth data and the IMU-derived estimates with the correction (black) and without the correction (red) for flexion/extension, internal/external rotation, and abduction/adduction are plotted for a representative time period.

Quantitative Comparisons
The above results illustrate very close qualitative agreement between the IMU-derived estimates of the joint angles and those measured directly using the embedded optical encoders. Next, we provide quantitative comparisons for the entire data set. To start, consider Figure 8 which shows the IMU-estimated flexion/extension angle versus that measured by the optical encoder for the duration of all five testing sessions (cumulatively about 30 min) of combined 3D rotations. The best fit line to this data yields a slope of 1.01, a y-intercept of 0.07 ( • ), and a correlation coefficient of r = 0.99. Thus, the estimates exhibit extremely high correlation with the truth data (and just slightly over predicting flexion/extension relative to the truth data). In addition, the root-mean square (RMS) error between the estimates and the truth data is 3.46 • or 2.96% relative to the 117 • range of motion. These results, and the analogous quantitative comparisons for all of the experiments, are summarized in Tables 1  and 2

below.
Sensors 2017, 17,1970 12 of 15 IMU-estimated flexion/extension angle versus that measured by the optical encoder for the duration of all five testing sessions (cumulatively about 30 min) of combined 3D rotations. The best fit line to this data yields a slope of 1.01, a y-intercept of 0.07 (°), and a correlation coefficient of r = 0.99. Thus, the estimates exhibit extremely high correlation with the truth data (and just slightly over predicting flexion/extension relative to the truth data). In addition, the root-mean square (RMS) error between the estimates and the truth data is 3.46° or 2.96% relative to the 117° range of motion. These results, and the analogous quantitative comparisons for all of the experiments, are summarized in Tables 1  and 2 below.  Table 1 reports the quantitative comparisons of IMU-estimated angles to those measured by the optical encoders for the three limiting cases of pure rotation about the FE, IE and AA axes. Reported are the range of motion (ROM) about each axis, the RMS error between the estimated and measured angles, and the correlation coefficient, slope and y-intercept of the associated linear fit. Analogous results are reported in Table 2 for the combined 3D rotation movements. Regardless of the trial (pure rotation about any one axis or combined rotation about all three axes), the IMU-derived angle  Table 1 reports the quantitative comparisons of IMU-estimated angles to those measured by the optical encoders for the three limiting cases of pure rotation about the FE, IE and AA axes. Reported are the range of motion (ROM) about each axis, the RMS error between the estimated and measured angles, and the correlation coefficient, slope and y-intercept of the associated linear fit. Analogous results are reported in Table 2 for the combined 3D rotation movements. Regardless of the trial (pure rotation about any one axis or combined rotation about all three axes), the IMU-derived angle estimates remain within 4 • of those measured by the optical encoders, have correlation coefficients exceeding 0.94 and slopes between 0.99 and 1.02. In short, the method yields estimates that replicate the truth data with high confidence. The ranges of motion for each angle are motivated by the behavior of the human knee for which flexion/extension typically has the greatest range of motion, followed by internal/external rotation, and abduction/adduction. However, the knee analog is exercised over far greater ranges of motion than observable on a healthy human knee for the purpose of this thorough validation study. Although the results presented above demonstrate excellent agreement with the truth data, there are possible remaining error sources in addition to those associated with the correction DCM C WT/WS . Referring to Equation (9), errors may arise from estimates of R T/WT and R WS/S that are provided by the chosen Kalman filtering method. As illustrated in Figures 3a and 6a, substantial drift arises about the vertical axis despite the filtering method (and the use of magnetometer data). Errors may also arise from estimates of R AT/T and R S/AS which, as already noted above, are likely responsible for the very slight off-axis rotations arising in Figure 5.

Extension to Measurements on the Human Knee
The above method yields highly accurate results when benchmarked on a coordinate measurement machine (CMM). It therefore it has considerable potential for accurately estimating the 3D rotations across the human knee. Doing so will also require further research to address additional considerations. First, unlike the CMM used herein, the knee joint has some laxity which may also influence the estimated 3D joint angles. Second, the stringent limits used in Equations (10)-(13), while fully appropriate for the CMM, must necessarily be relaxed for the human knee. Careful validation tests using human subjects will guide the selection of these new limits. Third, while the IMUs were rigidly fastened to the links of the CMM, the IMUs attached to the human shank and thigh may move slightly relative to the underlying skeleton due to soft tissue movements, a problem common to all motion capture methods. The accuracy of the estimated 3D rotations also depends on accurately establishing the orientation of the shank and thigh anatomical frames relative to their respective IMU frames as input to Equation (9). This motivates the need to study a variety of functional alignment movements to achieve this intermediate result.

Summary and Conclusions
This study contributes a new method to estimate the three-dimensional rotations across the knee using a pair of shank-and thigh-mounted IMUs. Central to this method is constructing a common world frame of reference for the two IMUs by exploiting a vector constraint equation for the medial-lateral axis of the knee. This constraint arises during the time periods when the knee behaves predominantly as a hinge joint, during which and between times the common world frame can be constructed despite sensor orientation drift errors. Estimates of knee flexion/extension, internal/external rotation, and abduction/adduction follow from an estimate of the DCM between the shank and thigh anatomical frames of reference. The method is carefully tested using a coordinate measuring machine that kinematically replicates the 3D rotations across the human knee joint. The embedded optical encoders yield high precision truth data over four qualitatively distinct knee movements spanning (1) pure flexion/extension, (2) pure internal/external rotation, (3) pure abduction/adduction, and (4) combinations of all three rotations. Over N = 50 trials, the RMS error and linear correlation coefficient between the IMU estimates and the optical encoder measurements is 3.90 • and r = 0.99 for pure flexion/extension, 1.83 • and r = 0.99 for pure internal/external rotation, and 0.12 • and r = 0.99 for pure abduction/adduction. During combined movements, the RMS errors and correlation coefficients remain largely unchanged (3.46 • and r = 0.99 for flexion/extension, 2.48 • and r = 0.99 for internal/external, and 1.69 • and r = 0.94 for abduction/adduction.) In summary, the IMU-derived estimates of 3D knee rotations replicate the truth data with high confidence. Our future work includes extending this methodology to trials on human subjects with direct comparison to the 3D knee rotations measured using optical motion capture.