How do I find the minimum-area ellipse that encloses a set of points?

Question

I have a set of points that resembles more of an ellipse than a circle. I implemented the optimization formulation below and the solution gives a circle. I tried with various initial values, still to no avail. Have I overlooked or mistaken something?

$$ \begin{align} \min_{r_{maj},r_{min}} \quad & \pi r_{maj}r_{min} \\ s.t. \quad & \left(\frac{x_i-x_c}{r_{maj}}\right)^2 + \left(\frac{y_i-y_c}{r_{min}}\right)^2 \leq 1 \quad \forall\;i \end{align} $$

$r_{maj},r_{min}$: major and minor radius
$(x_c,y_c)$: ellipse centre

the red cross is the 'ellipse' centre obtained from the optimization result.

Answer 1

Theory

The 1997 paper "Smallest Enclosing Ellipses -- Fast and Exact" by Gärtner and Schönherr addresses this question. The same authors provide a C++ implementation in their 1998 paper "Smallest Enclosing Ellipses -- An Exact and Generic Implementation in C++". The paper's 150 pages long, but, fortunately for us, the venerable CGAL library implements the algorithm as described here.

However, CGAL only provides the general equation for the ellipse, so some massaging is necessary to convert the general equation to the canonical equation and, from there, to get a parametric equation suitable for plotting.

All this is included in the implementation below.

Using WebPlotDigitizer to extract your data gives:

-2.024226110363392 5.01315585320002
1.9892328398384915 3.0400196291391692
-0.0269179004037694 1.980973446537659
-0.987886944818305 -0.9505049812548929
4.9851951547779265 -1.9398893695943187

Fitting this using the program below gives:

a = 2.47438
b = 5.42919
cx = 0.767976
cy = 0.792924
theta = 0.784877

We can then plot this with gnuplot

set parametric
plot "points" pt 7 ps 2, [0:2*pi] a*cos(t)*cos(theta) - b*sin(t)*sin(theta) + cx, a*cos(t)*sin(theta) + b*sin(t)*cos(theta) +
cy lw 2

to get

Implementation

The code below does this:

// Compile with clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O2 -g -DNDEBUG -Wall -Wextra -pedantic -march=native -frounding-math main.cpp -lgmpxx -lmpfr -lgmp

#include <CGAL/Cartesian.h>
#include <CGAL/Min_ellipse_2.h>
#include <CGAL/Min_ellipse_2_traits_2.h>
#include <CGAL/Exact_rational.h>

#include <cassert>
#include <cmath>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

typedef CGAL::Exact_rational             NT;
typedef CGAL::Cartesian<NT>              K;
typedef CGAL::Point_2<K>                 Point;
typedef CGAL::Min_ellipse_2_traits_2<K>  Traits;
typedef CGAL::Min_ellipse_2<Traits>      Min_ellipse;


struct EllipseCanonicalEquation {
  double semimajor; // Length of semi-major axis
  double semiminor; // Length of semi-minor axis
  double cx;        // x-coordinate of center
  double cy;        // y-coordinate of center
  double theta;     // Rotation angle
};

std::vector<Point> read_points_from_file(const std::string &filename){
  std::vector<Point> ret;

  std::ifstream fin(filename);
  float x,y;
  while(fin>>x>>y){
    std::cout<<x<<" "<<y<<std::endl;
    ret.emplace_back(x, y);
  }

  return ret;
}


