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
,

Quick Hull

Computational Geometry 2012. 9. 16. 02:07

파란선: 직전 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 = pts.size(); i-->1;) {
        if (pts[i].x > maxx) maxid = i, maxx = pts[i].x;
        if (pts[i].x < minx) minid = i, minx = pts[i].x;
    }
    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를 이용하거나 복사복은 이용하여야 한다.

 

 
728x90

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

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

3차의 Hermite spline을 쓰는 interpolation에서 입력 데이터 점에서의 기울기 값을 수정하면 입력 데이터의 monotonic 한 성질을 그대로 유지하는 보간을 할 수 있다. 주어진 데이터가 $\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\}$이고, $x/y$-좌표가 증가하는 순으로 정렬되었다고 하자. 구간 $[x_k, x_{k+1} ]$의 시작과 끝에서 값 ($p_0 = y_k $, $p_1=y_{k+1}$)과 기울기 ($m_0, m_1$) 가 주어지면 Cubic Hermite spline 보간 함수는 다음과 같이 매개변수 $t \in [0,1]$의 함수로 표현된다: 

$$f(t) = p_0 h_{00}(t) + m_1 h_{10}(t) + p_1 h_{01}(t) + m_1 h_{11}(t),            \quad \quad t= \frac {x-x_k}{x_{k+1}-x_k} $$

$$f(0) = p_0,\quad f'(0) = m_0,\quad f(1) = p_1,\quad f'(1)=m_1$$

Hermite 기저 함수의 정의는 아래와 같다:

$$ h_{00}(t) = 2t^3 - 3t^2+1, \quad h_{10}(t) = t^3 - 2t^2+1,\quad h_{01}(t)=-2t^3+3t^2,\quad h_{11}(t)=t^3-t^2$$

 

보간된 spline이 단조성(monotonicity)을 유지하도록 각 입력점에서 기울기 값($m_{k}$)을 조절해주어야 한다. 그럼 기울기 $m_{k}$는 어떻게 설정해야 인접 구간의 보간 함수와 smooth 하게 연결될 수 있을까? 입력점의 각 위치에서 기울기는 각각 왼쪽과 오른쪽 인접 입력점과의 평균 변화율($\delta_{k-1}, \delta_{k}$)을 구한 후 다시 이들의 평균값을 기울기로 사용하면 된다. (양 끝점에서는 오른쪽이나 왼쪽 구간의 평균 변화율로 대체)

$$ m_k = \frac {1}{2}(\delta_{k-1}+\delta_{k} )= \frac {1}{2}\left(\frac { y_k-y_{k-1}}{x_k - x_{k-1}}+\frac {y_{k+1}-y_k}{x_{k+1} + x_k}\right)$$

이것만으로는 단조성을 보증을 할 수 없으므로 다음 조건을 추가로 부여해야 됨이 관련 논문에서 증명이 되었다. 

$$1. \text{For}~~\delta_{k} = 0 \quad\Longrightarrow\quad m_k =m_{k+1} = 0$$

$$2. \text{For}~~\epsilon=\sqrt {\left(\frac {m_k }{\delta_k }\right)^2 + \left(\frac {m_{k+1}}{\delta_{k} }\right)^2 } > 3\\ \quad\Longrightarrow\quad m_{k} \leftarrow3\frac {m_{k}}{\epsilon}, ~m_{k+1} \leftarrow 3\frac {m_{k+1}}{\epsilon}$$

 

아래의 코드는 위키피디아에 제시된 알고리즘을 구현한 것이다

(입력 데이터가 한 방향으로 monotonic 할 때만 동작함)

참고: http://en.wikipedia.org/wiki/Monotone_cubic_interpolation

void estimateTangents(std::vector <double>& x, std::vector <double>& y, std::vector<double>& m);

void monotoneCubicSplineInterp(std::vector<CPoint>& points, std::vector<CPoint>& curve) {
    curve.clear();
    if (points.size() < 2) return ;
    //if (!isMonotonicPolyline(points)) return ;
    std::vector<double> xx(points.size()), yy(points.size());
    for (int i = points.size(); i-->0;) {
        xx[i] = points[i].x ;
        yy[i] = points[i].y;
    };
    //tangent table;
    std::vector<double> m(points.size(), 0);
    estimateTangents(xx, yy, m);
    //hermite spline;
#define STEPS 12
    curve.push_back(points[0]); //add first point;
    for (int i = 0; i < points.size() - 1; i++) {
        double h = xx[i + 1] - xx[i];
        for (int k = 1; k <= STEPS; k++) {
            double t = double(k) / STEPS;         // 0 < t <= 1;
            double t2 = t * t ;
            double t3 = t2 * t ;
            double h00 =  2. * t3 - 3. * t2     + 1.;
            double h10 =       t3 - 2. * t2 + t;
            double h01 = -2. * t3 + 3. * t2;
            double h11 =       t3      - t2;
            double yf = h00 * yy[i] + h10 * h * m[i] + h01 * yy[i + 1] + h11 * h * m[i + 1];
            curve.push_back(CPoint(int(xx[i] + t * h + 0.5), int(yf + 0.5)));
        }
    }
};

Natural cubic spline과의 비교:

 

 

 

 

 

 

 

 
 
728x90

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

Chain Hull  (2) 2012.09.16
Quick Hull  (2) 2012.09.16
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Curvature of a Curve  (0) 2012.08.07
Douglas-Peucker Algorithm  (0) 2012.05.03
Posted by helloktk
,

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가 접혀있는지를 체크하면 된다.

728x90

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

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