Comparison of Three Hand Pose Reconstruction Algorithms Using Inertial and Magnetic Measurement Units

The correct estimation of human hand kinematics has received a lot of attention in many research fields of neuroscience and robotics. Not surprisingly, many works have addressed hand pose reconstruction (HPR) problem and several technological solutions have been proposed. Among them, Inertial and Magnetic Measurement Unit (IMMU) based systems offer some elegant characteristics (including cost-effectiveness) that make these especially suited for wearable and ambulatory HPR. However, what still lacks is an exhaustive characterization of IMMU-based orientation tracking algorithms performance for hand tracking purposes. In this work, we have developed an experimental protocol to compare the performance of three of the most widely adopted HPR computational techniques, i.e. extended Kalman filter (EKF), Gauss-Newton with Complementary filter (CF) and Madgwick filter (MF), on the same dataset acquired through an IMMU-based sensing glove. The quality of the algorithms has been benchmarked against the ground truth measurement provided by an optical motion tracking system. Results suggest that performance of the three algorithms is similar, though the MF algorithm appears to be slightly more accurate in reconstructing the individual joint angles during static trials and to be the fastest one to run.


I. INTRODUCTION
The hand is one of the most important tools that humans use to interact with the environment. It not only represents the primary end-effector for object manipulation and grasping, but also the cognitive organ of the sense of touch and a major means of non-verbal communication. Understanding and studying how hands are used in tasks of daily living is important for various reasons, including human motion analysis, rehabilitation, neuro-prosthetics, advancedhuman -machine interaction, tele-robotics and robotics [1,2]. Indeed, accurate hand pose reconstruction (HPR) can provide crucial information not only to assess manipulative task performance in humans, but also for design and control of humanoids, e.g. to devise suitable guidelines for the development of artificial hand systems, including prostheses [3]. However, ideal HPR is difficult to achieve due to the complex biomechanics of the hand and a large number of degrees of freedom (DoFs) involved. Over the years, many different systems for monitoring and quantifying hand motion have been proposed, such as camera-based marker systems and instrumented gloves (equipped, for instance, with piezo-resistive, fibro-optic or inertial-magnetic sensors). Interested readers can refer to [4] for more details. All these systems rely on different sensing technologies and apply different mathematical models to decode the pose of the hand and the task being performed. This constitutes one of the major obstacles to compare study results obtained by different studies and to provide reliable benchmarking. Although optical systems are widely considered the gold standard in motion tracking, non-vision based wearable systems are more suitable (and commonly employed) for ambulatory and daily life monitoring, providing an ecological tool to study human behavior [4]. Among wearable solution, Inertial and Magnetic Measurement Unit (IMMU) based systems offer some elegant characteristics, including but not limited to a more effective integration with soft-robotic hands hand bodies [5]. Despite the fact that IMMU-based orientation tracking algorithms have been compared extensively in literature, no comparative studies have been carried out in the context of HPR. Note that IMMU are sometimes also referred to as IMU (Inertial measurement unit) or MARG (Magnetic, Angular Rate, and Gravity unit). In this paper, we will consistently use the acronym IMMU.
The aim of this study is to make a step further towards standardization and benchmarking of the three mostly widely used state of the art IMMU-based hand pose reconstruction algorithms. The three HPR algorithms considered in this study are: Extended Kalman Filter (EKF), Gauss-Newton method combined with the Complementary Filter (CF) and Madgwick Filter (MF). In the comparison, we focused on algorithms' repeatability; static hand pose reconstruction and dynamic hand pose reconstruction. We assessed accuracy of the pose reconstruction and computational complexity. Quality of the algorithms was benchmarked against the ground truth measurement provided by an optical motion tracking system.

