Hand-Eye Calibration

Hand-eye calibration estimates the fixed rigid transform between a robot’s hand (its gripper, flange or end-effector) and a sensor rigidly mounted on it, usually a camera. A robot that grasps what its camera sees needs this transform: an object located in the camera frame must be carried into the hand frame, and from there into the robot’s base frame, before the arm can move to it. Horaud and Dornaika (1995, §1) name two uses: mapping sensor-centred measurements into the robot workspace, and letting the robot move the sensor precisely, as in inspection, 3D reconstruction with a moving camera and visual servoing.

The transform cannot be measured with a ruler, because the camera’s optical centre and axes are inside the camera. It is estimated instead from motion: the robot moves the camera between several stations, the robot’s own kinematics report how the hand moved, images of a fixed calibration target report how the camera moved, and the two motions are linked through the unknown transform. This gives the homogeneous matrix equation $AX = XB$, which nearly every method solves (Horaud and Dornaika, 1995, §1). The camera’s intrinsic parameters come first, from target-based camera calibration; hand-eye calibration estimates the extrinsic link between two bodies. The same equation, with other sensors in place of the camera, appears whenever two rigidly connected sensors each measure their own motion.

Problem Definition

Four frames are involved: the robot base $b$, the gripper (hand) $g$, the camera $c$ and a calibration target $t$ fixed in the workspace. As in the camera pose article, $T_{ab}$ is the $4 \times 4$ rigid transform that maps coordinates in frame $b$ to frame $a$, so that $T_{ac} = T_{ab}T_{bc}$ and $T_{ba} = T_{ab}^{-1}$. The group of these transforms, $SE(3)$, is described in geometric transformations.

At each station $i$ the robot’s forward kinematics give the gripper-to-base transform $T_{bg}^{(i)}$, and the image of the target, through perspective-n-point pose estimation, gives the target-to-camera transform $T_{ct}^{(i)}$. The unknown is the camera-to-gripper transform $X = T_{gc}$. Because the target does not move in the base frame,

$$T_{bt} = T_{bg}^{(i)}\,X\,T_{ct}^{(i)} \quad \text{for every station } i .$$

Equating stations $i$ and $j$ and rearranging gives

$$\underbrace{\big(T_{bg}^{(j)}\big)^{-1} T_{bg}^{(i)}}_{A_{ij}}\; X \;=\; X\; \underbrace{T_{ct}^{(j)} \big(T_{ct}^{(i)}\big)^{-1}}_{B_{ij}},$$

where $A_{ij}$ is the motion of the gripper between the two stations, expressed as the map from gripper coordinates at station $i$ to gripper coordinates at station $j$, and $B_{ij}$ is the corresponding motion of the camera. This is the form documented for OpenCV’s calibrateHandEye in the eye-in-hand case (OpenCV 4.10.0, calib3d.hpp), and the form the camera calibration article states, with $A$ the robot motion and $B$ the camera motion. Horaud and Dornaika (1995, §2.1) name the matrices the other way round: their $A$ is the camera motion, obtained from two calibrations of the camera against the target as $A = A_2 A_1^{-1}$; their $B$ is the hand motion, $B = B_2^{-1}B_1$ from the hand-to-base transforms; and their $X$ maps hand coordinates to camera coordinates. Inverting both sides of one form gives the other, so the two describe the same problem with $X$ inverted; Andreff, Horaud and Espiau (2001, §2.2) and the industrial machine vision article use Horaud and Dornaika’s labels. Code that mixes the conventions returns the inverse of the transform it was asked for, or nonsense.

Writing each matrix as a rotation $R$ and a translation $\mathbf{t}$, $AX = XB$ splits into a rotation equation and a translation equation (Horaud and Dornaika, 1995, eqs. 15 and 16):

$$R_A R_X = R_X R_B, \qquad (R_A - I)\,\mathbf{t}_X = R_X \mathbf{t}_B - \mathbf{t}_A .$$

The rotation equation says $R_A = R_X R_B R_X^\top$: the hand and camera rotations are similar matrices, so they have the same rotation angle, and their axes are related by $\mathbf{n}_A = R_X \mathbf{n}_B$ (eqs. 17 and 18). One motion therefore tells us where one direction of the camera frame points in the hand frame. The translation equation is linear in $\mathbf{t}_X$ once $R_X$ is known, but $R_A - I$ is singular, since it annihilates the rotation axis $\mathbf{n}_A$: one motion leaves the component of $\mathbf{t}_X$ along its own axis undetermined.

Eye-in-hand and eye-to-hand. In the configuration above, the camera rides on the hand and looks at a fixed target. In the eye-to-hand configuration the camera is fixed in the workspace and the target is carried by the gripper. OpenCV documents the same equation for it, with the base-to-gripper transforms in place of the gripper-to-base transforms and $X$ the camera-to-base transform (OpenCV 4.10.0, calib3d.hpp). The visual servoing article uses the same two names for camera placement; there the constant camera-to-end-effector transform enters the controller’s Jacobian, and its control loop tolerates a coarse estimate.

Robot-world and hand-eye calibration. The station-by-station equation $T_{bt} = T_{bg}^{(i)} X T_{ct}^{(i)}$ contains a second unknown, the target’s pose in the base frame. Treating both as unknowns, without forming relative motions, gives an equation of the form $AX = YB$, as in the title of Zhuang, Roth and Sudhakar (1994, not read). Dornaika and Horaud (1998, §I) write it $AX = ZB$, with $X$ the gripper-to-camera transform and $Z$ the robot-to-world transform, the world frame being the frame of the calibration pattern; they also describe a second set-up with a fixed camera observing a tool on the gripper. OpenCV’s calibrateRobotWorldHandEye solves the same equation with its own labels: $A$ the world-to-camera transform, $X$ the base-to-world transform, $Z$ the gripper-to-camera transform and $B$ the base-to-gripper transform, from at least three measurements (OpenCV 4.10.0, calib3d.hpp).

Input and Output

• Input: for each of $n \ge 3$ stations, the robot’s reported gripper pose $T_{bg}^{(i)}$ and an image of a calibration target with known geometry; the camera’s intrinsic parameters and lens distortion. Methods based on relative motions use the $n-1$ consecutive motions, as Horaud and Dornaika (1995, eq. 13) write them, or all $n(n-1)/2$ pairs, as OpenCV 4.10.0’s implementations do.
• Output: the rotation and translation of $X$, for example as OpenCV’s R_cam2gripper and t_cam2gripper; in the robot-world variant also the target’s (or world’s) pose in the base frame. Good practice reports residuals as well, discussed under evaluation below.

