Project 3 / Camera Calibration and Fundamental Matrix Estimation with RANSAC

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:

  1. Part I: Camera Projection Matrix
  2. Part II: Fundamental Matrix Estimation
  3. Extra Credit: Normalization
  4. Part III: Fundamental Matrix with RANSAC

Part I: Camera Projection Matrix

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.

Part II: Fundamental Matrix Estimation

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.

Extra Credit: Normalization

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.

Part III: Fundamental Matrix with RANSAC

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: