Linearized Expressions of 3D Rotational Motion

1. Background

Rotational motion [1] in the three-dimensional physical world refers to the circular movement of an object around a fixed axis [2]. It is a basic concept and commonly used in physics, computer graphics, computer vision, robotics, aerospace, and other scientific fields [3,4,5,6]. More specifically, in many practical applications, rotational motion is not only an important means of constructing the motion of the target system [7,8] but also an important state parameter to be estimated [9,10,11]. However, rotation motion is typically formulated as a non-linear and non-convex process [11,12]. It not only makes some linear algorithms, such as Kalman filter [13,14], impossible to directly apply [15], but also makes solving the related problem more sophisticated and difficult [12,16].

In this paper, a new perspective is proposed that the rotation motion can be linearized without dropping any constraints and increasing any singularities. It is particularly worth pointing out that there are no special and advanced techniques, such as Lie groups and Lie algebras [17,18], for the methods used in this paper. Only some basic knowledge of linear algebra [19,20,21] is required to understand the ideas presented in this paper.

2. Rotation Motion

Mathematically, given a point $\mathit{x}\in {\mathbb{R}}^3$, if it is rotated about a fixed unit vector $\mathit{r}\in {\mathbb{S}}^2$ through the origin and the rotated angle is $\theta$, then the rotation motion can be formulated by multiplication of a matrix and a vector,

$$\mathbf{R}\mathit{x}\to \mathit{y}$$

(1)

where $\mathit{y}\in {\mathbb{R}}^3$ is the rotated point by $\mathit{x}$ and $\mathbf{R}$ is a $3\times 3$ rotation matrix. There are some properties regarding the rotation motion [22],

• Shape and Size: Rotation does not alter the shape or size of the object. More Specifically,
  • $\parallel \mathbf{R}\mathit{x}\parallel =\parallel \mathbf{R}\parallel ·\parallel \mathit{x}\parallel =\parallel \mathit{x}\parallel$. It means the length of the point $\mathit{x}$ will not be changed by rotation motion. Without loss of generality, the rotation motion in this paper is studied with unit norm vectors, i.e., $\mathit{x}\in {\mathbb{S}}^2$ on the unit sphere surface, which represents a unit direction vector in the 3D physical world.
  • $\angle \left(\mathbf{R}{\mathit{x}}_1,\mathbf{R}{\mathit{x}}_2\right)=\angle \left({\mathit{x}}_1,{\mathit{x}}_2\right)$ where ${\mathit{x}}_1,{\mathit{x}}_2\in {\mathbb{S}}^2$ are any two unit vectors. It means the angle (structure) of the object is unchanged.
• Axis of Rotation: Rotational motion occurs around a fixed axis. Given a vector $\mathit{x}\in {\mathbb{S}}^2$, if $\mathbf{R}\mathit{x}=\mathit{x}$, it means the after the rotation $\mathbf{R}$, the vector direction $\mathit{x}$ remains unchanged, and the vector $\mathit{x}$ must be parallel to the rotation axis $\mathit{r}$, which is from Euler’s rotation theorem and screw theory [23,24,25,26]. In addition,

  $$\mathbf{R}\mathit{x}=\mathit{x}\Rightarrow \left(\mathbf{R}-{\mathbf{I}}_{3\times 3}\right)\mathit{x}=0$$

  (2)

  where ${\mathbf{I}}_{3\times 3}$ is the 3 × 3 identity matrix. Algebraically, Eq.(2) means the rotation axis direction $\mathit{r}$ lies in the null space of $\mathbf{R}-{\mathbf{I}}_{3\times 3}$. Alternatively, let $\mathbf{R}\mathit{x}=\lambda \mathit{x}$, then the rotation axis direction $\mathit{r}$ is an eigenvector of $\mathbf{R}$ corresponding to the eigenvalue $\lambda =1$.

At first glance, the rotation motion is a linear process, which might be misled by the matrix representation Eq.(1). Unfortunately, $\mathbf{R}$ is not a regular matrix in ${\mathbb{R}}^{3\times 3}$. A rotation only has three degrees of freedom, since it can be determined by a rotation axis $\mathit{r}\in {\mathbb{S}}^2$ and a rotation angle $\theta \in [0,\pi ]$. Formally, $\mathbf{R}\in \mathbb{SO}\left(3\right)$, and a rotation can be represented by a matrix and non-linear constraints,

$$\mathbf{R}=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right],s.t.{\mathbf{R}}^T\mathbf{R}={\mathbf{I}}_{3\times 3},det\left(\mathbf{R}\right)=+1$$

(3)

According to the famous Rodrigues’ rotation formula [27,28], given the rotation axis $\mathit{r}={\left[\begin{array}{ccc}{r}_1& r_2& r_3\end{array}\right]}^T\in {\mathbb{S}}^2$ and rotation angle $\theta$, the matrix can be constructed by

$$\mathbf{R}=exp\left(\theta {\left[\mathit{r}\right]}_{\times }\right)={\mathbf{I}}_{3\times 3}+sin\left(\theta \right){\left[\mathit{r}\right]}_{\times }+\left(1-cos\left(\theta \right)\right){\left[\mathit{r}\right]}_{\times }^2,{\left[\mathit{r}\right]}_{\times }=\left[\begin{array}{ccc}0& -r_3& r_2\\ r_3& 0& -r_1\\ -r_2& r_1& 0\end{array}\right]$$

(4)

Moreover, a rotation matrix can be represented by a quaternion vector with unit length [29,30], such that $\mathit{q}={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T\in {\mathbb{S}}^3$. Specifically,

$$\mathbf{R}=\left[\begin{array}{ccc}{q}_0^2+q_1^2-q_2^2-q_3^2& 2\left(q_1{q}_2-q_0{q}_3\right)& 2\left(q_1{q}_3+q_0{q}_2\right)\\ 2\left(q_1{q}_2+q_0{q}_3\right)& q_0^2-q_1^2+q_2^2-q_3^2& 2\left(q_2{q}_3-q_0{q}_1\right)\\ 2\left(q_1{q}_3-q_0{q}_2\right)& 2\left(q_2{q}_3+q_0{q}_1\right)& q_0^2-q_1^2-q_2^2+q_3^2\end{array}\right]$$

(5)

It is worth noting that many different forms of parameter representation can represent the same rotation, and each form of representation has its appropriate application scenario [2,31,32,33]. For example, since a rotation can be determined by a fixed rotation axis $\mathit{r}\in {\mathbb{S}}^2$ and a rotation angle $\theta \in \left[0,\pi \right]$, then the rotation can be represented by axis–angle representation $\theta ·\mathit{r}$ [34], which geometrically is a point in a 3D ball of radius $\pi$ [28] and has been widely used in pose estimation problems [28,35,36]. Nonetheless, throughout most of this paper, quaternions are used to represent rotations [29,30] and thus to illustrate the proposed linearization theory.

Anyone who has ever used any other parametrization of the rotation group will, within hours of taking up the quaternion parametrization, lament his or her misspent youth [2,31].

3. Linear Expressions for Rotation Motion

Given $\mathit{a}=\left[\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right]\in {\mathbb{S}}^2$ and a rotation $\mathbf{R}\in \mathbb{SO}\left(3\right)$, the rotation motion can be formulated as

$$\mathbf{R}\mathit{a}=\mathit{b}\iff \angle \left(\mathbf{R}\mathit{a},\mathit{b}\right)=0\iff {\mathit{b}}^T\mathbf{R}\mathit{a}=1$$

(6)

where $\mathit{b}=\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]\in {\mathbb{S}}^2$ is rotated by $\mathit{a}$. This formulation is a special case of the next expression for rotational motion. Formally, given $\mathit{a},\mathit{b}\in {\mathbb{S}}^2$ and $\mathbf{R}\in \mathbb{SO}\left(3\right)$ whose quaternion expression is $\mathit{q}\in {\mathbb{S}}^3$,

