Geometry & Recognition

주어진 데이터를 interpolation을 할 때 3차 spline을 많이 사용한다. 그럼 왜 3차 일까? 데이터를 연결하는 spline은 되도록이면  부드럽게 연결되어야 한다. 곡선이 부드럽게 그려지기 위해서는 급격한 꺾임이 없어야 하는데 얼마나 급격히 꺾이는가는 곡률에 비례하고, 곡률은 그 지점에서 함수의 두 번 미분값에 비례한다.(곡선의 곡률: https://kipl.tistory.com/105) 따라서 어떤 함수 $f(x)$가 구간 $[a,b]$에서 얼마나 부드럽게 연결되는가에 대한 척도는 곡선의 각부분에서 곡률의 크기의 제곱을 를 더한  다음 (음수가 아닌) 적분이 제공할 수 있다.

$$\kappa [f]={\int }_a^b(f^{″}(x){)}^{2}dx$$

$a\le x\le b$ 구간에서 균일하게 샘플링된 데이터가 있고 이들을 보간하는  3차 spline이 $S(x)=\{S_i(x)=a_i{x}^3+b_i{x}^2+c_{i}x+d_i\}$라고 하자. 또 두 번 이상 미분가능한 임의의 함수 $f(x)$도 주어진 샘플링 데이터를 지나가는 보간함수라고 하자. 이 경우 natural cubic spline이 주어진 sampling 데이터를 지나가면서 가장 부드럽게 이어지는 곡선임을 보일 수 있다.

$$\kappa [f]\ge \kappa [S]\mathrm{\forall }f$$

Spline $S$와 함수 $f$의 차이를 $h(x)=f(x)-S(x)$라 하면 각 node에서 $h(x_i)=0$이다. 이제 

$$\kappa [f]={\int }_a^b(h^{″}(x)-S^{″}(x){)}^2=\kappa [h]+\kappa [S]-2{\int }_a^b{h}^{″}(x)S^{″}(s)dx$$

그런데, cubic spline의 경우 각 node에서 2차 미분이 연속이므로  

$$\begin{array}{rl}{\int }_a^b{h}^{″}(x)S^{″}(x)dx& =h^{\prime }(x)S^{″}(x){|}_a^b-{\int }_a^b{h}^{\prime }(x)S^{‴}(x)dx\\ & =h^{\prime }(b)S^{″}(b)-h^{\prime }(a)S^{″}(b)-\sum _i{\int }_{x_i}^{x_{i+1}}{h}^{\prime }(x)S_i^{‴}(x)dx\\ & =h^{\prime }(b)S^{″}(b)-h^{\prime }(a)S^{″}(a)-\sum _i{\int }_{x_i}^{x_{i+1}}{h}^{\prime }(x)(6a_i)dx\\ & =h^{\prime }(b)S^{″}(b)-h^{\prime }(a)S^{″}(a)-6\sum _i{a}_i[h(x_{i+1})-h(x_i)]\\ & =h^{\prime }(b)S^{″}(b)-h^{\prime }(a)S^{″}(a)\end{array}$$  $$\therefore \kappa [f]=\kappa [h]+\kappa [S]+h^{\prime }(b)S^{″}(b)-h^{\prime }(a)S^{″}(a)$$

임의의 곡선 $h$에 대해서 $\kappa [h]\ge 0$이므로, 양끝에서 2차 미분이 $S^{″}(a)=S^{″}(b)=0$으로 고정된 natural cubic spline은  

$$\kappa [f]\ge \kappa [S]$$

이므로 주어진 샘플링 데이터를 곡률척도가 가장 작게 즉, 가장 부드럽게 보간하는 곡선임을 알 수 있다.

구간 $[a=t_0,...t_{n-1}=b]$에서 일정한 간격(꼭 일정한 간격일 필요는 없지만 여기서는 1로 선택함)으로 샘플링된 데이터 $\{f_0,...f_{n-1}\}$이 있다. $n-1$ 개의 각 구간에서 데이터를 보간하는 삼차다항식 함수 모임 $$S(t)=\{S_j(t)|t_j\le t<t_{j+1},j=0,1,2...,n-2\}$$을 찾아보자. 전체 구간에서 연속적이고 충분히 부드럽게 연결되기 위해서는 우선 각 node에서 연속이어야 하고, 또한 1차 도함수와 2차 도함수도 연속이 되어야 한다.  물론 각 node에서는 샘플링된 데이터값을 가져야 한다. 

$$\begin{array}{rl}(a)& S(t_j)=f_j\\ (b)& S_{j+1}(t_{j+1})=S_j(t_{j+1})\\ (c)& S_{j+1}^{\prime }(t_{j+1})=S_j^{\prime }(t_{j+1})\\ (d)& S_{j+1}^{″}(t_{j+1})=S_j^{″}(t_{j+1})\end{array}$$

$n-1$개 구간에서 각각 정의된  3차 함수를 결정하려면 $4\times (n-1)$개의 조건이 필요하다. (a)에서 $n$개, (b), (c), (d)에서 $3\times (n-2)$개의 조건이 나오므로 총 $n+3\times (n-2)=4n-6$개의 조건만 생긴다. 삼차식을 완전히 결정하기 위해서는 2개의 추가적인 조건을 부여해야 하는데 보통 양끝에서 2차 도함수값을 0으로 하거나(natural boundary) 또는 양끝에서 도함수 값을 특정한 값으로 고정시킨다(clamped boundary). 여기서는 양끝의 2차 도함수를 0으로 한 natural cubic spline만 고려한다. 그리고 $S_j(t)$가 $n-2$개의 구간에 대해서만 정의되어 있지만, 끝점을 포하는 구간에서도 정의하자. 이 경우 $S_{n-1}(t),t\ge t_{n-1}$에 부여된 조건은 $t_{n-1}$에서 $S_{n-2}(t)$와 연속, 미분연속 그리고 2차 도함수가 0인 조건만 부여된다.

$j-$번째 구간의 삼차함수를 

$$S_j(t)=a_j+b_j(t-t_j)+c_j(t-t_j{)}^2+d_j(t-t_j{)}^3$$

로 쓰면 (a)에서 

$$S_j(t_j)=a_j=f_i,j=0,1,..n-2$$

(b)에서 

$$a_{j+1}=a_j+b_j+c_j+d_j,j=0,1,...,n-2$$

(c)에서 

$$b_{j+1}=b_j+2c_j+3d_j,j=0,1,...,n-2$$

(d)에서 

$$c_{j+1}=c_j+3d_j,j=0,1,...,n-2$$

이므로 (b)와 (c)에서 $d_j$을 소거하면

$$c_j=3(a_{j+1}-a_j)-b_{j+1}-2b_j$$

그리고 (a)에서 $$d_j=-2(a_{j+1}-a_j)+b_j+b_{j+1}$$ 이므로 $b_j$에 대해 정리하면 다음과 같은 점화식을 얻을 수 있다.

$$b_{j+1}+4b_j+b_{j-1}=3(a_{j+1}-a_{j-1})=3(f_{j+1}-f_j),j=1,2,...,n-2$$

물론 $j=0$일 때는 (note, $S_0^{″}(t_0)=0\to c_0=0$) $a_1=a_0+b_0+d_0,b_1=b_0+3d_0$이므로

$$b_1+2b_0=3(f_1-f_0)$$

$j=n-1$일 때 계수는 $S_{n-1}(t_{n-1})=f_{n-1}$이고 $S_{n-1}^{″}(t_{n-1})=c_{n-1}=0$이므로 $$2b_{n-1}+b_{n-2}=3(f_{n-1}-f_{n-2})$$

임을 알 수 있다. 따라서 계수를 구하는 과정은 $b_j(j=0,...,n-1)$을 구하는 것으로 결정된다. 이를 행렬로 표현하면

$$\left[\begin{array}{ccc}2& 1& \\ 1& 4& 1\\ & 1& 4& 1\\ & & & \cdots \\ & & & & 1& 4& 1\\ & & & & & 1& 2\end{array}\right]\left[\begin{array}{c}{b}_0\\ b_1\\ ⋮\\ b_{n-2}\\ b_{n-1}\end{array}\right]=\left[\begin{array}{c}3(f_1-f_0)\\ 3(f_2-f_0)\\ ⋮\\ 3(f_{n-1}-f_{n-3})\\ 3(f_{n-1}-f_{n-2})\end{array}\right]$$

와 같다. band 행렬은 upper triangle로 변환한 후 역치환과정을 거치면 쉽게 해를 구할 수 있다.

평면에서 주어진 점들을 보간하는 곡선은 이들 점을 표현하는 곡선의 매개변수를 일정한 간격으로 나누어서 샘플링된 결과로 취급하면, $x,y$ 성분에 대해서 각각 natural cubic spline를 구하여 얻을 수 있다. 그러나 곡선이 급격히 변할 때는 매개변수를 일정하게 잡는 것보다는 입력점의 사이거리를 기준으로 선택하는 것이 더 유연한 보간곡선을 만들어 준다.

struct Cubic {
    double a,b,c,d;  /* a + b*t + c*t^2 +d*t^3 */
    Cubic() {}
    Cubic(double a_, double b_, double c_, double d_) 
        : a(a_), b(b_), c(c_), d(d_) { }
    /* evaluate;*/
    double eval(double t) {
        return (((d*t) + c)*t + b)*t + a;
    }
};
std::vector<Cubic> calcNaturalCubic(std::vector<double>& x) {
    std::vector<double> gamma(x.size());
    std::vector<double> delta(x.size());
    std::vector<double> D(x.size());
    /* solve the banded equation:
    [2 1       ] [  D[0]]   [3(x[1]   - x[0])  ]
    |1 4 1     | |  D[1]|   |3(x[2]   - x[0])  |
    |  1 4 1   | | .    | = |      .           |
    |    ..... | | .    |   |      .           |
    |     1 4 1| | .    |   |3(x[N-1] - x[N-3])|
    [       1 2] [D[N-1]]   [3(x[N-1] - x[N-2])]
    ** make the banded matrix to an upper triangle;
    ** and then back sustitution. D[i] are the derivatives at the nodes.
    */
    const int n = x.size() - 1;  // note n != x.size()=N;
    // gamma;
    gamma[0] = 0.5;
    for (int i = 1; i < n; i++)
        gamma[i] = 1/(4-gamma[i-1]);
    gamma[n] = 1/(2-gamma[n-1]);
    // delta;
    delta[0] = 3*(x[1]-x[0])*gamma[0];
    for (int i = 1; i < n; i++) 
        delta[i] = (3*(x[i+1]-x[i-1])-delta[i-1])*gamma[i];
    delta[n] = (3*(x[n]-x[n-1])-delta[n-1])*gamma[n];
    // D;
    D[n] = delta[n];
    for (int i = n; i-->0;)
        D[i] = delta[i] - gamma[i]*D[i+1];

    /* compute the coefficients;*/
    std::vector<Cubic> Spline(n);
    for (int i = 0; i < n; i++)
        Spline[i] = Cubic(x[i], D[i], 3*(x[i+1]-x[i])-2*D[i]-D[i+1],
                                      2*(x[i]-x[i+1]) + D[i] + D[i+1])) ;
    return Spline;
}
std::vector<CPoint> NaturalCubicSpline(const std::vector<CPoint>& points) {
    if (points.size() < 2) return std::vector<CPoint> (); // null vector;
    std::vector<double> xp(points.size()), yp(points.size());
    for (int i = points.size(); i-->0;)
        xp[i] = points[i].x, yp[i] = points[i].y;
        
    std::vector<Cubic> splineX = calcNaturalCubic(xp);
    std::vector<Cubic> splineY = calcNaturalCubic(yp);
#define STEPS 12
    std::vector<CPoint> curve;
    curve.reserve(splineX.size() * STEPS + 1);
    curve.push_back(points.front());
    for (int i = 0; i < splineX.size(); i++) {
        for (int j = 1; j <= STEPS; j++) {
            double t = double(j) / STEPS;
            curve.push_back(CPoint(int(splineX[i].eval(t)+0.5), int(splineY[i].eval(t)+0.5)));
        }
    }
    return curve;
}

