Caliber: Camera Localization and Calibration Using Rigidity Constraints

Abstract

This article presents a camera calibration system, Caliber, and the underlying pose estimation problem it solves, which we call sensor localization with rigidity (SL-R). SL-R is a constraint-satisfaction-like problem that finds a set of poses satisfying certain constraints. These constraints include not only relative pose constraints such as those found in SLAM and motion estimation problems, but also rigidity constraints: the notion of objects that are rigidly attached to each other so that their relative pose is fixed over time even if that pose is not known a priori. We show that SL-R is NP-hard, but give an inference-based algorithm that works well in practice. SL-R enables Caliber, a tool to calibrate systems of cameras connected by rigid or actuated links, using image observations and information about known motions of the system. The user provides a model of the system in the form of a kinematic tree, and Caliber uses our SL-R algorithm to generate an estimate for the rigidity constraints, then performs nonlinear optimization to produce a solution that is locally least-squares optimal in terms of reprojection error. In this way, Caliber is able to calibrate a variety of setups that would have previously required special-purpose code to calibrate. We demonstrate Caliber in a number of different scenarios using both synthetic and experimental data.

Appendices

Appendix 1: Systems of AX = XB Equations

In this section we cover solutions to systems of $AX=XB$ equations of rigid transformations. These solutions are not new in an of themselves, but they are a subroutine that allows our algorithm to solve non-trivial instances of SL-R.

1.1 Overview

The problem of solving a system of $AX=XB$ equations of rigid transformations is defined as follows: given a set of pairs of 4x4 rigid transformation matrices $A_i$, $B_i$, determine a 4x4 rigid transformation X such that for all i, $A_i X = X B_i$. Systems of equations of this sort naturally arise from situations such as hand-eye calibration and non-overlapping views from rigidly attached cameras.

1.2 Separation of Rotation and Translation

We opt to separate the $AX=XB$ equation on full rigid transformations into its translational and rotational parts. This gives

$$\begin{aligned}&R_A R_X = R_X R_B \end{aligned}$$

(17)

$$\begin{aligned}&\left( R_A - I \right) \mathbf{t}_x = R_X \mathbf{t}_b - \mathbf{t}_a \end{aligned}$$

(18)

1.3 Screw Theory

It is also useful to observe that these equations describe a similarity transformation $B = X^{-1}AX$. Furthermore, (Chen 1991) gives a method of characterizing $AX=XB$ systems using screw theory. Briefly, screws are a way of representing rigid transformations. They consist of a screw axis, which is a line in space, a translation magnitude along that axis, and a rotation magnitude about that axis. Screws have the property that the translation and rotation magnitudes are invariant with rigid transformations, whereas the screw axis transforms with rigid transformations. Therefore, if we have one or more $AX=XB$ equations, we can think of this problem as one where we attempt to find the transformation X that takes the screw axes given by $A_1, A_2, \ldots$ to those given by $B_1, B_2, \ldots$

1.4 Cases with Rotation

Nonparallel (skew or intersecting) screw axes (0 undetermined dimensions) (Park and Martin 1994) gives a method to compute the unique answer in the case where the screw axes are nonparallel and the rotation magnitudes are not near the identity or a $180^\circ$ rotation (so the screw axes are well-defined). If the problem is overconstrained by having more than two $AX=XB$ equations, they give a least-squares solution. In particular, they compute the rotational matrix that transforms a $3 \times n$ matrix of screw axis direction vectors $\mathbf{\alpha }$ (corresponding to the As) to another one $\mathbf{\beta }$ (corresponding to the Bs) as

$$\begin{aligned} M&= \mathbf{\beta } \mathbf{\alpha }^T \end{aligned}$$

(19)

$$\begin{aligned} R&= \left( M^T M \right) ^{-\frac{1}{2}} M^T \end{aligned}$$

(20)

After the rotation is solved for, they find a least-squares solution for the translation. We use this solution directly for the fully-determined case, and our solutions for the other cases are based on this method.

