728x90

볼록 다각형(convex polygon)에서 한 방향으로(예를 들면 y-축방향) 가장 멀리 떨어진 두 꼭지점(antipodal points) 통과하는 평행한 두 직선을 생각해보자. 평행 직선의 방향은 이 경우에는 x-축(y-축과 수직)과 나란해야 한다 . 이제 두 평행 직선의 방향을 각각의 꼭지점을 축으로 해서 조금씩 회전시키면, 둘 중 하나는 볼록 다각형에 먼저 접하게 된다. 이 접하는 변을 기준으로 다시 antipodal point을  찾아 (당연히 접하지 않는 평행직선이 통과하는 꼭지점이 된다) 새로운 평행 직선을 구성한 후 동일한 과정을 반복한다. 이 과정을 N회 반복하면 볼록 다각형의 모든 변에 대해서 변을 공유하는 직선과 그 직선의 antipodal point를 찾을 수 있다. 

이 rotating calipers를 이용하면 convex hull의 bounding box, convex hull의 diameter, convex polygon 사이 거리 등의 계산에 응용할 수 있다.

double COS(CPoint A, CPoint B) {             // directional cosine;
    double norm2A = A.x * A.x + A.y * A.y;
    double norm2B = B.x * B.x + B.y * B.y;
    return (A.x * B.x + A.y * B.y) / sqrt(norm2A * norm2B);
}
int ccw(CPoint A, CPoint B, CPoint C) {     // cross(AB, BC); 반시계방향 > 0;
    return (B.x - A.x) * (C.y - B.y) - (B.y - A.y) * (C.x - B.x) ;
}
// find antipodal points along y-direction;
bool initCalipers(std::vector<CPoint> &Q, 
                 int &minid, double &minCos,
                 int &maxid, double &maxCos) {
    int N = Q.size();
    // ccw check (중복 // 일직선을 고려);
    int orient = 0;
    for (int i = 0; i < N; i++) {
        orient = ccw(Q[(i - 1 + N) % N], Q[i], Q[(i + 1) % N]);
        if (orient != 0) break;
    }
    if (orient == 0) return false;

    minid = maxid = 0;
    for (int i = 1; i < N; i++) {
        if (Q[i].y < Q[minid].y) minid = i;
        if (Q[i].y > Q[maxid].y) maxid = i;
    }

    CPoint n0 = CPoint(1, 0);
    if (orient < 0) n0 = -n0;  // convex hull이 CW인 경우;
    CPoint dir0 = Q[(minid + 1) % N] - Q[minid];
    CPoint dir1 = Q[(maxid + 1) % N] - Q[maxid];
    minCos = COS(dir0, n0);
    maxCos = COS(dir1, -n0);
    return true;
}
int rotatingCalipers(std::vector<CPoint> &Q,
                    int &i0, double &dirCos0,
                    int &i1, double &dirCos1) {
    int N = Q.size();
    int i0next = (i0 + 1) % N ;
    int i1next = (i1 + 1) % N ;
    CPoint dir0 = Q[i0next] - Q[i0];
    CPoint dir1 = Q[i1next] - Q[i1];
    if (dirCos0 > dirCos1) {          // i0을 포함하는 평행선이 다음 edge와 이루는 각이 작음;
        CPoint ndir0 = Q[(i0next + 1) % N] - Q[i0next]; // 각을 재는 기준선은 현재 다다음 edge;
        dirCos0 = COS(ndir0, dir0);    
        dirCos1 = COS(dir1, -dir0);    
        i0 = i0next;                  // i0는 다음 꼭지점으로 이동; i1 = antipodal point;
        return 0;
    } else {                          // i1을 포함하는 평행선이 다음 edge와 이루는 각이 작음
        CPoint ndir1 = Q[(i1next + 1) % N] - Q[i1next];
        dirCos1 = COS(ndir1, dir1);   
        dirCos0 = COS(dir0, -dir1);   
        i1 = i1next;                  // i1은 다음 꼭지점으로 이동; i0 = antipodal point;
        return 1;
    }
}
void runCalipers(std::vector<CPoint> &Q) {
    int i0, i1;
    double dirCos0, dirCos1;
    bool done = false;
    if (!initCalipers(Q, i0, dirCos0, i1, dirCos1)) return ;  // fails
    while (!done) {
    	int i0old = i0, i1old = i1;
        if (rotatingCalipers(Q, i0, dirCos0, i1, dirCos1) == 0) {
            // Q[i1] = antipodal point;
            // do something useful...
        } else {
            // Q[i0] = antipodal point;
            // do something useful...
        }
    }
};

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

