0%

凸包问题

发表于 更新于
from Peter Norvig's pytudes project

把一束钉子钉在一块木板上，然后用橡皮筋把它们围起来，让橡皮筋绷紧（taut)，就像这样：

image

橡皮筋已经把钉子集合的凸包画出来了。这是计算机图形学、机器人运动规划、地理信息系统、动物行为学等领域应用中的一个重要问题。更正式地说，我们说：

给定一个有限集, *P, 对于平面上的所有点, **P 的凸包是多边形 H, 满足:

  • P 中的每一点都在 H内或 H上.
  • H 上的每个顶点都是 P 中的一个点.
  • H 是凸的: 连接H的任意两个顶点的线段 ，要么是H的边，要么位于H内部.

在此notebook中，我们开发了一个算法来寻找凸包（并展示了如何使用matplotlib绘图的例子）。

首先，我们要弄清楚如何表示我们感兴趣的对象：

  • 点
  • 点集
  • 多边形：我们用顶点的有序列表来表示多边形

第一步：完成必要的import:

1
2
3
4
5
6
7
from __future__ import division, print_function

%matplotlib inline 
import matplotlib.pyplot as plt
import collections
import random
import math

点和点集

我将类Point定义为x和y坐标的named tuple，Points(n)定义为创建一组n个随机点的函数。

对于函数Points(n)有两个难点：

​ 1.第二个可选参数用来生成随机数种子

​ 2.由于matplotlib默认绘制在3×2矩形上，因此将从3×2矩形中均匀采样点（每条边上有一个0.05的小边框，以防止点碰到矩形的边）。

1
2
3
4
5
6
7
Point = collections.namedtuple('Point', 'x, y')

def Points(n, seed=42):
    "Generate n random points within a 3 x 2 box."
    random.seed((n, seed))
    b = 0.05
    return {Point(random.uniform(b, 3-b), random(b, 2-b)) for _ in range(n)}

可视化点和线段

现在来看看怎么可视化点；定义一个函数plot_points.我们想看到的是：

  • 点自身
  • 可选的，点之间的线段。一个可选式的参数允许你确定是不是要这个线，以及要什么颜色的。该参数使用由matplotlib定义的标准样式格式。例如”r”。表示没有线的红色圆点，“bs-”表示有线条的蓝色正方形，“go”表示点之间有虚线的绿色圆圈。线从点到点依次排列；如果要使线从最后一点向后一点闭合（以形成一个完整的多边形），请指定closed = True。 （为此，点的集合必须是一个列表； closed = False时，该集合可以是任何集合。）
    可选地，在这些点上的标签可以使我们彼此区分。 如果指定labels = True，则会得到标签（从0到n的整数）。
1
2
3
4
5
6
7
8
9
10
11
def plot_points(points, style='r.', labels=False, closed=False):
    """Plot a collection of points, Optionally change the line style, labe points with numbers,
    and /or for a closed polygon by closing the line from the last point to the first.
    """
    if labels:
        for (i, (x, y)) in enumerate(points):
            plt.text(x, y, ' '+str(i))
    if closed:
        points = points + [points[0]]
    plt.plot([p.x for p in points], [p.y for p in points], style, linewidth=2.5)
    plt.axis('scaled'); plt.axis('off')

凸性(Convexity)

我们想构造一个凸包，所以最好有一些确定多边形是否凸的方法。 让我们检查一个是：

1
2
3
octagon = [Point(-10, 0), Point(-7, -7), Point(0, -10), Point(+7, -7), 
           Point(+10, 0), Point(+7, +7), Point(0, +10), Point(-7, 7)]
plot_points(octagon, 'bs-', labels=True, closed=True)

img

如果你从左侧的点0开始，并沿八角形逆时针依次移动，从点到点跟随边缘，则可以看到在每个顶点处都在向左转。

现在让我们考虑一个非凸多边形：

img

吃豆人（pacman）多边形是非凸的； 你会看到一条从点3到点5的线经过了多边形。 你还可以看到，当你将逆时针方向从3移到4到5时，你将在4处右转。这导致了这个想法：如果我们在逆时针旋转多边形时没有向右转，则该多边形为凸形。

转向

现在我们如何确定从点A到B到C的转弯在B处左转还是右转（或直行）？ 考虑下图：

img

如果角度β大于角度α，则在B处向左转； 换句话说，如果β的tan值大于α：

​ $(C.y - B.y) / (C.x - B.x) > (B.y - A.y) / (B.x - A.x)$

但是，如果执行该计算，则每个分母为零时将需要特殊情况。 因此，将每边乘以分母：

​ $(B.x - A.x) * (C.y - B.y) > (B.y - A.y) * (C.x - B.x) $

（注意：这一步会让你非常紧张！通常，将不等式的两边都乘以负数会逆转不等式，在这里分母可能是负数。在这种情况下，它是可行的；基本上是因为我们正在做两次乘法 因此负数可以抵消，但是数学证明是棘手的，涉及向量代数中的一些概念，因此我在这里不再重复；取而代之，我将在下面提供良好的测试范围。）

1
2
3
4
def turn(A, B, C):
    "Is the turn from A->B->C a 'right', 'left', 'straight' turn?"
    diff =  (B.x - A.x) * (C.y - B.y)  -  (B.y - A.y) * (C.x - B.x) 
    return ('right' if diff < 0 else 'left' if diff > 0 else 'straight')

凸包算法概述

现在 我们有了寻找凸包策略的第一部分：

沿着点以某种顺序走一条路径(还不确定以什么顺序) 这条路上的任何没有标记为左转的点都不是凸包的一部分

什么是好的顺序？让我们看看，如果我们从最左边开始，一直走到最右边会发生什么。我们可以通过调用在点上排序的内置函数来实现这种排序（因为点是元组，排序按字典顺序排序：首先由它们的第一个X坐标，如果有关系，下一个是它们的Y坐标）。我们从11个随机点开始，我将定义一个函数来帮助绘制部分外壳(hull)：

1
2
3
def plot_partial_hull(points, hull_indexes=()):
    plot_points(points, labels=True)
    plot_points([points[i] for i in hull_indexes], 'bs-')
1
plot_partial_hull(sorted(Points(11)))

img

现在，我将按照从点0到1到2到3的顺序开始构建外壳：

1
plot_partial_hull(sorted(Points(11)), [0, 1, 9, 10])

img

1
plot_partial_hull(sorted(Points(11)), [10, 8, 4, 0])

ig

1
plot_partial_hull(sorted(Points(11)), [0, 1, 9, 10] + [10, 8, 4, 0])

img

算法的基本思想就是这样，但有一些边界情况需要担心：

  • 退化多边形：如果只有1/2(甚至0)个点会发生什么？这样的一组点应该被认为是凸的，因为没有办法画出这些线段外的点。

  • 共线(colinear)点：如果三个或更多的点是共线的，我们应该只保留两个“外部”点。不保留它们的理由是我们希望凸包是最小可能的点集。我们需要保留外面的，因为它们在外壳上标出了真正的角。当它是一个“straight”转向和当它是一个“right”转向时，我们可以通过拒绝一个点来实现。

  • 第一个和最后一个点：聪明(astute)的读者可能已经注意到，在A->B->C的转弯中，我们的算法只拒绝中间点B。这意味着，排序顺序中的第一个和最后一个点永远不会成为拒绝的候选点，因此总是会在外壳上。是这样吗？是的。第一个点是最左边的点，x值最低的点（如果有一样的，它是最左边的最低点）。这是一个极端的角落，所以它应该一直在外壳上。最后一个点同理。

凸包算法的实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def convex_hull(points):
    "Find the convex hull of a set of points"
    if len(points) <= 3:
        return points
    upper = half_hull(sorted(points))
    lower = half_hull(reversed(sorted(points)))
	return upper + lower[1:-1]

def half_hull(sorted_points):
    hull =[]
    for C in sorted_points:
        while len(hull) >= 2 and turn(hull[-2], hull[-1], C) != 'left':
            hull.pop()
        hull.append(C)
    return hull

可视化结果

1
2
3
4
5
def plot_convex_hull(points):
    hull = convex_hull(points)
    plot_points(points)
    plot_points(hull, 'bs-', closed=True)
    print(len(hull), 'of', len(points), 'points on hull')

n == 10,000的时候会不会很慢？

1
2
P10K = Points(10000)
%timeit convex_hull(P10K)

104 ms ± 3.69 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

1
plot_convex_hull(P10K)

img

26 of 10000 points on hull