Parallel screw axes (1 undetermined dimension) (Chen 1991) shows that in the case where the screw axes are parallel but not coincident, the rotation can be uniquely determined, but the translation has one undetermined dimension. For two such pairs of rigid transformations $A_1, B_1$ and $A_2, B_2$, we can construct three orthogonal directions using the (mutual) screw axis direction, the direction separating the two screw axes, and their cross product. We can then solve for the rotation as before.

Since $R_A$ always has at least one eigenvalue of unity, $\left( R_A - I \right)$ in the translation Eq. 18 has a one-dimensional null space, and the solution space has one undetermined dimension. Namely, this is the common screw axis direction. Since all the screws are parallel, this null space is the same for all of the equations. In this case, we attempt to keep the translation as close to the origin as possible, which we do by adding another row to the matrix $\mathbf{r}_a^T \mathbf{t}_x = 0$. Finally, we stack the equations produced by each pair of $AX=XB$ equations to produce a single least-squares solution for $\mathbf{t}_x$.

Note that this gives an alternative method for solving the skew screw axis case: choose one screw axis as one direction, the direction along the shortest line between two screw axes as another, and their cross product as the third. Since in this case the null space of $\left( R_A - I \right)$ for different equations are not the same, the translation has a unique solution.

Coincident screw axes (2 undetermined dimensions) In the coincident screw axis case, only one direction for the As and for the Bs can be determined. The requirement that $R_X$ rotate one direction onto another lets us solve for two rotational dimensions, but the third dimension is undetermined. In this case we choose the rotation axis to be along the (average) cross product of the screw axes of each $A_i, B_i$ pair, which minimizes the rotation angle. We then solve for the translation as per the parallel axis case.

1.5 Cases Without Rotation

In some cases we may not have rotations at all; here the screw axis has direction equal to the translation but its offset perpendicular to the axis direction is undefined. As such, we may be able to construct axis directions, but the origin position of the axes are undetermined. In terms of the separated equations, the rotation Eq. 17 reduces to $R_X = R_X$, which gives us no information. The left side of the translation equation 18 is zero, which means that $\mathbf{t}_x$ cannot be determined at all, but there may be some information about $R_X$.

Nonparallel translation (3 undetermined dimensions) If the translations are nonparallel, we can use two non-parallel translations and their cross product to construct three mutually orthogonal axis directions. This allows us to determine the rotation uniquely, as in the noncoincident screw axis cases. As noted above, the translation is completely undetermined.

Parallel translation case (4 undetermined dimensions) Similarly, the parallel translation case is analogous to the coincident screw axis case in that the rotation can be determined up to one dimension. As in the coincident screw axis case, the requirement that $R_X$ rotate one direction onto another lets us solve for two rotational dimensions, but the third dimension is undetermined. We can use the same method as the coincident screw axis case to determine a feasible rotation. As in the nonparallel translation case, the translation is completely undetermined.

No information case (6 undetermined dimensions)

If we have no information at all (no equations, or all As and Bs are the identity), then we simply guess X to be the identity.

1.6 Omitted Cases

(Chen 1991) addresses the cases where some or all of the rotations are $180^\circ$. We have not implemented these cases so far because they are rare in practical use. Typically, the As and Bs in the equations result from the relationship between two positions of a joint; this means skew screw axes for a ball joint, coincident screw axes for a revolute or cylindrical joint, and parallel translation for a prismatic joint. Therefore, it is rare to find cases with only $180^\circ$ rotations or mixed types of equations.

Appendix 2: Complexity

In this section we show that SL-R is NP-hard. Our proof is of weak NP-hardness due to the assumption that a transformation matrix may be specified in constant space. Furthermore, we show that the uniqueness variant of SL-R—that is, given an instance of SL-R for which a solution is guaranteed to exist, does there exist more than one solution?— is also NP-hard. Therefore, if P $\ne$ NP, there exists no polynomial-time algorithm that can solve SL-R, nor determine whether an instance of SL-R is ambiguous.

