A review of methods for mapping 3D rotations to one-dimensional planar rotations: analysis, comparison, and a novel efficient Quaternion-based approach

Abstract

In numerous applications, extracting a single rotation component (termed “planar rotation”) from a 3D rotation is of significant interest. In biomechanics, for example, the analysis of joint angles within anatomical planes offers better clinical interpretability than spatial rotations. Moreover, in parallel kinematics robotic machines, unwished rotations about an axis – termed “parasitic motions” – need to be excluded. However, due to the non-Abelian nature of spatial rotations, these components cannot be extracted by simple projections as in a vector space. Despite extensive discussion in the literature about the non-uniqueness and distortion of the results due to the nonlinearity of the SO(3) group, they continue to be used due to the absence of alternatives. This paper reviews the existing methods for planar-rotation extraction from 3D rotations, showing their similarities and differences as well as inconsistencies by mathematical analysis as well as two application cases, one of them from biomechanics (flexural knee angle in the sagittal plane). Moreover, a novel, simple, and efficient method based on a pseudo-projection of the Quaternion rotation vector is introduced, which circumvents the ambiguity and distortion problems of existing approaches. In this respect, a novel method for determining the orientation of a box from camera recordings based on a two-plane projection is also proposed, which yields more precise results than the existing Perspective 3-Point Problem from the literature. This paper focuses exclusively on the case of finite rotations, as infinitesimal rotations within a single plane are non-holonomic and, through integration, produce rotation components orthogonal to the plane.

1. Introduction

The orientation of rigid bodies serves as a fundamental characteristic in describing the state of single or multi-body systems. The extraction of a single rotation component (termed “planar rotation”) from a 3D rotation plays a crucial role in many applications. In biomechanics, for example, the analysis of joint angles within anatomic body planes offers better clinical interpretability than spatial rotations. Also, for example, in lower-mobility parallel kinematic robotic chains, often an unwished rotation – termed “parasitic motion” – is targeted to be excluded, and its determination from the 3D rotation is needed. However, due to the non-Abelian structure of spatial rotations, it is known that such “planar” components cannot be determined by simple projections as in a vector space, and well-known inconsistencies arise. Despite an extensive literature and discussion about the non-uniqueness and distortion of the results, such extraction techniques continue to be used due to the absence of alternatives, lacking a consistent standardization.

In this paper, the existing extraction techniques from the literature are reviewed, and a thorough analysis of their differences and similarities is performed. Hereby, the existing approaches can be grouped into two basic procedures. In the first group, the spatial rotation is decomposed into a sequence of three elementary rotations about follow-up orthogonal coordinate axes (Euler, roll-pitch-yaw, or Bryant angles), and one of them is identified as the sought-extracted rotation, where the differences lie in the chosen sequence of rotations. In the second group, a particular coordinate axis of the rotated coordinate frame is projected onto the target plane of the initial coordinate frame, assuming that the desired extracted angle is the angle between the projected coordinate axis and a reference coordinate axis defining the zero-degree angle. The paper shows on one side that the two approaches are fully isomorphic with respect to each other depending on the rotation sequence and projected axis used. On the other side, the paper shows that using these extraction methods leads to inconsistencies and ambiguities, which makes standardization difficult.

In response to these problems, the paper presents an alternative novel efficient method for extracting a rotational component in a plane. This method decomposes a spatial Quaternion into a Quaternion about an axis normal to the extraction plane and a Quaternion about an axis within this plane. As will be shown, unlike traditional methods, the Quaternion-based approach yields unique extractions of 1D components from 3D rotations. The Quaternion decomposition method was preliminarily proposed in a previous conference paper [23] and is here expanded to its full extent.

The differences between the methods are illustrated in a virtual example and an experimental case for the extraction of knee flexural angle in the sagittal plane from gait analysis, comparing it to a Vicon measurement. In this respect, also a novel method for determining the orientation of a box based on a two-plane projection from camera recordings is proposed, which yields more precise results than the existing solutions for the Perspective 3-Point Problem (P3p) from the literature used by default in most applications (as e.g. implemented in Matlab), and which is reviewed here for better reference.

2. Review of existing planar rotation extraction methods

Single rotation components (termed in the following “planar rotations”) are, in many instances, such as navigation [48, 56, 57], human motion analysis [5], robotics [10, 14, 51, 58], and folding of packing cartons [38], of higher interest than the spatial orientation itself. In the analysis of human motion, for example, the extraction of planar rotations is essential for the analysis of anatomical segment angles. Angles projected onto anatomical body planes are favored as more relevant for intuitive clinical interpretation than their three-dimensional counterparts [54]. The sagittal plane, which divides the body into a right and a left section, is particularly relevant as an extraction plane, as human locomotion predominantly takes place in this plane [46].

The extraction of planar rotations from spatial orientations is also of crucial importance in the field of robotics. In parallel kinematics machines with less than three rotational degrees of freedom, intended autonomous movements can be accompanied by undesired, constrained rotations about axes that do not correspond to any degree of freedom, so-called parasitic rotations. Although the reduced number of degrees of freedom offers advantages such as cost reduction due to fewer actuators, simpler kinematics, and easier control [36], it can also lead to problems in various applications due to unpredictability and lower accuracy [10, 34, 37, 44, 47, 50]. Thus, several approaches are considered in the literature for eliminating parasitic rotations by initially identifying them from known spatial rotations [10, 34, 35].

2.1. Problem statement

Consider a frame ${\mathcal K}_{B}$ visualized by a box that undergoes a spatial rotation with respect to a space-fixed “world” reference frame ${\mathcal K}_{W}$ (Figure 1). The world frame can also be moving, in which case the rotation is then a relative motion without restriction of generality. The task is to extract a “planar” rotation angle $\alpha ^{z}$ that describes as close as possible the “component” of the general rotation in a plane, here chosen as the $x,y$ plane of the world frame, again without restrictions of generality.

Figure 1. Example used to discuss the extraction methods.

Due to the non-commutativity of spatial rotations, 3D rotations are not isomorphic to vector spaces, and hence the sought “planar” rotation (a) cannot be determined by classical vector-space projection and (b) is not unique. For this reason, several approaches exist in the literature for estimating this rotation angle, with diverging results.

