Geometry & Recognition

Brute-force : $O(n^2)$,

Optimal: $O(n\log n)$

참고: people.csail.mit.edu/indyk/6.838-old/handouts/lec17.pdf

참고: arxiv.org/pdf/1911.01973.pdf

//  points should be sorted in order of increasing x-component.
//
#define DIST2(p, q) ((p).x-(q).x)*((p).x-(q).x) + ((p).y-(q).y)*((p).y-(q).y)
// 수정: index range -> [left, right]; ClosestPair(points, 0, points.size()-1, closest);
int ClosestPair(std::vector<CPoint>& points, int left, int right, int closest[2]) {
    if (right - left >= 3) {
        int lpair[2], rpair[2], min_dist2;
        int mid = left + (right - left) / 2;                    // half index;
        int ldist = ClosestPair(points, left, mid, lpair);     // left side;
        int rdist = ClosestPair(points, mid+1, right, rpair); // right side;
        if (ldist < rdist) {
            closest[0] = lpair[0]; closest[1] = lpair[1];
            min_dist2 = ldist;
        } else {
            closest[0] = rpair[0]; closest[1] = rpair[1];
            min_dist2 = rdist;
        }
        // find points which lies near center strip(2d-strip);
        // Note our distance is the squar of actual distance;
        int width = int(sqrt(double(min_dist2))+ .5);
        int ll = left;
        while (points[ll].x < (points[mid].x - width - 1)) ll++;
        int rr = idx2;
        while (points[rr].x > (points[mid + 1].x + width + 1)) rr--;
        for (int i = ll; i < rr; i++) {
            for (int j = i + 1; j <= rr; j++) {
                int dist2 = DIST2(points[i], points[j]);
                if (min_dist2 > dist2) {
                    min_dist2 = dist2;
                    closest[0] = i; closest[1] = j;
                }
            }
        }
        return min_dist2;
    } 
    else if (right == left + 2) {
        return ClosestPair3(points, left, closest);
    }
    else if (right == left + 1) {
        closest[0] = left; closest[1] = right;
        return DIST2(points[left], points[right]);
    }
    else return INT_MAX;
};

Bezier 곡선은 Bernstein 다항식을 이용하여 표현할 수 있지만, 차수가 높아질수록 다항식을 직접 계산하는 방법은 수치 오차가 커질 수 있다. 이에 비해 De Casteljau's Algorithm은 linear interpolation을 반복하여 Bezier 곡선 위의 점을 계산하므로 수치적으로 매우 안정적이다. 이 알고리즘은 Bezier 곡선을 정의하는 control points으로부터 주어진 매개변수에 해당하는 곡선 위의 한 점을 재귀적으로 계산하는 방법이다.

매개변수 $t(0 \le t \le 1)$가 주어졌을 때, 인접한 제어점들을 다음과 같이 선형 보간한다.

$$ P_i^{(1)}=(1-t)P_i+tP_{i+1}$$여기서 $P_i$는 원래의 제어점이다. 같은 과정을 새롭게 생성된 점들에 대해 반복하면

$$P_i^{(k)}=(1-t)P_i^{(k-1)}+tP_{i+1}^{(k-1)}    $$

을 얻는다. 이 과정을 하나의 점만 남을 때까지 반복하면 최종적으로 얻어지는 점

$$P_0^{(n)}$$이 바로 $t$에서의 Bezier 곡선 위의 점이다.

그렇다면 왜 단순한 linear interpolation을 반복하는 것만으로 Bezier 곡선을 정확하게 계산할 수 있을까? 그 이유는 Bezier 곡선이 linear interpolation을 재귀적으로 적용한 결과와 정확히 동일한 수학적 정의를 갖기 때문이다. 즉, De Casteljau 알고리즘은 Bezier 곡선을 근사하는 방법이 아니라 Bernstein 다항식으로 정의되는 Bezier 곡선을 다른 방식으로 계산하는 알고리즘이다.