Our proof has two parts. First, we construct a gadget that forces any solution to make a binary choice. Then, we use this gadget to construct a reduction from the NP-complete (unique) partition problem.

1.1 Binary Choice via XAX = B

We begin by inventing a SL-R gadget that forces the solution to make a binary choice. This gadget relies on equations of the form $XAX = B$, where all three quantities are rigid transformations and X is to be solved for.

In terms of SL-R such an equation consists of four absolute poses $i, j, k, \ell$ with relative pose constraints $P_{ik}$ and $P_{j\ell }$, and a rigidity constraint X on the relative poses $P_{ji}$ and $P_{k\ell }$.

Pictorially, such an equation looks like this:

(21)

Mathematical properties Inspiration taken from (Heath 2008).

Suppose we represent each relative pose as a 4x4 rigid transformation matrix

$$\begin{aligned} \begin{bmatrix} R&\mathbf{t} \\ 0&1 \end{bmatrix} \end{aligned}$$

(22)

where R is a 3x3 rotation matrix, and $\mathbf{t}$ is a 3-vector representing translation.

From a single equation we have

$$\begin{aligned}&XAX = B \end{aligned}$$

(23)

$$\begin{aligned}&AXAX = AB \end{aligned}$$

(24)

$$\begin{aligned}&AX = \left( AB \right) ^\frac{1}{2} \end{aligned}$$

(25)

Suppose we consider only the rotational part. For general real matrices, the square root has many possible values; however, since our solution is restricted to rotation matrices, this limits the possibilities. If $R_A R_B = I$, then I and any $180^\circ$ rotation is a possible square root. Otherwise, $R_A R_B$ has a defined axis of rotation, and $\left( R_A R_B \right) ^\frac{1}{2}$ must share that same axis, since applying any rotation twice preserves the axis of rotation. The magnitude of the rotation must be half that of $R_A R_B$, plus an integer multiple of $\pi$. Since a difference of $2 \pi$ leaves the rotation matrix unchanged, this produces two possible values for $\left( R_A R_B \right) ^\frac{1}{2}$ and therefore two possible values for $R_X$.

From here, we examine the translational part given the rotation. Denote the rotational and translational parts by $R_X$ and $\mathbf{t}_x$ respectively. Then solving for the unknown $\mathbf{t}_x$ reduces to

$$\begin{aligned}&R_X R_A \mathbf{t}_x + R_X \mathbf{t}_a + \mathbf{t}_x = \mathbf{t}_b \end{aligned}$$

(26)

$$\begin{aligned}&\left( R_X R_A + I \right) \mathbf{t}_x = \mathbf{t}_b - R_X \mathbf{t}_a \end{aligned}$$

(27)

This has a unique solution as long as $R_X R_A + I$ is invertible. To see when this is the case, consider this product for a unit vector x:

$$\begin{aligned} x^{-1} \left( R_X R_A + I \right) x = x^{-1} R_X R_A x + 1 \end{aligned}$$

(28)

The right side is only nonpositive when $R_X R_A$ is a $180^\circ$ rotation. Apart from this case, $R_X R_A + I$ is positive-definite and therefore invertible. Thus, unless $R_X R_A$ is a $180^\circ$ rotation, $\mathbf{t}_x$ has a unique solution given $R_X$.

The gadget Suppose we have a cycle of rigidity constraints X and relative pose constraints A and B that produces an $XAX = B$ equation. Let the rotational parts of A and B equal

$$\begin{aligned} R_A = R_B = \begin{bmatrix} 0&\quad -1&\quad 0 \\ 1&\quad 0&\quad 0 \\ 0&\quad 0&\quad 1 \end{bmatrix} \end{aligned}$$

(29)

From this we have

