Geometry & Recognition

Contrast Limited Adaptive Histogram Equalization (CLAHE). CLAHE는 영상의 평탄화 (histogram equalization) 과정에서 contrast가 과도하게 증폭이 되는 것을 방지하도록 설계된 adaptive algorithm의 일종이다. CLAHE algorithm은 전체 영상 영역을 바로 적용하지 않고 tile이라 불리는 작은 영역으로 영상을 분할하여 tile 별로 적용한다. 그리고 인접하는 tile 간에 bilinear 보간을 사용함으로써 tile 경계에서 급격한 변화가 생기는 것을 방지한다. gray 영상에 CLAHE 알고리즘을 적용하면 영상의 contrast을 개선해 준다. 그리고 컬러 영상에도 적용을 할 수 있는데, 보통 luminance 채널에만 적용한 결과(HSV로 변환 후 V에 대해서만 적용)가 RGB 각 채널별로 적용하는 것보다 좋은 결과를 준다.

아래는 RGB 영상의 각 채널에 CLAHE 알고리즘을 적용한 결과이다.

rgn: 히스토그램을 구하는 이미지 상의 영역;

tile: adaptive HE가 적용되는 영역

첫 행/열 tile의 폭/높이 = rgn 폭/높이의 절반

마지막 tile 크기 = 이미지의 나머지 폭/높이 

cliplimit = 각 영역의 히스토그램 슬롯당 들어가는 최대 픽셀 수

/* "Graphics Gems IV" 의 code를 수정함;
** 2023-01-13 updated: 이미지 사이즈가 rgn 사이즈의 정수배가 되지 않을 때 문제 해결;*/
int CLAHE ( BYTE* image, int width, int height,
            BYTE gMin/*0*/, BYTE gMax/*255*/, int rgn_nx/*16*/, int rgn_ny/*16*/,
            int bins/*256*/, double fcliplimit )
{
    std::vector<int> map_array(rgn_nx * rgn_ny * bins);
    int rgn_xsz		= width / rgn_nx;
    int rgn_ysz		= height / rgn_ny; 
    int rgn_area	= rgn_xsz * rgn_ysz;
    int clipLimit	= int( fcliplimit * ( rgn_xsz * rgn_ysz ) / bins );
    clipLimit = ( clipLimit < 1 ) ? 1 : clipLimit;
    //  
    /* calculate greylevel mappings for each contextual region */
    for (int iy = 0; iy < rgn_ny; iy++) {
        for (int ix = 0; ix < rgn_nx; ix++) {
            int *hist = &map_array[bins * ( iy * rgn_nx + ix )];
            int start = (iy * rgn_ysz) * width + ix * rgn_xsz;
            RgnHistogram (&image[start], width, rgn_xsz, rgn_ysz, hist, bins);
            ClipHistogram ( hist, bins, clipLimit );
            //convert hist to cuimulative histogram normalized to [gMin=0, gMax=255]
            MapHistogram ( hist, gMin, gMax, bins, rgn_area );
        }
    }

    /* bilinearInterp greylevel mappings to get CLAHE image */
    /* 첫열/행에 해당하는 타일은 rgn의 절반 크기만 처리하므로 ix(iy)는 0~rgn_nx(rgn_ny)까지임*/
    int szx, szy, ixl, ixr, iyt, iyb;
    BYTE *first_tile = &image[0]; // 같은 iy를 갖는 첫 타일의 주소;
    for (int iy = 0; iy <= rgn_ny; iy++ ) {
        if ( iy == 0 ) {				
            // 첫 타일은 절반만;
            szy = rgn_ysz / 2; iyt = iyb = 0;
        } else if ( iy == rgn_ny ) {
            // 마지막 타일은 나머지 전부;
            szy = height - ((rgn_ny-1)*rgn_ysz + rgn_ysz/2);
            iyt = rgn_ny - 1; iyb = iyt;
        } else {
            szy = rgn_ysz; iyt = iy - 1; iyb = iyt + 1;
        }
        BYTE *ptile = &first_tile[0]; //각 타일의 시작 주소;
        for (int ix = 0; ix <= rgn_nx; ix++ ) {
            if ( ix == 0 ) {	
                //첫 타일은 절반만;
                szx = rgn_xsz / 2; ixl = ixr = 0;
            } else if ( ix == rgn_nx ) {		
                //마지막 타일은 나머지 전부;
                szx = width - ((rgn_nx-1)*rgn_xsz + rgn_xsz/2);
                ixl = rgn_nx - 1; ixr = ixl;
            } else {
                szx = rgn_xsz; ixl = ix - 1; ixr = ixl + 1;
            }
            // cumulative histogram data;
            int *LU = &map_array[bins * ( iyt * rgn_nx + ixl )];
            int *RU = &map_array[bins * ( iyt * rgn_nx + ixr )];
            int *LB = &map_array[bins * ( iyb * rgn_nx + ixl )];
            int *RB = &map_array[bins * ( iyb * rgn_nx + ixr )];
            BilinearInterp (ptile, width, LU, RU, LB, RB, szx, szy );
            ptile += szx;			  
        }
        first_tile += szy * width;
    }
    return 0;						 
}

