找到形成凸多边形的最大点子集

Posted 2023-02-14

【中文标题】找到形成凸多边形的最大点子集

【英文标题】：Finding largest subset of points forming a convex polygon

【发布时间】：2014-03-13 18:19:48

【问题描述】：

我正在寻找一种算法来查找从给定点集形成凸多边形的最大点子集（我的意思是数量上的最大子集）。 我认为这可能可以使用 DP 解决，但我不确定。 是否可以在 O(n^3) 中做到这一点？ 其实我只需要最大子集的大小，所以它不需要有唯一的解决方案

编辑：

只是为了保持简单，

给定输入： 一组二维点

期望输出：形成凸多边形的最大点数，如示例中的输出为 5（ABHCD 是可能的凸多边形之一）

有一个类似的问题 spoj.com/problems/MPOLY 可以使用 O(N^3) 中的 DP 解决，我的问题是关于该问题的概括，它不必包含 (0,0)

【问题讨论】：

那不是一个圆圈吗？或者说，最大的点子集给定一组点？

我稍微编辑了我的问题，也许它会帮助理解我的问题

找到所有可能的凸多边形是一种选择吗？

就像找到所有可能的组合，然后在列表中搜索最大的。

我正在寻找一个多项式解决方案，因为 N 可能高达 250，不过感谢您的想法

【参考方案1】：

我想我通过扩展Andrew's algorithm for convex hulls 想出了一个算法。

最初，我想出了一个 O(N^4) 算法，但注意到它过于复杂并将其归结为 O(N^2) 算法。似乎代码中的某个地方可能存在一个错误，至少在垂直点的情况下会导致问题。这可能是一种特殊情况（需要更改几行代码），或者是算法存在较大缺陷的迹象。如果是后者，那么我倾向于说它是 NP 难的，并将算法作为启发式提供。我把所有的时间都花在了编码和调试上。无论如何，它在其他情况下似乎工作正常。但是，当正确的算法似乎是 O(2^N) 时，很难测试其正确性。

def maximal_convex_hull(points):
    # Expects a list of 2D points in the format (x,y)


    points = sorted(set(points)) # Remove duplicates and sort
    if len(points) <= 1: # End early
        return points

    def cross(o, a, b): # Cross product
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])


    # Use a queue to extend Andrew's algorithm in the following ways:
    # 1. Start a new convex hull for each successive point
    # 2. Keep track of the largest convex hull along the way
    Q = []
    size = 0
    maximal_hull = None
    for i,p in enumerate(points):
        remaining = len(points) - i + 1
        new_Q = []
        for lower, upper in Q: # Go through the queue
            # Build upper and lower hull similtanously, slightly less efficent
            while len(lower) >= 2 and cross(lower[-2], lower[-1], p) < 0:
                lower.pop()
            lower.append(p)
            while len(upper) >= 2 and cross(upper[-2], upper[-1], p) > 0:
                upper.pop()
            upper.append(p)

            new_size = len(lower) + len(upper)
            if new_size > size: # Check for a new highest maximal convex hull
                size = new_size
                maximal_hull = (lower, upper)


        # There is an odd bug? that when one or both if statements are removed
        #  xqg237.tsp (TSPLIB) gives slightly different solutions and is
        #   missing a number of points it should have.
        #  Seems like if the opposite should be true if anything since they
        #   were intended to be easy optimizations not required code.
            if remaining + new_size >= size: # Don't bother if not enough
                new_Q.append((lower, upper)) # Add an updated convex hulls
        if remaining > size: # Don't bother if not enough
            new_Q.append(([p], [p])) # Add a new convex hull
        Q = new_Q # Update to the new queue

    # Extract and merge the lower and upper to a maximium convex hull
    # Only one of the last points is ommited because it is repeated
    #    Asserts and related variables used for testing
    #    Could just have "return lower[:-1] + list(reversed(upper[:-1]))[:-1]"
    lower, upper = maximal_hull
    max_hull_set = set(lower) | set(lower) | set(upper)
    lower = lower
    upper = list(reversed(upper[:-1]))[:-1]
    max_convex_hull = lower + upper
    assert len(lower) + len(upper) == len(max_hull_set)
    assert max_hull_set == set(max_convex_hull)
    return max_convex_hull

编辑：此算法不正确，但它为正确算法提供了灵感，因此我将其留在这里。

【讨论】：

谢谢！我的想法是，通过扩展安德鲁的算法，我得到了一个使用动态编程的 O(N^4) 算法，我认为它没有任何缺陷，我仍然对你的 O(N^2) 有一些疑问，但我周末再去看看:)

Andrew 的算法是 O(N)（但需要排序，使其成为 O(N log N)），该算法运行 O(N) 版本的安德鲁算法。为 O(N^2) 运行时做准备。虽然，算法可能还远远不够。

