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.

  • 1
    $\begingroup$ I dont think that the problem, at least as stated, is convex. Are you solving for $r_{\min},r_{\max},x_c$ and $y_c$? Or just $r_{\min}$ and $r_{\max}$? Can you provide the data for those 5 poins? $\endgroup$ Commented Nov 21, 2021 at 8:03
  • 3
    $\begingroup$ Another approach to solve this problem is to use principal component analysis to get an approximation of the elliptical shape of the points. Then you can scale this ellipse until it contains all of your points. $\endgroup$ Commented Nov 21, 2021 at 11:17
  • $\begingroup$ What do you think the correct answer should look like? The circle looks correct to me. Perhaps you want to rotate the ellipse? $\endgroup$ Commented Nov 21, 2021 at 12:37
  • 1
    $\begingroup$ This is a well studied problem. See Mike Todd's book, Minimum-Volume Ellipsoids: Theory and Algorithms. $\endgroup$ Commented Nov 22, 2021 at 0:21

2 Answers 2



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


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

  // 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;
    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){
    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;
  • $\begingroup$ Thank you Richard! I was trying to reproduce some results using the formulation I got from a paper with a clogged brain without going deeper into the basics. $\endgroup$ Commented Nov 24, 2021 at 0:16

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.
  • 2
    $\begingroup$ This is a very good answer, but the fact that logdet is concave still means that at the end of the day you're minimizing a concave function. $\endgroup$ Commented Nov 21, 2021 at 21:24
  • $\begingroup$ Thank you Daniel for the introduction, I understood this further from Mike Todd's book too. $\endgroup$ Commented Nov 24, 2021 at 0:18
  • $\begingroup$ @RobertBasset I think it's okay. The bigger $\log \det A$ is, the smaller is the ellipse volume. Therefore we would like to maximize $\log \det A$, or equivalently, minimize $-\log \det A$. Here is a reference I found people.duke.edu/~psun/bio/MVEpaper_accept.pdf $\endgroup$
    – Mark
    Commented Aug 4, 2023 at 17:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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