기하 알고리즘 - Convex Hull, Rotating Calipers

Convex Hull

 2차원 평면 위의 점들을 모두 포함하는, 넓이가 최소인 볼록 다각형을 구해보자. 이는 볼록 껍질이라는 것을 구하면 해결된다.

 

 볼록 껍질을 구하는 것은 여러 방법이 있지만, 나는 반껍질 두 개를 구하는 방식을 선호한다. 한 쪽 반껍질을 구하고, 다른 반껍질은 반대로 구하면 되기 때문이다. (저번 게시글에서 볼록 껍질을 구할 때 각도 정렬이 필요하다고 했지만..) 이 방법은 각도 정렬을 하지 않아서 편하다.

 

 일단 반껍질을 구하기 전, x좌표, y좌표 순으로 증가하도록 점들을 정렬하자. 점들을 정렬했다면, 이제 볼록 껍질을 만들 수 있다. 점들을 스택에 넣으며, 마지막 세 점의 방향이 반시계 방향이 되도록 적절히 스택에서 pop을 하면 된다. 글로 쓰는 것은 내가 봐도 설명이 구리기 때문에, 아래 그림과 함께 설명하도록 하겠다.

처음은 그냥 스택에 넣는다. 초록색 점이 스택에 들어있는 점들이다.

스택의 마지막 세 점의 ccw가 양수가 아니기 때문에, 마지막에서 두 번째 점은 스택에서 제거한다.

마지막 세 점의 ccw가 양수이기 때문에, 스택에 넣어준다.

똑같이, 마지막 세 점의 ccw가 양수이기 때문에 스택에 넣어준다.

스택의 마지막 세 점의 ccw가 양수가 아니기 때문에, 마지막에서 두 번째 점은 스택에서 제거한다.

 

마찬가지로 마지막 세 점의 ccw가 양수가 아니기 때문에, 마지막에서 두 번째 점은 스택에서 제거한다.

 힘겨운 과정이 모두 끝났다. 스택에 남은 초록색 점들이 반껍질을 이루는 점들이다. 반대편 반껍질은 ccw가 음수인 경우에만 스택에 추가해주는 방식으로 해주면 된다. 점들의 순서가 반시계 방향인 것에 대한 이점이 많이 때문에, 반대편 반껍질은 순서를 뒤집어서 합쳐주면 된다.

 주의해야 할 점은, 각 반껍질의 양 끝 점이 겹치기 때문에 이를 빼주어야 한다. 다음은 이를 구현한 코드이다.

#include <bits/stdc++.h>
 
#define X first
#define Y second
#define endl '\n'
using namespace std;
 
typedef long long ll;
typedef pair<ll, ll> pll;
 
ll ccw(pll A, pll B, pll C) {
    return A.X * B.Y + B.X * C.Y + C.X * A.Y - A.Y * B.X - B.Y * C.X - C.Y * A.X;
}
 
void convexHull(vector<pll>& v){
    sort(v.begin(), v.end());
    vector<pll> L, R;
    for (pll i: v){
        while (L.size() >= 2 && ccw(L[L.size() - 2], L[L.size() - 1], i) <= 0) L.pop_back();
        L.emplace_back(i);
        while (R.size() >= 2 && ccw(R[R.size() - 2], R[R.size() - 1],i) >= 0) R.pop_back();
        R.emplace_back(i);
    }
    v.clear();
    for (int i=0;i<L.size()-1;i++){
        v.emplace_back(L[i]);
    }
    for (int i=R.size()-1;i>0;i--){
        v.emplace_back(R[i]);
    }
}
 
int main() {
    ios_base::sync_with_stdio(false);cin.tie(nullptr);
 
    int n; cin >> n;
    vector<pll> inp;
    for (int i=0;i<n;i++){
        ll x, y; cin >> x >> y;
        inp.emplace_back(x, y);
    }
    convexHull(inp);
    cout << inp.size() << endl;
    return 0;
}

 

 

BOJ 1708 - 볼록 껍질

 위 코드를 그대로 사용하면 된다.

 

