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

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.

Answer 2

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 3

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