Parametrization of an arbitrary conic/ellipse

Question

I have the coefficients for a conic (I actually know that it is an ellipse) in the form $$Ax^2 + Bxy + Cy^2 + Dx + Ey + F =0$$

Is there an efficient algorithm which returns the parametrization of the eclipse, i.e., $\langle x(t), y(t) \rangle?$

Answer 1

If all tools are at my disposal

Find the center. Take the partial derivatives.

$2Ax + By + D = 0\\ Bx + 2Cy + E = 0$

Solve this system of linear differential equations.

Lets call this $(x_0, y_0)$

$\begin{bmatrix} 2A & B\\B&2C \end{bmatrix}\begin {bmatrix} x_0\\y_0 \end{bmatrix} = \begin{bmatrix} -D\\-E\end{bmatrix}\\ \begin {bmatrix} x_0\\y_0 \end{bmatrix} = \frac {1}{B^2 - 4AC}\begin{bmatrix} 2CD-BE\\2AE - BD\end{bmatrix}$

$A(X-x_0)^2 + B(x-x_0)(y-y_0) + C(y-y_0)^2 = -F - Ax_0^2-Bx_0y_0 - C y_0^2$

Let $K^2 = -F - Ax_0^2-Bx_0y_0 - C y_0^2$

$\begin{bmatrix} (x-x_0)&(y-y_0)\end{bmatrix}\begin{bmatrix} A &\frac 12B\\\frac 12 B & C\end{bmatrix}\begin{bmatrix} (x-x_0)\\(y-y_0)\end{bmatrix} = K^2$

Diagonalize that matrix. The eigenvectors will be a rotation matrix. The eigenvalues will give you the coefficients. If it is an ellipse, the eigenvalues will both be positive.

Since you are trying to code this...

$\lambda^2 - (A+C) \lambda + AC-B^2 = 0\\ \lambda = \frac {(A+C)}{2} \pm \frac {\sqrt {(A-C)^2+4B^2}}{2}$

$\phi = \arctan {B}{A-\lambda_1} = \arctan \frac {C-\lambda_2}{B}$

(I hope I have that right)

$x = \frac {K}{\lambda_1} \cos (\theta-\phi) + x_0\\ y = \frac {K}{\lambda_2} \sin (\theta-\phi) + y_0$

Answer 2

Based on the comments to your question, it looks like the underlying problem that you’re trying to solve is to determine efficiently whether or not a point lies within the ellipse. A straightforward way to do this is to plug the coordinates of the point into the left-hand side of the general equation. The sign of the result will discriminate among the three possibilities. If you arrange for the leading coefficient $A$ to be positive (as will $C$ in that case), then the point is outside of the ellipse if the result is positive, inside if negative and of course on the ellipse if zero. Points very near the boundary might run into machine resolution limitations, so beware of that.

The reason this works is that, if you start from the unit circle $x^2+y^2=1$, for which the appropriate inequalities are obvious, any ellipse can be obtained from it via a combination of scaling, rotation and translation. This transformation preserves the interior and exterior, so the inequalities involving the transformed equation have the same senses. With positive scale factors, the coefficients $A$ and $C$ in the general equation will also be positive.

If you have a lot of points to test against a fixed ellipse and most of them are outside of the ellipse’s bounding box, it could be more efficient overall to do some preliminary range checks before evaluating the general conic expression. Finding the points on the ellipse where the partial derivative w/r $y$ vanishes will give you the $x$-coordinate bounds, while the points where the $x$-derivative vanishes will give you the $y$ bounds.

Another approach to finding the bounding box is to find the horizontal and vertical tangents to the ellipse. Working in homogeneous coordinates, the tangent lines to an ellipse given by the matrix $Q$ that go through a point $\mathbf p$ are captured by the degenerate conic $T=\mathcal M_{\mathbf p}^TQ^{-1}\mathcal M_{\mathbf p}$. Here, $\mathcal M_{\mathbf p}$ is the “cross-product matrix” of $\mathbf p$, i.e., $\mathcal M_{\mathbf p}\mathbf x=\mathbf p\times\mathbf x$. Instead of $Q^{-1}$, we can use the adjugate of $Q$ or any other convenient non-zero multiple of it instead.

Taking $\mathbf p = [0:1:0]$ will produce the vertical tangents to $Q$. For the general conic equation, then, we have $$Q = \begin{bmatrix}A&\frac B2&\frac D2\\\frac B2&C&\frac E2\\\frac D2&\frac E2&F \end{bmatrix}$$ and so $$T=\begin{bmatrix}0&0&-1\\0&0&0\\1&0&0\end{bmatrix} \begin{bmatrix} 4CF-E^2 & DE-4BF & 2BE-2CD \\ DE-4BF & 4AF-D^2 & 2BD-2AE \\ 2BE-2CD & 2BD-2AE & 4AC-4B^2 \end{bmatrix} \begin{bmatrix}0&0&1\\0&0&0\\-1&0&0\end{bmatrix} = \begin{bmatrix} 4AC-4B^2 & 0 & 2CD-2BE \\ 0&0&0 \\ 2CD-2BE & 0 & 4CF-E^2 \end{bmatrix}.$$ Now, you could use this directly for the range check by checking the sign of $(4AC-4B^2)\,x^2+(4CD-4BE)\,x+(4CF-E^2)$, but that’s not a whole lot better than simply using the original conic equation. However, we can tease out the individual lines by “splitting” this degenerate conic: Find an $\alpha$ for which $T+\alpha\mathcal M_{\mathbf p}$ has rank one. The two lines are then any row and column of the resulting matrix that have a common nonzero diagonal element. All of the $2\times2$ minors of $T+\alpha\mathcal M_{\mathbf p}$ are identically zero except for $$\alpha^2-4C^2(D^2-4AF)-4C(AE^2+4B^2F-2BDE).$$ Taking either root results in the equations $$2(AC-B^2)\,x = (BE-CD)\pm\sqrt{C^2(D^2-4AF)+C(AE^2+4B^2F-2BDE)}$$ for the tangents, which give you the left and right edges of the bounding box. The calculation for the top and bottom is similar, but with $\mathbf p=[1:0:0]$ instead, and again there’ll only be one $2\times2$ minor of $T+\alpha\mathcal M_{\mathbf p}$ that doesn’t vanish identically. The resulting equations for the horizontal tangents are $$2(AC-B^2)\,y=(BD-AE)\pm\sqrt{A^2(E^2-4CF)+A(CD^2+4B^2F-2BDE)}$$ from which we can read the max. and min. $y$-coordinates of the ellipse.