Closest point on ellipse to given point implementation in Octave
John Bryan
-
For two cases of ellipses
- Circle
- Ellipse centered at origin with the major and minor axes on the x and y axes.
the calculation to find the nearest point on the ellipse to a given point is implemented in Octave.
-
- Case 1: Circle.
- r = radius of circle
- point $D = (D_x,D_y) = (h,k)$, the center of circle
- point $P = (P_x,P_y) = (c,d)$, given point inside or outside of circle
- point $E = (E_x,E_y) = (x,y)$,closest point on circle to $P$
- vector $\vec{D} =$ vector from origin to point $D$.
- vector $\vec{P} =$ vector from origin to point $P$.
- vector $\vec{E} =$ vector from origin to point $E$.
Given $D$ and $P$, find $E$. From [1], $E$ will be the point on the circle that lies on the vector
connecting $D$ and $P$.
\begin{equation}
\vec{E} = \vec{D} + r{{\vec{P} - \vec{D}} \over {\Vert{\vec{P} - \vec{D}}\Vert}}.
\tag{eq. 1}
\end{equation}
$$
\begin{align}
E_x & = D_x + r{{P_x-D_x} \over \sqrt{{(P_x - D_x)}^2 + {(P_y - D_y)}^2}} \\
& = h + r{{c-h} \over \sqrt{{(c - h)}^2 + {(d - k)}^2}}.
\tag{eq. 2}
\end{align}
$$
$$
\begin{align}
E_y & = D_y + r{{P_y-D_y} \over \sqrt{{(P_x - D_x)}^2 + {(P_y - D_y)}^2}} \\
& = k + r{{d-k} \over \sqrt{{(c - h)}^2 + {(d - k)}^2}}.
\tag{eq. 3}
\end{align}
$$
- Case 2: Ellipse centered at origin with the major and minor axes on the x and y axes.
- a = horizontal semi-axis length
- b = vertical semi-axis length
- P = (c,d) = given point outside of ellipse
- E = (x,y) = closest point on ellipse to P
From [2], if one considers a bounding circle around the given point (c,d), which passes through the nearest point (x,y) on the ellipse, if one draws a line from (c,d) to (x,y), it must be perpendicular to the shared tangent of the circle and the ellipse. Any other point on the ellipse is outside of the bounding circle, and so if further away from the given point (c,d). So, at (x,y), the product of the gradient of tangent, $m_T$, and gradient of the line connecting (c,d) and (x,y), $m_l$, equals -1. $m_Tm_l=-1$.
\begin{equation}
\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1.
\tag{eq. 4}
\end{equation}
\begin{equation}
y^2=b^2(1-\frac{x^2}{a^2})
\tag{eq. 5}
\end{equation}
\begin{equation}
y= g(x) = \begin{cases}
b\sqrt{1-(x/a)^2} & \text{if y positive or zero} \\[2ex]
-b\sqrt{1-(x/a)^2} & \text{if y negative}
\end{cases}
\tag{eq. 6}
\end{equation}
\begin{equation}
\text{Let}\quad p(x)= b\sqrt{1-(x/a)^2}
\tag{eq. 7}
\end{equation}
\begin{equation}
y= g(x) = \begin{cases}
p(x) & \text{if y positive or zero} \\[2ex]
-p(x) & \text{if y negative}
\end{cases}
\tag{eq. 8}
\end{equation}
Differentiating eq.5 to find the gradient of the tangent, $m_T$, we have,
\begin{equation}
2yy^\prime = -2(b/a)^2x
\end{equation}
\begin{equation}
m_T = y^\prime = -\frac{b^2}{a^2}\frac{x}{y}
\tag{eq. 9}
\end{equation}
The equation for the gradient of line, $m_l$:
\begin{equation}
m_l = \frac{y-d}{x-c}
\tag{eq. 10}
\end{equation}
Condition for perpendicular lines - product of gradients = -1.
$$m_T \cdot m_l = -1.$$
$$-\frac{b^2x}{a^2y} \cdot \frac{y-d}{x-c} = -1.$$
$$b^2x(y-d) = a^2y(x-c)$$
$$b^2xy-b^2xd = a^2yx-a^2yc$$
$$b^2xy - a^2yx + a^2yc = b^2xd$$
$$y(b^2x-a^2x+a^2c)a=b^2xd$$
$$y(x(b^2-a^2)+a^2c)=(b^2d)x$$
$$y = h(x) = \frac{(b^2d)x}{(b^2-a^2)x+a^2c}
\tag{eq: 11}
$$
For $y=y$, equate h(x) from eq.11 with g(x) from eq.8,
\begin{equation}h(x)=g(x)=\begin{cases}
p(x) & \text{if y positive or zero} \\[2ex]
-p(x) & \text{if y negative}
\end{cases}
\tag{eq: 12}
\end{equation}
Find x for $$ f(x) = h(x) - g(x) = 0 $$
Use Newton iteration method, $x_{n+1}=x_n-(f(x_n)/f^\prime(x_n))$, to find x.
For the initial x value, use the x value from eq.2, with $r=max(a,b)$.
\begin{equation}
f(x) = h(x) - g(x)
= \begin{cases}
h(x) - p(x) & \text{if y positive or zero} \\[2ex]
h(x) + p(x) & \text{if y negative} \\[2ex]
\end{cases}
\tag{eq: 13}
\end{equation}
\begin{equation}
f^\prime(x) = h^\prime(x) - g^\prime(x)
= \begin{cases}
h^\prime(x) - p^\prime(x) & \text{if y positive or zero} \\[2ex]
h^\prime(x) + p^\prime(x) & \text{if y negative} \\[2ex]
\end{cases}
\tag{eq: 14}
\end{equation}
$$
\begin{align} h^\prime(x) & = \frac{b^2d((b^2-a^2)x+a^2c)-b^2d(b^2-a^2)x}{((b^2-a^2)x+a^2c)^2} \\
& = \frac{a^2b^2cd}{((b^2-a^2)x+a^2c)^2}
\end{align}
\tag{eq: 15}
$$
$$
\begin{align} p^\prime(x) &= (b)(\frac{1}{\sqrt{1-(x/a)^2}})(-\frac{2x}{a^2}) \\
& = \frac{-2bx}{a^2\sqrt{1-(x/a)^2}} \end{align}
\tag{eq: 16}
$$
Once x is calculated from the Newton iteration method, find y from eq.8.
-
Octave implementation.
-
Case 1 Results
-
Case 2 Results
-
References
-
Math StackExchange
-
Stack Overflow
