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

An answer to this question on the Scientific Computing Stack Exchange.

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

enter image description here
the red cross is the 'ellipse' centre obtained from the optimization result.

Answer

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

Smallest area ellipse of a set of points

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;
}