A. Kinematic Hand Model
In this article, a 20 DoF kinematic hand model is used as shown in Figure 1. In this model, each digit has a proximal, intermediate and distal phalanx. Each proximal interphalangeal (PIP) and distal interphalangeal (DIP) joint has a single flexion-extension axis (1 DoF), while the metacarpophalangeal (MCP) joint at the base of the digit has both flexion-extension and abduction-adduction axes (2 DoFs). The thumb is modeled with 4 DoFs corresponding to rotation and abduction/adduction at the carpometacarpal (CMC) joint (2 DoFs), flexion/extension at the MCP joint (1 DoF) and flexion/extension at the interphalangeal (IP) joint (1DoF). In this model, no pronation/supination of the fingers is considered and the position of the metacarpals is considered fixed. All the joints are assumed to be ideal hinge/universal joints whose axes are fixed relative to their associated links. This model is the same as proposed in [6].   Two optical markers were placed on each of the finger segments with three markers on each of the hand units. An additional marker was placed over the styloid ulnar process to create a stable, and always visible, hand base frame.

B. Hand Tracking Experimental Setup
Hand tracking data was acquired simultaneously with two independent systems. The IMMU-based motion tracking data was recorded with the PowerGlove system [7]. Optical motion tracking data was acquired with VICON system and used as a ground truth measurement. These two systems are described in more detail in the following subsections.

1) IMMU System
To capture the data used for the algorithm comparison, a PowerGlove system was used. The PowerGlove ( Figure 2) is an IMMU-based device that uses inertial (gyroscopes and accelerometers) and magnetic sensors placed on various hand and finger segments to accurately assess full 3D hand kinematics. An explanation of the anatomical acronyms used is given in TABLE I. The device acquires 100 samples per second from the accelerometers and magnetometers and 200 samples per second from the integrated gyroscopes. This device has been developed at the University of Twente and previously validated against an optical tracking system [7]. During this study eighteen inertial sensor units, each containing a 3D gyroscope and a 3D accelerometer (MPU-9150) were attached to the dorsal side of the right hand (3 sensors), to the metacarpal, proximal and distal phalanges of the thumb (3 sensors) and to the proximal, middle and distal phalanges of the digits (12 sensors). The sensors were secured in place using double-sided and single-sided medical tape. The 3 sensors on the dorsum of the hand, as well as all the sensors on the distal segments of the fingers were additionally equipped with 3D magnetometers. The data captured by the PowerGlove was used as dataset on which the performance of the three HPR algorithms have been evaluated.

2) Optical Tracking System
The optical motion tracking setup consisted of 6 infrared cameras (MX-13, VICON) arranged around a wooden table to create a relatively small capture volume appropriate for hand kinematic measurements. The arrangement and configuration of the cameras was optimized in terms of minimum coverage of the measurement volume and maximum marker visibility by adjusting the camera focus and IR-strobe intensity. Each subject's hand was equipped with thirty-six hemispheric 4mm retro-reflective markers and four spherical 9mm markers placed on top of the IMMUs (Figure 2), markers were placed on the IMMUs which ameliorated distortions due to e.g. tissue artifact. The kinematic data was recorded, reconstructed and labeled using VICON Nexus Software (version 1.8.5) with the video frame rate of 100 frames/second. To remove fluctuations of the markers and high-frequency drift, labeled trajectories were low-pass filtered (4th order Butterworth filter with 1.5Hz cutoff frequency for the static poses and 3Hz cut-off frequency for the dynamic pointing task) using the same software. The cut-off frequency was selected based on the fact that the performed hand poses were meant to be stationary, and on the earlier research pointing out that the bandwidth of most of the common activities of daily living is around 1Hz [8].

C. Hand Pose Reconstruction Algorithms
In this study, we tested three algorithms commonly used for HPR: Extended Kalman Filter, Gauss-Newton method combined with the Complementary Filter and Madgwick Filter, and comparatively evaluated them against the optical tracking system. Optical tracking systems are known for their high accuracy and are used as a gold standard for the HPR throughout the literature [1]. The performance of the chosen HPR algorithms was tested using a common dataset acquired with the experimental setup described above. See below for the description of the used HPR algorithms.