Bezier Curve Approximation of a Circle  (0) 2021.04.10
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
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2021.04.04 11:56 신고  댓글주소  수정/삭제  댓글쓰기

    혹시 로테이팅 캘리퍼스의 정당성 증명 자료가 있으신가요?

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

$$\underset{q \in P}{\text{max}}~|q - c| \le (1+\epsilon)R_P$$

$\epsilon=0.1$ (red=exact MEB)인 경우;

void ApproxMinBall(std::vector<CPoint> &pts, double epsilon, CfPt &center, double &radius) {
    double epsilonsq = epsilon * epsilon;
    int maxiters = int( 1 / epsilonsq );
    double maxdist2 = 0;
    int id = rand() % pts.size();
    center.x = pts[id].x; 
    center.y = pts[id].y;
    for (int iter = 1; iter < maxiters; iter++) { 
        // find the furthest point from (cx, cy);
        maxdist2 = 0;
        id = 0;
        for (int i = pts.size() - 1; i >= 0; --i) {
            double dx = pts[i].x - center.x, dy = pts[i].y - center.y;
            double dist2 = dx * dx + dy * dy;
            if (dist2 > maxdist2) {
                maxdist2 = dist2;
                id = i;
            }
        }
        // update center;
        center.x += (pts[id].x - center.x) / iter;
        center.y += (pts[id].y - center.y) / iter;
    }
    radius = sqrt(maxdist2);
}

kipl.tistory.com/300

 

Minimum Enclosing Circle

평면 위에서 주어진 점집합을 포함하는 가장 작은 원을 찾는 문제는 가장 단순하게는 임의의 두 점이 지름이 되는 원과 임의의 세 점이 만드는 외접원을 모두 조사하여 찾으면 된다. 따라서 O(n^4

kipl.tistory.com

 

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

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
단순 다각형 만들기  (0) 2021.01.25
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2021.03.20 21:51 신고  댓글주소  수정/삭제  댓글쓰기

    이게 말로만 듣던 최소 외접원인가요?

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

각 전경 픽셀에 대해서 모든 배경 픽셀까지 거리를 구한 다음 최솟값을 할당하면 된다. 픽셀 수가 n = w * h 개면 time complexity는 O(n^2)이다 (이미지 폭이나 높이의 4승이므로 연산량이 상당하다.) linear time Euclidean distance transform은 kipl.tistory.com/268.

ListPlot3D[ Reverse@ImageData[pic], PlotRange -> Full]

mathematica를 이용한 3D plot

int brute_force_edt(double *img, int w, int h) {
    int n = w * h;
    int *list = new int [n];
    int fg_cnt = 0;
    int bg_ptr = n - 1;               //배경 픽셀 위치는 끝에서 역으로 채움;
    //전경과 배경 픽셀 위치를 분리;
    for (int i = 0; i < n; ++i)
        if (img[i]) {                 // foreground;
            list[fg_cnt++] = i;
            img[i] = DBL_MAX;
        } else                       // background;
            list[bg_ptr--] = i;

    for (int i = 0; i < fg_cnt; ++i) {     // 전경 픽셀;
        int xi = list[i] % w, yi = list[i] / w;
        for (int j = fg_cnt; j < n; ++j) { // 배경 픽셀까지 거리를 계산해서 최소값을 할당;
            int xj = list[j] % w, yj = list[j] / w;
            double dx  = xi - xj, dy = yi - yj;  //
            double dst = dx * dx + dy * dy;
            if (img[list[i]] > dst)
                img[list[i]] = dst;
        }
    }
    delete [] list;
    return fg_cnt;
}

'Image Recognition > Fundamental' 카테고리의 다른 글

Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Posted by helloktk

댓글을 달아 주세요

728x90

평균이 $\lambda$인 Poisson 분포를 가지는 random number를 만드는 레시피는 [0,1]에서 추출된 균일한 random number 수열이 $u_1, u_2, u_3,...$로 주어질 때, $u_1  u_2 u_3....u_{k+1}\le  e^{-\lambda}$을 처음 만족하는 $k$를 찾으면 됨이 Knuth에 의해서 알려졌다.

int PoissonValue(double val) { 
    double L = exp(-val);
    int k = 0;
    double p = 1;
    do {
        k++;
        // generate uniform random number u in [0,1] and let p ← p × u.
        p *= double(rand()) / RAND_MAX;
    } while (p > L);
    return (k - 1);
}

이미지에서 각각의 픽셀 값은 영상신호를 받아들이는 ccd 센서에서 노이즈가 포함된 전기신호의 평균값이 구현된 것으로 볼 수 있다. 따라서 영상에 Poisson 노이즈를 추가하는 과정은 역으로 주어진 픽셀 값을 평균으로 가지는 가능한 분포를 찾으면 된다.(물론 이미지에 무관하게 노이즈를 추가할 수도 있다. 예를 들면 모든 픽셀에 대해서 같은 평균값을 갖는 포아송 노이즈를 주는 경우다)

void AddPoissonNoise(CRaster& raster, CRaster& noised) {
    if (raster.IsEmpty() || raster.GetBPP() != 8) return;
    CSize sz = raster.GetSize();
    noised = raster;
    srand(unsigned(time(0)));
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)noised.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            int a = PoissonValue(*p++);
            *q++ = a > 0xFF ? 0xFF: a;
        }
    }
}

 

