728x90

2차원 점집합의 median

'Computational Geometry' 카테고리의 다른 글

Bresenham's Line Algorithm  (0) 2021.04.07
Rotating Calipers  (3) 2021.03.31
Convex Hull Peeling  (0) 2021.03.29
Graham Scan  (0) 2021.03.28
Jarvis March  (0) 2021.03.26
Approximate Minimum Enclosing Circle  (1) 2021.03.18
Posted by helloktk

댓글을 달아 주세요

Graham Scan

Computational Geometry 2021. 3. 28. 10:35
728x90

Convex hull을 찾을 때, 중복점과 일직선상 위에 있는 점들에 대해서 적절한 처리를 하지 않으면 효율적으로 작성된 알고리즘이라도 정상적으로 동작하지 않을 수 있다. 그림은 Graham Scan 알고리즘을 이용해서 50,000개 점들의 convex hull을 구한 결과다. qsort()을 호출할 때 compare()에서 global 변수를 쓰지 않기 위해서 pivot point를 뺀 후 사용한다. 정렬 후에는 다시 더해준다.

inline int LeftSide(CPoint A, CPoint B, CPoint C) {
    return (B.x - A.x) * (C.y - A.y) - (C.x - A.x) * (B.y - A.y); // cross(AB, AC);
}
static int cross(CPoint *p1, CPoint *p2) {
    return p1->x * p2->y - p1->y * p2->x;
}
static int norm2(CPoint *p) {
    return p->x * p->x + p->y * p->y;
}
static int compare(const void *a, const void *b) {
    CPoint *p1 = (CPoint *)a, *p2 = (CPoint *)b;
    int ccw = cross(p2, p1);
    // 거리가 먼 곳 --> 가까운 곳;
    if (ccw == 0) return norm2(p2) > norm2(p1) ? 1 : -1;
    return ccw; // 반시계방향 정렬;
}
static void angularSort(std::vector<CPoint> & pts) {
    if ((pts.size()) < 2) return;
    int id = 0;
    for (int i = pts.size() - 1; i >= 1; --i) 
        if ((pts[id].y > pts[i].y) || (pts[id].y == pts[i].y && pts[i].x > pts[id].x)) 
            id = i;  //lowest--y-value; //largest--x-value;
    std::swap(pts[0], pts[id]);
    for (int i = pts.size() - 1; i >= 1; --i) pts[i] -= pts[0];
    qsort(&pts[1], pts.size() - 1, sizeof(CPoint), compare);
    for (int i = pts.size() - 1; i >= 1; --i) pts[i] += pts[0];
};
// AB와 BC가 나란할 때, AB 방향과 BC방향이 같은가(0) 반대면 (1): dot(AB, BC);
static int folded(CPoint A, CPoint B, CPoint C) {
    return (B.x - A.x) * (C.x - B.x) + (B.y - A.y) * (C.y - B.y) <= 0;
}
int GrahamScan(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    if (pts.size() < 3) return 0;
    angularSort(pts);
    hull.resize(pts.size());
    hull[0] = pts[0];
    int i = 1; 
    while (i < pts.size()) { 
        // peek out points degenerate with pts[0];
        if (pts[i] != hull[0]) {
            hull[1] = pts[i++]; break;  //  (i++)
        }
        i++;
    }
    if (i == pts.size()) return 0;
    int k = 2; // 새로 추가되는 hull index;
    while (i < pts.size()) {
        // (i)가 hull의 left에 있으면(1) hull의 꼭지점에 추가;
        int ccw = LeftSide(hull[k - 2], hull[k - 1], pts[i]); 
        if (ccw > 0) hull[k++] = pts[i++];
        // (i)가 hull의 직전 변과 collinear한 데 반대방향을 향하면 skip;
        else if (ccw == 0 && folded(hull[k - 2], hull[k - 1], pts[i])) ++i;
        // (i)가 hull의 right 또는 위에 있으면(직전 변과 같은 방향)
        // hull의 직전 꼭지점을 순차적으로 제거한 후 다시 검사;
        else --k;
        ASSERT(k >= 2);
    }
    hull.resize(k);
    return hull.size();
};

'Computational Geometry' 카테고리의 다른 글

Rotating Calipers  (3) 2021.03.31
Convex Hull Peeling  (0) 2021.03.29
Graham Scan  (0) 2021.03.28
Jarvis March  (0) 2021.03.26
Approximate Minimum Enclosing Circle  (1) 2021.03.18
Minimum Bounding Rectangle  (0) 2021.03.15
Posted by helloktk

댓글을 달아 주세요

728x90

기하 알고리즘에서 입력점의 개수가 많아지면 중복점이 생기고, 이 때문에 최적화에 주안점을 둔 알고리즘은 중복점에 대한 예외처리를 하지 않으면 정상적으로 동작하지 않는 경우가 종종 생긴다. 안정성과 구현의 간결함을 유지하고 싶으면 좀 더 brute force method에 가까운 알고리즘을 사용하는 것도 고려해 볼 수 있다. 다음은 Jarvis March 알고리즘을 써서 2차원 convex hull을 구하는 코드이다. (그림은 50,000개의 입력점을 사용함.)

BOOL ccw(CPoint a, CPoint b, CPoint c){	
    // cross(ac, ab) < 0: ab가 ac보다 오른쪽에 있음 = b가 ac변의 오른쪽에 있음.
    return ((c.x - a.x) * (b.y - a.y) - (c.y - a.y) * (b.x - a.x)) < 0;
}
int JarvisMarch(std::vector<CPoint> &pts, std::vector<CPoint> &hull) {
    int n = pts.size();
    if (n < 3) return 0;
    int left = 0;
    for (int i = 1; i < n; i++) //convex hull의 첫번째는 맨 왼쪽 점;
        if (pts[i].x < pts[left].x) left = i;
    hull.clear();
    int turns = 0; //debug 용도
    int prev = left;
    while (1) {
        hull.push_back(pts[prev]);
        int next = (prev + 1) % n;
        // prev-next 변이 모든 점의 오른쪽에 있게 하는 next를 찾으면 됨;
        // 없으면 next가 convex hull 위의 점임; 
        for (int i = 0; i < n; i++) {
            if (i == next || i == prev) continue;
            if (ccw(pts[prev], pts[i], pts[next]))  
                next = i;
        }
        prev = next;
        turns++;     //debug 용도;
        if (prev == left) break;
    } 
    TRACE("number of turns = %d\n", turns);
    remove_collinear(hull); // removes collinear vertices(optional);
    return hull.size();
}
더보기
int remove_collinear(std::vector<CPoint>& V) {
    int N = V.size();
    if (N < 3) return 0;
    // 먼저, 일직선상에 있지 않는 vertex을 찾음;
    int start; 
    for (start = 0; start < N; start++)
        if (ccw(V[(start - 1 + N) % N], V[start], V[(start + 1) % N]))
            break;

    std::vector<CPoint> H;
    H.push_back(V[start]);
    int prev = start;
    int curr = (prev + 1) % N;
    int next = (curr + 1) % N;
    while (curr != start) {
        if (ccw(V[prev], V[curr], V[next])) {
            H.push_back(V[curr]);
            prev = curr;
        } 
        curr = (curr + 1) % N;
        next = (curr + 1) % N;
    }
    std::swap(V, H);
    return V.size();
};

'Computational Geometry' 카테고리의 다른 글

Convex Hull Peeling  (0) 2021.03.29
Graham Scan  (0) 2021.03.28
Jarvis March  (0) 2021.03.26
Approximate Minimum Enclosing Circle  (1) 2021.03.18
Minimum Bounding Rectangle  (0) 2021.03.15
Minimum Enclosing Circle  (0) 2021.03.01
Posted by helloktk

댓글을 달아 주세요

728x90

주어진 점집합을 포함하는 최소 면적의 rectangle을 찾는 문제는 점집합의 convex hull을 포함하는 rectangle을 찾는 문제와 같다. convex hull의 한 변을 공유하면서 전체를 포함하는 rectangle 중 최소 면적인 것을 찾으면 된다. Rotating Caliper를 쓰면 O(n) step이 요구된다.

