The ellipse is known, but how to draw it?

Question

Let the equation of an ellipse in general be given by: $$ Ax^2+Bxy+Cy^2+Dx+Ey+F=0 $$ Everything real valued, and well known it seems. But I have one question: how can we efficiently draw such an ellipse?

My attempt

I have a Delphi Pascal Unit at my disposal that does the Art of Contouring. Simply calculate the isoline(s) of the function $\,f(x,y)=Ax^2+Bxy+Cy^2+Dx+Ey+F\,$ at level zero: $\,f(x,y)=0\,$. One such result is shown in the picture below. I've been able to reproduce in this way the outcome in an answer by Anders Kaseorg to the question why a Least square fit of ellipse worsens with increasing border thickness:

However, the Contouring method involves processing of all pixels in the canvas of a picture. Computers are fast nowadays and there is plenty of memory available, so it may sound old fashioned, but I find this method not very efficient. I wonder if there are better ways to do the drawing.

Answer 1

Starting with

$A x^2 + B xy + C y^2 + D x + E y + F = 0$

write the above in terms of the $2D$ vector $r = [x, y]^T$ as

$ r^T Q r + b^T r + c = 0 $

where,

$ Q = \begin{bmatrix} A && B/2 \\ B/2 && C \end{bmatrix} $

$ b = [ D , E ]^T $

$ c = F $

Since it is assumed that this quadratic equation represents an ellipse, then $Q$ is invertible. So the next step is to find the center of the ellipse, using the formula,

$ \text{Center} = r_0 = -\frac{1}{2} Q^{-1} b $

Now the equation, using the center $r_0$ becomes

$ (r - r_0)^T Q (r - r_0) - {r_0}^T Q {r_0} + c = 0 $

Dividing through by $( {r_0}^T Q {r_0} - c )$ , results in,

$ (r - r_0)^T E (r - r_0) = 1 $

where,

$E = \dfrac{1}{{r_0}^T Q {r_0} - c} Q = \begin{bmatrix} E_{11} && E_{12} \\ E_{12} && E_{22} \end{bmatrix} $

The next step is to diagonalize $E$ and this is always possible because $E$ is a symmetric matrix, and for that define $\theta$ as follows

$\theta = \frac{1}{2} \tan^{-1} \dfrac{2 E_{12} }{E_{11} - E_{22}} $

Then $E = R D R^T $

where $R$ is a rotation matrix by an angle $\theta$,

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

And $D$ is a diagonal matrix, whose diagonal entries are given by,

$D_{11} = E_{11} \cos^2 \theta + E_{22} \sin^2 \theta + 2 E_{12} \cos \theta \sin \theta $

$D_{22} = E_{11} \sin^2 \theta + E_{22} \cos^2 \theta - 2 E_{12} \cos \theta \sin \theta $

Now we have the diagonalization specification complete, so we can write,

$ (r - r_0)^T R D R^T (r - r_0) = 1 $

Define the vector $q =[q_x, q_y]^T = R^T (r - r_0) $, then

$ q^T D q = 1 $

And this reads $D_{11} q_x^2 + D_{22} q_y^2 = 1 $

which represents an ellipse in standard position in the $ q_x q_y $ coordinates, with semi-axis lengths of $\dfrac{1}{\sqrt{D_{11}}} $ and $\dfrac{1}{\sqrt{D_{22}}} $. Thus parametrically, we have

$(q_x, q_y)(t) = D^{-\frac{1}{2}} u(t) $

where $u(t) = (\cos t , \sin t ) $ is a rotating unit vector.

To obtain the ellipse in $r$ coordinates, we use

$ r = r_0 + R q = r_0 + R D^{-\frac{1}{2}} u $

The matrix $R D^{-\frac{1}{2}} = [v_1, v_2] $ is the matrix whose columns are the semi-axes vectors of the ellipse.

Answer 2

Any pair of conjugate diameters of an ellipse can be used as oblique cartesian axes, and finding them can be simpler than finding the axes of the ellipse.

First of all one must find the center of the ellipse:

$$ O=(x_O,y_O)= \left({2CD-BE\over B^2-4AC}, {2AE-BD\over B^2-4AC}\right). $$

Intersect now the ellipse with the line of equation $y=y_O$, passing through $O$ and parallel to $x$ axis: the intersection points $M$ and $N$ are the endpoints of the first diameter.

To find the other diameter just find the points of the ellipse where $dy/dx=0$, that is where $y=-(D+2Ax)/B$. These are the endpoints $P$ and $Q$ of the conjugate diameter.

Finally, a parametric equation of the ellipse can be written as:

$$ (x,y)=(M-O)\cos t+(P-O)\sin t. $$

Answer 3

On one hand, we have the general equation for an ellipse: $$ Ax^2+Bxy+Cy^2+Dx+Ey+F=0 $$ On the other hand, the following expression for an ellipse of inertia has been proposed elsewhere: $$ \frac{\sigma_{yy} (x-\mu_x)^2 - 2 \sigma_{xy} (x-\mu_x) (y-\mu_y) + \sigma_{xx} (y-\mu_y)^2}{ \sigma_{xx} \sigma_{yy} - \sigma_{xy}^2 } = 1 $$ The precise meaning of the coefficients is not important at this stage, as shown below.
Quoting from the accepted answer: First of all one must find the center of the ellipse.
For reasons of consistency (with the work called elsewhere) we put: $(x_O,y_O)=(\mu_x,\mu_x)$. Consider: $$ A(x-\mu_x)^2+B(x-\mu_x)(y-\mu_y)+C(y-\mu_y)^2 = \\ Ax^2-2A\mu_x x+A\mu_x^2 + Bxy-B\mu_x y-B\mu_y x+B\mu_x\mu_y + Cy^2-2C\mu_y y+C\mu_y^2 = \\ Ax^2+Bxy+Cy^2 + (-2A\mu_x-B\mu_y)x + (-2C\mu_y-B\mu_x)y + (A\mu_x^2+B\mu_x\mu_y+C\mu_y^2) $$ This result must be compatible with: $Ax^2+Bxy+Cy^2+Dx+Ey+F$, so: $$ \begin{cases} 2A\mu_x+B\mu_y = -D \\ B\mu_x+2C\mu_y = -E \end{cases} \quad \Longrightarrow \quad \begin{bmatrix} 2C & -B \\ -B & 2A \end{bmatrix} / (4AC-B^2) \begin{bmatrix} \mu_x \\ \mu_y \end{bmatrix} = \begin{bmatrix} -D \\ -E \end{bmatrix} $$ The solution is in concordance with the the answer by Intelligenti pauca: $$ \begin{bmatrix} \mu_x \\ \mu_y \end{bmatrix} = \begin{bmatrix} 2CD-BE \\ 2AE-BD \end{bmatrix} / (B^2-4AC) $$ Furthermore we have: $$ Ax^2+Bxy+Cy^2+Dx+Ey+F = 0 \\ A(x-\mu_x)^2+B(x-\mu_x)(y-\mu_y)+C(y-\mu_y)^2 - (A\mu_x^2+B\mu_x\mu_y+C\mu_y^2) + F = 0 \\ \frac{A(x-\mu_x)^2 - 2(-B/2)(x-\mu_x)(y-\mu_y) + C(y-\mu_y)^2}{AC-(-B/2)^2} = \frac{A\mu_x^2+B\mu_x\mu_y+C\mu_y^2-F}{AC-(-B/2)^2} $$ To be comparable with: $$ \frac{\sigma_{yy} (x-\mu_x)^2 - 2 \sigma_{xy} (x-\mu_x) (y-\mu_y) + \sigma_{xx} (y-\mu_y)^2}{ \sigma_{xx} \sigma_{yy} - \sigma_{xy}^2 } = 1 $$ We are almost there. Define: $$ \lambda = \frac{A\mu_x^2+B\mu_x\mu_y+C\mu_y^2-F}{AC-(-B/2)^2} $$ Then we finally have: $$ \frac{(\lambda A)(x-\mu_x)^2 - 2(-\lambda B/2)(x-\mu_x)(y-\mu_y) + (\lambda C)(y-\mu_y)^2}{(\lambda A)(\lambda C)-(-\lambda B/2)^2} = 1 $$ Cast into the ellipse of inertia format: $$ \sigma_{yy} = \lambda A \quad ; \quad \sigma_{xy} = -\lambda B/2 \quad ; \quad \sigma_{xx} = \lambda C $$ Now we can apply the method as described in Drawing an ellipse: $$ \begin{cases} x = \mu_x + a_x\cos(t) + a_y\sin(t) \\ y = \mu_y + b_x\cos(t) + b_y\sin(t) \end{cases} $$ With $(\mu_x,\mu_y)$ as defined above and: $$ (a_x,a_y) = (\sqrt{\sigma_{xx}},0) \\ (b_x,b_y) = \left(\frac{\sigma_{xy}}{\sqrt{\sigma_{xx}}},\sqrt{\frac{\sigma_{xx}\sigma_{yy}-\sigma_{xy}^2}{\sigma_{xx}}}\right) $$ Needless to say that the picture in the question can be exactly reproduced in this way.