function [F] = computeF_8point(corr, which_eight)
%Joshua Stough
%COMP 290-89, Pollefeys
%Assignment 2, Computing F
%Returns the fundamental matrix F, computed using the 
%normalized eight point algorithm.  If which_eight is
%not specified, all the correspondences are used.

s = size(corr);
if (nargin < 2) 
    which_eight = 1:s(2);
    n = s(2);
else
    n = 8;
end

%First, normalize the input coordinates, zero-mean and
%root 2 distance to the origin.
x = corr(1:2,which_eight);
xp = corr(3:4,which_eight);

mx = mean(x');
x(1,:) = x(1,:)-mx(1);
x(2,:) = x(2,:)-mx(2);
%centroid is origin. Now sum of squares should
%average 2.
s = mean(sum(x.^2));
s = sqrt(2/s);
x = x*s;
T = [[s 0 0]' [0 s 0]' [-s*mx(1) -s*mx(2) 1]']; 
%Now the average magnitude is sqrt(2);
%same for xp.
mx = mean(xp');
xp(1,:) = xp(1,:)-mx(1);
xp(2,:) = xp(2,:)-mx(2);
%centroid is origin. Now sum of squares should
%average 2.
s = mean(sum(xp.^2));
s = sqrt(2/s);
xp = xp*s;
TP = [[s 0 0]' [0 s 0]' [-s*mx(1) -s*mx(2) 1]']; 

%So, x and xp are corresponding points and they're
%normalized and zero-mean.  
%Now, build A.
A = zeros(n,9);
%build each column.
A(:,1) = (xp(1,:).*x(1,:))';
A(:,2) = (xp(1,:).*x(2,:))';
A(:,3) = xp(1,:)';
A(:,4) = (xp(2,:).*x(1,:))';
A(:,5) = (xp(2,:).*x(2,:))';
A(:,6) = xp(2,:)';
A(:,7) = x(1,:)';
A(:,8) = x(2,:)';
A(:,9) = ones(n,1);

%A is now built.  Find the constrained best 
%solution for F.
[U,S,V] = svd(A);

f = [V(1:3,9)'; V(4:6,9)'; V(7:9,9)'];
%We actually want the F that minimizes the Frobenius
%norm |f - F| subject to det F = 0.
[U,S,V] = svd(f);
S(3,3) = 0;
F = U*S*V';

%Denormalize.
F = TP'*F*T;

%That's it.