평면 상에 주어진 점들의 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();
};

**네이버 블로그 이전;

 

 

728x90

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

Minimum Enclosing Circle  (0) 2021.03.01
Creating Simple Polygons  (0) 2021.01.25
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Posted by helloktk
,

단순 다각형(simple polygon)의 무게중심(center of gravity or center of mass)은 다각형을 균일한 밀도의 판으로 생각했을 때 판의 무게중심과 같다. 가장 단순한 다각형인 삼각형의 무게중심은 세 꼭짓점의 산술평균으로 표현된다.

$$ \text{CoG} = \frac{1}{3} ({{\bf P} + {\bf Q} + {\bf R}}).$$

증명: 삼각형의 한 변 PQ에 나란한 띠로 삼각형을 분할하자. 그러면 각 띠의 무게중심은 띠의 기하학적 중심이므로 꼭지점 R와 변 PQ의 중심을 연결한 선분 RA상에 있어야 한다. 그리고 띠의 무게중심은 그 점에 모든 질량이 뭉친 것으로 생각할 수 있으므로 전체 삼각형의 무게중심은 선분 RA상의 어느 지점에 있어야 한다. 마찬가지 논리로  삼각형을 선분 QR에 나란한 띠로 분할나누면 이 띠들의 무게중심은 선분 PB상에 있음을 알 수 있다. 따라서 삼각형의 무게중심도 PB상의 한점이어야 하므로 선분 PB와 선분 RA의 교점이 무게중심이 된다. 

다각형은 삼각형으로 분할되므로 이 분할된 삼각형의 무게중심을 이용하면 쉽게 계산할 수 있다. 분할된  삼각형의 무게중심을 면적으로 가중치를 준 평균값이 다각형의 무게중심이 된다. 

 

실제 계산에서는 다각형을 삼각 분할하지 않고도 간단한 방법에 의해서 무게중심을 구할 수 있다. 원점과 다각형의 각 변의 꼭짓점을 이용해서 삼각형들을 구성하면 원래의 다각형을 겹치게(원점이 내부에 있으면 겹침이 없다) 분할할 수 있다. 분할된 삼각형으로 무게중심을 구할 때 겹치는 영역의 기여를 제거해야 한다. 그런데 다각형 밖의 영역을 분할하는 삼각형은 다각형 내부를 분할하는 삼각형과는 다른 orientation을 가지게 된다. 삼각형의 면적은 한 꼭짓점을 공유하는 두 변의 외적에 비례하므로, 반대의 orientation을 갖는 삼각형은 자동으로 반대 부호의 면적을 가지게 된다. 따라서 분할된 삼각형 면적 가중치를 외적으로 주어서 무게중심을 구하면 겹치는 영역이 자동으로 상쇄되는 효과를 얻을 수 있다.

노란색 영역에 있는 분할 삼각형은 음의 면적

$$\begin{align} \text{CoG} &= \frac{1}{\text{다각형 면적}} \sum (\text{삼각형 CoG})  (\text{면적})  \\ &= \frac{1}{\text{다각형 면적}} \sum \frac{1}{3} \left( {\bf V}_i + {\bf V}_{i+1} + {\bf O}\right ) \frac{ {\bf V}_{i} \times {\bf V}_{i+1} }{2} \\ &= \frac{1}{3}\frac{1}{   \text{다각형 면적} }\sum ( {\bf V}_{i} + {\bf V}_{i+1}) \frac{{\bf V}_{i} \times {\bf V}_{i+1}}{2} \end{align}$$ 

다각형의 면적($=\sum \frac{1}{2}({\bf V}_i \times {\bf V}_{i+1})$)을 구할 때 삼각형과 동일하게 orientation에 따라 부호를 포함하도록 설정하면 다각형의 면적 부호가 삼각형의 면적 부호로 상쇄되므로 다각형의 orientation에 무관하게 성립하는 공식이 된다.

CPoint polygon_centroid(CPoint V[], int N) {
    double cx = 0, cy = 0, area2 = 0;
    for(int i = 0, j = N - 1; i < N; j = i++) {
        double tri_area2 = V[i].x * V[j].y - V[i].y * V[j].x; // area * 2;
        cx += (V[i].x + V[j].x) * tri_area2;
        cy += (V[i].x + V[j].x) * tri_area2;
        area2 += tri_area2;                                   //total area * 2
    }
    cx /= 3 * area2;
    cy /= 3 * area2;
    return CPoint(int(cx + 0.5), int(cy + 0.5))
};

** 네이버 블로그에서 수정 이전;

728x90

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

Creating Simple Polygons  (0) 2021.01.25
단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Posted by helloktk
,

평면상의 다각형(모서리의 교차가 없는 단순 다각형)의 면적을 구하는 것은 단순하지 않을 것처럼 보이지만 계산식은 무척이나 간단하게 주어진다. 기본적인 아이디어는 다각형에 임의의 점을 찍으면 이 점과 이웃한 두 개의 꼭짓점으로 형성이 되는 삼각형의 합으로 다각형을 분할할 수 있다. 분할된 삼각형의 면적을 구하여 합산하면 다각형의 면적을 구할 수 있다.

 