void BoundingRect(std::vector<CPoint>& hull, CPoint vtx[4]) {//non-optimal.
    if (hull.size() < 3) return;
    double minArea = DBL_MAX;
    int bIndex = 0, tIndex = 0, lIndex = 0, rIndex = 0;
    for (int j = 0, i = hull.size() - 1; j < hull.size(); i = j++) {
        CPoint &A = hull[i], &B = hull[j];
        CPoint BA = B - A;  // convex hull의 한 기준변;
        CPoint BAnormal = CPoint(-BA.y, BA.x); //(B-A)에 반시계방향으로 수직인 벡터;
        double BAsq = DOT(BA, BA);
        //기준변과 대척점(antipodal point)을 찾는다: BAnormal방향으로 가장 멀리 떨어진 점
        int id = FindAntipodal(BAnormal, hull) ;
        CPoint QA = hull[id] - hull[i];
        //기준변과 대척점을 통과하는 직선과의 사이거리;
        double D1 = fabs(double(DOT(BAnormal, QA))) / sqrt(BAsq);
        //left_side_end_point;//기준변 반대방향으로 가장 멀리 떨어진 꼭지점;
        int id1 = FindAntipodal(-BA, hull);
        //right_side_end_point;//기준변 방향으로 가장 멀리 떨어진 꼭지점;
        int id2 = FindAntipodal(BA, hull);

        ///가장 왼쪽과 가장 오른쪽 꼭지점을 연결하는 선분
        CPoint H = hull[id1] - hull[id2];  
        //기준변 방향으로 정사영하면 두 평행선 사이거리를 줌;
        double D2 = fabs(double(DOT(BA,H))) / sqrt(BAsq);
        double area = D1 * D2; 
        if (area < minArea) {
            minArea = area;
            bIndex = i ;  tIndex = id;
            lIndex = id1; rIndex = id2;
        }
    }
    //directional vector for base_edge;
    CPoint dir = hull[(bIndex + 1) % hull.size()] - hull[bIndex] ; 
    vtx[0] = LineIntersect(hull[bIndex], dir, hull[lIndex]);
    vtx[1] = LineIntersect(hull[tIndex], dir, hull[lIndex]);
    vtx[2] = LineIntersect(hull[tIndex], dir, hull[rIndex]);
    vtx[3] = LineIntersect(hull[bIndex], dir, hull[rIndex]);
}

'Computational Geometry' 카테고리의 다른 글

Jarvis March  (0) 2021.03.26
Approximate Minimum Enclosing Circle  (1) 2021.03.18
Minimum Bounding Rectangle  (0) 2021.03.15
Minimum Enclosing Circle  (0) 2021.03.01
단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
Posted by helloktk

댓글을 달아 주세요

728x90

평면 상에 주어진 점들의 convex hull을 구하는 알고리즘 중의 하나인 Graham scan에서는 먼저 주어진 점들을 한 점을 기준으로 각도로 정렬하는 과정이 필요했다. 그러면 점들이 순차적으로 연결된 단순 다각형에서는 sorting과정이 없이 Graham scan 알고리즘을 바로 적용이 가능한 것일까? 이에 대한 counter example은 쉽게  찾을 수 있다. 단순 폴리곤에 대해서도 항상 각도에 대한 정렬 과정이 필요한 것인가? 답은 아니다. 정렬 과정이 없이도 단순 폴리곤에 대해서는 쉽게 convex hull을 구할 수 있는 알고리즘이 존재한다. 정렬 과정이 없이 단순 폴리곤의 convex hull을 찾는 알고리즘에 대해서 알아보자.

Melkman Algorithm;

우선 폴리곤의 이웃하는 세 꼭짓점을 잡아서, 반시계 방향의 삼각형을 구성한다. 이 삼각형을 deque에 넣는다 (bottom = top). 폴리곤을 순환하면서 새 점이 들어올 때 이미 만들어진 convex hull의 내부점이 아니면 이 점이 포함되도록 convex hull을 업데이트한다: Graham scan 알고리즘의 scanning 과정을 bottom을 기준으로 반시계 방향으로 convexity를 만족시킬 때까지 bottom을 제거하고, top을 기준으로는 시계방향으로 convexity를 만족시킬 때까지 top을 제거한다. 이 과정이 끝나면 새 점을 deque에 추가한다. 이 과정을 나머지 모든 점들에 대해서도 진행한다.

