11
$\begingroup$

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.

$\endgroup$
4
  • 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$ Nov 21 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$ Nov 21 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$ Nov 21 at 12:37
  • 1
    $\begingroup$ This is a well studied problem. See Mike Todd's book, Minimum-Volume Ellipsoids: Theory and Algorithms. $\endgroup$ Nov 22 at 0:21
21
$\begingroup$

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;
}
$\endgroup$
1
  • $\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$ Nov 24 at 0:16
10
$\begingroup$

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.
$\endgroup$
2
  • 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$ Nov 21 at 21:24
  • $\begingroup$ Thank you Daniel for the introduction, I understood this further from Mike Todd's book too. $\endgroup$ Nov 24 at 0:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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