728x90

$m$개의 control point $\{\mathbf{Q}_i\}$가 주어지면  $p$차의 B-Spline curve는 basis 함수를 $N_{i, p}(t)$을 써서

$$ \mathbf{B}(t) = \sum_{i=0}^{m-1} N_{i, p}(t) \mathbf{Q}_i $$

로 표현할 수 있다. 이를 이용해서 일정한 순서로 샘플링된 평면상의 $N$개의 입력점 $\{ \mathbf{P}_i \}$을 찾아보자. B-spline 곡선을 이용하면 이 문제는 control 점을 찾는 문제가 된다. 곡선이 입력점을 잘 표현하기 위해서는 곡선과 입력점과의 차이를 최소로 하는 control 점을 찾아야 한다:

$$    \mathbf{ Q}^* = \text{argmin}(L), \quad\quad L:= \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2 $$

여기서 $\{ t_k| k=0,1,...,N-1\}$는 입력점이 얻어지는 sample time으로 $0= t_0\le t_1\le...\le t_{N-1}= 1$로 rescale 할 수 있다. 

행렬을 이용해서 식을 좀 더 간결하게 표현할 수 있다. 

$$L = \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} N_{i, p}(t_k) \mathbf{Q}_i -  \mathbf{P}_k \right|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} A_{ki} \mathbf{Q}_i - \mathbf{P}_k \right|^2, \\ \quad A_{ki} = N_{i, p}(t_k) $$

로 쓸 수 있으므로, $\hat {Q}= (\mathbf{Q}_0, \mathbf{Q}_1,..., \mathbf{Q}_{m-1})^t$, $\hat {P} = (\mathbf{P}_0, \mathbf{P}_1,..., \mathbf{P}_{N-1})^t$인 벡터, 그리고 $ \hat {A} = (A_{ki})$인 행렬로 보면

$$ L= \left| \hat{A} \cdot \hat {Q} - \hat {P} \right|^2 =  (\hat{Q}^t \cdot \hat {A}^t -\hat{P}^t ) \cdot (\hat{A}\cdot \hat{Q} - \hat{P}).$$

위 식의 값을 최소로 하는 최소자승해는 $\hat {Q}^t$에 대한 미분 값이 0인 벡터를 찾으면 된다; $$ \frac {\partial L}{\partial \hat {Q}^t } = \hat {A}^t \cdot \hat {A} \cdot \hat {Q} - \hat {A}^{t} \cdot \hat {P} = 0.$$ 이 행렬 방정식의 해는

$$ \hat {Q}^* = ( \hat {A} ^t \cdot \hat {A})^{-1} \cdot ( \hat {A}^t \cdot \hat {P})$$ 로 표현된다. $\hat {A} ^t \cdot \hat {A}$가 (banded) real symmetric ($m\times m$) 행렬이므로 Cholesky decomposion을 사용하면 쉽게 해를 구할 수 있다. $\hat{A}$가 banded matrix 형태를 가지는 이유는 basis가 local support에서만 0이 아닌 값을 취하기 때문이다. 

open b-spline(cubic)

int BSplineFit_LS(std::vector<CPoint>& data,
                  int degree,             // cubic(3); 
                  int nc,                 // num of control points;
                  double X[], double Y[]) // estimated control points;
{
    // open b-spline;
    std::vector<double> knot((nc - 1) + degree + 2);
    for (int i = 0; i <= nc + degree; i++) knot[i] = i;
    
    int ndata = data.size();
    std::vector<double> t(ndata);                // parameter;
    double scale = (knot[nc] - knot[degree]) / (ndata - 1);
    for (int i = 0; i < ndata; i++) 
        t[i] = knot[degree] + scale * i;

    std::vector<double> A(ndata * nc);
    for (int i = 0; i < ndata; i++)
        for (int j = 0; j < nc; j++)
            A[i * nc + j] = Basis(j, degree, &knot[0], t[i]); //A(i,j)=N_j(t_i)

    // S = A^t * A; real-symmetric matrix;
    std::vector<double> Sx(nc * nc);
    std::vector<double> Sy(nc * nc);
    for (int i = 0; i < nc; i++) {
        for (int j = 0; j < nc; j++) {
            double s = 0;
            for (int k = 0; k < ndata; k++)
                s += A[k * nc + i] * A[k * nc + j];
            Sx[i * nc + j] = s;
        }
    }
    //copy;
    for (int i = 0; i < nc * nc; i++) Sy[i] = Sx[i];
    // X = A^t * P.x;  Y = A^t * P.y
    for (int i = 0; i < nc; i++) {
        double sx = 0, sy = 0;
        for (int k = 0; k < ndata; k++) {
            sx += A[k * nc + i] * data[k].x;
            sy += A[k * nc + i] * data[k].y;
        };
        X[i] = sx; Y[i] = sy;
    };
    // solve real symmetric linear system; S * x = X, S * y = Y;
    // solvps(S, X) destories the inputs;
    // ccmath-2.2.1 version;
    int res1 = solvps(&Sx[0], X, nc);
    int res2 = solvps(&Sy[0], Y, nc);
    return res1 == 0 && res2 == 0;
};

**네이버 블로그 이전;

'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

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

Fixed-point version: 선분 길이가 4096 픽셀 정도까지는 1-pixel 이내의 오차로 그려진다.

void DDALine(CPoint A, CPoint B) {
    const int mulfac = 4096; // 2^12 = 4096;
    int dx = B.x - A.x;
    int dy = B.y - A.y;
    int adx = dx < 0 ? -dx : dx;
    int ady = dy < 0 ? -dy : dy;
    int steps = adx < ady ? ady : adx;
    if (steps == 0) SetPixel(A.x, A.y); //single point;
    else {
        int curX = A.x * mulfac;      
        int curY = A.y * mulfac;
        dx = (dx * mulfac) / steps;
        dy = (dy * mulfac) / steps;
        while (steps-- >= 0) {
            SetPixel(curX / mulfac, curY / mulfac);
            curX += dx;
            curY += dy;
        }
    }
}

'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

댓글을 달아 주세요

B-Spline

Computational Geometry 2021. 4. 25. 09:43
728x90

$$ \mathbf{C}(u) = \sum_{i=0}^{n} N_{i,p} (u) \mathbf{P}_i$$

where $N_{i, p}(u)$'s are B-spline basis functions of degree $p$, $n+1$ control points $\{\mathbf{P}_i\}$, and a knot vector $\mathbf{U}=(u_0, u_1, ..., u_m)$. Note, $m=n + p + 1$.

B-spline의 basis 함수 $N_{i, p}(t)$는 유한한 구간에서만 0이 아닌 값을 갖는(local support) 함수로 다음과 같은  recursion relation을 가진다:

\begin{align} N_{i,0}(t) &= \left\{ \begin{matrix} 1 & \mathrm{if} \quad u_i \leq x \lt u_{i+1} \\ 0 & \mathrm{otherwise} \end{matrix} \right. \quad \quad (0 \le i \le n+p), \\ N_{i,p}(t) &= \frac{t - u_i}{u_{i+p} - u_i} N_{i,p-1}(t) + \frac{u_{i+p+1} - t}{u_{i+p+1} - u_{i+1}} N_{i+1,p-1}(t).\end{align}

$N_{i,p}(t)$의 local support는 $[u_i, u_{i+p+1}]$로 주어진다.

일반적으로 knot vector가 특별한 구조를 가지지 않으면 B-spline곡선은 Bezier 곡선의 경우와 달리 처음과 끝 control 점을 통과하지 않는다(open B-spline). 처음과 끝점을 통과하도록 만들기 위해서는 knot 벡터의 시작 $(p+1)$개 성분이 같은 값을 갖고, 마지막 $(p+1)$개 성분도 같은 값을 갖도록 조정하면 된다(clamped B-spline):

$$clamped: u[0]=u[1]=...=u[p],...,u[n+1]=...=u[n+p]=u[n+p+1]$$

knot는 균일할 필요는 없지만, 균일한 경우는 보통

$$ u_i = \left\{ \begin{array} {rcl} 0 & , & 0 \le i \le p, \\ \frac{i-p}{n+1-p}&, & p+1\le i \le n \\ 1&,& n+1\le i \le n+p+1\end{array}\right.$$

로 선택한다. 

ref: ae.sharif.edu/~aerocad/curve%20modelling-B_Spline.pdf

clamped(black), open(blue)

// Note, p = degree;
void BSpline(int p, std::vector<CPoint>& control, int resolution, CDC *pDC) {
    int n = control.size() - 1;  ASSERT(n > 0);
    std::vector<double> u(n + p + 2);                // knot vector;
    calculateKnots(&u[0], n, p);                     // uniform knot;
    double delta = (u[n + p + 1] - u[0]) / (resolution - 1);  // parameter increment;
    CPoint Q;
    for (double t = u[0] ; t <= u[n+p+1]; t += delta) {
        calculatePoint(&u[0], n, p, t, &control[0], &Q);
        if (i == 0) pDC->MoveTo(Q);
        pDC->LineTo(Q);
    }
    pDC->LineTo(control[n]);
}
// de Boor Cox's algorithm;
double Blend(int i, int p, double *u, double t) {  
    if (p == 0) { // termination condition of recursion;
        if ((u[i] <= t) && (t < u[i+1]))        return 1;
        else                                    return 0;
    }
    else {
        double coef1, coef2;
        if (u[i+p] == u[i]) {
            if (t == u[i])  coef1 = 1;
            else            coef1 = 0;
        } 
        else                coef1 = (t - u[i]) / (u[i+p] - u[i]); 
        
        if (u[i+p+1] == u[i+1]) {
           if (t == u[i+1]) coef2 = 1;
           else             coef2 = 0;
        } 
        else                coef2 = (u[i+p+1] - t) / (u[i+p+1] - u[i+1]);
        return coef1 * Blend(i, p-1, u, t) + coef2 * Blend(i+1, p-1, u, t);
    }
}
// clamped b-spline knot vector; uniform example;
// u[0]=u[1]= ... =u[p],.....,u=[n+p-2]=u[n+p-1]=u[n+p]=u[n+p+1]
void calculateKnots(double *u, int n, int p) {
    int j = 0; 
    while (j <= p)     u[j++] = 0;                 // multiplicity = p+1;
    while (j <= n)     u[j++] = j - p;             // m = n+p+1;
    while (j <= n+p+1) u[j++] = n - p + 1;         // multiplicity = p+1;
}
void calculatePoint(double *u, int n, int p, double t, CPoint *control, CPoint *output) {
    double x = 0, y = 0;
    for (int i = 0; i <= n; i++) {
        double b = Blend(i, p, u, t);
        x += control[i].x * b;
        y += control[i].y * b;
    }
    output->x = int(x + 0.5);
    output->y = int(y + 0.5);
}

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

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
Convexity in Bezier Curve  (0) 2021.04.22
Posted by helloktk

댓글을 달아 주세요

728x90

// test flatness of a bezier curve with control points {p1, c1, c2, p2};
bool FlatBezier(const double eps, 
                const CfPt& p1, const CfPt& c1, const CfPt& c2, const CfPt& p2) {
    CfPt U = c1 - (p1 * 2 + p2) / 3;
    CfPt V = c2 - (p1 + p2 * 2) / 3;
    double normU = hypot(U.x, U.y);
    double normV = hypot(V.x, V.y);
    return (3. * max(normU, normV) / 4.) <= eps;
}
// subdivision stage; De Casteljau's algorithm
void GetBezierPoints(const double eps, 
                     const CfPt &p1, const CfPt &c1, const CfPt &c2, const CfPt &p2,
                     std::vector<CfPt>& smoothed) {
    if (!FlatBezier(eps, p1, c1, c2, p2)) {
        // Subdivide the curve at t = 0.5
        // left;
        CfPt mid_segm = (p1  + c1 * 3 + c2 * 3 + p2) / 8.0;
        CfPt new_c1   = (p1 + c1) / 2.0;
        CfPt new_c2   = (p1 + c1 * 2 + c2) / 4.0;
        GetBezierPoints(eps, p1, new_c1, new_c2, mid_segm, smoothed);
        // right;
        new_c1  = (c1 + c2 * 2 + p2) / 4.0;
        new_c2  = (c2 + p2) / 2.0;
        GetBezierPoints(eps, mid_segm, new_c1, new_c2, p2, smoothed);
    } 
    else smoothed.push_back(p2); // flat enough;
        
}
// linear interpolation between a and b;
CfPt lerp(const double t, const CfPt &a, const CfPt &b) {
    return a * (1 - t) + b * t;
}
void SmoothPathWithBezier(std::vector<CfPt>& points, std::vector<CfPt>& smoothed) {
    int N = points.size();
    if (N < 3) return ; 
    smoothed.clear();
    CfPt cntl[4];
    CfPt *pts = &points[0];   // iterator;
    smoothed.push_back(pts[0]);    
    for (int i = 2; i < N; i++, pts++) {
        if (i == 2) {
            cntl[0] = pts[0];
            cntl[1] = lerp(2./ 3, pts[0], pts[1]);
        } else {
            //0번 control-pts= (0-1)의 중간점;
            //1번 control-pts= (0-1)의 5:1 내분점; 
            cntl[0] = lerp(1./ 2, pts[0], pts[1]);
            cntl[1] = lerp(5./ 6, pts[0], pts[1]);
        }
        if (i == N - 1) {
            //1:2 내분점;
            cntl[2] = lerp(1./ 3, pts[1], pts[2]);
            cntl[3] = pts[2];
        } else {
            //2번 control-pts;(1-2)의 1:5 내분점;
            //3번 control-pts;(1-2)의 중간점;
            cntl[2] = lerp(1./ 6, pts[1], pts[2]);
            cntl[3] = lerp(1./ 2, pts[1], pts[2]);
        }
        //입력 중에 처음 두개가 같거나, 나중 두개가 같으면 마지막이 control pt;
        if ((pts[0] == pts[1]) || (pts[1] == pts[2])) 
            smoothed.push_back(cntl[3]);
        else {
            double tolerance = 1.; // 1-pixel;
            GetBezierPoints(tolerance, cntl[0], cntl[1], cntl[2], cntl[3], smoothed);
        }
    }
}

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

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
Convexity in Bezier Curve  (0) 2021.04.22
De Casteljau's Algorithm  (0) 2021.04.22
Posted by helloktk

댓글을 달아 주세요

728x90

cubic Bezier 곡선은 네 개의 컨트롤 점($\mathbf {P}_0, 주어진다:

$$\mathbf {B}(t) = (1-t)^3 \mathbf {P}_0 + 3(1-t)^2 t \mathbf {P}_1 + 3 (1-t) t^2 + t^3 \mathbf {P}_3.$$

만약 곡선이 충분히 평탄하다면 굳이 삼차 다항식을 계수로 가지는 Bezier 곡선으로 표현할 필요가 없이$\mathbf {P}_0$와 잇는 직선으로 근사하는 것이 여러모로 유리할 것이다. 그러면 Bezier 곡선의 평탄성을 어떻게 예측할 것인가? 이는 cubic Bezier 곡선과 시작점과 끝점을 연결하는 직선 $\mathbf {L}(t)= (1-t) \mathbf {P}_0 + t \mathbf {P}_3$과의 차이를 비교하므로 판단할 수 있다.

 

직선도 cubic Bezier 곡선 표현으로 할 수 있다. 시작점이 $\mathbf {P}_0$이고 끝점이 $\mathbf {P}_3$일 때 이 둘 사이를 3 등분하는 두 내분점 $(2\mathbf {P}_0 + \mathbf{P}_3)/3$, $( \mathbf {P}_0 + 2\mathbf {P}_3)/3$을 control point에 포함시키면,

$$\mathbf {L}(t)= (1-t)^3  \mathbf {P}_0 + (1-t)^2 t (2 \mathbf {P}_0 +\mathbf {P}_3)+ (1-t) t^2  (\mathbf {P}_0 + 2 \mathbf {P}_3)+ t^3 \mathbf {P}_3$$

처럼 쓸 수 있다. Bezier 곡선과 직선의 차이를 구하면

\begin{align} \Delta (t) &= |\mathbf {B}(t) - \mathbf {L}(t)| \\   &= |(1-t)^2 t (3\mathbf {P}_1 -2 \mathbf {P}_0 -\mathbf{P}_3) + (1-t) t^2(3 \mathbf {P}_2 - \mathbf {P}_0 - 2 \mathbf {P}_3)| \\ &= 3 (1-t) t |(1-t) \mathbf {U} + t \mathbf {V}|, \\ \\ &\quad \mathbf {U}=\mathbf {P}_1 - \frac{1}{3}(2 \mathbf {P}_0 + \mathbf {P}_3),   \quad \mathbf {V}=\mathbf {P}_2 - \frac{1}{3} (\mathbf {P}_0 + 2 \mathbf {P}_3)  \end{align}

그런데 

$$ \text {max} \{3(1 - t) t \} = \frac {3}{4}, \quad 0\le t \le 1, \\ \text {max}\{ |(1-t) \mathbf {U} + t \mathbf {V}| \} =\text {max}\{ |\mathbf{U}| , |\mathbf {V}|\} $$

이므로 아래와 같은 벗어난 거리의 상한 값을 얻을 수 있다:

$$ \Delta (t) \le \frac {3}{4} \text{max} \{ |\mathbf {U} | , |\mathbf {V}| \}$$

이 상한값이 주어진 임계값보다 작으면 3차 Bezier 곡선은 두 끝점을 잇는 직선으로 근사할 수 있다. 좀 더 기하학적으로 표현하면

$$|\mathbf{U}| = \mathbf{P}_1\text{에서 왼쪽 3등분 내분점까지 거리}$$

$$|\mathbf {V}|= \mathbf {P}_2 \text{에서 오른쪽 3등분 내분점까지 거리}$$

이므로,

$$ \Delta(t) \le \frac {3}{4}\times( \text{control point와 3등분 내분점 간의 거리의 최댓값}).$$

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

B-Spline  (0) 2021.04.25
Bezier Curve Smoothing  (0) 2021.04.23
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
Posted by helloktk

댓글을 달아 주세요

728x90

Bezier 곡선은 control point $\{ \mathbf {P}_i\}$의 선형 결합으로 주어진다:

$$\mathbf {B}(t) = \sum_{i=0}^{n} B_{i, n} (t) \mathbf {P}_i , \quad B_{i, n}(t)=\left(\begin {array}{c} n \\ i \end {array} \right) t^i (1-t)^{n-1}.$$

Bernstein 다항식이 선형 결합의 가중치를 역할을 하는데 0과 1 사이의 양의 실수 값을 가진다. 그리고 이들 가중치의 합은 1이다:

$$ 0\le  B_{i, n}(t) \le 1, \quad i=0,1,2,... n ,    \quad    0\le t\le 1 \\\sum_{i=0}^{n}  B_{i, n}(t) = 1$$

따라서 Bezier 곡선은 control points로 만든 convex region에 있음을 알 수 있다. Bezier 곡선의 convexity 특성은 여러 가지 좋은 특성을 준다.  몇 가지만 나열하면, 첫째가 Bezier 곡선은 항상 컨트롤 포인트의 convex hull 내에 놓이게 되므로 곡선의 제어 가능성을 보장한다. 둘째, 교차 여부를 쉽게 확인할 수 있다. 셋째는 culling을 가능하게 한다.

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

Bezier Curve Smoothing  (0) 2021.04.23
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
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

곡선의 길이를 구하고자 하면 곡선을 따라 속력을 적분하면 된다. 곡선을 기술하는 벡터가 $\mathbf {B}(t) = [B_x(t), B_y(t)]$ 로 주어지면 곡선의 길이는

$$\text {Arc Length}(t_1, t_2) =\int_{t_1}^{t_2} | \mathbf {B}'(t)|dt=\int_{t_1}^{t_2} \sqrt{(B'_x(t))^2 + (B'_y(t))^2} dt.$$

그러나 일반적으로 이 적분은 닫힌 형태로 주어지지 않으므로 Gauss-Legendre quadrature와 같은 수치 적분에 의존해서 그 값을 계산하여야 한다.
Bezier curve의 경우는 기하학적인 방법을 이용하여서 좀 더 간편하게 그 길이를 구할 수 있다. Bezier curve의 길이에 대한 간단한 근사는 현의 길이($|\mathbf{Q}_3- \mathbf {Q}_0|$)로 근사하는 것이다. 그러나 현의 길이는 lower bound만을 준다. 다음으로 고려할 근사는  control 점 사이 거리의 합이다. 이것은 실제 길이의 upper bound를 준다. 따라서 근사적인 길이는 이 두 거리의 평균으로 잡으면 될 것이다:

$$ \text{Arc Length}\sim \frac{ \text {현의 길이} + \text {control point 사이 거리 합}}{2}.$$

그러나 이 근사가 유효하기 위해서는 현의 길이와 control point 간의 거리의 합

$$|\mathbf{Q}_3- \mathbf {Q}_2|+ |\mathbf {Q}_2- \mathbf {Q}_1| +|\mathbf {Q}_1- \mathbf {Q}_0|$$

의 차이가 충분히 작아야 한다. 만약에 이 두 길이의 차이기 정해진 임계값보다도 크면 Bezier curve를 둘로 분할해서 좌우 두 Bezier curve에 대해서 동일한 검사를 한다. 충분히 평탄한 경우는 위의 평균값으로 근사하고 그렇지 못한 경우는 충분히 평탄해질 때까지 분할하는 과정을 반복한다. 아래 코드는 3차 Bezier curve의 길이를 이 알고리즘을 이용해서 구하는 과정을 구현한 것이다.

#define DistL2(A, B) sqrt(((A).x-(B).x)*((A).x-(B).x) + ((A).y-(B).y)*((A).y-(B).y))
static void BezierSubDiv(CfPt Q[4], CfPt leftQ[4], CfPt rightQ[4]) {
    CfPt T[4][4];                      /* Triangle Matrix */
    for (int j = 0; j < 4; ++j) T[0][j] = Q[j];
    // de Casteljau divides and conquers at t = 0.5;
    for (int i = 1; i < 4; ++i)
        for (int j = 0 ; j <= 3 - i; ++j)
            T[i][j] = (T[i - 1][j] + T[i - 1][j + 1]) / 2;
            
    // left subdiv control points;
    for (int j = 0; j < 4; ++j) leftQ[j]  = T[j][0];
    // right subdiv control points;
    for (int j = 0; j < 4; ++j) rightQ[j] = T[3 - j][j];
}                                        
void BezierLength(CfPt Q[4], double tol, double *length) {
    CfPt leftQ[4], rightQ[4];                
    double controlLen = 0;  
    for (int i = 0; i < 3; ++i) controlLen += DistL2(Q[i], Q[i + 1]);
    double cordLen = DistL2(Q[0], Q[3]);
    // test flatness;
    if (fabs(controlLen - cordLen) > tol) {
        BezierSubDiv(Q, leftQ, rightQ);              /* divide */
        BezierLength(leftQ, tol, length);            /* left side */
        BezierLength(rightQ, tol, length);           /* right side */
        return;
    }
    *length = *length + (cordLen + controlLen) / 2 ;
    return ;
}

De Casteljau's algorithm

int main() {
    double k[4] = {0.5522847498, //touch at t = 0.5;
                   0.551915023,  //min-deviation;
                   0.551778477,  //match area;
                   0.551777131   //match length;
                   };
    for (int i = 0; i < 4; i++)
        printf("bezier length = %.8f\n", BezierLengthCircle(k[i]));
    printf("exact length  = %.8f\n", 2.0 * atan(1.));
    return 0;
}
double BezierLengthCircle(double k) {
    CfPt Q[4] = {{0, 1.}, {k, 1.}, {1., k}, {1., 0.}};
    double len = 0;
    double tol = 1.e-20;//
    BezierLength(Q, tol, &len);
    return len;
};

 

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

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
Bresenham's Line Algorithm  (0) 2021.04.07
Posted by helloktk

댓글을 달아 주세요

728x90

Bezier 곡선을 이용한 원의 근사처럼 타원을 근사해보도록 하자. 원점을 중심으로 하고 장축 반지름이 $a$, 단축 반지름이 $b$인 타원을 1사분에서 3차 Bezier curve을 이용해서 근사하려면 4개의 control point가 필요하다. 원의 경우처럼, 끝점에서 접선의 기울기를 같게 하는 조건을 부여하면 control point는 

$$\mathbf {P_0} = (0, b), \quad \mathbf {P}_1 = (ka, b), \quad \mathbf {P}_2 =( a, kb),\quad \mathbf {P}_3 = (a, 0)$$

로 잡을 수 있다. 따라서

$$\mathbf {B}(t) = (1-t^3) \mathbf {P}_0 + 3t(1-t)^2 \mathbf {P}_1 + 3t^2 (1-t) \mathbf {P}_2 + t^3 \mathbf {P}_3 = \left(\begin {array}{c} a x(t) \\ b y(t) \end {array}\right) $$

$$ x(t) = 3k (1-t)^2 t + 3 (1-t) t^2 + t^3 \\ y(t) = 3k t^2 (1-t) + 3t(1-t)^2 + (1-t)^3 $$

$t=1/2$일 때 $(a/\sqrt {2}, b/\sqrt {2})$을 통과하는 조건을 부여하면, 원과 마찬가지로

$$ k = \frac {4}{3}(\sqrt {2}-1)= 0.5522847498...$$

을 얻는다. Mahalanobis measure를 기준으로 거리를 측정하면  타원의 경우도 벗어남 에러가

$$ \Delta (t) = \sqrt { \frac {B_x^2(t)}{a^2} + \frac {B_y^2(t)}{b^2} }  -1 =\sqrt {x^2(t)+y^2(t)}-1$$

원의 경우와 같음을 쉽게 알 수 있다.

회전된 타원; 

2사분면: $(x, y) \rightarrow (-x, y)$

3사분면: $(x, y) \rightarrow (-x, -y)$

4사분면: $(x, y) \rightarrow (x, -y)$

void BezierEllipse(CDC *pDC, CPoint center, double a, double b) {
    const double k = 0.5522847498;
    CPen red(PS_SOLID, 3, RGB(0xFF, 0, 0));
    CPen *pOld = pDC->SelectObject(&red);
    CPoint P[4];  //control pts;
    P[0] = CPoint(center.x,                    int(center.y - b + 0.5));
    P[1] = CPoint(int(center.x + k * a + 0.5), int(center.y - b + 0.5));
    P[2] = CPoint(int(center.x + a + 0.5),     int(center.y - k * b + 0.5));
    P[3] = CPoint(int(center.x + a + 0.5),     center.y);
    pDC->PolyBezier(P, 4);
    pDC->SelectObject(pOld);
}
void BezierEllipse(CDC *pDC, CPoint center, double a, double b, double ang) {
    const double k = 0.5522847498;
    const double cosang = cos(ang);
    const double sinang = sin(ang);
    CPoint P[4];
    double x[4], y[4], xt[4], yt[4];
    //1사분면: 
    x[0] = 0,     y[0] = b;
    x[1] = k * a, y[1] = b;
    x[2] = a,     y[2] = k * b;
    x[3] = a,     y[3] = 0;
    CPen red(PS_SOLID, 3, RGB(0xFF, 0, 0));
    CPen *pOld = pDC->SelectObject(&red);    
    for (int i = 0; i < 4; i++) {
        xt[i] = x[i] * cosang - y[i] * sinang;
        yt[i] = x[i] * sinang + y[i] * cosang;
        P[i] = CPoint(int(center.x + xt[i] + 0.5), int(center.y - yt[i] + 0.5));
    }
    pDC->PolyBezier(P, 4);
    pDC->SelectObject(pOld);
    //4사분면;(x, -y);
    CPen blue(PS_SOLID, 3, RGB(0, 0, 0xFF));
    pOld = pDC->SelectObject(&blue);
    for (int i = 0; i < 4; i++) {
        xt[i] = x[i] * cosang - (-y[i]) * sinang;
        yt[i] = x[i] * sinang + (-y[i]) * cosang;
        P[i] = CPoint(int(center.x + xt[i] + 0.5), int(center.y - yt[i] + 0.5));
    }
    pDC->PolyBezier(P, 4);
    pDC->SelectObject(pOld);
    //3-사분면;(-x,-y)
    CPen green(PS_SOLID, 3, RGB(0, 0xFF, 0));
    pOld = pDC->SelectObject(&green);
    for (int i = 0; i < 4; i++) {
        xt[i] = (-x[i]) * cosang - (-y[i]) * sinang;
        yt[i] = (-x[i]) * sinang + (-y[i]) * cosang;
        P[i] = CPoint(int(center.x + xt[i] + 0.5), int(center.y - yt[i] + 0.5));
    }
    pDC->PolyBezier(P, 4);
    pDC->SelectObject(pOld);
    //2사분면;(-x,y);    
    CPen magenta(PS_SOLID, 3, RGB(0xFF, 0, 0xFF));
    pOld = pDC->SelectObject(&magenta);
    for (int i = 0; i < 4; i++) {
        xt[i] = (-x[i]) * cosang - y[i] * sinang;
        yt[i] = (-x[i]) * sinang + y[i] * cosang;
        P[i] = CPoint(int(center.x + xt[i] + 0.5), int(center.y - yt[i] + 0.5));
    }
    pDC->PolyBezier(P, 4);
    pDC->SelectObject(pOld);
}

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

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
Bresenham's Line Algorithm  (0) 2021.04.07
Rotating Calipers  (3) 2021.03.31
Posted by helloktk

댓글을 달아 주세요

728x90

한 개의 Bezier 곡선을 이용해서 원을 표현할 수 없음은 잘 알려진 사실이다. 그럼 Bezier 곡선을 이용해서 얼마나 원(호)을 잘 근사할 수 있을까? 3차 Bezier 곡선을 이용해서 원점에 중심을 둔 반지름 1인 원의 1 사분면 원호를 근사해보도록 하자. 3차 Bezier 곡선은 4개의 control point $\{ \mathbf {P}_i | i=0,1,2, 3\}$이 주어진 경우

$$ \mathbf {B}(t) = (1-t^3) \mathbf {P}_0 + 3t(1-t)^2 \mathbf {P}_1 + 3t^2 (1-t) \mathbf {P}_2 + t^3 \mathbf {P}_3$$

으로 표현된다. 원호 근사에 필요한 control point는 다음 조건을 부여하면 얻을 수 있다. (1) $t=0,1$일 때 원 위에 있어야 하므로

$$\mathbf {B}(t=0) = (0,1) \quad \rightarrow \quad \mathbf {P_0} = (0,1)$$

$$\mathbf {B}(t=1) = (1,0) \quad \rightarrow \quad \mathbf {P_3} = (1,0)$$

그리고 또 그 두 지점에서 접선의 기울기가 원의 기울기와 같아야 하므로

$$ {\mathbf {B}'}(t=0) \propto (1,0)\quad \text {and}\quad {\mathbf {B}'}(t=1)\propto (0,-1)$$

에서

$$ \mathbf {P}_1 = (k, 1), \quad \mathbf {P}_2 = (1, k)$$

처럼 잡을 수 있다. (2) Bezier 곡선의 중간지점이  원 위에 있도록 조건을 부여하면

$$ \mathbf {B}(t=1/2) = (1/\sqrt {2}, 1/\sqrt {2})$$

을 얻고, 이를 이용하면

$$k= \frac {4}{3} (\sqrt {2}-1) = 0.5522847498...$$

을 얻는다.

그럼 원에서 얼마나 벗어날까? 원 중심에서 거리를 차이를 구해보면 Bezier 곡선이 항상 원의 바깥으로 치우쳐 있음을 알 수 있다:

$$\Delta(t) = ||\mathbf {B}(t)||-1\ge 0$$ 

최대로 벗어나는 정도는 $t=(3\pm \sqrt{3})/6$일 때 $\Delta_\text {max}=\frac{1}{3}\sqrt{ \frac{71}{6}-2\sqrt{2}}-1=0.00027253...$이므로 대부분의 경우 크게 벗어남이 없는 원의 근사를 준다.

(2-2) $t=1/2$에서 Bezier 곡선이 원을 통과하는 조건 대신 원에서 벗어남을 최소로 하는 조건을 부여하면 더 좋은 근사를 얻을 수 있다:

$$k=\text{argmin}|\Delta|_\text{max}$$

이 경우  Bezier 곡선은 원의 바깥에 놓이지 않고 교차하게 된다. $t=1/2$에서 최솟값, $t=\frac{1}{2}\left(1\pm \frac{\sqrt{3k^2+20k -12}}{2-3k}\right) $일 때 최댓값을 가지는데, 두 값의 절대값이 같게 되려면(closed form이 없다)

$$ k = 0.551915023...$$

을 선택해야 하고, 이 때 벗어남의 최대값은 $|\Delta|_\text{max} =  0.00019607...$이므로 더 좋은 근사가 된다.

 

(2-3) 또 다른 제한조건은 없을까? Bezier 곡선이 만드는 면적이 사분원의 면적을 표현하도록 제한을 가하는 경우:

$$\frac{\pi}{4} = \int_{0}^{1} \frac{1}{2}\left( B_y B'_x - B_x B'_y\right) dt \\ = \int_0^1  \left(-\frac{3}{2} \left(3 k^2 (t-1)^2 t^2+k \left(-2 t^4+4 t^3-6 t^2+4 t-1\right)+2 (t-1) t\right) \right)dt \\= \frac{1}{2}+\frac{3k}{5}-\frac{3k^2}{20}$$ 에서 

$$ k = 2 - \sqrt{ \frac{22 - 5 \pi }{3}} =0.551778477...$$

이다. 이 경우 면적은 같으나 벗어남 오차는 $t=1/2$일 때

$$|\Delta |_\text{max} = \frac{1}{8} \sqrt{332-40 \sqrt{66-15 \pi }-30 \pi }-1= 0.00026849...$$

로 주어지는데, 중심을 지나는 경우보다는 벗어남이 작지만 최소는 아니다.

(2-4) Bezier 곡선의 길이가 원주가 되는 제한조건을 걸 수도 있다:

$$\frac{\pi}{2}  = \int_0^1\sqrt{ (B'_x)^2 + (B'_y)^2}dt.$$

그런데 우측 적분이 closed form으로 주어지지 않는다. 때문에 $k$ 값을 구하기 위해서 전적으로 numerical method에 의존해야되는데,  그 결과만 쓰면

$$k=0.551777131...$$

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

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
Bresenham's Line Algorithm  (0) 2021.04.07
Rotating Calipers  (3) 2021.03.31
Convex Hull Peeling  (0) 2021.03.29
Posted by helloktk

댓글을 달아 주세요

728x90

화면이나 그림에서 직선을 그리려면 직선의 방정식을 만족하는 모든 픽셀들을 찾으면 된다. 하지만 픽셀의 위치는 연속적인 값이 아니라  이산적이므로 그림에 그려진 직선을 구성하는 픽셀은 그 직선을 표현을 하는 식이 주는 값에 가장 가까운 정수를 찾아서 그린 것이다. 직선을 그리기 위해서 실수 연산을 한 후 정수로 반올림하는 과정이 들어가게 된다. 

Bresenham 알고리즘은 평면 위의 두 점 $(x_0, y_0), (x_1, y_1)$을 연결하는 직선을 정수 연산만을 사용하여 그릴 수 있음을 보여준다. 직선의 방정식은 $y = mx +b$의 형태로 표현이 되고, 두 점이 주어지는 경우 기울기는 $m = \frac {y_1 - y_0}{x_1 - x_0} = \frac {\Delta y}{\Delta x}$, $y-$절편은 $b =y_0 - m  x_0$로 얻을 수 있다. 직선의 기울기가 1보다 크면 $x$와 $x+1$에서 변할 때 종속변수 $y$값의 차이가 커지므로 선이 연속적으로 그려지지 않게 보인다.  이 경우 $y$를 독립변수로 사용하면 되고, 기울기가 음수면 시작과 끝을 바꾸면 된다. 따라서 기울기가 $0\le m \le 1$인 경우만 살펴보면 충분하다. 

이제, $(x_k, y_k)$ 위치에서 직선의 픽셀이 정해진 경우,  $x_{k+1}= x_k + 1$에서 직선식의 값 $y =m (x_k + 1) + b$는 어떤 정수 값으로 대체되어야 하는가? 기울기가 1보다 작으므로 $y$는 구간 $[y_k, y_k +1]$에서 값을 취한다. 즉,  오른쪽 옆이나 오른쪽 위 픽셀에 점을 찍어야 한다. 수식으로는 구간의 아래($d_{lower}$)와 위쪽($d_{upper}$)과의 차이를 비교하여 작은 쪽으로 결정하면 된다:

$$ d_{lower} = y - y_k = m(x_k + 1) +b - y_k$$

$$ d_{upper} = (y_k + 1) - y = y_k + 1 - m (x_k + 1) -b$$

이 과정을 실수연산 없이 판별하는 방법을 알아내자. $d_{lower}$와 $d_{upper}$의 차이는 

$$ d_{lower} - d_{upper} = 2m (x_k + 1) - 2 y_k + 2b -1$$

로 표현되고, 나눗셈을 하지 않기 위해서 $\Delta x$을 곱한 값을 discriminator로 사용하자:

$$p_k = \Delta x(d_{lower} - d_{upper} )= 2 \Delta y (x_k + 1) - 2 \Delta x y_k  + \Delta x (2b-1)$$

의 꼴로 표현되므로 직선의 정보($\Delta x, \Delta y$)와 현재 점을 찍은 픽셀($x_k, y_k$)을 알면 다음번에 점을 찍을 픽셀을 정수 연산만 가지고 판별할 수 있는 장점이 있다. $k+1$일 때 판별식이

$$ p_{k+1} = 2 \Delta y (x_k + 1)  - 2 \Delta x y_{k+1}  + \Delta x (2b-1)$$

이므로 다음과 같은 recursion relation을 만족한다:

$$p_{k+1} = p_k + 2\Delta y - 2\Delta x  (y_{k+1} - y_k) \\ p_0 = 2\Delta y - \Delta x$$

이를 이용하면 시작 픽셀 $(x_0, y_0)$에서 $p_0 = 2\Delta y - \Delta x$을 계산하여 0보다 작지 않으면 $(x_k + 1, y_k+1)$ 픽셀에 점을 찍고,  그다음 점 찍을 픽셀은  $p_{k+1} = p_k + 2\Delta y - 2\Delta x$을 이용해서 판별한다. $p_0 < 0$ 이면 ($x_k +1, y_k$) 픽셀에 점을 찍고,  그다음 점을 찍을 픽셀은 $p_{k+1}= p_k + 2\Delta y$을 계산하여 판별한다.

void BresenhamLine(int x1, int y1, int x2, int y2, DWORD color) {
    int dx = x2 - x1 ;
    int incx = 1;
    if (dx < 0) {
        dx = -dx;
        incx = -1;
    }
    int dy = y2 - y1;
    int incy = 1;
    if (dy < 0) {
        dy = -dy;
        incy = -1;
    }
    setPixel(x1, y1,color);
    if (dx > dy) {
        int p = (dy << 1) - dx;
        while (x1 != x2) {
            if (p >= 0) {
                y1 += incy;
                p -= (dx << 1);
            }
            x1 += incx;
            p += (dy << 1);
            setPixel(x1, y1, color);
        }
    } else {
        int p = (dx << 1) - dy;
        while (y1 != y2) {
            if (p >= 0) {
                x1 += incx;
                p -= (dy << 1);
            }
            y1 += incy;
            p += (dx << 1);
            setPixel(x1, y1, color);
        }
    }
};

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

Bezier Curve Approximation of an Ellipse  (0) 2021.04.11
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
Posted by helloktk

댓글을 달아 주세요

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

평면 위에서 주어진 점집합을 포함하는 가장 작은 원을 찾는 문제는 가장 단순하게는 임의의 두 점이 지름이 되는 원과 임의의 세 점이 만드는 외접원을 모두 조사하여 찾으면 된다. 따라서 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

댓글을 달아 주세요

728x90

기하 알고리즘을 특히 폴리곤 관련 알고리즘을 테스트하기 위해서는 폴리곤 데이터를 만들어야 한다. 여기서 2차원의 점들을 이용하여 간단히 simple polygon을 만드는 방법을 소개한다.

 

단순 폴리곤을 만들기 위해서는 점들을 정렬을 해야 한다. 각각의 점들을 x 축에 프로젝션을 시키면 x 축에 대해서 정렬이 된다(같은 x 값인 경우에는 y의 크기대로 정렬). 따라서 x 값으로 정렬된 점들을 쭉 이어나가면 하나의 poly-line을 만들 수 있다. poly-line의 끝을 잇는다고 해서 단순 폴리곤이 항상 만들어지지는 않는다. 이를 해결하기 위해서는 점들을  한 직선을 기준으로 위-아래로 분리하고, 직선의 윗부분은 x 값이 큰 순서대로 정렬을 하고, 아래 부분은 x 값이 작은 순서대로 정열을 하면 폴리 라인이 아닌 단순 폴리곤을 얻을 수 있다. 직선은 어떻게 잡을까? 이것은 가장 작은 x 값을 갖는 점과 가장 큰 x 값을 같은 점을 잇는 직선을 생각하면 편리하다.

// cross(BC, BA)
// A->B-C: 반시계(C가 AB 왼편: > 0), 시계(C가 AB오른편: < 0), 일직선(0);
int CCW(CPoint A, CPoint B, CPoint C) {
    C.x -= B.x; C.y -= B.y; 
    A.x -= B.x; A.y -= B.y;
    return C.x * A.y - C.y * A.x;
}
int comp(const void *A, const void *B) { //x->y의 크기 순서대로 정렬;
    int ax = ((CPoint *)A)->x - ((POINT *)B)->x ;
    if (ax > 0) return  1;
    if (ax < 0) return -1;
    int ay = ((CPoint *)A)->y - ((CPoint *)B)->y ;
    if (ay > 0) return  1;
    if (ay < 0) return -1;
    return 0;
}
void MakeSimplePolygon(std::vector<CPoint>& pts, std::vector<CPoint>& spoly) {
    if (pts.size() < 1) return;
    int maxid = 0, minid = 0;
    for (int i = 1; i < pts.size(); i++) {
        if (pts[i].x > pts[maxid].x) maxid = i;
        if (pts[i].x < pts[minid].x) minid = i;
    }
    std::vector<CPoint> LP, RP;
    for (int i = 0; i < pts.size(); i++) {
        if (i == maxid || i == minid) continue;
        // 기준선의 왼편(LP)/오른편(RP)에 놓인 점 분리;
        if (CCW(pts[minid], pts[maxid], pts[i]) >= 0) LP.push_back(pts[i]); 
        else RP.push_back(pts[i]);
    }
    if (LP.size() > 0) 
        qsort(&LP[0], LP.size(), sizeof(CPoint), comp);
    if (RP.size() > 0)    
        qsort(&RP[0], RP.size(), sizeof(CPoint), comp);
	
    spoly.clear();
    spoly.push_back(pts[minid]);
    for (int i = 0; i < LP.size(); i++) spoly.push_back(LP[i]);
    spoly.push_back(pts[maxid]);
    for (int i = RP.size() - 1; i >= 0; i--) spoly.push_back(RP[i]);
    // spoly = clockwise simple polygon;
};

다른 방법으로는 위와 같은 직선을 만들고, 아랫부분의 점들의 convex hull을 구성하고, 나머지 점들을 마찬가지의 방법으로 정렬을 하여 폴리곤을 만들거나, 아니면 한 점을 잡아서(맨 아래 오른쪽) 그 점을 기준으로 해서 나머지 점들을 각도에 대해서 정렬을 한 다음 그 정렬된 순서대로 폴리곤을 구성하면 된다.

기준점에 대한 각도를 정렬하여서 만든 예(동일한 각도가 생기는 경우에는 기준점에서의 거리로 비교)

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

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

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
단순 다각형의 면적(2D)  (0) 2021.01.23
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

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

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

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

 

실제 계산에서는 다각형을 삼각 분할하지 않고도 간단한 방법에 의해서 무게중심을 구할 수 있다. 원점과 다각형의 각 변의 꼭짓점을 이용해서 삼각형들을 구성하면 원래의 다각형을 겹치게(원점이 내부에 있으면 겹침이 없다) 분할할 수 있다. 분할된 삼각형으로 무게중심을 구할 때 겹치는 영역의 기여를 제거해야 한다. 그런데 다각형 밖의 영역을 분할하는 삼각형은 다각형 내부를 분할하는 삼각형과는 다른 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에 무관하게 성립하는 공식이 된다.

POINT polygon_centroid(POINT 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;
        cx += (V[i].x + V[j].x) * tri_area2;
        cy += (V[i].x + V[j].x) * tri_area2;
        area2 += tri_area2;
    }
    cx /= 3 * area2;
    cy /= 3 * area2;
    return POINT(cx, cy)
};

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

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

단순 다각형 만들기  (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
Point in Polygon  (2) 2020.12.14
Posted by helloktk

댓글을 달아 주세요

728x90

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

 

세 점 ${\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;
}

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

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

단순 다각형의 Convex hull  (0) 2021.01.24
단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Incremental Delaunay Triangulation  (1) 2020.12.01
Posted by helloktk

댓글을 달아 주세요

728x90

일직선에 있지 않은 세 점 $A$, $B$, $C$가 만드는 외접원 (circumcircle)에 주어진 점 $D$가 포함되는가를 판별하는 것은 기하 알고리즘에서 매우 중요하다. 직접 외접원의 중심과 반지름을 구해서 판단을 할 수 있지만, 이 경우 수치 계산의 오류가 들어갈 여지가 매우 많다 (외접원의 중심과 반지름을 구하는 식을 만들어 보면 쉽게 알 수 있다). 여기서는 이런 문제에 부딪히지 않고 보다 robust 한 판단을 주는 식을 구해보자.

이는 공간 상의 네 점 $A(a_x, a_y, a_z)$, $B(b_x, b_y, b_z)$, $C(c_x, c_y, c_z)$, $D(d_x, d_y, d_z)$만드는 사면체의 (부호를 가지는) 부피를 구하는 공식이 다음 행렬식으로 주어진다는 사실에서 출발한다.

사용자 삽입 이미지

이 식을 이용하면 평면에서 3점의 외접원에 주어진 점이 포함이 되는가를 판별하는 식 역시 3차원 공간에서의 체적으로 표현이 가능하다.

$z= x^2 + y^2$로 표현이 되는 포물면(paraboloid)의 한 점 $(a, b, a^2 + b^2)$에서 접평면의 방정식을 구하면,

$$z=2ax + 2bx - (a^2 + b^2)$$

로 주어진다. 이 평면을 $z$-방향으로 $r^2$만큼 이동시키면

$$z=2ax + 2by -(a^2 +b^2) + r^2$$

로 주어지고, 이 평면이 포물면을 절단하는 곡선을 바닥에 정사영하면 중심이 $(a,b)$이고 반지름이 $r^2$인 원이 된다:

$$(x-a)^2 + (y-b)^2 = r^2$$

3점을 포물면 위로 올리면 3차원의 점이 되고, 이 3점이 만드는 평면과 포물면의 절단 곡선을 바닥으로 정사영하면 원래의 외접원이 된다. 따라서 외접원의 내부점은 평면 아래에, 외부점은 평면 위에 놓인다. 원주점은 평면에 놓인다. 동일한 (평면에 놓인) 바닥 삼각형에 대해서 꼭지점(테스트하려는 점)의 상대적 위치가 정반대에 놓이므로 행렬식으로 표현된 이들 4점이 만드는 사면체의 부호가 달라진다. 즉, 외접원 내부인가 외부인가를 포물면 위로 올려진 외접원 3점과 테스트하려는 점이 구성하는 사면체 체적의 부호로 판별이 가능하다

 

다시 정리하면, 세 점 $A(a_x, a_y)$, $B(b_x, b_y)$, $C(c_x, c_y)$, $D(d_x, d_y)$가 주어진 경우에 $A, B, C$(반시계 방향)에 의해서 형성이 되는 원에 점 $D$가 포함이 되는지의 여부는 공간상의 네 점

$$A = (a_x, a_y, a_x^2+a_y^2),$$

$$B = (b_x, b_y, b_x^2+b_y^2),$$

$$C = (c_x, c_y, c_x^2+c_y^2),$$

$$D = (d_x, d_y, d_x^2+d_y^2)$$

가 형성하는 체적의 부호에 의해서 결정이 된다.

사용자 삽입 이미지


위의 행렬식은 좀 더 간단히 정리할 수 있다. 1, 2, 3행에서 4행을 행단위로 빼도 값이 변하지 않으므로 오른쪽 식을 얻고(4열에 대해서 전개 후), 1열에 $2d_x$를 곱한 것과 2열에 $2d_y$를 곱한 것을 각각 3열에서 빼면 아래줄의 식을 얻는다.

사용자 삽입 이미지

이 3x3 행렬식은 쉽게 계산이 된다. $ccw(A, B, C)$를 계산하여 전체의 부호를 고정시킨다.

// A, B, C가 만드는 외접원에 D가 들어가는가?;
// det = 0인 경우는 세 점이 일직선이나, D가 원주에 놓인 경우다;
// 마지막 판단을 수정하여서 원주에 놓인 경우도 해결할 수 있다.
int incircle(POINT A, POINT B, POINT C, POINT D) {
    double ccw = CCW(A, B, C);
    double adx = A.x - D.x, ady = A.y - D.y,
           bdx = B.x - D.x, bdy = B.y - D.y,
           cdx = C.x - D.x, cdy = C.y - D.y,
           bdxcdy = bdx * cdy, cdxbdy = cdx * bdy,
           cdxady = cdx * ady, adxcdy = adx * cdy,
           adxbdy = adx * bdy, bdxady = bdx * ady,
           alift = adx * adx + ady * ady,
           blift = bdx * bdx + bdy * bdy,
           clift = cdx * cdx + cdy * cdy;
    double det = alift * (bdxcdy - cdxbdy)
               + blift * (cdxady - adxcdy)
               + clift * (adxbdy - bdxady);

    if (ccw > 0) return det > 0; // CCW인 경우에..
    else         return det < 0; // CW인  경우에..     
};

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

단순 다각형의 무게중심  (0) 2021.01.24
단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
Point in Polygon  (2) 2020.12.14
Incremental Delaunay Triangulation  (1) 2020.12.01
Chain Hull  (2) 2012.09.16
Posted by helloktk

댓글을 달아 주세요

728x90

기하 관련 프로그래밍을 하다 보면 평면 상의 한 점이 다각형 내의 점인가 아니면 밖의 점인가를 판단해야 할 필요성이 많이 생긴다. 이 문제에 대한 답을 생각해보기로 하자. $N$개의 꼭짓점 $\{ x_i, y_i| i=0,...., N-1\}$으로 구성이 된 다각형을 놓고 생각을 하도록 하자. 평면상의 주어진 점 $(x_p, y_p)$에서 오른쪽으로 출발하는 반직선을 쭉 그으면 다각형과 교차하는 점들이 생기게 된다. 교차하는 상황을 잘 헤아려 보면 해당점이 다각형의 내부에 들어 있는 경우에는 항상 홀수번의 교차점이 생기고, 외부에 있는 경우에는 짝수번(0번 포함)의 교차점이 생김을 알 수 있다. 따라서 주어진 점에서 시작하는 반직선과 다각형의 교차점 개수의 홀짝 여부를 판별하면 내부/외부점인지에 대한 문제를 풀 수 있다. 여기서 한 가지 주의할 사항은, 주어짐 점에서 그은 반직선이 다각형의 변(edge)을 포함하는 경우는 제외해야 하고, 꼭짓점이 반직선에 걸리는 경우는 꼭짓점을 공유하는 두 변과 만나지만 한 번으로 카운트해야 한다.

BOOL is_inside_polygon(POINT p, POINT polygon[], int N){
    int counter = 0;
    POINT* q1 = &polygon[0];
    for (int i = 1; i <= N; i++) {
        POINT *q2 = &polygon[i % N];
        if (p.y > min(q1->y, q2->y)) {                                //y-intersect.;
            if (p.y <= max(q1->y, q2->y)){                            //y-intersect(=single);
                if (p.x <= max(q1->x, q2->x)){                        //x-intersect;
                    if (q1->y != q2->y) {                             //not parallel;
                        double xq = (p.y - q1->y) * (q2->x - q1->x) / (q2->y - q1->y) + q1->x;
                        // x-intersect;
                        if (q1->x == q2->x || p.x <= xq)
                            counter++;
                    }
                }
            }
        }
        q1 = q2; // move to next edge;
    }
    if (counter % 2 == 0) return FALSE;
    else return TRUE;
}

네이버 블로그에서 이전

참고: 

1. www.ics.uci.edu/~eppstein/161/960307.html

2.geomalgorithms.com/a03-_inclusion.html

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

단순 다각형의 면적(2D)  (0) 2021.01.23
삼각형 외접원의 Inclusion Test  (0) 2020.12.30
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
Posted by helloktk

댓글을 달아 주세요

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

    꼭짓점과 반직선이 만나는 경우에 예외 처리가 제일 짜증나는 부분인 것 같습니다.
    1번으로 따로 카운트하는 방법도 있지만,
    반직선의 기울기를 아주 살짝 위를 향하게 하는 방법도 있다고 합니다!

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

    혹은 모든 점에 대해 ccw를 이용하는 방법도 있다고는 하는데,
    솔직히 아무도 안 쓸 것 같네요

728x90

// 수퍼 삼각형은 입력 점집합을 모두 포함하도록 충분히 큰 삼각형을 선택한다;
// 삼각화한 후 최외각이 convex가 아니면 수퍼 삼각형을 더 크게 다시 잡아야 함.
// 마지막 수정: 2021.4.7, 코드 간결화;
#define REMOVE_FLAG    (-1)
//
// max num of triangles = 2 * (N - 2) - 1;
// max num of edegs     = 3 * (N - 2);
int inc_delaunay_triangulation(POINT P[], int N, TRIANGLE *V) {
    // 메모리는 입력 점집합의 크기가 (N+3)인 경우를 기준으로 할당되어야 한다;
    int emax = 3 * (N + 1) ;
    EDGE *E = (EDGE *)malloc(sizeof(EDGE) * emax) ;
    set_supertriangle(P, N) ;
    // First element of triangle index is supertriangle vertex;
    push_tri(&V[0], N, N + 1, N + 2);
    int nt = 1;
    for (int i = 0; i < N; i++) {
        int ne = 0;
        /* 추가점(P[i])를 외접원에 포함시키는 삼각형을 모두 찾아 제거한다.
        ** 이들 삼각형의 공유되는 내부 edge를 제외한 나머지 외부 edge와
        ** 추가점으로 새로운 삼각형을 만들어 리스트에 추가한다.
        */ 
        for (int j = 0; j < nt; j++) {
            TRIANGLE *t = &V[j];
            if (incircle(P[t->a], P[t->b], P[t->c], P[i])) {
                // 제거되는 삼각형의 edge는 백업한다; 내부 edge는 중복된다.
                push_edge(&E[ne++], t->a, t->b);
                push_edge(&E[ne++], t->b, t->c);
                push_edge(&E[ne++], t->c, t->a);
                // P[i]를 포함하는 외접원 삼각형은 제거;
                // swap(j, nt-1) to remove triangle at j
                V[j--] = V[--nt];  
            }
        }
        // 중복되어 있는 내부 edge를 제거;
        remove_double_edges(E, ne);
        
        // 외부 edge와 P[i]가 만드는 새 삼각형을 리스트에 추가;
        for (int j = 0; j < ne; j++) {
            if (is_removed_edge(&E[j])) continue;
            push_tri(&V[nt++], E[j].a, E[j].b, i);
        }
    }
    //done with super_triangle;
    free(E);
    return remove_supertriangle_relics(N, V, nt) ;
}
void remove_double_edges(EDGE *E, int ne) {
    for (int j = 0; j < ne - 1; j++)
    	for (int k = j + 1; k < ne; k++)
            if (is_double_edge(&E[j], &E[k])) {
            	pop_edge(&E[j]); pop_edge(&E[k]);
            }
}
void push_tri(TRIANGLE *V, int a, int b, int c) {
    V->a = a; V->b = b; V->c = c;
}
int is_removed_edge(const EDGE *E) {
    return (E->a == REMOVE_FLAG || E->b == REMOVE_FLAG);
};
int is_double_edge(const EDGE* E1, const EDGE* E2) {
    return ((E1->a == E2->a && E1->b == E2->b) ||
            (E1->a == E2->b && E1->b == E2->a)) ;
}
void push_edge(EDGE *E, int a, int b) {
    E->a = a; E->b = b;
}
void pop_edge(EDGE* E) {
    E->a = E->b = REMOVE_FLAG;
}

// 삼각형 중에서 수퍼 삼각형의 꼭지점(N,N+1,N+2)을 가지고 있는 것을 제거;
int remove_supertriangle_relics(int N, TRIANGLE *V, int nt) {
    for (int i = 0; i < nt; i++)
        if (V[i].a >= N || V[i].b >= N || V[i].c >= N)
            // swap(i, nt - 1) to remove triangle at j
            V[i--] = V[--nt];
    return nt;
}

// A,B,C가 만드는 외접원에 D가 들어가는가?;
int incircle(POINT A, POINT B, POINT C, POINT D) {
	POINT AD = A - D, BD = B - D, CD = C - D;
    double ccw = CCW(A, B, C);
    double det = norm2(AD) * cross(BD, CD)
               + norm2(BD) * cross(CD, AD)
               + norm2(CD) * cross(AD, BD)
    if (ccw > 0) return det > 0; //CCW인 경우에..
    else         return det < 0; //CW인  경우에..    
}

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

삼각형 외접원의 Inclusion Test  (0) 2020.12.30
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
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2020.12.02 20:40 신고  댓글주소  수정/삭제  댓글쓰기

    그 개어렵다는 델로네 삼각분할이군요...
    나중에 배울때 참고하겠습니다^^

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

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. For~~\delta [k] = 0 \quad\Longrightarrow\quad m [k]=m [k+1]=0$$

$$2. 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();
    int n = points.size() ;
    if (n < 2) return ;
    //if (!isMonotonicPolyline(points)) return ;
    std::vector<double> xx(n), yy(n);
    for (int i = 0; i < n; i++) {
        xx[i] = points[i].x ;
        yy[i] = points[i].y;
    };
    //tangent table;
    std::vector<double> m(n, 0);
    estimateTangents(xx, yy, m);
    //hermite spline;
#define STEPS 12
    curve.push_back(points[0]); //add first point;
    for (int i = 0; i < n - 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과의 비교:

 

 

 

 

 

 

 

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

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
Douglas-Peucker Algorithm  (0) 2012.05.03
Posted by helloktk

댓글을 달아 주세요

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

댓글을 달아 주세요