1) Extended Kalman Filter
EKF computes the state of the dynamic system in a twostep process. First, it estimates the current state variables, along with their uncertainties. Once the outcome of the next measurement (distorted by a measurement error) is observed, these estimates are updated using a weighted average (Kalman gain), with more weight being given to the estimate with the higher certainty. EKF uses multivariate Taylor Series expansion to linearize the model about current mean and covariance. In this study, a separate EKF has been used for each IMU sensor. The implemented filter uses a general state space model for system's dynamics and measurements. The state vector x consists of the single rotation quaternion (q) and tri-axial gyroscope bias vector ([ : The process model is based on the gyroscope integration and thus during the state update, state estimates are updated, based on the new gyroscope readings, according to: is the current quarternion estimate of the orientation, is the previous quarternion estimate of the orientation, is the quarternion multiplication operator, is the angular velocity measured at time k, is the sampling period, is the current gyroscope bias vector, is the previous gyroscope bias vector and is the change in the gyroscope bias.
The measurement model is constructed by stacking together the accelerometer and magnetometer measurement vectors. The a posteriori error state and corresponding covariance estimate, given the a priori state and covariance, is determined using the linearized matrix of the measurement function and the calculated Kalman gain. Finally, during the time update, error states are set to zero and the corresponding covariance matrix P is propagated according to discretized process model F and process noise Q.

2) Gauss-Newton combined with Complementary Filter
This algorithm exploits the quaternion-based method presented by Comotti et al. in [9]. It consists of two phases. In the first phase, accelerometer and magnetometer data are used to estimate the orientation. In the second phase, this estimation is fused with the orientation computed through integration of the gyroscope data.
The algorithm estimates the orientation using inertial and magnetic data, by minimizing a cost function. Exploiting the Gauss-Newton minimization procedure, the algorithm processes the measurement of gravity and Earth's magnetic flux to evaluate the current sensor orientation with respect to a fixed global reference frame. As suggested in [10], a compensation for magnetic distortion is performed.
The second step of the algorithm integrates the angular rate, measured by the gyroscope, and fuses the result with the output of the previous step using a complementary filter. It is worth pointing out that the gyroscope orientation is obtained using the last quaternion computed by the whole algorithm. Finally, the complementary filter merges together the outcome of the Gauss-Newton method and the angular rate integration. For further details refer to [11].

3) Madgwick Filter
The Madgwick Filter [10] is a constant gain filter adapted to estimate the orientation of a rigid body in quaternion form using data from the IMMU sensors. In this algorithm, the gyroscope measurements are first processed with a correction algorithm in order to minimize the effects due to the gyroscope bias and the drift error. They are then represented as quaternions, integrated and used to compute the absolute body orientation with the quaternion propagation starting from the orientation estimated at the previous timestamp. In the second step, the accelerometer and magnetometer measurements are fused together and converted to the quaternion representation through the gradient descent algorithm. The gradient descent algorithm aims to find a single, optimal sensor orientation solution that aligns the predefined direction of the magnetic field of the earth frame with the measured field in the sensor frame. The output of the gradient descent algorithm is then used to correct the orientation estimated by considering only the gyroscope measurements.

4) Filter Settings
In the simplified EKF implementation the process and measurement noise covariance matrices were initialized using the gyroscope and accelerometer noise values reported in the MPU 9150 datasheet [12], and gyroscope bias stability and magnetometer noise values reported in [13]. The Madgwick Filter was initiated with β = 0.041 and ζ = 0.02, as suggested by the author of the filter [2]. Note however that the initial β value was set to 2 for the first second of every trial, to ensure the convergence of the algorithm states from the initial conditions. Finally, the Complementary Filter was run with α = 0.85. The value of α was chosen experimentally, after some pilot experiments, as the value that allows for the attenuation of the gyroscope related drift without impacting the quality of the data.