세 점 ${\bf P, Q, R}$(이 순서대로 반시계방향으로 배열)이 만드는 삼각형의 면적은

$$\text {삼각형의 면적}=  \frac{1}{2} ({\bf R} - {\bf  Q})  \times ( {\bf P} - {\bf Q}); \quad  (\Rightarrow \text {시계 방향이면 면적이 음수})$$

로 주어지므로, 꼭짓점이 ${\bf P}_0(x_0, y_0), {\bf P}_1(x_1, y_1),....$(반시계 방향)으로 주어지는 $N$각형의 면적은 아래와 같이 주어진다. 

$$\begin{align} \text{다각형 면적} &= \sum \text{각 삼각형의 면적} \\ &= \frac{1}{2}\sum ({\bf P}_{i+1}-{\bf Q})\times ({\bf P}_{i}-{\bf Q})\quad \quad ({\bf Q}\text{는  임의의 점})\\ &= \frac{1}{2} \sum\left(  {\bf P}_{i+1} \times {\bf P}_{i} - {\bf P}_{i+1}\times {\bf Q} + {\bf P}_{i} \times {\bf Q}\right) \\ &=\frac{1}{2} \sum {\bf P}_{i+1} \times {\bf P}_{i} \\ &= \frac{1}{2}\sum \left( x_{i+1} y_{i} -x_{i}y_{i+1} \right) \end{align}$$

이 결과는 $\bf Q$에 무관하다. 다각형의 꼭짓점이 시계 방향으로 정렬이 된 경우는 면적이 음수로 나온다(윈도우 DC는 위-아래가 역전되어 있으므로 orientation이 반대로 보인다). 그리고 이 공식은 단순 다각형에만 적용이 되고 모서리의 교차가 있는 경우에는 적용이 되지 않는다.

double simple_polygon_area2D(POINT point[], int N) {
    double area = 0;
    for (int i = 0, j = N - 1; i < N; j = i++) 
        area += (point[i].x * point[j].y) - (point[i].y * point[j].x);
    area /= 2;
    // return area;  // signed area;
    return area < 0 ? -area: area;
}

**네이버 블로그 수정 이전;

728x90

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

단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Catmull-Rom Spline  (0) 2020.12.07
Posted by helloktk
,

세 점 $A$, $B$, $C$가 만드는 삼각형의 외접원의 반지름을 구해보자. 세 변을 나타내는 벡터는

$$\vec{a}= \vec{OC}-\vec{OB}, \quad \vec{b}=\vec{OA}-\vec{OC},\quad\vec{c}=\vec{OB}-\vec{OA}.$$

세 점이 일직선 위에 놓이지 않기 위해서는 

$$\vec{a}\times \vec{b} \ne  0$$

를 만족해야 한다. 

그림의 삼각형에서 $\angle(BOC)$가 $\angle(BAC)$의 두 배이므로,

$$\sin \theta =\frac{a/2}{R}.$$

이제 삼각형 면적은 세 변의 길이 ($a, b, c$)와 외접원의 반지름($R$)으로 표현된다:

$$\text{area}=\frac{1}{2}bc \sin \theta ~~\longrightarrow~~ \frac{abc}{4R}.$$

또한, 삼각형의 면적을 세 점이 일직선 위에 있는가를 테스트하는 식으로 표현하면 

$$\text{area}=\frac{1}{2} | \vec{a}\times \vec{b}|$$

처럼 쓸 수 있다. 따라서 삼각형의 외접원의 반지름은 아래처럼 주어진다.

$$R = \frac{abc}{4\cdot \text{area}} = \frac{| \vec{a}||\vec{b} ||\vec{c} |}{2|\vec{a}\times \vec{b}| }.$$

#define SQR(x) ((x) * (x))
double circumradius(CPoint A, CPoint B, CPoint C) {
    double ax = C.x - B.x, ay = C.y - B.y;
    double bx = A.x - C.x, by = A.y - C.y;
    double crossab = ax * by - ay * bx;
    if (crossab != 0) { 
        double a = sqrt(SQR(ax) + SQR(ay));
        double b = sqrt(SQR(bx) + SQR(by)); 
        double cx = B.x - A.x, cy = B.y - A.y;       
        double c = sqrt(SQR(cx) + SQR(cy));
        return (0.5 * a * b * c / fabs(crossab));
    } else 
        return 0;   //collinear
};
728x90
Posted by helloktk
,

Chain Hull

Computational Geometry 2012. 9. 16. 19:28

참고: 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) {
    int 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>& pts, std::vector<POINT>& chull) {
    int n = pts.size();
    if (n < 3) return 0;
    // copy to working array;
    chull.resize(n + 1);
    for (int i = 0; i < n; i++) chull[i] = pts[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;
};

 

 

 

 

728x90

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

Catmull-Rom Spline  (0) 2020.12.07
Incremental Delaunay Triangulation  (1) 2020.12.01
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
,