void ClipHistogram ( int* hist, int gLevels, int clipLimit ) {
    int excess = 0;
    for (int i = 0; i < gLevels; i++ ) { /* calculate total number of excess pixels */
        int diff =  hist[i] - clipLimit;
        if ( diff > 0 ) excess += diff;	 /* excess in current bin */
    };

    /* clip histogram and redistribute excess pixels in each bin */
    int incr = excess / gLevels;		/* average binincrement */
    int upper =  clipLimit - incr;		/* bins larger than upper set to cliplimit */
    for (int i = 0; i < gLevels; i++ ) {
        if ( hist[i] > clipLimit ) hist[i] = clipLimit; /* clip bin */
        else {
            if ( hist[i] > upper ) {			/* high bin count */
                excess -= hist[i] - upper;
                hist[i] = clipLimit;
            } else {							/* low bin count */
                excess -= incr;
                hist[i] += incr;
            }
        }
    }

    while ( excess ) { /* redistribute remaining excess  */
        int start = 0;
        while ( excess && start < gLevels) {
            int step = gLevels / excess;
            if ( step < 1 ) step = 1;		 /* step size at least 1 */
            for (int i = start; i < gLevels && excess; i += step) 
                if (hist[i] < clipLimit) { 
                    hist[i]++; excess--;
                }
            start++;
        }
    }
}
void BilinearInterp ( BYTE * image, int stride, 
                    int* LU, int* RU, int* LB,  int* RB,
                    int rgn_xsz, int rgn_ysz ) {
    int area = rgn_xsz * rgn_ysz; /* Normalization factor */
    for (int y = 0, yrev = rgn_ysz; y < rgn_ysz; y++, yrev--) {
        BYTE *line = &image[y * stride];
        for (int x = 0, xrev = rgn_xsz; x < rgn_xsz; x++, xrev-- ) {
            int v = line[x]; /* get histogram bin value */
            v = (yrev * ( xrev * LU[v] + x * RU[v] ) + y * ( xrev * LB[v] + x * RB[v] ) ) / area;
            line[x] = v > 255 ? 0xFF: v < 0 ? 0x00: v;
        }
        //image += stride - rgn_xsz;  // goto the beginning of rgn's next-line;
    }
}
void RgnHistogram ( BYTE* image, int stride, int rgn_xsz, int rgn_ysz, 
                        int* hist, int gLevels) {
    // image is pointing the beginnign of a given region;
    for (int i = 0; i < gLevels; i++ ) hist[i] = 0;
    for (int y = 0; y < rgn_ysz; y++ ) {
        BYTE *line = &image[y * stride];
        for (int x = 0; x < rgn_xsz; x++)
            hist[line[x]]++;
    }
}
// transform hist[] to a normalized cumulative histogram;
void MapHistogram ( int* hist, BYTE gMin, BYTE gMax, int gLevels, int rgn_area ){
    double scale = double( gMax - gMin ) / rgn_area;
    int iMin = int(gMin);
    for (int i = 0, sum = 0; i < gLevels; i++ ) {
        sum += hist[i];
        int a = int(iMin + sum * scale );
        hist[i] = a > gMax ? gMax: a < gMin ? gMin: a;
    }
}

