function [F_ans, inliers] = computeF_RANSAC(corre, how_close)
%Joshua Stough
%COMP 290-89, Pollefeys
%Assignment 2, Computing F
%Returns the fundamental matrix F, computed using the 
%RANSAC robust estimator and the 7-point computation
%of F.  Also returns, for the best case F, the 4xn 
%matrix of inliers (a column subset of the
%correspondences argument corr, also 4xm).
%I compute inliers using the symmetric epipolar
%distance, 10.10, p. 271.

global corr

corr = corre;

if (nargin < 2) how_close = .5; end

%Algorithm 3.5, p. 105, and Alg. 10.4, p. 275.
N = 1000;
sc = 0;
lgp = log(.01);

x = corr(1:2,:);
xp = corr(3:4,:);
s = size(corr);
x(3,:) = ones(1,s(2));
xp(3,:) = ones(1,s(2));

epsilon = 1;

while N > sc
    %Pick random sample of seven point correspondences.
    which_seven = round(1 + (s(2)-1)*rand(1,7));
    
    %compute the fundamental matrix or matrices.
    %This function returns a 3x3x3, with zeros in
    %the last two slices if there is only one solution.
    Fs = computeF_7point_ransac(which_seven);
    
    for i = 1:3
        F = Fs(:,:,i);
        %Now compute inliers and outliers.
        %the lines corresponding to the points.
        if sum(sum(F)) == 0
            disp('ERROR: F is zero matrix.');
        end
        
        lx = F*x;
        lxp = F'*xp;
        %normalize the lines.
        lxnorms = sqrt(sum(lx(1:2,:).^2));
        lxpnorms = sqrt(sum(lxp(1:2,:).^2));
        lxnorms(lxnorms == 0) = 1e-15;
        lxpnorms(lxpnorms == 0) = 1e-15;
        
        lx = lx./(ones(3,1)*lxnorms);
        lxp = lxp./(ones(3,1)*lxpnorms);
        
        %perpedicular distance from a point to the line it
        %is supposed to be on is the dot product.
        %square distances in image1 and two.
        dist_1 = abs(sum(lxp.*x));
        dist_2 = abs(sum(lx.*xp));
        
        %less than 1/2 pixel is good.
        ins = (dist_1 < how_close) & (dist_2 < how_close);
        %this means the correspondence errors on both sides
        %are less than 1/2 pixel.
        k = find(ins);
        
        temp = 1 - nnz(ins)/s(2);
        if temp < epsilon 
            inliers = corr(:, k');
            epsilon = temp;
            F_ans = F;
            disp(sprintf('%f %f %f',[epsilon sc N]));
        end
    end
    
    %Recompute N and inrement the sample count.
    N = lgp/log(1 - (1-epsilon)^7);
    sc = sc + 1;
end
