多角形と点の内外判定 - 輪郭をなぞるだけのブログ

理解するのに少し時間がかかったので、次回以降に考えなくて済むようにメモ。
なので中身はあまりない。

ある点 p が多角形 polygon (凸とは限らない) の中にあるかどうかを判定。
Crossing Number Algorithm を使う。(下記 reference 参照)
~~ほぼパクリだが~~下の実装例では、点 p から左に半直線を伸ばしたときに、多角形の辺といくつ交点を持つかを調べている。
(交点の個数の偶奇だけが問題なので、衝突が検出されるたびに boolean 変数を反転している。)

交点を持つとき、どのような位置関係になるかを見てみる。
点 p を通る半直線の下側に点 a が、上側に点 b があるときを考える。
半直線と線分 ab が交点を持つかどうかは、下図のように点 p と点 a を通る直線の下側に点 b があるかで判定できる。
これは、外積を用いた正負判定で簡単に書ける。

また、 $eps$ の部分は浮動小数点数 $d$ が与えられたとき、 $|d| \leq eps$ を満たす $d$ を $0$ にみなすということである。(多分)

実装例

using ld = long double;
using Point = complex<ld>;
using Polygon = vector<Point>;
constexpr ld eps = 1e-9;

// 外積
ld cross(const Point& a, const Point& b){
    return a.real() * b.imag() - a.imag() * b.real();
}

// 内積
ld dot(const Point& a, const Point& b){
    return a.real() * b.real() + a.imag() * b.imag();
}

// CCW : Counter-ClockWise
enum class CCW_RESULT{
    CCW = +1, CW = -1, BEHIND = +2, FRONT = -2, ON = 0
};

// a, b を通る直線に対して、c がどの位置にあるか
CCW_RESULT ccw_enum(Point a, Point b, Point c){
    b -= a, c -= a;
    if(cross(b, c) > eps)   return CCW_RESULT::CCW; // 上側
    if(cross(b, c) < -eps)  return CCW_RESULT::CW;  // 下側
    
    // cross(b, c) == 0
    if(dot(b, c) < 0)   return CCW_RESULT::BEHIND;  // 同一直線上の左側
    if(norm(b) < norm(c))   return CCW_RESULT::FRONT;   // 同一直線上の右側
    return CCW_RESULT::ON;  // 線分 ab 上にある
}

int ccw(Point a, Point b, Point c){
    return static_cast<int>(ccw_enum(a, b, c));
}

enum class CONTAIN_RESULT {
    OUT = 0, ON = 1, IN = 2
};

// polygon の隣り合う頂点は辺をなす
// つまり、線分 polygon[i] - polygon[i+1] が辺である
CONTAIN_RESULT contain_point_enum(const Polygon& polygon, const Point& p){
    int n = polygon.size();

    // 点 p が多角形 polygon の辺上にあるか
    for(int i = 0; i < n; ++i){
        int place = ccw(polygon[i], polygon[(i + 1) % n], p);
        if(place == static_cast<int>(CCW_RESULT::ON)){
            return CONTAIN_RESULT::ON;
        }
    }

    // Crossing Number Algorithm
    bool in = false;
    for(int i = 0; i < n; ++i){
        Point a = polygon[i] - p;
        Point b = polygon[(i + 1) % n] - p;
        if(a.imag() > b.imag()) swap(a, b);
        if(a.imag() < eps && eps < b.imag() && cross(a, b) < -eps){
            in = !in;
        }
    }

    return in ? CONTAIN_RESULT::IN : CONTAIN_RESULT::OUT;
}

int contain_point(const Polygon& polygon, const Point& p){
    return static_cast<int>(contain_point_enum(polygon, p));
}