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,:))