n이 양수일 때 n 과 n-1은 짝홀/홀짝의 짝이므로 2의 파워가 아니면 이진수로 표현할 때 한 곳에서 비트 차이가 난다. 따라서 AND연산을 하면 0이 아니다.  그런데 n이 2의 파워이면 최상위 비트만 1이고, n-1은 n의 최상위 비트를 제외한 하위비트 모두가 1로 구성되므로 AND 연산을 하면 0이다.

예:   n=6 = 110,  6-1=101 -->  6 & 5 = 100;

       n=8 = 1000, 8-1=0111 --> 8 & 7=0000

// n이 양수이고, n-1과 겹치는 비트가 없으면 된다;
bool IsPowerOf2(int n) {
    return n > 0 && (n & (n - 1)) == 0;
}

n>1 일 떄 2로 나누기를 계속할 때 1을 제외한 홀수가 나오면 안됨; 

bool IsPowerOf2(int n) {
    if (n == 1) return true;
    if (n == 0 || n & 1) return false;
    return IsPowerOf2(n >> 1);
} // n == 0;

x보다 크기 않은 2^i을 구한 후 x와 같은가 체크;

bool IsPowerof2(int x){
    int i = 0;
    while ((1 << i) < x) i++;
    if (x == (1 << i)) return true;
    return false;
}

2의 파워이면 log_2(n) = 자연수 이므로 ceil과 floor가 같다.

bool IsPowerOf2(int n) {
   if (n == 0) return false;
   return (ceil(log2(n)) == floor(log2(n)));
}

주어진 자연수보다 작지 않은 최소의 2^n;

int NextPowerOf2(int n) { //32-bit;
    n--;
    n |= n >> 1;
    n |= n >> 2;
    n |= n >> 4;
    n |= n >> 8;
    n |= n >> 16;
    n++;
    return n;    
} // NextPowerOf2(5) -> 8; NextPowerOf2(8) -> 8;

// labeling using depth-first search; non-recursive version;
int GetConnectedComponents(BYTE *image, int w, int h, int *table) {
    int label = 0;    // starting label = 1;
    std::vector<int> stack(w * h);
    // initialize the table;
    for (int k = npixs; k-- > 0;) 
        table[k] = image[k] ? -1: 0;  // Foreground = -1; Background = 0;

    for (int pos = w * h; pos-- > 0;) {
        if (table[pos] == -1) { // Foreground;
            ++label;            // assign next label;
            int top = -1;       // stack initialization;
            stack[++top] = pos;
            while (top >= 0) {
                int adj = stack[top--];
                int xx = adj % w;
                int yy = adj / w;
                if (table[adj] == -1) {// Foreground;
                    table[adj] = label;
                    // check 4-way connectivity;
                    if (xx + 1 < w) stack[++top] = adj + 1; //RIGHT;
                    if (yy + 1 < h) stack[++top] = adj + w; //BOTTOM;
                    if (yy > 0)     stack[++top] = adj - w; //TOP;
                    if (xx > 0)     stack[++top] = adj - 1; //LEFT;
                }
            }
        }
    }
    return label;    // total # of CCs;
};

이미지 처리 과정에서 미분에 해당하는 그래디언트 필드(gradient field: $g_x$, $g_y$ )를 이용하면 이미지 상의 특징인 corner, edge, ridge 등의 정보를 쉽게 얻을 수 있다. 이미지의 한 지점이 이러한 특징을 가지는 특징점이 되기 위해서는 그래디언트 필드의 값이 그 점 주변에서 (3x3나 5x5정도 크기의 window) 일정한 패턴을 유지해야 한다. 이 패턴을 찾기 위해서 그래디언트 필드에 PCA를 적용해보자. 수평과 수직방향의 그래디언트 field인 $g_x$와 $g_y$ 사이의 covariance 행렬은 다음 식으로 정의된다:

