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개의 교차편도함수값으로 결정된다.

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

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

$$ 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$$ 이 경우 행렬 $M$은 $x$나 $y$ 위치에 의존하지 않는 상수행렬이어야 하고, 행렬 \(P\)는 네 꼭짓점에서의 함수값과 미분값이 들어가도록 형태를 가질 것이다.

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

$$ p(t) = \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}$$구체적으로 계산을 하면

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

$$ 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$를 교환하면 $P$가 $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} = 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차원에 대해서도 적용하면 

$$P = 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) =  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}$$처럼 표현할 수 있다. 여기서

$$w(x) = 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) = w(x) Q w(y)^T =  \sum_{i=0}^3 \sum _{j = 0}^3 w_i(x) Q_{ij}  w_i(y)= \sum_{i=0}^3 w_i(x)\left( \sum_{j=0}^3 Q_{ij} w_j(y)\right)$$로 쓸 수 있고, 결국 bicubic interpolation은 Catmull-Rom cubic interpolation을 $x$, $y$ 방향으로 차례로 적용하는 separable interpolation과 동일하다. https://kipl.tistory.com/277. 이 결과를 이용하면 https://kipl.tistory.com/55 코드는 아래처럼 좀 더 간결한 형태로 정리할 수 있다.

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