For a unified comparison of existing methods, we adopt the well-known Quaternion description (also known as Euler-Rodrigues parameters) for the rotation matrix ${}^{}_{}{\textbf {R}}_{B}$ transforming vector components of the moving frame ${\mathcal K}_{B}$ with respect to the world frame ${\mathcal K}_{W}$ . By Euler’s theorem, the rotation from the world frame to the moving frame can be described by a single rotation about an axis $\underline{u}$ by an angle $\alpha$ (Figure 1), which is embodied by the 4-parametric Quaternion $Q$ as [3], see also Appendix A:

(1) \begin{align} Q = \begin{bmatrix} \;q_{0}\; \\[0.5em] \underline{q} \end{bmatrix} \quad \textrm{ with }\quad q_0 = \cos \frac{\displaystyle \alpha }{\displaystyle 2} \quad ; \quad \underline{q} = \begin{bmatrix} \;q_{x}\; \\[0.5em] q_{y} \\[0.5em] q_{z} \end{bmatrix} = \sin \frac{\displaystyle \alpha }{\displaystyle 2} \; \underline{u} \ \ . \end{align}

The elements of the rotation matrix are termed as:

(2) \begin{align}{}^{}_{}{\mathsf{R}}_{B} = \begin{bmatrix} \;\;\;\rho _{11} \;\;\;&\;\;\;\rho _{12} \;\;\;&\;\;\;\rho _{13} \;\;\; \\[0.5em] \;\;\;\rho _{21} \;\;\;&\;\;\;\rho _{22} \;\;\;&\;\;\;\rho _{23} \;\;\; \\[0.5em] \;\;\;\rho _{31} \;\;\;&\;\;\;\rho _{32} \;\;\;&\;\;\;\rho _{33} \;\;\; \\[0.5em] \end{bmatrix} \ \, \end{align}

where each column of the matrix is the projection of the corresponding axis of $\,{\mathcal K}_{B}$ onto $\,{\mathcal K}_{W}$ . The rotation matrix ${}^{}_{}\boldsymbol{\mathsf{R}}_{B}$ can then be expressed by the four components of the Quaternion $Q$ as [41], Appendix A:

(3)

[Note that the rotation matrix can be easily determined from a given Quaternion by Eq. (3), or, vice-versa, the Quaternion can be easily determined by inverse relationships for a given rotation matrix (Appendix A). Thus, both representations are fully equivalent, and the most suited one can be employed for a given application. For example, when integrating angular velocities, the Quaternion method is more immediate and natural (Appendix A), while, for given rotation matrices, the inverse relationship is assumed to be pre-computed for the rest of this paper.]

It is clear that if the spatial rotation takes place only about an axis orthogonal to the plane (here the $z$ axis), then the rotation axis $\underline{u}$ is also orthogonal to the plane and the sought rotation angle $\alpha ^{z}$ is identical to the full Quaternion rotation angle $\alpha$ . On the other hand, if the rotation vector $\underline{u}$ lies within the $x,y$ plane, then no rotation about the orthogonal axis $z$ takes place. The question is, what are the planar rotation angles for rotation axes having components both in the $x,y$ plane ( $\underline{q}^{xy}$ ) and the orthogonal direction $z$ ( $\underline{q}^{z}$ )? This is the basis of the differences between current planar rotation extraction methods.

Hereby, two groups of approaches can be identified: (1) decomposing the spatial rotation in a sequence of elementary rotations about coordinate axes and selecting one of them as the sought planar rotation angle $\alpha ^{z}$ (here termed Euler-angle decomposition methods), and (2) projecting an axis of the rotated frame on the plane and determining the angle between the projected axis and a reference axis in the plane as the sought planar rotation angle $\alpha ^{z}$ (here termed coordinate-axis projection methods). As will be shown, both methods are identical. In the sequel, both methods are reviewed and discussed.

2.2. Euler-angle decomposition methods

The Euler-angles approach views a spatial rotation as a sequence of three elementary rotations about the axes of a step-wise rotating orthogonal coordinate system under the condition that no two consecutive rotations share the same absolute axis. The extraction of planar rotations from spatial orientations via Euler angles is based on the assumption that a partial Euler rotation angle reflects the desired “planar” rotation angle [10, 12, 29, 34, 35]. The choice against computationally more robust approaches such as helical angles [53] or Quaternions [27] stems from the assumption that Euler angles are more intuitive for the extraction process. However, this approach is not unproblematic, since the sequence of Euler angle elementary rotations can significantly influence the value of the extracted “planar” angle [11, 15, 23], and the resulting angles can deviate considerably from the expectations [8, 21]. Despite these well-known problems, Euler angles are still the most common approach in the literature [5]. Efforts to circumvent these problems have led to attempts to establish predefined Euler sequences as a standard for specific applications, such as those proposed by the International Society of Biomechanics for human motion analysis [13, 55]. While this strategy improves the comparability of study results, it does not address the inherent lack of uniqueness of the extracted angles.

To ensure the uniqueness of angle projections, sequences where the first and last rotation occurs around the same local axis (“Proper Euler angles”) are not applied for this approach. This includes sequences such as $z$ - $y$ - $z$ . Instead, only rotation sequences that ensure that each rotation occurs around a different axis (“Tait-Bryan angles”) are used. This includes sequences like $z$ - $y$ - $x$ etc., leading to six potential rotation sequences for planar angle extractions. Often, the sequence is chosen such that the first rotation is about the axis that is normal to the desired projection plane, followed by rotations around the remaining axes [10, 11, 15, 23, 34, 35].

Each sequence can be interpreted in two distinct ways: extrinsically, with rotations about inertially-fixed axes ( ${\mathcal K}_{W}$ in the given example), and intrinsically, with rotations about step-wise moving body-fixed axes ( ${\mathcal K}_{B}$ in the given example). Both intrinsic [11, 15, 23] and extrinsic rotations [10, 34, 35] are utilized in the literature for planar angle extractions. As intrinsic Euler angles are mathematically equivalent to the reversed sequence of extrinsic rotations, the discussion will focus on intrinsic rotations, implicitly covering also extrinsically interpreted sequences.

The elementary rotations about follow-up $\;x$ -, $y$ -, and $z$ -axes, respectively, and defined by the unknown Euler angles $\varphi$ , $\Theta$ and the sought angle $\alpha ^{z}$ , respectively, are:

