BiLinear Interpolation of Image in C

구글을 뒤져서 다음의 코드를 찾아 냈다.

https://gist.github.com/folkertdev/6b930c7a7856e36dcad0a72a03e66716

bilinear interpolation in C

 

오랜동안 살아있는 링크로 알고 있지만, 혹시나 언제 지워질지 모르니 복사를 해 둔다.

/// Sources:
//
//  - https://chao-ji.github.io/jekyll/update/2018/07/19/BilinearResize.html
void bilinear_interpolation(float *data, uint32_t input_width,
                            uint32_t input_height, uint32_t output_width,
                            uint32_t output_height, float *output) {
    float x_ratio, y_ratio;

    if (output_width > 1) {
        x_ratio = ((float)input_width - 1.0) / ((float)output_width - 1.0);
    } else {
        x_ratio = 0;
    }

    if (output_height > 1) {
        y_ratio = ((float)input_height - 1.0) / ((float)output_height - 1.0);
    } else {
        y_ratio = 0;
    }

    for (int i = 0; i < output_height; i++) {
        for (int j = 0; j < output_width; j++) {
            float x_l = floor(x_ratio * (float)j);
            float y_l = floor(y_ratio * (float)i);
            float x_h = ceil(x_ratio * (float)j);
            float y_h = ceil(y_ratio * (float)i);

            float x_weight = (x_ratio * (float)j) - x_l;
            float y_weight = (y_ratio * (float)i) - y_l;

            float a = data[(int)y_l * input_width + (int)x_l];
            float b = data[(int)y_l * input_width + (int)x_h];
            float c = data[(int)y_h * input_width + (int)x_l];
            float d = data[(int)y_h * input_width + (int)x_h];

            float pixel = a * (1.0 - x_weight) * (1.0 - y_weight) +
                          b * x_weight * (1.0 - y_weight) +
                          c * y_weight * (1.0 - x_weight) +
                          d * x_weight * y_weight;

            output[i * output_width + j] = pixel;
        }
    }
}

이해 하기 쉽게 깔끔히 작성된 코드라서 분석을 해보기로 했다.

 

크기변경

분석~

새로운 이미지의 픽셀에 해당하는 원본의 픽셀 위치를 계산한다.

계산된 위치의 픽셀[1]과, 바로 오른쪽[2]/ 바로 아래[3] / 바로 오른쪽-아래픽셀[4], 총 4개의 픽셀에 거리에 따른 가중치를 적용해서 새로운 픽셀값을 계산해 낸다. 바로 옆에 붙어 있는 픽셀로만 계산하도록 되어 있다.( 비율과 상관없이 멀리 떨어진 픽셀은 신경도 안쓴다.)

첫 픽셀의 경우[마지막 픽셀도] 옆에 붙어 있는 픽셀을 참조하진 않지만 참조시켜도 가중치가 첫 픽셀에만 집중되기 때문에 결과는 같다. 하지만 마지막 픽셀의 경우 원본의 참조 범위가 넘어가기 때문에 참조하는 픽셀은 이웃한 픽셀없이 해당하는 픽셀로 한정 시키고 있다. 생성할 이미지가 원본에 비해 아주 클 경우 참조하는 픽셀은 1개만 되는 경우도 많다. 하지만 이런 경우가 테스트를 해보니 그리 많지 않아 최적화를 위해 if문으로 분기시킬 이유는 없었다.

이웃한 픽셀을 참조할지 말지는 xr값에 소수점이 있느냐 없느냐로 구분되고 이 소수점 부분이 가중치 계산에 이용된다.

 

변환

 

 

 

 

가중치가 0.5, 0.5일경우 균등하게 적용된다.

 

 

하지만.. 위에 코드는 한가지 색밖에 처리할 수 없다. RGBA같이 4개의 색을 처리하려면 코드를 좀 수정해줘야 한다.

위코드를 바탕으로 RGBA/32비트를 처리하는 코드로 재작성한 코드는 다음과 같다.

 

imageDsp.h

#include <inttypes.h>
struct IMAGE_RESAMPLE_PARAM {
    int32_t   srcW;
    int32_t   srcH;
    int32_t   srcStride; // in pixel
    uint32_t* srcPixels;
    