타원은 베지어 곡선을 이용해서 근사적으로 표현할 수 있음을 잘 알고 있다(https://kipl.tistory.com/313). 표준형 타원의 경우 1사분면에서 표현만 알면 대칭성에 의해서 나머지 사분면에서는 바로 구할 수 있으므로 거리를 구하는 점이 1사분면에 있는 경우만 알아보면 충분하다. 장축 반지름이 $a$, 단축 반지름이 $b$인 타원은 다음과 같이 베이지 곡선으로 표현된다.

$$\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}ax(t)\\ by(t)\end{array}\right)$$

$$x(t)=3k(1-t{)}^{2}t+3(1-t)t^2+t^3$$

$$y(t)=x(1-t),0\le t\le 1$$

$k$ 값은 $t=1/2$일 때 $(a/\sqrt{2},b/\sqrt{2})$을 통과하는 조건을 부여하여 정하면

$$k=\frac{4}{3}(\sqrt{2}-1)=0.5522847498...$$임을 알 수 있다. 

한 점 $\mathbf{P}=(p,q)$에서 베지어 곡선 위의 한 점 $\mathbf{B}(t)$까지의 거리가 최단거리가 되기 위해서는 다음 직교 조건

$$f(t)=(\mathbf{B}(t)-\mathbf{P})\cdot {\mathbf{B}}^{\prime }(t)=0$$

