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;
    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
end

    %         compute the projection point in the proper quadrant
    %         and mirror...
    
for i=1:n

    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