Variants

• Calibration without a target. Andreff et al. (2001) replace pose estimation against a calibration jig with structure from motion, so the camera motions are known only up to one global scale factor, which they add as an unknown to the linear formulation (§2.3, eq. 5; §3.1). They aim at calibration in the working site, for example for mobile and space robots that cannot carry a jig, with small motions; their conclusion reports a small loss of accuracy against classical calibration, good enough for visual servoing.
• Camera–LiDAR extrinsics. The rigid transform between a camera and a LiDAR is the same kind of unknown. The LiDAR article describes target-based methods, which make the LiDAR points on a checkerboard lie on the plane the camera sees, and targetless ones, such as maximizing mutual information between LiDAR reflectivity and image intensity. When each sensor can estimate its own motion, for example by visual odometry and by scan registration such as point cloud registration, the motion pairs satisfy $AX = XB$ as for a hand and a camera; Horaud and Dornaika (1995, §1) already count range finders among the sensors their formulations cover. For sensors on a ground vehicle, Andreff et al.’s planar-motion result applies in our reading: motions that all rotate about the vertical leave the offset along the vertical undetermined.
• Camera–IMU extrinsics. An inertial measurement unit measures angular velocity and specific force rather than poses, so its motions are not directly available as rigid transforms. Kalibr’s camera–IMU tool estimates the spatial and temporal parameters of a camera system with respect to an intrinsically calibrated IMU in a full batch optimization that models the system’s pose with splines; its documentation points to Furgale, Rehder and Siegwart (2013, not read) and to an earlier paper on continuous-time batch estimation for the approach (Kalibr wiki, accessed November 2026). It asks for a fixed target, motions that excite all IMU axes in rotation and translation, low motion blur, and timestamps from one clock with low jitter, and reports good results with a 20 Hz camera and a 200 Hz IMU; it warns that symmetric targets such as checkerboards can produce flips in the target pose estimates and recommends an AprilGrid. Temporal calibration is on by default. The camera calibration article summarizes the continuous-time approach, and visual-inertial odometry covers systems that estimate the camera–IMU extrinsics and time offset online.

Motion Conditions for a Unique Solution

One motion is never enough. Horaud and Dornaika (1995, §3) report that Shiu and Ahmad, Tsai and Lenz, and the other authors they summarize found at least three positions, that is two motions, necessary to determine $X$ uniquely. Andreff, Horaud and Espiau (2001, §1 and §3.3) state the condition more precisely: two calibration motions with non-parallel rotation axes are sufficient, which was shown algebraically by Tsai and Lenz and geometrically by Chen in a screw-motion analysis. Dornaika and Horaud (1998, §I) state the same condition for the robot-world problem: their methods, like the linear method of Zhuang et al., give a unique solution when the robot performs two motions with distinct rotation axes. OpenCV’s documentation states that a minimum of two motions with non-parallel rotation axes is necessary, so at least three poses, and strongly recommends many more.

The reason follows from the two equations. If every motion rotates about the same axis $\mathbf{n}_B$ in the camera frame, then any $X' = X\,T_k$, where $T_k$ is a rotation about $\mathbf{n}_B$ combined with a translation along it, satisfies every equation, because such a $T_k$ commutes with each camera motion. The rotation of $X$ about the common axis and its translation along that axis are left free. This is our derivation from the equations; the example below checks it numerically.

Andreff et al. (2001, §3.3, Table 1) work out what partial motions do determine, for camera motions that may be known only up to scale:

• Pure translations of the hand ($R_A = I$) produce pure translations of the camera, with $\mathbf{t}_A = R_X \mathbf{t}_B$ in our labels. Three independent pure translations give the rotation of $X$ (and the unknown scale) linearly, but the translation of $X$ cannot be observed (Proposition 1).
• Pure rotations about the wrist centre, at least two with non-parallel axes, give the rotation and the translation of $X$, the translation up to the unknown scale; when camera motions come from pose estimation against a known target, the scale is known and the translation is fully recovered (Proposition 3).
• Planar motions, as for a car-like robot whose rotations all share one axis, combined with a pure translation not parallel to that axis, give the rotation and the scale, but the translation only up to an unknown height along the normal to the plane of motion (Lemma 2).
• Two independent general motions with non-parallel axes determine the whole transform (Proposition 5).

Beyond uniqueness, the motions control accuracy. Andreff et al. (2001, §1) note that earlier methods used minimal rotation parameterizations, which are ill-defined for small rotation angles, so their calibration motions had to be as large as possible, and that Tsai and Lenz even gave a rule for choosing large motions. OpenCV 4.10.0’s implementation of the Tsai–Lenz method discards motions whose rotation vector, in the form $2\sin(\theta/2)\,\mathbf{n}$, has a norm below 0.3 or above 1.7. A source comment says these limits correspond to angles under 17° or over 120°, chosen by checking the method’s equation on simulated data; by our arithmetic a norm of 1.7 is about 116°.

Challenges

• Two sources of error. The camera motions carry the errors of camera calibration against the target, and the hand motions carry the errors of the robot’s kinematic model, whose parameters are never perfect (Horaud and Dornaika, 1995, §6). In their real data, Horaud and Dornaika judged camera calibration errors negligible next to robot calibration errors (§7).
• Error propagation through decoupling. Solving for rotation first and translation second passes the rotation error into the translation; Andreff et al. (2001, §2.2) give this as the motivation for simultaneous methods. Horaud and Dornaika (1995, §3) add that the geometric properties that allow the rotation equation to be linearized do not hold exactly in the presence of noise, so the linear problem can become ill-conditioned.
• Weak motions. Near-parallel axes and small rotations make the problem nearly degenerate, as the previous section describes; a fit with small residuals can still be far from the true transform.
• Units and weighting. Simultaneous methods mix rotation residuals, which have no unit, with translation residuals in metres or millimetres. Andreff et al. (2001, §3.3.4) cite Wei et al.’s observation that the accuracy of simultaneous estimation depends on the physical unit used for translation.
• Conventions. The direction of each transform, the labelling of $A$ and $B$, and quaternion order differ between papers and libraries, as the Problem Definition shows; rotation representations lists the library conventions for quaternions.
• Fixed intrinsics. Methods assume the camera’s intrinsic parameters do not change during calibration; Horaud and Dornaika (1995, §8) call this assumption not very realistic.

