Point in polygon

Point in polygon

어떠한 점이 다각형 내부에 있는지, 외부에 있는지 판정하는 알고리즘을 알아볼 것이다.

CCW(Counter Clock Wise)

CCW는 벡터가 회전하는 방향을 판정하는 알고리즘으로 벡터의 외적을 이용하면 쉽게 구할 수 있다.

세 점 $a, b, c$가 있을 때, 이 점들을 이용해 벡터 $\vec{ab}, \vec{bc}$를 만든다. 이 두 벡터의 외적값을 이용해 세 점에 대한 방향성을 판정한다.

우선 3차원에서의 점 $\vec{p} = (p_x, p_y, p_z), \vec{q} = (q_x, q_y, q_z)$가 있을 때 두 벡터의 외적 $\vec{p} \times \vec{q}$는 $\vec{p} \times \vec{q} = (p_yq_z - p_zq_y, p_xq_z - p_zq_x, p_xq_y - p_yq_x)$로 나타낼 수 있고, $\vec{p} \times \vec{q}$의 방향은 오른손 법칙을 따른다.

2차원에서의 벡터 외적은 3차원에서의 백터 외적에서 $z$좌표만 0으로 두면 된다. 두 벡터 $\vec{p}, \vec{q}$에 대해 $\vec{p} = (p_x, p_y, 0), \vec{q} = (q_x, q_y, 0)$로 나타낼 수 있고, $\vec{p} \times \vec{q} = (0, 0, p_xq_y - p_yq_x)$가 성립한다.

오른손 법칙에 따라 $\vec{p}, \vec{q}$에 대해 $p_xq_y - p_yq_x$가 양수이면, $\vec{q}$는 $\vec{p}$로 부터 반시계 방향에 있고, 반대로 음수이면 $\vec{q}$는 $\vec{p}$로 부터 시계방향에 있다.

세 개의 점 $a, b, c$가 주어졌을 때 벡터의 방향성 판단을 정리하면 다음과 같다.

1. $\vec{ab} \times \vec{bc} > 0$
   $\vec{ab}$에 비교했을 때, $\vec{bc}$의 방향은 반시계 방향이다. (CCW)
2. $\vec{ab} \times \vec{bc} < 0$
   $\vec{ab}$에 비교했을 때, $\vec{bc}$의 방향은 시계 방향이다. (CW)
3. $\vec{ab} \times \vec{bc} = 0$
   점 $a, b, c$는 일직선 상에 위치한다.

CCW에서 정방향은 반시계이기 때문에, 결과가 반시계일 경우 1, 시계는 -1, 직선은 0으로 표현한다.

c++구현은 다음과 같다.

#define X first
#define Y second

typedef long long ll;
typedef pair<ll, ll> pll;

int ccw(pll v1, pll v2){
    
    ll cross = v1.X * v2.Y - v2.X * v1.Y;
    
    if(cross > 0) return 1;
    if(cross < 0) return -1;
    return 0;
}

CCW는 여러 용도로 사용이 가능한데 그 중 하나가 두 선분의 교차 판정이다. 두 선분에 대한 방정식을 세워 판정할 수도 있지만, CCW를 이용해 더 간단하게 풀 수 있다.

일반적으로는 위와 같이 선분 ab와 cd가 있을 때, 벡터 ab와 bc의 ccw, 벡터 ab와 bd의 ccw 곱이 0이하라면(즉, 둘의 ccw가 다르거나 0이 있다면) 두 선분이 교차한다고 할 수 있다.

단, 몇가지 예외가 있는데,

위와 같은 경우 제대로 판정이 되지 않는다. 그래서 추가적으로 벡터 cd와 벡터 da의 ccw, 벡터 cd와 벡터 db의 ccw의 곱 또한 0이하인지 확인해야 한다.

ccw의 곱 2개가 모두 0인 경우도 선분이 만나지 않을 수 있다. 이 경우 두 선분의 겹치는 부분이 존재하는지 추가적으로 확인해야 한다.

ccw를 이용한 선분 교차 판정 구현은 다음과 같다.

pll p2v(pll p1, pll p2){
    return {p2.X - p1.X, p2.Y - p1.Y};
}

bool cross(pll a, pll b, pll c, pll d){
    int cc1 = ccw(p2v(a, b), p2v(b, d)) * ccw(p2v(a, b), p2v(b, c));
    int cc2 = ccw(p2v(c, d), p2v(d, a)) * ccw(p2v(c, d), p2v(d, b));
    
    if(!cc1 && !cc2){
        if(a > b) swap(a, b);
        if(c > d) swap(c, d);
        if(a > c){
            swap(a, c);
            swap(b, d);
        }
        return c <= b;
    }
    return cc1 <= 0 && cc2 <= 0;
}

Point in Concave polygon

