Find the equation of an ellipse using three points

conic sectionsgeometryrotationssoft-question

I came across this interesting problem yesterday and I am not quite able to find the equation of the ellipse after it has performed that roll. The original problem shows the ellipse to rotate till it is tangent to the x-axis at $(5,0)$ and that animation also shows the ellipse passing through two more points roughly speaking; $(0,7.9)$ and $(0.3.5)$.

To summarize, I require the equation of the rotated ellipse, knowing we have three points $(5,0)$, $(0,7.9)$ and $(0.3.5)$ and the fact that the ellipse is tangent to x-axis at $(5,0)$.

I have tried using $\frac{(x-x_0)^2}{a^2}+\frac{(y-y_0)^2}{b^2}=1$ as a standard equation (with 4 unknowns) ,but this just represents an ellipse with its axes parallel to coordinate axes so this was a dead end.

I am interested to understand how to write this ellipse's equation.

Best Answer

I will assume that $a$ and $b$ are known. What is unknown is the location of the center, and angle of rotation of the axes relative to the $xy$ axes.

The equation of this ellipse is therefore (in matrix-vector format) is

$ (r - C)^T Q (r - C) = 1$

where $ r = \begin{bmatrix} x \\ y \end{bmatrix} $ is a point on the ellipse, and $ C = \begin{bmatrix} C_x \\ C_y \end{bmatrix}$ is the center of the ellipse. Matrix $Q$ is $2 \times 2$, symmetric and positive definite. $Q$ can be written in factored form as

$ Q = R D R^T $

where $ D = \begin{bmatrix} \dfrac{1}{a^2} && 0 \\ 0 && \dfrac{1}{b^2} \end{bmatrix} $

and

$ R = \begin{bmatrix} \cos \theta && - \sin \theta \\ \sin \theta && \cos \theta \end{bmatrix} $

For a set of $ \theta, C_x , C_y $, the equation error of the ellipse equation at the point $r_i = \begin{bmatrix} x_i \\ y_i \end{bmatrix}$ is

$ e_i = (r_i - C)^T R D R^T (r_i - C) - 1 , i = 1, 2, 3$

Using the three points given (without the tangency condition), leads to the overall error function

$ E = \displaystyle \sum_{i=1}^3 e_i^2 = \sum_{i=1}^3 \left( (r_i - C)^T R D R^T (r_i - C) - 1 \right)^2 $

We want to minimize $E$, so taking the partial derivatives of $E$ with respect to the three variables $\theta, C_x, C_y$ gives

$\dfrac{\partial E}{\partial \theta} = 2 \sum_{i=1}^3 \left( (r_i - C)^T R D R^T (r_i - C) - 1 \right) \left( 2 (r_i - C)^T \dfrac{\partial R}{\partial \theta} D R^T (r_i - C) \right) = 0 $

$\dfrac{\partial E}{\partial C_x} = 2 \sum_{i=1}^3 \left( (r_i - C)^T R D R^T (r_i - C) - 1 \right) \left( - 2 e_1^T Q( r_i - C) \right) $

$\dfrac{\partial E}{\partial C_y} = 2 \sum_{i=1}^3 \left( (r_i - C)^T R D R^T (r_i - C) - 1 \right) \left( - 2 e_2^T Q( r_i - C) \right) $

where $e_1 = [1, 0]^T , e_2 = [0, 1]^T $. And

$ \dfrac{\partial R }{\partial \theta} = \begin{bmatrix} - \sin \theta && - \cos \theta \\ \cos \theta && - \sin \theta \end{bmatrix} $

These are three non-linear equations in $\theta, C_x, C_y$ of the form $\mathbf{f} = \mathbf{0}$ where $\mathbf{f}$ is a $3 \times 1$ vector of the right hand sides of the above three equations. One way that I know to solve this vector system for the unknowns is by the Newton-Raphson multivariate method. The method requires the computation of the Jacobian of the vector $\mathbf{f}$ which is defined as follows

$ J = [J_{ij}] $ where $J_{ij} = \dfrac{\partial \mathbf{f}_i }{ \partial \mathbf{x}_j } $

