I am trying to find the farthest and closest points of a ellipse without using any brute force type of coding. The processing power is limited so it should be as pinpoint as possible. I have tried a binary search method but it should still be faster than that.
All values of ellipse is known (Center Point(C), Rotation(r), focal points, long side's distance(a), short side's distance(b)) as well as the point outside of the ellipse(P).
Thanks!
Rewritten on 2016-04-29:
Let $\vec{c} = (x_c, y_c)$ be the center of the ellipse, $\vec{a} = (x_a, y_a)$ its semi-major axis, and $\vec{b} = (x_b, y_b)$ its semi-minor axis, $\lVert\vec{a}\rVert\ge\lVert\vec{b}\rVert$. Note that the axes are orthogonal, $\vec{a}\cdot\vec{b} = 0$.
To simplify the situation, let's switch to a coordinate system where the center of the ellipse is at origin, semi-major axis is parallel to the $x$ axis, and semi-minor axis to the $y$ axis. In this coordinate system, the point $\vec{p} = (x_p, y_p)$ is at $$\begin{cases} x = (\vec{p} - \vec{c}) \cdot \frac{\vec{a}}{\lVert\vec{a}\rVert} \\ y = (\vec{p} - \vec{c}) \cdot \frac{\vec{b}}{\lVert\vec{b}\rVert} \end{cases}$$ We only need to consider the positive quadrant of the ellipse, because we can negate both $x$ and $a$ if $x \lt 0$; and negate $y$ and $b$ if $y \lt 0$.)
Using $r_x = \lVert\vec{a}\rVert$ and $r_y = \lVert\vec{b}\rVert$, we can parametrize the positive quadrant of the ellipse $t = x / r_x$, i.e. $$\begin{cases} x_{+}(t) = r_x \; t \\ y_{+}(t) = r_y \; \sqrt{1 - t^2} \end{cases}$$
Point $t$ in the positive quadrant is closest to our point $(x, y)$ when the distance squared $s_N(t)$ reaches a minimum: $$s_N(t) = \left ( (x - x_{+}(t))^2 + (y - y_{+}(t))^2 \right ) = (x - t r_x)^2 + (y - r_y \sqrt{1 - t^2})^2$$ It reaches a minimum when $t = 0$, $t = 1$, or $ds_N(t)/dt = 0$. The roots of the last one are the same as the roots of the quartic equation $$f(t) = T_4 t^4 + T_3 t^3 + T_2 t^2 + T_1 t + T_0 = 0$$ where $$\begin{align} T_4 &= (r_x^2 - r_y^2)^2 \\ T_3 &= 2 x r_x ( r_y^2 - r_x^2 ) \\ T_2 &= x^2 r_x^2 + y^2 r_y^2 - (r_x^2 - r_y^2)^2 \\ T_1 &= 2 x r_x ( r_x^2 - r_y^2 ) \\ T_0 &= - x^2 r_x^2 \end{align}$$ Note that you can evaluate it as $f(t) = T_0 + t (T_1 + t(T_2 + t(T_3 + t T_4)))$, and that $0 \le t \le 1$. It typically has one or three roots within that range, so you probably want to use a method that also examines its derivative (which is, of course, $T_1 + t(2 T_2 + t(3 T_3 + t 4 T_4))$). Note that this is a simple polynomial, and thus very fast to evaluate. Compute all $s_N(t)$ for $t = 0$, $t = 1$, and for $f(t) = 0, 0 \lt t \lt 1$, and pick the $t$ that yields the smallest $s_N(t)$. That point $t$ on the positive quadrant of the ellipse is the nearest to our point $(x, y)$: $$\vec{p}_{nearest} = \vec{c} + \vec{a} t + \vec{b} \sqrt{1 - t^2} = \begin{cases} x_{nearest} = x_c + x_a t + x_b \sqrt{1 - t^2} \\ y_{nearest} = y_c + y_a t + y_b \sqrt{1 - t^2} \end{cases}$$
The point furthest on the ellipse must be in the opposite quadrant. The point $t$ in the quadrant opposite to point $(x,y)$ is $$\begin{cases} x_{-}(t) = - r_x \; t \\ y_{-}(t) = - r_y \; \sqrt{1 - t^2} \end{cases}$$ and the distance squared is $$s_F(t) = (x - x_{-}(t))^2 + (y - y_{-}(t))^2 = (x + t r_x)^2 + (y + r_y \sqrt{1 - t^2})^2$$ If $s_F(t)$ is a global maximum $0 \le t \le 1$, then the point on the ellipse furthest from point $(x,y)$ is $$\vec{p}_{furthest} = \vec{c} - \vec{a} t - \vec{b} \sqrt{1 - t^2} = \begin{cases} x_{furthest} = x_c - x_a t - x_b \sqrt{1 - t^2} \\ y_{furthest} = y_c - y_a t - y_b \sqrt{1 - t^2} \end{cases}$$ $s_F(t)$ reaches a maximum when $ds_F(t)/dt = 0$, i.e. when $$g(t) = T_4 t^4 + T_3 t^3 + T_2 t^2 + T_1 t + T_0 = 0$$ where $$\begin{align} T_4 &= (r_x^2 - r_y^2)^2 \\ T_3 &= 2 x r_x ( r_x^2 - r_y^2 ) \\ T_2 &= x^2 r_x^2 + y^2 r_y^2 - (r_x^2 - r_y^2)^2 \\ T_1 &= 2 x r_x ( r_y^2 - r_x^2 ) \\ T_0 &= - x^2 r_x^2 \end{align}$$ Note that the coefficients are the same as for the nearest case, except with $T_3$ and $T_1$ negated. Again, of all the values of $t$ where $g(t) = 0, 0 \lt t \lt 1$, and of $t = 0$ and $t = 1$, pick the $t$ where $s_F(t)$ is the largest, to find the point on the ellipse furthest from point $(x, y)$.
Implementation-wise, the only open question is how to solve the quartic equations. For these coefficients, there do exist analytical solutions that can be found using e.g. Maple, but they are computationally intensive (tested 3,500 to 10,000 clock cycles on Intel Core i5-4200U using double precision floating-point arithmetic).
Here is a test program, that takes a Xorshift* pseudorandom number generator seed, generates a test ellipse, and outputs an SVG image for easy verification:
Note that complex arithmetic is required for the intermediate values. To make it a bit faster, I split the real-only case out, so that complex intermediates are only used when really necessary. (The intermediates are not purely imaginary either at any point other than immediately after taking the square root of a negative real.)
To avoid the cases where the root has a small imaginary part due to numerical inaccuracy, I introduced the
IEPSmacro. In this particular case, it can be rather large, as the roots are tested anyway; as long as it is at least as large as the largest expected numerical inaccuracy in the imaginary part, the code should find the points correctly.Evaluating the quartic function to be solved, as well as its derivative, should only take a few clock cycles as they are simple polynoms. If limited accuracy is also desired, it may well be better to solve the quartic functions ($C_0 + C_1 t + C_2 t^2 - C_1 t^3 + t^4 = 0$ and $C_0 - C_1 t + C_2 t^2 + C_1 t^3 + t^4 = 0$) numerically.
I suspect Newton's method is not suitable, as it tends to overshoot. Personally, I think I'd split the range into a small number of spans, and use a binary search to find the root $t$ within each span, if the span spans zero, or bulges towards zero. (Since $t = x/\sqrt{x^2 + y^2}$ is the solution for the circle case, it is also an interesting point in the general case here.)