int LEFTSIDE(CPoint A, CPoint B, CPoint C){ // cross(AB, AC) > 0: C가 AB 대해서 왼쪽
	return ((B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x)) > 0;
}
int convexhull_spolygon(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    int N = pts.size();
    if (N < 3) return 0;
    std::vector<CPoint> Q(2 * N + 1);  // deque;
    int bot = N - 2, top = bot + 3;
    // |0|1|......|N-2|N-1|N+0|N+1|N+2|..............|2N|;
    //              bot    top    
    Q[bot] = Q[top] = pts[2];
    if (LEFTSIDE(pts[0], pts[1], pts[2])) { //2->0->1->2; Q에 ccw로 입력;
        Q[bot + 1] = pts[0]; Q[bot + 2] = pts[1];
    } else {
        Q[bot + 1] = pts[1]; Q[bot + 2] = pts[0];
    }
    int i = 2; // last touched index;
    while (++i < N) {
        // 기존의 convex_hull에 들어있으면 제외.
        if (LEFTSIDE(Q[bot + 0], Q[bot + 1], pts[i]) && 
            LEFTSIDE(Q[top - 1], Q[top + 0], pts[i]))
            continue;
        // bottom에 대해서 ccw 방향으로 체크(bot 증가 방향)
        // pts[i]가 (bot)->(bot+1)라인의 오른쪽에 있으면, 기존의 bottom을 제외;
        while (!LEFTSIDE(Q[bot], Q[bot + 1], pts[i]) && (bot < top)) bot++;
        Q[--bot] = pts[i];
        // 추가점에서 top->top-1의 ccw 방향으로 convexity 체크하여 만족하지 않은 경우
        // top을 감소
        while (!LEFTSIDE(Q[top - 1], Q[top], pts[i]) && (top > bot)) top-- ;
        Q[++top] = pts[i];
    }
    hull.resize(top - bot);
    for (int i = hull.size() - 1; i >= 0; --i) hull[i] = Q[i + bot];
    return hull.size();
};

**네이버 블로그 이전;

 

 

'Computational Geometry' 카테고리의 다른 글

Minimum Enclosing Circle  (0) 2021.03.01
단순 다각형 만들기  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Posted by helloktk

댓글을 달아 주세요

728x90

사진 속에 포함된 물체를 판별하기 위해서는 보통 이미지를 이진화시켜 binary 이미지를 만들고, conneted component labeling 등을 이용하여 영상에서 겹치지 않는 다수의 물체들을 분리해 낸다. 분리된 물체의 외각을 감싸는 convex hull을 구하기 위해서는 우선 경계점들을 추적하여 폴리곤을 형성한 다음 이들의 convex hull 알고리즘을 이용함으로써 구할 수 있다. convex hull 알고리즘 중에서 chain-hull이나 Melkman 알고리즘은 점들이 일정한 규칙에 의해서 정렬되어 있다는 사실을 이용하여 빠르게 convex hull을 찾을 수 있게 한다.

 

이미지를 읽는 raster-scan 방식은 물체의 경계점들을 일정하게 정렬된 순서로 읽게 되므로(chain-hull에서도 동일한 순서가 들어오지만, 전체 점들을 다 얻은 다음에 다시 분류작업을 해야 하므로, 읽은 즉시 바로 convex  hull을 구성하는 알고리즘은 제공하지 못한다) on-site convex hull의 구성이 가능하게 된다. 아래의  코드는 binary raster image를 읽으면서 경계점들의 convex hull을 찾는다.

#define bg     0x00             /* foreground color (black) */ 
#define fg     0xFF             /* background color (white) */
//cross(AB, AC) > 0 if C is lying on the left side of edge(AB);
#define CCW(A,B,C) ((B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x))

/* find a convex hull of foreground object in a binary image; */ 
/* in some case L[0] = R[0], or L[ll] = R[rr] if first line or last line of object 
** is composed of a single point (수정본:: 배경 체크 부분);
** convex hull이 중복점을 갖는 경우를 허용하지 않을 때는 comment에 따라서 처리해야 한다.
*/ 
int RasterHull(BYTE *image, int w, int h, CPoint H[/* 2xh */]) { 
    CPoint *L = &H[0]; //stack of left-side hull; 
    CPoint *R=  &H[h]; //stack of right side hull; 
    int rr = -1, ll = -1, x; 
    CPoint q; 
    for (int y = 0 ; y < h; y++) { 
        for (x = 0; x < w; x++)
            if (image[y * w + x] != bg) break; 
        if (x == w) continue; 
        q = CPoint(x, y);
        while (ll > 0) {
            // q가 이전 edge의 오른쪽이면 convex hull의 꼭지점이 됨;
            if (CCW(L[ll - 1], L[ll], q) < 0) break; 
            else ll--; 
        } 
        L[++ll] = q; 
        for (x = w - 1; x >= 0; x--) //x=-1 never occurs; 
            if (image[y * w + x] != bg) break; 
        // 행에 싱글 픽셀 있는 경우 왼쪽과 오른쪽이 점이 중복이 되는 경우가
        // 생긴다. 이 중복을 없애고 싶으면
        if (x <= q.x) continue ;
        q = CPoint(x, y); 
        while (rr > 0) { 
            // q가 이전 edge의 왼편에 있으면 convex hull의 꼭지점이 됨;
            if (CCW(R[rr - 1], R[rr], q) > 0) break; 
            else rr--; 
        } 
        R[++rr] = q; 
    } 
    // 왼쪽 처음에 점과 오른쪽 처음 점유 중복될 가능성이 있다(싱글 픽셀)
    // 왼쪽 마지막 점과 오른쪽 마지막 점이 중복이 되는 경우도 있다(싱글 픽셀) 
    // convex hull에서 중복점을 완전히 제외하고 싶으면 이를 체크하는 부분이 필요;
    for (int i = 0; i <= ll; i++) H[i] = L[i]; //left hull;
    // remove bottom degeneracy; 
    if (H[ll] == R[rr]) rr--;
    for (int j = rr; j >= 0; j--) H[i++] = R[j];//right hull; 
    //remove top degeneracy;
    if (i > 1 && H[0] == H[i - 1]) i--;
    return (i); 
};