(4) \begin{align} \textrm{Rot}[ x, \varphi ] & = \left[\begin{array}{ccc} 1 & 0 & 0 \\[3pt] 0 & \cos \varphi & -\sin \varphi \\[3pt] 0 & \sin \varphi &{+}\cos \varphi \end{array}\right]\!;\ \textrm{Rot}[ y, \Theta ] = \left[\begin{array}{ccc} {+}\cos \Theta & 0 & \sin \Theta \\[3pt] 0 & 1 & 0 \\[3pt] -\sin \Theta & 0 & \cos \Theta \end{array}\right]\!;\ \textrm{Rot}[ z, \alpha ^{z} ] & = \left[\begin{array}{ccc} \cos \alpha ^{z} & -\sin \alpha ^{z} & 0 \\[3pt] \sin \alpha ^{z} &{+}\cos \alpha ^{z} & 0 \\[3pt] 0 & 0 & 1 \end{array}\right] \end{align}

Due to the non-Abelian structure inherent in spatial rotations, each Euler sequence, represented by a non-commutative multiplication of the elementary rotations, leads to distinct definitions of the rotation matrix ${}^{}_{}\boldsymbol{\mathsf{R}}_{B}$ . For the six possible Tait-Bryan angles employed in the literature for extraction, the rotation matrices shown in Eq. (5) to Eq. (10) are obtained. To maintain clarity and facilitate further discussion, only the relevant elements of the resulting matrices for determining the sought angle $\alpha ^{z}$ are shown. Furthermore, the additional index on ${}^{}_{}\boldsymbol{\mathsf{R}}_{B}$ and $\alpha ^{z}$ serves to indicate the chosen Euler sequence and assists in distinguishing the extracted angles based on the extraction method.

(5) \begin{align} \boldsymbol{\mathsf{R}}^{\textrm{zyx}}_B & = \textrm{Rot}[ z, \alpha ^{z}_{\textrm{zyx}} ]\cdot \textrm{Rot}[ y, \Theta ]\cdot \textrm{Rot}[ x, \varphi ] = \left[ \begin{array}{@ c@ c@ c@{\hspace*{15pt}}} {+} \cos \alpha ^{z}_{\textrm{zyx}} \cos \Theta & * & * \\[3pt] {+}\sin \alpha ^{z}_{\textrm{zyx}}\cos \Theta & * & * \\[3pt] * & * & * \end{array}\right] \end{align}

(6) \begin{align} \boldsymbol{\mathsf{R}}^{\textrm{zxy}}_B & = \textrm{Rot}[ z, \alpha ^{z}_{\textrm{zxy}} ]\cdot \textrm{Rot}[ x, \varphi ]\cdot \textrm{Rot}[ y, \Theta ] = \left[ \begin{array}{@c@c@{}c@{}} * & -\sin \alpha ^{z}_{\textrm{zxy}}\cos \varphi &*\\[3pt] *& {+} \cos \alpha ^{z}_{\textrm{zxy}}\cos \varphi &*\\[3pt] * & * & * \end{array}\right] \end{align}

(7) \begin{align} \boldsymbol{\mathsf{R}}^{\textrm{xyz}}_B & = \textrm{Rot}[ x, \varphi ]\cdot \textrm{Rot}[ y, \Theta ]\cdot \textrm{Rot}[ z, \alpha ^{z}_{\textrm{xyz}} ] = \left[\begin{array}{@ c@{\quad}c@{}c@{}} {+}\cos \alpha ^{z}_{\textrm{xyz}}\cos \Theta & -\sin \alpha ^{z}_{\textrm{xyz}}\cos \Theta & * \\[3pt] * & * & *\\[3pt] * & * & * \end{array}\right] \end{align}

(8) \begin{align} \boldsymbol{\mathsf{R}}^{\textrm{yxz}}_B & = \textrm{Rot}[ y, \Theta ]\cdot \textrm{Rot}[ x, \varphi ]\cdot \textrm{Rot}[ z, \alpha ^{z}_{\textrm{yxz}} ] = \left[ \begin{array}{@ c@{\quad}c@{}c@{}} * & * & * \\[3pt] {+}\sin \alpha ^{z}_{\textrm{yxz}}\cos \varphi & {+}\cos \alpha ^{z}_{\textrm{yxz}}\cos \varphi &*\\[3pt] * & * & * \end{array}\right] \end{align}

(9) \begin{align} \boldsymbol{\mathsf{R}}^{\textrm{xzy}}_B & = \textrm{Rot}[ x, \varphi ]\cdot \textrm{Rot}[ z, \alpha ^{z}_{\textrm{xzy}} ]\cdot \textrm{Rot}[ y, \Theta ] = \left[ \begin{array}{@{}c@{}c@ c@{}} *&-\sin \alpha ^{z}_{\textrm{xzy}}&*\\[3pt] * & * & *\\[3pt] * & * & * \end{array}\right] \end{align}

(10) \begin{align}\boldsymbol{\mathsf{R}}^{\textrm{yzx}}_B & = \textrm{Rot}[ y, \Theta ]\,\cdot \textrm{Rot}[ z, \alpha ^{z}_{\textrm{yzx}} ]\cdot \textrm{Rot}[ x, \varphi ] = \left[ \begin{array}{@{}c@{}c@{}c@{}} * & * & *\\[3pt] {+}\sin \alpha ^{z}_{\textrm{yzx}}&*&*\\[3pt] * & * & * \end{array}\right] \end{align}

The sought “planar” rotation angle $\alpha ^{z}$ is then obtained by coefficient comparison between the rotation matrix Eq. (2) and Eq. (5) to Eq. (10), respectively, as:

(11) \begin{align} \alpha ^{\textrm{z}}_{\textrm{zyx}} &= \textrm{atan2} \big(\pm \rho _{21},\ \ \pm \rho _{11} \big) & &\text{[5, 9, 54, 55]} \end{align}

(12) \begin{align} \alpha ^{\textrm{z}}_{\textrm{zxy}} & = \textrm{atan2} \big ({-} \pm \,\rho _{12},\ \ \pm \rho _{22} \big ) & &\text{[36, 10, 11, 13, 12, 15, 29, 46, 54]} \end{align}

(13) \begin{align} \alpha ^{\textrm{z}}_{\textrm{xyz}} & = \textrm{atan2} \big ({-} \pm \,\rho _{12},\ \ \pm \rho _{11} \big ) && \text{[5, 8, 16, 54]} \end{align}

(14) \begin{align} \alpha ^{\textrm{z}}_{\textrm{yxz}} & = \textrm{atan2} \big (\pm \rho _{21},\ \ \pm \rho _{22} \big ) && \text{[8]} \end{align}

(15) \begin{align} \alpha ^{\textrm{z}}_{\textrm{xzy}} & = \textrm{asin} \big ({-}\rho _{12} \big) &&\text{[49]} \end{align}

(16) \begin{align} \alpha ^{\textrm{z}}_{\textrm{yzx}} & = \textrm{asin} \big (\rho _{21}\big ) &&\text{[8]}, \end{align}

where the $\pm$ sign corresponds to the sign of the cosine of the second elementary rotation, and the references next to the expressions show the applications from the literature for each extraction method. It is essential to be aware of potential singularities that may arise when the second rotation is close to $90^{\circ }$ , in which case the “planar” rotation angle $\alpha ^{z}$ cannot be determined.

2.3. Coordinate-axis projection methods

The second group of extraction methods assumes that the desired extracted “planar” angle can be viewed as the angle between a projected coordinate axis of the moving ${\mathcal K}_{B}$ onto the target plane and a reference coordinate direction in the plane for the zero angle definition [17, 33] (Figure 2). As will be shown here (to our knowledge for the first time), the procedure yields exactly the same results as the Euler-angle decomposition for a corresponding sequence.

Out of the three axes of the moving frame ${\mathcal K}_{B}$ , the axis normal to the projection plane in the non-rotated state (here the $z$ -axis) cannot be used, as it would be “invisible” when no rotation takes place or when the rotation is already a “planar” rotation about the plane normal. Therefore, any choice of the other two axes (here, the $x$ - and $y$ -axes) can be made. As the columns of ${}^{}_{}\boldsymbol{\mathsf{R}}_{B}$ are already the projection of axes of ${\mathcal K}_{B}$ onto ${\mathcal K}_{W}$ , the projected (moving) $x$ and $y$ axes, denoted as $\overline{\underline{x}}$ and $\overline{\underline{y}}$ , are in terms of ${}^{}_{}\boldsymbol{\mathsf{R}}_{B}$

(17) \begin{align} \overline{\underline{x}} = \begin{bmatrix} \;\overline{x}_{x}\;\\[0.5em] \overline{x}_{y} \end{bmatrix} = \begin{bmatrix} \;\rho _{11}\;\\[0.5em] \;\rho _{21}\; \end{bmatrix} \quad \textrm{ and }\quad \overline{\underline{y}} = \begin{bmatrix} \;\overline{y}_{x}\;\\[0.5em] \overline{y}_{y} \end{bmatrix} = \begin{bmatrix} \;\rho _{12}\;\\[0.5em] \;\rho _{22}\; \end{bmatrix} \ \ . \end{align}

The axes defining zero degrees are in the projection plane

(18) \begin{align} \overline{\underline{e}}_{x} = \begin{bmatrix} \;1 \;\\[0.5em] 0 \end{bmatrix} \quad \textrm{ and }\quad \overline{\underline{e}}_{y} = \begin{bmatrix} 0 \\[0.5em] \;1\; \end{bmatrix} \ \ . \end{align}

Figure 2. Visualization of the axis projection method.

The resulting angle extractions (termed as $\alpha ^{\textrm{z}}_{\textrm{Ax}}$ and $\alpha ^{\textrm{z}}_{\textrm{Ay}}$ , respectively) then result as:

(19) \begin{align} \alpha ^{\textrm{z}}_{\textrm{Ax}} ={+}\textrm{sgn}\left (\;\overline{x}_{y}\;\right )\,\textrm{acos}\bigg (\;\frac{\overline{\underline{e}}_{x}^{\textrm{T}}\;\overline{\underline{x}}}{\|\,\overline{\underline{e}}_{x}\,\|\cdot \|\,\overline{\underline{x}}\,\|}\;\bigg )=\textrm{sgn}\left (\;\rho _{21}\;\right )\,\textrm{acos}\left (\;\frac{\rho _{11}}{\sqrt{\rho _{11}^2+\rho _{21}^2}}\;\right )\;\;\text{[17]} \end{align}

(20) \begin{align} \alpha ^{\textrm{z}}_{\textrm{Ay}} = -\textrm{sgn}\left (\;\overline{y}_{x}\;\right )\, \textrm{acos}\bigg (\;\frac{\overline{\underline{e}}_{y}^{\textrm{T}}\;\overline{\underline{y}}}{\|\,\overline{\underline{e}}_{y}\,\|\cdot \|\,\overline{\underline{y}}\,\|}\;\bigg )= -\textrm{sgn}\left (\;\rho _{12}\;\right )\,\textrm{acos}\left (\;\frac{\rho _{22}}{\sqrt{\rho _{12}^2+\rho _{22}^2}}\;\right )\;\;\text{[33]} \end{align}

where the references next to the expressions show the applications from the literature for each extraction method.

Using the relationship

(21) \begin{align} \textrm{acos}\bigg (\;\frac{A}{\sqrt{A^2+B^2}}\;\bigg ) = \textrm{atan2}\big (\;\;\,B\;\;\,,\;\;\,A\;\;\,) \end{align}

one obtains

(22) \begin{eqnarray} \alpha ^{\textrm{z}}_{\textrm{Ax}} &=& \textrm{atan2}\big (\pm \rho _{21},\ \pm \rho _{11}\big ) \end{eqnarray}

(23) \begin{eqnarray} \alpha ^{\textrm{z}}_{\textrm{Ay}} &=& \textrm{atan2}\big (\mp \rho _{12},\ \pm \rho _{22} \big ) \end{eqnarray}

which are identical to the Euler $z$ - $y$ - $x$ and $z$ - $x$ - $y$ extraction angles, respectively, including singularities:

