理论

Gärtner 和 Schönherr 于1997 年发表的论文“最小封闭椭圆——快速且精确”解决了这个问题。相同的作者在他们 1998 年的论文“ Smallest Enclosure Ellipses -- An Exact and Generic Implementation in C++ ”中提供了 C++ 实现。这篇论文有 150 页长，但幸运的是，古老的 CGAL 库实现了这里描述的算法。

但是，CGAL 只提供椭圆的一般方程，因此需要进行一些处理才能将一般方程转换为规范方程，并从那里得到适合绘图的参数方程。

所有这些都包含在下面的实现中。

使用WebPlotDigitizer提取您的数据给出：

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

使用以下程序进行拟合给出：

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

然后我们可以用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

要得到

执行

下面的代码执行此操作：

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

通过制定约束，您只能生成半长轴和半短轴与坐标轴对齐的椭圆，但从您所附的图中可以清楚地看出，您可以沿着向量（或附近）。您可以尝试的另一个公式是寻找一个中心点和一个 2 x 2、对称、正定矩阵使得$(1, 1)^\top / \sqrt{2}$$\mathbf{p}_c$$\mathbf{A}$

对于所有。为了优化最小面积的椭圆，您需要根据找到该椭圆面积的一些公式。对称的 2 x 2 矩阵有 3 个自由度，而这个额外的自由度包括旋转。$i$$\mathbf{A}$

一些提示：

1. 它是 2 x 2，所以一切都可以根据迹线和行列式来完成。
2. 就矩阵的特征值而言，迹线和行列式是什么？它们与椭圆有何关系？
3. 函数是凸的。$f(\mathbf{A}) = -\log\det\mathbf{A}$