A Regularization-Patching Dual Quaternion Optimization Method for Solving the Hand-Eye Calibration Problem

Abstract

The hand-eye calibration problem is an important application problem in robot research. Based on the 2-norm of dual quaternion vectors, we propose a new dual quaternion optimization method for the hand-eye calibration problem. The dual quaternion optimization problem is decomposed to two quaternion optimization subproblems. The first quaternion optimization subproblem governs the rotation of the robot hand. It can be solved efficiently by the eigenvalue decomposition or singular value decomposition. If the optimal value of the first quaternion optimization subproblem is zero, then the system is rotationwise noiseless, i.e., there exists a “perfect” robot hand motion which meets all the testing poses rotationwise exactly. In this case, we apply the regularization technique for solving the second subproblem to minimize the distance of the translation. Otherwise we apply the patching technique to solve the second quaternion optimization subproblem. Then solving the second quaternion optimization subproblem turns out to be solving a quadratically constrained quadratic program. In this way, we give a complete description for the solution set of hand-eye calibration problems. This is new in the hand-eye calibration literature. The numerical results are also presented to show the efficiency of the proposed method.

1 Introduction

The hand-eye calibration problem is an important part of robot calibration, which has wide applications in aerospace, medical, automotive and industrial fields [10, 16]. The problem is to determine the homogeneous matrix between the robot gripper and a camera mounted rigidly on the gripper or between a robot base and the world coordinate system. In 1989, Shiu and Ahmad [32] and Tsai and Lenz [33] used one motion (two poses) to formulate the hand-eye calibration problem as solving a matrix equation

$$\begin{aligned} A X = X B, \end{aligned}$$

(1)

where X is the unknown homogeneous transformation matrix from the gripper (hand) to the camera (eye), A is the measurable homogeneous transformation matrix of the robot hand from its first to second position and B is the measurable homogeneous transformation matrix of the attached camera and also, from its first to second position. To allow the simultaneous estimation of both the transformations from the robot base frame to the world frame and from the robot hand frame to sensor frame, Zhuang et al. [41] derived another homogeneous transformation equation

$$\begin{aligned} A X = Z B, \end{aligned}$$

(2)

where X and Z are unknown homogeneous transformation matrices from the gripper to the camera and from the robot base to the world coordinate system, respectively, A is the transformation matrix from the robot base to the gripper and B is the transformation matrix from the world base to the camera. It is worth mentioning that there are other kinds of mathematical models for hand-eye calibration problem. In this paper, we focus on the models (1) and (2).

Over the years, many different methods and solutions are developed for the hand-eye calibration problem. Based on how the rotation and translation parameters are estimated, these approaches are broadly divided into two categories: separable solutions and simultaneous solutions. The separable solutions arise from solving the orientational component separately from the positional component. By using rotation matrix and translation vector to represent homogeneous transformation matrices, the hand-eye calibration equation is decomposed into rotation equation and position equation. The rotation parameters are first estimated. After that, the translation vectors could be estimated by solving a linear system. The different techniques that focus on the parametrization of rotation matrices include angle axis [32, 33, 35], Lie algebra [25], quaternions [3, 4, 14], Kronecker product [21, 31] and so on. The main drawback in these methods is that rotation estimation errors propagate to position estimation errors.

On the other hand, the simultaneous solutions arise from simultaneously solving the orientational component and the positional component. The rotation and translation parameters are solved either analytically or by means of numerical optimization. For analytical approaches, many techniques were proposed including quaternions [22], screw motion [2], Kronecker product [1], dual tensor [5], dual Lie algebra [6] and so on. The approaches based on numerical optimization include Levenberg–Marquardt algorithm [28, 42], gradient/Newton method [11], linear matrix inequality [12], alternative linear programming [40] and so on. For more details about solution methods for hand-eye calibration problem, one can refer to [10, 30] and references therein.

Among the solution methods for hand-eye calibration problem, the technique of dual quaternions was used to represent rigid motions by Daniilidis and Bayro-Corrochano [8]. Based on the dual quaternion parameterization, a simultaneous solution for the hand-eye problem was proposed by using the singular value decomposition [7, 8]. Compared with homogeneous transformation matrices, dual quaternions are considered to be a more efficient and compact way of representing the rotation and translation changes in a rigid body. Since the number of variables is reduced and the constraints are simple, the methods based on dual quaternions are computationally more efficient than the methods based on homogeneous transformation matrices [39]. In fact, it has been shown that the dual quaternion representation gives an efficient and robust way to estimate the solution of hand-eye calibration problem, as shown in [18,19,20, 23, 29, 34]. This also motivates us to use dual quaternions to formulate the hand-eye calibration problem as an optimization problem.

As far as we know, most existing methods for solving hand-eye calibration problems could handle the case when the rotation axes are not parallel. However, these methods may fail for the special case when the rotation axes are parallel, which is called the problem of pose singularity. For these methods, there is no description of the solution when the rotation axes are parallel, and it is difficult to judge whether the derived solutions are reliable or not.

In this paper, we propose a new dual quaternion optimization method for the hand-eye calibration problem based on the 2-norm of dual quaternion vectors, aiming to give a complete description of the solution set of the hand-eye calibration problem.

The theoretical base of dual quaternion optimization was established in [27], where a total order for dual numbers, the magnitude of a dual quaternion and the norm for dual quaternion vectors were proposed. Then, a two-stage solution scheme for equality constrained dual quaternion optimization problems was proposed in [26], with the hand-eye calibration problem and the simultaneous localization and mapping problem as application examples. It was shown in [26] that an equality constrained dual quaternion optimization problem could be solved by solving two quaternion optimization subproblems.

In the solution scheme of [26], the optimization solution set of the first quaternion optimization subproblem is designed as a constraint of the second quaternion optimization subproblem. This poses a challenge for implementing such a two-stage solution scheme in practice. In this paper, we propose a regularization-patching method to solve such a dual quaternion optimization problem arising from the hand-eye calibration problem. To apply the two-stage scheme of [26] to the hand-eye calibration problem, we may solve the first quaternion optimization subproblem efficiently by the eigenvalue decomposition or singular value decomposition. If the optimal value of the first subproblem is equal to zero, a regularization function is used to solve the second quaternion optimization subproblem. Otherwise, the solution of the second subproblem is determined by solving a patched quaternion optimization problem. In fact, the optimal value of the first subproblem is equal to zero if and only if there exists a “perfect” robot hand motion which meets all the testing poses exactly. In this case, we say that the hand-eye calibration system is rotationwise noiseless. The flowchart of proposed method is presented in Fig. 1. In this way, we give a complete description for the solution set of the hand-eye calibration problem. This is new in the hand-eye calibration literature and should be useful in applications.

Fig. 1

Flowchart of proposed method

In the next section, we present some preliminary knowledge on dual numbers, quaternions and dual quaternions. Based on dual quaternion optimization, the reformulations and analysis for hand-eye calibration equations $AX=XB$ and $AX=ZB$ are given in Sects. 3 and 4, respectively. In Sect. 5, we present the numerical results to show the efficiency of proposed methods. Some final remarks are made in Sect. 6.

Throughout the paper, the sets of real numbers, dual numbers, quaternion numbers and dual quaternion numbers are denoted by ${\mathbb {R}}$, ${\mathbb {D}}$, ${\mathbb {Q}}$ and $\mathbb{D}\mathbb{Q}$, respectively. The sets of n-dimensional real vectors, quaternion vectors and dual quaternion vectors are denoted by ${{\mathbb {R}}}^n$, ${{\mathbb {Q}}}^n$ and ${\mathbb{D}\mathbb{Q}}^n$, respectively. Scalars, vectors and matrices are denoted by lowercase letters, bold lowercase letters and capital letters, respectively.

2 Preliminaries

2.1 Dual Numbers

A dual number $q \in {\mathbb {D}}$ can be written as $q = q_{st} + q_\mathcal {I}\epsilon$, where $q_{st}, q_\mathcal {I}\in {\mathbb {R}}$ and $\epsilon$ is the infinitesimal unit satisfying $\epsilon ^2 = 0$. We call $q_{st}$ the standard part of q, and $q_\mathcal {I}$ the infinitesimal part of q. Dual numbers can be added in terms of components and multiplied by the formula

$$\begin{aligned} (p_{st} + p_\mathcal {I}\epsilon ) (q_{st} + q_\mathcal {I}\epsilon ) = p_{st} q_{st} +(p_{st} q_\mathcal {I}+ p_\mathcal {I}q_{st})\epsilon . \end{aligned}$$

The dual numbers form a commutative algebra of dimension two over the reals. The absolute value of $q = q_{st} + q_{\mathcal {I}} \epsilon \in {\mathbb {D}}$ is defined as

