## Vanishing points in presence of noise

Most self-calibration algorithms require a prior knowledge of the camera calibration matrix $K$; as an instance, you need it to normalize the image points as $\tilde{\mathbf{x}}=K^{-1} \mathbf{x}$ and therefore fit the essential matrix $E$.

With most commercial cameras it is safe to assume that the pixels are square and have no slant; in this case, $K= \begin{bmatrix} -fk & 0 & X_c \\ 0 & -fk & Y_c \\ 0 & 0 & 1 \end{bmatrix}$, where

• $X_c, Y_c$ are the coordinates of the principal point in pixel units; it is usually an acceptable first approximation to assume that the principal point is at the centre of the sensor;
• $k$ is the number of pixels per mm; it is usually safe to take the nominal value specified by the sensor vendor;
• $f$ is the focal length in mm. Taking the nominal value specified by the lens vendor could not be sufficient; even worse if, working with a zoom lens, one tries to read the focus scale.

So it would be nice to obtain at least $f$, or better the whole calibration matrix, from the images. The calibration matrix can be computed from the image of the absolute conic which, in the above conditions, is $\omega = \begin{bmatrix} \omega_1 & 0 & \omega_2 \\ 0 & \omega_1 & \omega_3 \\ \omega_2 & \omega_3 & \omega_4 \end{bmatrix}$; it is related to $K$ by $\omega \propto K^{-T}K^{-1}= \begin{bmatrix} \frac{1}{f^2 k^2} & 0 & -\frac{X_c}{f^2 k^2} \\ 0 & \frac{1}{f^2 k^2} & -\frac{Y_c}{f^2 k^2} \\ -\frac{X_c}{f^2 k^2} & -\frac{Y_c}{f^2 k^2} & 1+\frac{X_c^2+Y_c^2}{f^2 k^2} \end{bmatrix}$.

$\omega$, in turn, can be obtained from three couples of vanishing points $\mathbf{v}_1, \mathbf{v}_1'$; $\mathbf{v}_2, \mathbf{v}_2'$$\mathbf{v}_3, \mathbf{v}_3'$ corresponding to pencils of lines which are orthogonal in the object space; each couple contributes a constraint $\mathbf{v}^T \omega \mathbf{v}'=0$. Three constraints are sufficient because $\omega$ has four distinct non-null entries but, being the matrix of a conic section, it is homogeneous (defined up to a scale factor).

In an ideal world, things are so simple. All the lines of the pencil intersect exactly in the same point, whose coordinates are readily computed and can be used to obtain $\omega$.

In the real world, nothing is perfect and things are not so simple. Straight lines must be obtained by fitting them to points (edgels), which are noisy. As a result, each couple of lines in the pencil will intersect in a different point (though all the points should be very close).

Taking the centroid of the intersections may work, but it is possible to do better by realizing that we are facing a constrained minimum problem: we are looking for a pencil of straight lines, minimizing the sum of squared distances of each line from its own cloud of points, all lines constrained to intersect in the same point.

Let

• $N$ be the number of lines,
• $M_i, i \in 1 \cdots N$ be the number of points in the cloud of the i-th line,
• $x_{ij}, y_{ij}, j \in 1 \cdots M_i$ be the coordinates of the j-th point in the cloud of the i-th line,
• $\theta_i, d_i$ be the parameters of the polar representation of the i-th line,
• $x_I, y_I$ be the coordinates of the intersection;

then the problem can be written as

$min_{\{\theta_i\},\{d_i\}} \sum_{i=1}^N \sum_{j=1}^{M_i} (x_{ij} cos \theta_i + y_{ij} sin \theta_i - d_i)^2$

subject to

$x_I cos \theta_i + y_I sin \theta_i - d_i = 0, i \in 1 \cdots N$

The lagrangian is

$\Lambda (\{\theta_i\}, \{d_i\}, x_I, y_I)=\sum_{i=1}^N (\sum_{j=1}^{M_i} (x_{ij} cos \theta_i + y_{ij} sin \theta_i - d_i)^2 + \lambda_i (x_I cos \theta_i + y_I sin \theta_i - d_i))$

where $\lambda_i$ are the Lagrange multipliers.

This is a nonlinear problem and as such needs an initial approximation to start with. The unconstrained straight lines can be fitted eg with the method of Alciatore and Miranda, which is closed form and fast, and in this context possible outliers can be rejected eg with RANSAC; the intersection can be obtained as the centroid of all the intersections. It is a good idea to check that no intersection be too far from the centroid; this could be a symptom of something gone astray in the previous steps.

Unluckily, it is not possible to adapt to this problem the linear method I described here (reduction to the subspace of constraints) because the common intersection point is not known in advance.