Before i jump to the section on calibration of a single camera, it is worthwhile to present a refresher on basic concepts of orthogonality, properties of homogeneous linear equations, matrix decomposition and basics of Singular Value Decomposition.
What is an Orthogonal or Orthonormal Matrix?
If two vectors 
are perpendicular to each other, they are termed orthogonal, represented by 
. Their inner/dot product represented by 
is zero. For example vectors 
and 
are orthogonal to each other, since their dot product: 
is equal to zero.
Two vectors are orthonormal if they are orthogonal and are of unit length, i.e. norms of vectors are one. For example, vectors 
and 
are orthonormal because:
- Dot product:
,
- Norm of
, and
- Norm of
.
A similar concept is applied to matrices. A matrix is said to be orthogonal if:
- Matrix elements are real numbers (i.e. all floating point numbers that are not complex in nature),
- Matrix is a square matrix (i.e. same number of rows and columns), and
- Either column or row vectors are orthogonal.
In other words, if the product of a matrix 
and its transpose 
is equal to an identity matrix 
, i.e. 
. This is also true if the transpose of the matrix is the same as its inverse, i.e. 
.
If in such a matrix, column or row vectors are of unit length, then these matrices are called orthonormal matrices.
For example,
is an orthogonal and orthonormal matrix because,
- Dot product of column vectors:
,
- Norm of column vector
, and
- Norm of column vector
.
Also,
, and
Homogeneous Linear Equations
A set of equations where the constant terms are all zeros, are called homogeneous equations:
where, there are 
variables 
. This can be written in the matrix form:
or, 
Note, in the above system of equations, there always exits a solution: 
or, 
, called the “trivial solution”.
If the columns in the square matrix 
of size 
, are linearly independent, i.e. 
, then the trivial solution is the one and only one, unique solution. In other words, 
. This type of matrix is invertible and is referred to as “Regular Matrix”. On the other hand, if for this square matrix , 
, then the equations are not linearly independent, i.e. 
. Then, there exists infinite number of non-zero solutions, and is NOT invertible. Such a matrix is called the “Singular Matrix”.
What is Matrix Decomposition?
Matrix decomposition refers to reduction of the matrix to the simplest and most significant form-”basic building blocks”, without any loss of generality. For example, number 15 can be reduced or factorized into a multiple of 5*3, where both 5 and 3 are primes and can not be reduced further. A similar concept is applied to matrices, where a matrix can be represented by a multiple of basic matrices which can not be reduced further. One of the reasons why we decompose a matrix is to reduce a complex form into products of simpler forms, for example, a complex transformation applied to an image can be reduced to products of simpler transformations.
Singular Value Decomposition (SVD)
Singular Value Decomposition (SVD) is one of the many methods of factorization of a matrix with applications in least squares fitting of data, solving linear equations, computing the rank, range and null space of matrix, etc. A set of linear equations can be written in the form of 
, where is an 
matrix where rows 
are greater than the columns . Then we can factorize using SVD, where 
such that:
-
is a column orthogonal matrix of size ,
-
is a diagonal matrix with positive or zero elements of size , and
-
is an orthogonal matrix of size .
The above is termed as the Singular Value Decomposition (SVD) of . Let 
, where 
, then 
are called the singular values of . Columns of and are left and right singular vectors for the corresponding singular values respectively.
These singular values are not the same as the eigenvalues. For a matrix , matrix 
is normal with non-negative eigenvalues where singular values of are square roots of the eigenvalues of .
For a set of homogeneous linear equations 
, any vector 
which falls in the null space of is a solution. Any column where corresponding singular values approaches zero is a solution vector.
For a set of non-homogeneous linear equations 
, the goal is to find the solution with the smallest length 
, where ](http://s0.wp.com/latex.php?latex=x%3DV%5Bdiag%28+%5Csigma+_1+%5E+%7B-2%7D%2C+%5Csigma+_2+%5E+%7B-1%7D%2C+%5Ccdots+%2C+%5Csigma+_n+%5E%7B-1%7D+%29%5D%28U%5ETb%29&bg=ffffff&fg=000000&s=0 "x=V[diag( \sigma _1 ^ {-2}, \sigma _2 ^ {-1}, \cdots , \sigma _n ^{-1} )](U^Tb)")
. In this case, for every singular value 
where 
, 
is replaced by 0. If 
is not in the range of then there exists no vector to satisfy the equation. Thus it is not possible to obtain an exact solution.
If 
, i.e. is a square matrix, SVD is given by , where both and are orthogonal matrices. Thus, ](http://s0.wp.com/latex.php?latex=A%5E%7B-1%7D%3DV%5Bdiag%28+%5Csigma+_1+%5E+%7B-2%7D%2C+%5Csigma+_2+%5E+%7B-1%7D%2C+%5Ccdots+%2C+%5Csigma+_n+%5E%7B-1%7D+%29%5D%28U%5ET%29&bg=ffffff&fg=000000&s=0 "A^{-1}=V[diag( \sigma _1 ^ {-2}, \sigma _2 ^ {-1}, \cdots , \sigma _n ^{-1} )](U^T)")
.
The range of has the dimension of rank of . If is singular, then rank of will be less than , where 
. The columns of corresponding to singular values not equal to zero, are an orthonormal set for the range of , and columns of corresponding to singular values not equal to zero, are an orthonormal set for the null space of .
Calculating SVD
SVD can be calculated by first computing a bi-diagonal matrix of using Householder reflections, then computing SVD using an iterative method like QR algorithm for computation of eigenvalues [1]. In MATLAB there is a limit of 75 QR step iterations to get the singular values [2]. Sometimes, solution may not converge. LAPACK subroutines can be used for SVD : real double precision A-DGESVD, complex double precision A-ZGESVD, real single precision A-SGESVD, and complex single prevision A-CGESVD [2].
Camera Calibration
Here i am assuming you are familiar with concepts of camera calibration. For more details, please refer to Calibration slides by Bryan Morse at Brigham Young University. The work presented here is based on the Pinhole camera model.
If we know the 3D points in the world model, and the 2D projected pixel coordinates, then we can linearly solve (using SVD) for the calibration matrix, rotation matrix, translation vector, parameters for scaling in X and Y directions, and the optical center of the image (x,y). Using this we can re-project the 3D points and by comparing to the 2D points, figure out the average pixel error in the X and Y direction.
I have used the 3D points and 2D points provided by Dr. Boufama, at the University of Windsor. You can download these files from the links below:
- 3D Points
- 2D Points
Here is a visual representation of the points:
 3D Points |
 2D Points |
Here is the MATLAB Code:
% *************************************************************************
% Title: Function-Produce Calibration Matrix and Average Pixel Error using
% Singular value decomposition (SVD)
% Author: Siddhant Ahuja
% Created: February 2010
% Copyright Siddhant Ahuja, 2010
% ***Inputs***
% 1. File3D: 3-D File containing world coordinates
% 2. File2D: 2-D File containing pixel coordinates
% ***Outputs***
% calibMatrix: 3x4 Calibration Matrix
% rotationMatrix: 3x3 Rotation Matrix
% translationVector: 3x1 Translation vector
% alpha_u: Parameter for scaling in X-direction
% alpha_v: Parameter for scaling in Y-direction
% u_0: Optical Centre of the image (x value)
% v_0: Optical Centre of the image (y value)
% reprojMatrix: Reprojection of 3D points to 2D using Calibration Matrix
% avgError_u: Average Pixel Error in x direction
% avgError_v: Average Pixel Error in y direction
% timeTaken: Time taken by the code
% Example Usage of Function:
% [calibMatrix, rotationMatrix, translationVector, alpha_u, alpha_v, u_0,
% v_0, reprojMatrix, avgError_u, avgError_v, timeTaken]= funcCalibrate('3D.txt', '2D.txt');
% *************************************************************************
function [calibMatrix, rotationMatrix, translationVector, alpha_u, alpha_v, u_0, v_0, reprojMatrix, avgError_u, avgError_v, timeTaken] = funcCalibrate(File3D, File2D)
% Assuming no Radial Distortion is present.
% Initialize the timer to calculate the time consumed.
tic;
% Read 3-D file
[matrix3D]=textread(File3D);
% Delete First row of the 3D Matrix as it only contains number of data points
matrix3D(1,:)=[];
% Find the size of the 3D Matrix
[nR3D, nC3D]=size(matrix3D);
% Check to make sure matrix has 3 columns
if(nC3D~=3)
error('The matrix for 3D points does not have 3 columns.');
end
% Separate out the X values from the matrix
X = matrix3D(:,1);
% Separate out the X values from the matrix
Y = matrix3D(:,2);
% Separate out the X values from the matrix
Z = matrix3D(:,3);
% Plot the points in 3-D
figure;
plot3(matrix3D(:,1),matrix3D(:,2),matrix3D(:,3),'.');
% Read 2-D file
[matrix2D]=textread(File2D);
% Delete First row of the 2D Matrix as it only contains number of data points
matrix2D(1,:)=[];
% Find the size of the 2D Matrix
[nR2D, nC2D]=size(matrix2D);
% Check to make sure matrix has 2 columns
if(nC2D~=2)
error('The matrix for 2D points does not have 2 columns.');
end
% Separate out the u values from the matrix
u_values=matrix2D(:,1);
% Separate out the v values from the matrix
v_values=matrix2D(:,2);
% Plot the points in 2-D
figure;
plot(matrix2D(:,1),matrix2D(:,2),'.');
% Check to make sure number of rows of the 3D Matrix is the same as the
% number of rows of the 2D Matrix
if(nR3D~=nR2D)
error('Please make sure number of 3D and 2D points is the same.');
end
% ***Linear Solution of the Calibration Matrix***
% Let Calibtration Matrix be denoted by (M)
% Writing linear equations of the form AV=0, where A is a 2nx12 measurement
% matrix and V is a 12-element unknown vector.
% Create a Matrix of ones (o) with the same length as that of the u_values
% vector
o = ones(size(u_values));
% Create a Matrix of zeros (z) with the same length as that of the u_values
% vector
z = zeros(size(u_values));
% Populate the odd rows for A
AoddRows = [ X Y Z o z z z z -u_values.*X -u_values.*Y -u_values.*Z -u_values ];
% Populate the even rows for A
AevenRows = [ z z z z X Y Z o -v_values.*X -v_values.*Y -v_values.*Z -v_values ];
% Concatenate odd and even rows of A
A=[AoddRows; AevenRows];
% Compute SVD on the matrix
[U, S, V] = svd(A,0);
% Assuming no noise, since the elements of the diagonal matrix S are in descending order, to
% get the eigenvectors corresponding to the smallest eignevalue, we can
% just grab the last column of matrix
m = V(:,end);
% Construct the camera calibration matrix M
M = reshape(m,4,3)';
calibMatrix=M;
% Since the norm of Projection matrix M is equal to 1, we can calculate the
% absolute scale factor lambda
abs_lambda=sqrt(M(3,1)^2 + M(3,2)^2 + M(3,3)^2);
% Scale the Matrix with the scale factor
M = M / abs_lambda;
% In case the origin of the world frame is in front of the camera, we have
% s=lambda/abs_lambda, or From T_z=s*m_34, we can re-write s in the form of
% sign value of m_34. Here we assume it is in the front.
inFront=1;
if inFront
s = sign(M(3,4));
else
s = -sign(M(3,4));
end
% Thus, we can now calculate T_z or T(3)
T(3) = s*M(3,4);
% Create a 3x3 Rotation matrix and fill it with zeros
R = zeros(3,3);
% From equations, last row of the rotation matrix is the same as the first
% three elements of the last row of the calibration matrix
R(3,:)=s*M(3,1:3);
% Matrix M can be written as:
% M=( m1' )
% ( m2' m4)
% ( m3' )
% We can now calculate mi, where mi is a 3 element vector
m1 = M(1,1:3)';
m2 = M(2,1:3)';
m3 = M(3,1:3)';
m4 = M(1:3,4);
% Now, we can calculate the centres of projection u_0 and v_0
u_0 = m1'*m3;
v_0 = m2'*m3;
% Calculating the alpha values in u and v directions,
alpha_u=sqrt( m1'*m1 - u_0^2 );
alpha_v=sqrt( m2'*m2 - v_0^2 );
% We can now calculate the first and second rows of the rotation matrix
R(1,:) = s*(u_0*M(3,1:3) - M(1,1:3) ) / alpha_u;
R(2,:) = s*(v_0*M(3,1:3) - M(2,1:3) ) / alpha_v;
% We can also calculate the first and second elements of the Translation
% vector
T(1) = s*(u_0*M(3,4) - M(1,4) ) / alpha_u;
T(2) = s*(v_0*M(3,4) - M(2,4) ) / alpha_v;
T = T';
translationVector=T;
% The rotation matrix R obtained with this estimation procedure is not guaranteed to be orthogonal.
% Therefore we calculate the rotation matrix that is closest to the estimated matrix (in the
% Frobenius norm sense). Let R = UDV'T then R = UV'T.
[U,D,V] = svd(R);
R = U*V';
rotationMatrix=R;
% Reproject 3Dpoints to 2D points using the calibration matrix to calculate average pixel errors
newMatrix2D=zeros(nR3D,2);
for i=1:1:nR3D
num_u=M(1,1)*matrix3D(i,1) + M(1,2)*matrix3D(i,2) + M(1,3)*matrix3D(i,3) + M(1,4);
num_v=M(2,1)*matrix3D(i,1) + M(2,2)*matrix3D(i,2) + M(2,3)*matrix3D(i,3) + M(2,4);
den=M(3,1)*matrix3D(i,1) + M(3,2)*matrix3D(i,2) + M(3,3)*matrix3D(i,3) + M(3,4);
newMatrix2D(i,1)=num_u/den;
newMatrix2D(i,2)=num_v/den;
end
% Reprojected matrix
reprojMatrix=newMatrix2D;
% Calculate difference between the reprojectedMatrix and the original 2D
% Matrix
errorDiff=reprojMatrix-matrix2D;
% Calculate average pixel error in x direction
avgError_u=mean(errorDiff(:,1));
% Calculate average pixel error in y direction
avgError_v=mean(errorDiff(:,2));
% Stop the timer to calculate the time consumed.
timeTaken=toc;
To run the above code, please type in the following commands in MATLAB:
% Run this first
clc
clear all
[calibMatrix, rotationMatrix, translationVector, alpha_u, alpha_v, u_0, v_0, reprojMatrix, avgError_u, avgError_v, timeTaken]= funcCalibrate('3D.txt', '2D.txt');
Here are some of the results obtained by running the code:
- Calibration Matrix:
- Rotation Matrix:
- Translation Vector:
- Parameter for scaling in X-direction: 999.9978
- Parameter for scaling in Y-direction: 999.9978
- Optical center of the image: (256.0000, 256.0000)
- Average Pixel Error in X direction: -3.7259e-11
- Average Pixel Error in Y direction: -1.9027e-11
References
[1] Trefethen, Lloyd N.; Bau III, David (1997), Numerical linear algebra, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-361-9.
[2] http://www.mathworks.com/access/helpdesk/help/techdoc/ref/svd.html
So what is a filter?
. This property if called additive.
, where
is just a constant. This property is called Homogeneity, which implies if we double the strength of the input, the output will double up too. If we triple the input, the output’s strength gets tripled too.
,
