最小面积四边形算法

Question

有一些算法可以找到包含给定（凸）多边形的最小边界矩形。

有人知道找到最小面积边界四边形（任何四边形，而不仅仅是矩形）的算法吗？

我已经在互联网上搜索了几个小时，但是虽然我找到了一些关于这个问题的理论论文，但我没有找到一个实现......

编辑：Mathoverflow 的人向我指出了一篇带有数学解决方案的文章 (my post there)，但我没有找到实际的实现。我决定采用卡尔的蒙特卡洛方法，但我会在有时间的时候深入研究论文并在这里报告......

谢谢大家！

Answer 1

小气算法

从你的凸多边形开始。虽然点数超过 4 个，但找到合并成本最低的一对相邻点，然后合并它们，将多边形中的点数减少 1。

我所说的“合并”是指将两侧的两条边延伸到某一点，然后用该点替换这两条边。

“最便宜”是指合并为多边形增加最小面积的一对。

Before:                       After consolidating P and Q:

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

在 O(n log n) 中运行。但这只会产生一个近似值，我对此并不完全满意。算法在生成漂亮的正五边形方面做得越好，最后一次合并必须吃掉的区域就越多。

Answer 2

蒙特卡洛方法

感谢 cmets 对问题的澄清。我认为需要的不是数学上正确的结果，而是比其他形状的任何可比拟合更好的“拟合”。

与其在这个问题上投入大量的算法脑力，不如让计算机来处理它。生成 4 个随机点的组；检查通过凸连接这 4 个点形成的四边形是否与多边形相交，并计算四边形的面积。重复 100 万次，检索面积最小的四边形。

您可以应用一些约束来使您的点不是完全随机的；这可以显着提高收敛性。

---

蒙特卡洛，改进

我一直坚信，在飞机上随机扔 4 个点是一个非常低效的开始，即使对于蛮力解决方案也是如此。因此，以下细化：

• 对于每个试验，随机选择 p 个不同的顶点和 q 个多边形的不同边，这样 p + q em> = 4。
• 对于每一个 q 边，构造一条穿过该边端点的线。
• 对于每个 p 个顶点，构造一条穿过该顶点并随机分配斜率的线。
• 验证 4 条线确实形成了一个四边形，并且这个四边形包含（并且不相交！）多边形。如果这些测试失败，请不要继续进行此迭代。
• 如果这个四边形的面积是目前看到的所有面积中的最小值，请记住该四边形顶点的面积和坐标。
• 重复任意次数，并返回找到的“最佳”四边形。

与总是需要 8 个随机数（4 个点中的每一个的 x 和 y 坐标）相反，此解决方案只需要 (4 + p) 个随机数。此外，产生的线条并不是在平面上盲目地挣扎，而是每条都接触到多边形。这确保了四边形从一开始就至少非常接近多边形。

Answer 3

我认为围绕点的 2D OBB 是一个很好的起点。这可能会给出一个“好的”（但不是最好的）适合；如果您发现仍然需要更严格的界限，可以将拟合扩展到四边形。

首先，您应该计算输入点的凸包。这使您处理的分数更少，并且根本不会改变答案。

我会远离 Gottschalk 在其他地方的论文中引用的基于协方差的方法。这并不总是最适合，并且对于非常简单的输入（例如正方形）可能会出错。

在 2D 中，http://cgm.cs.mcgill.ca/~orm/maer.html 中描述的“旋转卡尺”方法应该给出最合适的 OBB。这个问题也很容易被认为是一维最小化问题：

1. 对于给定的角度，将 x 轴和 y 轴旋转该角度。这为您提供了两个正交向量。
2. 对于每个点，投影到两个轴上。跟踪每个轴上的最小和最大投影。
3. OBB 的面积为 (max1-min1)*(max2-min2)，您可以通过角度以及最小和最大投影轻松计算 OBB 的点。
4. 重复，对区间 [0, pi/2] 进行尽可能精细的采样，并跟踪最佳角度。

John Ratcliffe 的博客 (http://codesuppository.blogspot.com/2006/06/best-fit-oriented-bounding-box.html) 有一些代码可以在 3D 中执行这种采样方法；如果遇到困难，您应该能够使其适应您的 2D 案例。 警告：John 有时会在他的博客文章中发布轻微的 NSFW 图片，但该特定链接很好。

如果您在完成这项工作后仍然对结果不满意，您可以稍微修改一下方法；我能想到两个改进：

1. 您可以选择两个非正交轴，而不是像上面提到的那样选择正交轴。这会给你最合适的菱形。为此，您实际上需要在 [0, pi] x [0, pi] 上进行 2D 搜索。
2. 使用最适合的 OBB 作为更详细搜索的起点。例如，您可以从 OBB 中取出 4 个点，移动其中一个直到线接触到一个点，然后重复其他点。

希望对您有所帮助。抱歉信息超载:)

Answer 4

我认为定向边界框应该非常接近（尽管它实际上是一个矩形）。这是关于定向边界框的标准参考论文：Gottschalk's paper (PDF)

查看第 3 部分（安装 OBB）。

Answer 5

我相信最小面积边界四边形的每一边都将通过多边形的至少 2 个顶点。如果这个假设是正确的，那么对解决方案执行蛮力搜索应该很容易。

• 生成由 2 个顶点定义且不与多边形相交的唯一线集。
• 检查每组 4 条线，以确定哪条线产生最小面积边界四边形。

更新：边界四边形的每一边将通过至少 2 个顶点的假设是错误的，但一组相关的线可以为解决方案提供基础。至少，边界四边形的每一边都将通过一个顶点，该顶点用于定义由 2 个顶点定义且不与多边形相交的唯一线集。

Answer 6

以下是对蒙特卡洛算法的改进，也可能会导致直接答案。

引理：如果最优四边形的一条边仅与多边形的一个点接触，则它与该边的中点接触。

证明：将 X 和 Y 定义为接触点两侧的两条线段的长度。想象一下将边缘围绕单个触摸点旋转一个无穷小的角度 A。旋转会通过将四边形增加 XA/2 和将其减少 YA/2 来影响四边形的大小，反之亦然。如果 X != Y，则两个旋转方向之一导致较小的四边形。由于四边形是最小的，我们必须有 X=Y。

要使用这个事实，请注意，如果我们选择多边形与四边形接触的一些边和点，并且我们不连续选择两个点，我们可以通过选择通过每个点的边来唯一确定四边形使该点成为边缘的中点（如果不可能，则拒绝选择的配置）。在蒙特卡洛算法中，这简化为说我们不必为这条边选择斜率 - 它可以明确地确定。

我们仍然有两个相邻点被选中的案例 - 希望我已经启发了其他人来这里......

Answer 7

我有同样的问题要解决，我使用的代码实际上是用一个额外的矩形来实现 Jason Orendorff 的想法，它限制了过程并使结果更像方形。最后，这是一个很好的启发式方法，在我的情况下效果很好。我希望其他人也可以从这段代码中受益：

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;
    }
}

如果过程无法找到解决方案，可以通过例如从小的边界矩形开始并顺序增加它来进一步改进算法。 在实践中，我使用的边界矩形比最小边界矩形大 5%。