$$\begin{aligned} |q| = \left\{ \begin{aligned} |q_{st}| + \frac{ q_{st} q_{\mathcal {I}} }{|q_{st}|}\epsilon ,&\quad \textrm{if}\ q_{st} \not = 0, \\ |q_{\mathcal {I}}|\epsilon ,&\quad \mathrm{otherwise.} \end{aligned}\right. \end{aligned}$$

A total order “$\le$” for dual numbers was introduced in [27]. Given two dual numbers $p, q \in {\mathbb {D}}$, $p = p_{st} + p_\mathcal {I}\epsilon$, $q = q_{st} + q_\mathcal {I}\epsilon$, where $p_{st}, p_\mathcal {I}, q_{st}, q_\mathcal {I}\in {\mathbb {R}}$, we say that $p \le q$, if either $p_{st} < q_{st}$, or $p_{st} = q_{st}$ and $p_\mathcal {I}\le q_\mathcal {I}$. In particular, we say that p is positive, nonnegative, non-positive or negative, if $p > 0$, $p \ge 0$, $p \le 0$ or $p < 0$, respectively.

2.2 Quaternion Numbers

A quaternion number $q \in {\mathbb {Q}}$ has the form $q = q_0 + q_1 \textbf{i}+ q_2 \textbf{j}+ q_3 \textbf{k},$ where $q_0, q_1, q_2, q_3 \in {\mathbb {R}}$ and $\textbf{i}, \textbf{j}, \textbf{k}$ are three imaginary units of quaternions satisfying

$$\begin{aligned} \textbf{i}^2 = \textbf{j}^2 = \textbf{k}^2 = \textbf{i}\textbf{j}\textbf{k}=-1, \quad \textbf{i}\textbf{j}= -\textbf{j}\textbf{i}= \textbf{k}, \quad \textbf{j}\textbf{k}=-\textbf{k}\textbf{j}=\textbf{i}, \quad \textbf{k}\textbf{i}=-\textbf{i}\textbf{k}=\textbf{j}. \end{aligned}$$

The conjugate of q is the quaternion $q^*=q_0 - q_1 \textbf{i}- q_2 \textbf{j}- q_3 \textbf{k}$. The scalar part of q is $\text {Sc}(q) = \frac{1}{2}(q+q^*)=q_0$. Clearly, $\text {Sc}(q^*) = \text {Sc}(q)$ and $(p q)^* = q^* p^*$ for any $p, q \in {\mathbb {Q}}$. The multiplication of quaternions is associative and distributive over vector addition, but is not commutative. The magnitude of q is

$$\begin{aligned} |q| = \sqrt{q q^*} = \sqrt{q^* q}=\sqrt{q_0^2 + q_1^2 +q_2^2 + q_3^2}. \end{aligned}$$

The quaternion $q \in {\mathbb {Q}}$ is called a unit quaternion if $|q|=1$. It is well known [36] that the unit quaternion

$$\begin{aligned} q = \cos \left( \frac{\theta }{2}\right) + \sin \left( \frac{\theta }{2}\right) n_1 \textbf{i}+ \sin \left( \frac{\theta }{2}\right) n_2 \textbf{j}+ \sin \left( \frac{\theta }{2}\right) n_3 \textbf{k}, \end{aligned}$$

can be used to described the rotation around a unit axis $\textbf{n} = (n_1, n_2, n_3)^\top \in {\mathbb {R}}^3$ with an angle of $-\pi \le \theta \le \pi$. On the other hand, given a unit quaternion $q = q_0 + q_1 \textbf{i}+ q_2 \textbf{j}+ q_3 \textbf{k}\in {\mathbb {Q}}$, the rotation matrix R can be obtained by

$$\begin{aligned} R= \left( \begin{array}{ccc} q_0^2 + q_1^2-q_2^2-q_3^2 &{} 2(q_1 q_2 -q_0 q_3) &{} 2(q_1 q_3 + q_0 q_2) \\ 2(q_1 q_2 +q_0 q_3) &{} q_0^2 - q_1^2+q_2^2-q_3^2 &{} 2(q_2 q_3 - q_0 q_1) \\ 2(q_1 q_3 -q_0 q_2) &{} 2(q_2 q_3 + q_0 q_1) &{} q_0^2 - q_1^2 - q_2^2 + q_3^2 \end{array} \right) . \end{aligned}$$

(3)

For any $a = a_0 + a_1 \textbf{i}+ a_2 \textbf{j}+ a_3 \textbf{k}\in {\mathbb {Q}}$, denote $\overrightarrow{a}= (a_0, a_1, a_2, a_3)^\top \in {\mathbb {R}}^4$ and

$$\begin{aligned} M(a) = \left( \begin{array}{rrrr} a_0 &{} -a_1 &{} -a_2 &{} -a_3 \\ a_1 &{} a_0 &{} -a_3 &{} a_2 \\ a_2 &{} a_3 &{} a_0 &{} -a_1 \\ a_3 &{} -a_2 &{} a_1 &{} a_0 \end{array} \right) , \quad W(a) = \left( \begin{array}{rrrr} a_0 &{} -a_1 &{} -a_2 &{} -a_3 \\ a_1 &{} a_0 &{} a_3 &{} -a_2 \\ a_2 &{} -a_3 &{} a_0 &{} a_1 \\ a_3 &{} a_2 &{} -a_1 &{} a_0 \end{array} \right) . \end{aligned}$$

Clearly, $|a| = \Vert \overrightarrow{a} \Vert _2$. By direct calculations, we have the following propositions.

Proposition 2.1

For any $a = a_0 + a_1 \textbf{i}+ a_2 \textbf{j}+ a_3 \textbf{k}\in {\mathbb {Q}}$ and $b = b_0 + b_1 \textbf{i}+ b_2 \textbf{j}+ b_3 \textbf{k}\in {\mathbb {Q}}$, the following statements hold:

• (i) $\textrm{Sc}(r _1 a + r_2 b) = r_1 \textrm{Sc}(a) + r_2 \textrm{Sc}(b)$ for any $r_1, r_2 \in {\mathbb {R}}$.

• (ii) $\textrm{Sc}(a^* b) = \textrm{Sc}(a b^*) =\textrm{Sc}(b^* a) = \textrm{Sc}(b a^*) = \overrightarrow{a}^\top \overrightarrow{b}$.

• (iii) $M(a^*) = M(a)^\top$, $W(a^*) = W(a)^\top$.

• (iv) $\overrightarrow{ab} = M(a) \overrightarrow{b} = W(b) \overrightarrow{a}$.

• (v) $M(a)^\top M(a) = W(a)^\top W(a) = \Vert \overrightarrow{a} \Vert _2^2 I_{4\times 4}$, where $I_{4\times 4}$ is the identity matrix of size $4 \times 4$.

Proposition 2.2

If a and b are two quaternion numbers satisfying $\text {Sc}(a^* b) = 0$, then for any $q \in {\mathbb {Q}}$, we have $\text {Sc} (q^* a^* b q) = \text {Sc} (q^* b^* a q) = 0.$

Proof

Since $\text {Sc}(a^* b) = 0$, we have $a^* b + b^*a = 0$. According to Proposition 2.1, one can obtain that $\text {Sc} (q^* a^* b q) = \text {Sc} (q^* b^* a q) = \text {Sc} \left( q^* (\frac{a^* b + b^* a}{2}) q \right) = 0$ for any $q \in {\mathbb {Q}}.$ $\square$

Next we introduce the 2-norm for quaternion vectors, which can be found in [27]. Denote $\textbf{x}= (x_1, x_2, \cdots , x_n)^\top \in {{\mathbb {Q}}}^n$ for quaternion vectors. The 2-norm of $\textbf{x}\in {{\mathbb {Q}}}^n$ is defined as

$$\begin{aligned} \Vert \textbf{x}\Vert _2 = \sqrt{\sum _{i=1}^n |x_i|^2 } = \sqrt{\sum _{i=1}^n \Vert \overrightarrow{x_i} \Vert _2^2 }. \end{aligned}$$

The conjugate transpose of $\textbf{x}$ is defined as $\textbf{x}^* = (x_1^*, x_2^*, \cdots , x_n^*)$. More details about quaternions and quaternion vectors could be found in [37].

2.3 Dual Quaternion Numbers

A dual quaternion number $q \in \mathbb{D}\mathbb{Q}$ has the form $q = q_{st} + q_\mathcal {I}\epsilon$, where $q_{st}, q_\mathcal {I}\in {\mathbb {Q}}$. The conjugate of $q = q_{st} + q_\mathcal {I}\epsilon$ is $q^* = q_{st}^* + q_\mathcal {I}^* \epsilon$. The magnitude of a dual quaternion number $q = q_{st} + q_\mathcal {I}\epsilon$ is defined as

$$\begin{aligned} |q| = \left\{ \begin{aligned} |q_{st}| + \frac{\textrm{Sc}( q_{st}^* q_{\mathcal {I}})}{|q_{st}|}\epsilon ,&\quad \textrm{if}\ q_{st} \not = 0, \\ |q_{\mathcal {I}}|\epsilon ,&\quad \mathrm{otherwise.} \end{aligned}\right. \end{aligned}$$

The dual quaternion number $q \in \mathbb{D}\mathbb{Q}$ is called a unit dual quaternion if $|q|=1$. Note that $q = q_{st} + q_{\mathcal {I}} \epsilon \in \mathbb{D}\mathbb{Q}$ is a unit dual quaternion if and only if $q_{st}^* q_{st}= 1$ and $q_{st}^* q_{\mathcal {I}} + q_{\mathcal {I}}^* q_{st} = 0$. According to Proposition 2.2, we have the following result.

Corollary 2.3

If $q = q_{st} + q_{\mathcal {I}} \epsilon \in \mathbb{D}\mathbb{Q}$ is a unit dual quaternion, then $\text {Sc}(q_{st}^* q_{\mathcal {I}}) = \text {Sc}(q_{\mathcal {I}}^* q_{st}) = 0$, and for any $a \in {\mathbb {Q}}$, we have $\text {Sc}( a^* q_{st}^* q_{\mathcal {I}} a ) = \text {Sc}( a^* q_{\mathcal {I}}^* q_{st} a ) = 0.$

It has been shown that the 3D motion of a rigid body can be represented by a unit dual quaternion [7]. Consider a rigid motion in SE(3) represented by a $4 \times 4$ homogeneous transformation matrix

$$\begin{aligned} T = \left( \begin{array}{cc} R &{} \textbf{t} \\ \textbf{0}^\top &{} 1 \end{array}\right) , \end{aligned}$$

(4)

where $R \in {\mathbb {R}}^{3 \times 3}$ is the rotation matrix about an axis through the origin and $\textbf{t} \in {\mathbb {R}}^3$ is the translation vector. Let $q_{st} \in {\mathbb {Q}}$ be the unit quaternion encoding the rotation matrix R, and let $t \in {\mathbb {Q}}$ be the quaternion satisfying $\overrightarrow{t} = \left( \begin{array}{c} 0 \\ \textbf{t} \end{array} \right)$. Then the transformation matrix T is represented by the dual quaternion ${q = q_{st} + q_\mathcal {I}\epsilon },$ where $q_\mathcal {I}= \frac{1}{2} t q_{st}$. It is not difficult to check that $q \in \mathbb{D}\mathbb{Q}$ is a unit dual quaternion since

$$\begin{aligned} \textrm{Sc}(q_{st}^* q_\mathcal {I}) = \frac{1}{2} \textrm{Sc}(q_{st}^* t q_{st}) =0. \end{aligned}$$

On the other hand, given a unit dual quaternion $q = q_{st} + q_\mathcal {I}\epsilon \in \mathbb{D}\mathbb{Q}$, the corresponding homogeneous transformation matrix T can be obtained by (4), where the rotation matrix $R \in {\mathbb {R}}^{3 \times 3}$ can be derived from the unit quaternion $q_{st}$ according to (3) and the translation vector $\textbf{t} \in {\mathbb {R}}^{3 }$ can be derived from

$$\begin{aligned} \left( \begin{array}{c} 0 \\ \textbf{t} \end{array} \right) = \overrightarrow{2 q_\mathcal {I}q_{st}^*}. \end{aligned}$$

(5)

It follows that $\Vert \textbf{t}\Vert _2^2 = 4 q_\mathcal {I}q_{st}^* q_{st} q_\mathcal {I}^* = 4|q_\mathcal {I}|^2$. In other words, for a unit dual quaternion, the magnitude of its infinitesimal part is half of the length of the corresponding translation vector.

Denote $\textbf{x}= (x_1, x_2, \cdots , x_n)^\top \in {\mathbb{D}\mathbb{Q}}^n$ for dual quaternion vectors. We may also write

$$\begin{aligned} \textbf{x}= \textbf{x}_{st} + \textbf{x}_\mathcal {I}\epsilon , \end{aligned}$$

where $\textbf{x}_{st}, \textbf{x}_\mathcal {I}\in {{\mathbb {Q}}}^n$. The 2-norm of $\textbf{x}\in {\mathbb{D}\mathbb{Q}}^n$ is defined as

$$\begin{aligned} \Vert \textbf{x}\Vert _2 = \left\{ \begin{aligned}\sqrt{\sum _{i=1}^n |x_i|^2},&\quad \textrm{if}\ \textbf{x}_{st} \not = \textbf{0}, \\ \Vert \textbf{x}_\mathcal {I}\Vert _2\epsilon ,&\quad \mathrm{otherwise.} \end{aligned}\right. \end{aligned}$$

(6)

Denote by $\textbf{x}^*:= (x_1^*, x_2^*, \cdots , x_n^*)$ the conjugate transpose of $\textbf{x}\in \mathbb{D}\mathbb{Q}^n$. According to Proposition 6.3 of [27], it holds that

$$\begin{aligned} \Vert \textbf{x}\Vert _2 = \Vert \textbf{x}_{st} \Vert _2 + \frac{ \text {Sc}({\textbf{x}}_{st}^* {\textbf{x}_{\mathcal {I}}})}{ \Vert {\textbf{x}}_{st} \Vert _2} \epsilon , \end{aligned}$$

(7)

for any $\textbf{x}\in \mathbb{D}\mathbb{Q}^n$ with $\textbf{x}_{st} \ne \textbf{0}$.

3 Hand-Eye Calibration Equation $AX=XB$

The hand-eye calibration problem is to find the matrix X such that

$$\begin{aligned} A^{(i)} X =X B^{(i)} \end{aligned}$$

(8)

for $i=1,2,\ldots , n$, where X is transformation matrix from the gripper (hand) to the camera (eye), $A^{(i)}$ is the transformation matrix between the grippers of two different poses and $B^{(i)}$ the transformation matrix between the cameras of two different poses. The transformation matrices X, $A^{(i)}$ and $B^{(i)}$ are encoded with the unit dual quaternions

$$\begin{aligned} x=x_{st} + x_{\mathcal {I}} \epsilon , \quad {a}^{(i)}={a}^{(i)}_{st} + a^{(i)}_{\mathcal {I}} \epsilon , \quad {b}^{(i)} = b^{(i)}_{st} + b^{(i)}_{\mathcal {I}} \epsilon , \end{aligned}$$

for $i=1,2,\ldots ,n$. Let $\textbf{a} = \left( a^{(1)}, a^{(2)}, \ldots , a^{(n)} \right) ^\top \in \mathbb{D}\mathbb{Q}^n$ and $\textbf{b} = \big (b^{(1)}, b^{(2)}, \ldots , b^{(n)} \big )^\top \in \mathbb{D}\mathbb{Q}^n$. The hand-eye calibration problem (8) can be reformulated as the dual quaternion optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \Vert \textbf{a} x- x \textbf{b} \Vert _2 \\ \text {s.t.}&{} |x| = 1, \quad x \in \mathbb{D}\mathbb{Q}. \end{array} \end{aligned}$$

(9)

In the following, we assume that the optimal solution set of (9) is non-empty.

Denote $\textbf{f}(x)=\textbf{a} x- x \textbf{b} \in \mathbb{D}\mathbb{Q}^n$. According to (6) and (7), we have

$$\begin{aligned} \Vert \textbf{f}(x) \Vert _2 = \left\{ \begin{array}{ll} \Vert \textbf{f}_{st}(x) \Vert _2 + \frac{ \text {Sc}\left( \textbf{f}_{st}^*(x) \textbf{f}_{\mathcal {I}}(x) \right) }{ \Vert \textbf{f}_{st}(x) \Vert _2} \epsilon , &{} \text {if } \textbf{f}_{st}(x) \ne \textbf{0}, \\ \Vert \textbf{f}_{\mathcal {I}}(x) \Vert _2 \epsilon , &{} \text {otherwise}. \end{array} \right. \end{aligned}$$

Problem (9) can be divided to two different cases, which need to be handled very differently. One case is that the standard part of the optimal value of (9) is zero. Another case is that the standard part of the optimal value of (9) is positive. Physically, the standard part of the optimal value of (9) is zero if and only if there exists a “perfect” robot hand motion x, which meets all the n testing poses rotationwise exactly. In this case, we say that system is rotationwise noiseless. The following proposition provides a way to check whether the system is rotationwise noiseless or not.

Proposition 3.1

If ${\hat{x}}$ is an optimal solution of (9), the standard part ${\hat{x}}_{st}$ is an optimal solution of the quaternion optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \Vert \textbf{a}_{st} x_{st}- x_{st} \textbf{b}_{st} \Vert _2^2 \\ \mathrm{s.t.}&{} |x_{st}| = 1, \quad x_{st} \in {\mathbb {Q}}. \end{array} \end{aligned}$$

(10)

Hence, the system is rotationwise noiseless if and only if the optimal value of (10) is equal to zero.

Proof

According to the definition of total order for dual numbers, the result could be easily proved since $\textbf{f}_{st}(x) = \textbf{a}_{st} x_{st}- x_{st} \textbf{b}_{st}$. $\square$

Denote the optimal solution set of (10) by $X_{st}$. If the optimal value of (10) is equal to zero, we consider the regularized quaternion optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \Vert \textbf{f}_{\mathcal {I}} (x) \Vert _2^2 + \gamma (x_{st}^* x_{st} + x_{\mathcal {I}}^* x_\mathcal {I}) \\ \text {s.t.}&{} x_{st} \in X_{st}, \quad x_{st}^*x_{\mathcal {I}} + x_{\mathcal {I}}^*x_{st} = 0, \quad {x_{st} \in {\mathbb {Q}}, \quad x_{\mathcal {I}} \in {\mathbb {Q}}}, \end{array} \end{aligned}$$

(11)

where $\gamma$ is the parameter that balances the loss function and the regularization term. In fact, $x_{st} \in X_{st}$ implies $x_{st}^* x_{st} =1$, and $x_{\mathcal {I}}^* x_\mathcal {I}$ is proportional to the norm square of translation vector. By adding the regularization term, we try to find the best solution with minimal distance of translation. This explains the role of regularization.

If the optimal value of (10) is not equal to zero, we consider the quaternion optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \text {Sc}\left( \textbf{f}_{st}^*(x) \textbf{f}_{\mathcal {I}}(x) \right) \\ \text {s.t.}&{} x_{st} \in X_{st}, \quad x_{st}^*x_{\mathcal {I}} + x_{\mathcal {I}}^*x_{st} = 0, \quad {x_{st} \in {\mathbb {Q}}, \quad x_{\mathcal {I}} \in {\mathbb {Q}}}. \end{array} \end{aligned}$$

(12)

By using the matrix representation for quaternion numbers, problems (10), (11) and (12) could be solved efficiently. For $i=1,2,\ldots , n$, we have

$$\begin{aligned} \overrightarrow{ a^{(i)}_{st} x_{st}-x_{st}b_{st}^{(i)} }= \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] \overrightarrow{ x_{st}} \end{aligned}$$

and

$$\begin{aligned} \left| a^{(i)}_{st} x_{st}-x_{st}b_{st}^{(i)} \right| ^2 = \overrightarrow{ x_{st}}^\top \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] ^\top \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] \overrightarrow{ x_{st}}. \end{aligned}$$

Denote

$$\begin{aligned} L_{11} = \sum _{i=1}^n \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] ^\top \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] . \end{aligned}$$