Approaches

Separable closed-form methods

The early solutions decompose $AX = XB$ into the rotation equation and the translation equation and solve each by linear algebra (Horaud and Dornaika, 1995, §3). Shiu and Ahmad (1989) cast the rotation problem as eight linear equations in four unknowns, as Horaud and Dornaika summarize it. Tsai and Lenz (1989) represent $R_X$ by its axis $\mathbf{n}_X$ and angle $\theta_X$ and use the fact that, for a rotation taking $\mathbf{n}_B$ to $\mathbf{n}_A$, the vector $\mathbf{n} = \tan(\theta_X/2)\,\mathbf{n}_X$ satisfies a linear equation in which each motion contributes a rank-deficient $3 \times 3$ block, so that at least two motions are needed (Horaud and Dornaika, 1995, eq. 19). With rotation matrices acting on column vectors, the identity is

$$\mathbf{n}_A - \mathbf{n}_B = \mathbf{n} \times (\mathbf{n}_A + \mathbf{n}_B),$$

which follows from writing $R_X$ in the Cayley form $R_X = (I + [\mathbf{n}]_\times)(I - [\mathbf{n}]_\times)^{-1}$ (our derivation). Horaud and Dornaika’s eq. 19, in the HAL copy, prints the cross product in the opposite order, which corresponds to the transposed rotation. OpenCV 4.10.0, which follows Tsai’s paper through a port of Lazarevic’s MATLAB code according to its source comments, uses the motions’ vectors $\mathbf{P} = 2\sin(\theta/2)\,\mathbf{n}$ in place of unit axes, which weights each motion by its angle. It solves $[\mathbf{P}_A + \mathbf{P}_B]_\times\,\mathbf{P}' = \mathbf{P}_B - \mathbf{P}_A$ by least squares over all motions, recovers $\mathbf{P}_X = 2\mathbf{P}'/\sqrt{1 + \|\mathbf{P}'\|^2}$, and then solves the stacked translation equations, also by least squares; its comments label these equations 12, 14 and 15 of the paper.

Park and Martin (1994, not read) solve the problem on the Euclidean group, as their title says. In OpenCV’s implementation, each motion’s rotation vectors $\boldsymbol{\alpha}_i = \log R_{A_i}$ and $\boldsymbol{\beta}_i = \log R_{B_i}$ (axis times angle) are accumulated into $M = \sum_i \boldsymbol{\beta}_i \boldsymbol{\alpha}_i^\top$, and the rotation is

$$R_X = (M^\top M)^{-1/2} M^\top ,$$

followed by the same linear least squares for the translation. Because $\boldsymbol{\alpha}_i = R_X \boldsymbol{\beta}_i$ for exact data, this is the orthogonal Procrustes problem of aligning two sets of vectors, solved by a polar decomposition, as in the Kabsch algorithm for point sets (our reading of the formula). It also shows the degeneracy directly: if all axes are parallel, $M$ has rank one and $M^\top M$ cannot be inverted.

Quaternion methods

Horaud and Dornaika (1995) reduce both their formulations to the pair of equations $\mathbf{v}'_i = R\,\mathbf{v}_i$ and $(K_i - I)\,\mathbf{t} = R\,\mathbf{p}_i - \mathbf{p}'_i$ (eqs. 25 and 26), and minimize

$$f_1(R) = \sum_i \|\mathbf{v}'_i - R\,\mathbf{v}_i\|^2, \qquad f_2(R, \mathbf{t}) = \sum_i \|R\,\mathbf{p}_i - (K_i - I)\,\mathbf{t} - \mathbf{p}'_i\|^2$$

(eqs. 27 and 28), where for $AX = XB$ the vectors $\mathbf{v}'_i$ and $\mathbf{v}_i$ are the unit rotation axes of the two motions. With $R$ written as a unit quaternion $q$, $f_1 = q^\top \mathcal{A}\, q$ for a $4 \times 4$ symmetric matrix $\mathcal{A}$ built from the axes, so the optimal rotation is the eigenvector of $\mathcal{A}$ with the smallest eigenvalue, a closed form they attribute to Faugeras and Hébert; the translation then follows by linear least squares (§5.1). This is the eigenvector construction of Horn’s quaternion method for point sets in the Kabsch algorithm article, with the matrix set up so that the smallest eigenvalue, rather than Horn’s largest, gives the answer. Their second method minimizes $\lambda_1 f_1 + \lambda_2 f_2$ over rotation and translation together, with a penalty $\lambda(1 - q^\top q)^2$ (their $\lambda$, a third weight) that keeps $q$ near unit length, by Levenberg–Marquardt; all their experiments used $\lambda_1 = \lambda_2 = 1$ and $\lambda = 2 \times 10^6$ (§5.2, eq. 30).

Horaud and Dornaika also give a second formulation that avoids decomposing the camera’s $3 \times 4$ projection matrix into intrinsic and extrinsic parameters. With $Y$ the transform from the hand frame to the calibration frame at the first station, $X = A_1 Y$ in their notation, and substituting into $AX = XB$ gives

$$M_2\, Y = M_1\, Y B$$

for the projection matrices $M_1$ and $M_2$ of the two stations (§2.2, eq. 7). They show that this equation decomposes into the same two equations as the classical one (§4), so both solvers apply. Its advantage is that a perfect pinhole decomposition is no longer assumed, and they note that the formulation together with the classical one covers other camera-based sensors, including single scan-line cameras, stereo heads and range finders (§1).

Daniilidis (1999, not read) uses dual quaternions, as the title says, and OpenCV’s documentation lists the method among those that estimate rotation and translation simultaneously. In OpenCV 4.10.0’s implementation, each motion contributes six linear equations in the eight components of the dual quaternion of $X$, built from the vector parts of the two motions’ rotation quaternions $\mathbf{a}, \mathbf{b}$ and of their dual parts $\mathbf{a}', \mathbf{b}'$:

$$\begin{bmatrix} \mathbf{a} - \mathbf{b} & [\mathbf{a} + \mathbf{b}]_\times & \mathbf{0} & 0_{3 \times 3} \\ \mathbf{a}' - \mathbf{b}' & [\mathbf{a}' + \mathbf{b}']_\times & \mathbf{a} - \mathbf{b} & [\mathbf{a} + \mathbf{b}]_\times \end{bmatrix} \begin{bmatrix} q \\ q' \end{bmatrix} = \mathbf{0}.$$

The stacked matrix is decomposed by the SVD; the solution is a combination $\lambda_1 \mathbf{v}_7 + \lambda_2 \mathbf{v}_8$ of the right singular vectors of the two smallest singular values, and the two coefficients follow from the constraints that $q$ has unit norm and is orthogonal to $q'$, a quadratic equation. The code’s comments label these steps equations 31, 34 and 35. Andreff et al. (2001, §2.2) describe this line of work, citing a 1996 conference paper by Daniilidis and Bayro-Corrochano, as a dual-quaternion method that solves simultaneously for hand-eye rotation and translation.

Linear formulation with Kronecker products

Andreff et al. (2001, §3) keep the rotation as a full orthogonal matrix and write the rotation equation linearly in its nine entries with the Kronecker product: $(I_9 - R_{A} \otimes R_{B})\,\mathrm{vec}(R_X) = \mathbf{0}$, with $\mathrm{vec}$ stacking a matrix row by row (their eq. 26, which has this form in either labelling of $A$ and $B$). Stacked with the translation equations, this gives one linear system whose null space contains the solution; they argue that it handles small rotations, for which minimal parameterizations are ill-conditioned, and that it allows the algebraic analysis of partial calibration summarized above. A direct solution does not guarantee that the estimated $R_X$ is orthogonal, and correcting it afterwards cannot be carried consistently into the translation, so they prefer a two-step solution in the general case (§3.3.4). OpenCV’s CALIB_HAND_EYE_ANDREFF solves the stacked system by SVD and then normalizes the rotation, citing their 1999 conference paper.

Simultaneous and nonlinear refinement

Horaud and Dornaika’s stability analysis added uniform or Gaussian noise of 1% to 6% to the hand and camera motions and compared the Tsai–Lenz method, their closed form and their nonlinear method over 1,000 trials for each noise level (§6). The Tsai–Lenz and closed-form methods behaved almost the same, while the nonlinear method gave more accurate results in all the situations tested; they concluded that decoupling rotation and translation introduces a bias. With four motions, 6% rotational noise and 2% translational noise, the translation error was 4% for the nonlinear method and 6.5% for the other two. On real data, from 17, 7 and 6 robot positions, the results seemed to them to confirm that the nonlinear method estimates translation better, sometimes at the cost of a slightly less accurate rotation, and that the projection-matrix formulation is more accurate than the classical one (§7). They explain the rotation result by the robot’s imperfect kinematic calibration: the linear and closed-form methods estimate rotation independently of the robot’s translations, which carry those errors (§7).

For the robot-world problem, Dornaika and Horaud (1998, §II–III) describe Zhuang et al.’s linear quaternion solution, which first solves linearly for the components of two quaternions and then normalizes them, and which needs the scalar parts of two of the quaternions to be non-zero and is not feasible in some special configurations; they propose a closed form that minimizes a positive semi-definite quadratic form in the two unit quaternions jointly, with two Lagrange multipliers, so that one quaternion is an eigenvector of a $4 \times 4$ matrix built from the motions and the other follows from it (eqs. 9–15), and a nonlinear method that estimates both rotations and translations together. In their simulated and real experiments the nonlinear method performed better than the closed form, which performed slightly better than the linear method (§I).

Any of these can be followed by a refinement that minimizes a cost closer to the measurements, such as the reprojection error of the target corners in every image, with $X$ and the target’s pose in the base frame as unknowns. This is a nonlinear least-squares problem of the same kind as bundle adjustment, and the closed-form solvers provide its starting point. Which cost to refine, and how to weight rotation against translation, is a choice the sources read here leave open; the example prints two.

Implementations

OpenCV’s calibrateHandEye implements the Tsai–Lenz, Park–Martin, Horaud–Dornaika, Andreff–Horaud–Espiau and Daniilidis methods, with Tsai–Lenz as the default, takes the gripper-to-base and target-to-camera transforms of at least three stations, and returns the camera-to-gripper transform (OpenCV 4.10.0). Every method forms all pairs of stations. Its calibrateRobotWorldHandEye implements two methods for $AX = ZB$, which its documentation attributes to Shah (2013), using the Kronecker product, and to Li et al. (2010), using dual quaternions and the Kronecker product; neither paper was read for this article. Both functions are declared with CV_EXPORTS_W in 4.10.0, the macro that exports a function to the bindings; the camera calibration article notes that OpenCV 5.0.0 declares calibrateHandEye for C++ only. The OpenCV source also gives the Daniilidis paper’s year as 1998, while the journal issue (volume 18, number 3) is from 1999.

Example

The script below calibrates an eye-in-hand camera from rendered images of a checkerboard and asks how the answer degrades as the rotation axes of the robot’s motions are squeezed toward one direction. The camera has a focal length of 800 pixels and a $640 \times 480$ image, with known intrinsics and no distortion; the target has $9 \times 7$ inner corners on 25 mm squares. The true camera-to-gripper transform $X = T_{gc}$ and the target’s pose in the base frame are fixed by hand. Each trial uses ten stations: at station $i$ the camera is rotated from a fronto-parallel reference view by an angle of 15° to 45° (random sign) about an axis $\mathbf{u}_i$ drawn uniformly inside a cone of half-angle $s$ around the optical axis, and placed 0.5 m from the board centre, looking at it. The gripper pose that puts the camera there, $T_{bg}^{(i)} = T_{bt}\,(T_{ct}^{(i)})^{-1} X^{-1}$, plays the robot’s report and is exact. The conventions are those of the Problem Definition: $T_{ab}$ maps frame $b$ to frame $a$, $T_{ct}$ is target-to-camera (world-to-camera with the target as world), $X$ is camera-to-gripper, $A$ is the gripper motion and $B$ the camera motion over all 45 pairs of stations; SciPy’s quaternions are scalar-last and the dual-quaternion code uses scalar-first quaternions.

Each image is rendered by inverse homography with $3 \times 3$ supersampling and Gaussian pixel noise of standard deviation 0.01 on a 0–1 scale. The corners are located with scikit-image’s corner_subpix, started from the rounded true projections, which stands in for a corner detector. Each target pose comes from a planar perspective-n-point step: a homography by the normalized DLT, its decomposition, and Levenberg–Marquardt on the pixel residuals. Five estimators are compared: Tsai–Lenz and Park–Martin as OpenCV 4.10.0 implements them (without OpenCV’s discarding of small and large motions), the Daniilidis dual-quaternion solver following OpenCV’s construction, Horaud and Dornaika’s joint cost $f_1 + f_2$ with unit weights, and a reprojection refinement over $X$ and the target pose with the robot poses held exact. Both refinements start from the Park–Martin estimate. Our choices in the motion-space refinement are a rotation-vector parameterization instead of Horaud and Dornaika’s quaternion with a penalty term, and translations in metres. The spread of the motion axes is measured by the ratio of the second to the first singular value of the stacked rotation vectors of the robot motions, 0 when all axes are parallel; small motions, whose axes are poorly defined, get little weight in it.

Part 1 checks the solvers with exact corners: with well-spread axes all five must recover $X$, and with parallel axes ($s = 0$) the script does not solve, because the sources state that the solution is then not unique. It prints the rank of the Tsai–Lenz rotation system, the eigenvalue ratio of Park–Martin’s $M^\top M$ and the null-space dimension of the Daniilidis system, and checks that a second transform $X' = X\,T_k$, rotated by 30° about the camera’s optical axis and shifted 50 mm along it, satisfies every equation as well. Part 2 repeats the calibration from rendered images for ten trials at each spread and prints median errors. The script uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2; it renders and processes 500 images, takes about three minutes on the machine used for drafting, and printed the same output on two runs.

import numpy as np
from scipy.optimize import least_squares
from scipy.spatial.transform import Rotation as Rot
from skimage.feature import corner_subpix

# Conventions: T_ab maps coordinates in frame b to frame a. Frames: base b, gripper g, camera c, target t.
# Robot gives T_bg (gripper-to-base); PnP gives T_ct (target-to-camera). Unknown X = T_gc (camera-to-gripper).
K = np.array([[800.0, 0, 320], [0, 800, 240], [0, 0, 1]])
W, H, NX, NY, SQ = 640, 480, 9, 7, 0.025                    # 9 x 7 inner corners, 25 mm squares
gx, gy = np.meshgrid(np.arange(1, NX + 1) * SQ, np.arange(1, NY + 1) * SQ)
P_T = np.c_[gx.ravel(), gy.ravel(), np.zeros(gx.size)]      # corners in the target frame (z_t = 0)
CENTRE = np.array([(NX + 1) * SQ / 2, (NY + 1) * SQ / 2, 0.0])

def se3(R, t):
    T = np.eye(4); T[:3, :3] = R; T[:3, 3] = t
    return T

def inv(T):
    return se3(T[:3, :3].T, -T[:3, :3].T @ T[:3, 3])

def rotvec(R):
    return Rot.from_matrix(R).as_rotvec()

def expm(v):
    return Rot.from_rotvec(v).as_matrix()

def project(T_ct, P):
    x = (K @ (T_ct[:3, :3] @ P.T + T_ct[:3, 3:])).T
    return x[:, :2] / x[:, 2:]

def render(T_ct, rng, ss=3, noise=0.01):
    """Checkerboard image by inverse homography, ss x ss supersampling, Gaussian pixel noise."""
    Hi = np.linalg.inv(K @ np.c_[T_ct[:3, :2], T_ct[:3, 3]])      # pixel -> target plane
    img = np.full((H, W), 0.5)
    box = project(T_ct, np.array([[-SQ, -SQ, 0], [(NX + 2) * SQ, -SQ, 0],
                                  [-SQ, (NY + 2) * SQ, 0], [(NX + 2) * SQ, (NY + 2) * SQ, 0]]))
    u0, v0 = np.clip(np.floor(box.min(0)) - 2, 0, [W, H]).astype(int)
    u1, v1 = np.clip(np.ceil(box.max(0)) + 2, 0, [W, H]).astype(int)
    o = (np.arange(ss) + 0.5) / ss - 0.5
    u = np.arange(u0, u1)[None, :, None, None] + o[None, None, None, :]
    v = np.arange(v0, v1)[:, None, None, None] + o[None, None, :, None]
    u, v = np.broadcast_arrays(u, v)
    p = Hi @ np.stack([u.ravel(), v.ravel(), np.ones(u.size)])
    X, Y = p[0] / p[2], p[1] / p[2]
    i, j = np.floor(X / SQ), np.floor(Y / SQ)
    inside = (i >= 0) & (i <= NX) & (j >= 0) & (j <= NY)
    margin = (X > -SQ) & (X < (NX + 2) * SQ) & (Y > -SQ) & (Y < (NY + 2) * SQ)
    val = np.where(inside, np.where((i + j) % 2 == 0, 0.2, 0.8), np.where(margin, 0.8, 0.5))
    img[v0:v1, u0:u1] = val.reshape(u.shape).mean(axis=(2, 3))
    return img + noise * rng.standard_normal(img.shape)

def locate(img, uv_guess):
    """Sub-pixel corners, searched from the rounded true projections (a stand-in for a detector)."""
    rc = corner_subpix(img, np.round(uv_guess[:, ::-1]).astype(int), window_size=11)
    return rc[:, ::-1]

def pnp(uv):
    """Planar PnP: homography by normalized DLT, decomposition, then Levenberg-Marquardt on pixels."""
    m = (np.linalg.inv(K) @ np.c_[uv, np.ones(len(uv))].T).T[:, :2]
    rows = []
    for (x, y, _), (u, v) in zip(P_T, m):
        rows += [[x, y, 1, 0, 0, 0, -u * x, -u * y, -u], [0, 0, 0, x, y, 1, -v * x, -v * y, -v]]
    h = np.linalg.svd(np.array(rows))[2][-1].reshape(3, 3)
    h /= np.linalg.norm(h[:, 0]) * np.sign(h[2, 2])
    U, _, Vt = np.linalg.svd(np.c_[h[:, 0], h[:, 1], np.cross(h[:, 0], h[:, 1])])
    R0 = U @ np.diag([1, 1, np.linalg.det(U @ Vt)]) @ Vt
    f = lambda p: (project(se3(expm(p[:3]), p[3:]), P_T) - uv).ravel()
    p = least_squares(f, np.r_[rotvec(R0), h[:, 2]], method="lm", xtol=1e-15, ftol=1e-15).x
    return se3(expm(p[:3]), p[3:])

def motions(T_bg, T_ct):
    """All pairs i < j: A = gripper motion, B = camera motion, with A X = X B."""
    A, B = [], []
    for i in range(len(T_bg)):
        for j in range(i + 1, len(T_bg)):
            A.append(inv(T_bg[j]) @ T_bg[i])
            B.append(T_ct[j] @ inv(T_ct[i]))
    return A, B

def skew(v):
    return np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])

