给定一个凸多边形,如何找到定义面积最大的三角形的三个点。
相关问题:这个三角形的外接圆是否也能定义多边形的最小外接圆?
给定一个凸多边形,如何找到定义面积最大的三角形的三个点。
相关问题:这个三角形的外接圆是否也能定义多边形的最小外接圆?
是的,你可以比暴力算法做得更好。
通过暴力算法,我指的是检查所有点的三元组,并选择面积最大的那个。这需要 O(n3) 时间,但事实证明,不仅可以在 O(n2) 的时间内完成,还可以在 O(n) 时间内完成!
[更新:2017年发表的一篇论文举例说明了O(n)解决方案在特定类别的多边形中无效。有关详细信息,请参见this答案。但是O(n2)算法仍然是正确的。]
首先通过排序点/计算凸包(在O(n log n)时间内)来确保我们有按顺序排列出现在多边形中的凸多边形/凸壳。将点1、2、3、...、n称为(变量)点A、B和C,分别从1、2和3开始(按循环顺序)。我们将移动A、B和C,直到ABC成为最大面积的三角形。(这个想法类似于rotating calipers方法,在计算diameter (farthest pair)时使用。)
在固定A和B的情况下,递增C(例如,初始时,当A=1,B=2时,通过C = 3、C = 4等方式前进),只要三角形的面积增加,即只要Area(A,B,C)≤Area(A,B,C+1)。这个点C将是对于那些固定的A和B使Area(ABC)最大化的点。(换句话说,函数Area(ABC)作为C的函数是单峰的。)
接下来,只要它增加了面积,就递增B(而不改变A和C)。如果是这样,请再次按上述方式递增C。然后,如果可能的话,再次递增B等。这将给出以A为其中一个顶点的最大面积三角形。
(到目前为止的部分应该很容易证明,并且仅为每个A分别执行此操作将给出O(n²)算法。)
现在如果提高A可以改善面积,依此类推。
(这部分的正确性更加微妙:Dobkin和Snyder在他们的论文中给出了证明,但最近的一篇论文显示了一个反例。我还没有理解它。)虽然这有三个“嵌套”的循环,但是请注意,B和C总是“向前”前进,并且它们最多在总共2n次(同样,A最多n次)内前进,因此整个过程需要O(n)时间。
Python的代码片段(翻译成C应该很容易):
# Assume points have been sorted already, as 0...(n-1)
A = 0; B = 1; C = 2
bA= A; bB= B; bC= C #The "best" triple of points
while True: #loop A
while True: #loop B
while area(A, B, C) <= area(A, B, (C+1)%n): #loop C
C = (C+1)%n
if area(A, B, C) <= area(A, (B+1)%n, C):
B = (B+1)%n
continue
else:
break
if area(A, B, C) > area(bA, bB, bC):
bA = A; bB = B; bC = C
A = (A+1)%n
B = (A+1)%n
C = (A+2)%n
if A==0: break
回答你的相关问题:
三角形的外接圆不一定是多边形的最小外接圆。为了证明这一点,考虑一个非常扁平的等腰三角形,例如顶点分别在(0,0)、(10,0)和(5,1)处。最小外接圆的圆心为(5,0),半径为5,但是这个圆不与顶点(5,1)相切,因此它不是外接圆。(外接圆的圆心为(5,-12),半径为13)
编辑:
选择外接圆或包含多边形直径对踵点的圆中较小的那一个也不够,因为有可能构造出多边形,其某些点位于最大三角形的外接圆之外。例如考虑五边形,其顶点为:
(-5, 0)
(-4, -1)
( 5, 0)
( 4, 1)
(-4, 1)
最大三角形的顶点分别为(-4,-1)、(5,0)和(-4,1)。其外接圆不包括点(-5,0)。
来自http://www.wolframalpha.com/input/?i=triangle的信息: 三角形的面积 = sqrt((a+b-c)(a-b+c)(-a+b+c)*(a+b+c)) / 4 如果你将c连接到凸多边形的端点, 并且如果a和b会接触到你的凸多边形, 你可以在多边形周围迭代, 让a增长,b缩小,直到找到最大面积。 我会从中间开始,尝试每个方向寻找更大的面积。
我知道这是一个旧帖子,但通过参考上面的答案,我能够修改代码以最大化n边形的面积。
注意:使用Emgu OpenCV库找到了凸包。我正在使用CvInvoke.ContourArea()
方法来计算多边形的给定面积。这是用C#编写的。它还假定点按顺时针顺序排列。这可以在方法CvInvoke.ConvexHull()
中指定。
private PointF[] GetMaxPolygon(PointF[] convexHull, int vertices)
{
// validate inputs
if(convexHull.Length < vertices)
{
return convexHull;
}
int numVert = vertices;
// triangles are the smalles polygon with an area.
if (vertices < 3)
numVert = 3;
PointF[] best = new PointF[numVert]; // store the best found
PointF[] next = new PointF[numVert]; // test collection of points to compare
PointF[] current = new PointF[numVert]; // current search location.
int[] indexes = new int[numVert]; // map from output to convex hull input.
int[] nextIndices = new int[numVert];
//starting values 0,1,2,3...n
for(int i = 0; i < numVert; i++)
{
best[i] = convexHull[i];
next[i] = convexHull[i];
current[i] = convexHull[i];
}
// starting indexes 0,1,2,3... n
for(int i = 0; i < numVert; i++)
{
nextIndices[i] = i;
indexes[i] = i;
}
// starting areas are equal.
double currentArea = GetArea(current);
double nextArea = currentArea;
int exitCounter = 0;
while(true)
{
// equivelant to n-1 nested while loops
for(int i = numVert - 1; i > 0; i--)
{
while (exitCounter < convexHull.Length)
{
// get the latest area
currentArea = GetArea(current);
nextIndices[i] = (nextIndices[i] + 1) % convexHull.Length;
next[i] = convexHull[nextIndices[i]]; // set the test point
nextArea = GetArea(next);
if (currentArea <= nextArea) // compare.
{
indexes[i]= (indexes[i]+1) % convexHull.Length;
current[i] = convexHull[indexes[i]];
currentArea = GetArea(current);
exitCounter++; // avoid infinite loop.
}
else //stop moving that vertex
{
for(int j = 0; j< numVert; j++)
{
nextIndices[j] = indexes[j];
next[j] = convexHull[indexes[j]];//reset values.
}
break;
}
}
}
// store the best values so far. these will be the result.
if(GetArea(current)> GetArea(best))
{
for (int j = 0; j < numVert; j++)
{
best[j] = convexHull[indexes[j]];
}
}
// The first index is the counter. It should traverse 1 time around.
indexes[0] = (indexes[0] + 1) % convexHull.Length;
for(int i = 0; i < vertices-1;i++)
{
if(indexes[i] == indexes[i+1])// shift if equal.
{
indexes[i + 1] = (indexes[i + 1] + 1) % convexHull.Length;
}
}
//set new values for current and next.
for(int i = 0; i < numVert; i++)
{
current[i] = convexHull[indexes[i]];
next[i] = convexHull[indexes[i]];
}
// means first vertex finished traversing the whole convex hull.
if (indexes[0] == 0)
break;
}
return best;
}
使用的区域方法。这可能会根据需要进行最大化的变化。
private double GetArea(PointF[] points)
{
return CvInvoke.ContourArea( new Emgu.CV.Util.VectorOfPointF(points),false);
}