(13)

It follows that

$$\begin{aligned} \Vert \textbf{a}_{st} x_{st}- x_{st} \textbf{b}_{st} \Vert _2^2 = \sum _{i=1}^n \left| a^{(i)}_{st} x_{st}-x_{st}b_{st}^{(i)} \right| ^2 = \overrightarrow{x_{st}}^\top L_{11} \overrightarrow{x_{st}}. \end{aligned}$$

Denote the minimal eigenvalue of matrix $L_{11}$ by $\lambda _0$. As a result, problem (10) is equivalent to finding the unit eigenvectors corresponding to $\lambda _0$.

Similarly, for $i=1,2,\ldots , n$, we have

$$\begin{aligned} \overrightarrow{ a_{st}^{(i)}x_{\mathcal {I}} +a_{\mathcal {I}}^{(i)}x_{st} - x_{st}b_{\mathcal {I}}^{(i)} - x_{\mathcal {I}}b_{st}^{(i)} }= \left[ M\left( a^{(i)}_{st}\right) -W \left( b_{st}^{(i)} \right) \right] \overrightarrow{ x_{\mathcal {I}}} + \left[ M\left( a^{(i)}_{\mathcal {I}}\right) -W\left( b_{\mathcal {I}}^{(i)}\right) \right] \overrightarrow{ x_{st}}. \end{aligned}$$