$$ \Sigma = \left [ \begin {array}{cc} < g_x^2 > & < g_x g_y > \\    <g_x g_y> & <g_y^2 > \end {array}\right] =\left [ \begin {array}{cc} s_{xx} & s_{xy} \\ s_{xy} & s_{yy}\end {array}\right];$$

$<...> = \int_{W}(...) dxdy$는 픽셀 윈도에 대한 적분을 의미한다. $\Sigma$의 eigenvalue는 항상 음이 아닌 값을 갖게 되는데 (matrix 자체가 positive semi-definitive), 두 eigenvalue이 $λ_1$, $λ_2$면

$$λ_1 + λ_2 = s_{xx} + s_{yy} \ge 0, \quad \quad    λ_1  λ_2 = s_{xx} s_{yy} - s_{xy}^2 \ge0 $$

을 만족한다 (완전히 상수 이미지를 배제하면 0인 경우는 없다). eigenvalue $λ_1$, $λ_2$는 principal axis 방향으로 그래디언트 필드의 변동(분산)의 크기를 의미한다. edge나 ridge의 경우는 그 점 주변에서 잘 정의된 방향성을 가져야 하고, corner의 경우는 방향성이 없어야 한다. edge나 ridge처럼 일방향성의 그래디언트 특성을 갖거나 corner처럼 방향성이 없는 특성을 서로 구별할 수 있는 measure가 필요한데, $λ_1$과 $λ_2$를 이용하면 차원이 없는 measure을 만들 수 있다. 가장 간단한 차원이 없는 측도(dimensionless measure)는  eigenvalue의 기하평균과 산술평균의 비를 비교하는 것이다.

$$ Q = \frac { {λ_{1} λ_{2}} }{ \left( \frac {λ_{1}+λ_{2}}{2} \right)^2} = 4\frac { s_{xx} s_{yy} - s_{xy}^2}{(s_{xx} + s_{yy})^2};$$

기하평균은 산술평균보다도 항상 작으므로

$$ 0 \le Q \le 1 $$

의 범위를 갖는다. 그리고 $Q$의 complement로

$$P = 1-Q = \frac{(s_{xx}-s_{yy})^2 + 4 s_{xy}^2}{(s_{xx}+s_{yy})^2};$$를 정의할 수 있는 데 $0 \le P \le 1$이다. $Q$와 $P$의 의미는 무엇인가? 자세히 증명을 할 수 있지만 간단히 살펴보면 한 지점에서 $Q \rightarrow 1$이려면 $λ_{1} \approx λ_{2}$이어야 하고, 이는 두 주축이 동등하다는 의미이므로 그 점에서는 방향성이 없는 코너의 특성을 갖게 된다. 반대로 $Q \rightarrow 0$이면 강한 방향성을 갖게 되어서 edge(ridge) 특성을 갖게 된다.

실제적인 응용으로는 지문 인식에서 지문 영역을 알아내거나 (이 경우는 상당이 큰 윈도를 사용해야 한다) 또는 이미지 텍스쳐 특성을 파악하기 위해서는 이미지를 작은 블록으로 나누고 그 블록 내의 미분 연산자의 균일성을 파악할 필요가 있는데 이 차원이 없는 측도는 이미지의 상태에 상관없이 좋은 기준을 주게 된다.

참고 논문:

Image field categorization and edge/corner detection from gradient covariance
Ando, S.
Pattern Analysis and Machine Intelligence, IEEE Transactions on
Volume 22, Issue 2, Feb 2000 Page(s):179 - 190

** 네이버 블로그 이전;