def translation(A, B, R):
    """(R_A - I) t_X = R_X t_B - t_A, stacked over all motions."""
    M = np.vstack([a[:3, :3] - np.eye(3) for a in A])
    d = np.concatenate([R @ b[:3, 3] - a[:3, 3] for a, b in zip(A, B)])
    return np.linalg.lstsq(M, d, rcond=None)[0]

def tsai_lenz(A, B, rank_only=False):
    def P(T):                                                  # 2 sin(theta/2) * axis
        q = Rot.from_matrix(T[:3, :3]).as_quat()               # SciPy order (x, y, z, w)
        return 2 * q[:3] * np.sign(q[3])
    Pa, Pb = [P(a) for a in A], [P(b) for b in B]
    M = np.vstack([skew(pa + pb) for pa, pb in zip(Pa, Pb)])
    if rank_only:
        s = np.linalg.svd(M, compute_uv=False)
        return int((s > 1e-9 * s[0]).sum())
    p_ = np.linalg.lstsq(M, np.concatenate([pb - pa for pa, pb in zip(Pa, Pb)]), rcond=None)[0]
    p = 2 * p_ / np.sqrt(1 + p_ @ p_)                          # 2 sin(theta/2) * axis of R_X
    R = Rot.from_quat(np.r_[p / 2, np.sqrt(1 - p @ p / 4)]).as_matrix()
    return se3(R, translation(A, B, R))