'Image Recognition > Fundamental' 카테고리의 다른 글

Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Posted by helloktk

댓글을 달아 주세요

728x90

blue cell: obstacles, black cell = start( or end) 

$$\text{distance measure}:~\text{Manhattan distance}=|x_2 - x_1 | + |y_2 - y_1|$$

Grassfire 알고리즘을 써서 최단경로를 찾는 문제는 DFS로도 구현할 수 있으나 비효율적이고(왜 그런지는 쉽게 알 수 있다) BFS로 구현하는 것이 더 적합하다. 

void grassfire(CPoint q, int **map, int w, int h, int **dist) {
    for (int y = 0; y < h; y++) 
        for (int x = 0; x < w; x++) 
            dist[y][x] = INF;  //unvisited cells;

    std::queue<CPoint> Q;
    dist[q.y][q.x] = 0;     //start( or end) position: distance = 0;
    Q.push(q);
    while (!Q.empty()) {
        CPoint p = Q.front(); Q.pop();
        int distance = dist[p.y][p.x];
        if (p.x > 0     && map[p.y][p.x - 1] == 0 && dist[p.y][p.x - 1] == INF) {
            dist[p.y][p.x - 1] = distance + 1;
            Q.push(CPoint(p.x - 1, p.y)); 
        }
        if (p.x < w - 1 && map[p.y][p.x + 1] == 0 && dist[p.y][p.x + 1] == INF) {
            dist[p.y][p.x + 1] = distance + 1;
            Q.push(CPoint(p.x + 1, p.y));
        }
        if (p.y > 0     && map[p.y - 1][p.x] == 0 && dist[p.y - 1][p.x] == INF) {
            dist[p.y - 1][p.x] = distance + 1;
            Q.push(CPoint(p.x, p.y - 1));
        }
        if (p.y < h - 1 && map[p.y + 1][p.x] == 0 && dist[p.y + 1][p.x] == INF) {
            dist[p.y + 1][p.x] = distance + 1;
            Q.push(CPoint(p.x, p.y + 1));
        }
    }
}

'Image Recognition > Fundamental' 카테고리의 다른 글

Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Posted by helloktk

댓글을 달아 주세요

728x90

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

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;
}

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

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
단순 다각형의 무게중심  (0) 2021.01.24
Posted by helloktk

댓글을 달아 주세요