(24) \begin{align} \alpha ^{\textrm{z}}_{\textrm{Ax}} \equiv \alpha ^{z}_{\textrm{zyx}} \quad \textrm{ and }\quad \alpha ^{\textrm{z}}_{\textrm{Ay}} \equiv \alpha ^{z}_{\textrm{zxy}} \end{align}

3. Alternative novel Quaternion projection method

As an alternative to Euler-angle extraction, consider the decomposition of the total Quaternion (Figure 1)

(25) \begin{align} Q \;=\; Q^{\textrm{z}} \circ Q^{\textrm{xy}} \end{align}

where “ $\circ$ ” denotes Quaternion composition (Appendix A) and where the Quaternion “components”

(26) \begin{eqnarray} \,Q^{\textrm{z}} \;&=&\; \begin{bmatrix} \;\,q_{0}^{\textrm{z}}\;\,\\[0.5em] \underline{q}^{\textrm{z}} \end{bmatrix} \quad \textrm{ with }\quad \underline{q}^{\textrm{z}} \;=\; \begin{bmatrix} 0 \\[0.5em] 0 \\[0.5em] \;\,q_{z}^{\textrm{z}}\;\, \end{bmatrix} \end{eqnarray}

(27) \begin{eqnarray} Q^{\textrm{xy}} \;&=&\; \begin{bmatrix} \;q_{0}^{\textrm{xy}}\;\\[0.5em] \underline{q}^{\textrm{xy}} \end{bmatrix} \quad \textrm{ with }\quad \underline{q}^{\textrm{xy}} \;=\; \begin{bmatrix} \;q_{x}^{\textrm{xy}} \;\\[0.5em] q_{y}^{\textrm{xy}} \\[0.5em] 0 \end{bmatrix} \end{eqnarray}

are still unknown. The sought quantity is the extracted planar rotation $\alpha ^{\textrm{z}}_{\textrm{Q}}$ of the first component

(28) \begin{align} q_{0}^{\textrm{z}}= \cos{\left (\; \frac{\alpha ^{\textrm{z}}_{\textrm{Q}}}{2} \;\right ) } \quad \textrm{ and }\quad q_{z}^{\textrm{z}} = \sin{\left (\; \frac{\alpha ^{\textrm{z}}_{\textrm{Q}}}{2} \;\right ) } \end{align}

while the second component is an unknown rotation about an axis within the plane. Using the Quaternion inverse

(29) \begin{align} \left \{ Q^{\textrm{xy}} \right \}^{-1} =\; \overline{Q}^{\textrm{\,xy}} = \,\begin{bmatrix} q_{0}^{\textrm{xy}} \\[0.5em] \;-\underline{q}^{\textrm{xy}} \; \end{bmatrix} \end{align}

and the Quaternion composition rule leads to

(30) \begin{align} Q^{\textrm{z}} \;=\; Q \circ \left \{ Q^{\textrm{xy}} \right \}^{-1} =\; \begin{bmatrix} q_{0} q_{0}^{\textrm{xy}} + \underline{q} \cdot \underline{q}^{\textrm{xy}} \\[0.5em] \; - q_{0} \underline{q}^{\textrm{xy}} + q_{0}^{\textrm{xy}} \underline{q} - \underline{q} \times \underline{q}^{\textrm{xy}} \; \end{bmatrix} \end{align}

from where the following four scalar equations are obtained:

(31) \begin{eqnarray} \; q_{0}^{\textrm{z}} \; &=& \; q_{0}q_{0}^{\textrm{xy}} \; + q_{x}q_{x}^{\textrm{xy}} \;+ q_{y}q_{y}^{\textrm{xy}} \end{eqnarray}

(32) \begin{eqnarray} \;\;\, 0 \; &=& \; q_{x}q_{0}^{\textrm{xy}} \;- q_{0}q_{x}^{\textrm{xy}} \;\, + q_{z}q_{y}^{\textrm{xy}} \end{eqnarray}

(33) \begin{eqnarray} \;\; 0 \; &=& \; q_{y}q_{0}^{\textrm{xy}} \;- q_{z}q_{x}^{\textrm{xy}} \;- q_{0}q_{y}^{\textrm{xy}} \end{eqnarray}

(34) \begin{eqnarray} q_{z}^{\textrm{z}} \; &=& \; q_{z}q_{0}^{\textrm{xy}} \;\, + q_{y}q_{x}^{\textrm{xy}} - q_{x}q_{y}^{\textrm{xy}} \end{eqnarray}

Determining $q_{x}^{\textrm{xy}}$ and $q_{y}^{\textrm{xy}}$ from Eqs. (32) and (33) and inserting in Eq. (31) and Eq. (34) gives

(35) \begin{align} q_{x}^{\textrm{xy}} \;\; = \;\; \frac{\qquad \; q_{0}q_{x} + q_{y}q_{z}\qquad \; }{q_{0}^{2} + q_{z}^{2}}\,q_{0}^{\textrm{xy}} \qquad \qquad \qquad \quad \;\;\; \qquad \quad \;\;\; \left [\;\;=\quad \pm \frac{\;q_{0}q_{x} + q_{y}q_{z}\;}{\sqrt{\;q_{0}^{2} + q_{z}^{2}\;}}\;\;\right ] \end{align}

(36) \begin{align} q_{y}^{\textrm{xy}} \;\; = \;\; \frac{\qquad \; q_{0}q_{y} - q_{x}q_{z}\qquad \;}{q_{0}^{2} + q_{z}^{2}}\,q_{0}^{\textrm{xy}} \qquad \qquad \qquad \quad \;\;\; \qquad \quad \;\;\; \left [\;\;=\quad \pm \frac{\;q_{0}q_{y} - q_{x}q_{z}\;}{\sqrt{\;q_{0}^{2} + q_{z}^{2}\;}}\;\;\right ] \end{align}

(37) \begin{align} q_{0}^{\textrm{z}}\; \;\; = \;\; \frac{q_{0} \left ( \;q_{0}^{2} + q_{x}^{2} + q_{y}^{2} + q_{z}^{2} \;\right ) }{q_{0}^{2} + q_{z}^{2}}\,q_{0}^{\textrm{xy}} \;\; = \;\; \frac{q_{0}}{\;q_{0}^{2} + q_{z}^{2}\;}\,q_{0}^{\textrm{xy}} \qquad \qquad \; \left [\;\;=\quad \pm \frac{\qquad q_{0}\qquad }{\sqrt{\;q_{0}^{2} + q_{z}^{2}\;}}\, \,\right ] \end{align}