이를 이해하기 위해 먼저 1차 Bezier 곡선을 생각해 보자. 제어점이 $P_0$, $P_1$ 두 개뿐인 경우 곡선은 단순한 직선이며,

$$B(t)=(1-t)P_0+tP_1$$으로 표현된다. 이는 두 점 사이의 선형 보간식과 정확히 같다.

이제 2차 Bezier 곡선을 생각해 보자. 먼저 세 제어점 $P_0,P_1,P_2$에 대해

$$Q_0=(1-t)P_0+tP_1,\qquad Q_1=(1-t)P_1+tP_2$$를 계산한다. 이어서 다시 두 점 $Q_0,Q_1$를 선형 보간하면

$$B(t)=(1-t)Q_0+tQ_1$$을 얻는다. 이를 전개하면

\begin{aligned}B(t)&=(1-t)\big[(1-t)P_0+tP_1\big]+t\big[(1-t)P_1+tP_2\big] \\&=(1-t)^2P_0+2t(1-t)P_1+t^2P_2,\end{aligned}이 되는데, 이는 바로 2차 Bernstein 다항식으로 표현한 Bezier 곡선이다.

3차 이상의 경우에도 같은 원리가 성립한다. 선형 보간을 한 단계 수행할 때마다 Bernstein 다항식의 차수가 하나씩 증가하며, 이 과정을 반복하면 최종적으로

$$B(t)=\sum_{i=0}^{n}\binom{n}{i}(1-t)^{n-i}t^iP_i$$를 얻게 된다. 즉, De Casteljau 알고리즘으로 계산한 결과는 Bernstein 다항식으로 정의한 Bezier 곡선과 완전히 동일하면, 이 사실은 수학적 귀납법으로 쉽게 증명할 수 있다.

따라서 De Casteljau 알고리즘은 Bernstein 다항식을 근사하는 방법이 아니라, Bezier 곡선을 재귀적인 선형 보간만으로 정확하게 계산하는 알고리즘이다. Bernstein 다항식을 직접 계산하는 방법보다 계산량은 다소 많지만, 연산이 덧셈과 곱셈만으로 이루어져 큰 이항계수나 고차 다항식 계산에 따른 수치 오차를 효과적으로 줄일 수 있다. 또한 계산 과정에서 곡선을 두 개의 Bezier 곡선으로 자연스럽게 분할(subdivision)할 수 있어 CAD, 컴퓨터 그래픽스, 벡터 폰트 등 다양한 분야에서 널리 사용된다.

// De Casteljau's algorithm; recursive version; slow for larger deg;
double Bezier(int deg, double Q[], double t) {
   if (deg == 0) return Q[0];
   else if (deg == 1) return (1-t)*Q[0] + t*Q[1];
   else if (deg == 2) return (1-t)*((1-t)*Q[0] + t*Q[1]) + t*((1-t)*Q[1] + t*Q[2]);
   else return (1 - t) * Bezier(deg-1, &Q[0], t) + t * Bezier(deg-1, &Q[1], t);
}
// De Casteljau's algorithm(degree=n-1); 
// non-recursive. Bezier() modifies Q's;
double Bezier(int deg, double Q[], double t) {
    if (deg==0) return Q[0];
    else if (deg==1) return (1-t)*Q[0] + t*Q[1];
    else if (deg==2) return (1-t)*((1-t)*Q[0] + t*Q[1]) + t*((1-t)*Q[1] + t*Q[2]);
    
    for (int k = 0; k < deg; k++)
        for (int j = 0; j < (deg - k); j++)
            Q[j] = (1 - t) * Q[j] + t * Q[j + 1];
    return Q[0];
}
std::vector<CfPt> BezierCurve(const std::vector<CfPt> &cntls, const int segments) {
    std::vector<double> xp(cntls.size()), yp(cntls.size());
    std::vector<CfPt> curves(segments + 1);
    for (int i = 0; i <= segments; ++i) {
        double t = double(i) / segments;
        // clone control points; non-rec version modifies inputs;
        for (int k = cntls.size(); k-->0;) {
            xp[k] = cntls[k].x; yp[k] = cntls[k].y;
        }
        curves[i] = CfPt(Bezier(xp.size()-1, &xp[0], t), Bezier(yp.size()-1, &yp[0], t));
    }
    return curves;
};

