Geometry & Recognition

히스토그램처럼 균일한 간격으로 주어진 데이터에서 피크 값의 위치를 찾기 위해서는 먼저 노이즈 제거를 위해 몇 번의 smoothing 과정을 적용해야 한다 (구체적인 방법은 필요에 따라 다양하다). smoothing 과정을 거친 히스토그램에서 피크 위치를 찾는 것은 쉬운 작업이다. 그런데 경우에 따라 그 위치를 sub-order까지 구해야 필요가 생긴다. 예를 들면, 실수 값 데이터 시퀀스 대한 히스토그램을 만들려면 먼저 정수로 양자화시켜야 가능하다. 이 양자화된 정보를 담고 있는 히스토그램에서 피크의 위치를 찾으려 할 때 양자화 과정에서 잃어버린 정밀도를 복원하려면 interpolation을 써야 한다. 간단하게 parabolic 근사로 피크의 위치를 찾는 경우를 생각하자. 이 방법은 수학적으로 단순하지만 영상처리 알고리즘에서 많이 이용이 되고 있다. 데이터 시퀀스가 균일한 간격에서 주어졌으므로 계산은 $-1,0,1$의 세 군데 위치에서 중앙 근방에 피크가 나타나는지를 고려하면 된다. 세 점에서 히스토그램 값이 각각 $h_m,h_0,h_p$일 때 이들을 지나는 이차 곡선의 피크 위치를 찾자. 주어진 이차 곡선을

$$y=a(x-c{)}^2+b$$

꼴로 쓰면 $c$가 0에서 벗어난 정도를 나타낸다.

$$(-1,h_m)\to h_m=a(-1-c{)}^2+b$$

$$(0,h_0)\to h_0=a{c}^2+b$$

$$(+1,h_p)\to h_m=a(+1-c{)}^2+b$$

이므로

$$h_m-h_p=4ac,$$

$$h_m+h_p-2h_0=2a,$$

$$\therefore c=\frac{h_m-h_p}{2(h_m+h_p-2h_0)}$$

아래 코드는 피크의 위치가 중앙점에서 벗어난 정도를 준다.

bool parabolicInterpolate(double hm, double h0, double hp, double *c) {
      double a = hm + hp - 2 * h0;
      if (a >= 0) return false; // not a parabola(a==0), not convex up(a>0);
      *c = 0.5 * (hm - hp) / a;
      if (*c < -0.5 || *c > 0.5) return false; //  too far;
      return true;
}

3차의 Hermite spline을 쓰는 interpolation에서 입력 데이터 점에서의 기울기 값을 수정하면 입력 데이터의 monotonic 한 성질을 그대로 유지하는 보간을 할 수 있다. 주어진 데이터가 $\{(x_1,y_1),(x_2,y_2),\dots ,(x_n,y_n)\}$이고, $x/y$-좌표가 증가하는 순으로 정렬되었다고 하자. 구간 $[x_k,x_{k+1}]$의 시작과 끝에서 값 ($p_0=y_k$, $p_1=y_{k+1}$)과 기울기 ($m_0,m_1$) 가 주어지면 cubic Hermite spline 보간 함수는 다음과 같이 매개변수 $t\in [0,1]$의 함수로 표현된다: 

$$f(t)=p_0{h}_{00}(t)+m_1{h}_{10}(t)+p_1{h}_{01}(t)+m_1{h}_{11}(t),t=\frac{x-x_k}{x_{k+1}-x_k}$$

$$f(0)=p_0,f^{\prime }(0)=m_0,f(1)=p_1,f^{\prime }(1)=m_1$$

Hermite 기저 함수의 정의는 아래와 같다:

$$h_{00}(t)=2t^3-3t^2+1$$

$$h_{10}(t)=t^3-2t^2+t$$

$$h_{01}(t)=-2t^3+3t^2$$

$$h_{11}(t)=t^3-t^2$$

이를 보이기 위해 3차 보간함수를 

