function [H]=hfromnconics(Qs,Ss)
%
% hfromnconics solves a planar homography from N conic
% correspondences: Q_i <-> S_i
%
% INPUT
%  N real symmetric nonsingular 3x3 matrices:
%    Qs,Ss - 9xN matrices so that each column contains the conic
%    coefficients, Q_i=reshape(Q(:,i),3,3); S_i=reshape(S(:,i),3,3); 
%
% OUTPUT
%   H - the homography
%
% For details see the paper:
%
% Algorithms for computing a planar homography from conics in correspondence,
% J. Kannala, M. Salo and J. Heikkilä
% British Machine Vision Conference, 2006 (BMVC 2006)
%



N=size(Qs,2);

% compute all pairs of N correspondences
tmp=ones(N,N); tmp(1:N+1:N^2)=0;
[ytmp,xtmp]=find(tmp);
pairs=[xtmp ytmp];
%pairs=nchoosek(1:N,2);
Ncomb=size(pairs,1);

%normalize Frobenius norms to one
for k=1:N
  qnorms=sqrt(sum(Qs.^2,1)); 
  snorms=sqrt(sum(Ss.^2,1)); 
  Qs=Qs./repmat(qnorms,9,1);
  Ss=Ss./repmat(snorms,9,1);  
end

for k=1:N
  Qk=reshape(Qs(:,k),3,3); Sk=reshape(Ss(:,k),3,3);
  kk=det(Sk)/det(Qk);
  ds(k)=(abs(kk)^(1/3))*sign(kk);
end
Id=eye(3);
M=[];
for k=1:Ncomb
  l=pairs(k,1); m=pairs(k,2);
  Ql=reshape(Qs(:,l),3,3); Sl=reshape(Ss(:,l),3,3);
  Qm=reshape(Qs(:,m),3,3); Sm=reshape(Ss(:,m),3,3);
  Ak=ds(l)*inv(Sl)*Sm;
  Bk=ds(m)*inv(Ql)*Qm;
  bk=Bk(:);
  Mk=[Ak-bk(1)*Id -bk(2)*Id -bk(3)*Id;...
   -bk(4)*Id Ak-bk(5)*Id -bk(6)*Id;...
   -bk(7)*Id -bk(8)*Id Ak-bk(9)*Id];
  M=[M;Mk];
end
[U,S,V]=svd(M);
H=reshape(V(:,end),3,3);

%keyboard
