Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I have a set $X$ of points in $\mathbb{R}^2$ and I'm trying to find the smallest encompassing ellipse which main axes are parallel to the coordinate system's (to put it differently, its both centres share one coordinate). I need the gravitational centre and the vertical and horizontal radius.

Now, I managed to do this for a horizontally aligned rectangle, but that isn't much help (although it's a first approximation, as I can easily draw a rectangle around $X$). I also found some formulas on the internet, but they seem to be wrong as I keep getting $0$ for the radius.

Can this be done algebraically?

share|improve this question
    
Minimal with respect to what? Minimal area, perimeter, $\max\{a,b\}$ or what? –  enzotib Oct 5 '12 at 18:38
    
@enzotib: Minimal area. –  bitmask Oct 5 '12 at 20:07

2 Answers 2

I have an algorithm, but no closed formula for you.

Algorithm

Compute the convex hull of your points. The ellipsis you search for has to touch the corners of that hull in several points if it is to be the smallest one. Now iterate over all 4-element subsets of those corners, and solve the following equation:

$$\begin{pmatrix} x_1^2 & y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 & y_4 & 1 \end{pmatrix}\cdot\begin{pmatrix} a\\b\\c\\d\\e \end{pmatrix}=0 $$

The result won't be unique, but instead a one-dimensional space of solutions unless you hit a degenerate situation. Take any non-zero vector from that space to describe your conic secion. You may e.g. add a fifth equation $a=1$ as long as you ensure that you will recover from any degeneracies (i.e. unsolvable systems due to linear dependency).

The conic section described by the parameters you found will already be “parallel” to the coordinate axes, and it can be written as

$$ ax^2 + by^2 + cx + dy + e = 0 $$

Check whether this conic section is actually an ellipsis, and if it contains all other points of the convex hull. If it is and does, compute its area, and compare that against the best solution found so far.

Details

To perform all the computations you need, let's turn the above formula into something easier to manage.

\begin{align*} 0 &= ax^2 + by^2 + cx + dy + e \\ &= a\left(x+\frac{c}{2a}\right)^2 + b\left(y+\frac{d}{2b}\right)^2 - \left(\frac{c^2}{4a} + \frac{d^2}{4b} - e\right) \\ m_x &= -\frac{c}{2a} \\ m_y &= -\frac{d}{2b} \\ s &= \left(\frac{c^2}{4a} + \frac{d^2}{4b} - e\right) \end{align*}

$(m_x, m_y)$ describes the center of symmetry for your conic.

If it is to be an ellipsis, then $a$ and $b$ have to have the same non-zero sign (i.e. $ab>0$), in which case $s$ will have the same sign as well. Now divide your whole equation by $s$ to obtain

\begin{align*} \frac as(x-m_x)^2 + \frac bs(y-m_y)^2 + \frac ss &= 0 \\ \frac as(x-m_x)^2 + \frac bs(y-m_y)^2 &= 1 \\ a'(x-m_x)^2 + b'(y-m_y)^2 &= 1 \end{align*}

with

$$a' = \frac as \quad b' = \frac bs$$

denoting the inverse of the square lengths of the radii of your ellipsis.

You can turn the equality into an inequality $\ldots < 1$ to describe the interior of the ellipsis, and check whether the other points of the convex hull fall inside. You can use the term $a'\cdot b'$ to compare sizes, as it is indirectly proportional to the squared area of your ellipsis.

If you need the actual radii, you can compute them as

\begin{align*} r_x &= \frac1{\sqrt{a'}} = \sqrt{\frac sa} \\ r_y &= \frac1{\sqrt{b'}} = \sqrt{\frac sb} \\ 1 &= \left(\frac{x-m_x}{r_x}\right)^2 + \left(\frac{y-m_y}{r_y}\right)^2 \end{align*}

Alternatives

For users who read this question while investigating a different task:

  • For arbitrary center and fixed orientation, we used four points to define the conic, as seen above.
  • For arbitrary center and arbitrary orientation, we need five points, as five points uniquely define a general conic. Add a $xy$ monomial to the equation.
  • For a fixed center like the origin, and fixed orientation, use two points, as $c=d=0$.
  • For fixed center and arbitrary orientation, use three points.
share|improve this answer
    
Thanks. An algorithm will do fine. I'll give it a shot. –  bitmask Oct 5 '12 at 16:58

There is a very elegant algorithm to find an encompassing ellipse for points arbitrarily positioned in space via PCA (Principal Component Analysis) approach to find axes of the ellipse.

Suppose you have $n$ points of an ellipse stored in the matrix $E$:

$$ E = \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \\ \vdots & \vdots \\ x_n & y_n \end{bmatrix} $$

Let us also suppose that the center of the ellipse is in the origin $(0, 0)$ (if not, just subtract corresponding means from each of the coordinates, i.e. $x_i := x_i-\mu_x, y_i := y_i - \mu_y $).

Then, compute the covariance matrix $C$, in our case we can use the following formula as the resulting matrix will be proportional to the real covariance matrix:

$$ C = E^TE $$

Finally, find eigenvectors $v_i$ of $C$, in our particular case, $i = 1, 2$. The resulting eigenvectors will have the same direction as the axes of the ellipse you work on. Also, if you normalize a matrix $V = [v_1\;v_2]$, you may see it is actually a rotation matrix. So, if you want to rotate your ellipse in a way that all its points are parallel to coordinate axes, just do the multiplication: $EV$.

In order to find magnitudes of ellipse axes, use the rotated ellipse $EV$: find the biggest values of $x$ and $y$ in the set of points, and these values are the ones you need.

This algorithm will definitely work if the points lie in a 3D space, and it does not require points to lie exactly on the ellipse.

Sources on PCA, eigenvalues, eigenvectors:

  1. http://www.doc.ic.ac.uk/~dfg/ProbabilisticInference/IDAPILecture15.pdf
  2. http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.