Cutting a convex polygon by a line

概要

凸多角形 \(P = (v_1, v_2, ..., v_N)\) を 点\(p_0, p_1\) で切断し、その左側にできる凸多角形 \(P'\) を求める

計算量

\(O(N)\)

実装

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

def cross3(O, A, B):
    ox, oy = O; ax, ay = A; bx, by = B
    return (ax - ox) * (by - oy) - (bx - ox) * (ay - oy)
def cross_point(p0, p1, q0, q1):
    x0, y0 = p0; x1, y1 = p1
    x2, y2 = q0; x3, y3 = q1
    dx0 = x1 - x0; dy0 = y1 - y0
    dx1 = x3 - x2; dy1 = y3 - y2
    s = (y0-y2)*dx1 - (x0-x2)*dy1
    sm = dx0*dy1 - dy0*dx1
    if -EPS < sm < EPS:
        return None
    return x0 + s*dx0/sm, y0 + s*dy0/sm
EPS = 1e-9
def convex_cut(P, line):
    q0, q1 = line
    N = len(P)
    Q = []
    for i in range(N):
        p0 = P[i-1]; p1 = P[i]
        cv0 = cross3(q0, q1, p0)
        cv1 = cross3(q0, q1, p1)
        if cv0 * cv1 < EPS:
            v = cross_point(q0, q1, p0, p1)
            if v is not None:
                Q.append(v)
        if cv1 > -EPS:
            Q.append(p1)
    return Q

P = [(0, 0), (10, 1), (1, 10)]
line = ((5, -3), (7, 8))
for x, y in convex_cut(P, line):
    print("%.4f %.4f" % (x, y))
# =>
# 0.0000 0.0000
# 5.6481 0.5648
# 6.3846 4.6154
# 1.0000 10.0000

Verified

  • AOJ: "CGL_4_C: Convex Cut": source (Python3, 0.03sec)

  • AOJ: "2160: Voronoi Island": source (Python3, 0.04sec)