$$\begin{aligned} R_X&= R_A^{-1} \left( R_A R_B \right) ^{\frac{1}{2}} \in \left\{ I, I' \right\} \end{aligned}$$

(30)

$$\begin{aligned} I&= \begin{bmatrix} 1&\quad 0&\quad 0 \\ 0&\quad 1&\quad 0 \\ 0&\quad 0&\quad 1 \end{bmatrix} \end{aligned}$$

(31)

$$\begin{aligned} I'&= \begin{bmatrix} -1&\quad 0&\quad 0 \\ 0&\quad -1&\quad 0 \\ 0&\quad 0&\quad 1 \end{bmatrix} \end{aligned}$$

(32)

Now let us examine the translational parts of A, B, and X. Suppose

$$\begin{aligned} \mathbf{t}_a&= \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \end{aligned}$$

(33)

$$\begin{aligned} \mathbf{t}_b&= \begin{bmatrix} b \\ b \\ 0 \end{bmatrix} \end{aligned}$$

(34)

for some constant b. In the case where $R_X = I$, we have

$$\begin{aligned}&\left( R_X R_A + I \right) \mathbf{t}_x = \mathbf{t}_b - R_X \mathbf{t}_a \end{aligned}$$

(35)

$$\begin{aligned}&\left( R_A + I \right) \mathbf{t}_x = \mathbf{t}_b \end{aligned}$$

(36)

$$\begin{aligned}&\begin{bmatrix} 1&\quad -1&\quad 0 \\ 1&\quad 1&\quad 0 \\ 0&\quad 0&\quad 2 \end{bmatrix} \mathbf{t}_x = \begin{bmatrix} b \\ b \\ 0 \end{bmatrix} \end{aligned}$$

(37)

$$\begin{aligned}&\mathbf{t}_x = \begin{bmatrix} b \\ 0 \\ 0 \end{bmatrix} \end{aligned}$$

(38)

In the case where $R_X = I'$, we have

$$\begin{aligned}&\left( R_X R_A + I \right) \mathbf{t}_x = \mathbf{t}_b \end{aligned}$$

(39)

$$\begin{aligned}&\begin{bmatrix} 1&\quad 1&\quad 0 \\ -1&\quad 1&\quad 0 \\ 0&0&2 \end{bmatrix} \mathbf{t}_x = \begin{bmatrix} b \\ b \\ 0 \end{bmatrix} \end{aligned}$$

(40)

$$\begin{aligned}&\mathbf{t}_x = \begin{bmatrix} 0 \\ b \\ 0 \end{bmatrix} \end{aligned}$$

(41)

Therefore, in our gadget, X either represents a step “forwards” (+x direction) of size b, or a step to the “left” (+y direction) of size b followed by an about-face. If $b = 0$, then the choice is between standing still and an about-face. A pictorial representation is shown in Fig. 10. Also note that if b is an integer, all arithmetic here is integer, which shows that the NP-hardness of SL-R is not (solely) due to real arithmetic.

Fig. 10

Two-fold ambiguity for $XAX = B$ illustrated using a camera, a calibration target, and a revolute joint. A is a pose measurement taken by a revolute joint, B is a pose measurement taken by a camera and a checkerboard, and the two Xs are known to be equal. There are two distinct possible values of X that are consistent with the measurements. On the left, X has no rotational component, while on the right, X has a $180^\circ$ rotation

1.2 Reduction from the (Unique) Partition Problem

With this gadget in hand, we now show a polynomial reduction from the (unique) partition problem to the (unique) sensor localization with rigidity problem.

The partition problem Recall the partition problem (PP). The definition of the problem is as follows: given a set of positive integers $S = \left\{ x_1, x_2, \ldots x_n \right\}$, determine if there exists a subset of the integers $S_1$ such that

$$\begin{aligned} \sum _{x \in S_1} x = \sum _{x \not \in S_1} x \end{aligned}$$

(42)

The standard form of the problem is well-known to be NP-complete; (Dieudonné et al. 2010) shows that, given an instance of this problem for which it is guaranteed that at least one solution exists, the problem of determining whether there exists more than one solution (counting $S - S_1$ to be the same solution as $S_1$) is also NP-complete. They call this version the unique partition problem (UPP).

We likewise define the uniqueness version (USL-R) of SL-R to be: given an instance for SL-R which a solution is guaranteed to exist, does there exist more than one solution?

Integer representation We represent the integers of PP in SL-R as follows. For each $x_i \in S$, we introduce two rigidity constraints $X_i$ and $X_i'$. For both of these we construct an $XAX = B$ cycle of the form described above. For the one with $X_i$, we set $b = x_i$; for the one with $X_i'$, we set $b = 0$. This produces four options for the product $X_i X_i'$, corresponding to our choice of picking $R_X = I$ or $R_X = I'$ for each rotational part:

• 1. A step forward of size $x_i$ (choose $R_{X_i} = I$ and $R_{X_i'} = I$).

• 2. A step forward of size $x_i$ followed by an about-face (choose $R_{X_i} = I$ and $R_{X_i'} = I'$).

• 3. A step left of size $x_i$ (choose $R_{X_i} = I'$ and $R_{X_i'} = I'$).

• 4. A step left of size $x_i$ followed by an about-face (choose $R_{X_i} = I'$ and $R_{X_i'} = I$).

We construct each such cycle so that they share a common “root” absolute pose, so these are the only choices to be made here. Since we do constant work and produce a constant number of poses, relative pose constraints, and rigidity constraints per element of S, this part of the reduction is linear in time and problem size.

A Manhattan walk With our integer representation in place, we introduce one more cycle, which we will treat as two paths from the root to a destination. On one side we constrain the absolute pose of the the destination to equal

$$\begin{aligned}&P_\mathrm{dst} = \begin{bmatrix} 1&\quad 0&\quad 0&\quad \frac{1}{2} s \\ 0&\quad 1&\quad 0&\quad \frac{1}{2} s \\ 0&\quad 0&\quad 1&\quad 0 \\ 0&\quad 0&\quad 0&\quad 1 \end{bmatrix} \end{aligned}$$

(43)

$$\begin{aligned}&\text {where} \ s = \sum _i^n x_i \end{aligned}$$

(44)

On the other side we have a sequence of 2n poses with relative poses labeled such that they form the product

$$\begin{aligned} P_\mathrm{dst} = \prod _i^n X_i X_i' = X_1 X_1' X_2 X_2' \ldots X_n X_n' \end{aligned}$$

(45)

This part of the problem reduction is also linear in time and problem size.

The intuition is that we have constructed a problem where the goal is to reach a destination $\frac{s}{2}$ north and $\frac{s}{2}$ east of the starting point, using a set of steps of size $x_1, x_2, \ldots x_n$.

At each two-pose step towards the destination $X_i X_i'$ a solution must take one of the four options enumerated earlier. This constrains the translations to lie only on cardinal directions. Since the Manhattan distance to our target is s and there is only have a total of s distance in translations in the product of rigid transformations, any solution must step towards the destination in every step. The solution cannot do an about-face before the last step, since it would have to step away from the destination on the next step, and the solution cannot do an about-face on the last step, since it would be facing the wrong direction at the end. Therefore any solution must consist of only options 1 and 3, with forward steps being to the “east” and left steps to the “north”. Since the east steps total the same distance as the north steps in any solution, this forms a one-to-one correspondence to a solution to PP—the indices of the steps in which the solution moved east correspond to the indices of the integers in one half of the partition. To avoid the duplicate solution where $S_1$ is swapped with $S - S_1$, we require the first step to be east by adding a relative pose constraint.

An example is shown in Fig. 11.

Fig. 11

Pictorial example of a PP reduction to SL-R and its solution. We encode PP as a Manhattan walk where steps to the east represent one partition and steps to the north the other partition. We enforce the choice between these two directions using gadgets of the form $XAX = B$

Therefore PP reduces to SL-R, UPP reduces to USL-R, and both SL-R and USL-R are NP-hard.