円板と多角形の共通部分の面積

問題. Circles - Intersection of a Circle and a Polygon

平面上の単純多角形 と円板 が与えられる．このとき， の面積を求めよ．

単純多角形とは自己交差していない多角形で， は の頂点を反時計回りに訪れた頂点の列 で与えられ， の座標を で表す．
また，円板 は中心 で半径 としたとき を満たす領域である．

制約： ，　，　

解法

まず多角形の面積の求め方を考える．多角形 の面積は原点 と の隣接する2頂点 からなる三角形の符号付き面積の総和で計算できる．三角形 の符号付き面積はベクトル の外積の絶対値に等しく となる．このときの符号は から へ反時計回りなら正，時計回りなら負となる．したがって，多角形 の面積は
　　
となる．

次に多角形と円板の共通部分 の面積の求め方を考える．まず始めに議論を簡単にするために，円板の中心が原点になるように円板と多角形が平行移動されていると仮定する．このとき平行移動によって面積は変わらない．
共通部分 の面積は上の多角形のみの場合と同様に各三角形 と円板 の共通部分の符号付き面積を求めてその総和をとることによって求めることができる．三角形 と円板 の共通部分は2点 と が円板の内外のどこにあるのかの位置関係によって下図のように場合分けされる．青色の領域は三角形の符号付き面積で，緑色の領域は扇形の符号付き面積である．
　　　　

扇形の半径をなす 2 つのベクトルを としたとき中心角 は という関係から atan2() 関数 を用いて で求まる．このとき， の符号は から へ反時計回りなら正で，時計回りなら負となるので上で述べた三角形の符号付き面積の符号と一致する．したがって，このときの扇形の符号付き面積は となる．
後は端点を と となる線分と円板の交点などを求めて上の図のように と の位置関係によって場合分けを行い三角形の符号付き面積や扇形の符号付き面積を求めて総和をとれば答えが求まる．

計算時間：

ソースコードを表示

#include <bits/stdc++.h>

using Number = double;

constexpr Number EPS = 1e-10;
const Number PI = acos(-1.0);

inline int sign(Number x) { return (x < -EPS) ? -1 : (x > EPS) ? +1 : 0; }
inline bool equal(Number a, Number b) { return sign(a - b) == 0; }
inline bool le(Number a, Number b) { return sign(a - b) == -1; } // a < b
inline bool ge(Number a, Number b) { return sign(a - b) == 1; } // a > b

struct Point {
    Number x = 0.0, y = 0.0;

    explicit Point() {}
    Point(Number x, Number y) : x(x), y(y) {}

    // Arithmetic operator between points
    Point operator+(const Point &rhs) const {
        return Point(this->x + rhs.x, this->y + rhs.y);
    }
    Point operator-(const Point &rhs) const {
        return Point(this->x - rhs.x, this->y - rhs.y);
    }
    Point operator*(const Point &rhs) const { // cross product between points
        return Point(this->x * rhs.x - this->y * rhs.y,
                     this->x * rhs.y + this->y * rhs.x);
    }

    Point operator-() const {
        return Point(-this->x, -this->y);
    }
    Point& operator+=(const Point &rhs) {
        return *this = *this + rhs;
    }
    Point& operator-=(const Point &rhs) {
        return *this = *this - rhs;
    }

    Point operator*(Number rhs) const {
        return Point(this->x * rhs, this->y * rhs);
    }
    Point operator/(Number rhs) const {
        return Point(this->x / rhs, this->y / rhs);
    }

    bool operator==(const Point &rhs) const {
        return equal(this->x, rhs.x) && equal(this->y, rhs.y);
    }
    bool operator!=(const Point &rhs) const {
        return !equal(this->x, rhs.x) || !equal(this->y, rhs.y);
    }
    bool operator<(const Point &rhs) const {
        return le(this->x, rhs.x) || (equal(this->x, rhs.x) && le(this->y, rhs.y));
    }
    bool operator>(const Point &rhs) const {
        return ge(this->x, rhs.x) || (equal(this->x, rhs.x) && ge(this->y, rhs.y));
    }

    Number abs(void) const {
        return sqrt(this->x * this->x + this->y * this->y);
    }
    Number abs2(void) const {
        return this->x * this->x + this->y * this->y;
    }
    Number arg(void) const {
        return atan2(this->y, this->x);
    }
    Number dot(const Point &rhs) const {
        return this->x * rhs.x + this->y * rhs.y;
    }
    Point rotate(double angle) const {
        return Point(cos(angle) * this->x - sin(angle) * this->y,
                     sin(angle) * this->x + cos(angle) * this->y);
    }
};

std::ostream& operator<<(std::ostream &os, const Point &p) {
    return os << p.x << ' ' << p.y;
}

std::istream& operator>>(std::istream &is, Point &p) {
    return is >> p.x >> p.y;
}

inline Number dot(const Point &p1, const Point &p2) {
    return p1.x * p2.x + p1.y * p2.y;
}

inline Number abs_cross(const Point &p1, const Point &p2) {
    return p1.x * p2.y - p1.y * p2.x;
}

inline Number arg(const Point &p1, const Point &p2) {
    return atan2(abs_cross(p1, p2), dot(p1, p2));
}

