728x90

Brute-force : O(n^2),

Optimal: O(n log(n))

참고: people.csail.mit.edu/indyk/6.838-old/handouts/lec17.pdf

참고: arxiv.org/pdf/1911.01973.pdf

//  points should be sorted in order of increasing x-component.
//
#define DIST2(p, q) ((p).x-(q).x)*((p).x-(q).x) + ((p).y-(q).y)*((p).y-(q).y)
// 수정: index range -> [idx1, idx2]; ClosestPair(points, 0, points.size()-1, closest);
int ClosestPair(CPoint points[], int idx1, int idx2, int closest[2]) {
    if (idx2 - idx1 >= 3) {
        int lpair[2], rpair[2], min_dist2;
        int idx12 = idx1 + (idx2 - idx1) / 2;                    // half index;
        int dist1 = ClosestPair(points, idx1, idx12, lpair);     // left side;
        int dist2 = ClosestPair(points, idx12 + 1, idx2, rpair); // right side;
        if (dist1 < dist2) {
            closest[0] = lpair[0]; closest[1] = lpair[1];
            min_dist2 = dist1;
        } else {
            closest[0] = rpair[0]; closest[1] = rpair[1];
            min_dist2 = dist2;
        }
        // find points which lies near center strip(2d-strip);
        // Note our distance is the squar of actual distance;
        int width = int(sqrt(double(min_dist2))+0.5);
        int ll = idx1;
        while (points[ll].x < (points[idx12].x - width - 1)) ll++;
        int rr = idx2;
        while (points[rr].x > (points[idx12 + 1].x + width + 1)) rr--;
        for (int i = ll; i < rr; i++) {
            for (int j = i + 1; j <= rr; j++) {
                int dist2 = DIST2(points[i], points[j]);
                if (min_dist2 > dist2) {
                    min_dist2 = dist2;
                    closest[0] = i; closest[1] = j;
                }
            }
        }
        return min_dist2;
    } 
    else if (idx2 == idx1 + 2) {
        return ClosestPair3(points, idx1, closest);
    }
    else if (idx2 == idx1 + 1) {
        closest[0] = idx1; closest[1] = idx2;
        return DIST2(points[idx1], points[idx2]);
    }
    else return INT_MAX;
};

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

Data Fitting with B-Spline Curves  (0) 2021.04.30
Closest Pair of Points  (0) 2021.04.27
DDA Algorithm  (0) 2021.04.25
B-Spline  (0) 2021.04.25
Bezier Curve Smoothing  (0) 2021.04.23
Flatness of Cubic Bezier Curve  (0) 2021.04.23
Posted by helloktk

댓글을 달아 주세요

728x90

n 차 Bezier 곡선은 두 개의 (n - 1) 차 Bezier 곡선의 선형보간으로 표현할 수 있다. Bezier 곡선은 Bernstein 다항식을 이용해서도 표현할 수도 있지만, 높은 찻수의 곡선일 때는 De Casteljau's Algorithm을 이용하는 것이 수치적으로 보다 안정적인 결과를 준다.

// De Casteljau's algorithm; recursive version;
double Bezier(double t, double n, double Q[]){
   if (n == 1) return Q[0];
   else return (1 - t) * Bezier(t, n - 1, &Q[0]) + t * Bezier(t, n - 1, &Q[1]);
}
// De Casteljau's algorithm; non-recursive. calling of Bezier() modifies Q's;
double Bezier(double t, double n, double Q[]){
    for (int k = 1; k < n; k++)
        for (int j = 0; j < (n - k); j++)
            Q[j] = (1 - t) * Q[j] + t * Q[j + 1];
    return Q[0];
}
void BezierCurve(std::vector<CfPt> &cntls,
                 int npts, std::vector<CfPt> &curves) {
    int n = cntls.size();
    std::vector<double> xp(n);
    std::vector<double> yp(n);
    for (int k = 0; k < n; k++) {
        xp[k] = cntls[k].x;
        yp[k] = cntls[k].y;
    }
    curves.resize(npts);    
    for (int i = 0; i < npts; ++i) {
        double t = double(i) / (npts - 1);
        curves[i].x = Bezier(t, n, &xp[0]);
        curves[i].y = Bezier(t, n, &yp[0]);
    }
}

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

