Geometry & Recognition

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

한 개의 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...$$

두 개의 점이 주어지는 경우 이들을 잇는 가장 단순한 곡선은 직선이다. 또 직선은 두 점을 잇는 가장 짧은 거리를 갖는 곡선이기도 하다. 그러면 $N$개의 점이 주어지는 경우는 이들 모두 지나는 곡선을 어떻게 찾을까? 구간별로 직선으로 연결을 시켜서 하나의 곡선을 구성할 수 있으나 부드럽게 이어지는 형태의 곡선은 아니다. 곡선이 smooth 함은 그 곡선의 곡률이 잘 정의된다는 의미이다(곡률은 곡선의 이차 미분에 의존한다). 주어진 점들을 연결하는 충분히 짧고 smooth 한 곡선을 찾는 문제는 곡률에 의한 에너지를 최소화시키는 해를 구하는 문제로 환원할 수 있다. 이 문제의 해는 잘 알려진 cubic spline이다. 그러나 cubic spline은 지나는 점(컨트롤 점)의 변동에 대해서 안정적이 아니다. 컨트롤 점 중에서 한 점만 변해도 곡선이 전역적(global)으로 변하는 특성을 나타낸다. 컨트롤 점의 변동이 있을 때 그 점 주위에서만 변화하는 국소 특성을 갖는 곡선을 찾기 위해서는 smoothness 조건을 완화해야 한다. 그렇지만 충분히 부드러운 곡선이어야 하므로 적어도 일차의 미분 값의 연속성 정도는 보장되어야 한다. 이런 특성을 갖는 삼차 곡선이 Catmull-Rom spline이다. 이름은 Edwin Catmull와 Raphie Rom의 이름에서 유래한다.

특징: 

          주어진 점들 위에서 정의됨.

          국소적인 변동 특성(노이즈에 안정적).

          일차 미분의 연속성. 

두 점 $P_{0}$와 $P_{1}$을 지나는 일반적인 3차의 곡선을 적으면

$$ P(t) = a_{0}  + a_{1}  t + a_{2}  t^{2} + a_{3}  t^{3}$$

이 곡선이 $t=0$에서 $P_{0}$을 지나고, $t=1$에서 $P_{1}$을 지난다고 하더라도, 완전히 결정되지 않는다. 따라서 추가적인 조건을 주어야 하는데, 여기서는 $P_{0}$와 $P_{1}$에서의 미분 값이 주어진다고 가정한다. 즉,

$$P(0) = P_{0},\quad P(1) = P_{1},\quad  P'(0) = P_{0}',\quad P'(1) = P_{1}'$$

이 조건에서 계수들은

$$a_0 = P_0, \quad   a_1 = P_0', \\a_2 = 3(P_1- P_0) - 2P_0' - P_1', \\ a_3 = 2(P_0- P_1) + P_0' + P_1' $$

로 주어진다;

$$P(t)=\left[\begin{array}{cccc}1 & t &t^2& t^3\end {array}\right]\left [\begin {array}{rrrr}1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1\end {array}\right]\left [\begin {array}{r} P_0\\P_1\\P_0'\\P_1'\end {array}\right]$$

그러나  미분 값을 직접 정해서 넣는 것은 문제가 있다.

4 이상의 점들이 주어지는 경우에는 컨트롤상에서의 미분 값을 그 점 주위의 두 점을 있는 직선의 기울기 값으로 근사 시킬 수 있다:

$$P_i' \longrightarrow (P_{i+1}-P_{i-1})/2$$

이 미분 값을 사용하면 네 개의 컨트롤 점 $\{P_{i-1}, P_i, P_{i+1}, P_{i+2}\}$이 주어진 경우 $P_i (\leftarrow t=0)$와 $P_{i+1}(\leftarrow t=1)$ 사이를 보간하는 곡선은

$$\begin {align} P(t)&=\left [\begin {array}{cccc}1 & t &t^2& t^3\end {array}\right] \left[\begin{array}{rrrr}1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1\end{array}\right] \left[\begin{array}{c}P_i\\P_{i+1}\\(P_{i+1}-P_{i-1})/2\\(P_{i+2}-P_{i})/2\end{array}\right]\\ &=\frac{1}{2}\left[\begin{array}{cccc}1 & t &t^2& t^3\end{array}\right] \left [\begin {array}{rrrr}0&2&0&0\\-1&0&1&0\\2&-5&4&-1\\-1&3&-3&1\end {array}\right] \left [\begin {array}{c} P_{i-1}\\P_{i}\\P_{i+1}\\P_{i+2}\end {array}\right] \end {align}$$

로 표현된다. 다시 정리하면

$$ P(t) = \frac{1}{2}\left[ 2P_i + (-P_{i-1}+P_{i+1})t+ (2P_{i-1}-5P_i +4P_{i+1} - P_{i+2})t^2 \\ + (-P_{i-1} + 3P_i -3P_{i+1} +P_{i+2}) t^3 \right]$$

