function XYproj = ProjectPointsOntoHyperbola(XY,ParG) % % Projecting a given set of points onto a hyperbola % 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 hyperbola parameters % ParG = [Center(1:2), Axes(1:2), Angle]' % Center - the coordinates of the hyperbola's center % Axes - the axes (more precisely: % Axes(1) is the distance from the center to each vertex % Axes(2) is the vertical distance from each vertex to asymptotes % Angle - the angle of tilt of the hyperbola % Note: if hyperbola is defined by equation x^2/a^2-y^2/b^2=1, % then Center=(0,0),Axes=(a,b),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 Center = ParG(1:2); a = ParG(3); b = ParG(4); Angle = ParG(5); % hyperbola's parameters aa = a^2; bb = b^2; % squares of the axes D = aa + bb; % auxiliary parameter s2 = sqrt(2); if (a <= 0)||(b <= 0) disp(' Axes of the hyperbola must be positive') return end Q = [cos(Angle) -sin(Angle); sin(Angle) cos(Angle)]; % Q is the matrix for rotating the points and the hyperbola to the canonical system n = size(XY,1); % n is the number of points to be projected XYproj = zeros(n,2); % prepare an array of projections XY0 = [XY(:,1)-Center(1) XY(:,2)-Center(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 = abs(XY0(i,1)); v = abs(XY0(i,2)); % coordinates of the point au = a*u; bv = b*v; % products of coordinates and axes F0 = u*u/aa - v*v/bb - 1; % the F-value at the point (u,v) if F0 == 0 % point on the hyperbola already XYproj(i,:) = XY0(i,:); continue; end Finf = (sqrt(au) + sqrt(bv))^2*(au - bv)/D^2 - 1; % the F-value at the inflection point if F0 > 0 if Finf > 0 T = max(D*bv/(au+bv),D-au); if T <= 0 T = 0; else for iter=1:iterMax % loop of Newton's iterations to solve F(T)=0 F = (au/(D-T))^2 - (bv/T)^2 - 1; % value of F at the current T if (F >= 0)||(T >= bb) % gone too far, emergency stop break; end; Fder = 2*((au/(D-T))^2/(D-T) + (bv/T)^2/T); % 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 else T = bb; for iter=1:iterMax % loop of Newton's iterations to solve F(T)=0 F = (au/(D-T))^2 - (bv/T)^2 - 1; % value of F at the current T if (F <= 0)||(T <= 0) % gone too far, emergency stop break; end; Fder = 2*((au/(D-T))^2/(D-T) + (bv/T)^2/T); % 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 <= a+bb/a; % the point is not too far from the origin xprojx = a; % projects onto the vertex of the hyperbola yprojx = 0; xprojy = xprojx; yprojy = yprojx; else % the point is too far from the origin xprojx = u*aa/(aa+bb); % projection is above the horizontal axis yprojx = b*sqrt(max((xprojx/a)^2-1,0)); xprojy = xprojx; yprojy = yprojx; end else % the generic case xprojx = u*aa/(D-T); % first candidate for projection yprojx = b*sqrt(max((xprojx/a)^2-1,0)); yprojy = v*bb/T; % second candidate for projection xprojy = a*sqrt(1+(yprojy/b)^2); end else if Finf < 0 if au <= 0 T = 0; else T = min(D*au/(au+s2*bv),au/s2); end if T <= 0 T = 0; else for iter=1:iterMax % loop of Newton's iterations to solve F(T)=0 F = (au/T)^2 - (bv/(D-T))^2 - 1; % value of F at the current T if (F <= 0)||(T >= aa) % gone too far, emergency stop break; end; Fder = -2*((au/T)^2/T + (bv/(D-T))^2/(D-T)); % 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 else T = aa; for iter=1:iterMax % loop of Newton's iterations to solve F(T)=0 F = (au/T)^2 - (bv/(D-T))^2 - 1; % value of F at the current T if (F >= 0)||(T <= 0) % gone too far, emergency stop break; end; Fder = -2*((au/T)^2/T + (bv/(D-T))^2/(D-T)); % 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 vertical axis yprojx = v*bb/D; xprojx = a*sqrt(1+(yprojx/b)^2); xprojy = xprojx; yprojy = yprojx; else % the generic case xprojx = u*aa/T; % first candidate for projection yprojx = b*sqrt(max((xprojx/a)^2-1,0)); yprojy = v*bb/(D-T); % second candidate for projection xprojy = a*sqrt(1+(yprojy/b)^2); end 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,1)<0), XYproj(i,1) = -XYproj(i,1); end; if (XY0(i,2)<0), XYproj(i,2) = -XYproj(i,2); end; end % end the main loop over the points % rotate and shift back to the original system XYproj = XYproj*Q'; % rotation XYproj(:,1) = XYproj(:,1)+Center(1); % shifting XYproj(:,2) = XYproj(:,2)+Center(2); % shifting end % ProjectPointsOntoHyperbola