如何在蛮力搜索之外找到凸包中的最大三角形

Posted 2021-04-12

给定凸多边形，如何找到定义具有最大面积的三角形的3个点。

相关：该三角形的外接圆是否也会定义多边形的最小边界圆？

答案

是的，你可以做得比蛮力更好。

通过蛮力，我认为你的意思是检查所有三点，并选择一个最大面积的点。这在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前进到C = 3，C = 4，...），只要三角形的面积增加，即，只要面积（A，B，C）≤面积（A，B，C + 1）。这个点C将是那些固定的A和B最大化区域（ABC）的点。（换句话说，函数Area（ABC）是C的函数的单峰）。）

接下来，如果增加面积，则前进B（不改变A和C）。如果是这样，再次按上述方式推进C.然后如果可能的话再次前进B等。这将给出最大区域三角形，其中A作为顶点之一。

（到这里的部分应该很容易证明，并且简单地为每个A单独执行此操作将给出O（n2）算法。）

现在再次推进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
   if A==B: B = (B+1)%n
   if B==C: C = (C+1)%n
   if A==0: break

该算法在Dobkin和Snyder，On a general method for maximizing and minimizing among certain geometric problems，FOCS 1979中得到证明，上面的代码是他们的ALGOL-60代码的直接翻译。为休息时间结构道歉;应该可以将它们转换为更简单的while循环。

另一答案

回答你的相关问题：

三角形的外接圆不一定是多边形的最小边界圆。为了看到这一点，考虑一个非常平坦的等腰三角形，比如顶点为（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）处的点。

另一答案

根据this论文，有一类凸多边形，其中ShreevatsaR的答案引用的算法失败了。本文还提出了一种用于解决问题的O（n log n）分而治之算法。

显然，用于移动所有A的点B和C的更简单的O（n2）算法仍然有效。

另一答案

来自http://www.wolframalpha.com/input/?i=triangle三角形的面积= sqrt（（a + bc）（a-b + c）（ - a + b + c）*（a + b + c））/ 4如果使用c连接到终点你的凸多边形，如果a和b会碰到你的凸多边形，你可以围绕你的多边形迭代，允许a增长，b缩小，直到找到你的最大面积。我会从中间点开始尝试每个方向以获得更大的区域。

另一答案

我知道这是一个旧帖子，但通过参考answer above，我能够修改代码以最大化n边多边形的面积。

注意：凸壳是使用Emgu OpenCV library发现的。我正在使用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);
}