function XYproj = ProjectPointsOntoConicByWEPinv(XY,A)
%
%  The WEP method
%  version: solving an ordinary eigenvalue problem 
%           after inverting one of the two matrices
%
%   Projecting a given set of points onto a quadratic curve (conic)
%   and computing the distances from the points to the conic

%   The conic is defined by quadratic equation 
%      Ax^2 + 2Bxy + Cy^2 + 2Dx + 2Ey + F = 0

%   Input:  A = (A,B,C,D,E,F)' is the column vector of the parameters of the conic
%           XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)

%   Output:  RSS is the Residual Sum of Squares (the sum of squares of the distances)
%            XYproj is the array of coordinates of projections
%
%
%        Nikolai Chernov,  February 2012

n = size(XY,1);       %  n is the number of points to be projected
XYproj = zeros(n,2);  %  prepare an array of projections
xy = zeros(4,2);

C = [A(1) A(2) A(4);
     A(2) A(3) A(5);
     A(4) A(5) A(6)];   %    3x3 matrix of conic's parameters
 
 for i=1:n
     u = XY(i,1);  v = XY(i,2);   %  coordinates of the point to be projected
     
                  %  Next compute D, the auxiliary matrix 
                  %  it comes from the orthogonality condition
                  
     D = [2*A(2) A(3)-A(1) A(1)*v-A(2)*u+A(5);
         A(3)-A(1) -2*A(2) A(2)*v-A(3)*u-A(4);
         A(1)*v-A(2)*u+A(5) A(2)*v-A(3)*u-A(4) 2*A(4)*v-2*A(5)*u];

     if abs(det(C)) < abs(det(D))
         E = eig(inv(D)*C);

         %   next find a real generalized eigenvalue (just one is enough)

         jreal=0;
         for j=1:3
             if isreal(E(j))
                 jreal=j;
                 break;
             end
         end
         if jreal==0
             disp(' No real e-values...');
             return;
         end

         G = C-E(jreal)*D;    %  parameters of the degenerate conic
         %  in the family of conics C-t*D
     else
         E = eig(inv(C)*D);

         %   next find a real generalized eigenvalue (just one is enough)

         jreal=0;
         for j=1:3
             if isreal(E(j))
                 jreal=j;
                 break;
             end
         end
         if jreal==0
             disp(' No real e-values...');
             return;
         end

         G = D-E(jreal)*C;    %  parameters of the degenerate conic
         %  in the family of conics D-t*C
     end

     if det(G(1:2,1:2))>1.e-100
         det(G(1:2,1:2))
         disp(' Positive 2x2 determinant...'); 
         return;
     end
    
    [Q1,Q2] = GetLinesFromDegenerateConic(G);
    
    %   equation of the first  line is Q1(1)*x + Q1(2)*y + Q1(3) = 0
    %   equation of the second line is Q2(1)*x + Q2(2)*y + Q2(3) = 0
    
    %   Next find all the points of intersection of the above lines
    %   with the original conic (defined by the matrix C)
    
    %   The first line gives us up to two points
     
     k = 0;
     
     if abs(Q1(1)) > max(abs(Q1(2)),realmin)
         p = Q1(2)/Q1(1);  q = Q1(3)/Q1(1);
         a = A(1)*p*p-2*A(2)*p+A(3);
         b = A(2)*q-A(1)*p*q+A(4)*p-A(5);
         c = A(1)*q*q-2*A(4)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,2) = c/2/b;
             xy(k,1) = -p*xy(k,2) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,2) = (b + sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 else
                     k = k + 1;
                     xy(k,2) = (b - sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 end
             end
         end
     end

     if abs(Q1(2)) >= max(abs(Q1(1)),realmin)
         p = Q1(1)/Q1(2);  q = Q1(3)/Q1(2);
         a = A(3)*p*p-2*A(2)*p+A(1);
         b = A(2)*q-A(3)*p*q+A(5)*p-A(4);
         c = A(3)*q*q-2*A(5)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,1) = c/2/b;
             xy(k,2) = -p*xy(k,1) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,1) = (b + sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 else
                     k = k + 1;
                     xy(k,1) = (b - sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 end
             end
         end
     end

     %   The second line gives us up to two points

     if abs(Q2(1)) > max(abs(Q2(2)),realmin)
         p = Q2(2)/Q2(1);  q = Q2(3)/Q2(1);
         a = A(1)*p*p-2*A(2)*p+A(3);
         b = A(2)*q-A(1)*p*q+A(4)*p-A(5);
         c = A(1)*q*q-2*A(4)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,2) = c/2/b;
             xy(k,1) = -p*xy(k,2) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,2) = (b + sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 else
                     k = k + 1;
                     xy(k,2) = (b - sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 end
             end
         end
     end

     if abs(Q2(2)) >= max(abs(Q2(1)),realmin)
         p = Q2(1)/Q2(2);  q = Q2(3)/Q2(2);
         a = A(3)*p*p-2*A(2)*p+A(1);
         b = A(2)*q-A(3)*p*q+A(5)*p-A(4);
         c = A(3)*q*q-2*A(5)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,1) = c/2/b;
             xy(k,2) = -p*xy(k,1) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,1) = (b + sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 else
                     k = k + 1;
                     xy(k,1) = (b - sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 end
             end
         end
     end

     %       Select the intersection point closest to the point (u,v)
     if k > 0
         jmin = 1;
         dmin = realmax;
         for j=1:k
             if norm(xy(j,:) - XY(i,:))<dmin
                 dmin = norm(xy(j,:) - XY(i,:));
                 jmin = j;
             end
         end
         XYproj(i,:) = xy(jmin,:);   %    record the found projection
     else
         XYproj(i,:)=XY(i,:);
     end

     if ((~isfinite(XYproj(i,1)))||(~isfinite(XYproj(i,2))))
         XYproj(i,:) = XY(i,:);
     end
 end
 
end   %  ProjectPointsOntoConicByWEPinv