def park_martin(A, B, eig_only=False):
    M = sum(np.outer(rotvec(b[:3, :3]), rotvec(a[:3, :3])) for a, b in zip(A, B))
    w, V = np.linalg.eigh(M.T @ M)
    if eig_only:
        return w.min() / w.max()
    R = V @ np.diag(w ** -0.5) @ V.T @ M.T                     # (M^T M)^(-1/2) M^T
    return se3(R, translation(A, B, R))

def qmul(s, t):                                                # scalar-first quaternions
    return np.r_[s[0] * t[0] - s[1:] @ t[1:], s[0] * t[1:] + t[0] * s[1:] + np.cross(s[1:], t[1:])]

def dual_quat(T):
    q = Rot.from_matrix(T[:3, :3]).as_quat()[[3, 0, 1, 2]]
    q = q if q[0] >= 0 else -q
    return q, 0.5 * qmul(np.r_[0, T[:3, 3]], q)

def daniilidis(A, B, null_only=False):
    rows = []
    for a, b in zip(A, B):
        (qa, qa_), (qb, qb_) = dual_quat(a), dual_quat(b)
        u, u_, v, v_ = qa[1:], qa_[1:], qb[1:], qb_[1:]
        S = np.zeros((6, 8))
        S[:3, 0], S[:3, 1:4] = u - v, skew(u + v)
        S[3:, 0], S[3:, 1:4], S[3:, 4], S[3:, 5:] = u_ - v_, skew(u_ + v_), u - v, skew(u + v)
        rows.append(S)
    s, Vt = np.linalg.svd(np.vstack(rows))[1:]
    if null_only:
        return int((s < 1e-9 * s[0]).sum()) + 8 - len(s)
    v1, v2 = Vt[-2], Vt[-1]
    a_, b_, c_ = v1[:4] @ v1[4:], v1[:4] @ v2[4:] + v2[:4] @ v1[4:], v2[:4] @ v2[4:]
    best = None
    for r in np.roots([a_, b_, c_]).real:                      # unit-norm and orthogonality constraints
        val = r * r * v1[:4] @ v1[:4] + 2 * r * v1[:4] @ v2[:4] + v2[:4] @ v2[:4]
        if best is None or val > best[1]:
            best = (r, val)
    l2 = np.sqrt(1 / best[1]); l1 = best[0] * l2
    dq = l1 * v1 + l2 * v2
    q, q_ = dq[:4], dq[4:]
    t = 2 * qmul(q_, q * np.r_[1, -1, -1, -1])[1:]
    return se3(Rot.from_quat(q[[1, 2, 3, 0]]).as_matrix(), t)

