Minimum area quadrilateral algorithm

Question

There are a few algorithms around for finding the minimal bounding rectangle containing a given (convex) polygon.

Does anybody know about an algorithm for finding the minimal-area bounding quadrilateral (any quadrilateral, not just rectangles)?

I've searched the internet for several hours now, but while I found a few theoretical papers on the matter, I did not find a single implementation...

EDIT: People at Mathoverflow pointed me to an article with a mathematical solution (my post there), but for which I did not find an actual implementation. I decided to go with the Monte Carlo Method from Carl, but will dive into the paper and report here, when I have the time...

Thanks all!

Huppah

Source of the problem: I have digitized photos of parts of a product. All the parts are quadrilaterals. The product is flexible (and uneven) and the photos are taken from many different perspectives, so the digitized product parts become polygons. For further computational work on the digitized object, I need to auto-detect the best matching quadrilateral which is per definition of the project the minimum-area quadrilateral. As the photos are hires, a polygon has around 300 vertices. – Huppah

The smallest quadrilateral ? Or one that is small enough ? I suspect that some sort of detection of the edge points followed by clustering into 4 groups followed by line-of-best-fit for each group might work. – Albumenize

In real life "one small enough" will be sufficient. But because I want to know a solution and ideally not iterate through numbers of possibilites, my question is about the smallest quadrilateral... – Huppah

Is there necessarily a unique solution? Consider a regular hexagon with area 3. You can bound it with a square, rhombus, or trapezoid of area 4. – Glyptography

I've racked my brain a little more but I think it's time to throw in the towel. I suspect the folx over at mathoverflow.net will be delighted to jump on this problem, and that is probably the best advice I can offer! – Skirret

However, because you liked my Monte Carlo approach, I came up with a much improved variant of it that you may be interested in looking at. – Skirret

Answered on mathoverflow.net: mathoverflow.net/questions/11580/… – Glyptography

Answer 1

The Monte Carlo approach

Thanks for the clarifying comments on the problem. I've taken away that what's required is not a mathematically correct result but a "fit" that's better than any comparable fits for other shapes.

Rather than pouring a lot of algorithmic brain power at the problem, I'd let the computer worry about it. Generate groups of 4 random points; check that the quad formed by convexly joining those 4 points does not intersect the polygon, and compute the quad's area. Repeat 1 million times, retrieve the quad with the smallest area.

You can apply some constraints to make your points not completely random; this can dramatically improve convergence.

---

Monte Carlo, improved

I've been convinced that throwing 4 points randomly on the plane is a highly inefficient start even for a brute-force solution. Thus, the following refinement:

• For each trial, randomly select p distinct vertices and q distinct sides of the polygon such that p + q = 4.
• For each of the q sides, construct a line passing through that side's endpoints.
• For each of the p vertices, construct a line passing through that vertex and with a randomly assigned slope.
• Verify that the 4 lines indeed form a quadrilateral, and that this quadrilateral contains (and does not intersect!) the polygon. If these tests fail, don't pursue this iteration any further.
• If this quadrilateral's area is the minimum of all areas seen so far, remember the area and the coordinates of the quadrilateral's vertices.
• Repeat an arbitrary number of times, and return the "best" quadrilateral found.

As opposed to always requiring 8 random numbers (x and y coordinates for each of 4 points), this solution requires only (4 + p) random numbers. Also, the lines produced are not blindly floundering in the plane but are each touching the polygon. This ensures that the quadrilaterals are from the outset at least very close to the polygon.

Skirret

For real life Monte Carlo might and will work. But if it's possible to compute an exact solution, I would prefer it (also out of "scientific curiosity" I'd like to know)... – Huppah

Perhaps obvious, but there are some simple optimizations you can apply to the quadrilateral afterward to make sure you haven't left any low-hanging fruit. Each edge should touch the original polygon. A single side can be wiggled back and forth without changing the slope of any other sides. The same ideas could be applied after any approximation algorithm. – Glyptography

