§ 7. Getting practical

The goal is to estimate the camera matrix, given a number of correspondences mj,Mj⁢j=1⁢…⁢n

The geometric error associated to a camera estimate P⌃ is the distance between the measured image point mj and the re-projected point P⌃i⁢Mj:

where d() is the Euclidean distance between the homogeneous points.

The DLT solution is used as a starting point for the iterative minimization (e.g. Gauss-Newton)

The goal is to estimate the 3-D coordinates of a point M, given its projection mi and the camera matrix Pi for every view i=1⁢…⁢m.

The geometric error associated to a point estimate M⌃ in the i-th view is the distance between the measured image point mi and the re-projected point Pi⁢M⌃:

where d() is the Euclidean distance between the homogeneous points.

The linear solution is used as a starting point for the iterative minimization (e.g. Gauss-Newton).

The goal is to estimate F given a a number of point correspondences mℓi↔mri.

The geometric error associated to an estimate F⌃ is given by the distance of conjugate points from conjugate lines (note the symmetry):

where d() here is the Euclidean distance between a line and a point (in homogeneous coordinates).

The eight-point solution is used as a starting point for the iterative minimization (e.g. Gauss-Newton).

Note that F must be suitably parametrized, as it has only seven d.o.f.

The goal is to estimate H given a a number of point correspondences mℓi↔mri.

The geometric error associated to an estimate H⌃ is given by the symmetric distance between a point and its transformed conjugate:

where d() is the Euclidean distance between the homogeneous points. This also called the symmetric transfer error.

The linear solution is used as a starting point for the iterative minimization (e.g. Gauss-Newton).

If measurements are noisy, the projection equation will not be satisfied exactly by the camera matrices and structure computed in Sec. 5.3.2.

We wish to minimize the image distance between the re-projected point P⌃i⁢M⌃j and measured image points mij for every view in which the 3-D point appears:

where d() is the Euclidean distance between the homogeneous points.

As m and n increase, this becomes a very large minimization problem.

A solution is to alternate minimizing the re-projection error by varying P⌃i with minimizing the re-projection error by varying M⌃j.

See (Triggs et al.,2000) for a review and a more detailed discussion on bundle adjustment.

Least squares:

where θ are the regression coefficient (what is being estimated) and ri is the residual. M-estimators are based on the idea of replacing the squared residuals by another function of the residual, yielding

ρ is a symmetric function with a unique minimum at zero that grows sub-quadratically, called loss function.

Differentiating with respect to θ yields:

The M-estimate is obtained by solving this system of non-linear equations.

Given a model that requires a minimum of p data points to instantiate its free parameters θ, and a set of data points S containing outliers:

1. Randomly select a subset of p points of S and instantiate the model from this subset

2. Determine the set Si of data points that are within an error tolerance t of the model. Si is the consensus set of the sample.

3. If the size of Si is greater than a threshold T, re-estimate the model (possibly using least-squares) using Si (the set of inliers) and terminate.

4. If the size of Si is less than T, repeat from step 1.

5. Terminate after N trials and choose the largest consensus set found so far.

Three parameters need to be specified: t,T and N.

Both T and N are linked to the (unknown) fraction of outliers ϵ.

N should be large enough to have a high probability of selecting at least one sample containing all inliers. The probability to randomly select p inliers in N trials is:

By requiring that P must be near 1, N can be solved for given values of p and ϵ.

T should be equal to the expected number of inliers, which is given (in fraction) by 1-ϵ.

At each iteration, the largest consensus set found so fare gives a lower bound on the fraction of inliers, or, equivalently, an upper bound on the number of outliers. This can be used to adaptively adjust the number of trials N.

t is determined empirically, but in some cases it can be related to the probability that a point under the threshold is actually an inlier (Hartley and Zisserman,2003).

As pointed out by Stewart (1999), RANSAC can be viewed as a particular M-estimator.

The objective function that RANSAC maximizes is the number of data points having absolute residuals smaller that a predefined value t. This may be seen a minimising a binary loss function that is zero for small (absolute) residuals, and 1 for large absolute residuals, with a discontinuity at t.

By virtue of the prespecified inlier band, RANSAC can fit a model to data corrupted by substantially more than half outliers.

Another popular robust estimator is the Least Median of Squares. It is defined by:

It can tolerate up to 50% of outliers, as up to half of the data point can be arbitrarily far from the “true” estimate without changing the objective function value.

Since the median is not differentiable, a random sampling strategy similar to RANSAC is adopted. Instead of using the consensus, each sample of size p is scored by the median of the residuals of all the data points. The model with the least median (lowest score) is chosen.

A final weighted least-squares fitting is used.

With respect to RANSAC, LMedS can tolerate “only” 50% of outliers, but requires no setting of thresholds.

Following the development of Sec. 4.6 we know that for a camera P=K[R|t] the homography between a world plane at z=0 and the image is

where ri are the column of R.

Suppose that H is computed from correspondences between four or more known world points and their images, then some constraints can be obtained on the intrinsic parameters, thanks to the fact that the columns of R are orthonormal.

Writing H=h1,h2,h3, from the previous equation we derive:

where λ is an unknown scale factor.

The orthogonality r1T⁢r2=0 gives

or, equivalently (remember that ω=(KKT)-1)

Likewise, the condition on the norm r1T⁢r1=r2T⁢r2 gives

Introducing the Kronecker product as usual, we rewrite these two equations as:

As ω is a 3×3 symmetric matrix, its unique elements (the unknowns) are only six. This fact can be neatly taken into account using the v⁢e⁢c⁢h operator. The above equations are equivalent to ⁶

From a set of n images, we obtain a 2×n×6 coefficient matrix A by stacking up two equations for each image. The solution is the 1-dimensional right null-space of A.

At least five equations are needed (2⁢12 images). In practice, for a good calibration, one should use around 12 views).

K is obtained from the Cholesky factorization of ω, then R and t are recovered from:

Because of noise, the matrix R is not guaranteed to be orthogonal, hence we need to recover the closest orthogonal matrix.

Let R=Q⁢S be the polar decomposition of R. Then Q is the closest possible orthogonal matrix to R in Frobenius norm.

In this way we have obtained the camera matrix P by minimizing an algebraic distance which is not geometrically meaningful.

It is advisable to refine it with a (non-linear) minimization of a geometric error:

where P⌃i=K⌃[R⌃i|t⌃i] and the rotation has to be suitably parametrized with three parameters (see Rodrigues formula).

The linear solution is used as a starting point for the iterative minimization (e.g. Gauss-Newton).

The more relevant effect is the radial distortion, which is modeled as a non-linear transformation from the ideal (undistorted) pixel coordinates u,v to the observed (distorted) pixel coordinates u⌃,v⌃:

where rd2=u-u0au2+v-v0av2 and u0,v0 are the coordinates of the image centre.

This distortion makes a rectangle to appear as a “barrel” or a “cushion” (as in Fig. 19) depending whether the coefficient k1 is lower or greater than one.

Let us assume that the pinhole model is calibrated. The point m=u,v projected according to the pinhole model (undistorted) do not coincide with the measured points m⌃=u⌃,v⌃ because of the radial distortion.

We wish to recover k1 from Eq. (137). Each point gives two equation:

hence a least squares solution for k1 is readily obtained from n>1 points.

When calibrating a camera we are required to simultaneously estimate both the pinhole model's parameters and the radial distortion coefficient.

The pinhole calibration we have described so far assumed no radial distortion, and the radial distortion calibration assumed a calibrated pinhole camera.

The solution (a very common one in similar cases) is to alternate between the two estimation until convergence.

Namely: start assuming k=0, calibrate the pinhole model, then use that model to compute radial distortion. Once k1 is estimated, refine the pinhole model by solving Eq. (136) with the radial distortion in the projection, and continue until the image error is small enough.