D. Study Population & Design
The study design was approved by the internal review board of the faculty of Electrical Engineering, Mathematics and Computer Science of the University of Twente. The study protocol was executed at Laboratories of Roessingh Research and Development, Enschede, the Netherlands. Eight healthy right-handed adults (1 female, 7 males, age: 25.5±1.3 years, hand length: 19.4±1.1cm), studying at the University of Twente participated in this research. None of the subjects suffered from any movement disorders. Written informed consent was obtained from all subjects included in the study. Three types of experiments were performed, each addressing specific aspects of the HPR, as described below.
Repeatability: to investigate the repeatability of the estimation of the finger joint orientations, two hand postures were used: open hand posture and clenched hand posture. Thermoplastic material (ProtoPlast) was used to create two subject specific molds (one for each hand posture). A standardized protocol previously proposed by Wise et al. [14] has been partially adopted here to allow for easier comparison between this and similar studies. Starting with the open hand position, subjects were asked to place their hand flat in the mold and hold it still for 5 seconds, then to raise it up, flex all the fingers for 5 seconds, and return to the flat hand position. This flexion/flat cycle was repeated 10 times. Similarly, in the clenched hand position, subjects were asked to clench the mold for 5 seconds, then to release it and hold their hand still for 5 seconds. This clench/release cycle was repeated 10 times.
Static accuracy: to investigate the static accuracy of the algorithms, ten static hand poses were used. The included poses are shown in Figure 3 and are based on the GRASP taxonomy by Feix et al. [15] and the work of Ciotti et al. [4]. During the task subjects were asked to grab and hold real objects associated with each of the hand poses (see Figure 3). To avoid disturbances of the magnetic field no metal objects, except for the hammer head, were used, e.g. the key was 3D printed, the bag handle was textile. The hammer head was never located in the close vicinity of the IMMU sensors. Subjects performed each grasp 5 times maintaining the static hand pose for 5 seconds and placing their hand flat on the table for 5 seconds in between each two repetitions of the grasp. The order of the grasps was randomized across subjects.
Dynamic accuracy: to investigate algorithms' ability to keep track of rapid finger movements, a circular pointing task was included in the research protocol. This task was adopted from Noort et al. [7]. During this task, subjects made circular movements with a stretched index finger, keeping the hand still in a static posture such that the abduction/adduction angle of the index finger MCP joint was maximized. Five repetitions of this cyclical movement constituted a trial (one trial lasted 3 to 3.5 s). The task consisted of 10 trials.

E. Data Processing and Analysis 1) Optical Data
Labelled and filtered optical data was read into MATLAB using Biomechanical ToolKit. To correct for marker's flickering and occlusions, gaps in marker positions smaller than 20 frames were interpolated using a spline function. Only the frames in which all the finger markers were visible were analyzed further.
To calculate the flexion/extension and adduction/abduction angles, local reference system had to be established. Since the hand is not perfectly flat and curves in so called hand arch, two separate local coordinate systems were created (IM for index and middle fingers, and RL for ring and little fingers) and embedded at the sensors placed at the dorsum of the hand (Figure 2). Three markers placed at each of the hand sensors (H_IM_1, H_IM_2, H_IM_3 and H_RL_1, H_RL_2, H_RL_3, see Figure 2 for markers definitions) denoted its local reference frame as follows: • The x-axis is perpendicular to the plane containing H_IM_1, H_IM_2 and H_IM_3 markers (H_RL_1, H_RL_2 and H_RL_3 markers for RL coordinate system); • The y-axis is defined by the vector created by H_IM_1 and H_IM_2 markers (H_RL_1 and H_RL_2 markers for RL coordinate system); • The z-axis, is formed by the cross product of the x-and y-axes. To avoid magnetic field disturbance, extra care was taken for the used objects not to contain metallic elements. Thus, the key used in this study was 3D printed from a plastic material, and the metal needle from the pin, used for the pinch grasp, was removed. The only object that contained metal elements was a hammer. The head of the hammer was never located in the close vicinity of the IMU sensors.