Denote

$$\begin{aligned} L_{22} = \sum _{i=1}^n \left[ M \left( a^{(i)}_{\mathcal {I}}\right) -W \left( b_{\mathcal {I}}^{(i)}\right) \right] ^\top \left[ M\left( a^{(i)}_{\mathcal {I}}\right) -W\left( b_{\mathcal {I}}^{(i)}\right) \right] \end{aligned}$$

(14)

and

$$\begin{aligned} L_{12} = \sum _{i=1}^n \left[ M \left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] ^\top \left[ M\left( a^{(i)}_{\mathcal {I}}\right) -W\left( b_{\mathcal {I}}^{(i)}\right) \right] . \end{aligned}$$

(15)

It follows that

$$\begin{aligned}{} & {} \Vert \textbf{f}_{\mathcal {I}} (x) \Vert _2^2 = \sum _{i=1}^n \left| a_{st}^{(i)}x_{\mathcal {I}} +a_{\mathcal {I}}^{(i)}x_{st} - x_{st}b_{\mathcal {I}}^{(i)} - x_{\mathcal {I}}b_{st}^{(i)} \right| ^2 = \overrightarrow{x_{\mathcal {I}}}^\top L_{11} \overrightarrow{x_{\mathcal {I}}} \\{} & {} \quad + 2 \overrightarrow{x_{\mathcal {I}}}^\top L_{12} \overrightarrow{x_{st}} + \overrightarrow{x_{st}}^\top L_{22} \overrightarrow{x_{st}}. \end{aligned}$$

As a result, problem (11) is equivalent to the optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \overrightarrow{x_{\mathcal {I}}}^\top (L_{11} + \gamma I) \overrightarrow{x_{\mathcal {I}}} + 2 \overrightarrow{x_{\mathcal {I}}}^\top L_{12} \overrightarrow{x_{st}} + \overrightarrow{x_{st}}^\top (L_{22} + \gamma I) \overrightarrow{x_{st}} \\ \text {s.t.}&{} \overrightarrow{x_{st}} \in \overrightarrow{X_{st}}, \quad \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_{st}}\in {\mathbb {R}}^4, \quad \overrightarrow{x_{\mathcal {I}}} \in {\mathbb {R}}^4}, \end{array} \end{aligned}$$

(16)

where $\overrightarrow{X_{st}}$ is the set of all the unit eigenvectors corresponding to the minimal eigenvalue of matrix $L_{11}$. Once the set $\overrightarrow{X_{st}}$ is determined, problem (16) turns out to be a quadratically constrained quadratic program (QCQP).

To be specific, suppose that the dimension of the eigenspace of the minimal eigenvalue of $L_{11}$ is k. Let $Q \in {\mathbb {R}}^{4 \times k}$ be the matrix whose columns form an orthonormal basis of the eigenspace, i.e., $Q^\top Q = I_{k \times k}$. It is not difficult to see that $\overrightarrow{X_{st}} = \{ Q \textbf{y}: \textbf{y}^\top \textbf{y}=1, \textbf{y}\in {\mathbb {R}}^k \}$. Problem (16) can be rewritten as

$$\begin{aligned} \begin{array}{cl} \min &{} \overrightarrow{x_{\mathcal {I}}}^\top (L_{11} + \gamma I) \overrightarrow{x_{\mathcal {I}}} + 2 \overrightarrow{x_{\mathcal {I}}}^\top L_{12} Q \textbf{y}+ \textbf{y}^\top Q^\top (L_{22} + \gamma I) Q \textbf{y}\\ \text {s.t.}&{} {\textbf{y}}^\top \textbf{y}=1, \quad \textbf{y}^\top Q^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad {\textbf{y}\in {\mathbb {R}}^k, \quad \overrightarrow{ x_{\mathcal {I}}} \in {\mathbb {R}}^4}. \end{array} \end{aligned}$$

(17)

In particular, if the dimension of the eigenspace is one, i.e., $k=1$, the solution set $\overrightarrow{X_{st}} = \{ \textbf{q}, -\textbf{q}\}$, where $\textbf{q}\in {\mathbb {R}}^4$ is the normalized basis of the eigenspace. In this case, problem (17) could be solved efficiently by representing $\overrightarrow{ x_{\mathcal {I}}}$ in the orthogonal complement space of $\textbf{q}$.

In the following, we reformulate problem (12) as an optimization problem by using the matrix representation for quaternion numbers. According to Proposition 2.1, we have

$$\begin{aligned} \begin{aligned}&\text {Sc} \left( \left( a^{(i)}_{st} x_{st}-x_{st}b_{st}^{(i)} \right) ^* \left( a_{st}^{(i)}x_{\mathcal {I}} +a_{\mathcal {I}}^{(i)}x_{st} - x_{st}b_{\mathcal {I}}^{(i)} - x_{\mathcal {I}}b_{st}^{(i)} \right) \right) \\&= \overrightarrow{ x_{st}}^\top \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] ^\top \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] \overrightarrow{ x_{\mathcal {I}}} \\&\hspace{4cm} + \overrightarrow{ x_{st}}^\top \left[ M\left( a^{(i)}_{st}\right) -W\left( b_{st}^{(i)}\right) \right] ^\top \left[ M\left( a^{(i)}_{\mathcal {I}}\right) -W\left( b_{\mathcal {I}}^{(i)}\right) \right] \overrightarrow{ x_{st}} \end{aligned} \end{aligned}$$

for $i=1,2,\ldots , n$. It follows that

