%  General testing program
%  Calls various projection methods, forms histograms of bad cases

clear all;
warning off all;

Window = 1;      % size of rectangle containing all data points

n = 2;  %  number of data points

histCols = 20;   %  number of columns in histograms (maximum recordable precision)
histRows = 6;    %  number of rows in histograms (number of projection methods)

Emin = 10.^(-histCols);  %  min allowed error 

hist1=zeros(histRows,histCols); hist2=zeros(histRows,histCols);   %  files for histograms

N=0; N0=0; N1=0; N2=0; N3=0; Nmiss=0;    %  various counters of randomly generated conics

Times = zeros(histRows,1);      %  for timing

E1 = zeros(histRows,n);  E1 = zeros(histRows,n);
XYbest = zeros(n,2);  XYprojAll = zeros(n,2,histRows);

while(1)     %  main loop over random conics and points
    
    A = randn(6,1);   %  random conic parameters
    
    % the following line makes the conic nearly-parabolic (ill-conditioned)
    % this line can be commented out when generating typical conics
    
    %e=1e-13; A(3) = A(2)^2/A(1) + e*randn;   % make it a parabola (if desired)
    
    A = A/norm(A);    %  normalization (unnecessary)
 
	%  discard the conic if it is of the wrong type:
	   
    if (IsCrossingWindow(A,Window) == 0)  % discard the conic if it doed not
        Nmiss = Nmiss + 1;                % cross the data window
        continue;
    end;

    [ParG,code] = AtoG(A);   %  ParG is vector of geometric parameters of the conic
    if (code > 3)            %  discard degenerate conics
        N0 = N0 + 1;
        continue;
    end;
    N = N + 1;
    if (code==1), N1 = N1 + 1; end;   %  count ellipses
    if (code==2), N2 = N2 + 1; end;   %  count hyperbolas
    if (code==3), N3 = N3 + 1; end;   %  count parabolas

    XY = 2*Window*(rand(n,2)-0.5*ones(n,2));  %  random points

    Dmin = realmax*ones(n,1);
    
    tic
    XYproj = ProjectPointsOntoConicByWEPqz(XY,A);    %  1st projection method
    Times(1) = Times(1) + toc;
    E2(1,:) = ProjectionQuality(XY,XYproj,A);
    D = DistanceToConicApprx(XY,XYproj,A);
    for i=1:n
        if Dmin(i) > D(i)
            Dmin(i) = D(i);
            XYbest(i,:) = XYproj(i,:);
        end
    end
    XYprojAll(:,:,1) = XYproj;
    
    tic
    XYproj = ProjectPointsOntoConicByWEPinv(XY,A);    %  2nd projection method
    Times(2) = Times(2) + toc;
    E2(2,:) = ProjectionQuality(XY,XYproj,A);
    D = DistanceToConicApprx(XY,XYproj,A);
    for i=1:n
        if Dmin(i) > D(i)
            Dmin(i) = D(i);
            XYbest(i,:) = XYproj(i,:);
        end
    end
    XYprojAll(:,:,2) = XYproj;

    tic
    XYproj = ProjectPointsOntoConicByWEPeq3Schur(XY,A);    %  3rd projection method
    Times(3) = Times(3) + toc;
    E2(3,:) = ProjectionQuality(XY,XYproj,A);
    D = DistanceToConicApprx(XY,XYproj,A);
    for i=1:n
        if Dmin(i) > D(i)
            Dmin(i) = D(i);
            XYbest(i,:) = XYproj(i,:);
        end
    end
    XYprojAll(:,:,3) = XYproj;

    tic
    XYproj = ProjectPointsOntoConicByEberlyModified(XY,A);    %  4th projection method
    Times(4) = Times(4) + toc;
    E2(4,:) = ProjectionQuality(XY,XYproj,A);
    D = DistanceToConicApprx(XY,XYproj,A);
    for i=1:n
        if Dmin(i) > D(i)
            Dmin(i) = D(i);
            XYbest(i,:) = XYproj(i,:);
        end
    end
    XYprojAll(:,:,4) = XYproj;

    tic
    XYproj = ProjectPointsOntoConicByEberlyOrig(XY,A);    %  5th projection method
    Times(5) = Times(5) + toc;
    E2(5,:) = ProjectionQuality(XY,XYproj,A);
    D = DistanceToConicApprx(XY,XYproj,A);
    for i=1:n
        if Dmin(i) > D(i)
            Dmin(i) = D(i);
            XYbest(i,:) = XYproj(i,:);
        end
    end
    XYprojAll(:,:,5) = XYproj;

    tic
    XYproj = ProjectPointsOntoConicByQuarticEquation(XY,A);    %  6th projection method
    Times(6) = Times(6) + toc;
    E2(6,:) = ProjectionQuality(XY,XYproj,A);
    D = DistanceToConicApprx(XY,XYproj,A);
    for i=1:n
        if Dmin(i) > D(i)
            Dmin(i) = D(i);
            XYbest(i,:) = XYproj(i,:);
        end
    end
    XYprojAll(:,:,6) = XYproj;
	
	%   Record the results of all projection methods
    
    %   first, for each projected point compute the distance to the best projection
    
    for i=1:n
        for j=1:histRows
            E1(j,i) = norm(XYprojAll(i,:,j)-XYbest(i,:));
            E1(j,i) = max(E1(j,i),Emin);
            E2(j,i) = max(E2(j,i),Emin);
        end
    end
    
    %  next for each projected point record the errors computed in two ways: 
    %  E1 stores the distance to the best projected point
    %  E2 stores the quality of the projection (E in the article)
    %   update the histograms
    
    for i=1:n
        for j=1:histRows
            k = floor(-log10(E1(j,i)));
            if (k<1), k=1; end
            if (k<=histCols), hist1(j,k) = hist1(j,k) + 1;; end
            k = floor(-log10(E2(j,i)));
            if (k<1), k=1; end
            if (k<=histCols), hist2(j,k) = hist2(j,k) + 1;; end
        end
    end
    
    %  print out the results
    if mod(N*n,10000)==0
        fprintf(1,'\n  Histogram of distances to the best projection point:\n\n');
        for j=1:histRows
            for k=1:histCols
                fprintf(1,'  %5d',hist1(j,k));
            end
            fprintf(1,'\n');
        end
        fprintf(1,'\n  Histogram of qualities of projections:\n\n');
        for j=1:histRows
            for k=1:histCols
                fprintf(1,'  %5d',hist2(j,k));
            end
            fprintf(1,'\n');
        end
        fprintf(1,'\n  Timing:')
        for j=1:histRows
            fprintf(1,'  %8.5f',Times(j)/N);
        end
        fprintf(1,'\n\n conics used: %d  (Ell: %d  Hyp: %d  Par: %d)\n',N,N1,N2,N3)
        fprintf(1,'\n conics removed: %d degenerate and %d missing the window\n',N0,Nmiss)
    end
    if (N*n>=100000000)
        return;
    end
end