BOJ 4181 - Convex Hull

 세 점이 한 직선 위에 있어도 상관 없기 때문에, \(ccw == 0\) 인 경우에도 스택에 넣어주면 된다. 위 코드에서 ccw의 부등호만 수정해주면 된다.

 

BOJ 7420 - 맹독 방벽

 방벽의 길이를 최소화하기 위해서는, 점들의 볼록 껍질에서 길이가 \(L\)만큼 떨어진 곳에 방벽을 세워야 한다. 잘 생각해보면, (볼록 껍질의 둘레의 길이) + (반지름이 \(L\)인 원의 둘레의 길이) 가 답임을 알 수 있다. 다음은 이 문제의 정답 코드이다.

#include <bits/stdc++.h>
 
#define X first
#define Y second
#define endl '\n'
using namespace std;
 
typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> pll;
 
ld dist(pll A, pll B){
    ld dx = A.X - B.X;
    ld dy = A.Y - B.Y;
    return sqrt(dx * dx + dy * dy);
}
 
ll ccw(pll A, pll B, pll C) {
    return A.X * B.Y + B.X * C.Y + C.X * A.Y - A.Y * B.X - B.Y * C.X - C.Y * A.X;
}
 
void convexHull(vector<pll>& v){
    sort(v.begin(), v.end());
    vector<pll> L, R;
    for (pll i: v){
        while (L.size() >= 2 && ccw(i, L[L.size() - 1], L[L.size() - 2]) <= 0) L.pop_back();
        L.emplace_back(i);
        while (R.size() >= 2 && ccw(i, R[R.size() - 1], R[R.size() - 2]) >= 0) R.pop_back();
        R.emplace_back(i);
    }
    v.clear();
    for (int i=0;i<R.size()-1;i++){
        v.emplace_back(R[i]);
    }
    for (int i=L.size()-1;i>0;i--){
        v.emplace_back(L[i]);
    }
}
 
int main() {
    ios_base::sync_with_stdio(false);cin.tie(nullptr);
 
    int n, l; cin >> n >> l;
    vector<pll> inp;
    for (int i=0;i<n;i++){
        ll x, y; cin >> x >> y;
        inp.emplace_back(x, y);
    }
    convexHull(inp);
    inp.emplace_back(inp.front());
    ld ans = acos(-1) * 2 * l;
    for (int i=0;i<inp.size() - 1;i++){
        ans += dist(inp[i], inp[i + 1]);
    }
    cout << round(ans) << endl;
 
    return 0;
}

 

 

 

Rotating Calipers

캘리퍼스

 2차원 평면 위의 점들이 있다. 여기서 가장 먼 두 점의 거리를 구하려면 어떻게 해야 할까?

 가장 쉬운 방법은 당연하게도 \(O(N^2)\)로 구하는 것이다. 하지만 Rotating Calipers를 사용하면, \(O(N)\)에 구할 수 있다. 아이디어는 투 포인터를 기반으로 하고 있으니, 이를 알고 있다면 이해에 도움이 될 수 있다. (아마도)

 

 점들의 볼록 껍질을 생각해보자. 볼록 껍질 내부의 점은 거리가 가장 먼 두 점에 속할 수 없다.

 볼록 껍질 내부의 점 \(P\)가, 거리가 가장 먼 두 점 \(A, B\) 중 하나(편의 상 \(B=P\))라고 가정하자. 점 \(A\)를 중심으로 하고, 선분 \(\overline{\rm AP}\) 를 반지름으로 하는 원 밖에 볼록 껍질 위의 점은 항상 존재한다. 이는 \(\overline{\rm AP} < \overline{\rm AQ}\) 를 만족하는 점 \(Q\)가 항상 존재한다는 뜻이고, 앞서 한 가정과 모순된다.

 

 볼록 껍질 위의 점만을 확인해도 된다는 것을 알아냈다. 이제 한 점에서, 가장 먼 점을 찾을 것이다. 실제 캘리퍼스처럼, 평행선을 움직이며 찾아보자.

 

 먼저, 한 벡터를 고정하고, 이 벡터와 평행한 다른 벡터의 시점을 반시계방향으로 돌리며 찾는다. 앞에서, 볼록 껍질을 반시계 방향으로 구한 이유가 이것이다.