위 함수는 픽셀 값이 bg가 아닌 것은 모두 전경으로 처리하도록 써져 있는데, 만약 connected component를 구한 상태에서 특정 component의 convex hull을 구하고자 하는 경우에는 if-문안의 조건을 image[y * w + x] == index로 바꾸어야 한다 (여기서 index는 connected component label이다. 그리고 labeling 이미지는 일반적으로 정수 이미지이므로 그에 맞게 수정을 해야 한다.) 

 

data matrix 코드를 포함하는 이미지를 이진화시킨 후 connected component labeling을 하여 작은 blob을 제거하면 거의 코드만 남길 수 있다. 그 후 위의 알고리즘으로  convex hull(붉은색 라인)을 구성한 것을 보인 것이다.

C네이버 블로그에서 이전

Convex hull 알고리즘은 아래에서 찾을 수 있다.

Posted by helloktk

댓글을 달아 주세요

Chain Hull

Computational Geometry 2012. 9. 16. 19:28
728x90

참고: http://cm.bell-labs.com/who/clarkson/2dch.c에 구현된 알고리즘을 수정함:

  1. 점집합을 x-좌표의 값(같으면, y좌표의 값) 순으로 정렬한다. 시작점은 hull의 꼭짓점이 된다.
  2. x-좌표 순서대로 현재 점이 hull의 꼭짓점인가를 검사한다: hull의 꼭짓점이 되기 위해서는 직전까지 구성된 hull의 에지 왼편에 있어야 한다. 오른편에 있으면 직전까지 구성된 hull에서 현재 점이 왼편에 놓일 때까지 꼭짓점을 제거한다.
  3. 남는 점들을 1번 과정의 역순으로 정렬한다. 최대 x-좌표, 최소 y-좌표점이 다시 윗 hull의 시작점이 된다.
  4. x-좌표가 작아지는 순서대로 2번 과정을 수행한다.

/* A->B->C가 반시계방향으로 정렬되거나(< 0), 일직선상에 있으면(=) 참이다*/
static int ccw(POINT A, POINT B, POINT C);

더보기
/* A->B->C가 반시계방향으로 정렬되거나(< 0), 일직선상에 있으면(=) 참이다*/
static int ccw(POINT A, POINT B, POINT C) {
    double a = A.x - B.x, b = A.y - B.y, c = C.x - B.x, d = C.y - B.y;
    return (a * d - b * c) <= 0; /* true if points A,B,C counterclockwise */
}
static int cmpl(const void *a, const void *b) {
    POINT *A = (POINT* )a, *B = (POINT* )b;
    int v = (A)->x - (B)->x ;       //lower-hull은 x 증가: y감소 순으로 정렬;
    if (v > 0) return 1; 
    if (v < 0) return -1;
    v = (B)->y - (A)->y ;
    if (v > 0) return 1; 
    if (v < 0) return -1;
    return 0;
}
static int cmph(const void *a, const void *b) {
    return cmpl(b, a);      // upper-hull은 x감소, y-증가 순으로 정렬;
}
static int makeChain(POINT *V, int n, int (*cmp)(const void*, const void*)) {
    qsort(V, n, sizeof(POINT), cmp);
    // 꼭짓점(i(>j))가 현재까지 만들어진 hull(s까지)의 에지보다 더 오른쪽에 
    // 있으면(ccw(i,j,j-1)), 왼쪽에 위치할 때까지 기존 hull을 줄인다. 
    // 이 점을 hull에 추가.
    int s = 1;
    for (int i = 2; i < n; i++) {
        int j = s;
        while (j >= 1 && ccw(V[i], V[j], V[j - 1])) j--; // ccw의 =때문에 일직선상의 
                                                         // 마지막 점만 들어감;
        s = j + 1; // j = last index of hull;
        // (j + 1)에 새 hull-point을 추가;
        POINT t = V[i]; V[i] = V[s]; V[s] = t;         // swap: V[i] <-> V[j + 1]
    }
    return s;
}
int chainHull(std::vector<POINT>& P, std::vector<POINT>& chull) {
    int n = P.size();
    if (n < 3) return 0;
    // copy to working array;
    chull.resize(n + 1);
    for (int i = 0; i < n; i++) chull[i] = P[i];
    /* make lower hull */
    int u = makeChain(&chull[0], n, cmpl);      
    /* make upper hull */   
    chull[n] = chull[0];
    u += makeChain(&chull[u], n - u + 1, cmph); 
    chull.resize(u);
    return u;
};

 

 

 

 

