Code covered by the BSD License  

Highlights from
Fundamental Matrix Computation

image thumbnail

Fundamental Matrix Computation

by

 

This package, implements the 3 standard algorithms for the computation of the fundamental matrix.

compare_results(X1,X2,try_noises,EXPERIMENT_NUM,F_names)
% compare_results(X1,X2,try_noises,EXPERIMENT_NUM,F_names)
% Compares the F s computed by different approaches in case of different additive noises
% and different number of points used in the computation of F. 
% The criterion for comparison is the Residual.
% input:
%               X1              3xN     (optional) homogeneous coordinates of matched points in view 1
%               X2              3xN     (optional) homogeneous coordinates of matched points in view 2
%               try_noises      1x1     (optional) if enabled, adds additive noise to the coordinates
%               EXPERIMENT_NUM  1x1     (optional) the number of experiments to perform and average the criterion on
%               F_names         {n}     (optional) cell array of names of the functions to try compute F
%                                       with. Functions will be called using [F] = <name>(x1,x2);
% Notes: 
%- if no data is provided, synthetic data will be used. synthetic data generation does NOT avoid degenerate configuration 
%   and does not guarantee well spread points(although it is tried here to avoid such situations). 
%- the results in the book were computed without any outliers(perhaps with quantization error imposed by "pixel" matching).
%- this code compares the results of algorithms in presence of noise and outliers. The corresponding coordinates are quantized to pixel entities,
%   however, in practice, using techniques like line fitting, it is possible to get sub-pixel accuracy. In that case, one can use the flag 
%   QUANTIZATION_ENABLED and set it to 0.
% 
% Author: Omid Aghazadeh, KTH(Royal Institute of Technology), 2010/05/09
function compare_results(X1,X2,try_noises,EXPERIMENT_NUM,F_names)
global TOL_X TOL_FUN MAX_FUN_EVAL MAX_ITER;
if nargin<5, 
%     F_names ={'normalized_8_point','algebraic_l1','algebraic_l2','geometric_gold_l1','geometric_gold_l2','geometric_sampson_l1'}; 
    F_names = {'normalized_8_point','algebraic_l1','algebraic_l1_not_nomralized','geometric_gold_l1','geometric_sampson_l1','geometric_gold_l1_not_normalized'}; 
end
if nargin<4, EXPERIMENT_NUM = 50; end
if nargin<3, try_noises = 1; end
if nargin<2
    if exist('synthetic_model.mat','file')
        load ('synthetic_model.mat');
    else
        range_pic = [640;480];
        N = 50;
        thresh_ang = 0.25;
        thresh_en = 0.33;
        while 1
            [X1,X2,P1,P2] = make_synthetic_data(N,[10;10;10],range_pic); 
            cam_ang = P1(:)'*P2(:) / norm(P1(:)) /norm(P2(:));
            [u1,s1,v1]  = svd(X1(1:2,:));
            [u2,s2,v2]  = svd(X2(1:2,:));
            s1e = diag(s1) / N^2;
            s2e = diag(s2) / N^2;
            if sum(s1e<thresh_en) || sum(s2e<thresh_en)|| cam_ang < thresh_ang, continue, end % trying to avoid degenerate configuration
            save('synthetic_model.mat','X1','X2');
            break;
        end
        
    end
end
QUANTIZATION_ENABLED = 1;
clc; close all;
N_alg = length(F_names);
if ~try_noises
    NOISE_ENERGY = 0;
    NOISE_TYPE = 0;
else
    MIN_ENERGY_NOISE = 1; % noise energy definition
    MAX_ENERGY_NOISE = 3;
    NOISE_ENERGY_COUNT = 2;
    NOISE_ENERGY = linspace(MIN_ENERGY_NOISE,MAX_ENERGY_NOISE,NOISE_ENERGY_COUNT);
    PROB_SPURIOUS = 1e-1; SPURIOUS_MULTIPLIER = 10; % spurious noise ratio corresponds to the ration of outliers in the estimates
    NOISE_TYPE = 1:3; % 0: no noise, 1: gaussian, 2: uniform, 3: spurious