고정한 벡터는 빨간색으로 표현했다.

 맨 앞에서 언급한 투 포인터의 방식대로, 두 평행선 사이의 거리가 줄어들지 않는 한, 파란 벡터의 시점, 즉 비교할 벡터의 시점을 계속 옮겨주자. 이 점을 언제까지 돌려야 하는지는, ccw를 이용하여 확인할 수 있다. 파란 벡터의 시점을 \(v_i\)라고 하자. 고정한 벡터와, 시점이 \(v_i\)이고, 종점이 \(v_{i+1}\)인 벡터의 방향성을 확인해주면 된다. 

초록색 벡터의 시점을, 빨간색 벡터의 종점으로 옮긴다면?

 초록색 벡터(비교할 벡터)의 시점을, 빨간색 벡터(고정한 벡터)로 옮기자. 그리고 ccw를 구했을 때 이 값이 음수라면, 즉 시계 방향이라면 비교 벡터의 시점을 더 이상 옮겨줄 필요가 없어진다.

C를 B로 옮겨주어야 한다.

 이때 점 \(D\)는, \(D - C + B\)로 옮겨줄 수 있다. 이제 점 \(A, B, D\)의 방향성을 확인하기만 하면 된다.

두 점 사이의 거리를 구해준다.

 이제 두 점 사이의 거리를 구하고, 고정 벡터를 옮겨주자. 투 포인터로 진행하기 때문에, 파란 벡터의 시점은 계속하여 들고 가자.

 귀찮은 관계로 여기까지만 그리겠다.. 대충 이런 식으로 한 바퀴를 돌면 된다. 고정한 벡터와 비교할 벡터 모두 \(O(N)\)에 돌아간다. 다음은 이 과정을 구현한 코드이고, 코드를 보며 이해하는 것이 빠를 것이라고 생각한다.

pll operator + (const pll A, const pll B){
    return {A.X + B.X, A.Y + B.Y};
}
 
pll operator - (const pll A, const pll B){
    return {A.X - B.X, A.Y - B.Y};
}
 
int r = 0;
ll ans = 0;
for (int i = 0; i < n; i++){
    while (r < n * 2 && ccw(v[i], v[(i + 1) % n], v[(i + 1) % n] + v[(r + 1) % n] - v[r % n]) >= 0){
        ans = max(ans, dist(v[i], v[r % n]));
        r++;
    }
    ans = max(ans, dist(v[i], v[r % n]));
}

 

BOJ 2049 - 가장 먼 두 점

 주어진 점들의 볼록 껍질을 구하고, 여기에서 Rotating Calipers를 사용하여 가장 먼 두 점의 거리를 구하면 된다. 다음은 이 문제의 정답 코드이다.

#include <bits/stdc++.h>
 
#define X first
#define Y second
#define endl '\n'
using namespace std;
 
typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> pll;
 
pll operator + (const pll A, const pll B){
    return {A.X + B.X, A.Y + B.Y};
}
 
pll operator - (const pll A, const pll B){
    return {A.X - B.X, A.Y - B.Y};
}
 
ll dist(pll A, pll B){
    ll dx = A.X - B.X;
    ll dy = A.Y - B.Y;
    return dx * dx + dy * dy;
}
 
ll ccw(pll A, pll B, pll C) {
    return A.X * B.Y + B.X * C.Y + C.X * A.Y - A.Y * B.X - B.Y * C.X - C.Y * A.X;
}
 