(38) \begin{align} \;\, q_{z}^{\textrm{z}}\; \;\; = \;\; \frac{\;q_{z} \left ( q_{0}^{2} + q_{x}^{2} + q_{y}^{2} + q_{z}^{2} \;\right ) }{q_{0}^{2} + q_{z}^{2}}\,q_{0}^{\textrm{xy}} \;\; = \;\; \frac{q_{z}}{\; q_{0}^{2} + q_{z}^{2}\;}\,q_{0}^{\textrm{xy}} \qquad \qquad \;\! \left [\;\;=\quad \pm \frac{\qquad \, q_{z}\qquad }{\sqrt{\;q_{0}^{2} + q_{z}^{2}\;}}\;\right ] \end{align}

The expression in brackets follows from the normalization achieved by squaring and adding Eq. (37) and Eq. (38) and considering the normalizing condition $(q_{0}^{\textrm{z}})^{2} + (q_{z}^{\textrm{z}})^{2} = 1$ yielding

(39) \begin{align} q_{0}^{\textrm{xy}} \;=\; \pm \,\sqrt{\;q_{0}^{2} + q_{z}^{2}\;} \end{align}

Thus, one obtains by dividing the expressions in the brackets of Eq. (37) and Eq. (38)

(40) \begin{align} \frac{\;q_{z}^{\textrm{z}}\;}{\;q_{0}^{\textrm{z}}\;} \;=\; \frac{\;\pm q_z\;}{\;\pm q_{0}\;} \end{align}

in which all variables of the second Quaternion component are eliminated, and with Eq. (28) it holds

(41) \begin{align} \alpha ^{\textrm{z}}_{\textrm{Q}} = 2\cdot \textrm{atan2}\big (\pm q_{z},\ \pm q_{0} \big ) \ \ . \end{align}

It is remarkable that the planar rotation about the planar $z$ axis is just the $z$ component of the original Quaternion vector, scaled by the reciprocal of the real part of the original Quaternion. Thus, the extraction results as a simple scaled pseudo-projection of the original Quaternion vector to the axis orthogonal to the plane. A further remarkable property from Eq. (41) is that, although there are formally two solutions, these are separated by 360 $^{\circ }$ , thus actually the plus sign is unique for a full rotation. Even more remarkably, switching the order of components in Eq. (25), i.e., placing the $z$ rotation as second component of the decomposition, only leads to a change of sign in the cross product in Eq. (30), which again leads only to a change of the sign in the numerator of Eq. (35) and Eq. (36), respectively. Thus, the outcome for $\alpha ^{\textrm{z}}_{\textrm{Q}}$ remains the same, and one obtains the result that the decomposition Eq. (25) is commutative with respect to the planar rotation component. This makes the Quaternion extraction method a quasi-vector space with respect to the extracted angle $\alpha ^{z}$ . Finally, it can be observed that the solution of Eq. (37) and Eq. (38) becomes singular when both components $q_{0}$ and $q_{z}$ vanish at the same time. This corresponds to a rotation of 180 $^{\circ }$ about an axis in the projection plane. In this case, the space is rotated face-down, and if the rotation axis is skewed with respect to the $x$ and $y$ directions, then a spurious rotation about the $z$ appears to happen (Figure 3), although there is no such one. Such a case can however be excluded from the planar rotation extraction case, as it does not make sense to try to identify $z$ rotations when such a large rotation takes place only normal to the sought rotation axis.

Figure 3. Rotation of a rectangle ABCD by 180 $^{\circ }$ about its diagonal in the projection plane.

4. Comparison and discussion

The practical implications of the different planar angle extraction techniques of the literature, as well as the new Quaternion-based method, are illustrated in the following by two case examples: (1) a virtual rotation of a box generated by an arbitrary time development of the spatial angular velocity vector in body-fixed coordinates, and (2) an application from biomechanics, in which the sagittal knee angle is determined from experimental inertial measurement unit (IMU) measurements and compared to the gold-standard Vicon results from simultaneously infrared tracked retro-reflective markers.

4.1. Example 1: Virtual 3D rotation of a box

Let the kinematic input of a box be provided by the body-fixed angular velocity ${}^{B}_{}\underline{\omega }_{B}$

(42) \begin{align}{}^{B}_{}\underline{\omega }_{B} = \begin{bmatrix} \sin \big (\pi \tfrac{\tau }{\pi }\big ) \\[0.5em] \tfrac{3}{4}\sin \big (\tfrac{\pi }{6}\tau -\tfrac{\pi }{4}\big ) \\[0.5em] \;\tfrac{3}{2}\sin \big (\tfrac{2\pi }{3}\tau -\pi \big )\; \end{bmatrix} \end{align}

for the dimensionless time $0 \leq \tau \leq 6$ . The resulting time development $Q\left ( \tau \right )$ of the resulting Quaternion can be easily and precisely determined by Runge-Kutta integration of fourth order (Appendix A) with a step size of 2 ms. The progression of the box orientation over dimensionless time is shown in Figure 4, where the poses are “stretched” in translational direction for visual reasons to avoid overlaps. The upper panel displays the Quaternion vector with direction equal to the rotation vector $\underline{u}$ and magnitude equal to sine of half rotation angle, where Quaternion vectors close to the projection plane are colored gray. The second panel displays the total spatial rotation of the box. The other three panels display the rotation of the box when applying only the planar rotation resulting from the Quaternion extraction and the extraction by Euler sequences $z$ - $x$ - $y$ and $z$ - $y$ - $x$ , respectively. The resulting progression of the planar rotation angles $\alpha ^{z}$ is shown in Figure 5 for all extraction methods.

Figure 4. “Stretched” rotation points over time with the extractions $\alpha ^{z}_Q$ , $\alpha ^{z}_{zxy}$ , and $\alpha ^{z}_{zyx}$ .

Figure 5. Extracted planar angle $\alpha ^{z}$ for all extraction methods.

