function XYproj = ProjectPointsOntoParabola(XY,ParG) % Projecting a given set of points onto a parabola % and computing the distances from the points to the hyperbola % % This method is proposed by Nikolai Chernov and Hui Ma in their article % "Least squares fitting of quadratic curves and surfaces" % published in book % "Computer Vision", Editor S. R. Yoshida, Nova Science Publishers 2011; pp. 285-302 % % Input: XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2) % ParG is a vector 5x1 of the parabola parameters % ParG = [Vertex(1:2), p, Angle]' % Vertex - the coordinates of the parabola's vertex % p - distance from the focus to the directrix % Angle - the angle of tilt of the parabola % Note: if parabola is defined by equation y^2=2px, % then Vertex=(0,0), p is the factor in front of x, Angle=0 % % Output: XYproj is the array of coordinates of projections % % The algorithm is proven to converge and reaches an accuracy of 10-12 significant digit % It takes 6-8 iterations per point, on average % % Nikolai Chernov, March 2012 Vertex = ParG(1:2); p= ParG(3); Angle = ParG(4); % parabola's parameters if p <= 0 disp(' The p parameter of the parabola must be positive') return end n = size(XY,1); % n is the number of points to be projected XYproj = zeros(n,2); % prepare an array of projections pp=p*p; psqrt = sqrt(p); s2 = sqrt(2); Q = [cos(Angle) -sin(Angle); sin(Angle) cos(Angle)]; % Q is the matrix for rotating the points and the parabola to the canonical system XY0 = [XY(:,1)-Vertex(1) XY(:,2)-Vertex(2)]*Q; % points in canonical coordinates iterMax = 100; % Maximum number of Newton's ietrations. Usually, 6-8 are enough for i=1:n % Main loop over the points to be projected u = XY0(i,1); v = abs(XY0(i,2)); % coordinates of the point uu = u*u; vv = v*v; pu=p*u; if u <= p T = (v/p/s2)^(2/3); else T = min(v/2/sqrt(u-p)/psqrt,(v/p/2)^(2/3)); end if T <= 0 % safeguard against singular cases T = 0; else for iter=1:iterMax % loop of Newton's iterations to solve F(T)=0 F = (v/T)^2 - 2*p*(u-p+p*T); % value of F at the current T if F <= 0 % gone too far, emergency stop break; end; Fder = -2*((v/T)^2/T + pp); % derivative of F with respect to T Step = F/Fder; if T == T - Step % no progress, terminate iterations break; end T = T - Step; % Newton's iteration end % end of the loop of Newton's iterations end if T <= 0 % singularity: the point is on the horizontal axis if u <= p % point is to the left of the curvature center xprojx = 0; % project to the vertex (0,0) yprojx = 0; else % point is to the right of the curvature center xprojx = u-p; % project above the axis if xprojx <= 0 % safeguard against singular cases xprojx = 0; end yprojx = sqrt(2*p*xprojx); end xprojy = xprojx; yprojy = yprojx; else % the generic case xprojx = u + p*(T-1); % first candidate for projection if xprojx <= 0 xprojx = 0; end yprojx = sqrt(2*p*xprojx); yprojy = v/T; % second candidate for projection xprojy = yprojy^2/p/2; end % compute the projection of the point onto the hyperbola Fx = (xprojx-u)^2 + (yprojx-v)^2; % "quality" of first candidate Fy = (xprojy-u)^2 + (yprojy-v)^2; % "quality" of second candidate if Fx < Fy % the first candidate is better XYproj(i,1) = xprojx; XYproj(i,2) = yprojx; else % the second candidate is better XYproj(i,1) = xprojy; XYproj(i,2) = yprojy; end % end comparing the two candidates % compute the projection point in the proper quadrant if (XY0(i,2)<0), XYproj(i,2) = -XYproj(i,2); end; end % end the main loop XYproj = XYproj*Q'; % rotation XYproj(:,1) = XYproj(:,1)+Vertex(1); % shifting XYproj(:,2) = XYproj(:,2)+Vertex(2); % shifting end % ProjectPointsOntoParabola