Flatness of Cubic Bezier Curve  (0) 2021.04.23
Convexity in Bezier Curve  (0) 2021.04.22
De Casteljau's Algorithm  (0) 2021.04.22
Bezier Curve Arc Length  (0) 2021.04.21
Bezier Curve Approximation of an Ellipse  (0) 2021.04.11
Bezier Curve Approximation of a Circle  (0) 2021.04.10
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  (1) 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

댓글을 달아 주세요

728x90

FFT는 입력 데이터의 개수$(N)$가 2의 지수승으로 주어질 때  $O(N\log N)$의 연산만으로 빠르게 DFT을 수행하는 알고리즘이다. 주어진 $N$개의 data $\{F_0, F_1, ..., F_{N-1}\}$의  DFT $\{a_0, a_1,..., a_{N-1}\}$는

$$ a_k = \sum_{n = 0}^{N-1} F_n \exp\Big( -i \frac{2\pi nk }{N} \Big)=\sum_{n = 0}^{N-1} F_n (W^n )^k , \quad W = \exp\Big( -i \frac{2\pi}{N}\Big)$$ 로 표현된다. $N$이 2의 지수승으로 주어지면 이 합은 $n$이 짝수인 항과 홀수인 항으로 아래처럼 분리할 수 있다:

\begin{align} a_k &= \sum _{n=0}^{N/2-1} F_{2n} (W^{2n})^k+W^{k} \sum_{n=0}^{N/2-1} F_{2n+1}(W^{2n})^k\\ &=\text{(N/2개 even 항에 대한 DFT)} + W^k \times \text{(N/2개 odd 항에 대한 DFT)}\\ &= e_k + W^k o_k \end{align}

이는 원래의 문제가 처음 입력 데이터의 절반만 가지고 처리하는 두 개의 동일한 문제로 줄일 수 있음을 보여준다. 전형적인 divide and conquer 기법을 적용할 수 있는 구조이므로 매우 빠르게 연산을 수행할 수 있다(Cooley-Tukey algorithm). 더불어 주기성 때문에 $a_k$의 계산도 절반만 하면 된다: 

$$a_{k+ N/2} = e_k - W^k o_k, \quad k = 0,1,..., N/2-1.$$

역변환의 경우에는 twiddle factor $W$를 $W^*$로 바꾸면 되므로(전체 normalization factor $\frac{1}{N}$가 덧붙여진다) 코드에서 쉽게 수정할 수 있다: $\tt W_{im}  \rightarrow  - { W}_{im}$. 아래는 FFT를 재귀적으로 구현한 코드이다. 비재귀적 구현은 kipl.tistory.com/22에서 찾을 수 있다.

void split(int N, double* data) {
    double *tmp = new double [N / 2];
    for (int i = 0; i < N / 2; i++)   //copy odd elements to tmp buffers;
        tmp[i] = data[2 * i + 1];
    for (int i = 0; i < N / 2; i++)   // move even elements to lower half;
        data[i] = data[2 * i];            
    for (int i = 0; i < N / 2; i++)   // move odd elements to upper half;
        data[i + N / 2] = tmp[i];     
    delete [] tmp;
}
void fft(int N, double* re, double* im) {
    if (N < 2) return;
    else {
        split(N, re);                             // divide;
        split(N, im);        
        fft(N / 2, &re[    0], &im[    0]);       // conquer;
        fft(N / 2, &re[N / 2], &im[N / 2]);
        for (int k = 0; k < N / 2; k++) {
            double Ere = re[k];
            double Eim = im[k];
            double Ore  = re[k + N / 2];
            double Oim  = im[k + N / 2];
            double Wre =  cos(2. * PI * k / N);   // real part of twiddle factor:W
            double Wim =  -sin(2. * PI * k / N);  // imag part of twiddle factor:W
            double WOre = Wre * Ore - Wim * Oim;  // WO = W * O
            double WOim = Wre * Oim + Wim * Ore;
            re[k] = Ere + WOre;
            im[k] = Eim + WOim;
            re[k + N / 2] = Ere - WOre;
            im[k + N / 2] = Eim - WOim;
        }
    }
    return;
}

 6개의 주파수$\{f=2,5,9,11,21,29\text{Hz}\}$가 섞인 신호에서 [0,1]초 사이에 일정한 간격으로 sampling한  64개의 data을 FFT한 결과 (Nyquist frequency=32Hz);