• Circular pointing task
In order to investigate algorithms' ability to keep track of rapid finger movements, a circular pointing task was included in the research protocol. This task, which has previously been employed in the validation study of the PowerGlove [100], nicely compliments the assessment of the static reconstruction of the hand poses and allows us to see how different algorithms behave in more dynamic tracking settings. During this task, participants were asked to make circular movements with stretched index finger while the hand maintained a static posture such that the abduction/adduction angle of the MCP joint was maximized. Five repetitions of the cyclical movement constituted a trial. The task consisted of 10 trials.

Study Population
Thirteen healthy right-handed adults (2 females, 11 males, age: 25.08±1.65 years, hand length: 19.6±1 cm), studying at the University of Twente participated in this research. None of the participants suffered from any movement disorders.
Written informed consent was obtained from all participants included in the study. Next, for each finger, three vectors, one for each segment, were defined using the markers placed on the corresponding IMMU units. For each of the triangle marker clusters at the dorsum of the hand, three base vectors were constructed using the markers positions, the geometric center and the repeated cross-products. The angles between the vectors projected onto the local sagittal plane (IM for index and middle finger, RL for ring and little finger) described the flexion/extension movement. Similarly, the angles between the vectors projected onto the local coronal plane described the adduction/abduction of the MCP joints. Joint angles were measured for each of the static positions, dynamic circular pointing task as well as for the baseline flat hand position.
The PowerGlove data was read and processed using custom made MATLAB scripts. The data was first filtered using the 4th order Savitzky-Golay filter with a window length of 15 samples for the accelerometer and 31 samples for the gyroscope. Next, sensor calibration and anatomical segment calibration (see below) were applied.
As the MEMS based IMMUs often suffer from nonaccurate scaling, sensor axis misalignment and non-zero biases, IMMU sensor calibration was performed and used to correct the recorded PowerGlove data. The employed sensor calibration procedure was adapted from the work of Tedaldi et al. [16]. It is based on the assumption that in a static position, the norm of the measured acceleration is equal to the magnitude of the gravity acceleration plus a multi-source error factor (including bias, misalignment, noise, etc.). All those can be estimated via minimization over a set of static attitudes. After the calibration of the accelerometer, the gravity vector position measured by the accelerometer can be used as a reference to calibrate the gyroscope. Integrating the angular velocities between two consecutive static positions, we can estimate the gravity vector position in the new orientation. The gyroscope calibration is obtained by minimizing the errors between these estimates and the gravity reference given by the calibrated accelerometer. To this extend, to calibrate the IMMUs, the PowerGlove sensors were secured in a flat position within a tightened frame. The sensors were left running inside the frame for about 5 minutes for the temperature of the IMMUs to stabilize. The frame was then placed in about 40-50 different static positions to obtain the scale, misalignment and bias parameters of all accelerometers and gyroscopes. Those parameters were then used to calibrate the recorded data.
Magnetometer measurements, which often suffer from hard and soft iron artifacts were calibrated using the ellipsoid fitting method [17]. This method is based on the fact that the error model of the magnetic compass is an ellipsoid and thus a constraint least-square algorithm can be adopted to estimate the parameters of this ellipsoid by rotating the magnetic compass in various (random) orientations. Those parameters can then be used to compensate for the artifacts and transform the measurement surface into a sphere.
Furthermore, since during the experiments PowerGlove sensors were attached to the hand and finger segments, anatomical calibration procedure was performed to align sensor coordinate frames with the finger's segments' coordinate frames. For this reason, during each experimental session additional subject-specific calibration dataset was acquired. First, the subjects were asked to place the hand flat on the table, in such a way that the gravity vector was perpendicular to the palm of the hand, the fingers in neutral position with the phalanges aligned to each other, and the thumb in abduction. Secondly, the thumb was placed on the table with the dorsal side of the thumb and the nail positioned horizontally, in such a way that the gravity direction was perpendicular to the long axis of the thumb. The x-axes of the hand and finger segments, were determined from the accelerometer output during static postures of hand and thumb: (4) To find the z-axes of the finger segments, subjects were asked to flex the thumb couple of times in the IP joint, place the hand flat on the table and subsequently flex all the fingers in the MCP joints. A little book was used as a support for this task. The z-axes were then determined from the gyroscope output during the flexion/extension movements: The sensor to segment transformation was then defined as: (6) Finally, the hands were placed together and moved in an eight-shaped movement for five seconds. Angular velocity measured on various hand segments were assumed equal during this movement but measured in different coordinate frames. Relative orientation was deduced, to express the signals of all units on the hand (thumb, index/middle and ring/little) in a common reference frame.
Calibrated PowerGlove data was fed to the three algorithms. The algorithms returned the orientation of each of the finger segments in quaternion form. Obtained quaternions were transformed into Euler angles. Due to problems with magnetometer data calibration/quality, the analysis was restricted to flexion/extension angles that were calculated as a difference in roll angles of the two adjacent finger segments.