@Jason Orendorff: Right you are. Actually, I think the more effective course would be to combine my other solution with the Monte Carlo approach, so as to reap the benefits of your optimizations earlier. I've updated my answer to that effect. Thanks for your input! – Skirret

I guess you could now remove the element of randomness and you'd have an O(n^4) algorithm. Try every combination of vertices and sides. For each combination that makes a quadrilateral, write a formula for the area in terms of the p slope values and solve for the values that produce the minimum area. – Glyptography

Hmm. With p=1 that would be trivial. But with p=2 we'd already have to find valleys on a hilly plane, and with p=3 or p=4 we'd have 3, 4 unknowns... and I'm not clever enough to find bumps in a 5-dimensional solid. I'm not sure there are numeric solutions to that kinda stuff. – Skirret

But picking sides and vertices, in and of itself, is really only O(n^2). That's relatively harmless. – Skirret

If there's not a closed-form solution, just point a general-purpose solver at it; that'll be better than Monte Carlo at least. But I'm not convinced there's no closed-form solution. – Glyptography

OK, agreed. So in essence we're predigested the problem to a form that can hopefully be solved "mechanically", perhaps even with an exact result. This should warm Carsten's heart. – Skirret

But my guess is that he visited the smart people at MathOverflow and they laughed so hard at our "solutions" that he's stopped talking to us ;) – Skirret

Just noticed that any 4-tuple of angles (a, b, c, d), such that 0 <= a < b < c < d < 2π, b - a < π, c - b < π, d - c < π, and 2π+a - d < π, uniquely identifies a bounding quadrilateral whose sides make those four angles with the horizontal. Write a function that calculates the area of that quadrilateral given the tuple. How numerically well-behaved is this function? I don't know much about this. But I speculate you can point a solver at it and get excellent results. – Glyptography

@Carl, @All: I went indeed to MathOverFlow and was presented a paper with a solution, but had no time to go into the details of it yet. For the real-world app I decided earlier with the enhanced Monte Carlo method, because it's probably easier to implement. Thanks all for helping! – Huppah

P should always = 0. If you pick a line passing through a vertex, it cannot be a better fit than both of lines passing through the adjacent sides. – Packthread

That would be great news if it were true. Do you have any proof, reference or other information that leads to this statement, i.e. any way for me and others to verify it's not just your intuition? – Skirret

There are 2 triangles formed by the moving line , the two adjacent polygon edges, and the extended edges of the two hypothetical quad-sides that border the quad-side we're choosing now. Basically, you're picking a linear combination between two choices of triangle areas. – Packthread

Answer 2

A stingy algorithm

Start with your convex polygon. While there are more than 4 points, find the pair of neighboring points that's cheapest to consolidate, and consolidate them, reducing the number of points in your polygon by 1.

By "consolidate", I just mean extending the two edges on either side until they meet at a point, and taking that point to replace the two.

By "cheapest", I mean the pair for which consolidation adds the smallest amount of area to the polygon.

Before:                       After consolidating P and Q:

                                    P'
                                   /\
   P____Q                         /  \
   /    \                        /    \
  /      \                      /      \
 /        \                    /        \

Runs in O(n log n). But this produces only an approximation, and I'm not entirely happy with it. The better the algorithm does at producing a nice regular pentagon, the more area the last consolidation must eat up.

Glyptography

Are you sure that it's only an approximation? It's not obvious to me either way. – Quetzal

Well, take a regular octagon and distort it slightly so the algorithm is sure to consolidate points 1 and 2 first, then points 4 and 5. The algorithm produces a pointy kite shape instead of the desired near-square. – Glyptography

Ah, yes. That's a nice example. – Quetzal

This works quite well, with one implementation gotcha to be aware of. P' can be on the 'wrong side' of PQ in some cases. I added an explicit side of line check. I presume if you're using signed triangle area, a negative check would work also. .P---------Q .\ / . \ / Edit can't make the formatting work, but if the edges leading in to P&Q come together, the point P' crosses the line PQ – Jannet