열린 곡선으로 보간을 하는 경우 처음과 끝 구간에는 미분 값을 구할 수 없으므로 정의가 안되지만, 나머지 구간에서는 주어진 점들을 연결하는 충분히 부드러운 곡선이 된다. 양 끝 구간의 보간은 동일점을 반복함으로써 전체적인 구간에서 잘 정의되게 만들 수 있다(끝점에서 기울기를 그 점과 인접점과의 기울기/2로 잡는 것과 같은 효과임). 닫힌 곡선인 경우에는 모든 구간에서 잘 정의가 된다. Catmull-Rom spline은 Bezier나 B-Spline 곡선의 경우와 다르게 컨트롤 점의 convex hull 내부에서만 정의되는 특성을 따르지 않는다.

CPoint CRSpline(double t, CPoint p1, CPoint p2, CPoint p3, CPoint p4) {
    double tt = t * t ;
    double ttt = tt * t ;
    double x = 0.5 * ((-p1.x + 3 * p2.x - 3 * p3.x + p4.x) * ttt
        + (2 * p1.x - 5 * p2.x + 4 * p3.x - p4.x) * tt
        + (-p1.x + p3.x) * t
        + 2 * p2.x);
    double y = 0.5 * ((-p1.y + 3 * p2.y - 3 * p3.y + p4.y) * ttt
        + (2 * p1.y - 5 * p2.y + 4 * p3.y - p4.y) * tt
        + (-p1.y + p3.y) * t
        + 2 * p2.y);
    return CPoint(int(x + .5), int(y + .5)) ;
}

//open spline의 drawing;
void DrawCatmullRom(std::vector<CPoint>& Q, CDC* pDC) {  
#define STEPS (20)
    if (Q.size() < 4) return ;
    CPen red(PS_SOLID, 1, RGB(0xFF, 0, 0));
    CPen *pOld = pDC->SelectObject(&red);
    const int n = Q.size();
    for (int i = 0; i < n - 1; i++) {
        pDC->MoveTo(Q[i]);
        // i = 0 인 경우에는 처음 점을 반복, i = n - 2인 경우에는 마지막 점을 반복..
        int ip = max(i - 1, 0);
        int inn = min(i + 2, n - 1);
        for (int it = 1; it <= STEPS; it++)
            pDC->LineTo(CRSpline(double(it)/STEPS, Q[ip], Q[i], Q[i + 1], Q[inn]));
    };
    pDC->SelectObject(pOld);
}

**네이버 블로그 이전;

히스토그램은 영상에서 각 그레이 값에 해당하는 픽셀의 수를 주는 일종의 이산적 함수로 생각할 수 있다. histogram에서 피크의 위치는 histogram을 연속적인 함수로 모델링하거나 또는 여러 개의 그룹으로 분리를 할 때 중요한 정보를 제공한다. 영상으로부터 얻은 histogram은 대부분 이웃하는 그레이 값 사이에서 smooth 하게 변하지 않기 때문에 피크를 찾는 작업을 하기 전에 미리 mean filter나 gaussian filter와 같은 smoothing 과정을 거친 후 사용한다. 여기서는 low-pass filter 대신에 histogram의 bin 인덱스와 bin 값을 컨트롤 포인트로 사용해서 만든 Bezier 곡선을 이용해서 histogram을 smooth 한 곡선으로 근사하는 방법을 알아본다. 이 경우 Bezier 곡선은 255-차수의 곡선이 된다. 높은 차원의 Bezier 곡선 계산에 Berstein 함수를 사용하는 경우 truncation 등의 수치 에러 때문에 값이 불안정해지므로 De Casteljau's algorithm을 이용하여서 iterative 하게 값을 계산을 하면 된다.

// De Casteljau's algorithm (degree:=size(Q)-1)
double Bezier(int deg, double Q[], double t) {
    for (int k = 0; k < deg; k++)
        for (int j = 0; j < (deg - k); j++)
            Q[j] = (1 - t) * Q[j] + t * Q[j + 1];
    return Q[0];
}

void SmoothenHistogram (int hist[], int numLevels/* =256 */) {
    std::vector<double> p(numLevels);
    std::vector<int> hist2(numLevels);
    // cloning;
    for (int j = 0; j < numLevels; j++) hist2[j] = hist[j];
    for (int j = 0; j < numLevels; j++) {
        double t = double(j) / (numLevels - 1);
        // control points {p}; calling of Bezier() modifies p's;
        for (int i = 0; i < numLevels; i++) 
            p[i] = hist2[i];
        hist[j] = int(Bezier(numLevels-1, &p[0], t) + 0.5);
    }
};

data=violet,  smoothed=green;

Bezief curve 위키피디아: http://en.wikipedia.org/wiki/B%C3%A9zier_curve