내부 점 판정에 앞서 오목 다각형과 볼록 다각형의 차이에 대해 알아보면, 볼록 다각형은 내각이 전부 180도 이하인 다각형이고, 오목 다각형은 그렇지 않은 다각형이다.

오목 다각형의 문제는 다음과 같다.

N개의 점으로 이루어진 다각형과 M개의 점이 주어진다. 각각의 점들이 다각형의 내부에 있는지 판별하라. (단, 다각형의 경계 또한 내부이다.)

문제를 해결하는 방법은 판정할 점을 한 쪽 끝으로 반직선을 그어서(실제로 구현 시 판정할 점과 x좌표가 매우 큰 점을 잇는다.) 그 반직선과 다각형의 교점 개수를 확인한다.

교점이 홀수라면 다각형의 내부, 짝수라면 다각형의 외부에 존재한다.

N각형의 변은 N개 이므로 판정해야될 점이 M개라면 각 변에 대하여 선분 교차판정을 진행하여 $O(NM)$의 시간에 해결이 가능하다.

몇 가지 주의할 경우가 있는데 다음과 같다.

1. 판정할 점이 다각형의 경계 위에 있다면 교점의 개수와 관계없이 판정 결과가 정해진다. (경계도 내부로 치면 내부 아니면 외부)
2. 반직선이 다각형의 변과 일치하면(평행한데 만나면) 교차한 것으로 카운트하지 않는다.
3. 반직선이 다각형의 꼭짓점을 지나간다면(두 변이 만나는 점을 지난다면) 반직선 위에 있는 변만 교차한 것으로 카운트 한다.

2, 3번이 조금 번거로운데 쉽게 처리하는 방법으로 반직선을 구현할 때, 판정할 점의 좌표를 (x, y)라 하면 선분의 다른 끝점 좌표를 (INT_MAX, y+1)로 하여 약간의 기울기를 주는 것이다. 이렇게 하면 2, 3번이 모두 발생하지 않는다.

구현은 다음과 같다.

bool on_line(pll a, pll b, pll p){
    if(a > b) swap(a, b);
    return !ccw(p2v(a, b), p2v(b, p)) && a<=p && p<=b;
}

bool Point_In_Concave(vector<pll> C, pll P){
    int cnt = 0;
    for(int i=0; i<C.size();i++){
        if(on_line(C[i], C[(i+1)%C.size()], P)) return true;
        if(cross(C[i], C[(i+1)%C.size()], P, {INT_MAX, P.Y+1})) cnt++;
    }
    return cnt%2;
}

Point in Covex polygon

볼록 다각형 내부 점 판정 문제는 다음과 같다.

N개의 점으로 이루어진 볼록 다각형과 M개의 점이 주어진다. 각각의 점들이 다각형의 내부에 있는지 판별해라. (단, 다각형의 경계 또한 내부이다.)

볼록 다각형도 오목 다각형과 마찬가지로 내부 점 판정이 가능하지만, 좀 더 빠르게 판정할 수 있는 방법이 있다.

볼록 다각형의 한 점에서 시작해, 나머지 점들로 이어지는 반직선들을 긋는다. $n$각형에서 어떤 점이 다각형 내부에 위치하기 위해서는 1번과 $n-1$번 반직선 사이에 점이 존재해야 한다. 두 직선 사이에 점이 존재함을 확인했다면, 반직선으로 만들어지는 $n-2$개의 구간 중 어떤 구간에 점이 존재하는지 이분탐색으로 확인한다.

모든 비교는 CCW를 통해 할 수 있으며, 점이 $i$번 반직선과 $i+1$번 반직선 사이에 있음을 확인했으면, 세 점의 방향성을 판단하여 역시 CCW로 다각형 내부에 있는지 판정한다.

시간 복잡도는 다각형이 이미 만들어져 있다는 가정하에 $O(log \ n)$이고, 구현은 다음과 같다.

int bin_search(vector<pll>& P, int L, int R, pll point){
    if(L+1 == R) return L;
    int M = (L + R) / 2;
    if(ccw(p2v(P[0], P[M]), p2v(P[0], point)) > 0){
        return bin_search(P, M, R, point);
    } 
    else{
        return bin_search(P, L, M, point);
    }
}

bool Point_In_Convex(vector<pll>& P, pll S){
    if(ccw(p2v(P[0], P[1]), p2v(P[0], S))<0 
    || ccw(p2v(P[0], P[P.size()-1]), p2v(P[0], S))>0){ 
        return 0;
    } 
    int idx = bin_search(P, 1, P.size()-1, S);
    if(ccw(p2v(P[idx], S), p2v(S, P[idx+1]))>0) return 0;
    else return 1;
}

Reference

1. 공군 휴머니스트 air-wiki
2. https://anz1217.tistory.com/107