// Uses "Smallest Enclosing Ellipses -- An Exact and Generic Implementation in C++"
// under the hood.
EllipseCanonicalEquation get_min_area_ellipse_from_points(const std::vector<Point> &pts){
  // Compute minimum ellipse using randomization for speed
  Min_ellipse  me2(pts.data(), pts.data()+pts.size(), true);
  std::cout << "done." << std::endl;

  // If it's degenerate, the ellipse is a line or a point
  assert(!me2.is_degenerate());

  // Get coefficients for the equation
  // r*x^2 + s*y^2 + t*x*y + u*x + v*y + w = 0
  double r, s, t, u, v, w;
  me2.ellipse().double_coefficients(r, s, t, u, v, w);

  // Convert from CGAL's coefficients to Wikipedia's coefficients
  // A*x^2 + B*x*y + C*y^2 + D*x + E*y + F = 0
  const double A = r;
  const double B = t;
  const double C = s;
  const double D = u;
  const double E = v;
  const double F = w;

  // Get the canonical form parameters
  // Using equations from https://en.wikipedia.org/wiki/Ellipse#General_ellipse
  const auto a = -std::sqrt(2*(A*E*E+C*D*D-B*D*E+(B*B-4*A*C)*F)*((A+C)+std::sqrt((A-C)*(A-C)+B*B)))/(B*B-4*A*C);
  const auto b = -std::sqrt(2*(A*E*E+C*D*D-B*D*E+(B*B-4*A*C)*F)*((A+C)-std::sqrt((A-C)*(A-C)+B*B)))/(B*B-4*A*C);
  const auto cx = (2*C*D-B*E)/(B*B-4*A*C);
  const auto cy = (2*A*E-B*D)/(B*B-4*A*C);
  double theta;
  if(B!=0){
    theta = std::atan(1/B*(C-A-std::sqrt((A-C)*(A-C)+B*B)));
  } else if(A<C){
    theta = 0;
  } else { //A>C
    theta = M_PI;
  }

  return EllipseCanonicalEquation{a, b, cx, cy, theta};
}


int main(int argc, char** argv){
  if(argc!=2){
    std::cerr<<"Provide name of input containing a list of x,y points"<<std::endl;
    std::cerr<<"Syntax: "<<argv[0]<<" <Filename>"<<std::endl;
    return -1;
  }

  const auto pts = read_points_from_file(argv[1]);

  const auto eq = get_min_area_ellipse_from_points(pts);

  // Convert canonical equation for rotated ellipse to parametric based on:
  // https://math.stackexchange.com/a/2647450/14493
  std::cout << "Ellipse has the parametric equation " << std::endl;
  std::cout << "x(t) = a*cos(t)*cos(theta) - b*sin(t)*sin(theta) + cx"<<std::endl;
  std::cout << "y(t) = a*cos(t)*sin(theta) + b*sin(t)*cos(theta) + cy"<<std::endl;
  std::cout << "with" << std::endl;

  std::cout << "a = " << eq.semimajor << std::endl;
  std::cout << "b = " << eq.semiminor << std::endl;
  std::cout << "cx = " << eq.cx << std::endl;
  std::cout << "cy = " << eq.cy << std::endl;
  std::cout << "theta = " << eq.theta << std::endl;

  return 0;
}

Answer 2

With your formulation of the constraints, you can only produce ellipses where the semi-major and semi-minor axes are aligned with the coordinate axes, but it's clear from the figure that you attached that you can shrink the size of the ellipse slightly along the vector $(1, 1)^\top / \sqrt{2}$ (or thereabouts). Another formulation that you could try is to look for a center point $\mathbf{p}_c$ and a 2 x 2, symmetric, positive-definite matrix $\mathbf{A}$ such that

$$(\mathbf{p}_i - \mathbf{p}_c)^\top\cdot\mathbf{A}\cdot(\mathbf{p}_i - \mathbf{p}_c) \le 1$$

for all $i$. In order to optimize for the ellipse of smallest area, you'll want to find some formula for the area of that ellipse in terms of $\mathbf{A}$. Symmetric 2 x 2 matrices have 3 degrees of freedom, and it's that extra degree that includes a rotation.

Some hints:

1. It's 2 x 2 so everything can be done in terms of the trace and determinant.
2. What are the trace and determinant in terms of the eigenvalues of the matrix and how do those relate to the ellipse?
3. The function $f(\mathbf{A}) = -\log\det\mathbf{A}$ is convex.