function XYproj = ProjectPointsOntoConic3(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: 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]; E = eig(D,C); % find generalized eigenvalues of the pencil (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...'); coordinates of the point to be projected return; end G = D-E(jreal)*C; % parameters of the degenerate conic % in the family of conics D-t*C if det(G(1:2,1:2))>eps det(G(1:2,1:2)) disp(' Positive 2x2 determinant...'); coordinates of the point to be projected return; end [Q F] = eig(G); % eigendecomposition of the G matrix % The degenerate conic corresponding to the matrix G % consists of two lines. Next we find their equations [Fmax,Imax] = max(diag(F)); [Fmin,Imin] = min(diag(F)); Qmax = Q(1:3,Imax)*sqrt(abs(Fmax)); Qmin = Q(1:3,Imin)*sqrt(abs(Fmin)); Q1 = Qmax + Qmin; Q2 = Qmax - Qmin; % 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)),eps) 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); 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 if abs(Q1(2)) >= max(abs(Q1(1)),eps) 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); 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 % The second line gives us up to two points if abs(Q2(1)) > max(abs(Q2(2)),eps) 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); 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 if abs(Q2(2)) >= max(abs(Q2(1)),eps) 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); 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 % Select the intersection point closest to the point (u,v) jmin = 1; dmin = norm(xy(1,:) - XY(i,:)); for j=2:k if norm(xy(j,:) - XY(i,:))