Answer 3

Here's an algorithm to fit quadrangles to arbitrary masks using the technique from View Frustum Optimization To Maximize Object’s Image Area.

---

Here's the output -

---

import cv2
import numpy as np
import sympy


def appx_best_fit_ngon(mask_cv2, n: int = 4) -> list[(int, int)]:
    # convex hull of the input mask
    mask_cv2_gray = cv2.cvtColor(mask_cv2, cv2.COLOR_RGB2GRAY)
    contours, _ = cv2.findContours(
        mask_cv2_gray, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE
    )
    hull = cv2.convexHull(contours[0])
    hull = np.array(hull).reshape((len(hull), 2))

    # to sympy land
    hull = [sympy.Point(*pt) for pt in hull]

    # run until we cut down to n vertices
    while len(hull) > n:
        best_candidate = None

        # for all edges in hull ( <edge_idx_1>, <edge_idx_2> ) ->
        for edge_idx_1 in range(len(hull)):
            edge_idx_2 = (edge_idx_1 + 1) % len(hull)

            adj_idx_1 = (edge_idx_1 - 1) % len(hull)
            adj_idx_2 = (edge_idx_1 + 2) % len(hull)

            edge_pt_1 = sympy.Point(*hull[edge_idx_1])
            edge_pt_2 = sympy.Point(*hull[edge_idx_2])
            adj_pt_1 = sympy.Point(*hull[adj_idx_1])
            adj_pt_2 = sympy.Point(*hull[adj_idx_2])

            subpoly = sympy.Polygon(adj_pt_1, edge_pt_1, edge_pt_2, adj_pt_2)
            angle1 = subpoly.angles[edge_pt_1]
            angle2 = subpoly.angles[edge_pt_2]

            # we need to first make sure that the sum of the interior angles the edge
            # makes with the two adjacent edges is more than 180°
            if sympy.N(angle1 + angle2) <= sympy.pi:
                continue

            # find the new vertex if we delete this edge
            adj_edge_1 = sympy.Line(adj_pt_1, edge_pt_1)
            adj_edge_2 = sympy.Line(edge_pt_2, adj_pt_2)
            intersect = adj_edge_1.intersection(adj_edge_2)[0]

            # the area of the triangle we'll be adding
            area = sympy.N(sympy.Triangle(edge_pt_1, intersect, edge_pt_2).area)
            # should be the lowest
            if best_candidate and best_candidate[1] < area:
                continue

            # delete the edge and add the intersection of adjacent edges to the hull
            better_hull = list(hull)
            better_hull[edge_idx_1] = intersect
            del better_hull[edge_idx_2]
            best_candidate = (better_hull, area)

        if not best_candidate:
            raise ValueError("Could not find the best fit n-gon!")

        hull = best_candidate[0]

    # back to python land
    hull = [(int(x), int(y)) for x, y in hull]

    return hull

Here's the code I used to generate this image -

hull = appx_best_fit_ngon(mask_cv2)

for idx in range(len(hull)):
    next_idx = (idx + 1) % len(hull)
    cv2.line(mask_cv2, hull[idx], hull[next_idx], (0, 255, 0), 1)

for pt in hull:
    cv2.circle(mask_cv2, pt, 2, (255, 0, 0), 2)

Kimon

Answer 4

I think a 2D OBB around the points is a good starting place. That will probably give a "good" (but not best) fit; if you find you still need a tighter bound, you can extend the fit to quadrilaterals.

First off, you should compute the convex hull of your input points. That gives you fewer points to deal with, and doesn't change the answer at all.

I'd stay away from the covariance-based method that's referenced on Gottschalk's paper elsewhere. That doesn't always give the best fit, and can go very wrong for very simple input (e.g. a square).

