function XYproj = ProjectPointsOntoConicByWEPqz(XY,A) % % The WEP method % version: solving a generalized e-value problem with 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: XYproj is the array of coordinates of projections % % Nikolai Chernov, August 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); C0 = [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]; C = C0; if abs(det(C)) > abs(det(D)) G = C; C = D; D = G; end E = eig(C,D); % next find a real generalized eigenvalue (just one is enough) jreal=0; jreals=0; for j=1:3 if isreal(E(j)) jreal=j; jreals = jreals+1; end end if jreals==0 disp(' No real e-values...'); return; end if jreals<=2 x = E(jreal); end if jreals==3 x1 = E(1); x2 = E(2); x3 = E(3); x12 = abs(x1-x2); x13 = abs(x1-x3); x23 = abs(x2-x3); if ((x12<=x13)&&(x12<=x23)) x = x3; end; if ((x13<=x12)&&(x13<=x23)) x = x2; end; if ((x23<=x12)&&(x23<=x13)) x = x1; end; end G = C-x*D; % parameters of the degenerate conic % in the family of conics C-t*D 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,:))