Fant's Resampling

Fant's Resampling 알고리즘은 영상의 확대나 축소 과정에서 픽셀의 밝기값을 단순히 보간(interpolation)하는 대신, 픽셀 면적의 중첩(overlap)을 고려하여 새로운 픽셀 값을 계산하는 영상 리샘플링(resampling) 방법이다. 특히 영상을 축소할 때 발생하는 aliasing을 효과적으로 줄일 수 있어 고품질의 영상 크기 변환에 널리 사용된다.

기본 아이디어는 입력 영상의 각 픽셀을 일정한 면적을 갖는 직사각형으로 생각하고, 출력 영상의 각 픽셀과 입력 픽셀 사이에 겹치는 면적의 비율만큼 밝기값을 가중 평균하는 것이다. 즉, 출력 픽셀은 자신과 겹치는 입력 픽셀들의 평균 밝기를 나타내므로 영상의 전체 에너지가 잘 보존된다.

Fant 알고리즘은 이러한 면적 적분(area integration)을 직접 계산하지 않고, 입력 픽셀과 출력 픽셀의 경계를 순차적으로 이동시키면서 중첩 면적만을 누적하는 방식으로 계산한다. 따라서 계산량이 출력 픽셀 수에 비례하는 선형 시간(linear time)에 수행되며, 메모리 사용량도 매우 적다.

2차원 영상에서는 먼저 가로 방향으로 1차원 Resampling을 수행한 후, 이어서 세로 방향으로 같은 과정을 적용하는 separable 방식으로 구현하는 것이 일반적이다. 이 방법은 2차원 면적 적분을 직접 계산하는 것보다 계산량이 훨씬 적으면서도 동일한 결과를 얻을 수 있다.

Fant's Resampling은 구현이 간단하고 계산 효율이 높으며, 특히 축소 영상에서 aliasing을 효과적으로 억제할 수 있어 다양한 영상 처리 라이브러리와 그래픽 응용 프로그램에서 오랫동안 사용되어 온 대표적인 resampling 알고리즘이다.

