1 Introduction

Most of the applications for 3D User Interfaces in Augmented Reality (AR) and Virtual Reality (VR) will require the user to interact with immersive environments. One of the key goals of VR and AR systems is to foster the perception of “presence” in the user, which strives to provide the user with highly realistic visual and acoustic output [1,2,3,4]. In some VR environments, the perception of presence may be severely disrupted when the user is asked to interact using highly “unnatural” sequences of actions performed on an input device (e.g., joystick, gamepad, mouse or keyboard). Thus, a full hand motion tracking system would significantly help in circumventing these current limitations and improve 3D User Interfaces to become more realistic. Our proposed hand motion tracking interface seeks to determine the movement of the human hand by combining two different sources of information: orientation tracking using Inertial Measurement Units (IMUs) and position tracking using infrared cameras: OptiTrack V120 Trio as shown in Fig. 1.

Fig. 1.
figure 1

Position tracking of hand motion visualized in unity

Our research has yielded an algorithm to minimize the gyroscope drift problem which commonly occurs in miniature commercial-grade IMUs. IMUs contain gyroscopes, which provide angular velocity measurements that can be used to track orientation of the object that they are attached to. However, orientation tracking also requires the mathematical integration of angular velocity readings from the gyroscopes. Most low-cost gyroscopes may output an erroneous non-zero measurement level even when they are not moving, called the bias offset error. This bias offset error can produce an orientation tracking error called drift, which causes many problems in navigation and other applications that use low-cost gyroscopes to determine the orientation of moving parts [5, 6]. Several studies and algorithms have been proposed to eliminate gyroscope drift error in inertial measurement systems. Most of the studies use Kalman-based processes to solve the error in inertial sensors [7,8,9]. But this approach is complex and can be complicated to implement [10, 11].

To estimate the orientation, many studies use the idea of sensor fusion which is the method that combines more than one sensor to determine one measurement. For example, a sensor fusion approach may suggest using gyroscopes, accelerometers and magnetometers to determine orientation [12, 13]. The idea of sensor fusion is also the main method to determine the orientation in this research. This paper outlines the use of this algorithm to improve orientation tracking of human hand motion by using the combination of gyroscope, accelerometer and magnetometer measurements.

2 Methodology

2.1 Bias Offset Estimation

To minimize the drift, an algorithm to correct the gyroscope bias offset has been proposed. With the idea that the IMU sensor should provide the reading of zero when there is no input or movement taking place, the bias offset error could be determined and used to subtract from future gyroscope readings, to obtain unbiased gyroscope data as a result. The predicted bias offset error \( (\widehat{b}) \) of each gyroscope axis is determined during the sensor’s static periods by using a simple linear regression model [14] as shown in Eq. (1). The model coefficients \( (\beta_{1} \,{\text{and}}\,\beta_{0} ) \) can be found by using Eqs. (2) and (3), where \( b_{i} \) and \( t_{i} \) are the measured gyroscope bias and time at the ith sample, respectively.

$$ \widehat{b} = \beta_{0} + \beta_{1} t $$
(1)
$$ \beta_{1} = \frac{{\sum\nolimits_{i = 1}^{n} {\left( { t_{i} - \bar{t} } \right)\left( {b_{i} - \bar{b}} \right)} }}{{\sum\nolimits_{i = 1}^{n} {\left( { t_{i} - \bar{t} } \right)^{2} } }} $$
(2)
$$ \beta_{0} = \bar{b} - \beta_{1} \bar{t} $$
(3)