$$\angle \left(\mathbf{R}\mathit{a},\mathit{b}\right)=\tau \iff {\mathit{b}}^T\mathbf{R}\mathit{a}=cos\left(\tau \right)\iff {\mathit{q}}^T\mathbf{M}\mathit{q}=cos\left(\tau \right)$$

(7)

Here $\mathbf{M}$ is constructed by $\mathit{a}$ and $\mathit{b}$ [29,37]. Specifically,

$$\mathbf{M}=\left[\begin{array}{cccc}\hfill 0& -a_1& -a_2\hfill & -a_3\\ \hfill a_1& 0& a_3\hfill & -a_2\\ \hfill a_2& -a_3& 0\hfill & a_1\\ \hfill a_3& a_2& -a_1\hfill & 0\end{array}{]}^T\right[\begin{array}{cccc}\hfill 0& -b_1& -b_2\hfill & -b_3\\ \hfill b_1& 0& -b_3\hfill & b_2\\ \hfill b_2& b_3& 0\hfill & -b_1\\ \hfill b_3& -b_2& b_1\hfill & 0\end{array}]$$

(8)

When $\tau =0$, Eq.(7) degenerates back to Eq.(6).

By conducting eigendecomposition of the real symmetric matrix $\mathbf{M}=\mathbf{U}\Lambda {\mathbf{U}}^T$ [21,38], where and $\mathbf{U}=\left[\begin{array}{cccc}{\mathit{u}}_0& {\mathit{u}}_1& {\mathit{u}}_2& {\mathit{u}}_3\end{array}\right]\in {\mathbb{R}}^{4\times 4}$ is an orthogonal matrix $\mathbf{U}{\mathbf{U}}^T={\mathbf{I}}_{4\times 4}$. Observe

$$det\left(\mathbf{M}-\lambda {\mathbf{I}}_{4\times 4}\right)=0\Rightarrow {\left({\lambda }^2-1\right)}^2=0$$

(9)

Therefore, the eigenvalue of $\mathbf{M}$ is

$$\Lambda =\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right]$$

(10)

Accordingly, ${\mathit{q}}^T\mathbf{U}\Lambda {\mathbf{U}}^T\mathit{q}=cos\left(\tau \right)$. Let $\Gamma ={\left[{\Gamma }_0,{\Gamma }_1,{\Gamma }_2,{\Gamma }_3\right]}^T={\mathbf{U}}^T\mathit{q}$, then $\mathit{q}={\Gamma }_0{\mathit{u}}_0+{\Gamma }_1{\mathit{u}}_1+{\Gamma }_2{\mathit{u}}_2+{\Gamma }_3{\mathit{u}}_3$. Observe