    int32_t  dstW;
    int32_t  dstH;
    int32_t  dstStride; // in pixel
    uint32_t* dstPixels;
};

srcStride와 dstStride는 byte단위가 아니고 픽셀단위(byte/4)로 설정하게 만들었다. 32비트만 쓸 것이기에 이렇게 작성되었다.

 

imageDsp.cpp

#include <math.h>

bool ImageResampleFloat( const IMAGE_RESAMPLE_PARAM& param ) {
    
    const float x_ratio = param.dstW > 1 ? ( param.srcW - 1.0F ) / ( param.dstW - 1.0F ) : 0.0F;
    const float y_ratio = param.dstH > 1 ? ( param.srcH - 1.0F ) / ( param.dstH - 1.0F ) : 0.0F;

    for (int y = 0; y < param.dstH ; y++) {
        const float yr = y * y_ratio;
        const float y_l = floorf( yr );
        const float y_h = ceilf( yr );
        const float y_w = yr - y_l;

        uint32_t* srcPixelY1 = param.srcPixels + param.srcStride * (int)y_l;
        uint32_t* srcPixelY2 = param.srcPixels + param.srcStride * (int)y_h;
        uint8_t*  dstPixel   = (uint8_t*)( param.dstPixels + param.dstStride *  y );
       
        
        for (int x = 0; x < param.dstW ; x++ ) {
            const float xr = x * x_ratio;
            const float x_l = floorf( xr );
            const float x_h = ceil( xr );
            const float x_w = xr - x_l;
                        
            
            const float f4 = x_w * y_w;
            const float f1 = 1.0F - x_w - y_w + f4;
            const float f2 = x_w - f4;
            const float f3 = y_w - f4;
            
            const uint8_t* pxa = (uint8_t*)( srcPixelY1 + (int)x_l );
            const uint8_t* pxb = (uint8_t*)( srcPixelY1 + (int)x_h );
            const uint8_t* pxc = (uint8_t*)( srcPixelY2 + (int)x_l );
            const uint8_t* pxd = (uint8_t*)( srcPixelY2 + (int)x_h );
            
            dstPixel[0] = (uint8_t)( f1 * pxa[0] + f2 * pxb[0] + f3 * pxc[0] + f4 * pxd[0] );
            dstPixel[1] = (uint8_t)( f1 * pxa[1] + f2 * pxb[1] + f3 * pxc[1] + f4 * pxd[1] );
            dstPixel[2] = (uint8_t)( f1 * pxa[2] + f2 * pxb[2] + f3 * pxc[2] + f4 * pxd[2] );
            dstPixel[3] = (uint8_t)( f1 * pxa[3] + f2 * pxb[3] + f3 * pxc[3] + f4 * pxd[3] );
            
            dstPixel += 4;
        }
    }
    
    return true;
}

RGB/GA등 1~3개의 색을 처리하려면 위의 코드를 좀 수정해서 쓰면 될 것이다. 필자의 경우 모든이미지를 RGBA[32비트=8bits*4colors]로만 처리하기 때문에 특별히 수정할 일은 없었다.

 

최적화

이미지를 다루는 작업은 워낙 처리할 양이 많아서 최적화작업을 해주면 좋다.

1. 연산을 거의 전부 float으로 처리하기 때문에 int형을 이용한 계산을 하면 속도가 빨라 질 수 있다.

2. 이중루프를 이용하는 코드라 내부 루프의 내용에서 1Clock만 줄이더라도 1024x1024이미지의 경우 약 백만Clock을 줄일수가 있다. 내부루프에서 외부로 뺄 수 있는 것은 빼준다.

3. CPU를 이용한 코드라서 openmp를 이용하면 cpu수만큼 속도가 빨라질 수 있다.

4. xr변수의 경우 일정하게 증가하기 때문에 곱하기 연산에서 xr변수를 계속 증가시키는 방식을 적용한다.

    yr변수도 같은 방식으로 처리할 수 있지만 openmp를 사용하기 위해서 그냥 놔두었다.

5. 곱하기연산을 대체할 수 있는 것은 최대한 대체해 보았다.

 

 

위의 최적화를 적용한 코드는 다음과 같다.