$$\begin{aligned} \begin{aligned} \text {Sc} \left( \textbf{f}_{st}^*(x) \textbf{f}_{\mathcal {I}}(x) \right)&= \sum _{i=1}^n \text {Sc} \left( \left( a^{(i)}_{st} x_{st}-x_{st}b_{st}^{(i)} \right) ^* \left( a_{st}^{(i)}x_{\mathcal {I}} +a_{\mathcal {I}}^{(i)}x_{st} - x_{st}b_{\mathcal {I}}^{(i)} - x_{\mathcal {I}}b_{st}^{(i)} \right) \right) \\&= \overrightarrow{ x_{st}}^\top L_{11} \overrightarrow{ x_{\mathcal {I}}} + \overrightarrow{ x_{st}}^\top L_{12} \overrightarrow{ x_{st}}, \end{aligned} \end{aligned}$$

where $L_{11}$ and $L_{12}$ are given by (13) and (15), respectively. Note that $\overrightarrow{X_{st}}$ is the set of all unit eigenvectors corresponding to the minimal eigenvalue $\lambda _0$ of $L_{11}$. Under the constraints of (12), one can obtain that

$$\begin{aligned} \overrightarrow{ x_{st}}^\top L_{11} \overrightarrow{ x_{\mathcal {I}}} = \lambda _0 \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} =0, \end{aligned}$$

since $L_{11}$ is symmetric. It turns out that problem (12) is equivalent to the optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \overrightarrow{ x_{st}}^\top L_{12} \overrightarrow{ x_{st}} \\ \text {s.t.}&{} \overrightarrow{x_{st}} \in \overrightarrow{X_{st}}, \quad \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_{st}}\in {\mathbb {R}}^4, \quad \overrightarrow{x_{\mathcal {I}}} \in {\mathbb {R}}^4}. \end{array} \end{aligned}$$

(18)

Similarly, if Q is the matrix whose columns form an orthonormal basis of the eigenspace of $\lambda _0$, the optimal $\overrightarrow{x_{st}}$ can be derived by computing the unit eigenvectors corresponding to the minimal eigenvalue of $\text {Sym}\left( Q^\top L_{12} Q\right) = \left( Q^\top L_{12} Q + Q^\top L_{12}^\top Q \right) /2$. Since the objective function in (18) does not contain $\overrightarrow{x_\mathcal {I}}$, the optimal $\overrightarrow{x_\mathcal {I}}$ can be any vector which is orthogonal to the optimal $\overrightarrow{x_{st}}$. We may need to find a proper one via sewing a patch on the optimal set of $\overrightarrow{x_\mathcal {I}}$ once the optimal $\overrightarrow{x_{st}}$ is determined. Considering the continuity of the norm, it is naturally necessary to further search for $x_\mathcal {I}$ under the constrains of $\overrightarrow{x_{st}}^\top \overrightarrow{x_{\mathcal {I}}} = 0$, such that $\Vert \textbf{f}_\mathcal {I}(x)\Vert _2$ is reduced as much as possible, i.e.,

$$\begin{aligned} \begin{array}{cl} \min _{\overrightarrow{x_\mathcal {I}}} &{} \overrightarrow{x_{\mathcal {I}}}^\top L_{11} \overrightarrow{x_{\mathcal {I}}} + 2 \overrightarrow{x_{\mathcal {I}}}^\top L_{12} \overrightarrow{x_{st}} + \overrightarrow{x_{st}}^\top L_{22} \overrightarrow{x_{st}} \\ \text {s.t.}&{} \overrightarrow{x_{st}}^\top \overrightarrow{x_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_{\mathcal {I}}} \in {\mathbb {R}}^4}. \end{array} \end{aligned}$$

(19)

This explains the role of the patching.

Note that in this way, we give a complete description for the solution set of the hand-eye calibration problem. This is new in the hand-eye calibration literature and should be useful in applications.

To conclude, the solution method for hand-eye calibration equation $AX=XB$ is summarized in Algorithm 1.

Algorithm 1

Dual quaternion optimization for $AX = XB$

4 Hand-Eye Calibration Equation $AX=ZB$

In 1994, Zhuang et al. [41] generalized (1) to $AX=ZB$, where X is transformation matrix from the gripper to the camera, Z is the transformation matrix from the robot base to the world coordinate system, A is the transformation matrix from the robot base to the gripper and B is the transformation matrix from the world base to the camera. Given n measurements $\left( A^{(i)}, B^{(i)}\right) _{i=1}^n$, the problem is to find the best solution X and Z such that

$$\begin{aligned} A^{(i)} X =Z B^{(i)} \end{aligned}$$

(20)

for $i=1,2,\ldots , n$. The transformation matrices X, Z, $A^{(i)}$ and $B^{(i)}$ are encoded with the unit dual quaternions

$$\begin{aligned} x=x_{st} + x_{\mathcal {I}} \epsilon , \quad z=z_{st} + z_{\mathcal {I}} \epsilon , \quad {a}^{(i)}={a}^{(i)}_{st} + a^{(i)}_{\mathcal {I}} \epsilon , \quad {b}^{(i)} = b^{(i)}_{st} + b^{(i)}_{\mathcal {I}} \epsilon , \end{aligned}$$

for $i=1,2,\ldots ,n$. Let $\textbf{a} = \left( a^{(1)}, a^{(2)}, \ldots , a^{(n)} \right) ^\top \in \mathbb{D}\mathbb{Q}^n$ and $\textbf{b} = \big (b^{(1)}, b^{(2)}, \ldots , b^{(n)} \big )^\top \in \mathbb{D}\mathbb{Q}^n$. The hand-eye calibration problem (20) can be reformulated as the dual quaternion optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \Vert \textbf{a} x- z \textbf{b} \Vert _2 \\ \text {s.t.}&{} |x| = |z| = 1, \quad {x \in \mathbb{D}\mathbb{Q}, \quad z \in \mathbb{D}\mathbb{Q}}. \end{array} \end{aligned}$$

(21)

Similarly, we say that the system is rotationwise noiseless if and only if the standard part of the optimal value of (21) is zero.

Denote $\textbf{g}(x, z)=\textbf{a} x- z \textbf{b} \in \mathbb{D}\mathbb{Q}^n$. To solve problem (21), according to the definition of 2-norm for dual quaternion vectors, we first consider the quaternion optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \Vert \textbf{g}_{st}(x, z) \Vert _2^2 = \Vert \textbf{a}_{st} x_{st}- z_{st} \textbf{b}_{st} \Vert _2^2 \\ \text {s.t.}&{} |x_{st}| = |z_{st}| = 1, \quad {x_{st} \in {\mathbb {Q}}, \quad z_{st} \in {\mathbb {Q}}}. \end{array} \end{aligned}$$

(22)

Note that $a = a_{st} + a_{\mathcal {I}} \epsilon \in \mathbb{D}\mathbb{Q}$ is a unit dual quaternion if and only if $a_{st}^* a_{st}= 1$ and $\text {Sc}\left( a_{st}^* a_{\mathcal {I}} \right) = a_{st}^* a_{\mathcal {I}} + a_{\mathcal {I}}^* a_{st} = 0$. For $i=1,2,\ldots , n$, we have

$$\begin{aligned} \begin{aligned} \left| a^{(i)}_{st} x_{st}-z_{st}b_{st}^{(i)} \right| ^2&= \left( a^{(i)}_{st} x_{st}-z_{st}b_{st}^{(i)} \right) ^* \left( a^{(i)}_{st} x_{st}-z_{st}b_{st}^{(i)} \right) \\&= 2 - 2 \text {Sc} \left( x_{st}^* \left( a_{st}^{(i)}\right) ^* z_{st}b_{st}^{(i)} \right) \\&= 2 - 2 \overrightarrow{x_{st}}^\top M\left( a_{st}^{(i)}\right) ^\top W \left( b_{st}^{(i)}\right) \overrightarrow{z_{st}} \end{aligned} \end{aligned}$$

since x, z, $a^{(i)}$ and $b^{(i)}$ are unit dual quaternions. Denote

$$\begin{aligned} K_{11} = \sum _{i=1}^n M\left( a_{st}^{(i)}\right) ^\top W\left( b_{st}^{(i)}\right) . \end{aligned}$$

(23)

It follows that

$$\begin{aligned} \Vert \textbf{a}_{st} x_{st}- z_{st} \textbf{b}_{st} \Vert _2^2 = \sum _{i=1}^n \left| a^{(i)}_{st} x_{st}-z_{st}b_{st}^{(i)} \right| ^2 = 2n -2 \overrightarrow{x_{st}}^\top K_{11} \overrightarrow{z_{st}}. \end{aligned}$$

Then problem (22) is equivalent to the optimization problem

$$\begin{aligned} \begin{array}{cl} \max &{} \overrightarrow{x_{st}}^\top K_{11} \overrightarrow{z_{st}} \\ \text {s.t.}&{} \overrightarrow{x_{st}}^\top \overrightarrow{x_{st}} = \overrightarrow{z_{st}}^\top \overrightarrow{z_{st}} = 1, \quad {\overrightarrow{x_{st}} \in {\mathbb {R}}^4, \quad \overrightarrow{z_{st}} \in {\mathbb {R}}^4}. \end{array} \end{aligned}$$