$$f(t)=a_0+a_{1}t+a_2{t}^2+a_3{t}^3,0\le t\le 1$$으로 놓고, constraints를 부여하면

$$\begin{array}{rl}{p}_0& =f(0)=a_0\\ p_1& =f(1)=a_0+a_1+a_2+a_3\\ m_0& =f^{\prime }(0)=a_1\\ m_1& =f^{\prime }(1)=a_1+2a_2+3a_3\end{array}$$을 얻을 수 있다. 계수 $a_0,a_1,a_2,a_3$에 대해 풀면

$$f(t)=p_0+m_{0}t+(3p_1-3p_0-m_1-2m_0)t^2+(-2p_1+2p_0+m_0+m_1)t^3$$

$$=p_0(2t^3-3t^2+1)+p_1(-2t^3+3t^2)+m_0(t^3-2t^2+t)+m_1(t^3-t^2)$$ 처럼 표현됨을 알 수 있다.

보간 spline이 단조성(monotonicity)을 가지도록 하려면 각 입력점에서 기울기 값($m_k$)을 조절해주어야 하는데, 기울기 $m_k$를 어떻게 설정해야 인접 구간의 보간 함수와 smooth 하게 연결될 수 있을까? 주어진 입력점에서 기울기는 보통 왼쪽과 오른쪽 인접 입력점에서 각각 평균 변화율(${\delta }_{k-1},{\delta }_k$)을 구한 후 이들의 평균을 사용하면 된다. (양 끝점에서는 오른쪽이나 왼쪽 구간의 평균 변화율로 대체)

$$m_k=\frac{1}{2}({\delta }_{k-1}+{\delta }_k)=\frac{1}{2}(\frac{y_k-y_{k-1}}{x_k-x_{k-1}}+\frac{y_{k+1}-y_k}{x_{k+1}+x_k})$$

이것만으로는 단조성을 보증을 할 수 없으므로 다음 조건을 추가로 부여해야 됨이 관련 논문에서 증명이 되었다. 

$$1.\text{For}{\delta }_k=0⟹m_k=m_{k+1}=0$$

$$2.\text{For}ϵ=\sqrt{{\left(\frac{m_k}{{\delta }_k}\right)}^2+{\left(\frac{m_{k+1}}{{\delta }_k}\right)}^2}>3$$

$$⟹m_k\leftarrow 3\frac{m_k}{ϵ},m_{k+1}\leftarrow 3\frac{m_{k+1}}{ϵ}$$

아래의 코드는 위키피디아에 제시된 알고리즘을 구현한 것이다. (입력 데이터가 한 방향으로 monotonic 할 때만 동작함)

참고: http://en.wikipedia.org/wiki/Monotone_cubic_interpolation 

std::vector<double> estimateTangents(const std::vector<CfPt>& P) {
    // slopes;
    std::vector<double> delta(P.size());    
    for (int k = 0; k < delta.size()-1; k++)
        delta[k] = double(P[k+1].y - P[k].y) / (P[k+1].x - P[k].x) ;
        
    // average tangent at data points;
    std::vector<double> m(P.size());
    m.front() = delta.front();
    m.back() = delta[m.size()-2];
    for (int k = 1; k < m.size()-1; k++)
        m[k] = (delta[k-1] + delta[k]) / 2;  //error corrected.
        
    for (int k = 0; k < m.size()-1; k++) {
        if (delta[k] == 0) m[k] = m[k+1] = 0;
        else {
            double alpha = m[k] / delta[k];
            double beta  = m[k+1] / delta[k];
            double epsilon = hypot(alpha, beta);
            if (epsilon > 3) {
                double tau = 3 * delta[k] / epsilon ;
                m[k]   = tau * alpha;
                m[k+1] = tau * beta;
            }
        }
    }
    return m;
};