bool ImageResample( const IMAGE_RESAMPLE_PARAM& param ) {
    
    const int32_t x_ratio = param.dstW > 1 ? ( (param.srcW - 1) << 16 ) / (param.dstW - 1) : 0;
    const int32_t y_ratio = param.dstH > 1 ? ( (param.srcH - 1) << 16 ) / (param.dstH - 1) : 0;
    
#ifdef NDEBUG
    #pragma omp parallel for schedule( static )
#endif
    for (int y = 0; y < param.dstH ; y++) {
        const int32_t yr = y * y_ratio;
        const int32_t y_l = yr >> 16;
        const int32_t y_w = yr & 0xffff;
        const int32_t y_h = y_w ? y_l + 1 : y_l; // same to y_l or just below pixel
        assert( y_h < param.srcH );

        uint32_t* srcPixelY1 = param.srcPixels + param.srcStride * y_l;
        uint32_t* srcPixelY2 = param.srcPixels + param.srcStride * y_h;
        uint8_t*  dstPixel   = (uint8_t*)( param.dstPixels + param.dstStride *  y );
       
        int32_t xr = 0;
        
        for (int x = 0; x < param.dstW ; x++ ) {
            
            const int32_t x_l = xr >> 16;
            const int32_t x_w = xr & 0xffff; 
            const int32_t x_h = x_w ? x_l + 1 : x_l;
            assert( x_h < param.srcW );
                        
            const int32_t f4 = (int32_t)( ( (int64_t)x_w * (int64_t)y_w ) >> 16LL );
            const int32_t f3 = y_w - f4;
            const int32_t f2 = x_w - f4;
            const int32_t f1 = 65536 - y_w - f2; // 65536 - y_w - x_w  + f4;
                                
            const uint8_t* pxa = (uint8_t*)( srcPixelY1 + x_l );
            const uint8_t* pxb = (uint8_t*)( srcPixelY1 + x_h );
            const uint8_t* pxc = (uint8_t*)( srcPixelY2 + x_l );
            const uint8_t* pxd = (uint8_t*)( srcPixelY2 + x_h );
            
            dstPixel[0] = (uint8_t)( ( f1 * pxa[0] + f2 * pxb[0] + f3 * pxc[0] + f4 * pxd[0] ) >> 16 );
            dstPixel[1] = (uint8_t)( ( f1 * pxa[1] + f2 * pxb[1] + f3 * pxc[1] + f4 * pxd[1] ) >> 16 );
            dstPixel[2] = (uint8_t)( ( f1 * pxa[2] + f2 * pxb[2] + f3 * pxc[2] + f4 * pxd[2] ) >> 16 );
            dstPixel[3] = (uint8_t)( ( f1 * pxa[3] + f2 * pxb[3] + f3 * pxc[3] + f4 * pxd[3] ) >> 16 );
            dstPixel += 4;
           
            xr += x_ratio;
        }
    }
    
    return true;
}

openMP를 쓰지 않을 사람은 openMP부분만을 제거하면 된다.

 

float연산을 이용한 코드보다 2~3배정도 빨라지긴 했다.

 

컴파일에 필요한 flag와 option은 다음과 같다.

<Compile flag>
cpp ver : -std=c++17
openMP  : -fopenmp
simd    : -msse4.2 or -msse or 안드로이드라면 neon으로 적용하면 된다.

<Linker option>
math   : -lm
openMP : -lgomp

 

 

이것저것 만들면서 들었던 생각들

1. unrolling을 적용해 봤지만 내부 루프에 연산이 많아 너무 의미가 없어서 빼버렸다.

2. xr값의 변화를 incremental방식을 적용할 경우 int는 문제가 별로 없지만, float의 경우 마지막 픽셀 값이 정확히 떨어지지 않아서 x_h계산할 때 참조 메모리를 넘어가 버린다. [FLOAT epsilon]을 빼봤지만 차이가 없었다. 일정량을 빼버리기엔 계산하기 싫어서 그냥 포기했다.

2 - 1 . int연산의 경우 약간의 오차가 발생하긴 했지만 문제가 될 수준은 아니었다. x축의 마지막 픽셀에 오차가 좀 있다.

3. 어셈블러로 코딩할까 고민해 보았지만, 32bits/neon등 할 작업이 많아서 그냥 포기했다. compile된 코드를 보니 2배는 빠르게 만들 것같지만... 수정하려면 그 것도 머리 아프기도 하고..