(24)

Denote the maximal singular value of $K_{11}$ by $\sigma _1$, the set of optimal vector pairs of (24) by $\overrightarrow{\Omega _{st}}$. As a result, problem (22) aims at finding the unit singular vector pairs for $\sigma _1$, which can be solved efficiently by the singular value decomposition (SVD).

If the optimal value of (22) is equal to zero, i.e., $\sigma _1 = n$, consider the regularized optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \Vert \textbf{g}_{\mathcal {I}}(x, z) \Vert _2^2 + \gamma \left( {\Vert \overrightarrow{x_{st}}\Vert _2^2 + \Vert \overrightarrow{x_{\mathcal {I}}}\Vert _2^2 + \Vert \overrightarrow{z_{st}}\Vert _2^2 + \Vert \overrightarrow{z_{\mathcal {I}}}\Vert _2^2 }\right) \\ \text {s.t.}&{} \left( \overrightarrow{x_{st}}, \overrightarrow{z_{st}} \right) \in \overrightarrow{\Omega _{st}}, \quad \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad \overrightarrow{ z_{st}}^\top \overrightarrow{ z_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_{st}}, \overrightarrow{x_\mathcal {I}} \in {\mathbb {R}}^4, \quad \overrightarrow{z_{st}}, \overrightarrow{z_{\mathcal {I}}} \in {\mathbb {R}}^4}, \end{array}\nonumber \\ \end{aligned}$$

(25)

where $\gamma$ is the regularization parameter and

$$\begin{aligned} \Vert \textbf{g}_{\mathcal {I}}(x, z) \Vert _2^2 = \sum _{i=1}^n \left\| M\left( a_{st}^{(i)}\right) \overrightarrow{x_{\mathcal {I}}} + M\left( a_{\mathcal {I}}^{(i)}\right) \overrightarrow{x_{st}} - W\left( b_{\mathcal {I}}^{(i)}\right) \overrightarrow{z_{st}} - W\left( b_{st}^{(i)}\right) \overrightarrow{z_{\mathcal {I}}} \right\| _2^2. \end{aligned}$$

Once the set $\overrightarrow{\Omega _{st}}$ is determined, problem (25) could be also written as an QCQP. To be specific, suppose the singular value decomposition of matrix $K_{11}$ is $K_{11}=U \Sigma V^\top$, where $U, V \in {\mathbb {R}}^{4 \times 4}$ are orthogonal and $\Sigma \in {\mathbb {R}}^{4 \times 4}$ is diagonal. Let $Q_1 \in {\mathbb {R}}^{4 \times k}$ be the matrix whose columns are the columns of U corresponding to $\sigma _1$, and let $Q_2 \in {\mathbb {R}}^{4 \times k}$ be the matrix whose columns are the columns of V corresponding to $\sigma _1$. It is not difficult to see that $\overrightarrow{\Omega _{st}} = \left\{ \left( Q_1 \textbf{y}, Q_2 \textbf{y}\right) : \textbf{y}^\top \textbf{y}=1, \textbf{y}\in {\mathbb {R}}^k \right\}$. In fact, for any unit vectors $\textbf{y}_1$ and $\textbf{y}_2$, the value of objective function of (24) at the point $\left( \overrightarrow{x_{st}}, \overrightarrow{z_{st}} \right) = \left( Q_1 \textbf{y}_1, Q_2 \textbf{y}_2 \right)$ is

$$\begin{aligned} \overrightarrow{x_{st}}^\top K_{11} \overrightarrow{z_{st}} = \textbf{y}_1^\top Q_1^\top K_{11} Q_2 \textbf{y}_2 = \sigma _1 \textbf{y}_1^\top \textbf{y}_2 \le \sigma _1, \end{aligned}$$

according to the Cauchy–Schwarz inequality. Without loss of generality, we assume $\sigma _1 > 0$. Then the equality holds if and only if $\textbf{y}_1 = \textbf{y}_2$. As a result, problem (25) can be rewritten as an QCQP:

$$\begin{aligned}{} & {} \min \displaystyle \sum _{i=1}^n \left\| M\left( a_{st}^{(i)}\right) \overrightarrow{x_{\mathcal {I}}} + M\left( a_{\mathcal {I}}^{(i)}\right) Q_1 \textbf{y}- W\left( b_{\mathcal {I}}^{(i)}\right) Q_2 \textbf{y}- W\left( b_{st}^{(i)}\right) \overrightarrow{z_{\mathcal {I}}} \right\| _2^2 \nonumber \\{} & {} \qquad \quad \qquad + \gamma \left( { \Vert \overrightarrow{ x_{\mathcal {I}}}\Vert _2^2 + \Vert \overrightarrow{ z_{\mathcal {I}}}\Vert _2^2 + 2 \Vert \textbf{y}\Vert _2^2} \right) \nonumber \\{} & {} \text {s.t.}\, \textbf{y}^\top \textbf{y}= 1, \quad \textbf{y}^\top Q_1^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad \textbf{y}^\top Q_2^\top \overrightarrow{ z_{\mathcal {I}}} = 0, \quad {\textbf{y}\in {\mathbb {R}}^k, \quad \overrightarrow{x_{\mathcal {I}}} \in {\mathbb {R}}^4, \quad \overrightarrow{z_{\mathcal {I}}} \in {\mathbb {R}}^4}.\nonumber \\ \end{aligned}$$

(26)

In particular, when $k=1$, problem (26) could be solved efficiently by representing $\overrightarrow{ x_{\mathcal {I}}}$ and $\overrightarrow{ z_{\mathcal {I}}}$ in the corresponding orthogonal complement space of $Q_1$ and $Q_2$, respectively.

On the other hand, if the optimal value of (22) is not equal to zero, consider the optimization problem

$$\begin{aligned} \begin{array}{cl} \min &{} \text {Sc}\left( \textbf{g}_{st}^*(x, z) \textbf{g}_{\mathcal {I}}(x, z) \right) \\ \text {s.t.}&{} \left( \overrightarrow{x_{st}}, \overrightarrow{z_{st}} \right) \in \overrightarrow{\Omega _{st}}, \quad \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad \overrightarrow{ z_{st}}^\top \overrightarrow{ z_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_{st}}, \overrightarrow{x_\mathcal {I}} \in {\mathbb {R}}^4, \quad \overrightarrow{z_{st}}, \overrightarrow{z_{\mathcal {I}}} \in {\mathbb {R}}^4}. \end{array}\nonumber \\ \end{aligned}$$

(27)

According to Corollary 2.3, we have

$$\begin{aligned}{} & {} \text {Sc} \left( \left( a^{(i)}_{st} x_{st}\right) ^* a_{st}^{(i)}x_{\mathcal {I}} \right) = \text {Sc} \left( \left( a^{(i)}_{st} x_{st}\right) ^* a_{\mathcal {I}}^{(i)}x_{st} \right) = \text {Sc} \left( \left( z_{st} b_{st}^{(i)}\right) ^* z_{st}b_{\mathcal {I}}^{(i)} \right) \\{} & {} \quad =\text {Sc} \left( \left( z_{st} b_{st}^{(i)}\right) ^* z_{\mathcal {I}}b_{st}^{(i)} \right) = 0 \end{aligned}$$

since x, z, $a^{(i)}$ and $b^{(i)}$ are unit quaternions for $i=1,2,\ldots , n$. It follows that

$$\begin{aligned} \begin{aligned}&\text {Sc} \left( \left( a^{(i)}_{st} x_{st}-z_{st}b_{st}^{(i)} \right) ^* \left( a_{st}^{(i)}x_{\mathcal {I}} +a_{\mathcal {I}}^{(i)}x_{st} - z_{st}b_{\mathcal {I}}^{(i)} - z_{\mathcal {I}}b_{st}^{(i)} \right) \right) \\&= - \text {Sc} \left( \left( a^{(i)}_{st} x_{st}\right) ^* z_{st}b_{\mathcal {I}}^{(i)} + \left( a^{(i)}_{st} x_{st}\right) ^* z_{\mathcal {I}}b_{st}^{(i)} + \left( z_{st} b_{st}^{(i)}\right) ^* a_{st}^{(i)}x_{\mathcal {I}} + \left( z_{st} b_{st}^{(i)}\right) ^* a_{\mathcal {I}}^{(i)}x_{st} \right) \\&= - \overrightarrow{x_{st}}^\top \left[ M \left( a_{st}^{(i)}\right) ^\top W\left( b_{\mathcal {I}}^{(i)}\right) + M\left( a_{\mathcal {I}}^{(i)}\right) ^\top W\left( b_{st}^{(i)}\right) \right] \overrightarrow{z_{st}} \\&\quad - \overrightarrow{x_{st}}^\top M\left( a_{st}^{(i)}\right) ^\top W\left( b_{st}^{(i)}\right) \overrightarrow{z_{\mathcal {I}}} - \overrightarrow{x_{\mathcal {I}}}^\top M\left( a_{st}^{(i)}\right) ^\top W\left( b_{st}^{(i)}\right) \overrightarrow{z_{st}}. \end{aligned} \end{aligned}$$

