function [dataout]  = audio_abs_flat( S ) 
%function [dataout]  = audio_abs_flat( S ) 
%
% Given 6 points in the plane for which the distance from any of the 
% first 3 to any of the last 3 is known the possible solutions for the 
% structure are computed. 
%
% For example the first 3 can be microphones and the last 3 speakers. 
%
% Input matrix S is a 3x3 measurement matrix containing the squared distances. 
%
% Output matrix dataout contains the solutions as columns. There are up to 8 solutions which may be complex. 
%
% The code relies on buildM.c being compiled as mex. 
%
%
%
% %Example Code: 
%  X = [0 abs(rand(1)) 4+rand(1); 0 0 abs(rand(1)) ];
%  Y = randn( 2,3);
%  audio_observe33=inline('reshape(sum((X(:,[1:3 1:3 1:3])-Y(:,[1 1 1 2 2 2 3 3 3])).^2), 3, 3)','X','Y');
%  S = audio_observe33( X, Y);
%  Xest = audio_abs_flat( S ) ;
% %The correct answer is [X(:);Y(:)]
% 
%
% 
% A full description of the solver is given in my thesis. If you use the 
% solver please refer to: 
%
% @PHDTHESIS{stewenius-pdh-2005,
%    AUTHOR = {H. Stew\'enius},
%    TITLE = {Gr{\"o}bner Basis Methods for Minimal Problems in Computer Vision},
%    SCHOOL = {Lund University},
%    YEAR = 2005,
%    MONTH = APR,
%    PDF = {http://www.maths.lth.se/matematiklth/personal/stewe/THESIS/stewenius_phd_2005.pdf}
% }

%Copyright Henrik Stewenius 2006. 
%
% You are free to use the code for research. For commercial applications
% different licensing is possible. Use the code at own risk. 
%

if( not(exist('audio_abs_flat_helper')==3))
  if(not(exist('audio_abs_flat_helper.c')))
    error('You must download and mex audio_abs_flat_helper.c'); 
  else
    error(['You have not mexed (compiled) audio_abs_flat_helper.c'...
           'Mex by typing "mex audio_abs_flat_helper.c"'])
  end 
end 


%Renormalize to improve numerics. 
NORMS = norm(S, 'fro');
S = S/NORMS;

dataout =[];
S2 = [];
for i=1:3
   s0 = S(1,i);
   s1 = S(2,i);
   s2 = S(3,i);
   S2t =[ 1 s0 s1 s2 s1^2+s0^2-2*s0*s1 s1*s2-s0*s2-s0*s1+s0^2 s2^2-2*s0*s2+s0^2];
   S2 = [S2; S2t ]; 
end



%Computing the nullspace by using SVD. 
[U,D,V] = svd( S2);
X = V(:,end-3:end);

%Build the system matrix
M3 = audio_abs_flat_helper(X);
M3 = M3([1:30 end],:); 
%M3 = M3([4 5 7:25 26:29 30 31],:);
% keyboard

%Gauss-Jordan elimination- 
M4 = M3(:,1:26)\M3;


%Build the action matrix. 
M = [  -M4([ 11 12 14 17 21 22 24],end-8:end);
       1  0  0  0  0  0  0  0  0
       0 0  0  0  1  0  0  0  0];

%Solve the eigenproblem. 
[V,D] = eig( M);
%Backsubstitute to get the values for X and Y
for cnt=2:9
   if ( imag( D(cnt,cnt))==0 ) 
      xx = V(5:8,cnt)/V(end,cnt);
      Xest = X*xx/8; 
      
      xx0 = sqrt( -Xest(7) );
      xx1 = 0.5*Xest(6)/xx0;
      xx2 = sqrt( -Xest(5)-xx1^2);
%      disp( [xx0 xx1 xx2])
      if 1% (imag(xx0) == 0  ) & (imag(xx1) == 0  ) & (imag(xx2) == 0  )	 
	 Xvec =[ 0  xx0 xx1;0 0 xx2];
	 %%Generating the observation vectors
	 S2 = [];
	 
	 Yvec =[];
	 for i=1:3
	    Xtemp = [0 xx0 xx1; 0 0 xx2 ];
	    s0 = S(1,i);
	    s1 = S(2,i);
	    s2 = S(3,i);
	    A = [ 0 , 0 , 0  , 1  , -s0 , 1 ;
		  -2*xx0  , 0 , 0  , 0  , xx0^2+s0-s1 , 0 ;
		  -2*xx1  , -2*xx2 , 0  , 0  , xx1^2+xx2^2+s0-s2 , 0 ;
		  xx0^2+s0-s1 , 0 , 0  , 2*xx0  , -2*xx0*s0  , 0 ;
		  xx1^2+xx2^2+s0-s2 , 0 , -2*xx2 , 2*xx1  , -2*xx1*s0  , 0 ;
		  0 , xx0^2+s0-s1 , -2*xx0 , 0  , 0 , 0 ;
		  0 , xx1^2+xx2^2+s0-s2 , -2*xx1 , -2*xx2 , 0 , 0 ] ;
	    [UU,SS,VV]= svd( A ); 
	    temp= VV(:,end); 
	    %temp  = null( A ) ; 
	    Yvec = [Yvec temp(1:2)/temp(5) ];
	 end
	 dataout = [dataout  [Xvec(:) ;Yvec(:)]];
      end
   end
end

%Go back to the original problem. 
dataout = dataout*sqrt( NORMS);