A regularization-patching dual quaternion optimization method for solving the hand-eye calibration problem

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.


Introduction
The hand-eye calibration problem is an important part of robot calibration, which has wide applications in aerospace, medical, automotive and industrial fields [15,10].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 [29] and Tsai and Lenz [30] used one motion (two poses) to formulate the hand-eye calibration problem as solving a matrix equation 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, Roth and Sudhaker [38] derived another homogeneous transformation equation where X is the 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.It is worth mentioning that there are other kinds of mathematical models for hand-eye calibration problem.In this paper, we only 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 [29,30,32], Lie algebra [22], quaternions [3,4,14], Kronecker product [18,28] 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 [19], 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 [39,25], gradient/Newton method [11], linear matrix inequality [12], alternative linear programming [37], and pseudo-inverse [36].For more details about solution methods for hand-eye calibration problem, one can refer to [10,27] 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 [8,7].After that, many solution methods based on dual quaternions were proposed [26,16,20,31,17].It has been shown that the dual quaternion representation gives a stable way to estimate the solution.
The existing methods for the hand-eye calibration problem used to produce solutions in general cases i.e., the rotation axes are not parallel.There lacks a complete description of the solution set of the hand-eye calibration problem.
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 [24], 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 [23], with the hand-eye calibration problem and the simultaneous localization and mapping problem as application examples.It was shown in [23] that an equality constrained dual quaternion optimization problem could be solved by solving two quaternion optimization subproblems.
In the solution scheme of [23], 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 [23] 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 flow chart of proposed method is presented in Figure 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.
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 Sections 3 and 4, respectively.In Section 5, we present the numerical results to show the efficiency of proposed methods.Some final remarks are made in Section 6.
Throughout the paper, the sets of real numbers, dual numbers, quaternion numbers and dual quaternion numbers are denoted by R, D, Q and DQ, respectively.The sets of n-dimensional real vectors, quaternion vectors and dual quaternion vectors are denoted by R n , Q n and DQ n , respectively.Scalars, vectors and matrices are denoted by lowercase letters, bold lowercase letters and capital letters, respectively.

Dual numbers
A dual number q ∈ D can be written as q = q st + q I , where q st , q I ∈ R and is the infinitesimal unit satisfying 2 = 0. We call q st the standard part of q, and q I the infinitesimal part of q.Dual numbers can be added component-wise, and multiplied by the formula (p st + p I )(q st + q I ) = p st q st + (p st q I + p I q st ) .
The dual numbers form a commutative algebra of dimension two over the reals.The absolute value of q = q st + q I ∈ D is defined as A total order "≤" for dual numbers was introduced in [24].Given two dual numbers p, q ∈ D, p = p st + p I , q = q st + q I , where p st , p I , q st , q I ∈ R, we say that p ≤ q, if either p st < q st , or p st = q st and p I ≤ q I .In particular, we say that p is positive, nonnegative, nonpositive or negative, if p > 0, p ≥ 0, p ≤ 0 or p < 0, respectively.
Proposition 2.1.For any a = a 0 + a , where I 4×4 is the identity matrix of size 4 × 4.
Proposition 2.2.If a and b are two quaternion numbers satisfying Sc(a * b) = 0, then for any q ∈ Q, we have Sc(q * a * bq) = Sc(q * b * aq) = 0.
Next we introduce the 2-norm for quaternion vectors, which can be found in [24].Denote The conjugate transpose of x is defined as . More details about quaternions and quaternion vectors could be found in [34].