是的，但是我对那个算法的正确性有些怀疑，你能解释一下这个算法是如何工作的，因为我对 phython 不是很熟悉

就像这张图片中的imgur.com/DsJhF71，在插入C时没有弹出EFG，尽管最佳的上壳是AEGFD

该示例的重点是什么？您可能认为应该分别考虑上下船体。虽然这仍然会导致 O(N^2) 算法，只是具有更大的隐藏系数。

【参考方案2】：

这个问题已经在这些比赛中发表过：

SPOJ problem BOYSCOUT USACO DEC08 problem "Largest Fence"（页面上的最后一个问题）

其解决方案（O(N³) 算法）可以在此页面上找到："USACO DEC08 Problem 'fence' Analysis" Bruce Merry 和 Jacob Steinhardt。

下面是解释这个算法的尝试。还有here is my implementation 这个算法在 C++11 中。此代码适用于 NPython script，它生成满足这些要求的测试数据。

O(N²) 算法简化问题

简化问题：找到形成凸多边形并包含给定集合最左边的点的最大子集（或者如果有多个最左边的点，则该多边形应该包含其中的最上面的点）。

为了解决这个问题，我们可以通过有向线段连接每对点，然后（隐式）将每个线段视为一个图形节点，如图所示：

这里的父节点和对应的段是蓝色的。与其左子（红色）对应的线段从“父”线段的末端开始，它是相对于“父”线段方向左转最少的线段。其右孩子（绿色）对应的线段从“父”线段的起点开始，也是相对于“父”线段方向左转最少的线段。

任何在最左边点结束的段都对应一个“叶”节点，即它没有子节点。这种图的构造保证了没有循环，换句话说这个图是一个有向无环图。

现在要找到点的最大凸子集，只需找到该图中的最长路径即可。 DAG 中最长的路径可以通过动态规划算法在 O(E) 时间内找到，其中 E 是图中的边数。由于每个节点只有2条出边，每条对应一对点，这个问题可以在O(N²)时间内解决。

如果我们预处理输入数据，所有这些都可以实现，按角度对从同一点开始的有向线段进行排序。这允许在图中执行深度优先遍历。我们应该记住从每个处理过的节点开始的最长路径，这样每个图的边最多被访问一次，我们有 O(E) = O(N²) 时间算法（忽略预处理时间现在）。空间要求也是 O(N²)，因为我们需要为每对点存储排序后的方向，并且每个节点使用一个值（这也是一对点）。

这里是这个动态规划算法的 C++ 实现：

unsigned dpStep(ind_t i_left, ind_t from, ind_t to_ind)

    ind_t child = sorted[from][to_ind];

    while (child < i_left)
        child = sorted[from][++to_ind];

    if (child == i_left)
        return 1;

    if (auto v = memorize[from][child])
        return v;

    const auto x = max(dpStep(i_left, child, lefts[from][child]) + 1,
                       dpStep(i_left, from, static_cast<ind_t>(to_ind + 1)));

    memorize[from][child] = static_cast<ind_t>(x);
    return x;

输入参数是最左边的点和形成当前线段的一对点（其中线段的起点直接给出，而终点作为按角度数组排序的索引给出）。此代码片段中的“left”一词被稍微过度使用：它表示最左边的点 (i_left) 和包含图形左孩子的预处理数组 (lefts[][])。

O(N³) 算法

一般问题（最左边的点不固定）可以这样解决：

对左右方向的点进行排序 对点进行预处理以获得两个数组：(a) 每个点按相对于其他点的方向排序，(b) 线段末端的第一个数组中的位置，使得相对于“父级”方向的左转最少" 段。 将每个点作为最左边的点，通过“简化”算法找到最长的凸多边形。 定期从预处理数组中删除“最左侧”点左侧的所有点。 采用第 3 步中找到的最长路径。

第 4 步是可选的。它不会提高 O(N³) 时间复杂度。但它通过排除不需要的点略微提高了 DP 算法的速度。在我的测试中，这提供了 33% 的速度提升。

预处理有多种选择。我们可以计算每个角度（使用atan2）并对它们进行排序，正如 Bruce Merry 和 Jacob Steinhardt 在示例代码中所做的那样。这是一种最简单的方法，但atan2 可能过于昂贵，尤其是对于较小的点集。

或者我们可以排除 atan2 并排序切线而不是角度，就像在我的实现中所做的那样。这有点复杂但速度更快。