평면 위에서 주어진 점집합을 포함하는 가장 작은 원을 찾는 문제는 가장 단순하게는 임의의 두 점이 지름이 되는 원과 임의의 세 점이 만드는 외접원을 모두 조사하여 찾으면 된다. 따라서 $O(N^4)$의 복잡도를 가질 것이다. 그러나 이 문제는 점집합의 개수에 비례하는 시간 복잡도($O(N)$)를 가지는 알고리즘이 알려져 있고, 여기서는 재귀적 방법을 이용한 Welzl's algorithm을 구현한다.

Welzl's algorithm은 평면에 주어진 여러 점을 모두 포함하는 최소 크기의 원(Smallest Enclosing Circle)을 구하는 대표적인 무작위 재귀 알고리즘이다. 2차원에서는 최소 원이 경계에 있는 최대 3개의 점에 의해 결정된다는 성질을 이용한다. 알고리즘의 핵심 과정을 정리하면 

1. 입력 점들의 순서를 무작위로 섞고 한 점 $p$를 선택한다.
2. $p$를 제외한 점들에 대해 최소 포함 원을 재귀적으로 구한다.
3. $p$가 그 원의 내부에 있으면 현재 원을 그대로 사용한다.
4. $p$가 원 밖에 있다면, 최종 최소 원은 반드시 $p$를 경계에 포함해야 한다. 따라서 $p$를 경계점 집합 $R$에 추가하고 다시 계산한다.
5. $R$에 점이 3개가 되면 세 점을 지나는 원이 결정되므로 재귀를 종료한다. $R$이 2개라면 두 점을 지름의 양 끝으로 하는 원, 1개라면 그 점을 중심으로 하는 반지름 0의 원에서 시작할 수 있다.

즉, 모든 점으로 가능한 원을 조사하는 대신 원 밖으로 벗어나는 점만 최소 원의 경계를 결정하는 후보로 추가하는 것이 핵심이다. 무작위화 덕분에 평균적으로 경계점이 변경되는 경우가 드물어 기대 $O(N)$의 효율적인 계산이 가능하다. 이러한 이유로 Welzl 알고리즘은 최소 포함 원이나 3차원의 최소 포함 구(bounding sphere)를 계산하는 데 널리 알려진 방법이다.

int MinEncCir ( CfPt *P, int n, CfPt* boundary, int b, CfPt& center, double& rad ) {
    // exiting cases
    if ( b == 3 ) CalcCircle3 ( boundary, center, rad ); // a circle passing 3 points;
    else if ( ( n == 1 ) && ( b == 0 ) ) {
        rad = 0; b = 1;
        center = boundary[0] = P[0];
    } 
    else if ( ( n == 2 ) && ( b == 0 ) ) {
        boundary[0] = P[0]; boundary[1] = P[1]; b = 2;
        CalcCircle2 (boundary, center, rad );  // a circle with diagonal consisting of 2 points;
    }
    else if ( ( n == 0 ) && ( b == 2 ) ) {
        CalcCircle2 ( boundary, center, rad );
    }
    else if ( ( n == 1 ) && ( b == 1 ) ) {
        boundary[1] = P[0]; b = 2;
        CalcCircle2 ( boundary, center, rad );
    }
    else {// general case; ( b < 3 ) && ( n + b > 2 )
        // choose a random pivot;
        int k = rand() % n;
        if ( k != 0 ) SWAP( P[0], P[k] );
        int b1 = MinEncCir ( &P[1], n - 1, boundary, b, center, rad );
        if (!InCircle(P[0], center, rad) ) {
            // Now, P[0] belongs to the boundary.
            boundary[b++] = P[0];
            return MinEncCir ( &P[1], n - 1, boundary, b, center, rad );
        } else return b1;
    }
    return b;
}