如何在凸包中找到最大的三角形,除了蛮力搜索

25

给定一个凸多边形,如何找到定义面积最大的三角形的三个点。

相关问题:这个三角形的外接圆是否也能定义多边形的最小外接圆?


1
不,我正在处理 iPhone 游戏中的多边形形状的碰撞检测。最小外接圆将让我在执行更昂贵的多边形相交测试之前剔除一组潜在的相互碰撞形状。在此过程中,我正在学习计算几何算法并将它们转换为 Objective-C 语言。将来我可能只需要使用一个物理库,但我想从底层了解其工作原理。 - willc2
2
我发现在问题中过多地提供细节是不好的,因为人们往往会回答最有趣(暗示性)的问题,而不是明确陈述的问题,如果它“太简单”的话。作为新手,没有什么是简单的。 - willc2
请注意第二个问题的答案,而不是第一个。抱歉叫你出来,Stephan,这是我试图做两件事的错。 - willc2
5个回答

35

是的,你可以比暴力算法做得更好。

通过暴力算法,我指的是检查所有点的三元组,并选择面积最大的那个。这需要 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

这个算法在Dobkin和Snyder的论文"关于最大化和最小化某些几何问题的一般方法"中得到证明,发表于1979年的FOCS,上面的代码是他们ALGOL-60代码的直接翻译。抱歉使用了while-if-break结构;应该可以将它们转换成更简单的while循环。

1
这真的是O(n)吗? - Ninja420
@Ninja420:是的,没错。我已经在上面的回答中给出了证明;你也可以阅读链接的论文以获取更详细的信息。 - ShreevatsaR
相关问题,如果您能在这里提供详细说明就太好了,https://dev59.com/FnbZa4cB1Zd3GeqPBhJE - Ninja420
1
@Ninja420: 啊,我明白了。是的,当我写答案时,B和C每次只前进2n次的事实对我来说似乎很明显,但现在看起来有点微妙... 显然,我变老了!我会仔细考虑并在我再次清楚时,在另一个问题上发布答案。 - ShreevatsaR
2
@VikashB 不,那不是真的。例如,如果您查看论文中的图表(请参见[此答案](https://dev59.com/FnbZa4cB1Zd3GeqPBhJE#18443566)),似乎并不是多边形直径是三个“锚定局部最大值”之一的任何一条边。更一般地,想象一下正n边形,当n很大时,它几乎类似于一个圆。最大面积的三角形是等边三角形,而不是具有直径作为其中一条边的三角形。 - ShreevatsaR
显示剩余5条评论

6
根据这篇论文,存在一类凸多边形,ShreevatsaR的答案所引用的算法会失败。该论文还提出了一种O(n log n)分治算法来解决这个问题。
显然,对于所有A移动点B和C的简单O(n2)算法仍然有效。

由于这并没有回答问题本身,因此它可能应该作为对@ShreevatsaR答案的评论。 - B. Eckles
2
@B.Eckles 抱歉,这是我第一次贡献。不幸的是,我没有足够的声望来评论答案。由于链接的论文提出了解决问题的算法,所以我认为我的帖子可以作为一个答案。 - tomasf
我理解了! - B. Eckles
1
感谢您的纠正。我之前不知道这个。虽然我还没有尝试或理解反例,但知道它是有争议的很好。 - ShreevatsaR

4

回答你的相关问题:

三角形的外接圆不一定是多边形的最小外接圆。为了证明这一点,考虑一个非常扁平的等腰三角形,例如顶点分别在(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)。


如果我选择较小的一个,那么如何取:最大三角形的外接圆或者两个对踵点之间的圆心? - willc2
很棒,那是一个聪明的反例。我怀疑这样的反例会存在,但是没有想到...顺便说一下,我认为可以证明最小外接圆要么有一些直径对,要么是某个三角形的外接圆。 - ShreevatsaR
是的,我相信那是真的。 - Stephen Canon

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缩小,直到找到最大面积。 我会从中间开始,尝试每个方向寻找更大的面积。


0

我知道这是一个旧帖子,但通过参考上面的答案,我能够修改代码以最大化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);
}

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接