Fant's 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);
}