void convexHull(vector<pll>& v){
    sort(v.begin(), v.end());
    vector<pll> L, R;
    for (pll i: v){
        while (L.size() >= 2 && ccw(i, L[L.size() - 1], L[L.size() - 2]) <= 0) L.pop_back();
        L.emplace_back(i);
        while (R.size() >= 2 && ccw(i, R[R.size() - 1], R[R.size() - 2]) >= 0) R.pop_back();
        R.emplace_back(i);
    }
    v.clear();
    for (int i=0;i<R.size()-1;i++){
        v.emplace_back(R[i]);
    }
    for (int i=L.size()-1;i>0;i--){
        v.emplace_back(L[i]);
    }
}
 
int main() {
    ios_base::sync_with_stdio(false);cin.tie(nullptr);
 
    int n; cin >> n;
    vector<pll> v;
    for (int i=0;i<n;i++){
        ll x, y; cin >> x >> y;
        v.emplace_back(x, y);
    }
    convexHull(v);
    int r = 0;
    ll ans = 0;
    n = v.size();
    for (int i=0;i<n;i++){
        while (r < n * 2 && ccw(v[i], v[(i + 1) % n], v[(i + 1) % n] + v[(r + 1) % n] - v[r % n]) >= 0){
            ans = max(ans, dist(v[i], v[r % n]));
            r++;
        }
        ans = max(ans, dist(v[i], v[r % n]));
    }
    cout << ans << endl;
 
    return 0;
}

 

BOJ 1310 - 달리기 코스

 위 문제와 동일하다.

 

BOJ 9240 - 로버트 후드

 역시 동일한 문제다. 거리의 제곱을 출력하는 것이 아님에 유의.

 

BOJ 10254 - 고속도로

 최대 거리를 갖는 두 점의 좌표를 출력하는 문제다.

 

BOJ 15420 - Blowing Candles

 두 평행선을 그었을 때 모든 점이 이 사이에 들어갈 때, 이 평행선 사이의 거리의 최솟값을 구하는 문제다. 이것도 Rotating Calipers로 구할 수 있다. 점과 선분 사이의 거리의 최댓값을 찾자. 그리고 이들 중의 최솟값이 답이 된다.

 

BOJ 16525 - Escape, Polygon!

 다각형이 빠져나가려면, 나머지 한 직선의 위치가 중요하다.

이 경우에, 빨간 직선이 다각형의 우측에 존재한다면 다각형은 가둬진다. 그림의 경우는 다각형이 빠져나가는 경우다.

 불가능한 경우를 찾아보자. 한 직선을 고정하고, 나머지 두 직선의 위치를 찾는 문제가 된다. Rotating Calipers로 180도를 넘어가지 않는 선분의 최대 위치 \(r\)을 찾아주자. 그러면, 한 직선이 고정되었기 때문에, 나머지 두 직선은 \(r\)까지의 직선들 중 두 개를 뽑는 경우의 수와 같다. 이를 \(_{N}C_{3}\)에서 빼주며 진행하면 된다. 다음은 이 문제의 정답 코드이다. 파이썬으로 작성했다.

import queue
import sys, math
import heapq
 
from collections import deque
from bisect import bisect_left, bisect_right
 
input = sys.stdin.readline
 
pop = heapq.heappop
push = heapq.heappush
 
 
# input
def ip(): return int(input())
def sp(): return str(input().rstrip())
 
 
def mip(): return map(int, input().split())
def mfp(): return map(float, input().split())
def msp(): return map(str, input().split())
 
 
def lmip(): return list(map(int, input().split()))
def lmsp(): return list(map(str, input().split()))
 
 
############ Main! #############
 
 
def f(A, B, C):
    return A[0] * B[1] + B[0] * C[1] + C[0] * A[1] - A[1] * B[0] - B[1] * C[0] - C[1] * A[0]
 
 
def ccw(A, B, C, D):
    return f(A, B, [D[0] - C[0] + B[0], D[1] - C[1] + B[1]])
 
 
n = ip()
a = [lmip() for _ in range(n)]
 