std::vector<CfPt> monotoneCubicSplineInterp(const std::vector<CfPt>& points) {
    if (points.size() <= 2) return points;
    // tangent table;
    std::vector<double> m = estimateTangents(points);
    // hermite spline;
#define STEPS 12
    std::vector<CfPt> curve;
    curve.push_back(points.front()); //add first point;
    for (int i = 0; i < points.size()-1; i++) {
        double h = points[i+1].x - points[i].x;
        for (int k = 1; k <= STEPS; k++) {
            double t = double(k) / STEPS;         // 0 < t <= 1;
            double tt = t*t, ttt = tt*t;
            double h00 =  2 * ttt - 3 * tt     + 1;
            double h10 =      ttt - 2 * tt + t;
            double h01 = -2 * ttt + 3 * tt;
            double h11 =      ttt     - tt;
            double yf = h00 * points[i].y + h10 * h * m[i] 
                        + h01 * points[i+1].y + h11 * h * m[i+1];
            curve.push_back(CfPt(points[i].x + t * h, yf));
        }
    }
    return curve;
};

Natural cubic spline과의 비교:

https://kipl.tistory.com/569

이미지 처리에서 픽셀 좌표는 간격이 1인 2차원 그리드의 교차점으로 볼 수 있다. 이미지를 확대하거나 축소할 때, 픽셀과 픽셀 사이 구간에서 값이 필요한 경우가 생긴다. 간단하게는 가장 가까운 주변 픽셀의 값을 그대로 가져다 쓰거나, 또는 주변의 4 픽셀 값을 선형보간해서 사용할 수 있다. 2차원 그리드에서 선형보간은 bilinear interpolation이다. 이 보간법은 속도는 빠르지만 픽셀 값이 인접 그리드에서의 값으로 부드럽게 이어지지 않는 단점이 있다. 인접 그리드 경계에서 픽셀 값이 부드럽게 이어지기 위해서는 적어도 1차 미분이 연속인 보간법을 사용해야 하는데, 이러한 조건을 만족시키는 가장 낮은 찾수의 다항식 보간법이 bicubic interpolation이다.

4점 $(0,0),(0,1),(1,0),(1,1)$을 꼭짓점으로 하는 정사각형 내의 임의 지점 $D=\{(x,y)|0\le x\le 1,0\le y\le 1\}$ 에서 픽셀 값을 주는 보간곡면 $f(x,y)$을 주변의 16개 점 $\{(i,j)|-1\le i\le 2,-1\le j\le 2\}$에서 픽셀 값을 사용하는 bicubic interpolation을 이용해서 추정할 수 있다. 곡면 $f(x,y)$은 $x$와 $y$의 3차 함수로 다음과 같이 쓸 수 있다.

$$\begin{array}{r}f(x,y)=\sum _{i=0}^3\sum _{j=0}^3{a}_{ij}{x}^i{y}^j.\end{array}$$

문제는 16개 계수 $\{a_{ij}\}$를 어떻게 찾을 것인가? 이를 위해 16개의 조건이 필요한데 주변의 16개의 픽셀 값을 이용해서 만들 수 있다. 또한 픽셀 값이 인접 그리드 영역으로 smooth 하게 연결되기 위해서는 $f(x,y)$의 미분도 고려해야 한다. 

1. 4 꼭짓점에서 값 $f(x,y)$: $$\begin{array}{rl}f(0,0)& =a_{00};\\ f(1,0)& =a_{00}+a_{10}+a_{20}+a_{30};\\ f(0,1)& =a_{00}+a_{01}+a_{02}+a_{03};\\ f(1,1)& =a_{00}+a_{10}+a_{20}+a_{30}+a_{01}+a_{11}+a_{21}+a_{31}\\ & +a_{02}+a_{12}+a_{22}+a_{32}+a_{03}+a_{13}+a_{23}+a_{33};\end{array}$$
2. 4 꼭짓점에서 미분계수, $f_x$, $f_y$: $$\begin{array}{rl}{f}_x(0,0)& =a_{10};\\ f_x(1,0)& =a_{10}+2a_{20}+3a_{30};\\ f_x(0,1)& =a_{10}+a_{11}+a_{12}+a_{13};\\ f_x(1,1)& =a_{10}+2a_{20}+3a_{30}+a_{11}+2a_{21}+3a_{31}\\ & +a_{12}+2a_{22}+3a_{32}+a_{13}+2a_{23}+3a_{33};\\ f_y(0,0)& =a_{01}\\ f_y(1,0)& =a_{01}+a_{11}+a_{21}+a_{31};\\ f_y(0,1)& =a_{01}+2a_{02}+3a_{03};\\ f_y(1,1)& =a_{01}+a_{11}+a_{21}+a_{31}+2a_{02}+2a_{12}+2a_{22}+2a_{32}\\ & +3a_{03}+3a_{13}+3a_{23}+3a_{33};\end{array}$$ 
3. 4 꼭짓점에서 교차 미분계수, $f_{ij}$: $$\begin{array}{rl}{f}_{xy}(0,0)& =a_{11};\\ f_{xy}(1,0)& =a_{11}+2a_{21}+3a_{31};\\ f_{xy}(0,1)& =a_{11}+2a_{12}+3a_{13};\\ f_{xy}(1,1)& =a_{11}+2a_{21}+3a_{31}+2a_{12}+4a_{22}\\ & +6a_{32}+3a_{13}+6a_{23}+9a_{33};\end{array}$$

꼭짓점에서 $f(x,y)$의 미분계수는 꼭짓점과 인접 교차점 사이의 평균 변화율로 표현할 수 있다. 주변의 픽셀이 $p[i+1,j+1]=f(i,j)$로 주어진 경우;

$$\begin{array}{rl}{f}_{00}& =p[1,1],f_{10}=p[2,1],...\\ f_{x00}& =f_x(0,0)=\frac{1}{2}(f(1,0)-f(-1,0))=\frac{1}{2}(p[2,1]-p[0,1]),...\\ f_{xy00}& =f_{xy}(0,0)=\frac{1}{4}(f(1,1)-f(1,-1)-f(-1,1)+f(-1,-1))\\ & =\frac{1}{4}(p[2,2]-p[2,0]-p[0,2]+p[0,0]),...\end{array}$$

와 같다.

위의 식들은 16개의 미지수 {$a_{ij}$}를 가지는 16개의 연립 방정식이므로 이를 풀어서 계수를 구할 수 있다. 계수를 벡터 $\mathbf{v}$로 

$$\mathbf{v}=[a_{00},a_{10},a_{20},a_{30},a_{01},a_{11},a_{21},a_{31},a_{02},a_{12},a_{22},a_{32},a_{03},a_{13},a_{23},a_{33}{]}^t$$

로 나타내고, 각 꼭짓점에서의 픽셀 값, 미분 값을 $\mathbf{f}$ 벡터로 나타내면

$$\mathbf{f}=[f_{00},f_{10},f_{01},f_{11},f_{x00},f_{x10},f_{x01},f_{x11},f_{y00},f_{y10},f_{y01},f_{y11},f_{xy00},f_{xy10},f_{xy01},f_{xy11}{]}^t.$$

16개의 연립 방정식은 행렬식

$$\mathbf{f}=\mathbf{A}\cdot \mathbf{v}(해:\mathbf{v}={\mathbf{A}}^{-1}\cdot \mathbf{f})$$

으로 쓸 수 있다. 여기서,

$$\mathbf{A}=\left(\begin{array}{cccccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0& 0\\ 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 2& 3& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0\\ 0& 1& 2& 3& 0& 1& 2& 3& 0& 1& 2& 3& 0& 1& 2& 3\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 2& 0& 0& 0& 3& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1& 2& 2& 2& 2& 3& 3& 3& 3\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 2& 3& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 2& 0& 0& 0& 3& 0& 0\\ 0& 0& 0& 0& 0& 1& 2& 3& 0& 2& 4& 6& 0& 3& 6& 9\end{array}\right)$$

그리고 역행렬은

$${\mathbf{A}}^{-1}=\left(\begin{array}{cccccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -3& 3& 0& 0& -2& -1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 2& -2& 0& 0& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -3& 3& 0& 0& -2& -1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 2& -2& 0& 0& 1& 1& 0& 0\\ -3& 0& 3& 0& 0& 0& 0& 0& -2& 0& -1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& -3& 0& 3& 0& 0& 0& 0& 0& -2& 0& -1& 0\\ 9& -9& -9& 9& 6& 3& -6& -3& 6& -6& 3& -3& 4& 2& 2& 1\\ -6& 6& 6& -6& -3& -3& 3& 3& -4& 4& -2& 2& -2& -2& -1& -1\\ 2& 0& -2& 0& 0& 0& 0& 0& 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 2& 0& -2& 0& 0& 0& 0& 0& 1& 0& 1& 0\\ -6& 6& 6& -6& -4& -2& 4& 2& -3& 3& -3& 3& -2& -1& -2& -1\\ 4& -4& -4& 4& 2& 2& -2& -2& 2& -2& 2& -2& 1& 1& 1& 1\end{array}\right)$$

이 역행렬은 mathematica나 maple 등의 프로그램을 이용하면 쉽게 구할 수 있다. 이렇게 구한 $f(x,y)$는 $[0,1]\times [0,1]$ 영역에서 연속인 smooth 한 곡면을 형성한다. 또한 인접하는 그리드 영역의 곡면과 경계에서 smooth 하게 연결된다.

/* 픽셀의 위치가 grid의 교차점이 되도록 표현하였다*/

/* [ix, ix+1][iy, iy+1] 영역에서 보간을 하는 경우에, 주변 픽셀의 정보는 
** p[a][b]에 들어가고, p[a][b]= Image (ix - 1 + a, iy - 1 + b)
** 0 <= x < 1, 0 <= y < 1 
*/
double BicubicInterpolate(double p[4][4], double x, double y) { 
    double a00 = p[1][1];
    double a01 = -.5*p[1][0] + .5*p[1][2];
    double a02 = p[1][0] - 2.5*p[1][1] + 2*p[1][2] - .5*p[1][3];
    double a03 = -.5*p[1][0] + 1.5*p[1][1] - 1.5*p[1][2] + .5*p[1][3];
    double a10 = -.5*p[0][1] + .5*p[2][1];
    double a11 = .25*p[0][0] - .25*p[0][2] - .25*p[2][0] + .25*p[2][2];
    double a12 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + 
                  .5*p[2][0] - 1.25*p[2][1] + p[2][2] - .25*p[2][3];
    double a13 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - 
                 .25*p[2][0] + .75*p[2][1] - .75*p[2][2] + .25*p[2][3];
    double a20 = p[0][1] - 2.5*p[1][1] + 2*p[2][1] - .5*p[3][1];
    double a21 = -.5*p[0][0] + .5*p[0][2] + 1.25*p[1][0] - 1.25*p[1][2] - 
                 p[2][0] + p[2][2] + .25*p[3][0] - .25*p[3][2];
    double a22 = p[0][0] - 2.5*p[0][1] + 2*p[0][2] - .5*p[0][3] - 2.5*p[1][0] + 
                 6.25*p[1][1] - 5*p[1][2] + 1.25*p[1][3] + 2*p[2][0] - 5*p[2][1] + 
                 4*p[2][2] - p[2][3] - .5*p[3][0] + 1.25*p[3][1] - p[3][2] + .25*p[3][3];
    double a23 = -.5*p[0][0] + 1.5*p[0][1] - 1.5*p[0][2] + .5*p[0][3] + 1.25*p[1][0] - 
                 3.75*p[1][1] + 3.75*p[1][2] - 1.25*p[1][3] - p[2][0] + 3*p[2][1] - 
                 3*p[2][2] + p[2][3] + .25*p[3][0] - .75*p[3][1] + .75*p[3][2] - .25*p[3][3];
    double a30 = -.5*p[0][1] + 1.5*p[1][1] - 1.5*p[2][1] + .5*p[3][1];
    double a31 = .25*p[0][0] - .25*p[0][2] - .75*p[1][0] + .75*p[1][2] + 
                 .75*p[2][0] - .75*p[2][2] - .25*p[3][0] + .25*p[3][2];
    double a32 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + 1.5*p[1][0] - 
                 3.75*p[1][1] + 3*p[1][2] - .75*p[1][3] - 1.5*p[2][0] + 3.75*p[2][1] - 
                 3*p[2][2] + .75*p[2][3] + .5*p[3][0] - 1.25*p[3][1] + p[3][2] - .25*p[3][3];
    double a33 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .75*p[1][0] + 
                 2.25*p[1][1] - 2.25*p[1][2] + .75*p[1][3] + .75*p[2][0] - 2.25*p[2][1] + 
                 2.25*p[2][2] - .75*p[2][3] - .25*p[3][0] + .75*p[3][1] - 
                 .75*p[3][2] + .25*p[3][3];
    double x2 = x * x;
    double x3 = x2 * x;
    return a00 + (a01 + (a02  + a03 * y) * y) * y +
          (a10 + (a11 + (a12  + a13 * y) * y) * y) * x  +
          (a20 + (a21 + (a22  + a23 * y) * y) * y) * x2 +
          (a30 + (a31 + (a32  + a33 * y) * y) * y) * x3;
};

* 원본 4x4 이미지(RGB):  nearest-neighbor interploation으로 256x256 크기로 만듦;

* bicubic interpolation 결과 (256x256):

* bilinear interpolation 결과 (256x256):

코드 구현 일부:  bicubic interpolation은 주변의 16개 픽셀 정보가 필요한데, 가장자리 픽셀인 경우는 이를 충족시킬 수 없으므로 이를 해결하기 위해 영역 밖은 가장자리 픽셀이 반복된 것으로 처리하는 것이 가장 쉽다. 그리고 소스 영상이 큰 경우에는 문제가 되지 않지만 예제처럼 작은 소스 영상의 경우는 가장자리 부근에서 단순 채움을 사용하면 왜곡이 발생한다. 이를 해소하기 위해서는  소스 픽셀 위치가 픽셀의 중간을 나타낸 것으로 처리하면 된다. 
더보기

     for (int y = 0; y < dstHeight; y++) {
        BYTE *pdst = (BYTE *)dst.GetLinePtr(y);
        double src_fy = y * scaley - 0.5;
        int src_iy = int(floor(src_fy));         
        double ry = src_fy - src_iy;              
        for (int x = 0; x < dstWidth; x++) {
            double src_fx = x * scalex - 0.5;
            int src_ix = int(floor(src_fx));
            double rx = src_fx - src_ix;
            // Fill interpolation matrices for R,G,B;
            for (int by = 0; by < 4; by++) {
                int yy = by + src_iy - 1;
                CLAMP(yy, 0, srcHeight - 1);
                BYTE *psrc = (BYTE *)src.GetLinePtr(yy);
                for (int ax = 0; ax < 4; ax++) {
                    int xx = ax + src_ix - 1;
                    CLAMP(xx, 0, srcWidth - 1);
                    B[ax][by] = psrc[3 * xx + 0];
                    G[ax][by] = psrc[3 * xx + 1];
                    R[ax][by] = psrc[3 * xx + 2];
                }
            }
            int bb = (int)BicubicInterpolate(B, rx, ry);
            int gg = (int)BicubicInterpolate(G, rx, ry);
            int rr = (int)BicubicInterpolate(R, rx, ry);
            *pdst++ = bb < 0 ? 0: bb > 255 ? 255: bb;
            *pdst++ = gg < 0 ? 0: gg > 255 ? 255: gg;
            *pdst++ = rr < 0 ? 0: rr > 255 ? 255: rr;
        }
    };