아날로그 신호로부터 디지털 데이터로 바꿀 때 보통 시간당(sound) 또는 공간적 거리당(image) 일정한 횟수로  데이터를 수집한다. 이 비율을 sampling rate이라 한다. 그런데 sampling rate을 바꾸어야 할 경우가 있다. 1초마다 한 번씩 데이터를 모았는데 실제로 0.5초마다 수집된 데이터가 필요가 한 경우에는 기존의 수집된 데이터와 데이터 사이에서 원 아날로그 신호 값을 알아내야 한다. sampling rate을 바꿀 때 기존의 데이터를 이용해서 원 아날로그 신호에서 수집한 것처럼 새로운 데이터를 만들어 내는 과정이 resampling이다. 이미지에서는 원본을  확대하는 up sampling이나 축소하는 down sampling 등이 resampling의 일종이다. resampling을 하기 위해서는 기존 데이터를 이용해서 데이터와 데이터 사이에서의 값을 추정하는 interpolation이 필요하게 된다. 많이 사용하는 interpolation 방법으로는 nearest neighbor interpolation, linear interpolation, cubic interpolation, spline interpolation 등이 있다.  여기서는 cubic interpolation을 이용한 이미지 resampling을 살펴본다.

$x$축 위의 4 점 $x=-1$, $x=0$, $x=1$, $x=2$에서 값이 $p_0$, $p_1$, $p_2$, $p_3$으로 주어진 경우 $0 \le x \le1$ 구간에서 cubic interpolation 함수는 

$$ f(x)=ax^3 +b x^2 + cx +d$$

라고 쓰면, 우선 $x=0$에서 $p_1$, $x=1$에서 $p_2$를 지나므로

$$f(0) = d = p_1, \quad\quad f(1) = a + b+ c+d = p_2$$

를 만족해야 하고, 양끝에서 도함수 값은 주변 입력점의 평균변화율로 주면, 다음식을 만족한다.

$$  f'(0) = c = \frac{p_2 - p_0}{2}, \quad \quad f'(1) = 3a + 2b + c =\frac{p_3 - p_1}{2}$$

이를 이용해서 $a,b,c,d$를 구하면,

$$a= \frac{1}{2}( -p_0 +3 p_1 -3p_2 +p_3), ~b = \frac{1}{2}(2p_0 -5 p_1 + 4p_2 -p_3),~ c=\frac{1}{2} (-p_0 +p_2), ~d=p_1. $$

따라서 cubic interpolation 함수는

\begin{align} f(x) &=\frac {1}{2}\big(- x^3 + 2x^2 -x\big) p_0+\frac {1}{2}\big(3x^3 - 5x^2 +2 \big) p_1\\&+\frac {1}{2}\big( -3x^3 + 4 x^2 + x\big) p_2 +\frac {1}{2}\big( x^3 - x^2 \big) p_3 \end{align}

로 쓰인다. 이는 다음처럼 정의된 kernel 함수를(Catmull-Rom filter) convolution한 형태로도 표현할 수 있다.

$$ K(x) = \left\{ \begin{array}{ll} \frac{1}{2} (3|x|^3 - 5|x|^2 +2) &  |x| <1 \\  \frac{1}{2}(-|x|^3+5|x|^2 -8|x|+4 ) & 1 \le |x|<2 \\ 0 & \text{otherwise} \end{array} \right. $$

$$  f(x) = p_0 K(-1-x) + p_1 K(-x) + p_2 K(1-x) + p_3 K(2-x),\quad (0\le x <1)$$

여기서는 좀 더 간단히 표현하기 위해서 $$ f(x) = m_0 p_0 + m_1 p_1 + m_2 p_2 + m_3 p_3 $$ 로 표현한다. 물론 $m_i$는 $0\le x <1$의 함수이다. $f(x)$는 그림에서 보는 것처럼 $p_1$과 $p_2$를 통과한다.

구간 끝에서 도함수를 $$f'(0) = \frac{p_2 - p_0}{2}, ~f'(1) = \frac{p_3 - p_1}{2}$$로 설정했기 때문에 인접 구간(eg. $[p_2,p_3]$)에서의 interpolation과 경계에서 같은  도함수를 가지게 되어 곡선이 꺽임없이 부드럽게 연결이 되는 특성을 가진다(예를 들면, $p_2$와 $p_3$을 보간하는 cubic interpolation함수는 $p_2$에서 미분값이 $(p_3 - p_1)/2$인데, 이는 $p_1,p_2$구간의 interpolation  함수가 $p_2$에서 가지는 미분값과 같다).

이를 2차원 공간으로 확장하면, 평면 위의 격자점 $\{(x, y)| x, y=-1,0,1, 2\}$에서 16개의 값 $\{ p_{ij} |i, j=0,1,2, 3\}$가 주어진 경우 정사각형 영역 $ 0\le x, y \le 1$에서 bicubic interpolation은 먼저 4개의 행 $y=-1,0,1,2$에서 각각 $x$ 방향의 cubic interpolation 된 값을 구한다:

\begin{align} q_0 =f(x, y=-1)= m_0 p_{00} + m_1 p_{10} + m_2 p_{20} + m_3 p_{30} \\ q_1 =f(x, y=0)= m_0 p_{01} + m_1 p_{11} + m_2 p_{21} + m_3 p_{31} \\ q_2 =f(x, y=1)= m_0 p_{02} + m_1 p_{12} + m_2 p_{22} + m_3 p_{32} \\ q_3 =f(x, y=2)= m_0 p_{03} + m_1 p_{13} + m_2 p_{23} + m_3 p_{33} \end{align}

4개 행에서 구해서 값을 이용해서 다시 $y$ 방향으로 cubic interpolation 한 값을 구하면 bicubic interpolation이 완성된다: $$ q =f(x,y)= q_0 n_0 + q_1 n_1 + q_2 n_2 + q_3 n_3$$ 여기서 $n_i$는 $m_i$에서 $x$를 $y$로 치환한 값이다.

원본 이미지를 확대/축소 또는 변형 등의 변환 과정을 거쳐서 출력 이미지를 만들 때 원본 이미지의 픽셀 값을 resampling 하는 과정이 들어온다. 이때 원본 이미지의 픽셀과 픽셀 사이의 값이 필요한 경우가 발생하여 적절한 interpolation이 필요한데 많이 쓰이는 interpolation 중의 하나가 bicubic interpolation이다. 아래는 bicubic interpolation을 이용한 이미지 resampling을 수행하는 코드다.

interpolation에서 실수(float) 연산을 하지 않고 정수 연산만 사용하면서도 일정한 정밀도로 소수점 아래 부분의 정보를 유지하기 위해 나누기 전에 미리 256을 곱한 후에 나눈다(256 대신에 적당히 큰 수를 선택해도 되는데 2의 지수승을 잡으면 곱셈/나눗셈을 shift 연산으로 대체할 수 있는 장점이 있다). 이렇게 하면 나눈 몫에 $\tt 0xFF$ 비트 마스크를 적용할 때 남는 부분이 소수점 아랫부분을 표현한다. 정밀도는 $1/256$ 까지다. 중간 과정에서 소수점 이하 부분끼리 곱셈을 할 때는 항상 $256$으로 나누어서 $255$ 이하가 되도록 만들어야 한다. 최종 결과에서도 다시 $256$으로 나누기를 해주어야 된다. bicubic인 경우는 $x/y$ 양방향으로 적용하므로 $256\times 256$으로 나누어야 하고 cubic interpolation 계수에서 $1/2$이 들어오므로 추가로 $4$만큼 더 나누어 주어야 한다(코드의 마지막 결과에서 shift 연산 "$\tt >> 18$"이 들어온 이유다).

bicubic interpolation을 적용할 때 $\tt y=0$이나 $\tt x=0$에서는 이전 행이나 열이 없으므로 자신을 반복하는 방식으로 처리해 주어야 하고, 또 $\tt y=rows-2, rows-1$이나 $\tt x=cols-2, cols-1$일 때도 비슷한 처리가 있어야 한다.

int resample_bicubic ( BYTE *src, int cols, int rows,
                       BYTE *des, int newcols, int newrows ) {
    if (cols < 4 || rows < 4)
        return resample_bilinear(src, cols, rows, des, newcols, newrows);

    int ixn = cols - 4;       
    BYTE *pa, *pb, *pc, *pd;
    for (int j = 0; j < newrows; j++) {
        int yy = ( ( j * rows ) << 8 ) / newrows;
        int yp = yy >> 8;                        // src pixel y-position;
        int dy = yy & 0xFF;
        int dy2 = (dy * dy) >> 8;
        int dy3 = (dy2 * dy) >> 8;
        int n0 = -dy3 + 2 * dy2 - dy;
        int n1 = 3 * dy3 - 5 * dy2 + 2 * 256;
        int n2 = -3 * dy3 + 4 * dy2 + dy;
        int n3 = dy3 - dy2;
        //
        pb = src + yp * cols;                  //current line;
        if (yp == 0) pa = pb;
        else         pa = pb - cols;           //previous line;
        
        if (yy < rows - 2) {
            pc = pb + cols;                    //next line;
            pd = pc + cols;                    //next-next line;
        } else if (yp < rows - 1) {
            pc = pb + cols;        
            pd = pc;
        } else 
            pd = pc = pb;

        for (int i = 0; i < newcols; i++) {
            int xx = ( ( i * cols ) << 8 ) / newcols;
            int xp = xx >> 8;                    // src pixel x-position;
            int dx = xx & 0xFF;
            int dx2 = ( dx * dx ) >> 8;
            int dx3 = ( dx2 * dx ) >> 8;
            int m0 = -dx3 + 2 * dx2 - dx;
            int m1 = 3 * dx3 - 5 * dx2 + 2 * 256;
            int m2 = -3 * dx3 + 4 * dx2 + dx;
            int m3 = dx3 - dx2;
            int p = (xp == 0) ? 0 : (xp < ixn) ? xp - 1: ixn;    // p+3 <= ixn+3=cols-1;
            int a = ((m0 * pa[p] + m1 * pa[p + 1] + m2 * pa[p + 2] + m3 * pa[p + 3]) * n0 +
                     (m0 * pb[p] + m1 * pb[p + 1] + m2 * pb[p + 2] + m3 * pb[p + 3]) * n1 +
                     (m0 * pc[p] + m1 * pc[p + 1] + m2 * pc[p + 2] + m3 * pc[p + 3]) * n2 +
                     (m0 * pd[p] + m1 * pd[p + 1] + m2 * pd[p + 2] + m3 * pd[p + 3]) * n3) >> 18;
            *des++ = (a > 0xFF) ? 0xFF: (a < 0) ? 0: a;
        }
    }
	return 1;
}

bicubic interpolation을 하기 위해서는 4점이 필요하므로 폭이나 높이가 이보다 작은 경우는 bilinear interpolation을 사용한다. 다음 코드는 fixed-point bilinear interpolation을 구현한 코드다.

int resample_bilinear(BYTE *src, int cols, int rows,
                      BYTE *des, int newcols, int newrows ) {
    for (int j = 0; j < newrows; j++ ) {
        int yy = ((j * rows ) << 8) / newrows;
        int yp = yy >> 8;   // integer part; src y-position;
        int dy = yy & 0xFF; // fractional part;
        BYTE *curr = src + yp * cols;
        BYTE *next = curr + cols;
        for (int i = 0; i < newcols; i++) {
            int xx = ((i * cols ) << 8) / newcols;
            int xp = xx >> 8;       //src x-position;
            int dx = xx & 0xFF;
            int p00 = curr[xp];
            int p10 = curr[xp + 1];
            int p01 = next[xp];
            int p11 = next[xp + 1];
            int val = ((p11 * dx + p01 * (256 - dx)) * dy
                    + (p10 * dx + p00 * (256 - dx)) * (256 - dy)) >> 16;
            *des++ = val > 255 ? 0xFF: val < 0 ? 0 : val;
        }
    }
    return 1;
}

kipl.tistory.com/55