Geometry & Recognition

// 배열첨자(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);
}

방 안에 뜨거운 물체가 놓이면 그 열은 시간이 지나면서 사방으로 퍼져나간다. 열원에서 주위로 열이 시간에 따라 전달이 되는 프로세스를 기술하는 방정식이 열(전도)방정식이다. 이것을 이미지 필터링에서도 적용이 가능하다. 이미지에 나타난 각각의 노이즈를 열원으로 보면 주변으로 열의 전달은 열원의 온도를 낮추게 되는데 이를 이미지에서 노이즈를 제거하는 filtering 과정으로 생각할 수 있다. 그러나 열방정식은 방향성이 없는 등방성 확산(isotropic diffusion process)을 기술한다. 이미지에 열전도 방정식을 바로 적용하면 이미지의 블러링과 유사한 효과를 만든다. filtering를 하더라도 에지를 보존해야 이미지의 정보 손실이 적게 된다. 에지 영역에서는 gradient의 공분산 행렬($\Sigma$)의 고윳값의 차이가 크고, 노이즈나 corner 근방에서는 두 고윳값의 크기가 비슷하다. 따라서, 에지를 보존하기 위해서는 공간으로 확산하는 역할을 하는 이미지의 Laplace 항 ($\nabla^2 I = I_{xx} + I_{yy}$)을 gradient의 공분산 행렬의 고유 벡터와 고유치, 그리고 이미지의 hessian 행렬 등을 조합하여서 불변인 형태를 만들어야 한다. 여러 가능한 후보군 중에 하나의 예를 들면

$$ I(t+dt) - I(t) = f_1(s) \text{Tr}({\bf  U U^T H}) + f_2(s) \text{Tr}({\bf V V^T H}) $$

$$ {\bf H}=\begin{pmatrix} I_{xx} & I_{xy } \\ I_{xy} & I_{yy} \end{pmatrix}, \quad \nabla^2 I = \text{Tr} ({\bf H}) $$

여기서, $\bf U, V$는 그라디엔트 공분산 행렬의 두 고유 벡터이고, $s=\lambda_1 + \lambda_2$는 고유값의 합이다. $f_1$과 $f_2$의 함수 형태를 고유치가 큰 방향에서는 억제되고, 고유치가 작은 방향에서는 덜 억제가 되도록 적절히 선택하는 경우에 에지를 보존하는 필터링의 효과를 얻을 수 있다.

아래의 예는 이 식을 써서 300번 iteration을 한 결과이다. 원본 이미지는 첨부한 논문에서 추출하였다.

https://kipl.tistory.com/563