This project is about camera calibration and fundamental matrix estimation. I successfully finished all three parts of the project and also the extra credit (Normalization). The outline of this report is as follows:
I established a homogenous set of equations using the corresponding 2D and 3D points and then used the least square with SVD to solve the projection matrix. The calculate_projection_matrix
function is written as follows:
% solve the projection matrix
function M = calculate_projection_matrix( Points_2D, Points_3D )
for i = 1:20
A(2*i-1,:)=[Points_3D(i,:),1,zeros(1,4),-Points_2D(i,1)*Points_3D(i,:),-Points_2D(i,1)];
A(2*i,:)=[zeros(1,4),Points_3D(i,:),1,-Points_2D(i,2)*Points_3D(i,:),-Points_2D(i,2)];
end
[U,S,V] = svd(A);
M = V(:,end);
M = reshape(M,[],3)';
end
According to the property of the projection matrix, I can compute the camera center easily. The compute_camera_center
function is written as follows:
% compute the camera center
function [ Center ] = compute_camera_center( M )
Center = -inv(M(:,1:3))*M(:,4);
end
The running results in Matlab is shown as follows:
I can see the results match the correct answer (the residual is very small and the projection matrix is scaled verision of the given M matrix.), which verifies our implementation. I show the following figures for 2D projected points estimation and the camera center. From the following figures, I can see that our implementation works well.
I first implemented the fundamental matrix estimation without normalization of the input points. I used the 8 points algorithm to perform the estimation. The estimate_fundamental_matrix
function is written as follows:
% estimate fundamental matrix without normalization
function [ F_matrix ] = estimate_fundamental_matrix(Points_a,Points_b)
num = size(Points_a,1);
for i = 1:num
A(i,:) = [Points_a(i,1)*Points_b(i,1),Points_b(i,1)*Points_a(i,2),...
Points_b(i,1),Points_a(i,1)*Points_b(i,2),Points_b(i,2)*...
Points_a(i,2),Points_b(i,2),Points_a(i,1),Points_a(i,2),1];
end
[~, ~, V] = svd(A);
f = V(:, end);
F_matrix = reshape(f, [3 3])';
[U, S, V] =svd(F_matrix);
S(3,3) = 0;
F_matrix = U*S*V';
end
The results of epipolar lines are shown as follows:
From the results, I can see the epipolar lines are basically accurate, but there is still a liitle offset from the detected points. I improve this estimation using normalization discussed in the next section.
The normalization is mainly used to improve the poor numerical conditioning of the eight points algorithm. To implement this part, I add a new function (normalize.m
) to perform the normalization. The normalize
function is shown as follows:
% normalize the input 2D points
function [newpts, T] = normalize( pts )
centroid = sum(pts(:,1:2)) / size(pts,1);
diff = pts(:,1:2) - repmat(centroid,size(pts,1),1);
avgdist = sum(sqrt(diff(:,1).^2 + diff( :,2).^2)) / size(pts,1);
scale = sqrt(2) / avgdist;
T = diag([scale scale 1]) * [eye(3,2) [-centroid 1]'];
newpts = pts;
newpts(:,3) = 1;
newpts = newpts * T';
newpts = newpts ./ repmat( newpts(:,3), 1, 3);
newpts(:,3) = [];
end
Using the normalize
function, I reimplement the eight points algorithm, which is basically normalized eight points algorithm. The new estimate_fundamental_matrix
function is written as follows:
% estimate fundamental matrix with normalization
function [ F_matrix ] = estimate_fundamental_matrix(Points_a,Points_b)
[Points_a,T_a] = normalize( Points_a );
[Points_b,T_b] = normalize( Points_b );
num = size(Points_a,1);
for i = 1:num
A(i,:) = [Points_a(i,1)*Points_b(i,1),Points_b(i,1)*Points_a(i,2),...
Points_b(i,1),Points_a(i,1)*Points_b(i,2),Points_b(i,2)*...
Points_a(i,2),Points_b(i,2),Points_a(i,1),Points_a(i,2),1];
end
[~, ~, V] = svd(A);
f = V(:, end);
F_matrix = reshape(f, [3 3])';
[U, S, V] =svd(F_matrix);
S(3,3) = 0;
% S = S / S(1,1);
F_matrix = U*S*V';
F_matrix = T_b' * F_matrix * T_a;
end
The results of new epipolar lines in part II are shown as follows:
From the above results, I can see the estimation of epipolar lines are more accurate than the original ones. (You will be able to notice the difference with larger resolution.)
In the following implementation, I will use the normalized eight points instead of the original eight points algorithm.
I used the normalized eight points algorithm along with RANSAC to estimate the "best" fundamental matrix. I considered to use the property of x^T*F*x'=0 to perform the RANSAC algorithm. The ransac_fundamental_matrix
function is shown as follows:
% compute the fundamental matrix via RANSAC (used in the final code)
function [ Best_Fmatrix, inliers_a, inliers_b] = ransac_fundamental_matrix(matches_a, matches_b)
matches_a_h = [matches_a,ones(size(matches_a,1),1)];
matches_b_h = [matches_b,ones(size(matches_b,1),1)];
% dThresh = 46;
num = size(matches_a,1);
indi_f_old = inf;
% RANSAC
for i = 1:30000
randnum = randperm(num);
Best_Fmatrix_m = estimate_fundamental_matrix(matches_a(randnum(1:8),:), matches_b(randnum(1:8),:));
% use the property of x^T*F*x'=0 to do RANSAC
va = abs(diag(matches_b_h*Best_Fmatrix_m*matches_a_h'));
inliers = find(va<0.05);
inlierCount = size(inliers,1);
%%
if inlierCount > 0
indi_f = 1-inlierCount/num;
if indi_f < indi_f_old
Best_Fmatrix = Best_Fmatrix_m;
inlierValue_best = va(inliers);
inliers_best = inliers;
indi_f_old = indi_f;
end
end
end
[~,ind] = sort(inlierValue_best,'ascend');
inliers_a = matches_a(inliers_best(ind(1:30)),:);
inliers_b = matches_b(inliers_best(ind(1:30)),:);
end
Apart from the above implementation, I also tried to use the Euclidean distance between epipolar lines and detected points to do RANSAC. The implementation is shown as follows:
% compute the fundamental matrix via RANSAC (Euclidean distance version, not used in final code)
% it is just for report.
function [ Best_Fmatrix, inliers_a, inliers_b] = ransac_fundamental_matrix(matches_a, matches_b)
matches_a_h = [matches_a,ones(size(matches_a,1),1)];
matches_b_h = [matches_b,ones(size(matches_b,1),1)];
dThresh = 46;
num = size(matches_a,1);
indi_f_old = inf;
% RANSAC
for i = 1:30000
randnum = randperm(num);
Best_Fmatrix_m = estimate_fundamental_matrix(matches_a(randnum(1:8),:), matches_b(randnum(1:8),:));
% use L2 distance between epipolar line and the points to do RANSAC
L1_m = Best_Fmatrix_m*matches_a_h';
L1 = L1_m ./ repmat( (sqrt( L1_m(1,:).^2 + L1_m(2,:).^2 )), 3, 1 );
dist1 = abs(dot( matches_b_h', L1 ));
L2_m = Best_Fmatrix_m'*matches_b_h';
L2 = L2_m ./ repmat( (sqrt( L2_m(1,:).^2 + L2_m(2,:).^2 )), 3, 1 );
dist2 = abs(dot( L2, matches_a_h' ));
inliers = find( dist1 < dThresh & dist2 < dThresh );
inlierCount = size(inliers,2);
if inlierCount > 0
indi_f = sum( dist1(inliers).^2 + dist2(inliers).^2 ) / inlierCount;
if indi_f < indi_f_old
Best_Fmatrix = Best_Fmatrix_m;
dist1_best = dist1;
dist2_best = dist2;
inliers_best = inliers;
indi_f_old = indi_f;
end
end
end
% L2 distance
confidence = dist1_best(inliers_best) + dist2_best(inliers_best);
[~,ind] = sort(confidence,'ascend');
inliers_a = matches_a(ind(1:min(20,size(ind,2))),:);
inliers_b = matches_b(ind(1:min(20,size(ind,2))),:);
end
The results for the "Mount Rushmore" pairs are shown as follows:
The results for the "Notre Dame" pairs are shown as follows:
The results for the "Gaudi" pairs are shown as follows:
The results for the "Woodruff Dorm" pairs are shown as follows: