Intersection between a convex polygon and a line

概要

凸多角形 \(P = (v_1, v_2, ..., v_N)\) と 点 \(p_0, p_1\) を通る直線の交差判定と交点を求める

具体的計算

\(h(i) = (v_{(i-1 \text{ mod } N)+1} - p_0) \times (p_1 - p_0)\) と定義する。

この関数は、頂点が直線の左右どちらかにあるかによって符号が変化し、直線と頂点の距離に比例する。

関数 \(h\) の最大値をとる \(x\) を \(x_1\), 最小値をとる \(x\) を \(x_0\) とする時、 この関数は以下の性質を満たす

  • 区間 \([x_0+1, x_0+N-1]\) について 関数 \(h\) は上に凸な関数

  • 区間 \([x_1+1, x_1+N-1]\) について 関数 \(h\) は下に凸な関数

この最大値と最小値は2回の三分探索で求める。 1回目の三分探索では 区間 \([1, N]\) で定義される以下の関数 \(g\) 上で最小値をとる \(x_0\) を求める。(\(h(1) = h(N)\) の場合や \(h(1)\) が極値の場合は特殊ケースとして処理)

\(\displaystyle g(x) = \min\left(h(x), \frac{h(N) - h(1)}{N-1}(x-1) + h(1)\right)\)

2回目の三分探索では 区間 \([x_0+1, x_0+N-1]\) 上の関数 \(h\) で \(x_1\) を求める。

この最大値、最小値を求めた後、交点を含む辺を見つけるために \(h(i)\) の符号が変化する箇所を二分探索で求める。

計算量

\(O(\log N)\)

実装

凸多角形の頂点は反時計回りのリストとして持つとする。

def find_zero(x0, x1, f):
    v0 = f(x0)
    if v0 == 0:
        return x0+1
    left = x0; right = x1+1
    while left+1 < right:
        mid = (left + right) >> 1
        if v0 * f(mid) >= 0:
            left = mid
        else:
            right = mid
    return right

def binary_search(f, L):
    left = 0; right = L
    while left+1 < right:
        mid = (left + right) >> 1
        if f(mid) < 0:
            left = mid
        else:
            right = mid
    return right

def line_polygon_intersection(p0, p1, qs):
    x0, y0 = p0; x1, y1 = p1
    dx = x1 - x0; dy = y1 - y0
    h = lambda p: (p[0] - x0)*dy - (p[1] - y0)*dx
    L = len(qs)

    i0 = i1 = -1

    v0 = h(qs[0])
    v1 = h(qs[L-1])
    if v0 == v1:
        v2 = h(qs[1])
        # assert v0 != v2
        if v0 < v2:
            i0 = L-1
        else:
            i1 = L-1
    else:
        v2 = h(qs[1])
        if v1 > v0 <= v2:
            i0 = 0
        elif v1 < v0 >= v2:
            i1 = 0
        else:
            g = lambda x: min((v1 - v0)*x/(L-1) + v0, h(qs[x]))
            i0 = binary_search(lambda x: g(x+1) - g(x), L-1)

    if i1 == -1:
        B = i0 - L
        k = binary_search(lambda x: h(qs[B+x]) - h(qs[B+x+1]), L)
        i1 = (i0 + k) % L
    else:
        B = i1 - L
        k = binary_search(lambda x: h(qs[B+x+1]) - h(qs[B+x]), L)
        i0 = (i1 + k) % L

    if h(qs[i0]) * h(qs[i1]) > 0:
        # a line and a polygon are disjoint
        return []

    # a vertex to the left side of a line: i0
    # a vertex to the right side of a line: i1

    f = lambda i: h(qs[i-L])
    k0 = find_zero(i1, i0 if i1 < i0 else i0+L, f) % L
    k1 = find_zero(i0, i1 if i0 < i1 else i1+L, f) % L
    # vertices to the left side of a line: k0, k0+1, ..., k1-2, k1-1
    # vertices to the right side of a line: k1, k1+1, ..., k0-2, k0-1
    if k0 == k1:
        return [k0]
    return [k0, k1]

# e_i is a line segment between v_{i-1} and v_i

p0, p1 = (0, 0), (1, 1)
qs = [(0, 2), (2, 0), (4, 0), (4, 2), (2, 4), (0, 4)]
print(line_polygon_intersection(p0, p1, qs))
# => "[4, 1]": cross points is on either e_4 or e_1

p0, p1 = (0, -1), (4, 0)
qs = [(0, 2), (2, 0), (4, 0), (4, 2), (2, 4), (0, 4)]
# => "[3]": cross point is on e_3

参考

  • Chazelle, Bernard, and David P. Dobkin. "Intersection of convex objects in two and three dimensions." Journal of the ACM (JACM) 34.1 (1987): 1-27.

戻る