Denote

$$\begin{aligned} K_{12} = \sum _{i=1}^n M\left( a_{st}^{(i)}\right) ^\top W\left( b_{\mathcal {I}}^{(i)}\right) \end{aligned}$$

(28)

and

$$\begin{aligned} K_{21} = \sum _{i=1}^n M\left( a_{\mathcal {I}}^{(i)}\right) ^\top W\left( b_{st}^{(i)}\right) . \end{aligned}$$

(29)

By simple computation, one can obtain that

$$\begin{aligned} \begin{aligned} \text {Sc}\left( \textbf{g}_{st}^*(x, z) \textbf{g}_{\mathcal {I}}(x, z) \right)&= \sum _{i=1}^n \text {Sc}\left( \left( a^{(i)}_{st} x_{st}-z_{st}b_{st}^{(i)} \right) ^* \left( a_{st}^{(i)}x_{\mathcal {I}} +a_{\mathcal {I}}^{(i)}x_{st} - z_{st}b_{\mathcal {I}}^{(i)} - z_{\mathcal {I}}b_{st}^{(i)} \right) \right) \\&= - \left[ \overrightarrow{ x_{st}}^\top (K_{12} + K_{21} ) \overrightarrow{ z_{st}} + \overrightarrow{ x_{st}}^\top K_{11} \overrightarrow{ z_{\mathcal {I}}} + \overrightarrow{ x_{\mathcal {I}}}^\top K_{11} \overrightarrow{ z_{st}} \right] , \end{aligned} \end{aligned}$$

where $K_{11}$, $K_{12}$ and $K_{21}$ are given by (23), (28) and (29), respectively. Under the constraints of problem (27), $\overrightarrow{x_{st}}$ and $\overrightarrow{z_{st}}$ are left-singular and right-singular vectors corresponding to the maximal singular value $\sigma _1$ for $K_{11}$, which means

$$\begin{aligned} K_{11} \overrightarrow{z_{st}} = \sigma _1 \overrightarrow{x_{st}} \quad \text {and} \quad K_{11}^\top \overrightarrow{x_{st}} = \sigma _1 \overrightarrow{z_{st}}. \end{aligned}$$

Then we have $\overrightarrow{ x_{st}}^\top K_{11} \overrightarrow{ z_{\mathcal {I}}} = \sigma _1 \overrightarrow{ z_{st}}^\top \overrightarrow{ z_{\mathcal {I}}} = 0$ and $\overrightarrow{ x_{\mathcal {I}}}^\top K_{11} \overrightarrow{ z_{st}} = \sigma _1 \overrightarrow{ x_{\mathcal {I}}}^\top \overrightarrow{ x_{st}} = 0$ under the constraints of problem (27). As a result, problem (27) is equivalent to the optimization

$$\begin{aligned} \begin{array}{cl} \max &{} \overrightarrow{ x_{st}}^\top (K_{12} + K_{21} ) \overrightarrow{ z_{st}} \\ \text {s.t.}&{} \left( \overrightarrow{x_{st}}, \overrightarrow{z_{st}} \right) \in \overrightarrow{\Omega _{st}}, \quad \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad \overrightarrow{ z_{st}}^\top \overrightarrow{ z_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_{st}}, \overrightarrow{x_\mathcal {I}} \in {\mathbb {R}}^4, \quad \overrightarrow{z_{st}}, \overrightarrow{z_{\mathcal {I}}} \in {\mathbb {R}}^4}. \end{array}\nonumber \\ \end{aligned}$$

(30)

Similarly, given the singular value decomposition $K_{11} = U \Sigma V^\top$, let $Q_1$ be the matrix whose columns are the columns of U corresponding to $\sigma _1$, and let $Q_2$ be the matrix whose columns are the columns of V corresponding to $\sigma _1$. The optimal $\overrightarrow{x_{st}}$ and $\overrightarrow{z_{st}}$ can be derived by computing the unit eigenvectors corresponding to the maximal eigenvalue of $\text {Sym}\left( Q_1^\top (K_{12}+K_{21}) Q_2 \right)$. Since the objective function in (30) does not contain $\overrightarrow{ x_{\mathcal {I}}}$ and $\overrightarrow{ z_{\mathcal {I}}}$, the optimal $\overrightarrow{x_{\mathcal {I}}}$ can be any vector which is orthogonal to the optimal $\overrightarrow{x_{st}}$, and the optimal $\overrightarrow{z_{\mathcal {I}}}$ can be any vector which is orthogonal to the optimal $\overrightarrow{z_{st}}$. Considering the continuity of the norm, once the optimal $\overrightarrow{x_{st}}$ and $\overrightarrow{z_{st}}$ are determined, we try to find the best one in the optimal set of $\overrightarrow{x_\mathcal {I}}$ and $\overrightarrow{z_\mathcal {I}}$ such that the patching function $\Vert \textbf{g}_{\mathcal {I}} (x, z) \Vert _2^2$ is minimized, i.e.,

$$\begin{aligned} \begin{array}{cl} \max _{\overrightarrow{x_{\mathcal {I}}}, \overrightarrow{z_{\mathcal {I}}}} &{} \displaystyle \sum _{i=1}^n \left\| M\left( a_{st}^{(i)}\right) \overrightarrow{x_{\mathcal {I}}} + M\left( a_{\mathcal {I}}^{(i)}\right) \overrightarrow{x_{st}} - W\left( b_{\mathcal {I}}^{(i)}\right) \overrightarrow{z_{st}} - W\left( b_{st}^{(i)}\right) \overrightarrow{z_{\mathcal {I}}} \right\| _2^2 \\ \text {s.t.}&{} \overrightarrow{ x_{st}}^\top \overrightarrow{ x_{\mathcal {I}}} = 0, \quad \overrightarrow{ z_{st}}^\top \overrightarrow{ z_{\mathcal {I}}} = 0, \quad {\overrightarrow{x_\mathcal {I}} \in {\mathbb {R}}^4, \quad \overrightarrow{z_{\mathcal {I}}} \in {\mathbb {R}}^4}. \end{array}\nonumber \\ \end{aligned}$$

(31)

To conclude, the solution method for hand-eye calibration equation $AX=ZB$ is summarized in Algorithm 2.

Algorithm 2

Dual quaternion optimization for $AX = ZB$

5 Numerical Experiments

In this section, we report a set of synthetic experiments to show the efficiency of proposed methods for hand-eye calibration problem. All the codes are written in MATLAB R2017a. The numerical experiments were done on a desktop with an Intel Core i5-2430 M CPU dual-core processor running at 2.4GHz and 6GB of RAM.

In the implementation of our proposed methods, we use GloptiPoly [13] to construct SDP relaxations of QCQPs, and call the interior point optimizer in MOSEK [24] to solve SDPs. Further, GloptiPoly can also recover the solution to the original problem and certify its optimality. We set the regularization parameter $\gamma = 2 \times 10^{-6}$. For hand-eye calibration model $AX=XB$, we compare our method with the direct estimation proposed by Tsai et al. [33] (denoted by “Tsai89”), the Kronecker method proposed by Andreff et al. [1] (denoted by “Andreff99”), the classic dual quaternion method proposed by Daniilidis [7] (denoted by “Daniilidis99”), the improved dual quaternion method proposed by Malti et al. [23] (denoted by “Malti10”) and the dual quaternion method using polynomial optimization proposed by Heller et al. [12] (denoted by “Heller14”).

For hand-eye calibration model $AX=ZB$, we compare our method with the quaternion method proposed by Zhuang et al. [41] (denoted by “Zhuang94”), the quaternion method proposed by Dornaika et al. [9] (denoted by “Dornaika98”), the classic dual quaternion method proposed by Li et al. [18] (denoted by “Li10”), the dual quaternion method using polynomial optimization proposed by Heller et al. [12] (denoted by “Heller14”) and the dual quaternion method proposed by Li et al. [20] (denoted by “Li18”).

Numerical experiments are carried out as follows. First, the original homogeneous transformation matrices ${\hat{X}}$ and ${\hat{Z}}$ in (2) are given by

$$\begin{aligned} {\hat{X}} = \left( \begin{array}{rrrr} 0.9995 &{} -0.0100 &{} 0.0297 &{} 9.190 \\ 0.0116 &{} 0.9986 &{} -0.0523 &{} 5.397 \\ -0.0291 &{} 0.0526 &{} 0.9982 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1.0000 \\ \end{array} \right) , \end{aligned}$$

(32)

and

$$\begin{aligned} {\hat{Z}} = \left( \begin{array}{rrrr} 0.2790 &{} -0.0981 &{} -0.9553 &{} 164.226 \\ -0.5439 &{} 0.8037 &{} -0.2414 &{} 301.638 \\ 0.7914 &{} 0.5869 &{} 0.1709 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1.0000 \end{array} \right) . \end{aligned}$$