In 2D, the "rotating calipers" approach described at http://cgm.cs.mcgill.ca/~orm/maer.html should give the exact best-fit OBB. The problem is also easy to think about as a 1D minimization problem:

1. For a given angle, rotate the x and y axes by this angle. This gives you two orthogonal vectors.
2. For each point, project onto both axes. Keep track of the min and max projection on each axis.
3. The area of the OBB is (max1-min1)*(max2-min2), and you can easily compute the points of the OBB from the angle and the min and max projections.
4. Repeat, sampling the interval [0, pi/2] as finely as you want, and keeping track of the best angle.

John Ratcliffe's blog (http://codesuppository.blogspot.com/2006/06/best-fit-oriented-bounding-box.html) has some code that does this sampling approach in 3D; you should be able to adapt it to your 2D case if you get stuck. Warning: John sometimes posts mildly NSFW pictures on his blog posts, but that particular link is fine.

If you're still not happy with the results after getting this working, you could modify the approach a bit; there are two improvements that I can think of:

1. Instead of picking orthogonal axes as mentioned above, you could pick two non-orthogonal axes. This would give you the best-fitting rhombus. To go about this, you'd essentially do a 2D search over [0, pi] x [0, pi].
2. Use the best-fit OBB as the starting point for a more detailed search. For example, you could take the 4 points from the OBB, move one of them until the line touches a point, and then repeat for the other points.

Hope that helps. Sorry the information overload :)

Delgado

Answer 5

I think the oriented bounding box should be pretty close (though it's actually a rectangle). Here is the standard reference paper on oriented bounding boxes: Gottschalk's paper (PDF)

Look at section 3 (Fitting an OBB).

Biannual

Ah, I think I found the weakness there: The paper qualifies OBBs as rectanguloids (in 3-space), which would be equivalent to rectangles in 2D. But this would be a restriction from the OP's requirement that the bounding quadrilaterals in fact be general quadrilaterals. I suspect oriented rectangles would be considerably easier to compute, so the algorithms given in the paper would not be general enough. – Skirret

Answer 6

I believe that each side of the minimal-area bounding quadrilateral will pass through at least 2 of the vertices of your polygon. If this assumption is true, it should be easy to perform a brute force search for the solution.

• Generate the unique set of lines which are defined by 2 of your vertices and do not intersect the polygon.
• Examine each set of 4 lines to determine which produces the minimal-area bounding quadrilateral.

UPDATED: The assumption that each side of the bounding quad will pass through at least 2 vertices is false, but a related set of lines may provided a basis for a solution. At the very least, each side of the bounding quad will pass through one of the vertices used to define the unique set of the lines which are defined by 2 vertices and do not intersect the polygon.

Haun

Upon further reflection, my assumption may not even hold true for a regular pentagon, so this probably isn't the answer. – Haun

Answer 7

Here's an observation that leads to an improvement to the Monte Carlo algorithm, and may lead to a direct answer as well.

Lemma: If an edge of an optimal quadrilateral touches the polygon at only one point, it touches at the midpoint of that edge.

Proof: Define X and Y as the lengths of the two segments on either side of the touching point. Imagine rotating the edge about the single touching point by an infinitesimal angle A. The rotation affects the size of the quadrilateral by increasing it by XA/2 and decreasing it by YA/2, or vice versa. If X != Y, then one of the two rotation directions leads to a smaller quadrilateral. Since the quadrilateral is minimal, we must have X=Y.

To use this fact, note that if we pick some edges and points where the polygon touches the quadrilateral, and we don't pick two points in a row, we can uniquely determine the quadrilateral by picking the edge through each point which makes that point the midpoint of the edge (and if that isn't possible, reject the configuration that was picked). In the Monte Carlo algorithm, this reduces to saying that we don't have to pick the slope for this edge - it can be determined explicitly.

We still have the case where two adjacent points were picked - hopefully I've inspired someone else to pick up here...

Istle

Answer 8

