Find projection point on elipse

Question

Given the ellipse

$$ E(x, \ y) = ax^2 + bxy + cy^2 + dx + ey + f = 0 $$

With $4ac - b^2 = 1$

How can I compute the projection $(x_p, \ y_p)$ of the point $(x_0, y_0)$ in this ellipse?

$$ E(x_p, \ y_p) = 0 $$

Motivation: There is already a direct algorithm to fit ellipses from a set of datapoints as shown here where I can get the coefficients.

In linear interpolation, there's the R-Squared indicator bellow that tells us how good the model is. With $R_2 = 1$, all the points are in the line.

$$ R_2 = \dfrac{\left[n\sum x_iy_i - \left(\sum x_i\right)\left(\sum y_i\right) \right]^2}{\left[n\sum x_i^2 - \left(\sum x_i\right)^2\right]\left[n\sum y_i^2 - \left(\sum y_i\right)^2\right]} $$

But for the ellipse problem I don't know the formula of $R_2$. Then I want to compute the mean distance between the points and the ellipse to know how well my datapoints fits an ellipse.

To compute the distance of each point I need the projection point.

Answer 1

From the algebraic equation

$ a x^2 + b xy + c y^2 + d x + e y + f = 0 $

Obtain the vector parametric equation

$ r = C + E_1 \cos(t) + E_2 \sin(t) $

where $r =[x,y]^T$ is the coordinate vector, $C$ is the center, $E_1, E_2$ are the semi-axes vectors (minor and major), they are all $2 \times 1$ vectors.

The derivative of $r$ with respect to $t$ gives a tangent vector to the ellipse.

$ \dot{r} = - E_1 \sin(t) + E_2 \cos(t) $

Now we want $r$ such that this tangent vector is perpendicular to the vector $r - P$, where $P=(x_0, y_0) $ is the point to be projected onto the ellipse.

Hence, we now have

$ (r - P) \cdot \dot{r} = 0 $

which in expanded form is

$ \left(C - P + E_1 \cos(t) + E_2 \sin(t) \right) \cdot \left(-E_1 \sin(t) + E_2 \cos(t) \right) = 0$

And this leads directly to

$ C_1 \cos(t) + C_2 \sin(t) + C_3 \cos(2t) + C_4 \sin(2t) = 0 \hspace{30pt}(*) $

where

$C_1 = (C-P)\cdot E_2 , \\C_2 = -(C-P) \cdot E_1,\\ C_3 = E_1 \cdot E_2 ,\\ C_4 = \frac{1}{2} \left(E_2 \cdot E_2 - E_1 \cdot E_1 \right) $

Note that since $E_1$ is normally perpendicular to $E_2 $ as they represent the semi-minor and semi-major axes , then $C_3 = 0 $

Equation $(*)$ can be solved exactly for its roots $\{t_i\}$ which are a maximum of $4$ in number. To calculate the projection of $P$ we have to compute

$ P_i = C + E_1 \cos(t_i) + E_2 \sin(t_i) $

And then we find the distances $\| P_i - P \| $ for all the possible points $P_i$, and choose the minimum of these. The minimizing $P_i$ is the projection of $P$.

Answer 2

This answer is not the one that I expected but it solves my problem in a non-elegant way.

I will divide this answer in three parts:

1. Transform the ellipse with center $(x_c, \ y_c)$ and rotated counter clockwise with angle $\varphi$

   $$E(x, y) = ax^2+bxy+cy^2 + dx + ey+f = 0 \ \ \ \ \ \ \ \ \ \ \ \ \ (1)$$

   into another ellipse

   $$F(u , v) = \left(\dfrac{u}{A}\right)^2 + \left(\dfrac{v}{B}\right)^2 - 1 =0 \ \ \ \ \ \ \ \ \ \ \ \ \ (2)$$

   With $A > B$

2. Scale into a new ellipse of major axis equal 1

   $$G(a, b) = a^2 + \left(\dfrac{b}{H}\right)^2 - 1 \ \ \ \ \ \ \ \ \ \ \ \ \ (3)$$

3. Compute the new point $(u_0, \ v_0)$ using the translation and rotations operations of point $(x_0, y_0)$

4. Compute $(a_0, \ b_0)$ by scaling $(u_0, \ v_0)$

5. Compute $(a, \ b)$ such we get the minimal distance from $(a_0, \ b_0)$ to $(a, \ b)$

6. Scale, rotate and translate back

Part 1: Ellipse transformation

We can describe the variable $x$ and $y$ using