The predicted bias is subtracted from the measured gyroscope reading, yielding the unbiased angular rate \( (\omega_{B} ) \). To avoid gimbal lock problems, the orientation representation used in this project is quaternion notation, placing the real part in the forth component. The quaternion rate \( (\dot{q}) \) can be calculated by using the unbiased angular rate \( (\omega_{B} ) \), where \( \widehat{q}_{0} \) is the previous quaternion estimation of the orientation (For the first iteration of an algorithm, \( \widehat{q}_{0} \) is set to [0, 0, 0, 1], indicating no rotation or rotation by angle of zero). An unbiased angular rate vector \( \omega_{B} \) which is in \( {\mathbb{R}}^{3} \) can be mapped into quaternion space \( ({\mathbb{R}}^{4} ) \) by augmenting zero as its real part, this quaternion is called a pure quaternion [15]. Since the real part of \( \omega_{B} \) in \( {\mathbb{R}}^{4} \) is zero, the quaternion multiplication can be written in matrix form as shown in Eq. (4).

$$ \dot{q} = \frac{1}{2}\widehat{q}_{0} \otimes \omega_{B} = \frac{1}{2}\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & {\omega_{Bz} } \\ { - \omega_{Bz} } & 0 \\ \end{array} } & {\begin{array}{*{20}c} { - \omega_{By} } & {\omega_{Bx} } \\ {\omega_{Bx} } & {\omega_{By} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\omega_{By} } & { - \omega_{Bx} } \\ { - \omega_{Bx} } & { - \omega_{By} } \\ \end{array} } & {\begin{array}{*{20}c} 0 & {\omega_{Bz} } \\ { - \omega_{Bz} } & 0 \\ \end{array} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\widehat{q}_{0x} } \\ {\widehat{q}_{0y} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\widehat{q}_{0z} } \\ {\widehat{q}_{0w} } \\ \end{array} } \\ \end{array} } \right] $$
(4)
$$ q_{G} = \exp \left( {\left( {\Delta t} \right)\dot{q} \otimes \widehat{q}_{0}^{*} } \right) \otimes \widehat{q}_{0} $$
(5)

The quaternion representation of the orientation \( (q_{G} ) \) can be obtained by integrating the quaternion rate using Eq. (5), where \( \Delta t \) is the sampling period of the sensor’s data. This quaternion \( (q_{G} ) \) can be used to rotate points or vectors in three-dimensional space. In other words, it can be used to describe the orientation of the sensor module or the object that the IMU is attached to. To justify the orientation estimation using this method, we also proposed the method to compare this quaternion result with other referencing measurements from the accelerometer and magnetometer integrated in the same IMU.

2.2 Quaternion Correction

Using the Gravity Vector.

Each of the three accelerometers in the IMU provides the value of acceleration applied to the sensor module in its sensing direction. When the sensor is static, the accelerometers are assumed to measure only the acceleration due to gravity, which is always directed towards the Earth’s center and has only a vertical component in the Earth frame. When the IMU module is in an oblique orientation (not at right angle or parallel with the horizon), the acceleration due to gravity is decomposed into each orthogonal axis in the sensor frame. These components can be represented as the inclination of the sensor module in the Earth frame. The triple-product quaternion operator shown in Eq. (6) is used to rotate the coordinate frame [15] with respect to vector \( \vec{v} \) into another coordinate frame represented by vector \( \vec{w} \) since \( q \) represents the orientation or the inclination of the sensor with respect to the Earth frame.

$$ \vec{w} = q^{ *} \otimes \vec{v} \otimes q $$
(6)

Accordingly, we can use the quaternion calculated from the integration of gyroscope reading from the first part \( (q_{G} ) \) to rotate the referencing frame of the gravity vector from the Earth frame into the sensor frame in order to obtain the gravity vector referencing in the sensor frame as shown in Eq. (7). This calculated gravity vector \( (\vec{a}\left( {q_{G} } \right)) \) will further be used to compare with the measured gravity vector from the accelerometer readings.

$$ \vec{a}\left( {q_{G} } \right) = q_{G}^{ *} \otimes A_{int} \otimes q_{G} $$
(7)

