Day 038-Q2 | Sutherland-Hodgman 凸多角形クリッピング

問題

凸多角形 $P$（頂点 $N$ 個、反時計回り）を凸多角形 $C$（頂点 $M$ 個、反時計回り）でクリッピングした $P \cap C$ の面積を求めよ。Sutherland-Hodgman アルゴリズムを使って実装すること。

制約

パラメータ	範囲
$N, M$	$3 \le N, M \le 10^4$
座標	$-10^9 \le x, y \le 10^9$
誤差許容	絶対誤差 $10^{-6}$
制限	時間 3 sec / メモリ 256MB

入出力例

入力例 1

出力例 1

4.000000

正方形 [0,4]² を [1,3]² でクリップ → 面積 2×2 = 4.0

概念図: Sutherland-Hodgman クリッピング

$C$ の各辺が定義する半平面で $P$ を順番にカット。各辺処理で頂点リストが更新される。

ヒント（段階的開示）

ヒント1: 方向性

Sutherland-Hodgman はクリップ多角形 $C$ の各辺について順番に対象多角形をカットする。各辺は「左側（内側）」の半平面を定義し、外側の点を除いて辺との交点を追加する。

ヒント2: アプローチ

$C$ の各辺 $(E_1, E_2)$ について、対象多角形の各辺 $(A, B)$ を処理:
1. A 内・B 内 → B を出力
2. A 内・B 外 → 交点を出力
3. A 外・B 内 → 交点と B を出力
4. A 外・B 外 → 何も出力しない

ヒント3: コード骨格

def clip_by_edge(polygon, a, b):
    result = []
    n = len(polygon)
    for i in range(n):
        A = polygon[i]
        B = polygon[(i+1) % n]
        iA = is_inside(A, a, b)
        iB = is_inside(B, a, b)
        if iA:
            result.append(A)  # A は出力
            if not iB:
                result.append(line_intersect(A, B, a, b))
        else:
            if iB:
                result.append(line_intersect(A, B, a, b))
    return result

# 内側判定: cross(a, b, p) >= -EPS
# 面積: shoelace / 2

模範解答 (Python)

import sys
input = sys.stdin.readline

EPS = 1e-9

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

def is_inside(p, a, b):
    return cross(a, b, p) >= -EPS

def line_intersect(p1, p2, p3, p4):
    d1 = (p2[0]-p1[0], p2[1]-p1[1])
    d2 = (p4[0]-p3[0], p4[1]-p3[1])
    denom = d1[0]*d2[1] - d1[1]*d2[0]
    if abs(denom) < EPS:
        return p1
    t = ((p3[0]-p1[0])*d2[1] - (p3[1]-p1[1])*d2[0]) / denom
    return (p1[0] + t*d1[0], p1[1] + t*d1[1])

def clip_by_edge(polygon, a, b):
    if not polygon:
        return []
    result = []
    n = len(polygon)
    for i in range(n):
        A = polygon[i]
        B = polygon[(i+1) % n]
        iA = is_inside(A, a, b)
        iB = is_inside(B, a, b)
        if iA:
            result.append(A)
            if not iB:
                result.append(line_intersect(A, B, a, b))
        else:
            if iB:
                result.append(line_intersect(A, B, a, b))
    return result

def polygon_area(poly):
    n = len(poly)
    s = 0.0
    for i in range(n):
        x1, y1 = poly[i]
        x2, y2 = poly[(i+1) % n]
        s += x1*y2 - x2*y1
    return abs(s) / 2.0

def sutherland_hodgman(subject, clip):
    output = subject[:]
    n = len(clip)
    for i in range(n):
        if not output:
            return []
        a = clip[i]
        b = clip[(i+1) % n]
        output = clip_by_edge(output, a, b)
    return output

def read_polygon():
    n = int(input())
    pts = []
    for _ in range(n):
        x, y = map(float, input().split())
        pts.append((x, y))
    return pts

def solve():
    P = read_polygon()
    C = read_polygon()
    result = sutherland_hodgman(P, C)
    print(f"{polygon_area(result):.6f}")

solve()

Step-by-Step 解説

1内外判定
辺 $(A, B)$ の左側に点 $P$ があるかを外積 $\text{cross}(A, B, P) \ge 0$ で判定。反時計回りの多角形では「左 = 内側」。EPS で境界上の点も内側扱い。

2直線交点計算
パラメトリック表現 $P_1 + t(P_2 - P_1)$ で $t$ を求め交点座標を計算。分母が EPS 未満なら平行として $p_1$ を返す（安全な fallback）。

3clip_by_edge
各辺 $(A, B)$ について 4 ケースを判別し、次の多角形の頂点列を生成。$A$ を出力する（先頭を B の前ループで追加する仕組み）。

4全辺でループ
$C$ の $M$ 辺について `clip_by_edge` を繰り返す。各辺後に頂点数は最大 $+1$ 増える。全体 $O(NM)$。

5靴ひも公式（Shoelace）
$\text{Area} = \frac{1}{2}|\sum_i (x_i y_{i+1} - x_{i+1} y_i)|$ で面積を算出。

計算量

クリッピング: $O(NM)$
面積計算: $O(N+M)$
合計: $O(NM)$ — $N, M \le 10^4$ なら $10^8$ でギリギリ

よくあるミス

ミス	原因	正しい書き方
辺の向きと内側定義がずれる	CW/CCW を混同	反時計回りでは cross > 0 が内側
平行辺でゼロ除算	denom ≈ 0 未処理	`if abs(denom) < EPS: return p1`
空多角形チェック漏れ	中間で空になる	各 clip_by_edge 前後で `if not output: return []`
浮動小数点で境界上を外側判定	厳密な >= 0 のため	`>= -EPS` で許容

次のステップ

発展問題: 凸多角形 $K$ 個の共通部分の面積 $O(K^2 N)$
応用: 半平面交差（Half-plane Intersection）による凸領域の高速構築 $O(N \log N)$
最適化: 辺の角度ソート + 線形スキャンで $O(N + M)$ クリッピング

自己評価

完全理解おおむね理解要復習

Day 038-Q2 — Sutherland-Hodgman 凸多角形クリッピング

問題

制約