'Computational Geometry' 카테고리의 다른 글

Point in Polygon  (2) 2020.12.14
Incremental Delaunay Triangulation  (1) 2020.12.01
Chain Hull  (2) 2012.09.16
Quick Hull  (2) 2012.09.16
Monotone Cubic Interpolation  (0) 2012.09.09
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2020.12.30 15:10 신고  댓글주소  수정/삭제  댓글쓰기

    C++ STL에는 계산 기하학에서 유용하게 활용할 수 있는 <complex> 클래스가 있어요!
    어려워질 수록 struct로 정의하는 것이 나아보이기는 하지만,
    대부분의 경우에 간단하게 구현할 수 있어 추천드립니다.

  2. hgmhc 2020.12.30 15:11 신고  댓글주소  수정/삭제  댓글쓰기

    chain hull을 해주셨으니...
    rotating calipers를 막연히 기다려보겠습니다...

Quick Hull

Computational Geometry 2012. 9. 16. 02:07
728x90

파란선: 직전 convex hull을 표시함

Input = a set S of n points 
    Assume that there are at least 2 points in the input set S of points
QuickHull (S) { 
    // Find convex hull from the set S of n points
    Convex Hull := {} 
    가장 왼쪽(A)과 오른쪽(B)에 있는 두 점을 찾아 convec hull에 추가;
    선분 AB을 기준으로 나머지 (n-2)개 점들을 두 그룹 S1과 S2로 나눔;
        S1 = 선분 AB의 오른쪽에 놓인 점; 
        S2 = 선분 BA의 오른쪽에 놓인 점;
    FindHull (S1, A, B) 
    FindHull (S2, B, A) 
};
FindHull (Sk, P, Q) { 
   If Sk has no point, 
        then  return. 
   선분 PQ에서 가장 멀리 떨어진 점(C)을 찾아 convex hull에서 P와 Q의 사이에 추가;
   세점 P, C, Q는 나머지 Sk의 나머지 점들을 세 그룹S0, S1, and S2 
        S0 = 삼각형 PCQ에 포함되는 점;
        S1 = 선분 PC의 오른쪽에 놓인 점;
        S2 = 선분 CQ의 오른쪽에 놓인 점; 
    FindHull(S1, P, C) 
    FindHull(S2, C, Q) 
}
Output = convex hull

참고: http://www.cse.yorku.ca/~aaw/Hang/quick_hull/Algorithm.html

 

// P가 line(A, B)의 완전히 왼쪽에 있는가?

static bool leftSide(CPoint P, CPoint A, CPoint B) ;

더보기
static bool leftSide(CPoint P, CPoint A, CPoint B) {
   int ax = P.x - A.x,  ay = P.y - A.y;
   int bx = B.x - A.x,  by = B.y - A.y;
   return (ax * by - ay * bx) >= 0;
};

double dist2ToLine(CPoint P, CPoint A, CPoint B);

더보기
double dist2ToLine(CPoint P, CPoint A, CPoint B) {
    double dx = (B.x - A.x), dy = (B.y - A.y);
    double den = (dx * dx + dy * dy);
    if (den == 0.) return 0; //degenerate;
    double du = (P.x - A.x), dv = (P.y - A.y);
    double dist = (du * dy - dv * dx);
    return dist * dist / den; //distance-square!
};

#define SWAP_POINT(A, B)

static void findHull(CPoint *S, int n, CPoint A, CPoint B, std::vector <CPoint>& hull) ;

