正如 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] 的各个点,灰色虚线是边界框并使用外切圆。绿色椭圆是找到解。