def refine_motion(A, B, X0):
    """Horaud-Dornaika joint cost: sum ||n_A - R n_B||^2 + ||R t_B - (R_A - I) t - t_A||^2 (metres)."""
    nA = [rotvec(a[:3, :3]) / np.linalg.norm(rotvec(a[:3, :3])) for a in A]
    nB = [rotvec(b[:3, :3]) / np.linalg.norm(rotvec(b[:3, :3])) for b in B]
    def f(p):
        R, t = expm(p[:3]), p[3:]
        return np.concatenate([np.r_[na - R @ nb, R @ b[:3, 3] - (a[:3, :3] - np.eye(3)) @ t - a[:3, 3]]
                               for na, nb, a, b in zip(nA, nB, A, B)])
    p = least_squares(f, np.r_[rotvec(X0[:3, :3]), X0[:3, 3]], method="lm").x
    return se3(expm(p[:3]), p[3:])

def refine_reproj(T_bg, uvs, X0, T_ct0):
    """Reprojection error over X and the target pose T_bt, robot poses taken as exact."""
    Y0 = T_bg[0] @ X0 @ T_ct0
    def f(p):
        X, Y = se3(expm(p[:3]), p[3:6]), se3(expm(p[6:9]), p[9:])
        return np.concatenate([(project(inv(Tg @ X) @ Y, P_T) - uv).ravel() for Tg, uv in zip(T_bg, uvs)])
    p = least_squares(f, np.r_[rotvec(X0[:3, :3]), X0[:3, 3], rotvec(Y0[:3, :3]), Y0[:3, 3]], method="lm").x
    return se3(expm(p[:3]), p[3:6])