If the estimated orientation of the IMU, represented by the quaternion \( q_{G} \) had no error, \( \vec{a}\left( {q_{G} } \right) \) from Eq. (7) would match the measured gravity vector from accelerometer readings. Otherwise, the error between these two vectors can be used to correct the orientation estimation by calculating the angular difference between these two vectors. The difference in angle of the calculated gravity vector \( (\vec{a}\left( {q_{G} } \right)) \) and measured gravity vector from accelerometer readings \( (\vec{a}_{0} ) \) can be calculated and then represented as a rotation in quaternion form, as in Eq. (8). Then, the estimated quaternion using gravity vector correction \( (\widehat{q}_{GA} ) \) is calculated by multiplying the quaternion from Eq. (5) with the difference in quaternion \( (\Delta q_{A} ) \) between the calculated gravity vector and the measured gravity vector from accelerometer readings, from Eq. (8). The estimated quaternion \( (\widehat{q}_{GA} ) \) is described as in Eq. (11).

$$ \Delta q_{A} = {\mathcal{H}}\left( {\vec{q}_{Av} ,q_{Aw} } \right) $$
(8)
$$ \vec{q}_{Av} = \vec{a}_{0} \times \vec{a}\left( {q_{G} } \right) $$
(9)
$$ q_{Aw} = \left\| {\vec{a}_{0} } \right\|\left\| {\vec{a}\left( {q_{G} } \right)} \right\| { + }\vec{a}_{0} \cdot \vec{a}\left( {q_{G} } \right) $$
(10)
$$ \widehat{q}_{GA} = q_{G} \otimes \Delta q_{A} $$
(11)

This resulting quaternion \( \widehat{q}_{GA} \), which estimates the orientation of the sensor module by using gravity vector correction, will be used in the next step of the algorithm. Note that the quaternion result has to be normalized before proceeding with any calculations.

Using the Magnetic North Vector.

Similar to quaternion correction using gravity vector, this part of the algorithm pursues a similar approach in using the other sensing source as the referencing vector. Unlike the gravity vector, the directions of the magnetic north vector at different locations on Earth are different. However, at the particular location, the magnetic north vector is always pointing towards the North magnetic pole. Therefore, each of the three magnetometers in the IMU should provide the value of magnetic field strength applied to the sensor module in its sensing direction. Whether the sensor is static or moving, the magnetometers are assumed to measure the direction of the magnetic North vector, which is always directed towards the North magnetic pole in the Earth’s coordinate frame. When the IMU module is deviated from magnetic north vector direction (not parallel with the direction of the Earth magnetic field), the magnetic North vector is decomposed into each orthogonal axis in the sensor frame. These components can also be represented as the inclination of the sensor module in the Earth frame. Thus, we can use the quaternion calculated from the integration of gyroscope reading from the first part \( (q_{G} ) \) to rotate the referencing magnetic North vector in the Earth frame \( (M_{int} ) \) in order to obtain the magnetic North vector in the sensor frame as shown in Eq. (12). This calculated magnetic North vector \( (\vec{m}\left( {q_{G} } \right)) \) will further be used to compare with the measured magnetic North vector from the magnetometer readings.

$$ \vec{m}\left( {q_{G} } \right) = q_{G}^{ *} \otimes M_{int} \otimes q_{G} $$
(12)

If the estimated orientation of the IMU, represented by the quaternion \( q_{G} \) had no error, \( \vec{m}\left( {q_{G} } \right) \) from Eq. (12) would match the measured magnetic North vector from magnetometer readings. Otherwise, the error between these two vectors can be used to correct the orientation estimation by calculating the angular difference between these two vectors. The difference in angle of the calculated magnetic North vector \( (\vec{m}\left( {q_{G} } \right)) \) and measured magnetic North vector from magnetometer readings \( (\vec{m}_{0} ) \) can be calculated and then represented as a rotation in quaternion form, as in Eq. (13). Then, the estimated quaternion using gravity vector correction \( (\widehat{q}_{GM} ) \) is calculated by multiplying the quaternion from Eq. (5) with the difference in quaternion \( (\Delta q_{M} ) \) between the calculated magnetic North vector and the measured magnetic North vector from magnetometer readings, from Eq. (13). The estimated quaternion \( (\widehat{q}_{GM} ) \) is described as in Eq. (16).