$$ x = u \cdot \cos \varphi + v \cdot \sin \varphi + x_c \ \ \ \ \ \ \ \ \ \ \ \ \ (4) $$ $$ y = -u \cdot \sin \varphi + v \cdot \cos \varphi + y_c \ \ \ \ \ \ \ \ \ \ \ \ \ (5) $$

With $(1)$, $(2)$, $(4)$ and $(5)$ we can solve it. The computations are a little bit hard but we get

The center

$$ x_c = be-2cd $$ $$ y_c = bd-2ae $$

The angle

$$ \varphi = \dfrac{1}{2}\arctan\left(\dfrac{b}{c-a}\right) $$

The major and minor axis

$$ \kappa = \sqrt{(a-c)^2+b^2} $$ $$ \eta = \left(f+bde\right) -\left(ae^2+ cd^2\right) $$ $$ A = \sqrt{\dfrac{-\eta}{a+c-\kappa}} $$ $$ B = \sqrt{\dfrac{-\eta}{a+c+\kappa}} $$

Part 2: Scalling ellipse

$$ u = A \cdot a \ \ \ \ \ \ \ \ \ v = A \cdot b \ \ \ \ \ \ \ \ \ H = \dfrac{B}{A} $$

Part 3: Point transformation

To make a translation we only subtract the center $(x_c, \ y_c)$ and to make a clockwise rotation we put

$$ \begin{bmatrix} u_0 \\ v_0 \end{bmatrix} = \begin{bmatrix} \cos \varphi & -\sin \varphi \\ \sin \varphi & \cos \varphi \end{bmatrix}\begin{bmatrix} x_0 - x_c \\ y_0 - y_c \end{bmatrix} $$

Part 4: Scale point

Then we have

$$ a_0 = \dfrac{u_0}{A} \ \ \ \ \ \ \ \ \ \ \ \ \ b_0 = \dfrac{v_0}{A} $$

Part 5: Compute distance

Let's find the distance of an arbitrary point $(a_0, b_0)$ to the ellipse $(3)$

The distance square is given by

$$ D^2 = (a-a_0)^2 + (b-b_0)^2 $$

Option 3.1: Describe $a$ and $b$ in terms of $\theta$

Let

$$ a = \cos \theta \ \ \ \ \ \ \ \ \ \ \ \ \ \ b = H \cdot \sin \theta $$

Then

$$ D^2 = (\cos \theta - a_0)^2 + (H\sin \theta-b_0)^2 $$

Getting the minimal of $D^2$ we get

$$ \sin \theta (a_0-\cos \theta) - H\cos \theta (b_0-H\sin \theta) = 0 $$ $$ - (1-H^2) \cdot \sin \theta \cos \theta+ a_0 \cdot \sin \theta - Hb_0 \cdot \cos \theta = 0 $$

Then solve for $\theta$ with a numerical method.

Option 3.2: Use Lagrange multiplier

We have the minimizing function $D^2$ with the constraint of $(3)$. That means we have to find the solution for

$$ \nabla D^2 + \lambda \nabla G = \vec{0} $$ $$ \begin{bmatrix} 2(a-a_0) \\ 2(b-b_0) \end{bmatrix} + \lambda \begin{bmatrix} 2a \\ \dfrac{2b}{H^2} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} $$

Gives us

$$ a = \dfrac{u_0}{1+\lambda} \ \ \ \ \ \ \ \ \ \ \ \ b = \dfrac{v_0}{1+\frac{\lambda}{H^2}} $$

And in $(2)$ gives us

$$ u_0^2 \left(H+\dfrac{\lambda}{H}\right)^2 + v_0^2 \left(1+\lambda\right)^2 - \left(1+\lambda\right)^2 \left(H+\dfrac{\lambda}{H}\right)^2 = 0 $$

Then solve for $\lambda$ with a numerical method

Part 6: Scale back

Once we have $\theta$ or $\lambda$, we can compute the point $a$ and $b$

$$ a = \cos \theta \ \ \ \ \ \ \ \ \ \ \ \ \ \ b = H \sin \theta $$

or

$$ a = \dfrac{a_0}{1+\lambda} \ \ \ \ \ \ \ \ \ \ \ \ \ \ b = \dfrac{b_0}{1+\dfrac{\lambda}{H^2}}$$

To get the projection point we have to do

$$ \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} \cos \varphi & \sin \varphi \\ -\sin \varphi & \cos \varphi \end{bmatrix} \begin{bmatrix} A \cdot a \\ A \cdot b \end{bmatrix}+ \begin{bmatrix} x_c \\ y_c \end{bmatrix} $$