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 point가 만드는 convex region 내부에 있음을 의미한다. Bezier 곡선의 convexity 성질은 여러 좋은 특성을 준다.  몇 가지만 나열하면, 첫째가 Bezier 곡선은 항상 컨트롤 포인트의 convex hull 내에 놓이게 되므로 곡선의 제어 가능성을 보장한다. 둘째는 교차 여부를 쉽게 확인할 수 있다. 또한 culling을 가능하게 한다.

 

 
728x90

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

Bezier Smoothing  (0) 2021.04.23
Flatness of Cubic Bezier Curve  (0) 2021.04.23
De Casteljau's Algorithm  (0) 2021.04.22
Arc Length of Bezier Curves  (0) 2021.04.21
Bezier Curve Approximation of an Ellipse  (0) 2021.04.11
Posted by helloktk
,

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

// De Casteljau's algorithm; recursive version; slow for larger n;
double Bezier(double t, int 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, int 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 segments, std::vector<CfPt> &curves) {
    std::vector<double> xp(cntls.size()), xp(cntls.size());
    curves.resize(segments + 1);
    for (int i = 0; i <= segments; ++i) {
        double t = double(i) / segments;
        // clone control points; non-rec version modifies inputs;
        for (int k = cntls.size(); k-->0;) {
            xp[k] = cntls[k].x; yp[k] = cntls[k].y;
        }
        curves[i] = CfPt(Bezier(t, xp.size(), &xp[0]), Bezier(t, yp.size(), &yp[0]));
    }
}

 
 
 
 
 
 
728x90

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

Flatness of Cubic Bezier Curve  (0) 2021.04.23
Convexity of Bezier Curve  (0) 2021.04.22
Arc Length of Bezier Curves  (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
,

곡선의 길이를 구하고자 하면 곡선을 따라 속력을 적분하면 된다. 곡선을 기술하는 벡터가 $\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 {cord length} + \sum\text {distance between control points}}{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

다음은 4분원의 길이를 위의 알고리즘을 이용해서 구하는 예제이다. 네 개의 control point는

$$ Q_1=(0,1), ~~Q_2 = (k, 1),~~Q_3=(1, k),~~Q_4=(1,0)$$

이고, $k$ 값은 Bezier curve의 중간이 원에 접하는 경우, 원과 벗어남이 가장 작은 경우, 면적을 같게 하는 경우, 길이를 같게 하는 경우 각각에 대해 계산을 한다. 

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

 

 
728x90

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

Convexity of Bezier Curve  (0) 2021.04.22
De Casteljau's Algorithm  (0) 2021.04.22
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
,

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);
}
728x90

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

De Casteljau's Algorithm  (0) 2021.04.22
Arc Length of Bezier Curves  (0) 2021.04.21
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
,

한 개의 Bezier 곡선을 이용해서 원을 표현할 수 없음은 잘 알려진 사실이다. 그럼 Bezier 곡선을 이용해서 얼마나 원(호)을 잘 근사할 수 있을까? 원점에 중심을 둔 반지름 1인 원의 1 사분면 원호가 3차 Bezier 곡선으로 얼마나 잘 표현되는지 알아보자. 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$)과 끝점($t=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)$$

에서 나머지 두 control point는 다음과 같이 쓸 수 있다:

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

그럼 $k$ 값은 어떻게 정할 수 있을까?

(2-1) 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...$$

728x90

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

Arc Length of Bezier Curves  (0) 2021.04.21
Bezier Curve Approximation of an Ellipse  (0) 2021.04.11
Bresenham's Line Algorithm  (0) 2021.04.07
Rotating Calipers  (3) 2021.03.31
Convex Hull Peeling  (0) 2021.03.29
Posted by helloktk
,