$$f(t)=a^{2}x(t)x^{\prime }(t)-b^{2}y(t)y^{\prime }(t)-ap{x}^{\prime }(t)-bq{y}^{\prime }(t)=0$$

을 만족하는 근을 찾으면 된다. 풀어야 할 방정식이 5차이므로 직접적으로 근을 찾는 것은 불가능하므로 Newton-Raphson 방법을 사용하도록하자. 앞선 정확한 거리 계산에서와는 달리 초기 $t$ 값 설정에 민감하게 의존하지 않는다.

//cubic bezier_x: P0(0, 1), P1(k, 1),P2(1, k), P3(1, 0);
double B(double t) { //t in [0:1]
    const double k = 4.*(sqrt(2.)-1)/3;
    return t*(3*k + t*(3 - 6*k + (-2 + 3*k)*t));
}
// derivative of B(t);
double DB(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 3*k + t*(6 - 12*k + (-6 + 9*k)*t);
}
// derivative of DB(t);
double D2B(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 6 - 12*k + (-12 + 18*k)*t;
}
// ellipse radii=(a, b);
double dist2EllipseBezier3(double p, double q, double a, double b,
                           double& xt, double& yt) {
    if (a == b) return dist2Circle(p, q, a, xt, yt);
    double x = fabs(p), y = fabs(q);
    const double eps = 0.001 / max(a, b);
    double t = 0.5;  // mid
    while (1) {
        // Newton-Raphson;
        double f = a*a*B(t)*DB(t)-b*b*B(1-t)*DB(1-t)-a*x*DB(t)+b*y*DB(1-t);
        if (f == 0) break;
        double df = a*a*(DB(t)*DB(t)+B(t)*D2B(t))+b*b*(DB(1-t)*DB(1-t)+B(1-t)*D2B(1-t))
                    -a*x*D2B(t)-b*y*D2B(1-t);
        double dt = f / df;
        t -= dt;
        t = max(0, min(1, t));
        if (abs(dt) < eps ) break;
    }
 
    xt = a * B(t); yt = b * B(1-t);
    xt = p >= 0? xt: -xt;
    yt = q >= 0? yt: -yt;
    return hypot(p - xt, q - yt);
}