더보기
static void findHull(CPoint *S, int n, CPoint A, CPoint B, std::vector<CPoint>& hull) {
    if (n < 1) return ;
    //Find furtherst point 
    double maxdist = 0;
    for (int i = 0; i < n; i++) {
        double dist = dist2ToLine(S[i], A, B);
        if (dist > maxdist) {
            SWAP_POINT(S[0], S[i]);
            maxdist = dist;
        }
    }
    if (maxdist == 0.) return ; //collinear case;
    hull.push_back(S[0]);
    int j = 1;
    for (int i = 1; i < n; i++) {
        if (leftSide(S[i], A, S[0])) {
            SWAP_POINT(S[i], S[j]);
            j++;
        };
    };
    int k = j;
    for (int i = j; i < n; i++) {
        if (leftSide(S[i], S[0], B)) {
            SWAP_POINT(S[i], S[k]);
            k++;
        }
    };
    //1,...,j-1;
    findHull(&S[1], j-1, A, S[0], hull);
    //j,...,k-1;
    findHull(&S[j], k-j, S[0], B, hull);
};

void findMaxMin(std::vector<CPoint>& pts) ;

더보기
void findMaxMin(std::vector<CPoint>& pts) {
    int minx = pts[0].x, maxx = minx;
    int minid = 0, maxid = 0;
    for (int i = 0; i < pts.size(); i++) {
        if (pts[i].x > maxx)
        if (pts[i].x < minx)
    }
    SWAP_POINT(pts[0], pts[minid]); if (maxid == 0) maxid = minid;
    SWAP_POINT(pts[1], pts[maxid]);
};

void hullToPolyline(std::vector<CPoint>& hull);

void quickHull(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    if (pts.size() < 3) return;
    hull.clear();
    //Find left and right most points, say A & B, and add A & B to convex hull ;
    findMaxMin(pts);
    hull.push_back(pts[0]);
    hull.push_back(pts[1]);
    int j = 2;
    for (int i = 2; i < pts.size(); i++) {
        if (leftSide(pts[i], pts[0], pts[1])) {
            SWAP_POINT(pts[i], pts[j]);
            j++;
        }
    }
    //2,3,...,j-1;
    findHull(&pts[2], j-2, pts[0], pts[1], hull);
    //j,j+1,...,pts.size()-1;
    if (j < pts.size()) // in order to avoid dereferencing &pts[pts.size()];
        findHull(&pts[j], pts.size()-j, pts[1], pts[0], hull);
    hullToPolyline(hull);
};

출력 hull은 단순히 convex hull을 구성하는 정렬이 안된 점들만 주므로, hull의 에지를 구하고 싶으면 추가적인 수정이 필요함.

더보기
static int cmph(const void *a, const void *b) {
    CPoint *A = (CPoint *)a , *B = (CPoint *)b ;
    int v = (A->x - B->x) ;
    if (v > 0 ) return 1;
    if (v < 0 ) return -1;
    v = B->y - A->y ;
    if (v > 0 ) return 1;
    if (v < 0 ) return -1;
    return 0;
};
static int cmpl(const void * a, const void *b) {
    return cmph(b, a);
};
void hullToPolyline(std::vector<CPoint>& hull) {
    CPoint A = hull[0];
    CPoint B = hull[1];
    // 가장 멀리 떨어진 두 점(hull[0], hull[1])을 연결하는 직선을 기준으로 프로젝션을 구하여서 순차적으로 
    int k = 2;
    for (int i = 2; i < hull.size(); i++) {
        if (leftSide(hull[i], A, B)) {
            SWAP_POINT(hull[i], hull[k]); k++;
        };
    };
    // k-1; last index of hull left side of line(A,B);
    // upper part reordering: 
    qsort(&hull[0], k, sizeof(CPoint), cmph);
    //lower part reordering;
    if (k < hull.size())
        qsort(&hull[k], hull.size()-k, sizeof(CPoint), cmpl);
    }
};

또한 입력점들의 순서를 그대로 유지하고 싶으면, double pointer를 이용하거나 복사복은 이용하여야 한다.

 

'Computational Geometry' 카테고리의 다른 글

Incremental Delaunay Triangulation  (1) 2020.12.01
Chain Hull  (2) 2012.09.16
Quick Hull  (2) 2012.09.16
Monotone Cubic Interpolation  (0) 2012.09.09
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Curvature of Curve  (0) 2012.08.07
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2020.12.30 15:12 신고  댓글주소  수정/삭제  댓글쓰기

    quick hull이 평균 O(n log n)이였던가요?

728x90

Convex hull은 집합으로 주어진 점이나 영역을 포함하는 가장 작은 볼록 집합이다. 알려진 알고리즘으로는 

1. Jarvis March (Gift Wrap)

2. Graham scan

3. Quick hull (Chain hull)

