How to find largest triangle in convex hull aside from brute force search

Question

Given a convex polygon, how do I find the 3 points that define a triangle with the greatest area.

Related: Is it true that the circumcircle of that triangle would also define the minimum bounding circle of the polygon?

Answer 1

Yes, you can do significantly better than brute-force.

By brute-force I assume you mean checking all triples of points, and picking the one with maximum area. This runs in O(n³) time, but it turns out that it is possible to do it in not just O(n²) but in O(n) time!

[Update: A paper published in 2017 showed by example that the O(n) solution doesn't work for a specific class of polygons. See this answer for more details. But the O(n²) algorithm is still correct.]

By first sorting the points / computing the convex hull (in O(n log n) time) if necessary, we can assume we have the convex polygon/hull with the points cyclically sorted in the order they appear in the polygon. Call the points 1, 2, 3, … , n. Let (variable) points A, B, and C, start as 1, 2, and 3 respectively (in the cyclic order). We will move A, B, C until ABC is the maximum-area triangle. (The idea is similar to the rotating calipers method, as used when computing the diameter (farthest pair).)

With A and B fixed, advance C (e.g. initially, with A=1, B=2, C is advanced through C=3, C=4, …) as long as the area of the triangle increases, i.e., as long as Area(A,B,C) ≤ Area(A,B,C+1). This point C will be the one that maximizes Area(ABC) for those fixed A and B. (In other words, the function Area(ABC) is unimodal as a function of C.)

Next, advance B (without changing A and C) if that increases the area. If so, again advance C as above. Then advance B again if possible, etc. This will give the maximum area triangle with A as one of the vertices.

(The part up to here should be easy to prove, and simply doing this separately for each A would give an O(n²) algorithm.)

Now advance A again, if it improves the area, and so on._{(The correctness of this part is more subtle: Dobkin and Snyder gave a proof in their paper, but a recent paper shows a counterexample. I have not understood it yet.)}

Although this has three "nested" loops, note that B and C always advance "forward", and they advance at most 2n times in total (similarly A advances at most n times), so the whole thing runs in O(n) time.

Code fragment, in Python (translation to C should be straightforward):

 # 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

This algorithm is proved in Dobkin and Snyder, On a general method for maximizing and minimizing among certain geometric problems, FOCS 1979, and the code above is a direct translation of their ALGOL-60 code. Apologies for the while-if-break constructions; it ought to be possible to transform them into simpler while loops.

Answer 2

According to this paper, there is a class of convex polygons in which the algorithm cited by ShreevatsaR's answer fails. The paper also proposes a O(n log n) divide and conquer algorithm for solving the problem.

Apparently, the simpler O(n²) algorithm in which you move points B and C for all A is still valid.

Answer 3

answering your related question:

The circumcircle of the triangle is not necessarily the minimum bounding circle of the polygon. To see this, consider a very flat isosceles triangle, say with vertices at (0,0), (10,0) and (5,1). The minimum bounding circle has center (5,0) and radius 5, but this circle doesn't touch the vertex at (5,1), so it's not the circumcircle. (The circumcircle has center (5,-12) and radius 13)

edit:

Choosing the smaller of the circumcircle or the circle containing the antipodal points of the diameter of the polygon also doesn't suffice, because it is possible to construct polygons that have points outside the circumcircle of the maximal triangle. Consider the pentagon with vertices at:

(-5,  0)
(-4, -1)
( 5,  0)
( 4,  1)
(-4,  1)

The maximal triangle has vertices at (-4,-1), (5, 0), and (-4, 1). Its circumcircle does not include the point at (-5, 0).

Answer 4

from http://www.wolframalpha.com/input/?i=triangle The area of the triangle = sqrt((a+b-c)(a-b+c)(-a+b+c)*(a+b+c)) / 4 If you use c connected to the end points of your convex polygon and if a and b would touch your convex polygon you could iterate around your polygon allowing a to grow and b to shrink until you find your maximum area. I would start mid point and try each direction for a larger area.

Answer 5

I know this is an old post but by referencing the answer above I was able to modify the code to maximize the area for an n-sided polygon.

Note: The convex hull was found using Emgu OpenCV library. I'm using CvInvoke.ContourArea() method to calculate the given area of a polygon. This is written in C#. It also assumes that the points are arranged in clockwise order. This can be specified in the method 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;
}

The area method used. This could change depending on what is needed to maximize.

private double GetArea(PointF[] points)
{
    return CvInvoke.ContourArea( new Emgu.CV.Util.VectorOfPointF(points),false);
}