还有一些有意思的例子再[这里][[https://github.com/Qasak/pytudes/blob/master/ipynb/Convex%20Hull.ipynb][https://github.com/Qasak/pytudes/blob/master/ipynb/Convex%Hull.ipynb]

测试

到目前为止，一切看起来都很好！但如果我们能通过一套测试，我会更加自信：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def tests():
    # Tests of `turn`
    assert turn(octagon[0], octagon[1], octagon[2]) == 'left'
    assert turn(octagon[2], octagon[3], octagon[4]) == 'left'
    assert turn(octagon[1], octagon[0], octagon[7]) == 'right'
    assert turn(octagon[5], octagon[6], octagon[7]) == 'left'
    assert turn(octagon[2], octagon[1], octagon[0]) == 'right'
    assert turn(pacman[1], pacman[2], pacman[3]) == 'left'
    assert turn(pacman[3], pacman[4], pacman[5]) == 'right'
    assert turn(Point(0, 0), Point(0, 1), Point(0, 2)) == 'straight'
    assert turn(Point(2, 1), Point(3, 1), Point(4, 1)) == 'straight'
    assert turn(Point(2, 1), Point(4, 1), Point(3, 1)) == 'straight'
    assert turn(Point(0, 0), Point(1, 1), Point(2, 2)) == 'straight'
    assert turn(Point(0, 0), Point(-1, -1), Point(2, 2)) == 'straight'
    # More tests of `turn`, covering negative denominator
    A, B = Point(-2, -2), Point(0, 0)
    assert turn(A, B, Point(1, 3)) == 'left'
    assert turn(A, B, Point(2, 2)) == 'straight'
    assert turn(A, B, Point(3, 1)) == 'right'
    assert turn(A, B, Point(-1, 1)) == 'left'
    assert turn(A, B, Point(-1, -4)) == 'right'
    assert turn(A, B, Point(-1, -1)) == 'straight'
    assert turn(B, A, Point(-3, -4)) == 'left'
    assert turn(B, A, Point(-4, -3)) == 'right'
    assert turn(B, A, Point(-1, -1)) == 'straight'
    assert turn(B, A, Point(-3, -3)) == 'straight'
    
    # Tests of convex_hull
    assert convex_hull(octagon)== octagon
    assert convex_hull(circle) == convex_hull(donut)
    assert convex_hull(circle) == convex_hull(convex_hull(circle))
    for n in (0, 1, 2, 3):
        assert convex_hull(Points(n)) == Points(n)
    collinear = {Point(x, 0) for x in range(100)}
    assert convex_hull(collinear) == [min(collinear), max(collinear)]
    assert convex_hull(collinear | {P}) == [min(collinear), max(collinear), P]
    P = Point(5, 5)
    grid1 = {Point(x, y) for x in range(10) for y in range(10)}
    assert convex_hull(grid1) == [Point(0, 0), Point(9, 0), Point(9, 9), Point(0, 9)]

    return 'tests pass'
        
tests()

外壳上有多少点？

Points（N）在外壳上的点数似乎随着N的增加而缓慢增加。有多慢？让我们试试看。我们将对60次随机试验中的Points（N）取外壳上的点的平均数：

1
2
3
4
5
def average_hull_size(N, trials=60):
    """Compute the average hull size of N random points
    (averaged over the given number of random trials)."""
    return sum(len(convex_hull(Points(N, seed=trials+i)))
               for i in range(trials)) / trials

我们将对N的几个值进行此操作，取2的幂(We’ll do this for several values of N, taken as powers of 2:)：

1
2
3
4
5
6
hull_sizes = [average_hull_size(2**e) 
              for e in range(14)]

print('   N   Hull Size')
for e in range(14): 
    print('{:4}: {:4.1f}'.format(2**e, hull_sizes[e]))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
   N   Hull Size
   1:  1.0
   2:  2.0
   4:  3.7
   8:  5.1
  16:  7.1
  32:  8.6
  64: 11.0
 128: 12.6
 256: 14.6
 512: 16.4
1024: 18.1
2048: 19.8
4096: 21.6
8192: 23.2
1
2
3
4
5
6
def plot_hull_sizes(hull_sizes):
    plt.plot(hull_sizes, 'bo-')
    plt.ylabel('Hull size')
    plt.xlabel('log_2(number of points)')

plot_hull_sizes(hull_sizes)

img

这意味着我们可以通过计算直线的斜率（slope）和截距(intercept)来定义estimated_hull_size。（我不需要做线性回归；我只需要画一条从hull_sizes的第一点到最后一点的直线。）

1
2
3
4
5
6
7
def estimated_hull_size(N):
    slope = (hull_sizes[-1] - hull_sizes[0]) / (len(hull_sizes) - 1)
    return hull_sizes[0] + slope * math.log(N, 2)

plt.plot([estimated_hull_size(2**e) 
          for e in range(len(hull_sizes))], 
         'r--');

结束语和引伸阅读

凸壳问题是算法设计中的一个有趣的问题。这里的算法叫做安德鲁单调链 Andrew’s Monotone Chain. 它是格雷姆扫描的变种 Graham Scan.在这里阅读更多 Tamassia or Wikipedia.