$$ \Delta q_{M} = {\mathcal{H}}\left( {\vec{q}_{Mv} ,q_{Mw} } \right) $$
(13)
$$ \vec{q}_{Mv} = \vec{m}_{0} \times \vec{m}\left( {q_{G} } \right) $$
(14)
$$ q_{Mw} = \left\| {\vec{m}_{0} } \right\|\left\| {\vec{m}\left( {q_{G} } \right)} \right\| + \vec{m}_{0} \cdot \vec{m}\left( {q_{G} } \right) $$
(15)
$$ \widehat{q}_{GM} = q_{G} \otimes \Delta q_{M} $$
(16)

The resulting quaternion \( \widehat{q}_{GM} \) estimates the orientation of the sensor module by using magnetic North vector correction and will proceed to the next step of the algorithm. Note that the quaternion result has to be normalized before proceeding with any calculations.

From Eq. (11), \( \widehat{q}_{GA} \) is the product of quaternion multiplication between \( q_{G} \) and \( \Delta q_{A} \) which represents the combination of both rotations. Since \( \Delta q_{A} \) is the rotation that represents angular difference in quaternion form between the calculated gravity vector and the measured gravity vector from accelerometer readings, the quaternion product of \( q_{G} \otimes \Delta q_{A} \) then represents the estimated orientation that is fully corrected based on the measurement of accelerometer. Similarly, \( \widehat{q}_{GM} \) also represents the estimated orientation that is fully corrected based on the measurement of magnetometer.

2.3 Quaternion Interpolation

With the combination of bias offset estimation and quaternion correction using gravity vector, the algorithm can be used to correct the drift in rotation and improve the orientation estimates [16]. However, when a rapid movement of the hand occurs, the accelerometer measures not only just acceleration due to gravity but also the acceleration due to linear motion. Thus, the accelerometer readings can no longer be used as a reliable source of reference. In contrast, the magnetometer readings are not affected by the linear acceleration due to rapid movement of the hand (However, the magnetic field perceived by the sensor may not be uniform around the user, due to the presence of ferromagnetic objects nearby). The periods of time when the hand is more likely to have rapid translation and linear acceleration inclusion are the moments when the orientation estimates should rely more on the orientation estimates using magnetometer \( (\widehat{q}_{GM} ) \) than the orientation estimates using accelerometer \( (\widehat{q}_{GA} ) \). An appropriate method to partially or selectively combine two different rotations by using a parametric weight is the method used in computer graphic animation called Spherical Linear Interpolation (SLERP) or great arc interpolation [17].

$$ SLERP(q_{0} ,q_{1,} h) = \frac{{q_{0} \sin \left( {\left( {1 - h} \right)\Omega } \right) + q_{1} { \sin }\left( {h\Omega } \right)}}{{{ \sin }\left(\Omega \right)}} $$
(17)
$$ { \cos }\left(\Omega \right) = q_{0} \cdot q_{1} $$
(18)

Spherical Linear Interpolation of two quaternions as shown in Eq. (17) is the method in animation to calculate the quaternion or the rotation representation in the intermediate stages between the starting \( (q_{0} ) \) and the ending \( (q_{1} ) \) of rotation sequence. The control parameter \( (h) \) that interpolates the rotation ranges from 0 to 1, where the control parameter that tends towards zero interpolates the rotation close to the starting rotation, and the control parameter that tends towards one interpolates the rotation close to the ending rotation. The control parameter divides the interpolated rotation linearly. Evidently, the Spherical Linear Interpolation preserves the normalization of the quaternion and all the interpolated rotations lie on the same great arc.