end
N_noise_energy = length(NOISE_ENERGY);
N_noise_type= length(NOISE_TYPE);
TOL_X = 1e-4; TOL_FUN = 1e-4; MAX_FUN_EVAL = 1e4; MAX_ITER = 1e3; % parameters of the iterative minimization (LM) method
TRY_UB = 25; %size(X1,2)
% TRY_UB = round(2/3*size(X1,2));
RESIDUALS = zeros(N_noise_type,TRY_UB-7,N_alg,N_noise_energy,EXPERIMENT_NUM); RUN_TIMES = RESIDUALS;
Noise_Names = {'no noise','zm gaussian','uniform', 'spurious'}; 
Colors = {'r-','go-','g.-','bo-','b.-','bx-'};
% Colors = {'r-','g.-','b.-','b--'};
F_names_disp = F_names; 
for i =1:N_alg, F_names_disp{i}(F_names{i} == '_') = '-'; end;
for noise_type_i = 1:N_noise_type;
	noise_type = NOISE_TYPE(noise_type_i);
    for noise_energy_i = 1:N_noise_energy
        noise_energy = NOISE_ENERGY(noise_energy_i);
        for N = 8: TRY_UB
            for experiment = 1 : EXPERIMENT_NUM
                ixr = randperm(size(X1,2));
                x1 = X1(:,ixr(1:N));
                x2 = X2(:,ixr(1:N));

                x2n = x2; x1n=x1;
                switch noise_type
                    case 0 % no noise
                    case 1 % gaussian, zero mean!
                        x1n(1:2,:) = x1n(1:2,:) + randn(2,N)*noise_energy ;
                        x2n(1:2,:) = x2n(1:2,:) + randn(2,N)*noise_energy ;
                    case 2 % uniform, noise_param ~ energy
                        x1n(1:2,:) = x1n(1:2,:) + (rand(2,N)-0.5)*2*noise_energy;
                        x2n(1:2,:) = x2n(1:2,:) + (rand(2,N)-0.5)*2*noise_energy;
                    case 3 % spurious, noise_param ~ probability
                        r = rand(1,N);
                        ixnoise = find(r<PROB_SPURIOUS);
                        x1n(1:2,ixnoise) = x1n(1:2,ixnoise)+ (rand(2,length(ixnoise))-0.5)*2*noise_energy*SPURIOUS_MULTIPLIER;
                        x2n(1:2,ixnoise) = x2n(1:2,ixnoise)+ (rand(2,length(ixnoise))-0.5)*2*noise_energy*SPURIOUS_MULTIPLIER;
                end

                for algi = 1: N_alg
                    tic;
                    if QUANTIZATION_ENABLED
                        eval(sprintf('F=%s(round(x1n),round(x2n));',F_names{algi})); % compute the F using noisy points
                    else
                        eval(sprintf('F=%s((x1n),(x2n));',F_names{algi})); % compute the F using noisy points
                    end
                    
                    t = toc;

                    residual = get_residual(X1,X2,F); % testing F on the noiseless data, all points
                    RESIDUALS(noise_type_i,N-7,algi,noise_energy_i,experiment) = residual;
                    RUN_TIMES(noise_type_i,N-7,algi,noise_energy_i,experiment) = t;
                end
            end
        end
        
        figure; set(gcf,'Name',['Residuals/Noise Model ' Noise_Names{noise_type+1} ', noise energy=' num2str(noise_energy) ]); hold on;
        for i =1:N_alg
            plot(8:TRY_UB,log(mean(RESIDUALS(noise_type_i,:,i,noise_energy_i,:),5)),Colors{i});  % averaging over the experiments and plotting the data
        end
        title(sprintf('log Residual vs Number of Points, noise model: %s, noise energy: %f, Quantization enabled: %d',Noise_Names{noise_type+1},noise_energy,QUANTIZATION_ENABLED)); legend(F_names_disp); 

%         figure; set(gcf,'Name',['Run Time/Noise Model ' Noise_Names{noise_type+1} ', noise energy=' num2str(noise_energy) ]); hold on;
%         for i =1:N_alg
%             plot(8:TRY_UB,log(mean(RUN_TIMES(noise_type_i,:,i,noise_energy_i,:),5)),Colors{i});  % averaging over the experiments and plotting the data
%         end
%         title(sprintf('log runtime vs Number of Points, noise model: %s, noise energy: %f',Noise_Names{noise_type+1},noise_energy)); legend(F_names_disp); 
        drawnow;
    end
    
    
end
end

function res = get_residual(x1,x2,F)
d_x2_Fx1 = get_D_point_line(x2,F*x1); d_x1_Fx2 = get_D_point_line(x1,F'*x2);
res = mean(d_x2_Fx1.^2 + d_x1_Fx2.^2);
end

% function d = get_D_point_line(X1,Lp2)
% computes the distance of points X1 (homogeneous coordinates) to the lines
% Lp2. 
function d = get_D_point_line(X1,Lp2)
if sum(X1(3,:)~=1), error('the points are assumed to be normalized (last dimension = 1)'); end
d = abs(sum(X1 .* Lp2,1))./sqrt(sum(Lp2(1:2,:).^2));
end

function F= algebraic_l1(x1,x2), F = det_F_algebraic(x1,x2,1); end % sample functions to show the usage of algorithms
function F= algebraic_l2(x1,x2), F = det_F_algebraic(x1,x2,2); end
function F= algebraic_l1_not_nomralized(x1,x2), F = det_F_algebraic(x1,x2,1,0); end
function F= geometric_sampson_l1(x1,x2), F = det_F_gold(x1,x2,1,1); end
function F= geometric_gold_l1(x1,x2), F = det_F_gold(x1,x2,1); end
function F= geometric_gold_l1_not_normalized(x1,x2), F = det_F_gold(x1,x2,1,0,0); end
function F= geometric_gold_l2(x1,x2), F = det_F_gold(x1,x2,2); end
function F= normalized_8_point(x1,x2), F = det_F_normalized_8point(x1,x2); end

Contact us