Calculating camera calibration matrix with Scilab

Question

I'm not entirely sure whether this question belongs here or in DSP but I think this is the proper site.

I'm following these videos (first video, second video) to calibrate a camera for photogrammetry using Scilab. The method is pretty simple but I can't get proper results. I'm not sure whether I have misunderstood something in the calibration or if I'm using Scilab wrong.

The calibration method uses a calibration object of known geometry. I have 3d printed a cube with markers for that purpose. An image of the object is taken and a number of points on it are selected. The world coordinates x_w, y_w, z_w of those points are determined as well as pixel coordinates u and v (horizontal and vertical).

The coordinates of the points are put in a matrix A like this: $$ A=\begin{bmatrix} x_{w1} & y_{w1} & z_{w1} & 1 & 0 & 0 & 0 & 0 & -u_1x_{w1} & -u_1y_{w1} & -u_1z_{w1} & -u_1 \\ 0 & 0 & 0 & 0 & x_{w1} & y_{w1} & z_{w1} & 1 & -v_1x_{w1} & -v_1y_{w1} & -v_1z_{w1} & -v_1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ x_{wn} & y_{wn} & z_{wn} & 1 & 0 & 0 & 0 & 0 & -u_nx_{wn} & -u_ny_{wn} & -u_nz_{wn} & -u_n \\ 0 & 0 & 0 & 0 & x_{wn} & y_{wn} & z_{wn} & 1 & -v_nx_{wn} & -v_ny_{wn} & -v_nz_{wn} & -v_n \\ \end{bmatrix} $$

We want to solve projection matrix P which we can do by solving the eigenvalues and eigenvectors of A^TA. The projection matrix in vector form p is the eigenvector corresponding to the smallest eigenvalue. Then the projection matrix can be formed by rearranging the elements of p.

The calibration matrix is obtained by applying QR factorization to the first three columns of the projection matrix. The resulting upper triangle matrix is the calibration matrix K. It contains the values which we are interested in.

$$ K=\begin{bmatrix} f_x & 0 & o_x \\ 0 & f_y & o_y \\ 0 & 0 & 1 \\ \end{bmatrix} $$

The f_x and f_y are the effective focal lengths in x and y directions and should be about the same if the pixels are square. The o_x and o_y are the coordinates of the principal point which is where the optical axis intersects the image plane. It is usually near the center of the image so the coordinates should be in the order of hundreds or few thousands depending on the resolution of the image.

Here is the Scilab script I made to do the calculations. I have used only six points here which is way too little for accurate results but the algorithm should still give approximate results.

//Coordinates
//Pixels horizontal
u1=1196;
u2=1776;
u3=2553;
u4=1241;
u5=1753;
u6=2459;

//Pixels vertical
v1=1003;
v2=1246;
v3=1084;
v4=1859;
v5=2277;
v6=1947;

//Millimeters X
xw1=0;
xw2=0;
xw3=50;
xw4=0;
xw5=0;
xw6=50;

//Millimeters Y
yw1=50;
yw2=0;
yw3=0;
yw4=50;
yw5=0;
yw6=0;

//Millimeters Z
zw1=50;
zw2=50;
zw3=50;
zw4=0;
zw5=0;
zw6=0;

//The matrix
A=[xw1 yw1 zw1 1 0 0 0 0 -u1*xw1 -u1*yw1 -u1*zw1 -u1;
0 0 0 0 xw1 yw1 zw1 1 -v1*xw1 -v1*yw1 -v1*zw1 -v1;
xw2 yw2 zw2 1 0 0 0 0 -u2*xw2 -u2*yw2 -u2*zw2 -u2;
0 0 0 0 xw2 yw2 zw2 1 -v2*xw2 -v2*yw2 -v2*zw2 -v2;
xw3 yw3 zw3 1 0 0 0 0 -u3*xw3 -u3*yw3 -u3*zw3 -u3;
0 0 0 0 xw3 yw3 zw3 1 -v3*xw3 -v3*yw3 -v3*zw3 -v3;
xw4 yw4 zw4 1 0 0 0 0 -u4*xw4 -u4*yw4 -u4*zw4 -u4;
0 0 0 0 xw4 yw4 zw4 1 -v4*xw4 -v4*yw4 -v4*zw4 -v4;
xw5 yw5 zw5 1 0 0 0 0 -u5*xw5 -u5*yw5 -u5*zw5 -u5;
0 0 0 0 xw5 yw5 zw5 1 -v5*xw5 -v5*yw5 -v5*zw5 -v5;
xw6 yw6 zw6 1 0 0 0 0 -u6*xw6 -u6*yw6 -u6*zw6 -u6;
0 0 0 0 xw6 yw6 zw6 1 -v6*xw6 -v6*yw6 -v6*zw6 -v6];

//Calculate eigenvectors and eigenvalues
[evec, eval]=spec(A'*A);

//Select the eigenvector corresponding to the smallest eigenvalue
p=evec(:,1);

//Rearrange elements to form the projection matrix
P=[p(1) p(2) p(3) p(4);p(5) p(6) p(7) p(8);p(9) p(10) p(11) p(12)];

//Apply QR factorization to the first three columns of the projection matrix
[R,K]=qr(P(1:3,1:3));

//Extract the values from the calibration matrix
fx=K(1,1);
fy=K(2,2);
ox=K(1,3);
oy=K(2,3);

With a calibrated camera I could take two images of an object with the two camera positions side by side with distance called baseline b. Pixel coordinates for points on the object are determined on both images. Then the following equations can be used to solve the world coordinates of the points. (Video)

$$x=\frac{b(u_l-o_x)}{u_l-u_r}$$ $$y=\frac{bf_x(v_l-o_y)}{f_y(u_l-u_r)}$$ $$z=\frac{bf_x}{u_l-u_r}$$

But that doesn't give me correct coordinates because my calibration values are wrong. The X coordinates are in the correct order of magnitude but Y and Z are way off.

Can anyone spot what I've done wrong?

Answer 1

I found the answer. The videos I was looking at didn't mention a very important detail. The QR factorization needs to be applied to the inverse of the first three columns of P. The resulting K is also inversed.

Finally K needs to be normalized like this: $$K=\frac{1}{K_{3,3}}K$$

This video explained the calibration in greater detail.