Dual quaternion numbers
A dual quaternion number q ∈ DQ has the form q = q st + q I , where q st , q I ∈ Q.The conjugate of q = q st + q I is q * = q * st + q * I .The magnitude of a dual quaternion number q = q st + q I is defined as The dual quaternion number q ∈ DQ is called a unit dual quaternion if |q| = 1.Note that q = q st + q I ∈ DQ is a unit dual quaternion if and only if q * st q st = 1 and q * st q I + q * I q st = 0.According to Proposition 2.2, we have the following result.
Corollary 2.3.If q = q st + q I ∈ DQ is a unit dual quaternion, then Sc(q * st q I ) = Sc(q * I q st ) = 0 and for any a ∈ Q, we have Sc(a * q * st q I a) = Sc(a * q * 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 × 4 homogeneous transformation matrix where R ∈ R 3×3 is the rotation matrix about an axis through the origin and t ∈ R 3 is the translation vector.Let q st ∈ Q be the unit quaternion encoding the rotation matrix R and let t ∈ Q be the quaternion satisfying . Then the transformation matrix T is represented by the dual quaternion q = q st + q I , where On the other hand, given a unit dual quaternion q = q st + q I ∈ DQ, the corresponding homogeneous transformation matrix T can be obtained by ( 4), where the rotation matrix R ∈ R 3×3 can be derived from the unit quaternion q st according to (3) and the translation vector t ∈ R 3 can be derived from 0 It follows that t 2 2 = 4q I q * st q st q * I = 4|q 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 for dual quaternion vectors.We may also write where x st , x I ∈ Q n .The 2-norm of x ∈ DQ n is defined as Denote by According to Proposition 6.3 of [24], it holds that for any x ∈ DQ n with x st = 0.
3 Hand-Eye Calibration Equation AX = XB The hand-eye calibration problem is to find the matrix X such that for i = 1, 2, . . ., 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 for i = 1, 2, . . ., n.Let a = a (1) , a (2) , . . ., a (n) ∈ DQ n and b = b (1) , b (2) , . . ., b Denote f (x) = ax − xb ∈ DQ n .According to ( 6) and ( 7), we have 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 if the system is rotationwise noiseless or not.
Proposition 3.1.If x is an optimal solution of (9), the standard part xst is an optimal solution of the quaternion optimization problem 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 Denote the optimal set of (10) by X st .If the optimal value of ( 10) is equal to zero, we consider the regularized quaternion optimization problem where γ is the parameter that balances the loss function and the regularization term.In fact, x st ∈ X st implies x * st x st = 1, and x * I x 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 min Sc (f By using the matrix representation for quaternion numbers, problems (10), ( 11) and ( 12) could be solved efficiently.For i = 1, 2, . . ., n, we have and a Denote It follows that Denote the minimal eigenvalue of matrix L 11 by λ 0 .As a result, problem ( 10) is equivalent to finding the unit eigenvectors corrresponding to λ 0 .Similarly, for i = 1, 2, . . ., n, we have Denote and It follows that As a result, problem (11) is equivalent to the optimization problem where − → X st is the set of all the unit eigenvectors corrresponding to the minimal eigenvalue of matrix L 11 .Once the set − → X st is determinated, 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 ∈ R 4×k be the matrix whose columns form an orthonormal basis of the eigenspace, i.e., Q Q = I k×k .It is not difficult to see that − → X st = {Qy : y y = 1, y ∈ R k }.Problem ( 16) can be rewritten as In particular, if the dimension of the eigenspace is one, i.e., k = 1, the solution set − → X st = {q, −q}, where q ∈ R 4 is the normalized basis of the eigenspace.In this case, problem (17) could be solved efficiently by representing − → x I in the orthogonal complement space of 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 Sc a where L 11 and L 12 are given by ( 13) and ( 15) respectively.Note that − → X st is the set of all unit eigenvectors corresponding to the minimal eigenvalue λ 0 of L 11 .Under the constraints of ( 12), one can obtain that − → Similarly, if Q is the matrix whose columns form an orthonormal basis of the eigenspace of λ 0 , the optimal − → x st can be derived by computing the unit eigenvectors corresponding to the minimal eigenvalue of Sym Q Since the objective function in (18) does not contain − → x I , the optimal − → x I can be any vector which is orthogonal to the optimal − → x st .
We may need to find a proper one via sewing a patch on the optimal set of − → x I , while ensuring that f st (x) 2 is minimized.Considering the continuity of the norm, it is naturally necessary to further search for x I under the constrains of − → x st − → x I = 0, such that f I (x) 2 is reduced as much as possible, i.e., min − → 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 Require: Motions A (i) , B (i) n i=1 , regularization parameter γ.Ensure: The hand-eye transformation matrix X.

6:
Compute x st by finding the unit eigenvector corresponding to the minimal eigenvalue of Sym Q L 12 Q .

7:
Compute x I by solving (19) with the optimal x st .
8: end if , where R is computed from x st by (3) and t is computed from x st and x I by using (5).

Hand-Eye Calibration Equation AX = ZB
In 1994, Zhuang, Roth and Sudhaker [38] 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 A (i) , B (i) n i=1 , the problem is to find the best solution X and Z such that for i = 1, 2, . . ., n.The transformation matrices X, Z, A (i) and B (i) are encoded with the unit dual quaternions for i = 1, 2, . . ., n.Let a = a (1) , a (2) , . . ., a (n) ∈ DQ n and b = b (1) , b (2) , . . ., b (n) ∈ DQ n .The hand-eye calibration problem (20) can be reformulated as the dual quaternion optimization problem min ax − zb 2 s.t.
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 g(x, z) = ax − zb ∈ DQ n .To solve problem (21), according to the definition of 2-norm for dual quaternion vectors, we first consider the quaternion optimization problem Note that a = a st + a I ∈ DQ is a unit dual quaternion if and only if a * st a st = 1 and Sc (a * st a I ) = a * st a I + a * I a st = 0.For i = 1, 2, . . ., n, we have since x, z, a (i) and b (i) are unit dual quaternions.Denote It follows that Then problem ( 22) is equivalent to the optimization problem Denote the maximal singular value of K 11 by σ 1 , the set of optimal vector pairs of (24) by − → Ω st .As a result, problem (22) is to find the unit singular vector pairs for σ 1 , which can be solved efficiently by the singular value decomposition (SVD).
If the optimal value of ( 22) is equal to zero, i.e., σ 1 = n, consider the regularized optimization problem min g where γ is the regularization parameter and .
Once the set − → Ω 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 ΣV , where U, V ∈ R 4×4 are orthogonal and Σ ∈ R 4×4 is diagonal.Let Q 1 ∈ R 4×k be the matrix whose columns are the columns of U corresponding to σ 1 , and let Q 2 ∈ R 4×k be the matrix whose columns are the columns of V corresponding to σ 1 .It is not difficult to see that In fact, for any unit vectors y 1 and y 2 , the value of objective function of ( 24) at the point ( − → according to the Cauchy-Schwarz inequality.Without loss of generality, we assume σ 1 > 0. Then the equality holds if and only if y 1 = y 2 .As a result, problem (25) can be rewritten as an QCQP: In particular, when k = 1, problem (26) could be solved efficiently by representing − → x I and − → z 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 min Sc (g According to Corollary 2.3, we have Sc a since x, z, a (i) and b (i) are unit quaternions for i = 1, 2, . . ., n.It follows that Sc a and By simple computation, one can obtain that where K 11 , K 12 and K 21 are given by ( 23), ( 28) and ( 29), respectively.Under the constraints of problem (27), − → x st and − → z st are left-singular and right-singular vectors corresponding to the maximal singular value σ 1 for K 11 , which means Then we have − → x st = 0 under the constraints of problem (27).As a result, problem ( 27) is equivalent to the optimization Similarly, given the singular value decomposition K 11 = U ΣV , let Q 1 be the matrix whose columns are the columns of U corresponding to σ 1 , and let Q 2 be the matrix whose columns are the columns of V corresponding to σ 1 .The optimal − → x st and − → z st can be derived by computing the unit eigenvectors corresponding to the maximal eigenvalue of Sym Q 1 (K 12 + K 21 )Q 2 .Since the objective funcion in (30) does not contain − → x I and − → z I , the optimal − → x I can be any vector which is orthogonal to the optimal − → x st , and the optimal − → z I can be any vector which is orthogonal to the optimal − → z st .Considering the continuity of the norm, once the optimal − → x st and − → z st are determined, we try to find the best one in the optimal set of − → x I and − → z I such that the patching 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 Require: Measurements A (i) , B (i) n i=1 , regularization parameter γ.Ensure: The hand-eye transformation matrix X and robort-word transformation matrix Z.
Compute x st , z st , x I and z I by solving QCQP (26).

6:
Compute x st and z st by finding the unit eigenvector corresponding to the maximal eigenvalue of Sym Q 1 (K 12 + K 21 )Q 2 .

7:
Compute x I and z I by solving (31) with the optimal x st and z st .
8: end if 9: Compute X and Z from the dual quaternions x = x st + x I and z = z st + z I respectively.

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-2430M 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 use MOSEK [21] as SDP solver.Further, GloptiPoly can also recover the solution to the original problem and certify its optimality.We set the regularization parameter γ = 2 × 10 −6 .For hand-eye calibration model AX = XB, we compare our method with the direct estimation proposed by Tsai et al. [30] (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. [20] (denoted by "Malti10"), and the dual quaternion method using polynomial optimization proposed by Heller et al. [12] (denoted by "Heller14").
Numerical experiments are carried out as follows.First, the original homogeneous transformation matrices X and Ẑ in (2) Second, we generate n transformation matrices A (i) , i = 1, 2, . . ., n.Then the transformation matrix B (i) is computed by

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

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 han-eye calibration of SCARA robots [31].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,35].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 X and Ẑ 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 n = (0, 0, 1) .For A (1) , A (2) , A (3) , A (4) , the rotation angles are θ 1 = π 6 , θ 2 = π 3 , θ 3 = − π 6 and θ 4 − π 3 , while their translation vectors are randomly generated given by  3 and 4, respectively.

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 A (i) , B (i) 4 i=1 are the same with that in Subsection 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 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 Figures 2 and 3, respectively.

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, ∞-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 γ to improve the efficiency of the method?
4. Can we extend this method to the simultaneous localization and mapping problem?

Figure 1 :
Figure 1: The flow chart of proposed method.

Table 1 :
As described above, we have four measurements A (i) , B(i)for equation AX = ZB, and six motions Ã(s) , B(s) for equation AX = XB.The numerical results for AX = XB and AX = ZB with non-parallel rotation axis are reported in Tables1 and 2, respectively.The proposed Algorithms 1 and 2 show the best behavior in terms of estimation error.Note that the first three methods in Table1and the first three methods in Table2get 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 Table1, and Algorithm 2 may need more computation time to get the solution when compared with the first three methods in Table2.Numerical results for AX = XB with non-parallel rotation axis

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

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

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