void ProjectPointsOntoParabola(Mnx1 X, Mnx1 Y, M5x1 ParG, Mnx1 &Xproj, Mnx1 &Yproj)
/* <-------- Input --------> <------ Output -------->
Projecting a given set of points onto a parabola
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
It is an adaptation of Eberly's method of projecting points onto ellipses
Input: two n-vectors X and Y containing the coordinates of n points
ParG is a 5x1 vector of the ellipse parameters:
ParG(0) is the x-coordinate of the parabola's vertex
ParG(1) is the y-coordinate of the parabola's vertex
ParG(2) is the distance from the focus to the directrix (often denoted by p)
ParG(3) is the angle of tilt of the ellipse
Note: the user needs to make sure that ParG(2) > 0.
Note: the equation of the hyperbola in the canonical coordinates is
y^2 = 2px
Output: two n-vectors Xproj and Yproj containing the coordinates of n projections of the points
The algorithm is proven to converge and reaches an accuracy of 14-15 significant digit
It takes 7-8 iterations per point, on average.
Nikolai Chernov, August 2012
*/
{
int n=X.rows(),i,iter,iterMax;
const reals s2=sqrt(Two),TwoThirds=Two/Three;
reals p,pp,psqrt,c,s,u,v,uu,vv,T,Tnew,F,xprojx,yprojx,xprojy,yprojy,Fx,Fy;
Mnx1 X0(n,1), Y0(n,1); // column vectors of length n (matrices n times 1)
p = ParG(2); pp = p*p; psqrt = sqrt(p);
s = sin(ParG(3)); c = cos(ParG(3)); // rotation parameters
X0 = (X.array()-ParG(0))*c + (Y.array()-ParG(1))*s; // rotating x-coordinates
Y0 = (Y.array()-ParG(1))*c - (X.array()-ParG(0))*s; // rotating y-coordinates
iterMax = 100; // Maximum number of Newton's ietrations. Usually, 5-7 are enough
for (i=0; i