I have the same problem to solve and the code I am using is actually implementing the idea of Jason Orendorff with one additional rectangle, which is bounding the process and makes the result more square like. At the end it is a good heuristic which works well in my case. I hope somebody else could also benefit from this code:

import java.util.ArrayList;
import java.util.List;

import org.opencv.core.Point;
import org.opencv.core.Rect;


public class Example {
    
    
    public static Point[] getMinimalQuadrilateral(Point[] convexPolygon, Rect boundingRec) {
        
        if (convexPolygon.length <= 4) {
            throw new IllegalStateException();
        }
        
        //Create list with all entries
        List<ListItem<Point>> all_init_list = new ArrayList<ListItem<Point>>();
        for (int i = 0; i < convexPolygon.length; i++) {
            ListItem<Point> cur = new ListItem<Point>();
            cur.value = convexPolygon[i];
            all_init_list.add(cur);
        }
        
        //Link the list
        for (int i = 0; i < all_init_list.size() - 1; i++) {
            all_init_list.get(i).next = all_init_list.get(i + 1);
        }
        //Make it cyclic
        all_init_list.get(all_init_list.size() - 1).next = all_init_list.get(0);
        
        
        int countOfPoints = all_init_list.size();
        ListItem<Point> start = all_init_list.get(0);
        
        while (countOfPoints > 4) {
            
            //System.out.println("countOfPoints=" + countOfPoints);
            
            double minTriangleArea = Double.MAX_VALUE;
            ListItem<Point> best = null;
            ListItem<Point> best_intersection = new ListItem<Point>();
            ListItem<Point> cur = start;
            do {
                Point p1 = cur.value;
                Point p2 = cur.next.value;
                Point p3 = cur.next.next.value;
                Point p4 = cur.next.next.next.value;
                
                //Do work
                Point intersection = findIntersection(p1, p2, p4, p3);
                if (intersection != null && boundingRec.contains(intersection)) {
                    double cur_area = triangleArea(p2, intersection, p3);
                    if (cur_area < minTriangleArea) {
                        minTriangleArea = cur_area;
                        best = cur;
                        best_intersection.value = intersection;
                        //System.out.println("minTriangleArea=" + minTriangleArea);
                    }
                }
                
                cur = cur.next;
            } while (cur != start);
            
            //If there is best than remove 2 points and put their intersection instead
            if (best == null) {
                break;
            }
            best_intersection.next = best.next.next.next;
            best.next = best_intersection;
            countOfPoints--;
            start = best;
        }
        
        //Compose result
        Point[] result = new Point[countOfPoints];
        while (countOfPoints > 0) {
            result[countOfPoints - 1] = start.value;
            start = start.next;
            countOfPoints--;
        }
        
        return result;
    }
    
    
    public static double triangleArea(Point A, Point B, Point C) {
        double area = (A.x * (B.y - C.y) + B.x * (C.y - A.y) + C.x * (A.y - B.y)) / 2.0;
        return Math.abs(area);
    }
    
    
    public static Point findIntersection(Point l1s, Point l1e, Point l2s, Point l2e) {
        
        double a1 = l1e.y - l1s.y;
        double b1 = l1s.x - l1e.x;
        double c1 = a1 * l1s.x + b1 * l1s.y;
        
        double a2 = l2e.y - l2s.y;
        double b2 = l2s.x - l2e.x;
        double c2 = a2 * l2s.x + b2 * l2s.y;
        
        double delta = a1 * b2 - a2 * b1;
        if (delta == 0) {
            return null;
        }
        
        return new Point((b2 * c1 - b1 * c2) / delta, (a1 * c2 - a2 * c1) / delta);
    }
    
    
    private static final class ListItem<T> {
        public T value;
        public ListItem<T> next;
    }
}

The algorithm could be further improved by for example starting with small bounding rectangle and sequentially increasing it, if the process is not able to find a solution. In practice I am using bounding rectangle with 5% bigger size that the minimal bounding rectangle.

Swarts