where $\mathbf{x}_1 = \theta , \mathbf{x}_2 = C_x , \mathbf{x}_3 = C_y $.

This computation is very lengthy if done exactly. Instead, it is possible to compute the Jacobian by numerical approximation (i.e. numerical differentiation). It is then trivially simple and easy to compute the Jacobian.

The Newton-Raphson iteration is as follows

$ \mathbf{x}_{k+1} = \mathbf{x}_k - J^{-1}_k \mathbf{f}_k $

The only catch for this method is the need for a "good" initial guess $\mathbf{x}_0$ of the solution. Sometimes, from the nature of the problem, it is possible to obtain such a good initial guess.

I've implemented the above procedure with the points $P_1 = (5, 0), P_2 = (0, 7.9), P_3 = (0, 3.8) $. The initial guess of the parameters that I used was $ \mathbf{x_0}_1 = \theta = - \dfrac{\pi}{6}$, $\mathbf{x_0}_2 = C_x = 4 $, $ \mathbf{x_0}_3 = C_y = 5 $. The Newton-Raphson iteration converged, and the resulting ellipse is shown below. It is quite close to the ellipse in your linked problem (the rolling ellipse problem).

enter image description here


In the second version of the identification procedure, we have only two points, and the angle that normal vector to the ellipse makes with the positive $x$ axis at the first point. In this case we have the equation error for the two points as before, defined as

$ e_i = (r_i - C)^T R D R^T (r_i - C) - 1 , i = 1, 2$

For the third error, we have the normal vector at $r_1$ given by

$ n_1 = R D R^T (r_1 - C) $

And we have the given direction of where the normal vector should be pointing as $n_0$ (assumed a unit vector). The error related to the direction, is defined as the squared difference between $\cos$ the angle between $n_0$ and $n_1$ and $1$. This can be computed using the dot product. So our error function to be minimized is

$ E = e_1^2 + e_2^2 + ( \cos(\phi (n_0, n_1)) - 1 )^2 $

where $\phi(n_0, n_1) = \cos^{-1} \left( \dfrac{n_0 \cdot n_1 }{\| n_1 \|} \right) $

This means that

$ E = e_1^2 + e_2^2 + \left( \dfrac{ n_0^T Q (r_1 - C) }{\| n_1 \|} - 1 \right)^2 $

We can simplify the last term, by redefining the error function to be

$ E = e_1^2 + e_2^2 + ( n_0^T Q (r_1 - C) - \sqrt{(r_1 - C)^T Q^2 (r_1 - C) })^2 $

The derivative of the first two terms of $E$ with respect to $\theta, C_x, C_y$ is the same as before. For the last term $T_3$, we have

$T_3 = K^2 $

with $ K = n_0^T Q (r_1 - C) - \sqrt{(r_1 - C)^T Q^2 (r_1 - C) } $

$ \dfrac{\partial T_3}{\partial \theta} = (2K) \bigg( 2 n_0^T \left(\dfrac{\partial R}{\partial \theta } D R^T + R D \dfrac{\partial R^T }{\partial \theta} \right) (r_1 - C) - \dfrac{ (r_1 - C)^T \dfrac{ \partial R}{\partial \theta} D^2 R^T ( r_1 - C) }{\sqrt{(r_1 - C)^T Q^2 (r_1 - C) }} \bigg)$

and

$\nabla_C T_3 = (2 K) \bigg( -2 Q n_0 + \dfrac{ Q^2 (r_1 - C) }{\sqrt{(r_1 - C)^T Q^2 (r_1 - C) }} \bigg) $

Now we have our non-linear equation which are the derivatives of $E$ with respect to $\theta, C_x, C_y$.

The rest is the same as before. Using the Newton-Raphson method, and with a good initial guess we arrive at the solution. I've implemented the above equations within the Newton-Raphson recursion, and the two points $(5, 0)$ and $(0, 7.8)$ with horizontal tangency at the first point, and I got convergence to the ellipse shown below which is almost identical to the one we obtained before (the one with three points and not tangency information).

enter image description here