r = 1
ans = n * (n - 1) * (n - 2) // 6
for i in range(n):
    while ccw(a[i], a[(i + 1) % n], a[r % n], a[(r + 1) % n]) > 0:
        r += 1
    if r - i <= 1:
        continue
    ans -= max(0, (n - (r - i)) * (n - (r - i) - 1) // 2)
 
print(ans)
 
 

 

BOJ 19586 - 울타리

 벡터의 외적은, sin이다. 그렇다면 cos은?

 직사각형의 한 변을 고정하고, 투 포인터로 각각 L, U, R을 찾자. 외적은 180도, 내적은 90도를 확인할 수 있기 때문에, 이를 잘 활용해보자. 다음은 이 문제의 정답 코드이다.

#include <bits/stdc++.h>
 
#define X first
#define Y second
#define endl '\n'
using namespace std;
 
typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> pll;
 
pll operator + (const pll A, const pll B){
    return {A.X + B.X, A.Y + B.Y};
}
 
pll operator - (const pll A, const pll B){
    return {A.X - B.X, A.Y - B.Y};
}
 
ll dot(const pll A, const pll B){
    return A.X * B.X + B.Y * A.Y;
}
 
ll cross(const pll A, const pll B){
    return A.X * B.Y - B.X * A.Y;
}
 
ld dist(pll A, pll B){
    ll dx = A.X - B.X;
    ll dy = A.Y - B.Y;
    return sqrt(dx * dx + dy * dy);
}
 
ll ccw(pll A, pll B, pll C) {
    return A.X * B.Y + B.X * C.Y + C.X * A.Y - A.Y * B.X - B.Y * C.X - C.Y * A.X;
}
 
ld dist2(const pll &x, const pll &y, const pll &z){
    pll t = x + z;
    return 1.0 * abs(cross((y - x), (t - x))) / dist(x, t);
}
 
void convexHull(vector<pll>& v){
    sort(v.begin(), v.end());
    vector<pll> L, R;
    for (pll i: v){
        while (L.size() >= 2 && ccw(i, L[L.size() - 1], L[L.size() - 2]) <= 0) L.pop_back();
        L.emplace_back(i);
        while (R.size() >= 2 && ccw(i, R[R.size() - 1], R[R.size() - 2]) >= 0) R.pop_back();
        R.emplace_back(i);
    }
    v.clear();
    for (int i=0;i<R.size()-1;i++){
        v.emplace_back(R[i]);
    }
    for (int i=L.size()-1;i>0;i--){
        v.emplace_back(L[i]);
    }
}
 
int main() {
    ios_base::sync_with_stdio(false);cin.tie(nullptr);
 
    int n; cin >> n;
    vector<pll> v;
    for (int i=0;i<n;i++){
        ll x, y; cin >> x >> y;
        v.emplace_back(x, y);
    }
    convexHull(v);
    n = v.size();
 
    cout << fixed << setprecision(13);
    if (n == 1) {
        cout << 0;
        return 0;
    }
    if (n == 2){
        cout << 2.0 * dist(v[0], v[1]);
        return 0;
    }
 
    int r1 = 1;
    int r2 = 1;
    int r3 = 1;
    ld ans = 1e16;
    for (int i=0;i<n;i++){
        if (r1 % n == i) r1++;
        if (r2 % n == i) r2++;
        if (r3 % n == i) r3++;
 
        while (ccw(v[i], v[(i + 1) % n], v[(r1 + 1) % n] - v[r1 % n] + v[(i + 1) % n]) > 0 && dot(v[(i + 1) % n] - v[i], v[(r1 + 1) % n] - v[r1 % n]) > 0) r1++;
        while (ccw(v[i], v[(i + 1) % n], v[(r2 + 1) % n] - v[r2 % n] + v[(i + 1) % n]) > 0) r2++;
        while (ccw(v[i], v[(i + 1) % n], v[(r3 + 1) % n] - v[r3 % n] + v[(i + 1) % n]) > 0 || dot(v[(i + 1) % n] - v[i], v[(r3 + 1) % n] - v[r3 % n]) < 0) r3++;
        pll p = v[(i + 1) % n] - v[i];
        ans = min(ans, (dist2(v[i], v[r2 % n], p) + dist2(v[r1 % n], v[r3 % n], {p.Y, -p.X})) * 2);
    }
    cout << ans;
 
 
    return 0;
}

 

BOJ 10466 - Mailing Origami

 위 문제와 동일하다. 다만 이 문제는 직사각형의 최소 넓이를 묻고 있다.

 

BOJ 13310 - 먼 별

 시간이 지남에 따라, 가장 멀리 떨어진 두 별의 거리는 점점 감소하다가 다시 증가한다. 삼분 탐색으로 열심히 찾으면 된다.

과거 파이썬의 흔적

 다음은 이 문제의 정답 코드이다. Rotating Calipers를 go() 함수로 따로 뺐다.

#include <bits/stdc++.h>
 
#define X first
#define Y second
#define endl '\n'
using namespace std;
 
typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> pll;
 
pll operator + (const pll A, const pll B){
    return {A.X + B.X, A.Y + B.Y};
}
 
pll operator - (const pll A, const pll B){
    return {A.X - B.X, A.Y - B.Y};
}
 
ll dist(pll A, pll B){
    ll dx = A.X - B.X;
    ll dy = A.Y - B.Y;
    return dx * dx + dy * dy;
}
 
ll ccw(pll A, pll B, pll C) {
    return A.X * B.Y + B.X * C.Y + C.X * A.Y - A.Y * B.X - B.Y * C.X - C.Y * A.X;
}
 
 
void convexHull(vector<pll>& v){
    sort(v.begin(), v.end());
    vector<pll> L, R;
    for (pll i: v){
        while (L.size() >= 2 && ccw(i, L[L.size() - 1], L[L.size() - 2]) <= 0) L.pop_back();
        L.emplace_back(i);
        while (R.size() >= 2 && ccw(i, R[R.size() - 1], R[R.size() - 2]) >= 0) R.pop_back();
        R.emplace_back(i);
    }
    v.clear();
    for (int i=0;i<R.size()-1;i++){
        v.emplace_back(R[i]);
    }
    for (int i=L.size()-1;i>0;i--){
        v.emplace_back(L[i]);
    }
}
 
vector<pll> d;
 
ll go(int t, vector<pll>& v){
    vector<pll> p;
    for (int i=0;i<v.size();i++){
        p.emplace_back(v[i].X + d[i].X * t, v[i].Y + d[i].Y * t);
    }
    convexHull(p);
    int r = 0;
    int n = p.size();
 
    ll ret = 0;
    for (int i=0;i<n;i++){
        while (r < n * 2 && ccw(p[i], p[(i + 1) % n], p[(r + 1) % n] - p[r % n] + p[(i + 1) % n]) >= 0){
            ret = max(ret, dist(p[i], p[r % n]));
            r++;
        }
        ret = max(ret, dist(p[i], p[r % n]));
    }
    return ret;
}
 
int main() {
    ios_base::sync_with_stdio(false);cin.tie(nullptr);
 
    int n, t; cin >> n >> t;
    vector<pll> v;
    for (int i=0;i<n;i++){
        ll x, y, dx, dy; cin >> x >> y >> dx >> dy;
        v.emplace_back(x, y);
        d.emplace_back(dx, dy);
    }
 
    int lo = 0;
    int hi = t;
    while (lo + 3 <= hi){
        int l = (lo * 2 + hi) / 3;
        int r = (lo + hi * 2) / 3;
        if (go(l, v) <= go(r, v)) hi = r;
        else lo = l;
    }
    ll ans = go(lo, v);
    int idx = lo;
    for (int i=lo;i<=hi;i++){
        ll now = go(i, v);
        if (now < ans){
            idx = i;
            ans = now;
        }
    }
    cout << idx << endl << ans << endl;
 
    return 0;
}

 

 

 

마치며

 다음 글에서는 다각형 내부 점 판정을 다룰 생각이다. 근데 이거 아직도 구현 안해봐서 좀 걸릴 것 같다.

 

 수정 1: 볼록 껍질 코드를 본문에 맞게 수정했습니다.