这两种选择都需要 O(N² log N) 时间（除非我们使用一些分布排序）。第三种选择是使用众所周知的计算几何方法（排列和对偶）来获得 O(N²) 预处理。但我认为它对我们的问题没有用：由于常数因子大，对于较小的点集可能太慢了，对于较大的点集，它可能会提高一些速度，但太微不足道了，因为第 3 步将大大超过第 2 步。而且实现起来要困难得多。

其他一些结果：我尝试实现迭代 DP 而不是递归；这并没有明显改变速度。我还尝试一次执行两个 DP 搜索，将每个搜索的步骤交错（类似于在软件中实现的光纤或超线程）；这可能会有所帮助，因为 DP 主要由具有高延迟和阻止 CPU 吞吐量的充分利用的内存访问组成；结果不是很令人印象深刻，只有大约 11% 的速度提升，很可能使用真正的超线程会快得多。

【讨论】：

很抱歉打扰，但我无法理解给定 URL cerberus.delos.com:790/TESTDATA/DEC08.fence.htm 上的一件事它被写为“船体中的第一个和最后两个点足以检查下一个点是否'不打破凸性条件'。你能解释一下吗？为什么不需要所有积分？点是否按特定顺序排列？

@Naman：看起来很清楚。我只能试着用不同的语言来解释它。在每个 DP 步骤中，有 3 组点：（a）提到的 4 个点（第一对/最后一对），（b）已经在船体中的点（它们已经被 DP 算法选择，因此它们满足凸性条件并且它们是最优的), (c) 所有其他点。在每一步中，算法独立地尝试组 (c) 中的每个点，并检查该点相对于组 (a) 的点的凸性条件。如果点适合，则不需要检查其相对于（b）组点的凸性，保证满足。

不知何故我仍然无法理解为什么它是有保证的。我能解释一个案例吗？考虑第一对（2,2），（3,1），最后一对（8,2），（7,1）和船体中已经存在的点（6,6），（7,5）。现在一个新的点 (5,3) 来了。它将从第一对和最后一对成功凸性条件，但当我们针对包括 (b) 组中的所有点考虑它时，它不会成功。我知道我在理解上犯了一些小错误，但如果您能纠正我，我将不胜感激。

@Naman：现在我明白了你问题的最后一部分。是的，积分已经安排好了。在您的示例中，该对中的第一个点（也是序列中的第一个）是（3,1），最后一个点是（7,1）。所以 (5,3) 不能插入 (7,1) 之后。如果你这样做，你会得到 (8,2), (7,1), (5,3), (3,1), (2,2) 这不是凸的。

所以你的意思是在给定算法之前我需要对给定点进行排序/排序？如果是这样订购如何（按x，y或顺时针？）。如果我问愚蠢的问题，我很抱歉。如果你能帮我提供一些有详细解释的链接，那也很有帮助。

【参考方案3】：

这是我的一个动态规划 O(N^4) 算法，其想法来自 Nucleman 发布的 Andrew's Algorithm，我仍在寻找 O(N^3) 算法

upper_hull(most left point, previous point, current point)

    ret = 0
    if (current point != most left point)   /* at anytime we can decide to end the upper hull and build lower_hull from current point ending at most left point */
        ret = min(ret,lower_hull(most left point, -1, current point)) 
    for all point on the right of current point /* here we try all possible next point that still satisfy the condition of convex polygon */
        if (cross(previous point,current point,next point) >= 0) max(ret,1+upper_hull(most left point, current point, next point))
    return ret;


lower_hull(most left point, previous point, current point)

    if (current point == most left point) return 0;
    ret = -INF /* it might be impossible to build a convex hull here, so if it does it will return -infinity */
    for all point on the left of current point and not on the left of most left point
        if (cross(previous point,current point,next point) >= 0) max(ret,1+lower_hull(most left point, current point, next point))
    return ret;

先按x轴排序，然后按y轴排序，然后尝试所有点作为最左边的点运行upper_hull(p,-1,p)，请告诉我这个算法是否有任何缺陷

【讨论】：

一种可能性是预先计算所有的叉积结果 (O(N^3)) 并根据结果是正数还是负数将它们分成两个图表，然后使用深度优先搜索开始在最左边的点上找到最长的最短路径。似乎它应该在 O(E) 中是可行的，因为边缘是三元组 (a,b),(b,c)，所以是 O(N^3)。但是，您必须对 O(N) 个点（对于每个最左边的点）执行此操作，从而导致总体时间复杂度为 O(N^4)。因此，我现在相当确定你不会比 O(N^4) 更好。

感谢您的赏金，很高兴能为您提供帮助。

【参考方案4】：

您可以使用 delaunay 三角剖分并删除最长的边和顶点。我使用类似的算法来找到凹壳。您可以找到基于人口数据@http://www.phpdevpad.de/geofence 的 jan 示例。我还写了一个php类[e-mail]。