3) Outcome Parameters and Statistical Analysis
To assess the repeatability of the estimated finger joint angles in both flat and clenched positions, following the original protocol of Wise et al. [14], two different measures were calculated for each DoF (joint angle). Those measures were: • range: refers to the difference between the highest and lowest joint angle value measured during the flat hand/clenched hand phase, in degrees • standard deviation of the joint angle values measured during the flat hand/clenched hand phase, in degrees Those measures were then averaged across all joints, separately for each subject and each algorithm. To quantify the computational complexity of the compared algorithms, the repeatability datasets were also used to measure the time needed to obtain orientation estimation of each finger segment, separately for each subject and each algorithm.
To compare the performance of the chosen algorithms against the optical tracking system, the joint angles measured during the static hand poses experiment were analyzed. The algorithm performance was evaluated in terms of the estimation errors as compared to the ground truth results obtained with the VICON optical tracking system. Absolute DoF estimation errors, that is the absolute differences between joint angle values returned by the algorithm and the ground truth measurement provided by the VICON system: (7) and pose estimation errors, that is the mean of the DoF absolute estimation errors computed separately for each pose: (8) were considered and averaged over all of the trials [18]. These metrics were computed separately for each subject and each algorithm. Due to problems with optical markers' visibility, bag grasp and hammer grasp were excluded from the analysis. Likewise, due to limited visibility of the markers, thumb data was not analyzed. In case of some subjects and grasps, optical data of the ring and/or little finger was also missing. Those trials were thus not used for the calculation of the pose estimation error, and were included only in the calculation of the DoF absolute errors.
In order to evaluate the algorithm's performance during dynamic finger movements, the circular pointing task data was analyzed. The estimated flexion/extension angle of the MCP joint of the index finger was compared to the measurement output of the optical tracking system. Posthoc cross-correlation of the calculated joint angles was used to synchronize the PowerGlove data with the optical data. The accuracy of the algorithm was then estimated using root square mean (RMS) difference between both measurement systems averaged over all trials. This metric was computed separately for each subject and algorithm.
All statistical analyses were performed with SPSS (Statistical Package for Social Science, IBM, Armonk, New York, USA) with significance level set at α = 0.05. At first the data was checked for normality using Shapiro-Wilk test. All the data, apart from the computational times data, proved to be distributed in the Gaussian fashion, allowing for the parametric statistics to be used. To assess the effect of the algorithm on the joint angles repeatability, static pose estimation errors and dynamic trial RMS differences, repeated measures ANOVAs were conducted separately for each metric of interest. To assess the differences in computational complexity, nonparametric Friedman test was employed. In cases where the assumption of sphericity was not met, Greenhouse-Geisser correction was used. Post-hoc analyses were carried out with Bonferroni correction for multiple comparisons (3 comparisons -significance level α = 0.017). To support the statistical interferences drawn based on the p-values, subsequent Bayesian statistical analysis was carried out using JASP [19]. Bayesian analysis allowed for quantification of the evidence in favor of the null hypothesis (H0: no difference between the algorithms). To this extend Bayesian repeated measures ANOVAs were used and followed by post-hoc Bayesian paired t-tests. The prior was set to follow Cauchy's distribution with the width equal to 0.707 [20]. To assure the robustness of the statistical test to the choice of the prior's width, the same analysis was repeated for prior widths ranging from 0 to 1.5. Since the choice of prior width did not influence the results, only the Bayesian Factors associated with Cauchy's prior width of 0.707 are reported here.

