# Calculating camera calibration matrix with Scilab

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 xw, yw, zw 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 ATA. 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 fx and fy are the effective focal lengths in x and y directions and should be about the same if the pixels are square. The ox and oy 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?

• K should be normalized such that K(3,3) == 1. Also seems to me that K(1,2) is not close to zero (at least given these inputs), which seems to be referred to as $s == \alpha \cot \theta$ in other texts i find. Commented Jun 22, 2023 at 0:32
• I noticed that K(3,3) is usually shown as one but mine is not. However it isn't usually mentioned at all. I guess that normalization would be (1/K(3,3))K. But that doesn't give any good results either. Commented Jun 22, 2023 at 15:55
• I believe there is an arbitrary scaling here, thus the normalization. Everything else looks correct, assuming your camera holds for the pinhole model and whatever other assumptions are used for deriving this linear camera model (which is far outside my area of expertise). Even just using trying to use the projection directly P*[xw; yw; zw; 1] results didn't look great to me. Commented Jun 22, 2023 at 22:50

## 1 Answer

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.