$$s(t)=\sum_{i=0}^{5} (i+1) \cos(2\pi f_i t) $$

int fft_test_main() {
    const double PI = 4. * atan(1.);
    const int samples = 64;
    double re[samples], im[samples];
    const int nfreqs = 6;
    double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
    for ( int i = 0; i < samples; i++ ) {
        double signal = 0;
        for ( int k = 0; k < nfreqs; k++ ) 
            signal += (k + 1) * cos ( 2.0 * PI * freq[k] * i / samples );
        re[i] = signal;
        im[i] = 0.;
    }
    fft ( samples, &re[0], &im[0] );
    for ( int i = 0; i < samples; i ++ ) TRACE ( "%d, %f, %f\n", i, re[i], im[i] );
    return 0;
}

FFT결과:

real part (imaginary part = 0)

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

Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.

이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram(ch)과 평균을 계산할 수 있는 1차 culmuative histogram(cxh)을 입력으로 사용한다. 

$$ th= \text {argmax} \left( \sigma^2_B = \sum_{j=0}^{nclass-1} \omega_j m_j^2 \right)$$

 

* Otsu 알고리즘을 이용한 이미지 이진화 코드: kipl.tistory.com/17

* 좋은 결과를 얻으려면 히스토그램에 적당한 필터를 적용해서 smooth하게 만드는 과정이 필요하다.

// 0 <= start < n;
double histo_partition(int nclass, double cxh[], int ch[], int n, int start, int th[]) {
    if (nclass < 1) return 0;
    if (nclass == 1) {
        int ws; double ms;
        if (start == 0) {
            ws = ch[n - 1];
            ms = cxh[n - 1] / ws;
        } else {
            ws = ch[n - 1] - ch[start - 1];             // start부터 끝까지 돗수;
            ms = (cxh[n - 1] - cxh[start - 1]) / ws;    // start부터 끝까지 평균값;
        }
        th[0] = n - 1;
        return ws * ms * ms;                            // weighted mean;
    }

    double gain_max = -1;
    int *tt = new int [nclass - 1];
    for (int j = start; j < n; j++) {
        int wj; double mj;
        if (start == 0) {
            wj = ch[j]; 
            mj = cxh[j];                    //mj = cxh[j] / wj;
        }
        else {
            wj = ch[j] - ch[start - 1];     //start부터 j까지 돗수;
            mj = (cxh[j] - cxh[start - 1]); //mj = (cxh[j] - cxh[start - 1]) / wj;
        }
        if (wj == 0) continue;
        mj /= wj;                           //start부터 j까지 평균;
        double gain = wj * mj * mj + histo_partition(nclass - 1, cxh, ch, n, j + 1, tt);
        if (gain > gain_max) {
            th[0] = j;
            for (int k = nclass - 1; k > 0; k--) th[k] = tt[k - 1];
            gain_max = gain;
        }
    }
    delete [] tt;
    return gain_max;
};

trimodal 분포의 분리;

class0: 0~th[0]

class1: (th[0]+1)~th[1],

class2: (th[1]+1)~th[2]=255

th[0]=103, th[1]=172  (th[2]=255)
th[0]=88, th[1]=176, th[2]=255

 

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

Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Kuwahara Filter  (2) 2020.12.28
Moving Average을 이용한 Thresholding  (0) 2020.11.26
Posted by helloktk

댓글을 달아 주세요

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

