tangent lines from a point to a convex polygon

概要

ある点 \(w\) から凸多角形 \(P = (v_1, v_2, ..., v_N)\) への接線を引いた時に、その接線を通る凸多角形の点 \(v_i\) を求める。

この時、点 \(w\) は凸多角形 \(P\) の外側に存在するとする。

具体的計算

凸多角形 \(P\) 内の適当な点 \(s\) を選ぶ。

\(h(i)\) を 点 \(w\) と \(s\) を通る直線と点 \(w\) と \(v_i\) を通る直線がなす角 (\(h(i) \in [-\pi, \pi]\)) と定義する。

関数 \(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\) を求める。

この求めた \(x_0\) と \(x_1\) について、点 \(w\) と \(x_0\) を通る直線 と 点 \(w\) と \(x_1\) を通る直線が 点 \(w\) から 凸多角形 \(P\) への接線になる。

計算量

\(O(\log N)\)

実装

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

from math import atan2

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 convex_polygon_tangent(p0, qs):
    L = len(qs)
    d = L//3
    gx = (qs[0][0] + qs[d][0] + qs[2*d][0]) / 3
    gy = (qs[0][1] + qs[d][1] + qs[2*d][1]) / 3

    x0, y0 = p0
    dx = gx - x0; dy = gy - y0
    sv = - x0*dy + y0*dx; cv = - x0*dx - y0*dy
    h = lambda p: atan2(p[0]*dy - p[1]*dx + sv, p[0]*dx + p[1]*dy + cv)

    i0 = i1 = -1

    EPS = 1e-12

    v0 = h(qs[0])
    v1 = h(qs[L-1])
    if abs(v0 - v1) < EPS:
        v2 = h(qs[1])
        # assert v2 != v0
        if v0 < v2:
            i0 = L-1
        else:
            i1 = L-1
    else:
        v2 = h(qs[1])
        if v2 >= v0 < v1:
            i0 = 0
        elif v2 <= v0 > v1:
            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[x+B]) - h(qs[x+B+1]), L)
        i1 = (i0 + k) % L
    else:
        B = i1-L
        k = binary_search(lambda x: h(qs[x+B+1]) - h(qs[x+B]), L)
        i0 = (i1 + k) % L
    # assert i0 != -1 != i1
    return i0, i1

参考

  • 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.

戻る