// Counter-Clockwise predicate (a, b, c)
enum CCW {
    COUNTER_CLOCKWISE = 1,     // counter clockwise
    CLOCKWISE         = -1,    // clockwise
    ONLINE_FRONT      = 2,     // a--b--c on line or (a == b and b != c)
    ONLINE_BACK       = -2,    // c--a--b on line
    ON_SEGMENT        = 0,     // a--c--b on line or (a != b and b == c) or (a == b == c)
    OTHER             = -3,
};

CCW ccw(const Point &a, Point b, Point c) {
    b -= a;  c -= a;
    if (sign(abs_cross(b, c)) == 1) return COUNTER_CLOCKWISE;
    if (sign(abs_cross(b, c)) == -1) return CLOCKWISE;
    if (sign(dot(b, c)) == -1)       return ONLINE_BACK;
    if (sign(b.abs2() - c.abs2()) == -1)   return ONLINE_FRONT;
    return ON_SEGMENT;
}

class Line : public std::array<Point, 2> {
public:
    explicit Line() {}
    Line(const Point &p1, const Point &p2) {
        (*this)[0] = p1;
        (*this)[1] = p2;
    }
};

inline CCW ccw(const Line &l, const Point &p) {
    return ccw(l[0], l[1], p);
}

class Segment : public Line {
public:
    explicit Segment() {}
    Segment(const Point &p1, const Point &p2) : Line(p1, p2) {}
};

class Circle : public Point {
public:
    Number r;
    explicit Circle() {}
    Circle(const Point &p, Number r = 0.0) : Point(p), r(r) {}
};

Point Projection(const Line &l, const Point &p) {
    Point dir = l[1] - l[0];
    Number t = dot(p - l[0], dir) / dir.abs2();
    return l[0] + dir * t;
}

inline bool IsIntersect(const Circle &c, const Point &p) {
    return (c - p).abs() <= c.r + EPS;
}

inline bool IsIntersect(const Circle &c, const Segment &s) {
    return IsIntersect(c, s[0]) || IsIntersect(c, s[1]) ||
        (IsIntersect(c, Projection(s, c))
         && ccw(s[0], Projection(s, c), s[1]) == CCW::ONLINE_FRONT);
}

std::vector<Point> CrossPoint(const Circle &c, const Segment &s) {
    if (!IsIntersect(c, s)) return std::vector<Point>();

    const Point mid = Projection(s, c), e = (s[1] - s[0]) / (s[1] - s[0]).abs();
    const Number len = sqrt(c.r * c.r - (mid - c).abs2());

    std::vector<Point> ps;
    const Point p1 = mid + e * len, p2 = mid - e * len;

    const CCW ccw1 = ccw(s[0], p1, s[1]), ccw2 = ccw(s[0], p2, s[1]);
    if (ccw1 == CCW::ONLINE_FRONT || p1 == s[1]) ps.emplace_back(p1);
    if (ccw2 == CCW::ONLINE_FRONT || p2 == s[1]) ps.emplace_back(p2);

    if (ps.size() == 2 && ccw(s[0], ps.back(), ps.front()) == CCW::ONLINE_FRONT)
        std::swap(ps.front(), ps.back());
    return ps;
}

class Polygon : public std::vector<Point> {
public:
    explicit Polygon() {}
    explicit Polygon(int size) : std::vector<Point>(size){}
};

// 円板 c と多角形 poly の共通部分の面積を返す
Number intersectionArea(const Circle &c, const Polygon &poly) {
    Number area = 0.0;

    const int n = poly.size();
    for (int i = 0; i < n; ++i) {
        const Point &p1 = poly[i] - c, &p2 = poly[(i + 1) % n] - c;

        // p1 と p2 が同一直線上にある場合は面積は 0 なのでスキップ
        if (abs(ccw(c, p1, p2)) != 1) continue;

        if ((p1.abs() < c.r - EPS) && (p2.abs() < c.r - EPS)) {
            area += 0.5 * abs_cross(p1, p2);
        }
        else if (p1.abs() < c.r - EPS) {
            const std::vector<Point> ps = CrossPoint(c, Segment(p1, p2));
            area += 0.5 * abs_cross(p1, ps.front());
            area += 0.5 * c.r * c.r * arg(ps.front(), p2);
        }
        else if (p2.abs() < c.r - EPS) {
            const std::vector<Point> ps = CrossPoint(c, Segment(p1, p2));
            area += 0.5 * c.r * c.r * arg(p1, ps.front());
            area += 0.5 * abs_cross(ps.front(), p2);
        }
        else {
            const std::vector<Point> ps = CrossPoint(c, Segment(p1, p2));
            if (ps.size() == 0) area += 0.5 * c.r * c.r * arg(p1, p2);
            else {
                area += 0.5 * c.r * c.r * arg(p1, ps.front());
                area += 0.5 * abs_cross(ps.front(), ps.back());
                area += 0.5 * c.r * c.r * arg(ps.back(), p2);
            }
        }
    }

    return area;
}


int main() {
    std::cout << std::fixed << std::setprecision(10);
    std::cin.tie(0); std::ios::sync_with_stdio(false);

    int n, r;
    std::cin >> n >> r;

    Circle c(Point(0, 0), r);
    Polygon poly(n);
    for (auto &p : poly) std::cin >> p;

    std::cout << intersectionArea(c, poly) << std::endl;

    return 0;
}

なかなか手強かった．
平面上の面積を測る道具に プラニメータ というのがあるのを最近知ったので面積繋がりでメモ．時間があれば原理を調べたい．