// 배열첨자(dj)와 픽셀의 실제위치(srcindex, dstindex)를 따로 분리하여서
// 열방향을 따라서 작업을 하더라도 메모리 복사가 필요없이 처리가 가능하도록 하였음.
BOOL 
resampleRGB(BYTE *src, BYTE* dst, int length,      // 
               int stride,                         // RGB = 3(행방향), 행의 총 BYTE수(열방향) 
               double outpos[length + 1])          // forward map: src -> dst; 
{
    std::vector<double> inpos(length+1);
    // si = src index;
    // dj = dst index;
    // Find inverse map; dst--> src;
    for (int si = 0, dj = 0; dj < length; dj++) {
        while (si < (length - 1) && outpos[si + 1] < dj) si++; // max(si)=length-1;
        if (si < length - 1)
            inpos[dj] = si + (dj - outpos[si]) / (outpos[si + 1] - outpos[si]);
        else // max(inpos[dj])=lenght-1;
            inpos[dj] = si + 1 ;
    }
    inpos[length] = length;

    double inseg  = 1;            //첫 입력 픽셀은 완전히 사용가능;
    double outseg = inpos[1];       //첫 출력을 위해서 필요한 입력의 양은 
                                    //inpos[1]-inpos[0]이나 inpos[0]=0으로 함;
    double sizefac = outseg ;
    int srcindex = 0;
    int dstindex = 0;
    //첫번째 픽셀;
    int b = src[0], g = src[1], r = src[2];
    srcindex += stride;
    //두번째 픽셀;
    int nextb = src[srcindex+0], 
        nextg = src[srcindex+1], 
        nextr = src[srcindex+2];
    srcindex += stride;
    //
    double bsum = 0, gsum = 0, rsum = 0;
    for (int dj = 1; dj < length; ) {
        // linear interpolation;
        double bintensity = inseg * b + (1 - inseg) * nextb ;
        double gintensity = inseg * g + (1 - inseg) * nextg ;
        double rintensity = inseg * r + (1 - inseg) * nextr ;
        //
        if (inseg < outseg) {//output cycle;
            // accumulation of weighted contrib;
            bsum += inseg * bintensity ; 
            gsum += inseg * gintensity ; 
            rsum += inseg * rintensity ; 
            //기존에 현재의 픽셀을 추가하므로, 출력을 위한 입력픽셀의 소요량이 줄어든다.
            outseg -= inseg ;
            
            //현재의 현재의 픽셀값을 갱신한다.
            b = nextb; g = nextg; r = nextr;
            //새로 들어올 입력픽셀을 준비한다.
            inseg = 1.0 ;
            // srcindex < endindex-2;
            nextb = src[srcindex+0];
            nextg = src[srcindex+1];
            nextr = src[srcindex+2];
            srcindex += stride ;
        } else { //input cycle;
            // accumulation of weighted contrib;
            bsum += outseg * bintensity ;
            gsum += outseg * gintensity ;
            rsum += outseg * rintensity ;
            //hack;
            if (sizefac == 0) sizefac = 1;
            //set dst pixel if inpos[dj]>=0; dj=1,2,3,....;
            //src가 dst 내부로 들어가는 경우임;
            if (inpos[dj - 1] >= 0) {
                //x = 0, 1, 2...
                // 출력픽셀을 만드는데 sizefac 만큼 입력픽셀이 들어갔으므로 나누어 주어야 한다.
                dst[dstindex+0] = (BYTE)min(bsum / sizefac, 0xFF);
                dst[dstindex+1] = (BYTE)min(gsum / sizefac, 0xFF);
                dst[dstindex+2] = (BYTE)min(rsum / sizefac, 0xFF);
            }
            dstindex += stride ;

            // reset accumulator for next output;
            bsum = gsum = rsum = 0; 
            // source has been consumed outseg fraction; 
            // 현재의 입력픽셀이 다음 출력픽셀에 기여할 수 있는 남아 있는 양;
            inseg -= outseg ;
            // set new outseg; = 다음 출력픽셀을 완성하는데 필요한 입력픽셀의 양;
            outseg = inpos[dj + 1] - inpos[dj] ;
            // 출력 픽셀을 완성하는데 필요한 입력픽셀의 양(outseg는 다음 출력시까지 계속 변하므로 
            // 처음 세팅할 때 기억해 두어야 한다;
            sizefac = outseg ;
            dj++ ;
        }
    }
    return TRUE;
}

void test_fant(CRaster& raster, CRaster& out) {
    CSize sz = raster.GetSize() ;
    ASSERT(raster.GetBPP() == 24) ;
    out.SetDimensions(sz, raster.GetBPP()) ;
    std::vector<double> hpos(sz.cx + 1);
    std::vector<double> vpos(sz.cy + 1);
    // horizontal test ;
    for (int x = 0; x < sz.cx; x++) hpos[x] = 4 * double(x) / 3 ;   
    for (int y = 0; y < sz.cy; y++) {
        BYTE *psrc = (BYTE *)raster.GetLinePtr(y) ;
        BYTE *pdst = (BYTE *)out.GetLinePtr(y) ;
        resample(psrc, pdst, sz.cx, 3, hpos) ;
    }
    //vertical test;
    CRaster tmp ;
    tmp.SetDimensions(sz, raster.GetBPP());
    for (int x = 0; x < sz.cx; x++){
        //for each x;
        double starty = (double(x) * sz.cy) / (3 * sz.cx) ;
        double endy   = sz.cy - (double(x) * sz.cy) / (3 * sz.cx) ;
        for (int y = 0; y < sz.cy; y++)
            vpos[y] = starty + (endy - starty) * y / sz.cy ;
        BYTE *psrc = (BYTE *)out.GetLinePtr(0) + 3 * x ;
        BYTE *pdst = (BYTE *)tmp.GetLinePtr(0) + 3 * x ;
        resample(psrc, pdst, sz.cy, -(int)raster.GetBytesPerLine(), vpos);
    }
    out.SwapBuffer(tmp);
}