A. Calibration
The sensor calibration corrected sensor scaling, misalignment and bias errors in the accelerometers, gyroscopes and magnetometers. Sensor to segment calibration measurements were carried out for each subject, allowing transformation of sensor data from sensor frame to segment frame.

B. Repeatability Results
To compare repeatability of the algorithms' outcomes, the range (mean difference between highest and lowest joint angle value obtained in each trial) and standard deviation of those joint angles have been calculated across n=10 trials, separately for each subject and each algorithm. Withinsubjects comparisons allowed us to look at the differences between algorithms while accounting for intra-subject variability.  To look for significant differences between algorithms, repeated measures ANOVA and Bayesian repeated measures ANOVA were used. In the flat hand scenario, the repeatability of the algorithms was found to be significantly different between the algorithms (Figure 4) in terms of range (F(2,14) = 8.885; p=0.003, BF=10.64) but not in terms of standard deviation. Post-hoc pairwise comparison with Bonferroni adjustment (α = 0.017) revealed that the repeatability range of the CF was significantly smaller (p<0.005) than the same metric calculated for the two other algorithms. This result was further confirmed by the Bayesian analysis, where for range CF<Madgwick: BF=10.62; CF<EKF: BF=13.84; Madgwick≠EKF: BF=0.33. No significant differences were found in the clenched hand scenario, with the Bayesian analysis providing substantial evidence in favour of the null result (no differences between algorithms) for the clenched hand repeatability range (BF=0.25) ( Figure 5).

C. Static Hand Poses
The estimation accuracy of the algorithms was evaluated in terms of the estimation errors, as compared with the optical tracking system. Both individual DoF (joint angle) estimation errors and pose estimation errors (average of all joint angle estimation errors for a given grasp, see Methods section) were considered ( Figure 6). Statistical analysis revealed significant effect of the hand pose on the pose estimation error (F(7,14)=2.74, p=0.017, BF>3000), with glass and open hand poses associated with the lowest reconstruction accuracy. No significant effect of the algorithm was found, with Bayesian analysis providing support for the null hypothesis (no differences between the algorithms, BF=0.27). Similarly, no significant interaction between the pose and the algorithm has been observed (BF=0.012).

D. Circular Pointing Task
To assess algorithms' performance in more dynamic conditions, a circular pointing task was employed (see Methods). Index finger MCP flexion/extension angle was calculated using the three algorithms, as well as optical tracking data. RMS errors between the output of the algorithms and the VICON data were calculated and averaged over all trials (Figure 8). The CF algorithm was on average associated with 8.34±2.64 degrees RMSE, Madgwick algorithm with 8.74±3.04 degrees RMSE and EKF with 8.46±2.14. No statistical difference between the algorithms was found (BF=0.44). Figure 8. A) Representative reconstruction of the index finger's MCP joint flexion/extension angle during the circular pointing task. Data from one trial and one subject. Madgwick algorithm's reconstruction is plotted in blue, CF's in red and EKF's in yellow. Reference data from the Vicon system is plotted as a dashed black line. B) Average RMSE, as compared to the optical tracking data, for each of the algorithms during the circular pointing task. Error bars represent standard deviation

