I'm trying to implement Zhang's algorithm for camera calibration in MATlab. The steps to do that should be pretty straight forward and easy to implement, especially using MATLab. However, I've been stuck at a point where I'm getting complex numbers where I should be getting real numbers.
If any of you has done this before, you know exactly what has to be done.The first step would be to take images of a checkerboard and calculate the homography mapping between the model plane and the image plane using the below equation(I'm basically going to include most of my code because it's very small):
for i=1:n
p=[X(i,:) z -x(i,1).*X(i,1) -x(i,1).*X(i,2) -x(i,1)];
q=[z X(i,:) -x(i,2).*X(i,1) -x(i,2).*X(i,2) -x(i,2)];
E=[E;p;q];
end
where X holds the homogeneous coordinates of the checkerboard corners in the world plane and x holds the homogeneous coordinates in the image plane.
Next step would be solving the equation EH = 0 to get H:
[u d v]=svd(E);
H= v(:,end);
H= reshape(H,3,3)';
I'm not going to go through the details of what each variable is in the remaining code because I'm assuming that I'll be getting help from people who have already done this (it's a pretty popular algorithm), so I'll just put the code as is with some commentary:
vij = @(i,j,H) [ H(i,1)*H(j,1)
H(i,1)*H(j,2) + H(i,2)*H(j,1)
H(i,2)*H(j,2)
H(i,3)*H(j,1) + H(i,1)*H(j,3)
H(i,3)*H(j,2) + H(i,2)*H(j,3)
H(i,3)*H(j,3) ];
G = [ vij(1,2,H)'; (vij(1,1,H)-vij(2,2,H))' ];
V = [ V; G ];
of course in the above chunk of code we're inside a loop that runs for as many times as there are images
Up next is the most straight forward part of the whole algorithm where you just can't go wrong, where we calculate the intrinsic parameters:
[u1,d1,v1] = svd( V );
b = v1(:,end);
v0 = ( b(2)*b(4)-b(1)*b(5) ) / ( b(1)*b(3)-b(2)^2 );
lambda = b(6) - ( b(4)^2 + v0*(b(2)*b(4)-b(1)*b(5)) ) / b(1);
b(1)
alpha = sqrt( lambda / b(1) );
beta = sqrt( lambda*b(1) / (b(1)*b(3)-b(2)^2) );
gamma = -b(2)*alpha^2*beta / lambda;
u0 = gamma*v0 / beta - b(4)*alpha^2 / lambda;
A = [ alpha gamma u0;
0 beta v0;
0 0 1 ]
Then, we need to account for the difference in square sizes between the image plane and the world plane so we build a normalization matrix:
N = [ 2/width 0 -1
0 2/height -1
0 0 1 ];
Where height and width are those of the image in pixels.
Finally, the intrinsics matrix A can be obtained by the below line:
A = N\A;
Now, to make sure that I'm getting accurate results, I compared my results with the calibration tool from caltech ( Bouget ), and for some reason I'm getting u0 and v0 correctly but the rest ( alpha,beta,gamma ) are way off, in fact they're complex numbers.
And I realize that it's because lambda could be negative. Except that, it shouldn't be! and that's the issue!
I've looked at literally hundreds of implementations online and they are all very close to mine but all of them got ALL their intrinsics right. What's mindboggling to me is that my u0 and v0 are EXACTLY right, and the others are just complex jibberish.
Can you please help? I would really appreciate it!!
