백준 2389 세상의 중심에서... C++ / Welzl algorithm

https://www.acmicpc.net/problem/2389

 

최대 100개의 점을 모두 포함하는 원의 중심과 반지름을 구하는 문제입니다. 2차원에서 원은 3개의 점으로 구할 수 있습니다. 2개의 점을 지나는 최소 원은 선분을 지름으로 하는 원입니다. 

(이와 같은 논리로 d차원에서는 affine independent한 d+1개의 점이 있으면, 그 점을 지나는 boundary는 유일하게 결정됩니다.)

(참고 : https://freshrimpsushi.github.io/ko/posts/2385/)

 

100개의 점을 대상으로 100C1, 100C2, 100C3 케이스에 대해 원을 구성하고, 100개의 점들과 포함유무를 확인하는 브루트포스 방식도 가능합니다. (대충 복잡도 계산해보면 1600만 정도 돌아야 함)

이러한 방식의 복잡도가 O(n^4) 이지만 (n : point 개수), 모든 d 차원에서 O(n) 으로 최소인 원 (최소 내포 디스크) 을 구할 수 있는 방식이 Welzl (벨츨) 알고리즘입니다. 

 

Welzl 알고리즘의 동작 과정은 아래와 같습니다. 

P : 아직 처리하지 않은 점들의 집합

R : 현재 최소 구(sphere) 경계에 반드시 포함되어야 하는 점들의 집합 (최대 d+1개)

 

1) 기저 케이스 계산

P가 비어있거나 (더 이상 탐색할 거 없음) R이 d+1개의 점을 가질 때, R에 의해 결정되는 구를 계산한다. 

즉, 2차원의 경우, R에 3개의 점이 있다면, 3개의 점을 이용해 3개의 점을 지나는 최소 원을 구한다. 

 

2) 재귀 처리

P에서 임의의 점 p를 제거한다. 그리고 P - {p}, R을 대상으로 재귀적으로 후보가 되는 최소 구 D를 계산한다. 

p가 D 내부에 있다면 D는 P의 모든 점을 포함하는 최소 구가 된다. 

p가 D 외부에 있다면 P - {p}, R + {p} 를 대상으로 다시 D를 계산한다. (즉, 기저 케이스 계산 시, 포함 안 된 p를 반드시 포함하도록 하여 계산) 

 

Welzl에서는 속도 향상을 위해 랜덤하게 p를 선택합니다. 그리고 문제와 같은 2차원이 아닌 3,4,5,... 차원에 대해서도 기저 케이스 계산만 잘 한다면 동일한 방식으로 최소 내포 디스크를 구할 수 있습니다. 

 

알고리즘 상에서는 p를 랜덤하게 골랐지만, 문제 풀 때는 그냥 순차적으로 p를 골랐습니다. 

코드는 아래와 같습니다. (소수점 계산 오차 eps = 1e-9 정도로 잡았습니다.)

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>

using namespace std;

int N;
double x, y;

struct Point {
    double x;
    double y;
};

struct Circle {
    Point p;
    double r;
};

vector<Point> points;
double eps = 1e-9;

void init() {
    cin >> N;
    for(int i = 0; i < N; i++) {
        cin >> x >> y;
        points.push_back({x, y});
    }
}

double dist(Point &p1, Point &p2) {
    double d = 0;
    d += (p1.x - p2.x) * (p1.x - p2.x);
    d += (p1.y - p2.y) * (p1.y - p2.y);
    return sqrt(d);
}

Circle circleFromTwoPoints(Point &p1, Point &p2) {
    Circle c;
    c.p.x = (p1.x + p2.x) / 2;
    c.p.y = (p1.y + p2.y) / 2;
    c.r = dist(p1, p2) / 2;
    return c;
}

Circle circleFromThreePoints(Point &p1, Point &p2, Point &p3) {
    // (x - x1)^2 + (y - y1)^2 = (x - x2)^2 + (y - y2)^2 = (x - x3)^2 + (y - y3)^2
    // 위 식 이용해서 행렬식 구하면 외심 구할 수 있음 
    double det = 2 * (p1.x * (p2.y - p3.y) + p2.x * (p3.y - p1.y) + p3.x * (p1.y - p2.y));
    double x = (p1.x * p1.x + p1.y * p1.y) * (p2.y - p3.y) 
             + (p2.x * p2.x + p2.y * p2.y) * (p3.y - p1.y)
             + (p3.x * p3.x + p3.y * p3.y) * (p1.y - p2.y);
    x /= det;
    double y = (p1.x * p1.x + p1.y * p1.y) * (p3.x - p2.x) 
             + (p2.x * p2.x + p2.y * p2.y) * (p1.x - p3.x)
             + (p3.x * p3.x + p3.y * p3.y) * (p2.x - p1.x);
    y /= det;

    Circle c;
    Point p = {x, y};
    c.p = p;
    c.r = dist(p, p1);
    return c;
}

