function Fd = dual_fundmatrix(x1, x2)

%DUAL_FUNDMATRIX computes the dual fundamental matrix using 3 correspondences.
%
%   dualF = dual_fundmatrix(x1, x2)
%
%   The function computes the dual fundamental matrix, Fd, using 3 pairs
%   of matching points.
%
%   Input arguments:
%   - x1 and x2 must both be 3-by-3 matrix where the corresponding columns
%     contain a pair of matching points.
%
%   Output arguments:
%   - Fd would be a cell array of 3-by-3 matrix containing the estimtaed
%     dual fundamental matrix.  As there may be 1 or 3 possible solutions for
%     the dual fundamental matrix, each of them would be stored in the cell
%     array.
%
%   Reference: 
%   [1] R. Hartley and A. Zisserman, "Multiple View Geometry
%       in Computer Vision", Cambridge University Press 2000,
%       Chapter 19.
%
%   SEE ALSO trifocal_6pts, recons_3views
%
%Created in January 20 2004 by Du Huynh.
%Revised in March 2004 by Du Huynh
%
%Copyright Du Huynh
%The University of Western Australia
%School of Computer Science and Software Engineering

% check input arguments
if ~all(size(x1)==[3,3]) | ~all(size(x2) == [3,3])
   error('dual_fundmatrix: both x1 and x2 must be 3-by-3 matrices');
end

% This step below is similar to the code in fundmatrix_nonlin.
% The only difference is that we have a reduced fundamental matrix here,
% i.e. F has the form [0 p q; r 0 s; t u 0] where p+q+r+s+t+u=0
%
% See [1], page 491.
A = [];
for i=1:3
   A(i,:) = [x2(1,i)*x1(2,i) x2(1,i)*x1(3,i) ...
         x2(2,i)*x1(1,i) x2(2,i)*x1(3,i) ...
         x2(3,i)*x1(1,i)];
end
[v1,v2,s] = getnull(A);
v1(6) = -sum(v1(1:5));
v2(6) = -sum(v2(1:5));
FF = [0 v1(1) v1(2); v1(3) 0 v1(4); v1(5) v1(6) 0];
GG = [0 v2(1) v2(2); v2(3) 0 v2(4); v2(5) v2(6) 0];

% Step (v): We have Fhat = alpha*FF + (1-alpha)*GG
% The following coefficients are generated by Matlab's symbolic
% maths toolbox
coeff3 = ...
   FF(2,1)*FF(1,3)*FF(3,2)+GG(2,1)*GG(1,3)*FF(3,2)-GG(2,1)*GG(1,3)*GG(3,2)-...
   FF(2,1)*GG(1,3)*FF(3,2)-FF(2,1)*FF(1,3)*GG(3,2)-FF(3,1)*FF(1,2)*GG(2,3)+...
   FF(3,1)*FF(1,2)*FF(2,3)-GG(3,1)*FF(1,2)*FF(2,3)-GG(2,1)*FF(1,3)*FF(3,2)+...
   FF(2,1)*GG(1,3)*GG(3,2)+FF(3,1)*GG(1,2)*GG(2,3)-FF(3,1)*GG(1,2)*FF(2,3)+...
   GG(2,1)*FF(1,3)*GG(3,2)-GG(3,1)*GG(1,2)*GG(2,3)+GG(3,1)*GG(1,2)*FF(2,3)+...
   GG(3,1)*FF(1,2)*GG(2,3);
coeff2 = ...
   -2*GG(3,1)*GG(1,2)*FF(2,3)-2*GG(2,1)*GG(1,3)*FF(3,2)+...
   FF(2,1)*FF(1,3)*GG(3,2)+3*GG(2,1)*GG(1,3)*GG(3,2)+FF(3,1)*GG(1,2)*FF(2,3)+...
   FF(2,1)*GG(1,3)*FF(3,2)+GG(3,1)*FF(1,2)*FF(2,3)-2*FF(3,1)*GG(1,2)*GG(2,3)-...
   2*FF(2,1)*GG(1,3)*GG(3,2)+GG(2,1)*FF(1,3)*FF(3,2)+...
   3*GG(3,1)*GG(1,2)*GG(2,3)+FF(3,1)*FF(1,2)*GG(2,3)-...
   2*GG(3,1)*FF(1,2)*GG(2,3)-2*GG(2,1)*FF(1,3)*GG(3,2);
coeff1 = ...
   -3*GG(2,1)*GG(1,3)*GG(3,2)+GG(2,1)*GG(1,3)*FF(3,2)+GG(3,1)*FF(1,2)*GG(2,3)+...
   GG(2,1)*FF(1,3)*GG(3,2)+FF(2,1)*GG(1,3)*GG(3,2)-3*GG(3,1)*GG(1,2)*GG(2,3)+...
   FF(3,1)*GG(1,2)*GG(2,3)+GG(3,1)*GG(1,2)*FF(2,3);
coeff0 = ...
   GG(3,1)*GG(1,2)*GG(2,3)+GG(2,1)*GG(1,3)*GG(3,2);
alphas = roots([coeff3 coeff2 coeff1 coeff0]);
% There may be 1 or 3 real solutions.  We discard those complex solutions.
solns = alphas( find(imag(alphas) == 0) );
for i=1:length(solns)
   Fd{i} = solns(i)*FF + (1-solns(i))*GG;
end

% -------------------------------

% this function is designed for dual_fundmatrix.m.  Matrix A should be 3-by-5
function [v1,v2,s] = getnull(A)

% do column normalisation
for i=1:size(A,2)
   scale(i) = 1/norm(A(:,i));
end
Scale = diag(scale);
[U,S,V] = svd(A*Scale,0);
s = diag(S);
v1 = Scale*V(:,end-1);
v2 = Scale*V(:,end);
v1 = v1 / norm(v1,2);
v2 = v2 / norm(v2,2);
return