Approximate Convex Hull

주어진 점집합의 Convex Hull을 구하는 대부분의 정확한 알고리즘은 점들을 일정한 순서로 정렬하는 과정을 포함하므로 일반적으로 \(O(n\log n)\)의 시간 복잡도를 갖는다. 영상에서 추출한 점집합은 매우 많은 점을 포함할 뿐만 아니라, 중복점이나 동일한 직선 위에 놓인 점과 같은 degeneracy가 자주 발생하므로, 정확한 Convex Hull 알고리즘을 구현할 때에는 이러한 특수한 경우를 적절히 처리해야 한다.

한편, 영상처리에서는 Convex Hull이 물체의 외곽을 대략적으로 표현하는 용도로 사용되는 경우가 많아, 정확한 Convex Hull 대신 근사적인 결과만으로도 충분한 경우가 많다. Convex Hull 알고리즘에서 가장 많은 계산량을 차지하는 단계는 점들을 정렬 과정이다. 이 단계를 근사적으로 수행하면 계산량을 크게 줄일 수 있다. 이를 위해 먼저 점집합을 수평(또는 수직) 방향으로 일정한 폭을 갖는 여러 개의 bin으로 분할한다. 이 과정은 모든 점을 한 번만 순회하면 되므로 \(O(n)\)의 시간에 수행할 수 있다.

다음으로, 각 bin에서 bin을 나눈 방향에 수직인 방향으로 가장 바깥쪽에 위치한 두 점을 대표점으로 선택한다. 예를 들어, 수평 방향으로 bin을 나누었다면 각 bin에서 가장 위쪽(\(y\)=max)과 가장 아래쪽(\(y\)=min)의 점을 대표점으로 선택한다. 반대로 수직 방향으로 bin을 나누었다면 가장 왼쪽(\(x\)=min)과 가장 오른쪽(\(x\)=max)의 점을 대표점으로 선택한다.

이렇게 선택된 대표점들만을 이용하여 Convex Hull을 계산하면 원래 점집합의 근사적인 Convex Hull을 얻을 수 있다. 대표점들은 bin의 순서대로 이미 정렬되어 있으므로 추가적인 정렬 과정이 필요하지 않으며, Convex Hull 계산은 대표점의 개수(즉, bin의 개수)를 \(m\)이라 할 때 \(O(m)\)의 시간에 수행된다. 따라서 전체 알고리즘의 시간 복잡도는 \(O(n+m)\)이 된다.

구간(bin)의 개수를 증가시키면 대표점이 원래 점집합의 외곽을 더욱 세밀하게 반영하게 되므로, 계산된 Convex Hull은 점차 실제 Convex Hull에 가까워지는 경향이 있다. 반대로 bin의 개수를 줄이면 계산 속도는 빨라지지만 근사 오차는 증가한다. 따라서 bin의 개수는 계산 시간과 정확도 사이의 trade-off를 결정하는 중요한 매개변수가 된다. 즉, 이 알고리즘은 점들의 정확한 정렬 대신 공간을 일정한 구간으로 binning 하여, 각 구간을 대표하는 극점만을 이용하여 Convex Hull을 계산함으로써 계산량을 크게 줄이는 근사 Convex Hull 알고리즘이라 할 수 있다.

 

아래 그림은 500,000개 점의 근사적인 convex hull을 보여준다(bins=128)

 

// cross(P1-P0, P2-P0);
// = := colinear;
// < := <p0,p1,p2> CW order;
// > := <p0,p1,p2> CCW order;
static
int ccw(const CPoint &P0, const CPoint &P1, const CPoint &P2 ) {
    return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
}
std::vector<CPoint> approximateHull(const std::vector<CPoint>& Q, int bins) {
    if (Q.size() < 3) return std::vector<CPoint> ();
    if (bins < 0) bins = max(1, Q.size()/10);
    int Lmin = 0, Lmax = 0, Rmin = 0, Rmax = 0;
    int left = Q[0].x, right = left;
    for (int i = 1; i < Q.size(); i++) {
        if (Q[i].x <= left) {
            if (Q[i].x < left) {
                Lmin = Lmax = i;
                left = Q[i].x;
            } 
            else 
                if (Q[i].y < Q[Lmin].y) Lmin = i;
                else if (Q[i].y > Q[Lmax].y) Lmax = i;
        }
        if (Q[i].x >= right) {
            if (Q[i].x > right) {
                Rmin = Rmax = i;
                right = Q[i].x;
            }
            else 
                if (Q[i].y < Q[Rmin].y) Rmin = i;
                else if (Q[i].y > Q[Rmax].y) Rmax = i;
        }
    }
    if (left == right) return std::vector<CPoint> (); //vertical line;
    //
    std::vector<int> slotYhigh(bins + 2, -1);
    std::vector<int> slotYlow(bins + 2, -1);
    slotYlow.front() = Lmin; 
    slotYhigh.front() = Lmax;
    slotYlow.back()  = Rmin;
    slotYhigh.back()  = Rmax;
    //
    const CPoint &A = Q[Lmin];
    const CPoint &B = Q[Rmin];
    const CPoint &C = Q[Lmax];
    const CPoint &D = Q[Rmax];
    for (int i = 0; i < Q.size(); i++) {
        if (Q[i].x == right || Q[i].x == left) continue;
        if (ccw(A, B, Q[i]) < 0) { // below line(A,B);
            int b = 1 + (bins * (Q[i].x - left)) / (right - left);
            if (slotYlow[b]==-1) slotYlow[b] = i;
            else if (Q[i].y < Q[slotYlow[b]].y) slotYlow[b] = i;
            continue ;
        }
        if (ccw(C, D, Q[i]) > 0) {// above line(C,D);
            int b = 1 + (bins * (Q[i].x - left)) / (right - left);
            if (slotYhigh[b]==-1) slotYhigh[b] = i;
            else if (Q[i].y > Q[slotYhigh[b]].y) slotYhigh[b] = i;
            continue;
        }
    }
    std::vector<CPoint> hull(Q.size());
    int s = -1;
    for (int i = 0; i <= bins + 1; i++) {
        if (slotYlow[i]==-1) continue;
        while (s > 0) 
            if (ccw(hull[s-1], hull[s], Q[slotYlow[i]]) > 0) break;
            else s--;
        
        hull[++s] = Q[slotYlow[i]];
    };
    if (Rmax != Rmin) hull[++s] = Q[Rmax];
    //
    int s0 = s ;
    for (int i = bins; i >= 0; i--) {
        if (slotYhigh[i]==-1) continue;
        while (s > s0) 
            if (ccw(hull[s-1], hull[s], Q[slotYhigh[i]]) > 0) break;
            else s--;
        
        hull[++s] = Q[slotYhigh[i]];
    }
    if (hull[s] == hull[0]) s--;
    hull.resize(s+1);
    return hull;
}