James Gregson's Website

Mar 09, 2019

Tsai Camera Calibration

Projection of known 3D points into a pinhole camera can be modeled by the following equation:

\begin{align*} \begin{bmatrix} \hat{u}_c \\ \hat{v}_c \\ \hat{w}_c \end{bmatrix} = \begin{bmatrix} f & 0 & u_0 \\ 0 & f & v_0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} r_{xx} & r_{xy} & r_{xz} & t_x \\ r_{yx} & r_{yy} & r_{yz} & t_y \\ r_{zx} & r_{zy} & r_{zz} & t_z \end{bmatrix} \begin{bmatrix} x_w \\ y_w \\ z_w \\ 1 \end{bmatrix} \\ = K E x \end{align*}

The matrix \(K\) contains the camera intrinsics and the matrix \(E\) contains the camera extrinsics. The final image coordinates are obtained by perspective division: \(u_c = \hat{u}_c/\hat{w}_c\) and \(v_c = \hat{v}_c/\hat{w}_c\).

Estimating \(K\) and \(E\) jointly can be difficult because the equations are non-linear and may converge slowly or to a poor solution. Tsai's algorithm provides a way to obtain initial estimates of the parameters for subsequent refinement. It requires a known calibration target that can be either planar or three-dimensional.

The key step is in estimating the rotation and x-y translation of the extrinsics \(E\) in such a way that the unknown focal length is eliminated. Once the rotation is known, the focal-length and z translation can be computed.

This process requires the optical center \([u_0, v_0]\) which is also unknown. However, given reasonable quality optics this is often near the image center so the center can be used as an initial estimate.

This is the starting point for the estimation.

Estimating (most of) the Extrinsics

This is divided into two cases, one for planar targets and one for 3D targets.

Case 1: 3D Non-planar Target:

For calibration targets with 8+ points that span 3 dimensions. Points in the camera frame prior to projection are:

\begin{equation*} \begin{bmatrix} x_c \\ y_c \\ z_c \end{bmatrix} = \begin{bmatrix} r_{xx} x_w + r_{xy} y_w + r_{xz} z_w + t_x \\ r_{yx} x_w + r_{yy} y_w + r_{yz} z_w + t_y \\ r_{zx} x_w + r_{zy} y_w + r_{zz} z_w + t_z \end{bmatrix} \end{equation*}

The image points relative to the optical center are then:

\begin{equation*} \begin{bmatrix} u_c - u_0 \\ v_c - v_0 \end{bmatrix} = \begin{bmatrix} u_c' \\ v_c' \end{bmatrix} = \begin{bmatrix} (f x_c + u_0 z_c)/z_c - u_0 \\ (f y_c + v_0 z_c)/z_c - v_0 \end{bmatrix} = \begin{bmatrix} f x_c/z_c \\ f y_c/z_c \end{bmatrix} \end{equation*}

The ratio of these two quantities gives the equation:

\begin{equation*} \frac{u_c'}{v_c'} = \frac{x_c}{y_c} \end{equation*}

or, after cross-multiplying and collecting terms to one side:

\begin{align*} v_c' x_c - u_c' y_c = 0 \\ [ v_c' x_w, v_c' y_w, v_c' z_w, -u_c' x_w, -u_c' y_w, -u_c' z_w, 1, -1] \begin{bmatrix} r_{xx} \\ r_{xy} \\ r_{xz} \\ r_{yx} \\ r_{yy} \\ r_{yz} \\ t_x \\ t_y \end{bmatrix} = 0 \end{align*}

Each point produces a single equation like the above which can be stacked to form a homogeneous system. With a zero right-hand-side, scalar multiples of the solution still satisfy the equations and a trivial all-zero solution exists. Calling the system matrix \(M\), the solution to this equation system is the eigenvector associated with the smallest eigenvalue of \(M^T M\). The computed solution will be proportional to the desired solution. It is also possible to fix one value, e.g. \(t_y=1\) and solve an inhomogeneous system but this may give problems if \(t_y \neq 0\) for the camera configuration.

The recovered \(r_x = [ r_{xx}, r_{xy}, r_{xz} ]\) and \(r_x = [ r_{yx}, r_{yy}, r_{yz} ]\) should be unit length. The scale of the solution can be recovered as the mean length of these vectors and used to inversely scale the translations \(t_x, t_y\).

Case 2: 2D Planar Targets:

For 2D targets, arbitrarily set \(z_w=0\), leading to:

\begin{equation*} \begin{bmatrix} x_c \\ y_c \\ z_c \end{bmatrix} = \begin{bmatrix} r_{xx} x_w + r_{xy} y_w + t_x \\ r_{yx} x_w + r_{yy} y_w + t_y \\ r_{zx} x_w + r_{zy} y_w + t_z \end{bmatrix} \end{equation*}

Use the same process as before, a homogeneous system can be formed. It is identical to the previous one, except the \(r_{xz}\) and \(r_{yz}\) terms are dropped since \(z_w=0\). Each equation ends up being:

\begin{equation*} [ v_c' x_w, v_c' y_w, -u_c' x_w, -u_c' y_w, 1, -1] \begin{bmatrix} r_{xx} \\ r_{xy} \\ r_{yx} \\ r_{yy} \\ t_x \\ t_y \end{bmatrix} = 0 \end{equation*}

Solving this gives \(r_{xx}, r_{xy}, r_{zx}, r_{zy}, t_x, t_y\) up to a scalar \(k\). It does not provide \(r_{xz}, r_{yz}\). In addition, we know that the rows of the rotation submatrix are unit-norm and orthogonal, giving three equations:

\begin{align*} r_{xx}^2 + r_{xy}^2 + r_{xz}^2 = k^2 \\ r_{yz}^2 + r_{yy}^2 + r_{yz}^2 = k^2 \\ r_{xx} r_{yx} + r_{xy} r_{yy} + r_{xz} r_{yz} = 0 \end{align*}

After some algebraic manipulations the factor \(k^2\) can be found as:

\begin{equation*} k^2 = \frac{1}{2}\left( r_{xx}^2 + r_{xy}^2 + r_{yx}^2 + r_{yy}^2 + \sqrt{\left( (r_{xx} - r_{yy})^2 + (r_{xy} + r_{yx})^2\right)\left( (r_{xx}+r_{yy})^2 + (r_{xy} - r_{yx})^2 \right) } \right) \end{equation*}

This allows \(r_{xz}, r_{yz}\) to be found as:

\begin{align*} r_{xz} = \pm\sqrt{ k^2 - r_{xx}^2 - r_{xy}^2 } \\ r_{yz} = \pm\sqrt{ k^2 - r_{yz}^2 - r_{yy}^2 } \end{align*}

Estimate Remaining Intrinsics and Focal Length

Estimating the focal length and z-translation follows a similar process. Using the definitions above we have:

\begin{align*} u_c' = f \frac{x_c}{z_c} \\ v_c' = f \frac{y_c}{z_c} \end{align*}

Cross-multiplying and collecting gives:

\begin{align*} (r_{xx} x_w + r_{xy} y_w + r_{xz} z_w + t_x) f - u_c' t_z = u_c' (r_{zx} x_w + r_{zy} y_w + r_{zz} z_w) \\ (r_{yx} x_w + r_{yy} y_w + r_{yz} z_w + t_y) f - v_c' t_z = v_c' (r_{zx} x_w + r_{zy} y_w + r_{zz} z_w) \end{align*}

Solving this system gives the desired \(f, t_x\) values.