$$\left\{\begin{array}{cc}\hfill {\Gamma }^T\Gamma ={\mathit{q}}^T\mathbf{U}{\mathbf{U}}^T\mathit{q}={\Gamma }_0^2+{\Gamma }_1^2+{\Gamma }_2^2+{\Gamma }_3^2& =1\hfill \\ \hfill {\mathit{q}}^T\mathbf{U}\Lambda {\mathbf{U}}^T\mathit{q}={\Gamma }^T\Lambda \Gamma ={\Gamma }_0^2+{\Gamma }_1^2-{\Gamma }_2^2-{\Gamma }_3^2& =cos\left(\tau \right)\hfill \end{array}\right.\Rightarrow \left\{\begin{array}{c}\hfill {\Gamma }_0^2+{\Gamma }_1^2={\displaystyle \frac{1}{2}}\left(1+cos\left(\tau \right)\right)\\ \hfill {\Gamma }_2^2+{\Gamma }_3^2={\displaystyle \frac{1}{2}}\left(1-cos\left(\tau \right)\right)\end{array}\right.$$

(11)

Therefore, the possible rotation $\mathit{q}$ can be reformulated linearly

$$\begin{array}{cc}\hfill \mathit{q}=& {\Gamma }_0{\mathit{u}}_0+{\Gamma }_1{\mathit{u}}_1+{\Gamma }_2{\mathit{u}}_2+{\Gamma }_3{\mathit{u}}_3\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill =& \frac{\sqrt{2\left(1+cos\left(\tau \right)\right)}}{2}\left(cos\left(\varphi \right){\mathit{u}}_0+sin\left(\varphi \right){\mathit{u}}_1\right)+\frac{\sqrt{2\left(1-cos\left(\tau \right)\right)}}{2}\left(cos\left(\psi \right){\mathit{u}}_2+sin\left(\psi \right){\mathit{u}}_3\right)\hfill \end{array}$$

(13)

where $\varphi ,\psi \in [-\pi ,\pi ]$. This result clearly shows that the rotational motion can be expressed in linear forms.

Geometrically, if $\angle \left(\mathbf{R}\mathit{a},\mathit{b}\right)=\tau$, then all possible rotation $\mathit{q}$ that satisfy the constraint will be in a torus in ${\mathbb{S}}^3$. For better understanding, we draw some special cases, i.e., $\tau =\{0,\frac{\pi }{9},\frac{3\pi }{9},\frac{6\pi }{9},\frac{8\pi }{9}\}$ via stereographic projection [39,40], as in Figure 1. It is worth noting that when $\tau =\{0,\pi \}$, i.e., $\mathbf{R}\mathit{a}=\pm \mathit{b}$, the torus in ${\mathbb{S}}^3$ is degenerated into a one-dimensional circle, more specifically, a great circle [41] in ${\mathbb{S}}^3$. Formally,

$$\begin{array}{cccc}\hfill & ➀:\mathrm{when}\tau =0\hfill & \hfill \Rightarrow & \mathit{q}=cos\left(\varphi \right){\mathit{u}}_0+sin\left(\varphi \right){\mathit{u}}_1\hfill \\ \hfill & ➁:\mathrm{when}\tau =\pi \hfill & \hfill \Rightarrow & \mathit{q}=cos\left(\psi \right){\mathit{u}}_2+sin\left(\psi \right){\mathit{u}}_3\hfill \end{array}$$

(14)

In addition, when $\tau =\frac{\pi }{2}$, the torus is becoming to the Clifford torus in ${\mathbb{S}}^3$ [42,43], since it has the equal diameters in cross directions. Specifically,

$$\mathit{q}=\frac{1}{2}\left(cos\left(\varphi \right){\mathit{u}}_0+sin\left(\varphi \right){\mathit{u}}_1\right)+\frac{1}{2}\left(cos\left(\psi \right){\mathit{u}}_2+sin\left(\psi \right){\mathit{u}}_3\right)$$

(15)

These special cases will be discussed in depth in subsequent sections.

4. Special Case I: Great Circle in ${\mathbb{S}}^3$

If $\tau =0$, which means $\mathbf{R}\mathit{a}=\mathit{b}$, it is one of the most common cases in practical applications. It is worth noting that there are two facts regarding $\mathbf{R}\mathit{a}=\mathit{b}$,

• Not all rotations in $\mathbb{SO}\left(3\right)$ can rotate $\mathit{a}$ to $\mathit{b}$.
• There is more than one rotation that satisfies the given rotation motion.

Then according to Eq. (14) all rotations that satisfy $\mathbf{R}\mathit{a}=\mathit{b}$ should be in

$$\mathit{q}=cos\left(\varphi \right){\mathit{u}}_0+sin\left(\varphi \right){\mathit{u}}_1$$

(16)

where ${\mathit{u}}_0$ and ${\mathit{u}}_{\mathbf{1}}$ are the basis vectors in ${\mathbb{S}}^3$, and they depend on $\mathbf{M}$, i.e., $\mathit{a}$ and $\mathit{b}$. Geometrically speaking, $\mathit{q}$ is a unit vector in ${\mathbb{R}}^4$, which is a linear combination of the two vectors ${\mathit{u}}_0$ and ${\mathit{u}}_1$. In addition, ${\mathit{q}}^T\mathit{q}=1$, which means $\mathit{q}$ is on the surface of the unit quaternion sphere. Therefore, the region where all possible rotations satisfy a specific rotation motion is a great circle in ${\mathbb{S}}^3$ [41]. Moreover, since ${\mathbf{U}}^T\mathbf{U}={\mathbf{I}}_{4\times 4}$, there will be

$$\left\{\begin{array}{c}\hfill \left(cos\left(\varphi \right){\mathit{u}}_0^T+sin\left(\varphi \right){\mathit{u}}_1^T\right){\mathit{u}}_2=0\\ \hfill \left(cos\left(\varphi \right){\mathit{u}}_0^T+sin\left(\varphi \right){\mathit{u}}_1^T\right){\mathit{u}}_3=0\end{array}\right.\Rightarrow \left\{\begin{array}{c}\hfill {\mathit{q}}^T{\mathit{u}}_2=0\\ \hfill {\mathit{q}}^T{\mathit{u}}_3=0\end{array}\right.$$

(17)

Up until this point, the rotation motion is linearly expressed, which is closely related to $\mathbf{U}$.

In addition,

$$\mathbf{M}{\mathit{u}}_0=1·{\mathit{u}}_0;\mathbf{M}{\mathit{u}}_1=1·{\mathit{u}}_1;\mathbf{M}{\mathit{u}}_2=-1·{\mathit{u}}_2;\mathbf{M}{\mathit{u}}_3=-1·{\mathit{u}}_3$$

(18)

It is easy to confirm [44,45],

$$\begin{array}{cc}\hfill \mathbf{M}\left[\begin{array}{cc}0& -{\left(\mathit{a}-\mathit{b}\right)}^T\\ -\left(\mathit{a}-\mathit{b}\right)& {\left[\mathit{a}+\mathit{b}\right]}_{\times }^T\end{array}\right]=& \left[\begin{array}{cc}0& -{\left(\mathit{a}-\mathit{b}\right)}^T\\ -\left(\mathit{a}-\mathit{b}\right)& {\left[\mathit{a}+\mathit{b}\right]}_{\times }^T\end{array}\right]\left[\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right]\stackrel{·}{=}-1·{\Omega }^{-}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \mathbf{M}\left[\begin{array}{cc}0& -{\left(\mathit{a}+\mathit{b}\right)}^T\\ -\left(\mathit{a}+\mathit{b}\right)& {\left[\mathit{a}-\mathit{b}\right]}_{\times }^T\end{array}\right]=& \left[\begin{array}{cc}0& -{\left(\mathit{a}+\mathit{b}\right)}^T\\ -\left(\mathit{a}+\mathit{b}\right)& {\left[\mathit{a}-\mathit{b}\right]}_{\times }^T\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\stackrel{·}{=}+1·{\Omega }^{+}\hfill \end{array}$$

(20)

Therefore, any row of ${\Omega }^{-}$ can be the eigenvector of $\mathbf{M}$ with eigenvalue $-1$. Similarly, any row of ${\Omega }^{+}$ can be the eigenvector of $\mathbf{M}$ with eigenvalue 1.

Accordingly, (1) any row of ${\Omega }^{-}$ is a linear combination of ${\mathit{u}}_0$ and ${\mathit{u}}_1$. (2) any row of ${\Omega }^{+}$ is a linear combination of ${\mathit{u}}_2$ and ${\mathit{u}}_3$. Consequently, $\mathrm{rank}\left({\Omega }^{-}\right)=\mathrm{rank}\left({\Omega }^{+}\right)=2$; ${\Omega }^{-}$ and ${\Omega }^{+}$ are not orthogonal matrices. Therefore, ${\Omega }^{-}\mathit{q}=\mathbf{0}$, more specifically,

$$\left[\begin{array}{cccc}0& a_1-b_1& a_2-b_2& a_3-b_3\\ -\left(a_1-b_1\right)& 0& -\left(a_3+b_3\right)& a_2+b_2\\ -\left(a_2-b_2\right)& a_3+b_3& 0& -\left(a_1+b_1\right)\\ -\left(a_3-b_3\right)& -\left(a_2+b_2\right)& a_1+b_1& 0\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]=\mathbf{0}\iff \left[\begin{array}{cc}0& {\left(\mathit{a}-\mathit{b}\right)}^T\\ -\left(\mathit{a}-\mathit{b}\right)& {\left[\mathit{a}+\mathit{b}\right]}_{\times }\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]=\mathbf{0}$$

(21)

Now the linear expression of $\mathbf{R}\mathit{a}=\mathit{b}$ is constructed.

Furthermore, if more than one rotation motion constraints are given, such that

$$\left\{\begin{array}{cc}\hfill \mathbf{R}{\mathit{a}}_1& ={\mathit{b}}_1\hfill \\ \hfill \mathbf{R}{\mathit{a}}_2& ={\mathit{b}}_2\hfill \\ \hfill \cdots \\ \hfill \mathbf{R}{\mathit{a}}_N& ={\mathit{b}}_N\hfill \end{array}\right.\stackrel{\mathrm{linearization}}{\to }\left\{\begin{array}{c}\hfill {\Omega }_1^{-}\mathit{q}=\mathbf{0}\\ \hfill {\Omega }_2^{-}\mathit{q}=\mathbf{0}\\ \hfill \cdots \\ \hfill {\Omega }_N^{-}\mathit{q}=\mathbf{0}\end{array}\right.⟹\left[\begin{array}{c}{\Omega }_1^{-}\\ {\Omega }_2^{-}\\ \cdots \\ {\Omega }_N^{-}\end{array}\right]\mathit{q}=\mathbf{0}$$

(22)

Eventually, the corresponding linear system is constructed. It must be acknowledged that there is still a constraint of length 1, i.e., ${\mathit{q}}^T\mathit{q}=1$. However, the scale is not particularly important for homogeneous linear equations.

4.1. Alternative Way to Obtain Great Circle in ${\mathbb{S}}^3$

In this part, a geometrical approach to obtain a great circle in ${\mathbb{S}}^3$ is presented. Given the rotation motion $\mathbf{R}\mathit{a}=\mathit{b}$, the possible rotation can be formulated as [9,34]

$$\mathbf{R}={\mathbf{R}}_{\mathit{b}}\left(\alpha \right){\mathbf{R}}_a^b$$

(23)

where ${\mathbf{R}}_{\mathit{b}}\left(\alpha \right)$ is a rotation around rotated axis is $\mathit{b}$ and rotated angle is $\alpha$; ${\mathbf{R}}_{\mathit{a}}^{\mathit{b}}$ is a rotation that can rotate $\mathit{a}$ to $\mathit{b}$. To use quaternion to represent rotation [30],

$$\begin{array}{cc}\hfill {\mathbf{R}}_{\mathit{b}}\left(\alpha \right)& :=\left[cos\left(\frac{\alpha }{2}\right),sin\left(\frac{\alpha }{2}\right)\mathit{b}\right]\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathbf{R}}_{\mathit{a}}^{\mathit{b}}& :=\left[q_0,{\left[q_1,q_2,q_3\right]}^T\right]=\left[q_0,{\mathit{q}}_{123}\right]\hfill \end{array}$$

(25)

The formula for the multiplication of two quaternions is [29]

$$(r_1,\overrightarrow{v_1})(r_2,\overrightarrow{v_2})=(r_1{r}_2-\overrightarrow{v_1}\overrightarrow{v_2},r_1\overrightarrow{v_2}+r_2\overrightarrow{v_1}+\overrightarrow{v_1}\times \overrightarrow{v_2})$$

(26)

Then,

$$\begin{array}{c}\hfill \mathbf{R}:=\left[cos\left(\frac{\alpha }{2}\right)q_0-sin\left(\frac{\alpha }{2}\right){\mathit{b}}^T{\mathit{q}}_{123},q_{0}sin\left(\frac{\alpha }{2}\right)\mathit{b}+cos\left(\frac{\alpha }{2}\right){\mathit{q}}_{123}+sin\left(\frac{\alpha }{2}\right)\mathit{b}\times {\mathit{q}}_{123}\right]\end{array}$$

(27)

Therefore,

$$\begin{array}{c}\hfill \mathit{q}=\left[\begin{array}{c}\hfill q_0\\ \hfill {\mathit{q}}_{123}\end{array}\right]cos\left(\frac{\alpha }{2}\right)+\left[\begin{array}{c}\hfill -{\mathit{b}}^T{\mathit{q}}_{123}\\ \hfill q_0\mathit{b}+\mathit{b}\times {\mathit{q}}_{123}\end{array}\right]sin\left(\frac{\alpha }{2}\right)\stackrel{·}{=}{\mathit{q}}_{b1}cos\left(\frac{\alpha }{2}\right)+{\mathit{q}}_{b2}sin\left(\frac{\alpha }{2}\right)\end{array}$$

(28)

Observe

$$\begin{array}{cc}\hfill & {\mathit{q}}_{b1}^T{\mathit{q}}_{b1}=q_0^2+q_1^2+q_2^2+q_3^2=1\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & \begin{array}{cc}\hfill {\mathit{q}}_{b2}^T{\mathit{q}}_{b2}=& {\left({\mathit{b}}^T{\mathit{q}}_{123}\right)}^2+{\left(q_0\mathit{b}\right)}^2+{\left(\mathit{b}\times {\mathit{q}}_{123}\right)}^2={\mathit{q}}_{123}^T{\mathit{q}}_{123}{\left(cos\angle (\mathit{b},{\mathit{q}}_{123})\right)}^2+q_0^2+{\mathit{q}}_{123}^T{\mathit{q}}_{123}{\left(sin\angle (\mathit{b},{\mathit{q}}_{123})\right)}^2\hfill \\ \hfill =& q_0^2+q_1^2+q_2^2+q_3^2=1\hfill \end{array}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & {\mathit{q}}_{b1}^T{\mathit{q}}_{b2}=-q_0{\mathit{b}}^T{\mathit{q}}_{123}+{\mathit{q}}_{123}\left(q_0\mathit{b}+\mathit{b}\times {\mathit{q}}_{123}\right)=0\hfill \end{array}$$

(31)

Therefore, ${\mathit{q}}_{b1}$ and ${\mathit{q}}_{b2}$ are two orthonormal bases in ${\mathbb{R}}^4$, which can be considered as ${\mathit{u}}_0$ and ${\mathit{u}}_1$ in Eq. (16). Consequently, $\mathit{q}$ lies in the great circle in ${\mathbb{S}}^3$.

5. Special Case II: Clifford Torus in ${\mathbb{S}}^3$

5.1. Alternative Way to Obtain Clifford Torus in ${\mathbb{S}}^3$

In the previous text, we provide a general derivation for the connection between rotation constraint, i.e., ${\mathit{b}}^T\mathbf{R}\mathit{a}=0$ and Clifford torus. Here we show an alternative way to obtain Clifford torus from rotation constraint.

Given a rotation constraint ${\mathit{b}}^T\mathbf{R}\mathit{a}=0$, the possible rotation $\mathbf{R}$ can be (see Figure 2)

$$\mathbf{R}\left(\alpha ,\beta \right)={\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0,\alpha ,\beta \in \left[-\pi ,\pi \right]$$

(32)

where ${\mathit{v}}_2=\mathit{b}$, ${\mathit{v}}_1^T{\mathit{v}}_1=1$ and ${\mathit{v}}_1^T{\mathit{v}}_2=0$; ${\mathbf{R}}_0$ is a rotation that rotates $\mathit{a}$ to ${\mathit{v}}_1$; ${\mathbf{R}}_{v_1}\left(\alpha \right)=exp\left(\alpha {\left[{\mathit{v}}_1\right]}_{\times }\right)$ and ${\mathbf{R}}_{v_2}\left(\beta \right)=exp\left(\beta {\left[{\mathit{v}}_2\right]}_{\times }\right)$ are rotations around directions ${\mathit{b}}_1$ and $\mathit{b}$, respectively. The rotation motion can be as follows,

$${\mathit{v}}_1={\mathbf{R}}_0\mathit{a},⟶{\mathit{v}}_1={\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0\mathit{a},⟶{\mathit{b}}^T\left({\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0\mathit{a}\right)=0.$$

(33)

To use quaternion to represent rotation [30]

$$\begin{array}{cc}\hfill {\mathbf{R}}_{v_2}\left(\beta \right):=& {\mathit{q}}_{v_2}\left(\beta \right)={\left[\begin{array}{cc}cos\left({\displaystyle \frac{\beta }{2}}\right)& {\mathit{v}}_2^{T}sin\left({\displaystyle \frac{\beta }{2}}\right)\end{array}\right]}^T\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\mathbf{R}}_{v_1}\left(\alpha \right):=& {\mathit{q}}_{v_1}\left(\alpha \right)={\left[\begin{array}{cc}cos\left({\displaystyle \frac{\alpha }{2}}\right)& {\mathit{v}}_1^{T}sin\left({\displaystyle \frac{\alpha }{2}}\right)\end{array}\right]}^T\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\mathbf{R}}_0:=& {\left[\begin{array}{cc}{q}_0& {\mathit{q}}_{123}^T\end{array}\right]}^T\hfill \end{array}$$

(36)

The formula for the multiplication of two quaternions is [29]

$$(r_1,\overrightarrow{v_1})(r_2,\overrightarrow{v_2})=(r_1{r}_2-\overrightarrow{v_1}\overrightarrow{v_2},r_1\overrightarrow{v_2}+r_2\overrightarrow{v_1}+\overrightarrow{v_1}\times \overrightarrow{v_2})$$

(37)

Define ${\mathit{v}}_3={\mathit{v}}_1\times {\mathit{v}}_2$, then

$${\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right):={\left[\begin{array}{cc}cos\left({\displaystyle \frac{\beta }{2}}\right)cos\left({\displaystyle \frac{\alpha }{2}}\right)& {\mathit{v}}_2^{T}sin\left({\displaystyle \frac{\beta }{2}}\right)cos\left({\displaystyle \frac{\alpha }{2}}\right)+{\mathit{v}}_1^{T}sin\left({\displaystyle \frac{\alpha }{2}}\right)cos\left({\displaystyle \frac{\beta }{2}}\right)-{\mathit{v}}_3^{T}sin\left({\displaystyle \frac{\beta }{2}}\right)sin\left({\displaystyle \frac{\alpha }{2}}\right)\end{array}\right]}^T$$

(38)

In addition, let $w_1={\displaystyle \frac{1}{2}}(\alpha -\beta )$ and $w_2={\displaystyle \frac{1}{2}}(\alpha +\beta )$,

$$\begin{array}{cc}\hfill cos\left({\displaystyle \frac{\beta }{2}}\right)cos\left({\displaystyle \frac{\alpha }{2}}\right)=& {\displaystyle \frac{1}{2}}\left(cos\left({\displaystyle \frac{\alpha -\beta }{2}}\right)+cos\left({\displaystyle \frac{\alpha +\beta }{2}}\right)\right)={\displaystyle \frac{1}{2}}\left(cos\left(w_1\right)+cos\left(w_2\right)\right)\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill sin\left({\displaystyle \frac{\beta }{2}}\right)cos\left({\displaystyle \frac{\alpha }{2}}\right)=& {\displaystyle \frac{1}{2}}\left(sin\left({\displaystyle \frac{\alpha +\beta }{2}}\right)-sin\left({\displaystyle \frac{\alpha -\beta }{2}}\right)\right)={\displaystyle \frac{1}{2}}\left(sin\left(w_2\right)-sin\left(w_1\right)\right)\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill sin\left({\displaystyle \frac{\alpha }{2}}\right)cos\left({\displaystyle \frac{\beta }{2}}\right)=& {\displaystyle \frac{1}{2}}\left(sin\left({\displaystyle \frac{\alpha +\beta }{2}}\right)+sin\left({\displaystyle \frac{\alpha -\beta }{2}}\right)\right)={\displaystyle \frac{1}{2}}\left(sin\left(w_2\right)+sin\left(w_1\right)\right)\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill sin\left({\displaystyle \frac{\beta }{2}}\right)sin\left({\displaystyle \frac{\alpha }{2}}\right)=& {\displaystyle \frac{1}{2}}\left(cos\left({\displaystyle \frac{\alpha -\beta }{2}}\right)-cos\left({\displaystyle \frac{\alpha +\beta }{2}}\right)\right)={\displaystyle \frac{1}{2}}\left(cos\left(w_1\right)-cos\left(w_2\right)\right)\hfill \end{array}$$

(42)

Then,

$$\begin{array}{cc}\hfill {\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right):=& {\displaystyle \frac{1}{2}}{\left[\begin{array}{cc}cos\left(w_1\right)+cos\left(w_2\right)& sin\left(w_1\right)\left({\mathit{v}}_1^T-{\mathit{v}}_2^T\right)-cos\left(w_1\right){\mathit{v}}_3^T+sin\left(w_2\right)\left({\mathit{v}}_1^T+{\mathit{v}}_2^T\right)+cos\left(w_2\right){\mathit{v}}_3^T\end{array}\right]}^T\hfill \end{array}$$

(43)

Therefore,

$$\begin{array}{cc}\hfill & \mathbf{R}\left(\alpha ,\beta \right)={\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill :=& {\displaystyle \frac{1}{2}}\left[\begin{array}{c}\left(cos\left(w_1\right)+cos\left(w_2\right)\right)q_0-\left(sin\left(w_1\right)\left({\mathit{v}}_1^T-{\mathit{v}}_2^T\right)-cos\left(w_1\right){\mathit{v}}_3^T+sin\left(w_2\right)\left({\mathit{v}}_1^T+{\mathit{v}}_2^T\right)+cos\left(w_2\right){\mathit{v}}_3^T\right){\mathit{q}}_{123}\\ \begin{array}{c}\hfill \left(cos\left(w_1\right)+cos\left(w_2\right)\right){\mathit{q}}_{123}+\left(sin\left(w_1\right)\left({\mathit{v}}_1-{\mathit{v}}_2\right)-cos\left(w_1\right){\mathit{v}}_3+sin\left(w_2\right)\left({\mathit{v}}_1+{\mathit{v}}_2\right)+cos\left(w_2\right){\mathit{v}}_3\right)q_0+\\ \hfill \left(sin\left(w_1\right)\left({\mathit{v}}_1-{\mathit{v}}_2\right)-cos\left(w_1\right){\mathit{v}}_3+sin\left(w_2\right)\left({\mathit{v}}_1+{\mathit{v}}_2\right)+cos\left(w_2\right){\mathit{v}}_3\right)\times {\mathit{q}}_{123}\end{array}\end{array}\right]\hfill \end{array}$$

(45)

$$$$

(46)

$$\begin{array}{cc}\hfill =& {\displaystyle \frac{1}{2}}\left[\begin{array}{cccc}{\mathit{k}}_0& {\mathit{k}}_1& {\mathit{k}}_2& {\mathit{k}}_3\end{array}\right]\left[\begin{array}{c}cos\left(w_1\right)\\ sin\left(w_1\right)\\ cos\left(w_2\right)\\ sin\left(w_2\right)\end{array}\right],w_1={\displaystyle \frac{1}{2}}(\alpha -\beta ),\mathrm{and},w_2={\displaystyle \frac{1}{2}}(\alpha +\beta )\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill =& {\displaystyle \frac{1}{2}}\left({\mathit{k}}_{0}cos\left(w_1\right)+{\mathit{k}}_{1}sin\left(w_1\right)+{\mathit{k}}_{2}cos\left(w_2\right)+{\mathit{k}}_{3}sin\left(w_2\right)\right)\hfill \end{array}$$

(48)

where

$$\begin{array}{cc}\hfill {\mathit{k}}_0=\left[\begin{array}{c}{q}_0+{\mathit{v}}_3^T{\mathit{q}}_{123}\\ {\mathit{q}}_{123}-q_0{\mathit{v}}_3-{\mathit{v}}_3\times {\mathit{q}}_{123}\end{array}\right]& {\mathit{k}}_1=\left[\begin{array}{c}-{\left({\mathit{v}}_1-{\mathit{v}}_2\right)}^T{\mathit{q}}_{123}\\ q_0\left({\mathit{v}}_1-{\mathit{v}}_2\right)+\left({\mathit{v}}_1-{\mathit{v}}_2\right)\times {\mathit{q}}_{123}\end{array}\right]\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\mathit{k}}_2=\left[\begin{array}{c}{q}_0-{\mathit{v}}_3^T{\mathit{q}}_{123}\\ {\mathit{q}}_{123}+q_0{\mathit{v}}_3+{\mathit{v}}_3\times {\mathit{q}}_{123}\end{array}\right]& {\mathit{k}}_3=\left[\begin{array}{c}-{\left({\mathit{v}}_1+{\mathit{v}}_2\right)}^T{\mathit{q}}_{123}\\ q_0\left({\mathit{v}}_1+{\mathit{v}}_2\right)+\left({\mathit{v}}_1+{\mathit{v}}_2\right)\times {\mathit{q}}_{123}\end{array}\right]\hfill \end{array}$$

(50)

It is easy to verify

$${\left[\begin{array}{cccc}{\mathit{k}}_0& {\mathit{k}}_1& {\mathit{k}}_2& {\mathit{k}}_3\end{array}\right]}^T\left[\begin{array}{cccc}{\mathit{k}}_0& {\mathit{k}}_1& {\mathit{k}}_2& {\mathit{k}}_3\end{array}\right]=\left[\begin{array}{cccc}2& 0& 0& 0\\ 0& 2& 0& 0\\ 0& 0& 2& 0\\ 0& 0& 0& 2\end{array}\right]$$

(51)

Geometrically, Eq. (44) describes a Clifford torus [42,46] in ${\mathbb{S}}^3$.

5.2. Intersections of Two Different Clifford Tori

Given two rotation constraints ${\mathit{b}}_1^T\mathbf{R}{\mathit{a}}_1=0$ and ${\mathit{b}}_2^T\mathbf{R}{\mathit{a}}_2=0$, the possible rotation can be considered in the intersections of two Clifford tori (see Figure 3). Here, a solution to solve the intersection line is presented. Specifically,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\hfill {\mathit{b}}_1^T\mathbf{R}{\mathit{a}}_1=0\\ \hfill {\mathit{b}}_2^T\mathbf{R}{\mathit{a}}_2=0\end{array}\right.\Rightarrow \left\{\begin{array}{cc}\hfill \mathbf{R}\left(\alpha ,\beta \right)& ={\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0\hfill \\ \hfill {\mathit{b}}_2^T\mathbf{R}{\mathit{a}}_2& =0\hfill \end{array}\right.\Rightarrow {\mathit{b}}_2^T{\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0{\mathit{a}}_2=0\end{array}$$

(52)

where ${\mathit{v}}_2={\mathit{b}}_1$ and ${\mathit{v}}_1\in {\mathbb{S}}^2$ is an arbitrary vector and satisfies ${\mathit{v}}_1^T{\mathit{v}}_2=0$; ${\mathbf{R}}_0$ is a specific arbitrary rotation that can rotate $\mathit{a}$ to ${\mathit{v}}_1$; ${\mathbf{R}}_{v_1}\left(\alpha \right)=exp\left(\alpha {\left[{\mathit{v}}_1\right]}_{\times }\right)$ and ${\mathbf{R}}_{v_2}\left(\beta \right)=exp\left(\beta {\left[{\mathit{v}}_2\right]}_{\times }\right)$.

According to Rodrigues’ rotation formula [3], we have

$$\begin{array}{c}\hfill {\mathbf{R}}_{v1}\left(\alpha \right)=exp\left(\alpha {\left[{\mathit{v}}_1\right]}_{\times }\right)=\mathbf{I}+sin\left(\alpha \right){\left[{\mathit{v}}_1\right]}_{\times }+\left(1-cos\left(\alpha \right)\right){\left[{\mathit{v}}_1\right]}_{\times }^2\end{array}$$

(53)

$$\begin{array}{c}\hfill {\mathbf{R}}_{v2}\left(\beta \right)=exp\left(\beta {\left[{\mathit{v}}_2\right]}_{\times }\right)=\mathbf{I}+sin\left(\beta \right){\left[{\mathit{v}}_2\right]}_{\times }+\left(1-cos\left(\beta \right)\right){\left[{\mathit{v}}_2\right]}_{\times }^2\end{array}$$

(54)

Therefore,

$$\begin{array}{c}\hfill {\mathit{b}}_2^T\left(\mathbf{I}+sin\left(\beta \right){\left[{\mathit{v}}_2\right]}_{\times }+\left(1-cos\left(\beta \right)\right){\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(\mathbf{I}+sin\left(\alpha \right){\left[{\mathit{v}}_1\right]}_{\times }+\left(1-cos\left(\alpha \right)\right){\left[{\mathit{v}}_1\right]}_{\times }^2\right){\mathbf{R}}_0{\mathit{a}}_2=0\end{array}$$

(55)

Clearly, Eq.(55) shows a one-dimensional curve constructed by $\alpha$ and $\beta$ (see Figure 4).

Specifically, given a $\alpha \in [-\pi ,\pi ]$, $\beta$ would be computed by Eq. (55). Intuitively, one should solve the following equation

$$sin\left(\beta \right)A+cos\left(\beta \right)B+C=0$$

(56)

where $A={\mathit{b}}_2^T{\left[{\mathit{v}}_2\right]}_{\times }{\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0{\mathit{a}}_2$, $B=-{\mathit{b}}_2^T{\left[{\mathit{v}}_2\right]}_{\times }^2{\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0{\mathit{a}}_2$ and $C={\mathit{b}}_2^T{\mathit{v}}_2{\mathit{v}}_2^T{\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0{\mathit{a}}_2$. According to the tangent half-angle formula, when $\beta \ne \pm \pi$,

$${\displaystyle \frac{2tan\left(\frac{\beta }{2}\right)}{1+{tan}^2\left(\frac{\beta }{2}\right)}}A+{\displaystyle \frac{1-{tan}^2\left(\frac{\beta }{2}\right)}{1+{tan}^2\left(\frac{\beta }{2}\right)}}B+C=0\stackrel{y\triangleq tan\left(\frac{\beta }{2}\right)}{\to }2Ay+(1-y^2)B+(1+y^2)C=0\Rightarrow (C-B)y^2+2Ay+B+C=0$$

(57)

Accordingly, $\beta$ can be computed and rotation $\mathbf{R}\left(\alpha ,\beta \right)$ can be computed as well. Formally,

$$f:\alpha \in [\pi ,\pi ]\stackrel{\beta \mathrm{by}\mathrm{Eq}.()}{\to }\mathbf{R}\left(\alpha ,\beta \right)$$

(58)

5.3. Intersections of Three Different Clifford Tori

Given three rotation constraints ${\mathit{b}}_1^T\mathbf{R}{\mathit{a}}_1=0$, ${\mathit{b}}_2^T\mathbf{R}{\mathit{a}}_2=0$ and ${\mathit{b}}_3^T\mathbf{R}{\mathit{a}}_3=0$, the possible rotation can be considered in the intersections of three Clifford tori. Note that it belongs to the famous E3Q3 problem [47], and there are many solutions have been proposed to solve it [48,49,50,51,52]. Here this paper provides a more straight linear way. Specifically,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\hfill {\mathit{b}}_1^T\mathbf{R}{\mathit{a}}_1=0\\ \hfill {\mathit{b}}_2^T\mathbf{R}{\mathit{a}}_2=0\\ \hfill {\mathit{b}}_3^T\mathbf{R}{\mathit{a}}_3=0\end{array}\right.\Rightarrow \left\{\begin{array}{cc}\hfill \mathbf{R}\left(\alpha ,\beta \right)& ={\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0\hfill \\ \hfill {\mathit{b}}_2^T\mathbf{R}{\mathit{a}}_2& =0\hfill \\ \hfill {\mathit{b}}_3^T\mathbf{R}{\mathit{a}}_3& =0\hfill \end{array}\right.\Rightarrow \left\{\begin{array}{c}\hfill {\mathit{b}}_2^T{\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0{\mathit{a}}_2=0\\ \hfill {\mathit{b}}_3^T{\mathbf{R}}_{v_2}\left(\beta \right){\mathbf{R}}_{v_1}\left(\alpha \right){\mathbf{R}}_0{\mathit{a}}_3=0\end{array}\right.\end{array}$$

(59)

Therefore,

$$\begin{array}{c}\hfill {\mathit{b}}_2^T\left(\mathbf{I}+sin\left(\beta \right){\left[{\mathit{v}}_2\right]}_{\times }+\left(1-cos\left(\beta \right)\right){\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(\mathbf{I}+sin\left(\alpha \right){\left[{\mathit{v}}_1\right]}_{\times }+\left(1-cos\left(\alpha \right)\right){\left[{\mathit{v}}_1\right]}_{\times }^2\right){\mathbf{R}}_0{\mathit{a}}_2=0\end{array}$$

(60)

$$\begin{array}{c}\hfill {\mathit{b}}_3^T\left(\mathbf{I}+sin\left(\beta \right){\left[{\mathit{v}}_2\right]}_{\times }+\left(1-cos\left(\beta \right)\right){\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(\mathbf{I}+sin\left(\alpha \right){\left[{\mathit{v}}_1\right]}_{\times }+\left(1-cos\left(\alpha \right)\right){\left[{\mathit{v}}_1\right]}_{\times }^2\right){\mathbf{R}}_0{\mathit{a}}_3=0\end{array}$$

(61)

There are two unknown variables $\alpha$ and $\beta$ for two constraints, and the intersections can be computed.

According to tangent half-angle formula (a.k.a. Weierstrass substitution formulas),

$$sin\left(\alpha \right)={\displaystyle \frac{2tan\left(\frac{\alpha }{2}\right)}{1+{tan}^2\left(\frac{\alpha }{2}\right)}},cos\left(\alpha \right)={\displaystyle \frac{1-{tan}^2\left(\frac{\alpha }{2}\right)}{1+{tan}^2\left(\frac{\alpha }{2}\right)}},sin\left(\beta \right)={\displaystyle \frac{2tan\left(\frac{\beta }{2}\right)}{1+{tan}^2\left(\frac{\beta }{2}\right)}},cos\left(\beta \right)={\displaystyle \frac{1-{tan}^2\left(\frac{\beta }{2}\right)}{1+{tan}^2\left(\frac{\beta }{2}\right)}}$$

(62)

Let $x=tan\left(\frac{\alpha }{2}\right)$, and $y=tan\left(\frac{\beta }{2}\right)$, and when $\alpha \ne \pm \pi$, $\beta \ne \pm \pi$,

$$\begin{array}{c}\hfill {\mathit{b}}_2^T\left((1+y^2)\mathbf{I}+2y{\left[{\mathit{v}}_2\right]}_{\times }+2y^2{\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left((1+x^2)\mathbf{I}+2x{\left[{\mathit{v}}_1\right]}_{\times }+2x^2{\left[{\mathit{v}}_1\right]}_{\times }^2\right){\mathbf{R}}_0{\mathit{a}}_2=0\end{array}$$

(63)

$$\begin{array}{c}\hfill {\mathit{b}}_3^T\left((1+y^2)\mathbf{I}+2y{\left[{\mathit{v}}_2\right]}_{\times }+2y^2{\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left((1+x^2)\mathbf{I}+2x{\left[{\mathit{v}}_1\right]}_{\times }+2x^2{\left[{\mathit{v}}_1\right]}_{\times }^2\right){\mathbf{R}}_0{\mathit{a}}_3=0\end{array}$$

(64)

The system of equations can be reformulated as

$$\mathcal{F}(x,y)=\left[\begin{array}{ccc}{y}^2& y& 1\end{array}\right]\mathbf{A}\left[\begin{array}{c}{x}^2\\ x\\ 1\end{array}\right]=0,\mathcal{G}(x,y)=\left[\begin{array}{ccc}{y}^2& y& 1\end{array}\right]\mathbf{B}\left[\begin{array}{c}{x}^2\\ x\\ 1\end{array}\right]=0$$

(65)

where $\mathbf{A}$ and $\mathbf{B}\in {\mathbb{R}}^{3\times 3}$, Specifically,

$$\mathbf{A}=\left[\begin{array}{ccc}{A}_{1,1}& A_{1,2}& A_{1,3}\\ A_{2,1}& A_{2,2}& A_{2,3}\\ A_{3,1}& A_{3,2}& A_{3,3}\end{array}\right]={\mathit{b}}_2^T*\begin{array}{ccc}\hfill \left[l\right]\left(\mathbf{I}+2{\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(\mathbf{I}+2{\left[{\mathit{v}}_1\right]}_{\times }^2\right)& \left(\mathbf{I}+2{\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(2{\left[{\mathit{v}}_1\right]}_{\times }\right)& \mathbf{I}+2{\left[{\mathit{v}}_2\right]}_{\times }^2\hfill \\ \hfill \left(2{\left[{\mathit{v}}_2\right]}_{\times }\right)\left(\mathbf{I}+2{\left[{\mathit{v}}_1\right]}_{\times }^2\right)& \left(2{\left[{\mathit{v}}_2\right]}_{\times }\right)\left(2{\left[{\mathit{v}}_1\right]}_{\times }\right)& 2{\left[{\mathit{v}}_2\right]}_{\times }\hfill \\ \hfill \mathbf{I}+2{\left[{\mathit{v}}_1\right]}_{\times }^2& 2{\left[{\mathit{v}}_1\right]}_{\times }& \mathbf{I}\hfill \end{array}*{\mathbf{R}}_0{\mathit{a}}_2$$

(66)

$$\mathbf{B}=\left[\begin{array}{ccc}{B}_{1,1}& B_{1,2}& B_{1,3}\\ B_{2,1}& B_{2,2}& B_{2,3}\\ B_{3,1}& B_{3,2}& B_{3,3}\end{array}\right]={\mathit{b}}_3^T*\begin{array}{ccc}\hfill \left[l\right]\left(\mathbf{I}+2{\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(\mathbf{I}+2{\left[{\mathit{v}}_1\right]}_{\times }^2\right)& \left(\mathbf{I}+2{\left[{\mathit{v}}_2\right]}_{\times }^2\right)\left(2{\left[{\mathit{v}}_1\right]}_{\times }\right)& \mathbf{I}+2{\left[{\mathit{v}}_2\right]}_{\times }^2\hfill \\ \hfill \left(2{\left[{\mathit{v}}_2\right]}_{\times }\right)\left(\mathbf{I}+2{\left[{\mathit{v}}_1\right]}_{\times }^2\right)& \left(2{\left[{\mathit{v}}_2\right]}_{\times }\right)\left(2{\left[{\mathit{v}}_1\right]}_{\times }\right)& 2{\left[{\mathit{v}}_2\right]}_{\times }\hfill \\ \hfill \mathbf{I}+2{\left[{\mathit{v}}_1\right]}_{\times }^2& 2{\left[{\mathit{v}}_1\right]}_{\times }& \mathbf{I}\hfill \end{array}*{\mathbf{R}}_0{\mathit{a}}_3$$

(67)

and * is to calculate the matrix item by item. Geometrically, if only considering real solutions, it describes the intersections of two cures computed by observations in $z=0$ plane (see Figure 5).

Algebraically, it is a polynomial system containing two equations in two variables, which can be generally solved by resultant-based method [53,54,55]. We donate

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}{y}^2& y& 1\end{array}\right]\mathbf{A}\left[\begin{array}{c}{x}^2\\ x\\ 1\end{array}\right]=\mathcal{A}\left(y\right)\left[\begin{array}{c}{x}^2\\ x\\ 1\end{array}\right]=0,\left[\begin{array}{ccc}{y}^2& y& 1\end{array}\right]\mathbf{B}\left[\begin{array}{c}{x}^2\\ x\\ 1\end{array}\right]=\mathcal{B}\left(y\right)\left[\begin{array}{c}{x}^2\\ x\\ 1\end{array}\right]=0\end{array}$$

(68)

Both sides of the equations are multiplied by x,

$$\mathcal{A}\left(y\right)\left[\begin{array}{c}{x}^3\\ x^2\\ x\end{array}\right]=0,\mathcal{B}\left(y\right)\left[\begin{array}{c}{x}^3\\ x^2\\ x\end{array}\right]=0$$

(69)

By stacking equations,

$$\left[\begin{array}{cc}0& \mathcal{A}\left(y\right)\\ 0& \mathcal{B}\left(y\right)\\ \mathcal{A}\left(y\right)& 0\\ \mathcal{B}\left(y\right)& 0\end{array}\right]\left[\begin{array}{c}{x}^3\\ x^2\\ x\\ 1\end{array}\right]=\mathcal{H}\left(y\right)\left[\begin{array}{c}{x}^3\\ x^2\\ x\\ 1\end{array}\right]=0$$

(70)

To solve y, it is to solve $\left|\begin{array}{c}\mathcal{H}\left(y\right)\end{array}\right|=0$. This determinant is also known as a hidden-variable resultant in computer vision filed [56,57], and it is an univariate polynomial with degree 8 of the hidden variable y.

$$\begin{array}{c}\hfill \left|\begin{array}{c}\mathcal{H}\left(y\right)\end{array}\right|=a_8{y}^8+a_7{y}^7+a_6{y}^6+a_5{y}^5+a_4{y}^4+a_3{y}^3+a_2{y}^2+a_{1}y+a_0=0\end{array}$$

(71)

In fact, the problem is not very complex, one can directly calculate the explicit expression of each term for the resultant-based method using symbolic math. Practically, one can obtain explicit coefficients ${\left\{a_i\right\}}_{i=0}^8$ by the function resultant in Matlab [55] or the function sympy.polys.polytools.resultant in SymPy1, and this paper will not list all details to avoid long-and-tedious expressions for readability. Nonetheless, the solution of y can be directly computed by solving roots of a univariate polynomial [58,59]. After solving y, it is easy to calculate x.

5.3.1. Linear Solution to Solve Intersections of Three Different Clifford Tori

Alternatively, Eq. (70) can be rewritten as a quadratic eigenvalue problem, i.e., $\left({\lambda }^2{\mathbf{H}}_2+\lambda {\mathbf{H}}_1+{\mathbf{H}}_0\right)\mathit{v}=0$ [56,57,60].

$$\mathcal{H}\left(y\right)\left[\begin{array}{c}{x}^3\\ x^2\\ x\\ 1\end{array}\right]=0\Rightarrow \left(y^2{\mathbf{H}}_2+y{\mathbf{H}}_1+{\mathbf{H}}_0\right)\left[\begin{array}{c}{x}^3\\ x^2\\ x\\ 1\end{array}\right]=0$$

(72)

where ${\mathbf{H}}_0$, ${\mathbf{H}}_1$ and ${\mathbf{H}}_2\in {\mathbb{R}}^{4\times 4}$. Specifically,

$${\mathbf{H}}_2=\left[\begin{array}{cccc}0& A_{1,1}& A_{1,2}& A_{1,3}\\ 0& B_{1,1}& B_{1,2}& B_{1,3}\\ A_{1,1}& A_{1,2}& A_{1,3}& 0\\ B_{1,1}& B_{1,2}& B_{1,3}& 0\end{array}\right],{\mathbf{H}}_1=\left[\begin{array}{cccc}0& A_{2,1}& A_{2,2}& A_{2,3}\\ 0& B_{2,1}& B_{2,2}& B_{2,3}\\ A_{2,1}& A_{2,2}& A_{2,3}& 0\\ B_{2,1}& B_{2,2}& B_{2,3}& 0\end{array}\right],{\mathbf{H}}_0=\left[\begin{array}{cccc}0& A_{3,1}& A_{3,2}& A_{3,3}\\ 0& B_{3,1}& B_{3,2}& B_{3,3}\\ A_{3,1}& A_{3,2}& A_{3,3}& 0\\ B_{3,1}& B_{3,2}& B_{3,3}& 0\end{array}\right]$$

(73)

Notably, there is no x and y in ${\mathbf{H}}_2$, ${\mathbf{H}}_1$ and ${\mathbf{H}}_0$. Accordingly, it is easy to solve the system of polynomial equations by linear way [60].

Specifically, the to-be-solved eigenvector of this problem is $\mathit{v}={\left[\begin{array}{cccc}{x}^3& x^2& x& 1\end{array}\right]}^T$ and corresponding eigenvalue is $\lambda =y$. Furthermore, the quadratic eigenvalue problem can be transformed into a generalized eigenvalue problem, i.e., $\mathbf{A}x=\lambda \mathbf{B}x$ [61], which has been extensively studied [57,61,62,63]. Specifically,

$$\left({\lambda }^2{\mathbf{H}}_2+\lambda {\mathbf{H}}_1+{\mathbf{H}}_0\right)\mathit{v}=0\Rightarrow \left[\begin{array}{cc}\mathbf{0}& \mathbf{I}\\ -{\mathbf{H}}_0& -{\mathbf{H}}_1\end{array}\right]\left[\begin{array}{c}\mathit{v}\\ \lambda \mathit{v}\end{array}\right]=\lambda \left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\mathbf{H}}_2\end{array}\right]\left[\begin{array}{c}\mathit{v}\\ \lambda \mathit{v}\end{array}\right]$$

(74)

Therefore, there will be 8 eigenvalues of this linear problem. The easiest way to solve this problem is using the Matlab function built-in polyeig2. SciPy also provides a function scipy.linalg.eig3 to solve generalized eigenvalue problems.