bool isThreePointsInLine(Point &p1, Point &p2, Point &p3) {
    // 삼각형 넓이 0이면 일직선 상에 있는 것 
    // 외적 이용하면 넓이 나옴 
    double x1 = p2.x - p1.x;
    double y1 = p2.y - p1.y;
    double x2 = p3.x - p1.x;
    double y2 = p3.y - p1.y;

    // 부동소수점 오차 정도는 허용 가정 
    double s = (x1 * y2 - x2 * y1) / 2;
    return fabs(s) < eps;
}

Circle circleFromPoints(vector<Point> &R) {
    Circle c;

    // base cases
    if(R.size() == 0) {
        c.p = {0, 0};
        c.r = 0;
    }
    else if(R.size() == 1) {
        c.p = R[0];
        c.r = 0;
    }
    else if(R.size() == 2) {
        // 같은 점인 경우 
        if(R[0].x == R[1].x && R[0].y == R[1].y) {
            c.p = R[0];
            c.r = 0;
        } else {
            // 두 점 중심
            c = circleFromTwoPoints(R[0], R[1]);
        }
    }
    else {
        // 점 3개로 원 만들기
        // 3개 점 일치 케이스 
        if(R[0].x == R[1].x && R[0].y == R[1].y && R[0].x == R[2].x && R[0].y == R[2].y) {
            c.p = R[0];
            c.r = 0;
        } 
        // 2개 점 일치 케이스 
        else if(R[0].x == R[1].x && R[0].y == R[1].y) {
            c = circleFromTwoPoints(R[0], R[2]);
        }
        else if(R[0].x == R[2].x && R[0].y == R[2].y) {
            c = circleFromTwoPoints(R[0], R[1]);
        }
        else if(R[1].x == R[2].x && R[1].y == R[2].y) {
            c = circleFromTwoPoints(R[0], R[1]);
        }
        // 일직선 상에 있는 경우 
        else if(isThreePointsInLine(R[0], R[1], R[2])) {
            // 거리 가장 먼 걸로 처리 
            Circle c1 = circleFromTwoPoints(R[0], R[1]);
            Circle c2 = circleFromTwoPoints(R[0], R[2]);
            Circle c3 = circleFromTwoPoints(R[1], R[2]);
            if(c1.r >= max(c2.r, c3.r)) c = c1;
            if(c2.r >= max(c1.r, c3.r)) c = c2;
            if(c3.r >= max(c1.r, c2.r)) c = c3;
        }
        // 3개로 원 만드는 경우 
        else {
            c = circleFromThreePoints(R[0], R[1], R[2]);
        }
    }

    return c;
}

bool isInside(Circle &c, Point &p) {
    // c 안에 p 들어가는지 체크 
    // 오차 범위 : eps 
    return dist(c.p, p) < (c.r + eps);
}

Circle Welzl(vector<Point> &P, int idx, vector<Point> R) {
    // P : 전체 포인트 집합
    // idx : 제외할 포인트의 index
    // R : bound circle 계산할 포인트 집합 
    if(idx == 0 || R.size() == 3) {
        // idx == 0 : 더 이상 탐색할 거 없음 
        // R.size() == 3 : 3개 점으로 원 만들기 (2-dim)
        return circleFromPoints(R);
    }

    // idx - 1번째 제외
    Point p = P[idx - 1];
    Circle D = Welzl(P, idx - 1, R);
    // 제외되었던 것도 circle 안에 포함 
    if(isInside(D, p))
        return D;

    // 제외한 포인트 포함해서 다시 계산 
    R.push_back(p);
    return Welzl(P, idx - 1, R);
}

int main() {
    ios_base::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);

    init();
    vector<Point> R;
    Circle result = Welzl(points, points.size(), R);
    cout << result.p.x << ' ' << result.p.y << ' ' << result.r;
}