One can see that at the beginning, all planar extraction methods behave similarly, but deviate substantially as rotation progresses. A closer analysis reveals the causes for similarity and divergence: At the beginning, rotations are small and thus close to infinitesimal relationships, whereas rotations are commutative and thus lie in a vector space. Thus, the extraction process corresponds approximately to a proper vector projection, and all extraction methods yield close to equal results. For the other regions, results of the Euler extraction methods diverge substantially, and even between two sequences with equal locations of the extracted planar rotation in the Euler sequence (e.g. $x$ - $y$ - $z$ vs. $y$ - $x$ - $z$ ). This shows that the supposedly intuitive Euler-angle extraction leads to the unintuitive result that the extracted angle $\alpha ^{z}$ changes significantly depending on the sequence of other two rotations, here $x$ and $y$ . Moreover, the Euler extraction methods produce inconsistent results at poses A, B, and C, where the rotation takes place about an axis within the projection plane and thus no $z$ -rotation would be expected (Figure 6). Thus, the Euler extraction methods are intuitive only for small rotation angles but become erratic and inconsistent for large rotations depending on the rotation sequence, in particular also on the choice of the sequence of the other (uninvolved) two. This is not the case for the Quaternion extraction method, which yields consistent and unique results (zero planar rotation angle when rotation takes place about an axis in the projection plane) and a smoother rotation trajectory.

Figure 6. Rotated box with its rotation axis entirely within the projection plane.

4.2. Example 2: Extraction of sagittal knee angle from 3D IMU measurements

As a second example, the measurement of the knee sagittal angle for human walking along a track using IMUs is regarded. The results are validated with respect to the current gold-standard technique for capturing human motion within the lab based on retro-reflective markers attached to the subject’s skin at specific anatomic reference points and subsequently recorded by a synchronized set of cameras within the room [2, 7, 40]. As known, marker-based measurements excel with precision and robustness, but have certain limitations such as the necessity for specialized laboratories, restricted mobility, high costs, and the need for specially trained personnel [4, 52]. On the other hand, IMUs are at advantage regarding low cost, suitability for outdoor measurements, and easiness and speed to attach. However, they also feature a couple of shortcomings, such as integration drift due to sensor errors, undefined orientation with respect to limbs, and unknown (varying) tilt offsets between IMU and anatomical axes.

The experimental setup is shown in Figure 7. Subjects were outfitted with reflective markers by a standardized plug-in gait protocol [18]. A movement measurement system (Vicon MX, Vicon Motion Systems Ltd UK, 100 Hz) comprising eight cameras, supplemented by a synchronized video camera (Basler 602 fc, 50 Hz), was utilized. All raw kinematic data of the Vicon reference system underwent filtering using a Woltring filter. A video camera positioned laterally to the walking track with its image plane parallel to the sagittal plane served to determine the constant offset between the relative sensor angle and the actual knee angle. In addition, two IMUs (Myon Aktos-T, myon AG, Schwarzenberg, Switzerland, 2000Hz) were placed in a box and affixed using double-sided adhesive tape to the right leg at the shank adjacent to the musculus tibialis anterior, and laterally to the thigh. Two further calibrated DSLR cameras (EOS 2000D, Canon) with known absolute orientations were positioned to determine the IMU box orientation at the start. In order to validate the proposed method against a gold-standard reference, the knee angles of ten test subjects were assessed during five trials of walking along a designated track, leading to 60 trials in total, employing both retro-reflective markers and the method delineated in this paper. An analog signal was employed to concurrently trigger both DSLR cameras, IMU sensor recordings, the Vicon system, and the laterally placed video camera, ensuring synchronization across all systems and marking the beginning of the individual trial.

Figure 7. Illustration of the experimental process for knee sagittal angle measurement featuring the same subject twice in the same scene: left scene: 20 sec standstill phase for initial target orientation; right scene: dynamic assessment during walking motion.

Five steps were performed to optimize the precision of IMU-based tracking of the sagittal joint angle:

1. 1) determination of precise initial IMU orientation using two CCD cameras (Appendix B)

2. 2) identification of constant error offset in angular velocity (see below)

3. 3) integration of IMU-measured angular velocity using fourth-order Runge-Kutta (Appendix A)

4. 4) extraction of planar angles of shank and thigh rotation according to methods in Sections 2 and 3

5. 5) offset shift of relative sagittal knee angle to anatomic angle (see below)

Step 1: Initial IMU box orientation

A new method was developed to determine the precise orientation of the IMU box with respect to the track world frame (Appendix B). The method uses a Newton iteration to determine the pose, which matches the projections of a checkerboard face on two camera planes with known orientation (Figure 8). Recordings from the triggered CCD cameras during a standstill phase for 20 s were utilized to determine the initial orientation of the sensors. The recorded average accuracy was $0.37^{\circ } \pm 0.15^{\circ }$ compared to classical methods such as in MATLAB with an accuracy of $2.03^{\circ } \pm 0.76^{\circ }$ .

Figure 8. Left: exploded view of the target with an inertial measurement unit including a gyroscope. Middle: targets positioned on the observed leg, highlighting the relative sensor angle ( $\alpha$ ) and knee angle ( $\beta$ ). Right: lateral video camera frame used for a one-time knee angle estimation.

Step 2: Compensation of constant sensor error offset

The main cause of sensor error was found to be a constant offset for the measured angular velocity. This constant error offset $\underline{\delta }$ was determined for each sensor as the average value during the standstill phase of 20 s. Then, the value was subtracted for the actual sensor measurement ${}^{T}_{}\underline{\omega }_{}^{\ast }(t)$ to obtain the corrected body-fixed angular velocity ${}^{T}_{}\underline{\omega }_{}(t)$ relative to the sensor coordinate frame ${\mathcal K}_{T}$ as:

(43) \begin{align}{}^{T}_{}\underline{\omega }_{}(t) ={}^{T}_{}\underline{\omega }_{}^{\ast }(t) - \underline{\delta } \end{align}

Step 3: Integration of IMU-measured angular velocity

Using the initial IMU orientation obtained from step 1, a fourth-order Runge-Kutta integration method for determining the resulting Quaternions for the IMU sensors of shank $Q_{S}(t)$ and thigh $Q_{T}(t)$ including normalization of the resulting Quaternions after each integration step was applied (Appendix A).

Step 4: Extraction of planar rotation angles

