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!

Answer 1

I have implementation in Kotlin that worked perfectly for my case:

import androidx.compose.ui.geometry.Offset
import kotlin.math.abs
import kotlin.math.sqrt

object GetBoundaryBorder {

fun getPolygonsHull(polygons: List<List<Offset>>, n: Int = 4): List<Offset> {
    val allPoints = polygons.flatten()
    val initialHull = findConvexHull(allPoints)
    return findPolygonHull(initialHull, n)
}

fun getPolygonHull(polygons: List<Offset>, n: Int = 4): List<Offset> {
    val initialHull = findConvexHull(polygons)
    return findPolygonHull(initialHull, n)
}

private fun findConvexHull(points: List<Offset>): List<Offset> {
    if (points.size < 3) return points

    val start = points.minWith(compareBy({ it.y }, { it.x }))
    val sortedPoints = points.filter { it != start }.sortedWith(compareBy(
        { point ->
            val dx = point.x - start.x
            val dy = point.y - start.y
            Math.atan2(dy.toDouble(), dx.toDouble())
        },
        { point ->
            val dx = point.x - start.x
            val dy = point.y - start.y
            dx * dx + dy * dy
        }
    ))

    val stack = mutableListOf(start)
    for (point in sortedPoints) {
        while (stack.size >= 2 && !isCounterClockwise(
                stack[stack.size - 2],
                stack[stack.size - 1],
                point
            )
        ) {
            stack.removeAt(stack.size - 1)
        }
        stack.add(point)
    }

    return stack
}

private fun isCounterClockwise(p1: Offset, p2: Offset, p3: Offset): Boolean {
    val crossProduct = (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x)
    return crossProduct > 0
}

private fun findPolygonHull(hullOffsets: List<Offset>, numberOfSides: Int = 4): List<Offset> {
    // Handle special cases
    when (hullOffsets.size) {
        0, 1 -> return hullOffsets
        2 -> return expandLineToQuad(hullOffsets[0], hullOffsets[1])
        3 -> return expandTriangleToQuad(hullOffsets)
    }

    var currentHull = hullOffsets
    while (currentHull.size > numberOfSides) {
        var bestCandidate: Pair<List<Offset>, Double>? = null

        for (edgeIdx1 in currentHull.indices) {
            val edgeIdx2 = (edgeIdx1 + 1) % currentHull.size
            val adjIdx1 = (edgeIdx1 - 1 + currentHull.size) % currentHull.size
            val adjIdx2 = (edgeIdx1 + 2) % currentHull.size

            val edgePt1 = currentHull[edgeIdx1]
            val edgePt2 = currentHull[edgeIdx2]
            val adjPt1 = currentHull[adjIdx1]
            val adjPt2 = currentHull[adjIdx2]

            val intersection = lineIntersectionBorder(adjPt1, edgePt1, edgePt2, adjPt2) ?: continue
            val area = triangleAreaBorder(edgePt1, intersection, edgePt2)

            if (bestCandidate != null && bestCandidate.second < area) continue

            val betterHull = currentHull.toMutableList()
            betterHull[edgeIdx1] = intersection
            betterHull.removeAt(edgeIdx2)
            bestCandidate = Pair(betterHull, area)
        }

        bestCandidate?.let {
            currentHull = it.first.toMutableList()
        } ?: break  // If we can't find a valid candidate, break instead of throwing exception
    }

    return currentHull
}

private fun expandTriangleToQuad(triangle: List<Offset>): List<Offset> {
    // Find the longest edge
    val edges = listOf(
        Pair(0, 1),
        Pair(1, 2),
        Pair(2, 0)
    )

    // Calculate edge lengths and find longest
    val edgeLengths = edges.map { (i, j) ->
        Triple(i, j, distance(triangle[i], triangle[j]))
    }
    val longestEdge = edgeLengths.maxBy { it.third }

    // Find the point not on the longest edge
    val oppositePointIndex = (0..2).first { it != longestEdge.first && it != longestEdge.second }
    val oppositePoint = triangle[oppositePointIndex]

    // Get points of the longest edge
    val edgePoint1 = triangle[longestEdge.first]
    val edgePoint2 = triangle[longestEdge.second]

    // Calculate vector from longest edge to opposite point
    val edgeVector = Offset(
        edgePoint2.x - edgePoint1.x,
        edgePoint2.y - edgePoint1.y
    )

    // Calculate midpoint of longest edge
    val midpoint = Offset(
        (edgePoint1.x + edgePoint2.x) / 2,
        (edgePoint1.y + edgePoint2.y) / 2
    )

    // Calculate vector from midpoint to opposite point
    val toOppositeVector = Offset(
        oppositePoint.x - midpoint.x,
        oppositePoint.y - midpoint.y
    )

    // Calculate the length of this vector
    val oppositeLength = sqrt(toOppositeVector.x * toOppositeVector.x + toOppositeVector.y * toOppositeVector.y)

    // Create a point on the opposite side with the same distance
    val fourthPoint = Offset(
        midpoint.x - (toOppositeVector.x / oppositeLength) * oppositeLength * 1.5f,
        midpoint.y - (toOppositeVector.y / oppositeLength) * oppositeLength * 1.5f
    )

    // Create result in correct order
    return sortPointsClockwise(
        listOf(
            edgePoint1,
            edgePoint2,
            oppositePoint,
            fourthPoint
        )
    )
}

private fun distance(p1: Offset, p2: Offset): Float {
    val dx = p2.x - p1.x
    val dy = p2.y - p1.y
    return sqrt(dx * dx + dy * dy)
}

private fun sortPointsClockwise(points: List<Offset>): List<Offset> {
    // Calculate centroid
    val centroid = Offset(
        points.sumOf { it.x.toDouble() }.toFloat() / points.size,
        points.sumOf { it.y.toDouble() }.toFloat() / points.size
    )

    // Sort points by their angle from centroid
    return points.sortedBy { point ->
        -kotlin.math.atan2(
            (point.y - centroid.y).toDouble(),
            (point.x - centroid.x).toDouble()
        )
    }
}

private fun expandLineToQuad(p1: Offset, p2: Offset): List<Offset> {
    val dx = p2.x - p1.x
    val dy = p2.y - p1.y
    val length = sqrt(dx * dx + dy * dy)
    val perpX = -dy / length * 10  // Perpendicular vector, scaled
    val perpY = dx / length * 10

    return listOf(
        p1,
        Offset(p1.x + perpX, p1.y + perpY),
        p2,
        Offset(p2.x + perpX, p2.y + perpY)
    )
}

private fun lineIntersectionBorder(p1: Offset, p2: Offset, p3: Offset, p4: Offset): Offset? {
    val denom = (p1.x - p2.x) * (p3.y - p4.y) - (p1.y - p2.y) * (p3.x - p4.x)
    if (denom.toDouble() == 0.0) return null

    val x = ((p1.x * p2.y - p1.y * p2.x) * (p3.x - p4.x) - (p1.x - p2.x) * (p3.x * p4.y - p3.y * p4.x)) / denom
    val y = ((p1.x * p2.y - p1.y * p2.x) * (p3.y - p4.y) - (p1.y - p2.y) * (p3.x * p4.y - p3.y * p4.x)) / denom

    return Offset(x, y)
}

private fun triangleAreaBorder(p1: Offset, p2: Offset, p3: Offset): Double {
    return abs((p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y)) / 2.0)
}}

Answer 2

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)

Answer 3

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.

Answer 4

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...

Answer 5

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.

Answer 6

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 :)

Answer 7

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.

Answer 8

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).

Answer 9

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.