function XYproj = ProjectPointsOntoConic4(XY,A)

%   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:  XYproj is the array of coordinates of projections
%
%      The problem of projecting a given point (u,v) onto a conic
%      is known to reduce to a polynomial equation of degree 4 (a quartic equation)

%      Here we introduce variable t by the following equations:
%           u-x = t*m     v-y = t*n
%      where (m,n) is the normal vector to the conic at the point (x,y)
%      The above equation is based on the orthogonality of the projection
%
%      Then we construct a quartic equation for the variable t
%      The quartic equation is solved by MATLAB function "roots"

%      This gives up to four solutions
%      Each solution gives a point (x,y) on the conic
%      The one closest to (u,v) is the desired projection

%        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);
R0 = zeros(4,1);

for i=1:n
     u = XY(i,1);  v = XY(i,2);   %  coordinates of the point to be projected
     
     L0 = [1; -(A(1)+A(3)); A(1)*A(3)-A(2)*A(2)];
     Lx = [u; A(2)*v-A(3)*u+A(4); A(2)*A(5)-A(3)*A(4)];
     Ly = [v; A(2)*u-A(1)*v+A(5); A(2)*A(4)-A(1)*A(5)];
     
     M = A(1)*Lx*Lx' + 2*A(2)*Lx*Ly' + A(3)*Ly*Ly' + 2*A(4)*Lx*L0' + 2*A(5)*Ly*L0' + A(6)*L0*L0';
     
%      solve the quartic equation

     R = roots([M(1,1), M(1,2)+M(2,1), M(1,3)+M(2,2)+M(3,1), M(2,3)+M(3,2), M(3,3)]);
    
%      extrct real roots

     k = 0;
     for j=1:size(R,1)
         if isreal(R(j))
             k = k + 1;
             R0(k) = R(j);
         end
     end
     U = [R0(1:k).*R0(1:k) R0(1:k) ones(k,1)];
     
%       compute the coordinates of up to four points on the conic

     xy = [(U*Lx)./(U*L0) (U*Ly)./(U*L0)];
     
%       select the one closest to the point (u,v)
     
     jmin = 1;  dmin = norm(xy(1,:) - XY(i,:));
     for j=2: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

end

end   %  ProjectPointsOntoConic4