4. Monotone chain

5.......

여기서는 소개하는 convex hull 알고리즘은 가장 기초적인 알고리즘이다. 두 점에 의해서 만들어지는 선분이 convex hull의 일부이면 나머지 점들은 항상 왼편(설정에 따라 오른편으로 해도 된다)에 있게 된다는 사실을 이용한다. N개의 점이 있을 떄 가능한 선분의 수는 N(N-1) (방향까지 고려함)이고, 나머지 점에 대해서 테스트를 해야 하므로 총 비교 횟수는 N(N-1)(N-2) ~ O(N^3) 정도이다. 알고리즘의 구현이 간단하고, 단순 비교연산이므로 매우 robust해서 점집합의 크기가 작을 때 효과적이다.

 

static int ccw(CPoint A, CPoint B, CPoint P) ; // cross(AB, AP) > 0 if <ABP> is ccw;

더보기
static int ccw(CPoint A, CPoint B, CPoint P) { //cross(AB, AP) > 0 if <ABP> is ccw;
    return ((B.x - A.x) * (P.y - A.y) - (B.y - A.y) * (P.x - A.x));
}
int BruteForceConvexHull(std::vector<CPoint>& pts, std::vector<CPoint>& hull) {
    // O(n^3) ;
    if (pts.size() < 3) return 0;
    hull.clear() ;    
    std::vector<int> head; //hull을 구성하는 에지의 head index;
    std::vector<int> tail; //hull을 구성하는 에지의 tail index;
    for (int i = 0; i < pts.size(); i++) {
        for (int j = 0; j < pts.size(); j++) {
            if (pts[j] == pts[i]) continue;     //에지는 중복점이 아닌 경우만(점을 비교해야 함)
            bool onHull = true;
            for (int k = 0; k < pts.size(); k++) {
                if (pts[k] == pts[i] || pts[k] == pts[j]) continue;
                if (ccw(pts[i], pts[j], pts[k]) < 0) {  // cw;
                    onHull = false;                     // (i,j) is not an edge;
                    break ;
                }
            }
            if (onHull) { //(i->j) edge is on the hull;
                tail.push_back(i); head.push_back(j);
            }
        }
    }
    if (head.size() < 1) return 0;
    // hull을 구성하는 에지정보에서 단순폴리곤 만들기(bug 수정);
    int currIdx = 0;
    hull.push_back(pts[tail[currIdx]]);
    hull.push_back(pts[head[currIdx]]);
    while (head[currIdx] != tail[0]) {
        for (int j = 0; j < head.size(); j++)
            if (tail[j] == head[currIdx]) {
                currIdx = j;
                hull.push_back(pts[head[currIdx]]);
                break;
            }
    }
    return hull.size();
};

convex hull의 edge가 중간에 (collinear) 점을 포함하지 않도록 하기 위해서는, ccw == 0일 때 ABC가 접혀있는지를 체크하면 된다.

'Computational Geometry' 카테고리의 다른 글

Quick Hull  (2) 2012.09.16
Monotone Cubic Interpolation  (0) 2012.09.09
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Curvature of Curve  (0) 2012.08.07
Douglas-Peucker Algorithm  (0) 2012.05.03
Inside Quadrilateral, Inside Triangle  (0) 2012.01.18
Posted by helloktk

댓글을 달아 주세요

728x90

It is well known that the convex hull of a set of n points in the (Euclidean) plane can be found by an algorithm having worst-case complexity O(n log n). In this note we give a short linear algorithm for finding the convex hull in the, case that the (ordered) set of points from the vertices of a simple (i.e., non-self-intersecting) polygon.

ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/81/887/CS-TR-81-887.pdf

http://deepblue.lib.umich.edu/bitstream/2027.42/7588/5/bam3301.0001.001.pdf
We present space-efficient algorithms for computing the convex hull of a simple polygonal line in-place, in linear time. It turns out that the problem is as hard as stable partition, i.e., if there were a truly simple solution then stable partition would also have a truly simple solution, and vice versa. Nevertheless, we present a simple self-contained solution that uses O(log n) space, and indicate how to improve it to O(1) space with the same techniques used for stable partition. If the points inside the convex hull can be discarded, then there is a truly simple solution that uses a single call to stable partition, and even that call can be spared if only extreme points are desired (and not their order). If the polygonal line is closed, then the problem admits a very simple solution which does not call for stable partitioning at all.

http://springerlink.metapress.com/content/5clwpcjenldf5a78/fulltext.pdf

Posted by helloktk

댓글을 달아 주세요