function XYproj = ProjectPointsOntoConicByQuarticEquation(XY,A)
% Quartic equation method
% Projecting a given set of points onto a quadratic curve (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,:))