In Eq. (19), the Spherical Linear Interpolation method was adapted to determine the final estimated orientation \( (\widehat{q}_{OUT} ) \) between the orientation estimate using magnetometer \( (\widehat{q}_{GM} ) \) and the orientation estimate using accelerometer \( (\widehat{q}_{GA} ) \). The angle of the interpolated arc subtended by \( \widehat{q}_{GM} \) and \( \widehat{q}_{GA} \) is denoted by \( \Omega \), where \( \cos \left(\Omega \right) \) equals to the 4-dimensional inner product of \( \widehat{q}_{GM} \) and \( \widehat{q}_{GA} \). The control parameter of quaternion interpolation \( (\alpha ) \) is the parameter that indicates the stillness of the sensor module. When \( \alpha \) equals to 1, it indicates a static period of the sensor module. Thus, the final estimated quaternion orientation is solely estimated using accelerometer data \( (\widehat{q}_{GA} ) \) since the accelerometer is assumed to measure only the acceleration due to gravity when the sensor is static. When there is a movement, the accelerometer measurement includes acceleration due to linear motion, resulting in the drop of stillness. Consequently, the control parameter of quaternion interpolation \( (\alpha ) \) is close to zero. In this case, the final estimated quaternion is the interpolated orientation estimate that tends towards the orientation estimates using magnetometer \( (\widehat{q}_{GM} ) \).

$$ \widehat{q}_{OUT} = \frac{{\widehat{q}_{GM} \sin \left( {\left( {1 - \alpha } \right)\varOmega } \right) + \widehat{q}_{GA} \sin \left( {\alpha \varOmega } \right)}}{{{ \sin }\left( \varOmega \right)}} $$
(19)
$$ \cos \left(\Omega \right) = \widehat{q}_{GM} \cdot \widehat{q}_{GA} $$
(20)

This algorithm to estimate the orientation in quaternion form is also visually described in the block diagram shown in Fig. 2. The final orientation estimates \( (\widehat{q}_{OUT} ) \) is also fed back for the calculation of the quaternion estimates in the initial stage \( (q_{G} ) \) for the next iteration.

Fig. 2.
figure 2

Block diagram of the proposed drift correction algorithm using both (a) gravity vector and (b) Magnetic North vector correction

3 Implementation

In this work, we used two 3-Space Sensor modules, from Yost Labs. This module is an ultra-miniature, high-precision, high-reliability, low-cost Inertial Measurement Unit (IMU). The module has the dimensions of 23 mm × 23 mm x 2.2 mm (0.9 × 0.9 × 0.086 in.) and consists of tri-axial gyroscope, accelerometer, and magnetometer. The coordinate axes of the sensor’s frame are superimposed on the image of the sensor’s module shown in Fig. 3.

Fig. 3.
figure 3

Coordinate axes superimposed on the image of YEI 3-Space Embedded IMU module

These Yost Labs 3-Space Sensor modules were attached on the back of the hand and on the tip of the index finger as shown in Fig. 4. The gyroscope, accelerometer and magnetometer data were recorded and stored in a text file while a sequence of hand motions were being performed. The orientation correction algorithm was applied to the recorded data within MATLAB. In order to be able to effectively apply this orientation correction algorithm to the hand motion tracking system, the algorithm are supposed to work in real-time. Therefore, during the implementation of the algorithm in MATLAB, the programming code was aware of only one sample of data at a time in order to imitate real-time performance.

Fig. 4.
figure 4

The glove with IMUs attached on the back of the hand and on the tip of index finger

The bias offset error must be estimated when the IMUs detect periods of stillness. A flowchart showing the implementation of the orientation correction algorithm is displayed in Fig. 5. It shows that the algorithm begins with the calculation of the maximum value in a 50-sample window average of gyroscope data (from the 3 axes) called gyroMaxAvg.

Fig. 5.
figure 5

Flowchart showing the implementation of orientation correction algorithm for one iteration using both gravity and magnetic North vector corrections.