(33)

Second, we generate n transformation matrices $A^{(i)}$, $i=1,2,\ldots , n$. Then the transformation matrix $B^{(i)}$ is computed by $B^{(i)} = {\hat{Z}}^{-1} A^{(i)} {\hat{X}}$ for $i=1,2,\ldots , n$. We use the methods Zhuang94, Dornaika98, Li10, Heller14, Li18 and Algorithm 2 to solve the hand-eye calibration equation $AX=ZB$ with the given matrices $\left( A^{(i)}, B^{(i)}\right) _{i=1}^n$. For hand-eye calibration equation $AX=XB$, we construct $\frac{n(n-1)}{2}$ pairs of matrices $\left( \left( A^{(i)}\right) ^{-1} A^{(j)}, \left( B^{(i)}\right) ^{-1} B^{(j)} \right) _{i<j}$, denoted by $\left( \tilde{A}^{(s)}, \tilde{B}^{(s)}\right) _{s=1}^{n(n-1)/2}$. Then the methods Tsai89, Andreff99, Daniilidis99, Malti10, Heller14 and Algorithm 1 are used to solve the hand-eye calibration equation $AX=XB$ with the given matrices $\left( \tilde{A}^{(s)}, \tilde{B}^{(s)}\right) _{s=1}^{n(n-1)/2}$. The estimation errors are computed by

$$\begin{aligned} e_X = \Vert X- {\hat{X}} \Vert _2, \quad e_Z = \Vert Z- {\hat{Z}} \Vert _2. \end{aligned}$$

It is worth mentioning that the computation time for our algorithms consists of the time of constructing SDP relaxations, the time of solving SDPs and the time of checking the optimality of the solutions. This is the reason why the proposed algorithms may need more time to get the solution when compared with other direct methods, such as Tsai89, Andreff99, Daniilidis99 for hand-eye calibration model $AX=XB$ and Zhuang94, Dornaika98 and Li10 for hand-eye calibration model $AX=ZB$. Here we mainly focus on the quality of the solutions and the comparison of the computation time within the same type of methods. In fact, the proposed algorithms could be improved by using other efficient ways of solving the resulting QCQP, such as copositive relaxations or doubly nonnegative relaxations as in [15, 17].

5.1 Measurements with Non-parallel Rotation Axis

Four measurements of A with non-parallel rotation axis are given by

$$\begin{aligned} A^{(1)} = \left( \begin{array}{rrrr} 0.1752 &{} -0.6574 &{} 0.7329 &{} -10.5536 \\ 0.6325 &{} -0.4954 &{} -0.5954 &{} -30.5304 \\ 0.7545 &{} 0.5679 &{} 0.3290 &{} 50.4851 \\ 0 &{} 0 &{} 0 &{} 1.0000 \end{array} \right) , \\ A^{(2)} = \left( \begin{array}{rrrr} -0.0745 &{} 0.9661 &{} 0.2471 &{} -20.4123 \\ 0.8573 &{} -0.0645 &{} 0.5108 &{} -50.8904 \\ 0.5094 &{} 0.2499 &{} -0.8234 &{} 80.8685 \\ 0 &{} 0 &{} 0 &{} 1.0000 \end{array} \right) , \\ A^{(3)} = \left( \begin{array}{rrrr} -0.1456 &{} -0.6867 &{} 0.7122 &{} -20.5519 \\ 0.8252 &{} -0.4814 &{} -0.2955 &{} -30.6491 \\ 0.5458 &{} 0.5447 &{} 0.6367 &{} 60.4312 \\ 0 &{} 0 &{} 0 &{} 1.0000 \end{array} \right) , \\ A^{(4)} = \left( \begin{array}{rrrr} -0.1434 &{} -0.5250 &{} 0.8389 &{} -10.5892 \\ 0.8158 &{} -0.5427 &{} -0.2001 &{} -50.6730 \\ 0.5603 &{} 0.6557 &{} 0.5061 &{} 80.4641 \\ 0 &{} 0 &{} 0 &{} 1.0000 \end{array} \right) . \end{aligned}$$

As described above, we have four measurements $\left( A^{(i)}, B^{(i)}\right)$ for equation $AX=ZB$, and six motions $\left( \tilde{A}^{(s)}, \tilde{B}^{(s)}\right)$ for equation $AX=XB$. The numerical results for $AX=XB$ and $AX=ZB$ with non-parallel rotation axis are reported in Tables 1 and 2, respectively. Algorithms 1 and 2 show the best behavior in terms of estimation error. Note that the first three methods in Table 1 and the first three methods in Table 2 get the solution via solving linear equations, while the other methods need to call SDP solvers to get the solution. That explains why Algorithm 1 may need more computation time to get the solution when compared with the first three methods in Table 1, and Algorithm 2 may need more computation time to get the solution when compared with the first three methods in Table 2.

Table 1 Numerical results for $AX=XB$ with non-parallel rotation axis

Table 2 Numerical results for $AX=ZB$ with non-parallel rotation axis

5.2 Measurements with Parallel Rotation Axis

In this subsection, we test our algorithms for the case that all the axes of measurements are parallel, which is often the situation for the hand-eye calibration of SCARA robots [34]. In this case, it has been shown that the problem is not well defined and there exists a 1D manifold of equivalent solutions with identical algebraic error [2, 38]. To evaluate the quality of solutions, we try to find the solution such that the third component of its translation vector is equal to zero, and then compare it with the real solution ${\hat{X}}$ and ${\hat{Z}}$ given by (32) and (33), respectively.

Four measurements of A are generated with the same rotation axis, but with different angles. Without loss of generality, the normalized rotation axis is $\textbf{n} = ( 0, 0, 1)^\top$. For $A^{(1)}, A^{(2)}, A^{(3)}, A^{(4)}$, the rotation angles are $\theta _1 = \frac{\pi }{6}$, $\theta _2=\frac{\pi }{3}$, $\theta _3 = -\frac{\pi }{6}$ and $\theta _4 = -\frac{\pi }{3}$, while their translation vectors are randomly generated given by

$$\begin{aligned} \textbf{t}_1 = (-10.9865, 12.3788, -27.2571)^\top , \quad \textbf{t}_2 = (38.8986, 84.6736, -93.8814)^\top , \\ \textbf{t}_3 = (-75.7189, -53.6187, 28.5794)^\top , \quad \textbf{t}_4 = (-52.8133, 93.3732, -70.1666)^\top , \end{aligned}$$

respectively. The numerical results for $AX=XB$ and $AX=ZB$ with parallel rotation axis are reported in Tables 3 and 4, respectively. As we can see, the transformation matrices could be also estimated by the proposed algorithms with the minimal error. In this case, the linear equations used in Tsai89, Andreff99, Zhuang94 and Dornaika98 have infinite solutions, and the SVDs used in Malti10 and Li18 have infinite unit singular vectors corresponding to the minimal singular value. As expected, the particular solutions derived by these methods could result in larger estimation errors when compared with other methods.

Table 3 Numerical results for $AX=XB$ with parallel rotation axis

Table 4 Numerical results for $AX=ZB$ with parallel rotation axis

Fig. 2

Robustness testing for $AX=XB$

5.3 Measurement Estimation with Noise

In practice, the measurement of B is typically estimated using visual processing. Since visual estimation is noisy, this set of experiment aims comparing the robustness of the different methods to disturbances in the measurement of B.

The four measurements $\left( A^{(i)}, B^{(i)}\right) _{i=1}^4$ are the same as those in Sect. 5.1. The rotation and translation of $B^{(i)}$ are disturbed by adding zero mean Gaussian noise with increasing standard deviation. Note that the motions $\tilde{B}^{(s)}$ are also disturbed when adding noise to the measurements $B^{(i)}$. The standard deviation of the additive noise increases from 0 to 0.02 in steps of 0.002. For each standard deviation, the average errors of $e_X$ and $e_Z$ are recorded after 10 runs of each method. The robustness testing for $AX=XB$ and $AX=ZB$ with noisy measurements of B are plotted in Figs. 2 and 3, respectively. For $AX=XB$, Algorithm 1 may be not as robust as Daniilidis99 and Heller14. However, Algorithm 2 shows the competitive robustness to the noise when compared with other methods for $AX=ZB$.

Fig. 3

Robustness testing for $AX=ZB$

6 Final Remarks

In this paper, we establish a new dual quaternion optimization method for the hand-eye calibration problem based on the 2-norm of dual quaternion vectors. A two-stage method is also proposed by using the techniques of regularization and patching. However, there are still some problems that need further study. We have the following final remarks.

• 1. Can we use some other norms for dual quaternion vectors, e.g., 1-norm, $\infty$-norm, instead of 2-norm in this method?

• 2. We may also consider some other hand-eye calibration models, such as multi-camera hand-eye calibration.

• 3. How can we choose the regularization parameter $\gamma$ to improve the efficiency of the method?

• 4. Can we extend this method to the simultaneous localization and mapping problem?