projPointOnEllipse

PROJPOINTONELLIPSE Project a point orthogonally onto an ellipse.

   PROJ = projPointOnEllipse(PT, ELLI)
   Computes the (orthogonal) projection of point PT onto the ellipse given
   by ELLI.


   Example
     % create an ellipse and a point
     elli = [50 50 40 20 30];
     pt = [60 10];
     % display reference figures
     figure; hold on; drawEllipse(elli, 'b');
     axis equal; axis([0 100 0 100]);
     drawPoint(pt, 'bo');
     % compute projected point
     proj = projPointOnEllipse(pt, elli);
     drawEdge([pt proj], 'b');
     drawPoint(proj, 'k*');

   See also
     ellipses2d, distancePointEllipse, projPointOnLine, ellipsePoint
     drawEllipse

0001 function proj = projPointOnEllipse(point, elli)
0002 %PROJPOINTONELLIPSE Project a point orthogonally onto an ellipse.
0003 %
0004 %   PROJ = projPointOnEllipse(PT, ELLI)
0005 %   Computes the (orthogonal) projection of point PT onto the ellipse given
0006 %   by ELLI.
0007 %
0008 %
0009 %   Example
0010 %     % create an ellipse and a point
0011 %     elli = [50 50 40 20 30];
0012 %     pt = [60 10];
0013 %     % display reference figures
0014 %     figure; hold on; drawEllipse(elli, 'b');
0015 %     axis equal; axis([0 100 0 100]);
0016 %     drawPoint(pt, 'bo');
0017 %     % compute projected point
0018 %     proj = projPointOnEllipse(pt, elli);
0019 %     drawEdge([pt proj], 'b');
0020 %     drawPoint(proj, 'k*');
0021 %
0022 %   See also
0023 %     ellipses2d, distancePointEllipse, projPointOnLine, ellipsePoint
0024 %     drawEllipse
0025 %
0026 
0027 % ------
0028 % Author: David Legland
0029 % E-mail: david.legland@inrae.fr
0030 % Created: 2022-07-17, using Matlab 9.12.0.1884302 (R2022a)
0031 % Copyright 2022-2024 INRAE - BIA Research Unit - BIBS Platform (Nantes)
0032 
0033 % defaults
0034 nMaxIters = 4;
0035 
0036 % compute transform to centered axis-aligned ellipse
0037 center = elli(1:2);
0038 theta = deg2rad(elli(5));
0039 transfo = createRotation(-theta) * createTranslation(-center);
0040 pointT = transformPoint(point, transfo);
0041 
0042 % retrieve ellipse semi axis lengths
0043 a = elli(3);
0044 b = elli(4);
0045 
0046 % keep absolute values
0047 px = abs(pointT(:,1));
0048 py = abs(pointT(:,2));
0049 
0050 % initial guess of solution
0051 tx = 0.707;
0052 ty = 0.707;
0053 
0054 % iterate
0055 for i = 1:nMaxIters
0056     x = a * tx;
0057     y = b * ty;
0058 
0059     ex = (a*a - b*b) * power(tx, 3) / a;
0060     ey = (b*b - a*a) * power(ty, 3) / b;
0061 
0062     rx = x - ex;
0063     ry = y - ey;
0064 
0065     qx = px - ex;
0066     qy = py - ey;
0067 
0068     r = hypot(ry, rx);
0069     q = hypot(qy, qx);
0070 
0071     tx = min(1, max(0, (qx .* r ./ q + ex) / a));
0072     ty = min(1, max(0, (qy .* r ./ q + ey) / b));
0073     t = hypot(ty, tx);
0074     tx = tx ./ t;
0075     ty = ty ./ t;
0076 
0077 end
0078 
0079 % computes coordinates of projection
0080 projX = a * tx;
0081 projY = b * ty;
0082 
0083 % fix sign
0084 projX(pointT(:,1) < 0) = -projX(pointT(:,1) < 0);
0085 projY(pointT(:,2) < 0) = -projY(pointT(:,2) < 0);
0086 
0087 % project pack to basis of original ellipse
0088 proj = transformPoint([projX projY], inv(transfo));