영상에서 추출한 객체(object)의 경계선은 객체의 모양에 대한 많은 정보를 담고 있지만, 어떤 경우에는 불필요하게 많은 정보가 될 수도 있다. 이 경우 원래의 경계를 충분히 닮은 다각형으로 간략화하여서 불필요한 정보를 제거할 수 있다. 이 다각형으로의  근사가  물체의 경계를 얼마나 제대로 표현할 수 있는지에 대한 기준이 있어야 한다. N개의 점 $\{(x_i, y_i) | i = 1... N\}$으로 표현된 경계가 있다고 하자. 가장 단순한 근사는 가장 멀리 떨어진 두 점(A, B)으로 만들어진 선분일 것이다(선분의 길이가 대략적인 물체의 크기를 주고, 선분 방향이 물체의 orientation을 알려준다). 이 선분을 기준으로 다시 가장 멀리 떨어진 점(C)을 찾아서 일정 거리 이상 떨어져 있으면 꼭짓점으로 추가하여 두 개의 선분으로 분할한다. 두 개의 선분에 대해서 동일한 분할 작업을 재귀적으로 반복하므로 좀 더 정밀하게 물체를 표현하는 다각형을 찾을 수 있다. 선분에서 가장 멀리 떨어진 지점이 기준 거리 이내이면 분할을 멈춘다. 닫힌 다각형 근사일 때는 시작점 = 끝점으로 조건을 부여하면 된다.

** open polyline으로 우리나라 경계를 단순화시키는 예: top-right 쪽이 시작과 끝이므로 simplication 과정에서 변화지 않는다.

// distance square from (x, y) to the line segment AB ;

double pt2SegDist2(int x, int y, CPoint A, CPoint B) ;

더보기
double pt2SegDist2(int x, int y, CPoint A, CPoint B) {
    double dx = B.x - A.x, dy = B.y - A.y ;
    double lenAB2 = dx * dx + dy * dy ;
    double du = x - A.x, dv = y - A.y ;
    if (lenAB2 == 0.) return du * du + dv * dv;
    double dot = dx * du + dy * dv ; 
    if (dot <= 0.) return du * du + dv * dv;
    else if (dot >= lenAB2) {
        du = x - B.x; dv = y - B.y;
        return return du * du + dv * dv;
    } else {
        double slash = du * dy - dv * dx ;
        return slash * slash / lenAB2;
    }
};

//  재귀 호출을 이용한 Douglas-Peucker 알고리즘;

void DouglasPeucker(double tol, CPoint vtx[], int istart, int iend, int mark[]) {
    if (iend <= istart + 1) return; //종료 조건;
    // vtx[istart] to vtx[iend]을 잇는 선분을 기준으로 분할점을 찾음;
    int ibreak = istart;			// 선분에서 가장 먼 꼭짓점;
    double maxdist2 = 0.0;
    double tolsq = tol * tol;		
    for (int i = istart + 1; i < iend; i++) {
        double dist2 = pt2SegDist2(vtx[i].x, vtx[i].y, vtx[istart], vtx[iend]);
        if (dist2 <=  maxdist2) continue;
        ibreak = i;
        maxdist2 = dist2;
    } 
    if (maxdist2 > tolsq) {		// 가장 먼 꼭짓점까지 거리가 임계값을 넘으면 => 분할;
        mark[ibreak] = 1;		// vtx[ibreak]를 마킹;
        // vtx[ibreak] 기준으로 좌/우 polyline을 분할시도;
        DouglasPeucker(tol, vtx, istart, ibreak, mark);
        DouglasPeucker(tol, vtx, ibreak, iend, mark);
    }
    return;
};

// driver routine;

int polySimplify(double tol, std::vector <CPoint>& V, std::vector <CPoint>& out );

더보기
int polySimplify(double tol, std::vector <CPoint>& V, std::vector <CPoint>& out) {
    if (V.size() < 2) return 0;    
    std::vector<int> mark(V.size(), 0); // 꼭짓점 마킹용 버퍼;
    // DP-알고리즘 호출 전에 너무 가까이 붙어있는 꼭짓점을 제거하면 보다 빠르게 동작하게 만들 수 있다;
    mark[0] = mark[V.size() - 1] = 1;   // 처음과 끝은 marking;
    DouglasPeucker( tol, &V[0], 0, V.size() - 1, &mark[0] );
    // save marked vertices;
    out.clear();    
    for (int i = 0; i < V.size(); i++)
        if (mark[i]) out.push_back(V[i]);
    return out.size();
}
Posted by helloktk

