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

Posted 2023-02-22

【中文标题】除了蛮力搜索之外，如何在凸包中找到最大的三角形

【英文标题】：How to find largest triangle in convex hull aside from brute force search

【发布时间】：2009-10-25 16:46:20

【问题描述】：

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

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

【问题讨论】：

不，我正在为 iPhone 游戏进行多边形形状的碰撞检测。最小边界圆可以让我在进行更昂贵的多边形-多边形相交测试之前剔除一组可能发生碰撞的形状。在这个过程中，我正在学习计算几何算法并将它们翻译成 Objective-C。将来我可能只会使用物理库，但我想知道它是如何从头开始工作的。

我发现在问题中加入太多细节是不好的，因为如果“太简单”，人们往往会回答最有趣（隐含）的问题，而不是陈述的问题。作为一个新手，没有什么是简单的。

注意次要问题的答案，而不是主要问题的答案。很抱歉把你叫出来斯蒂芬，这是我的错，我尝试了两次。

【参考方案1】：

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

通过蛮力，我假设您的意思是检查所有三重点，并选择具有最大面积的点。这在 O(n³) 时间内运行，但事实证明不仅可以在 O(n²) 时间内完成但在 O(n) 时间内！

[更新： 2017 年发表的一篇论文通过示例表明 O(n) 解决方案不适用于特定类别的多边形。有关更多详细信息，请参阅this 答案。但是 O(n²) 算法还是正确的。]

如果有必要，首先对点进行排序/计算凸包（在 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)。对于那些固定的 A 和 B，该点 C 将最大化 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
   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-if-break 结构道歉；应该可以将它们转换为更简单的 while 循环。

【讨论】：

这真的是 O(n) 吗？

@Ninja420：是的。我已经在上面的答案中给出了证明；您还可以阅读链接的论文以获取更详细的内容。

相关问题，如果你能在这里提供详细的描述就好了，***.com/questions/18423040/…

@Ninja420：啊，我明白了。是的，当我写下答案时，B 和 C 各自只前进 2n 次这一事实对我来说一定是显而易见的，但现在看起来很微妙……显然我变老了！当我再次清楚时，我会考虑并在另一个问题上发布答案。

@VikashB 不，那不是真的。例如。如果您查看论文中的图（请参阅this answer），则多边形的直径似乎不是三个“锚定局部最大值”中的任何一个的边。更一般地，想象一个用于大 n 的规则 n 边形，它几乎类似于一个圆。面积最大的三角形是等边三角形，而不是直径作为边之一的三角形。

【参考方案2】：

根据this 论文，有一类凸多边形，其中 ShreevatsaR 的答案引用的算法失败了。论文还提出了一种O(n log n)分治算法来解决这个问题。

显然，对于 all A 移动点 B 和 C 的更简单的 O(n²) 算法仍然有效。

【讨论】：

由于这不是对问题本身的回应，它可能应该作为@ShreevatsaR 对答案的评论

@B.Eckles 抱歉，这是我第一次投稿。不幸的是，我没有足够的声誉来评论答案。而且由于链接的论文提出了解决问题的算法，我认为我的帖子有资格作为答案。

对我来说很有意义！

感谢指正。我不知何故错过了这个早些时候。我还没有尝试或理解反例，但很高兴知道它是有争议的。

【参考方案3】：

回答您的相关问题：

三角形的外接圆不一定是多边形的最小外接圆。要看到这一点，请考虑一个非常平坦的等腰三角形，例如顶点位于 (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) 处的点。

【讨论】：

如果我取较小的那个：（最大三角形的）外圆或以对映点为中心的圆？

酷，这是一个聪明的反例。我怀疑这样的存在，但没有想出一个......顺便说一句，我认为有可能证明最小边界圆要么有一些直径，要么是 some i> 三角形。

是的，我相信这是真的。

【参考方案4】：

来自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 缩小，直到找到最大面积。 我会从中间点开始，尝试每个方向以获得更大的区域。

【参考方案5】：

我知道这是一篇旧帖子，但通过引用 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);