E. Computational Complexity
The computational complexity of the algorithms was assessed by measuring the CPU time needed to compute the segment orientations from the repeatability tasks. The computer system used was a MacBook Pro, 2.66GHz Intel Core 2 Duo processor, 8GB RAM, OS El Capitan. Results of the Friedman test indicated significant statistical differences between computation times obtained with different algorithms (χ 2 (2,574)=472.11, p<0.0001, BF> 3000). Posthoc analysis proved Madgwick algorithm to be the fastest (on average 3.18±1.91s), EKF algorithm to be the second best (on average 4.36±0.42s), and CF algorithm to be the slowest one (on average 17.47±3.34s). Those results were corroborated by the Bayesian post-hoc testing, with BF always above 10.15.

IV. DISCUSSION
This work provides an overview of the performance of the three most commonly employed IMMU-based HPR algorithms. Outcomes of the algorithms are comparable, however MF outperforms the other algorithms in terms of accuracy in individual DOF reconstruction and reduced computational complexity. An overview of the main results is given in Table II.
It is worth mentioning that the outcomes presented here are dependent on the choice of the tuning parameters of the filters. If different parameters will be chosen, the results of the comparison might also be different. Furthermore, it needs to be remembered that the EKF implementation used in this study was very much simplified when compared to the more elaborate EKF implementations proposed in other HPR studies. With only seven states included in the states vector (see Methods section), the used implementation did not take full advantage of EKF's probabilistic nature. The original implementation of Kortier et al. [21], with larger state vector, biomechanical constraints update and hierarchical overlay of multiple Kalman Filters would likely perform better than its simplified counterpart. However, it would also likely increase the computation load, which may become an issue depending on the context in which the algorithm is being used. For instance, if HPR is used for therapeutic assessment, feedback based training or when integrated as sensory pathway in a humanoid, real-time HPR processing, and hence low computational load, would typically be required. To compare the algorithms, an optical tracking system has been used here as the ground truth reference. Internal consistency of the HPR in an optical and IMMU-based systems has been investigated in [7] with inaccuracies attributed to the skin movement artifacts. However, in this study optical markers were placed directly on top of the IMMU sensors, so skin movement artifact affected both measurement systems in the same way.

V. CONCLUSIONS
Objective and accurate analysis of hand kinematics plays a critical role in many fields including clinical practice, robotics, rehabilitation and neuroprosthetics. In light of many different strategies used to reconstruct the pose of the hand, the aim of this study was to compare three different HPR algorithms. The three algorithms compared in this study obtained very similar scores on most of the considered metrics, suggesting no gross differences in their HPR performance. Yet, the MF algorithm proved to be slightly more accurate in reconstructing the individual joint angles during static trials (its static absolute DoF error was significantly lower than those of the CF and EKF algorithms) and was the fastest one to run. In future studies it may also be interesting to look into internal consistency of reconstruction (ICC), and kinematic hand synergies as another potential comparison metrics. As muscle synergies have recently been shown to be affected by the pre-processing pipeline [22], investigating robustness of the kinematic synergies to the choice of the sensor fusion algorithm could be an important issue to tackle. Hand synergies could also be used to study minimal sensing approaches [1,18], where the posture of the hand is deducted from an incomplete and imperfect subset of data. In other words, it can be used to examine the predictability of a subset of joint movements in respect of the movements of other joints. This is of fundamental interest in the human hand tracking, e.g. in the rehabilitation-assessment field, where the need for wearability imposes constraints to the number of sensors. Towards this goal, dynamic accuracy of the here discussed algorithms will be evaluated in dailyliving tasks, trying to minimize marker occlusions that may occur when acquiring optical-based ground truth data. Finally, to generalize the results of this work, testing these techniques with other hardware platforms is also envisioned.