Official

N - 共通部分 / Intersection Editorial by KoD


以下では、ベクトル \(\bm{a} = (x_1, y_1), \bm{b} = (x_2, y_2)\) に対し \(c(\bm{a}, \bm{b}) = x_1y_2 - y_1x_2 \) と定めます。

\(S\) の頂点の位置ベクトルを \(\bm{p_i}\) とおきます。位置ベクトルが \(\bm{r}\) である点が \(S\) の周または内部に含まれることは、次が成り立つことと同値です。

  • 任意の \(1 \leq i \leq N\) に対し \(c(\bm{p_{i + 1} - p_i}, \bm{r} - \bm{p_i}) \geq 0\)(ただし、\(p_{N + 1} = p_1\) とする)

\(T\) の頂点の位置ベクトルを \(\bm{q_i}\) とおくと、\(S, T\) の共通部分は、\(S\) の周または内部のうち、\(c(\bm{q_{i + 1}} - \bm{q_i}, \bm{r} - \bm{q_i}) < 0\) である部分を取り除く処理を繰り返すことで求めることができます。この処理を行う方法を説明します。

\((c_i, d_i), (c_{i + 1}, d_{i + 1})\) を通る直線 \(l\)\(S\) の周の共有点は高々 \(2\) 個です。共有点が高々 \(1\) 個のとき、処理を行っても \(S\) は変化しないか、処理を行うことで \(S\) の面積が \(0\) になるかのいずれかです。共有点が \(2\) 個のとき、\(l\) によって \(S\) は二つに分けられ、そのうち一方が \(c(\bm{q_{i + 1}} - \bm{q_i}, \bm{r} - \bm{q_i}) \geq 0\) を満たし、もう一方は満たしません。\(c(\bm{q_{i + 1}} - \bm{q_i}, \bm{r} - \bm{q_i}) \geq 0\) を満たす部分の頂点および \((c_i, d_i), (c_{i + 1}, d_{i + 1})\) の和集合が処理後の凸多角形の頂点となります。

上の処理は \(O(N)\) で行うことができ、処理ごとに頂点は定数個しか増えないので、全体の計算量は \(O((N + M)M)\) となります。

実装例 (C++) :

#include <bits/stdc++.h>
using namespace std;

using F = double;
constexpr F EPS = 1e-10;

struct Point {
    F x, y;
    Point() : x(), y() {}
    Point(const F& x_, const F& y_) : x(x_), y(y_) {}
    Point operator-(const Point& other) const { return Point(x - other.x, y - other.y); }
    Point operator*(const F& other) const { return Point(x * other, y * other); }
    friend F cross(const Point& p, const Point& q) { return p.x * q.y - p.y * q.x; }
};

int main() {
    int n, m;
    cin >> n >> m;
    vector<Point> s(n);
    for (auto& [x, y] : s) {
        int a, b;
        cin >> a >> b;
        x = a, y = b;
    }
    vector<Point> t(m);
    for (auto& [x, y] : t) {
        int a, b;
        cin >> a >> b;
        x = a, y = b;
    }
    for (int i = 0; i < m; ++i) {
        const Point pivot = t[i];
        const Point vect = t[i + 1 == m ? 0 : i + 1] - t[i];
        const int len = s.size();
        vector<Point> next;
        for (int j = 0; j < len; ++j) {
            const Point p = s[j], q = s[j + 1 == len ? 0 : j + 1];
            const F p_c = cross(vect, p - pivot), q_c = cross(vect, q - pivot);
            if (p_c + EPS > 0) {
                next.push_back(p);
            }
            if (((p_c > EPS) and (q_c < -EPS)) or ((p_c < -EPS) and (q_c > EPS))) {
                next.push_back(p - (q - p) * (p_c / (q_c - p_c)));
            }
        }
        s = move(next);
    }
    const int len = s.size();
    F ans = 0;
    for (int i = 0; i < len; ++i) {
        ans = ans + cross(s[i], s[i + 1 == len ? 0 : i + 1]);
    }
    cout << fixed << setprecision(20);
    cout << ans / 2 << '\n';
}

posted:
last update: