Find projection point on elipse 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.
 A: 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:

*

*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$


*Scale into a new ellipse of major axis equal 1
$$G(a, b) = a^2 + \left(\dfrac{b}{H}\right)^2 - 1  \ \ \ \ \ \ \ \ \ \ \ \ \ (3)$$


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


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


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


*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}
$$
A: 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$.
