在c ++中查找包含一组点的最小面积椭圆
Posted
技术标签:
【中文标题】在c ++中查找包含一组点的最小面积椭圆【英文标题】：Find Minimum area ellipse enclosing a set of points in c++ 【发布时间】：2017-03-20 05:17:04 【问题描述】：我有一组二维点。我需要找到一个包围所有点的最小面积椭圆。有人可以给出如何解决这个问题的想法。对于一个圆圈来说，这很简单。中心与点之间的最大距离。但是对于椭圆来说，它非常复杂，我不知道。我必须在 C++ 中实现它。
【问题讨论】：
如果这个问题没有封闭形式的解决方案，它似乎很适合某种启发式搜索技术。 椭圆必须以原点为中心吗？椭圆的轴必须平行于坐标轴吗？ （这里的任何否定答案都会使问题变得非常复杂。） 我已经重新标记了您的问题（当您明确表示需要 C++ 时，为什么要标记 JAVA？） JAVA 在不知不觉中被标记， 椭圆的中心是否一定在 (0,0) 并且轴不旋转？如果没有，在一般情况下，您拥有给出正确解决方案的 MVEE（最小体积封闭椭圆）。 【参考方案1】：不确定我是否可以证明这一点，但在我看来，最佳解决方案的特征在于（至少）3 个点相切，而所有其他点都在椭圆内（想想看！）。因此，如果不出意外，您应该能够通过检查所有 ~n^3 个三元组点并检查它们是否定义了解决方案来强制它。应该可以通过删除必须严格位于任何周围椭圆内的所有点来改进这一点，但我不确定如何做到这一点。也许通过 x 和 y 坐标对点进行排序，然后做一些花哨的事情。
不是一个完整的解决方案，但它是一个开始。
编辑： 不幸的是，3 个点不足以定义一个椭圆。但也许如果你把它限制在与 3 点相切的最小区域的椭圆上？
【讨论】：
【参考方案2】：正如 Rory Daulton 建议的那样，您需要清楚地指定解决方案的约束条件，并且消除任何会使事情变得非常复杂。对于初学者来说，现在假设：
这是2D问题 椭圆是轴对齐的 中心是任意的而不是(0,0)
我会用approximation search（二分搜索和线性搜索之间的混合）将其作为标准genere and test问题来攻击它以加快速度（但你也可以从一开始就尝试蛮力，所以你看看它是否有效）。
计算解的约束
要限制搜索，您需要找到椭圆的大致放置位置和大小。为此，您可以使用外切圆作为您的点。很明显，椭圆面积将小于或等于圆，并且放置将在附近。圆圈不一定是最小的，因此我们可以使用以下示例：
-
找到点的边界框
让圆以该边界框为中心，半径为从其中心到任何点的最大距离。
这将是O(n) 复杂度，其中n 是您的积分数。
搜索“所有”可能的省略号并记住最佳解决方案
所以我们需要找到椭圆中心(x0,y0) 和半轴rx,ry 而area = M_PI*rx*ry 是最小的。对于近似搜索，每个变量的因子为O(log(m))，每次迭代都需要测试有效性，即O(n)，因此最终复杂度为O(n.log^4(m))，其中m是每个搜索参数可能变化的平均数量（取决于准确性和搜索约束）。通过简单的粗略搜索，它将是 O(n.m^4)，这真的很可怕，尤其是对于 m 可能非常大的浮点数。
为了加快速度，我们知道椭圆的面积将小于或等于找到的圆的面积，因此我们可以忽略所有较大的椭圆。 rx,ry 的约束可以从边界框的纵横比 +/- 一些保留值中得出。
这里是使用上面链接中的 approx 类的简单 C++ 小示例：
//---------------------------------------------------------------------------
// input points
const int n=15; // number of random points to test
float pnt[n][2];
// debug bounding box
float box_x0,box_y0,box_x1,box_y1;
// debug outscribed circle
float circle_x,circle_y,circle_r;
// solution ellipse
float ellipse_x,ellipse_y,ellipse_rx,ellipse_ry;
//---------------------------------------------------------------------------
void compute(float x0,float y0,float x1,float y1) // cal with bounding box where you want your points will be generated
int i;
float x,y;
// generate n random 2D points inside defined area
Randomize();
for (i=0;i<n;i++)
pnt[i][0]=x0+(x1-x0)*Random();
pnt[i][1]=y0+(y1-y0)*Random();
// compute bounding box
x0=pnt[0][0]; x1=x0;
y0=pnt[0][1]; y1=y0;
for (i=0;i<n;i++)
x=pnt[i][0]; if (x0>x) x0=x; if (x1<x) x1=x;
y=pnt[i][1]; if (y0>y) y0=y; if (y1<y) y1=y;
box_x0=x0; box_x1=x1;
box_y0=y0; box_y1=y1;
// "outscribed" circle
circle_x=0.5*(x0+x1);
circle_y=0.5*(y0+y1);
circle_r=0.0;
for (i=0;i<n;i++)
x=pnt[i][0]-circle_x; x*=x;
y=pnt[i][1]-circle_y; y*=y; x+=y;
if (circle_r<x) circle_r=x;
circle_r=sqrt(circle_r);
// smallest area ellipse
int N;
double m,e,step,area;
approx ax,ay,aa,ab;
N=3; // number of recursions each one improves accuracy with factor 10
area=circle_r*circle_r; // solution will not be bigger that this
step=((x1-x0)+(y1-y0))*0.05; // initial position/size step for the search as 1/10 of avg bounding box size
for (ax.init( x0, x1,step,N,&e);!ax.done;ax.step()) // search x0
for (ay.init( y0, y1,step,N,&e);!ay.done;ay.step()) // search y0
for (aa.init(0.5*(x1-x0),2.0*circle_r,step,N,&e);!aa.done;aa.step()) // search rx
for (ab.init(0.5*(y1-y0),2.0*circle_r,step,N,&e);!ab.done;ab.step()) // search ry
e=aa.a*ab.a;
// is ellipse outscribed?
if (aa.a>=ab.a)
m=aa.a/ab.a; // convert to circle of radius rx
for (i=0;i<n;i++)
x=(pnt[i][0]-ax.a); x*=x;
y=(pnt[i][1]-ay.a)*m; y*=y;
// throw away this ellipse if not
if (x+y>aa.a*aa.a) e=2.0*area; break;
else
m=ab.a/aa.a; // convert to circle of radius ry
for (i=0;i<n;i++)
x=(pnt[i][0]-ax.a)*m; x*=x;
y=(pnt[i][1]-ay.a); y*=y;
// throw away this ellipse if not
if (x+y>ab.a*ab.a) e=2.0*area; break;
ellipse_x =ax.aa;
ellipse_y =ay.aa;
ellipse_rx=aa.aa;
ellipse_ry=ab.aa;
//---------------------------------------------------------------------------
即使是这个只有 15 个点的简单示例也需要大约 2 秒的时间来计算。您可以通过添加启发式方法来提高性能，例如仅低于 circle_r^2 的测试区域等，或者使用一些数学规则更好地选择解决方案区域。如果您使用蛮力而不是近似搜索，预计计算时间甚至可能是几分钟或更长时间，因此 O(scary)...
请注意，此示例不适用于点的任何纵横比，因为我将 rx,ry 的上限硬编码为 2.0*circle_r，这可能还不够。相反，您可以根据点的纵横比和/或条件来计算上限 rx*ry<=circle_r^2...
还有其他（“更快”）方法，例如可以使用 CCD 的变化（循环坐标下降）。但是这样的方法通常不能保证会找到最优解，或者根本无法保证......
这里是示例输出的概述：
这些点是来自pnt[n] 的各个点，灰色虚线是边界框并使用外切圆。绿色椭圆是找到解。
【讨论】：
感谢您的帮助。由于我对 C++ 非常陌生，因此我需要一些时间来理解代码。您能否建议我必须使用哪些库才能使代码正常工作。 @Abdul onlymath.h 和链接答案中的 approx 类，您可以直接复制到您的源。【参考方案3】：
这些并没有为您提供 C++ 代码，但它们包含对您需要做的有效算法的深入讨论。
https://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf
http://www.stsci.edu/~RAB/Backup%20Oct%2022%202011/f_3_CalculationForWFIRSTML/Gaertner%20&%20Schoenherr.pdf
【讨论】：
虽然链接到提供答案的信息是一种有用的输入，但 *** 更希望您的答案本身能够提供帮助。这个想法是，答案至少应该包含解决方案的最小摘要，以便在链接失效的情况下它仍然具有一些价值。【参考方案4】：此处的其他答案提供近似方案或仅提供链接。我们可以做得更好。
Gärtner 和 Schönherr (1997) 的论文“Smallest Enclosing Ellipses -- Fast and Exact”解决了您的问题。同一作者在其 1998 年的论文“Smallest Enclosing Ellipses -- An Exact and Generic Implementation in C++”中提供了 C++ 实现。该算法在 CGAL here 中以更实用的形式实现。
但是，CGAL 只提供了椭圆的一般方程，因此我们使用一些变换来获得适合绘图的参数方程。
所有这些都包含在下面的实现中。
使用WebPlotDigitizer 提取数据，同时为轴的长度选择任意值，但保留它们的纵横比，给出：
-1.1314123177813773 4.316368664322679
1.345680085331649 5.1848164974519015
2.2148682495160603 3.9139687117291504
0.9938150357523803 3.2732678860664475
-0.24524315569075128 3.0455750009876343
-1.4493153715482157 2.4049282977126376
0.356472958558844 0.0699802473037554
2.8166270295895384 0.9211630387547896
3.7889384901038987 -0.8484766720657362
1.3457654169794182 -1.6996053411290646
2.9287101489353287 -3.1919219373444463
0.8080480385572635 -3.990389523169913
0.46847074625686425 -4.008682890214516
-1.6521060324734327 -4.8415723146209455
使用下面的程序拟合它给出：
a = 3.36286
b = 5.51152
cx = 0.474112
cy = -0.239756
theta = -0.0979706
然后我们可以用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 EllipseCanonicalEquationa, 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;
【讨论】：
以上是关于在c ++中查找包含一组点的最小面积椭圆的主要内容，如果未能解决你的问题，请参考以下文章
我需要一种算法，可以适应任何大小的n个矩形，在较大的矩形中最小化其面积