def errors(Xe, Xt):
    return np.degrees(np.linalg.norm(rotvec(Xe[:3, :3] @ Xt[:3, :3].T))), 1000 * np.linalg.norm(Xe[:3, 3] - Xt[:3, 3])

X_TRUE = se3(Rot.from_euler("xyz", [10, -20, 90], degrees=True).as_matrix(), [0.05, -0.03, 0.10])
T_BT = se3(Rot.from_euler("z", 30, degrees=True).as_matrix(), [0.6, 0.1, 0.0])
R_REF = Rot.from_rotvec([np.pi, 0, 0]).as_matrix()            # reference view: camera faces the board

def stations(spread_deg, seed, n=10, dist=0.5):
    """Camera rotations theta_i about axes u_i inside a cone of half-angle `spread` around the optical axis."""
    rng = np.random.default_rng(seed)
    U, phi, mag, sgn = rng.random(n), rng.uniform(0, 2 * np.pi, n), rng.uniform(15, 45, n), rng.choice([-1, 1], n)
    ca = 1 - U * (1 - np.cos(np.radians(spread_deg)))
    sa = np.sqrt(1 - ca ** 2)
    T_ct, T_bg = [], []
    for k in range(n):
        u = np.array([sa[k] * np.cos(phi[k]), sa[k] * np.sin(phi[k]), ca[k]])
        R_tc = R_REF @ expm(np.radians(sgn[k] * mag[k]) * u)   # camera-to-target rotation
        Tct = inv(se3(R_tc, CENTRE - dist * R_tc[:, 2]))       # camera looks at the board centre
        T_ct.append(Tct)
        T_bg.append(T_BT @ inv(Tct) @ inv(X_TRUE))             # forward kinematics that put it there
    return T_bg, T_ct

def axis_spread(A):
    """Second over first singular value of the stacked rotation vectors: 0 when all axes are parallel."""
    sv = np.linalg.svd(np.array([rotvec(a[:3, :3]) for a in A]), compute_uv=False)
    return sv[1] / sv[0]

NAMES = ["Tsai-Lenz", "Park-Martin", "Daniilidis", "refine motion", "refine reproj"]

def solve_all(T_bg, T_ct, uvs):
    A, B = motions(T_bg, T_ct)
    Xp = park_martin(A, B)
    return [tsai_lenz(A, B), Xp, daniilidis(A, B), refine_motion(A, B, Xp), refine_reproj(T_bg, uvs, Xp, T_ct[0])]

print("Part 1: exact corners")
for s in [90, 0]:
    T_bg, T_true = stations(s, seed=0)
    uvs = [project(T, P_T) for T in T_true]
    T_ct = [pnp(uv) for uv in uvs]
    A, B = motions(T_bg, T_ct)
    print(f"spread {s:2d} deg, axis spread ratio of the robot motions {axis_spread(A):.1e}")
    if s > 0:
        for name, Xe in zip(NAMES, solve_all(T_bg, T_ct, uvs)):
            print(f"  {name:14s} rotation error {errors(Xe, X_TRUE)[0]:.1e} deg, translation error {errors(Xe, X_TRUE)[1]:.1e} mm")
    else:
        print(f"  Tsai-Lenz rotation system rank {tsai_lenz(A, B, rank_only=True)} of 3;"
              f" Park-Martin min/max eigenvalue of M^T M {park_martin(A, B, eig_only=True):.1e};"
              f" Daniilidis null-space dimension {daniilidis(A, B, null_only=True)} (2 when unique)")
        X2 = X_TRUE @ se3(Rot.from_euler("z", 30, degrees=True).as_matrix(), [0, 0, 0.05])
        res = max(np.abs(a @ X2 - X2 @ b).max() for a, b in zip(A, B))
        print(f"  X' = X Rot_z(30 deg) Trans_z(50 mm) also satisfies every A X' = X' B: max residual {res:.1e}")

print("\nPart 2: rendered images, 10 stations, median over 10 trials")
rows, corner_err = {}, []
for s in [90, 30, 10, 3, 1, 0]:
    errs, spreads = [], []
    for trial in range(10):
        T_bg, T_true = stations(s, seed=trial)
        if axis_spread(motions(T_bg, T_bg)[0]) < 1e-9:        # commanded robot axes parallel
            break
        rng = np.random.default_rng(1000 + trial)
        uvs = []
        for T in T_true:
            uv = locate(render(T, rng), project(T, P_T))
            corner_err.append(np.linalg.norm(uv - project(T, P_T), axis=1))
            uvs.append(uv)
        T_ct = [pnp(uv) for uv in uvs]
        spreads.append(axis_spread(motions(T_bg, T_ct)[0]))
        errs.append([errors(Xe, X_TRUE) for Xe in solve_all(T_bg, T_ct, uvs)])
    rows[s] = (np.median(spreads), np.median(np.array(errs), axis=0)) if errs else None
print(f"corner location error: median {np.median(np.concatenate(corner_err)):.3f} px")
for k, (label, unit) in enumerate([("rotation error", "deg"), ("translation error", "mm")]):
    print(f"\n{label} ({unit})  ratio" + "".join(f"{n:>15s}" for n in NAMES))
    for s, r in rows.items():
        if r is None:
            print(f"spread {s:2d} deg        0   parallel axes: X not unique, not solved")
            continue
        print(f"spread {s:2d} deg  {r[0]:7.4f}" + "".join(f"{v:15.3f}" for v in r[1][:, k]))

Output:

Part 1: exact corners
spread 90 deg, axis spread ratio of the robot motions 7.7e-01
  Tsai-Lenz      rotation error 5.2e-14 deg, translation error 3.7e-13 mm
  Park-Martin    rotation error 1.1e-14 deg, translation error 2.5e-13 mm
  Daniilidis     rotation error 7.6e-14 deg, translation error 3.4e-13 mm
  refine motion  rotation error 4.9e-15 deg, translation error 2.4e-13 mm
  refine reproj  rotation error 3.2e-14 deg, translation error 1.4e-13 mm
spread  0 deg, axis spread ratio of the robot motions 5.9e-17
  Tsai-Lenz rotation system rank 2 of 3; Park-Martin min/max eigenvalue of M^T M -1.1e-16; Daniilidis null-space dimension 4 (2 when unique)
  X' = X Rot_z(30 deg) Trans_z(50 mm) also satisfies every A X' = X' B: max residual 1.5e-15

