Geometry & Recognition

단위 정사각형의 네 꼭짓점에서 함수값이 주어지고, 꼭짓점의 미분값을 주변 격자점으로부터 추정할 수 있을 때 내부의 임의의 점에서 함수값을 bicubic interpolation으로 구해 보자. Bicubic 보간함수 $x, y$ 각각에 대해 3차인 다항식

$$ f(x, y) = \sum_{i=0}^3 \sum_{j = 0}^3 a_{ij} x^i y^j =\begin {bmatrix}1 & x & x^2 & x^3 \end {bmatrix} \begin {bmatrix} a_{00} & a_{01} & a_{02} & a_{03} \\ a_{10} & a_{11} & a_{12} & a_{13} \\ a_{20}& a_{21}& a_{22} & a_{23}\\ a_{30} & a_{31} & a_{32} & a_{33}\end {bmatrix} \begin {bmatrix} 1 \\ y \\ y^2 \\ y^3  \end {bmatrix}$$으로 나타낼 수 있다. 16개의 계수 $a_{ij}$는 4 꼭짓점에서 4개의 함숫값, 8개의 1차 편도함숫값, 그리고 4개의 교차편도함수값으로 결정된다. https://kipl.tistory.com/55

그런데 이 bicubic interpolation을 $x$와 $y$ 방향을 동일한 방식으로 취급하므로, 보간식을 다음과 같이 두 방향의 대칭적인 구조가 드러나는 형태로 쓸 수 있다.

$$ f(x, y)  = \bf (M^T X) ^T P (M^T Y) = X^T M P M^T Y  $$

$$ \bf X= \begin{bmatrix} 1 & x & x^2 & x^3\end {bmatrix}^T,\quad Y = \begin {bmatrix} 1 & y & y^2 & y^3\end {bmatrix}^T$$ 이 경우 행렬 $\bf M$은 $x$나 $y$ 위치에 의존하지 않는 상수행렬이어야 하고, 행렬 \(\bf P\)는 네 꼭짓점에서의 함수값과 미분값이 들어가도록 형태를 가질 것이다.

행렬 $\bf M$과 $\bf P$의 형태를 추측하기 위해서 먼저 1차원 cubic interpolation을 고려하자. $0 \le t \le 1$ 구간의 양끝에서 값과 미분값이 주어진 경우 구간에서 3차 보간함수는 다음과 같이 표현된다.

$$ p(t) = \bf \begin{bmatrix} 1 & t & t^2 & t^3 \end {bmatrix} M \begin {bmatrix}  p(0) \\p(1) \\p'(0) \\ p'(1)\end {bmatrix}, \qquad M = \begin {bmatrix} 1& 0&0&0 \\ 0 & 0 & 1 & 0 \\ -3 & 3 & -2 & -1\\ 2 & -2 & 1 & 1 \end {bmatrix}$$구체적으로 계산을 하면

$$\bf \begin{bmatrix}1&t&t^2&t^3\end{bmatrix}M = \begin{bmatrix}2t^3-3t^2+1\\ -2t^3+3t^2\\ t^3-2t^2+t \\ t^3-t^2\end{bmatrix}^T$$가 되며, 네 다항식은 cubic Hermite basis function임을 알 수 있다.

2차원에서도 동일한 $x$와 $y$ 방향에 동일한 Hermite interpolation을 적용하면 행렬 $\bf P$는 4개의 꼭짓점에서 함숫값과 1차 도함수 및 교차도함수값이 $x$와 $y$의 교환에 대해 대칭적인 꼴을 가져야 한다는 것을 알 수 있다. 1차원 cubic interpolation에서 $\bf X^TM$과 곱해지는 열벡터가 $\begin {bmatrix} p(0)& p(1)& p'(0)& p'(1)\end {bmatrix}$이므로 동일하게 $\bf M^T Y$곱해지는 행벡터도 같은 $x$대신 $y$로 바뀐 형태가 되어야 하므로 행렬 $\bf P$는 

$$\bf P = \begin{bmatrix} f(0,0)& f(0,1) & f_y(0,0) & f_y(0,1) \\  f(1,0) &  f(1,1) & f_y(1,0) & f_y(1,1) \\ f_x(0,0) & f_x(0,1) & f_{xy}(0,0) & f_{xy}(0,1) \\ f_x(1,0) & f_x(1,1) & f_{xy}(1,0) & f_{xy} (1,1)\end {bmatrix}$$처럼 표현되어야 한다. 이 경우 $x$와 $y$를 교환하면 $\bf P$가 $\bf P^T$로 바뀌면서 보간식의 동일한 구조가 유지됨을 알 수 있다.