From the resulting Quaternions of shank $Q_{S}(t)$ and thigh $Q_{T}(t)$ , respectively, the corresponding sagittal angles $\alpha ^{\textrm{z}}_S(t)$ and $\alpha ^{\textrm{z}}_T(t)$ of shank and thigh, respectively, were determined for all methods described in the previous sections, specifically Eq. (11) to Eq. (16) for Euler-angle decompositions and Eq. (41) for Quaternion decomposition.

Step 5: Offset shift of relative sagittal knee angle to anatomic angle

The IMU relative planar angle $\alpha ^{\textrm{z}}(t)$ at the knee was computed as the difference between the planar angles of thigh and shank

(44) \begin{align} \alpha ^{\textrm{z}}(t) = \alpha ^{\textrm{z}}_S(t) - \alpha ^{\textrm{z}}_T(t) \end{align}

Between the IMU planar relative angle $\alpha ^{\textrm{z}}(t)$ and the “true” anatomical sagittal knee angle $\beta ^{\textrm{z}}(t)$ , a constant offset $\beta ^{\textrm{z}}_0$ (Figure 8) was determined such that:

(45) \begin{align} \beta ^{\textrm{z}}(t) = \alpha ^{\textrm{z}}(t) - \beta ^{\textrm{z}}_0 \end{align}

The constant offset $\beta ^{\textrm{z}}_0$ was estimated by an experienced orthopedic professional visually assessing the sagittal knee angle at an arbitrarily chosen time point within the gait cycle at the laterally positioned camera and comparing this to the Vicon value at this point, which reached an average accuracy of $2.1^{\circ }$ compared to the Vicon measurement.

Resulting sagittal knee values for all planar rotation extraction methods and comparison to Vicon

The sagittal knee angle for the selected gait cycle was measured using Vicon and the methods presented in this section. Vicon measurements and sensor data were time-normalized to 100% for the described gait cycle. The IMU-sampling rate and the Runge-Kutta step size were 1 ms and 2 ms, respectively. The root-mean-square error (RMSE) was calculated for each trial by comparing the calculated sagittal knee angle with the result from Vicon. Figure 9 displays the resulting average sagittal knee angle over a trial for one person, while the RMSEs for all sagittal angle extraction methods across all 60 trials are displayed in Table I.

Table I. Root mean square error of all planar rotation extraction methods.

Figure 9. Progression of the sagittal knee angle for one trial of a single participant during a gait cycle.

Clearly, the Quaternion extraction exhibits the lowest deviation of $2.8^{\circ }$ , closely approaching the gold standard. Conversely, all Euler-angle decompositions present a significantly higher RMSE, with a particularly pronounced deviation in Euler sequences where the extracted angle about the global $z$ -axis is the middle angle of the sequence ( $x$ - $z$ - $y$ and $y$ - $z$ - $x$ ).

From Figure 9, one can observe that the Quaternion extraction follows best the Vicon measurement, up to some wobbling effects, which are attributed to gross soft-tissue motions around heel strike (yellow area), and a slight overshoot at maximal knee flexion (dark gray area), which is attributed to skin stretching and muscle bulging at these stages. On the other hand, one can see strong deviations in the Euler extraction methods, which can reach an order of magnitude of 10 $^{\circ }$ in the worst case. Also, the wobbling effects are more pronounced in the Euler extraction methods, which is due to their sensibility to slight tilt rotations of the IMU sensor on the skin, which are better compensated by the Quaternion extraction method. Similar results were obtained when determining first the 3D relative orientation between the two IMUs and then extracting the planar sagittal knee angle from this, as shown in Figure 10.

Figure 10. Progression of the sagittal knee angle determined from relative IMU orientation.

Altogether, the newly introduced Quaternion extraction is visibly more robust and precise for knee flexural angle determination than the traditional Euler extraction methods, and, due to its simplicity and uniqueness, a promising technology for also other sagittal anatomical angle extractions (e.g., hip, elbow, and torso) from IMU measurements for applications in biomechanics.

Future research might focus on enhancing sensor placement, developing software-based detection and compensation mechanisms for soft-tissue artifact effects, and determining the sagittal plane during free walking or running from measurement data.

5. Conclusion

A review of existing methods for extracting “planar” rotation components from 3D rotations, such as those needed, for example, for the determination of sagittal anatomical joint rotations from IMU measurements in biomechanics or of parasitic motions in parallel kinematics manipulators, is presented. The existing methods can be grouped into two categories: one using a decomposition of the 3D rotation into a sequence of elementary follow-up rotations about orthogonal coordinate axes for different rotation sequences, identifying one of them as the sought planar rotation (termed Euler-angle decomposition methods), and one using a projection of one of the coordinate axes of the rotated frame onto the target plane and determining the planar rotation as the angle between the projected moving axis and an axis of the base frame (termed coordinate decomposition methods). Both groups are shown to be equivalent. Moreover, it is shown that the existing methods are highly sensitive with respect to the choice of the sequence of partial rotations, and that they lead to ambiguities when the rotations are not small and inconsistencies when a rotation about an axis in the plane takes place. A novel approach to resolve these difficulties is presented based on the decomposition of the overall rotation Quaternion in a component normal to the target plane and a component in the plane. Three remarkable properties of the new Quaternion planar rotation extraction method are proven: (1) that it corresponds to a scaled projection of the full Quaternion vector onto the target planar rotation axis (divided by the original Quaternion real part); (2) that it is unique; and (3) that it is commutative with respect to the location of the planar rotation in the decomposition of the full rotation concerning the value of the extracted planar rotation angle, that is, if it placed at the beginning (inertially-fixed) or at the end (body-fixed). This makes the Quaternion extraction method a simple and effective quasi-vector space projection with respect to the extracted planar rotation angle, a property that has not been shown or known before. Two case examples, one for a virtual rotation of a box driven by a time-dependent angular velocity input, the second involving an experimental setup for measuring the sagittal knee angle by IMUs compared to the gold-standard Vicon, illustrate the inconsistencies of the existing approaches and validate the superiority of new Quaternion decomposition approach. Moreover, a new method for determining the orientation of a box with checkerboard face by two CCD cameras is presented, which yields higher accuracy by a factor of 5 with respect to existing P3Pp methods.