If this value is less than a pre-defined threshold value (gyroThreshold), that IMU is considered to be in the static period (this condition has to be true for 25 consecutive samples). Then, the new predicted bias offset error will be calculated. Otherwise, the previous predicted bias offset error will be used to calculate the unbiased angular rate instead. The quaternion \( q_{G} \) is calculated by integrating the unbiased angular rate. The quaternion correction using gravity-vector and magnetic North vector are implemented simultaneously. The output quaternion \( (\widehat{\varvec{q}}_{{\varvec{OUT}}} ) \) is calculated by using quaternion interpolation between estimated quaternion correction using gravity vector \( (\widehat{\varvec{q}}_{{\varvec{GA}}} ) \) and estimated quaternion correction using magnetic North vector \( (\widehat{\varvec{q}}_{{\varvec{G}M}} ) \). In our experiments, the resulting orientations were exported as a text file and visualized by the animation of a 3D hand model in Unity. The motion of the hand model with corrected orientation was validated and compared with snapshots of the hand motion sequence.

4 Results and Discussion

Figures 6(a) and 7(a) show the orientation in quaternion for IMUs attached on the back of the hand and the index finger. They were recorded directly from the sensor modules for a duration of 22.348 s. These recordings are the quaternion output from the Kalman-based orientation filtering inside the Yost Lab’s 3-space sensor module.

Fig. 6.
figure 6

The plots of estimated quaternion results (a) with On-board Kalman-based Orientation Filtering and (b) with gravity and magnetic North vectors correction for IMUs attached on the back of the hand

Fig. 7.
figure 7

The plots of estimated quaternion results (a) with On-board Kalman-based Orientation Filtering and (b) with gravity and magnetic North vectors correction for IMUs attached on the index finger

The resulting estimated orientation after gravity-vector and magnetic North vector correction for IMUs attached on the back of the hand and the index finger are shown in Figs. 6(b) and 7(b), respectively. Comparing the quaternion results from the orientation correction algorithm using gravity and magnetic North vector and the quaternion recording from the on-board Kalman-based filtering in Fig. 6, it was found that the result from the algorithm using gravity and magnetic North vector correction can reduce the amount of drift that distinctly occurred in Kalman-based filtering estimation during the 2nd to 5th seconds and 17th to 21st seconds. Similarly in Fig. 7, the quaternion result of the IMU attached on the index finger from the orientation correction algorithm using gravity and magnetic North vector and the quaternion recording from the on-board Kalman-based filtering are compared. This figure shows that the result from the algorithm using gravity and magnetic North vector correction can also reduce the amount of drift in the quaternion result using Kalman-based filtering method that occurred during the 7th to 11th seconds, 12th to 14th seconds and 17th to 21st seconds of the recording.

To verify these results, the quaternion outputs obtained from the gravity and magnetic North correction algorithm were exported as text files and were applied to the rotation of the 3D hand model in Unity. Figure 8 shows the comparison between a sequence of actual photographs of the hand orientation and 3D visualizations of estimated hand orientation after gravity-vector and magnetic North vector correction. Comparing the results to the actual hand orientation, leads us to confirm that the estimated quaternion after gravity-vector and magnetic North vector correction provides an acceptable result of orientation tracking using the IMUs. The 9-stage sequences of 3D hand orientation in Fig. 8 are the visualization of the quaternion results in Figs. 6(b) and 7(b). The numbers (1–9) to the left of the stages in Fig. 8 correspond with the (1–9) intervals labelled at the bottom of Figs. 6(b) and 7(b).

Fig. 8.
figure 8

Comparison between a sequence of actual hand orientations and 3D visualizations of estimated hand orientation after gravity-vector and magnetic North vector correction.

5 Conclusion

We believe that our approach to correct the drift in the gyroscope measurements and compensate the orientation estimation using the gravity vector and magnetic North vector will be useful in orientation tracking of human hand motion. Since the precise data of the hand or the fingers orientation can be obtained, this can lead to the improvement of 3D user interfaces to become more realistic. Thus, it could also enhance the experiences of the natural interaction of humans with a 3D virtual environment.