이제 미분값을 중앙차분으로 근사하면 1차원 cubic interpolation에서 도함수는 $p(0), p(1)$ 뿐만 아니라 $p(-1), p(2)$에서 값도 필요하다. 구체적으로 

$$ \begin{bmatrix} p(0) \\ p(1) \\p'(0) \\p'(1)\end {bmatrix} = \bf D \begin {bmatrix} p(-1) \\ p(0) \\ p(1) \\ p(2) \end {bmatrix}, \qquad D= \begin {bmatrix}   0 & 1 & 0 & 0  \\ 0 & 0 & 1 & 0 \\ -\frac {1}{2} & 0 & \frac {1}{2} & 0 \\ 0 & -\frac {1}{2} & 0 &\frac {1}{2}  \end {bmatrix}$$이므로 이를 2차원에 대해서도 적용하면 

$$\bf P = \bf D Q D^T, \qquad    Q= \begin{bmatrix} f(-1,-1)&f(-1,0)&f(-1,1)&f(-1,2)\\ f(0,-1)&f(0,0)&f(0,1)&f(0,2)\\ f(1,-1)&f(1,0)&f(1,1)&f(1,2)\\ f(2,-1)&f(2,0)&f(2,1)&f(2,2) \end {bmatrix}$$처럼 쓸 수 있고, 최종적으로 보간함수를

$$ f(x, y) = \bf  X^T C Q C^T Y,   \qquad C= MD = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -\frac12 & 0 & \frac12 & 0\\ 1 & -\frac52 & 2 & -\frac12\\ -\frac12 & \frac32 & -\frac32 & \frac12  \end{bmatrix}$$처럼 표현할 수 있다. 여기서

$$\bf w(x) =\bf  X^TC = \begin{bmatrix} -\frac {1}{2} x + x^2 -\frac {1}{2} x^3\\ 1-\frac {5}{2}x^2 +\frac {3}{2} x^3 \\ \frac {1}{2} x + 2x^2 -\frac {3}{2} x^3 \\ -\frac {1}{2} x^2 + \frac {1}{2} x^3\end {bmatrix}^T$$인데, 각 성분은 Catmull-Rom cubic interpolation basis 함수임을 알 수 있다.

따라서 꼭짓점에서의 미분값을 중앙차분으로 근사한 경우 bicubic Hermite interpolation은

$$ f(x, y) = {\bf w(x) Q w(y)^T} =  \sum_{i=0}^3 \sum _{j = 0}^3 w_i(x) Q_{ij}  w_j(y)= \sum_{j=0}^3 w_j(y) \left( \sum_{i=0}^3 w_i(x) Q_{ij}\right)  $$로 쓸 수 있고, 결국 bicubic interpolation은 Catmull-Rom cubic interpolation을 $x$, $y$ 방향으로 차례로 적용하는 separable interpolation과 동일하다. https://kipl.tistory.com/277.

$y=-1, 0, 1, 2\,$일 때 $0 \le x \le 1$ 위치에서 삼차보간된 값이 각각 $j=0,1,2,3\,$에 해당하는

\begin{align}q_j (x) &= \sum_{i = 0} ^3 w_i (x) Q_{ij} \\ &= w_0(x)f(-1,j-1)+w_1(x)f(0, j-1)+w_2(x)f(1, j-1)+w_3(x) f(2, j-1)\end{align}이고, 이 4개의 $y\,$방향 값을 이용해서 $0\le y\le 1$에서 삼차보간한 결과가

$$ \sum_{j=0}^3 w_j (y) q_j(x) = f(x, y)$$이다. 따라서 16개의 $a_{ij}$를 구해서 구현한 bicubic interpolation 코드는 (https://kipl.tistory.com/55) 아래처럼 좀 더 간결한 형태로 정리할 수 있다.

// 등간격으로 샘플링된 값 p0, p1, p2, p3를 이용해서 중앙 구간에서 값을 구하는
// Catmull-Rom interpolation
double CubicInterpolate(double p0, double p1, double p2, double p3, double t) {
    double a0 = p1;
    double a1 = 0.5 * (p2 - p0);
    double a2 = p0 - 2.5*p1 + 2.0*p2 - 0.5*p3;
    double a3 = 0.5 * (-p0 + 3.0*p1 - 3.0*p2 + p3);
    return ((a3*t + a2)*t + a1)*t + a0;
}
// p[i][j] = f(i-1,j-1)
double BicubicInterpolate(const double p[4][4], double x, double y) {
    double q[4];
    for (int j = 0; j < 4; ++j)
        q[j] = CubicInterpolate(p[0][j], p[1][j], p[2][j], p[3][j], x);
    return CubicInterpolate(q[0], q[1], q[2], q[3], y);
}