Part 2: rendered images, 10 stations, median over 10 trials
corner location error: median 0.038 px

rotation error (deg)  ratio      Tsai-Lenz    Park-Martin     Daniilidis  refine motion  refine reproj
spread 90 deg   0.7788          0.011          0.009          0.011          0.010          0.005
spread 30 deg   0.3084          0.023          0.022          0.023          0.030          0.009
spread 10 deg   0.1010          0.094          0.094          0.093          0.087          0.059
spread  3 deg   0.0302          0.290          0.280          0.279          0.777          0.130
spread  1 deg   0.0101          0.824          0.972          0.984          1.293          0.832
spread  0 deg        0   parallel axes: X not unique, not solved

translation error (mm)  ratio      Tsai-Lenz    Park-Martin     Daniilidis  refine motion  refine reproj
spread 90 deg   0.7788          0.125          0.117          0.145          0.116          0.032
spread 30 deg   0.3084          0.166          0.166          0.174          0.160          0.043
spread 10 deg   0.1010          0.965          0.965          1.083          0.963          0.077
spread  3 deg   0.0302          2.369          2.368          2.503          2.371          0.171
spread  1 deg   0.0101          8.453          8.453          9.217          8.422          0.188
spread  0 deg        0   parallel axes: X not unique, not solved

With exact corners and well-spread axes, every method returns $X$ to rounding error, so the implementations and the convention for $A$, $B$ and $X$ agree. With parallel axes the degeneracy appears in each formulation in the way the theory predicts: the Tsai–Lenz system loses one rank, Park–Martin’s $M^\top M$ becomes singular, the Daniilidis system has a four-dimensional null space instead of two, and a transform 30° and 50 mm away from the true one satisfies all 45 equations to rounding. No amount of image accuracy fixes this, since the commanded robot motions themselves carry no information about the free rotation and translation.

With rendered images, the corners are located to a few hundredths of a pixel, and with axes spread over a hemisphere every method finds the rotation within about 0.01° and the translation within about 0.15 mm. As the spread ratio falls from 0.1 to 0.01, the errors of the closed-form methods grow roughly in proportion to its inverse, to about 1° and 8–9 mm: the near-degenerate direction is estimated from the small differences between nearly parallel axes, and noise in those differences is amplified. The three closed-form methods behave almost alike throughout, as Horaud and Dornaika found for Tsai–Lenz and their closed form, and the dual-quaternion solver, which estimates rotation and translation together, is no better here than the separable ones.

The two refinements differ. Horaud and Dornaika’s joint cost brings no clear improvement over Park–Martin in this set-up: its translation errors are almost the same, and its rotation is worse at spread ratios of 0.31, 0.03 and 0.01; in our reading, the cost gives every motion’s unit axis the same weight, including the axes of small motions, which are dominated by noise, whereas Tsai–Lenz and Park–Martin weight each axis by its rotation angle. This differs from Horaud and Dornaika’s simulations, where noise was added to the motions’ parameters directly and the nonlinear method was the most accurate. The reprojection refinement, which works on the corner measurements and estimates the target’s pose together with $X$, has the smallest errors at almost every spread, and its translation error stays below 0.2 mm even at the smallest spread, where the closed-form methods are off by 8 mm or more; its rotation error there is no better than that of Tsai–Lenz. In our reading, it gains because it fits each corner measurement once, under the pixel noise the images actually carry, while the motion-based methods treat 45 pairwise motions built from only ten poses as independent equations. The example holds the robot’s kinematics exact; with real robots, Horaud and Dornaika found the robot’s errors larger than the camera’s, and a reprojection refinement that trusts the robot’s poses inherits those errors.

Evaluation Metrics

In simulation the true transform is known, and errors are measured directly. Horaud and Dornaika (1995, §6) report the root-mean-square difference between estimated and nominal rotation matrices over 1,000 trials and the root-mean-square translation error relative to the nominal translation length, as functions of the noise level and of the number of motions; the example uses the geodesic rotation angle and the translation distance, whose definitions the camera pose article discusses.

With real data the true transform is unknown, so methods are judged by consistency. Horaud and Dornaika’s tables give, for each method, the sum of squares of the absolute rotation error and of the relative translation error (§7), that is, how well the estimated transform satisfies the measured motion equations. Natural further checks, which are our suggestions rather than practice reported in the sources read, are the reprojection error of target corners predicted through the chain $T_{bg}^{(i)} X\, T_{ct}^{(i)}$, residuals on stations held out from the fit, and the spread of the motion axes, since small residuals from weak motions do not imply an accurate transform.

Applications

• Grasping and manipulation. A camera on the wrist locates an object, and hand-eye calibration carries its pose into the robot’s frame, as in robotic manipulation and robotic grasping; robot perception lists it among the transforms a robot must know and keep current.
• Industrial robot guidance. In factory cells the calibrated transform turns part poses measured in images into robot coordinates, as the industrial machine vision article describes.
• Visual servoing. The constant camera-to-end-effector transform enters the relation between camera velocity and joint velocity in eye-in-hand control; visual servoing describes when a coarse value suffices.
• Moving-camera measurement. Horaud and Dornaika (1995, §1) list inspection of complex 3D parts and 3D reconstruction with a camera moved by a robot, where known camera poses come from the robot’s kinematics through $X$.
• Multi-sensor platforms. The same estimation links cameras, LiDARs and IMUs on vehicles, as outlined under Variants.

Open Problems

• Calibration in place. Calibrating without a jig, with small motions in a limited workspace, is the goal of Andreff et al. (2001); doing it reliably depends on tracking enough points through the motion, which their conclusion names as a difficulty.
• Joint robot and hand-eye calibration. Robot kinematic errors dominate in practice, and Horaud and Dornaika (1995, §8) name calibrating the robot together with the hand-eye transform, and dropping the assumption of fixed intrinsics, as further work.
• Choosing motions. How to choose stations that make the problem well conditioned, within a robot’s reachable and collision-free workspace and with the target in view, follows from the conditions above only qualitatively; Andreff et al. (2001) close by naming the motions that give the highest numerical accuracy as a question never answered.
• Time and continuous motion. For sensors that are not triggered together, such as cameras and IMUs, the time offset must be estimated with the transform, as Kalibr does, and the motion must excite every axis.