댓글을 달아 주세요

728x90

이미지에서 물체를 인식할 때 물체의 윤곽 정보는 매우 중요한 역할을 한다. 이 윤곽 정보는 보통의 경우 에지를 검출한 후 적절한 영상처리를 거쳐서 단일 픽셀의 선으로 표현된다. 그러나 이렇게 만든 윤곽선은 불필요하게 많은 점으로 구성이 되어 있는 경우가 많으므로, 실제 결과를 사용할 때 적당한 길이의 직선으로 연결된 폴리 라인(polyline)으로 근사를 하여 사용한다. 이 경우 윤곽선의 연접한 점들을 어디까지 직선으로 판단할 것인가에 대한 문제가 생긴다. 

이에 대한 한 가지 접근법은 먼저 윤곽선의 처음과 끝점을 잇는 선분을 만들고 사이의 점들이 이 선분에서 벗어난 정도를 조사한다. 최대로 벗어난 윤곽선의 점이 폴리 라인의 새 꼭짓점이 되고 그 꼭짓점을 기준으로 이전 선분을 두 개로 분할한다. 다시 분할된 두 선분에 대해서 같은 작업을 시행하여 더 이상 작은 선분으로의 분할이 없을 때까지 반복하는 것이다. 이 방법은 잘 알려진 Douglas-Peucker 알고리즘이다. 

두 번째 방법은 윤곽선의 처음 점에서 시작하는 직선의 끝점의 위치를 윤곽선을 따라 순차적으로 옮겨가면 선분을 만들고, 사이의 점들이 이 선분에서 너무 멀리 떨어지지 않는가를 체크한다. 선분에서 너무 벗어난 점이 있으면, 바로 직전의 윤곽선 점을 꼭짓점으로 분류한 후 이 꼭짓점에서 시작하는 선분을 앞에서와 같은 방식으로 만들어서 다음 꼭짓점 후보를 찾는다. DP 알고리즘이 분할에 기반한 것이라면, 이 방법은 점들의 합병(merging)에 기반한 알고리즘이다. 두 알고리즘 모두 재귀적으로 구현할 수 있다. 아래 그림은 DP-알고리즘(Green Line)과 merge-알고리즘 (Blue Line)을 동시에 비교한 것이다.


거리-tolerance를 10씩 증가시키면(점점 단순해진다) 두 알고리즘을 적용하여 본 결과:

// mark는 새로운 꼭지점의 후보를 기록한다(1=꼭짓점, 0=아님)
// mark[0] = mark[end_idx] = 1;
// end_idx = size(v) - 1: not changing;
void 
MergeSimplication(double tol, CfPt v[], int start_idx, const int end_idx, int mark[]) {
    if (end_idx <= start_idx + 1) return ; //at least 3 ;
    double distmax = 0 ;
    for (int end = start_idx + 2; end <= end_idx; end++) {
        double distmax = 0;
        for (int i = start_idx + 1 ; i < end; i++) {
            double dist = SegmentDistance(v[start_idx], v[end], v[i]);
            if (distmax < dist) distmax = dist;
        }
        if (distmax > tol) break; 
    } 
    // 정상적으로 끝나면 end==end_idx+1;
    if (end <= end_idx) {
        //end에서 tolerance을 넘었으므로 새로운 vertex은 end-1;
        mark[end - 1] = 1; //new_vertex;
        MergeSimplication(tol, v, end - 1, end_idx, mark);
    }
    return;
};

// SegmentDistance()는 한 점에서 선분까지 거리를 구하는 함수;
// distantance tolerance 이내에 근접한 이웃점들은 미리 제외시킨 후 알고리즘을 적용한다.

Posted by helloktk

댓글을 달아 주세요