Lanczos kernel:  a sinc function windowed by the central lobe of a second, longer, sinc function(wikipedia)

$$K(x) = \left\{ \begin {array}{ll}   \text {sinc}(\pi x) \text {sinc}\Big( \frac {\pi x}{a}\Big), & |x| < a \\ 0 ,& |x| \ge a\end {array} \right. = \text{sinc}(\pi x) \text{sinc}\Big( \frac{\pi x}{a}\Big)\text{Box}(|x|<a) $$

$$H(\omega)=\frac{1}{2\pi^2\sqrt{2\pi}}\Big( -(2\pi-3\omega)\text{Si}(2\pi-3\omega)+(4\pi -3\omega)\text{Si}(4\pi-3\omega)\\+(2\pi +3\omega)\text{Si}(2\pi+3\omega)+(4\pi+3\omega)\text{Si}(4\pi+3\omega)\Big).$$

여기서 sine intergral은 다음과 같이 정의된다:

$$\text{Si}(x)=\int_0^x \frac{\sin(t)}{t} dt.$$

$x=0$ 근방에서는 Tayor 전개식을 이용하는 것이 더 효율적인다.

\

Lanczos kernel을 사용해서 일정한 간격으로 sampling 된 데이터 $\{f_i \}$를 보간하려고 할 때 보간함수 $F(x)$는 Lanczos kernel과 데이터의 convolution으로 표현할 수 있다.

$$ F(x) = \sum_{i \in \text{window}} f_i K(x - i) $$

Interpolation 함수는 영상의 resampling 과정에서 사용할 수 있다. 주어진 이미지를 확대하거나 축소하려면 기존 이미지의 픽셀과 픽셀 사이 구간에서 픽셀 값을 알아내야 하는데, 이때 interpolation 함수를 이용하여 필요한 데이터를 얻는 과정인 resampling 한다.

 

다음 코드는 Lanczos3 kernel을 사용해서 이미지를 확대하거나 축소할 때 행에 대해서 필요한 resampling을 수행한다. 여기서는 kernel의 중심이 소스 이미지의 픽셀과 픽셀의 가운데에 놓이도록 설정하였다($-0.5$의 의미). 이는 작은 영상을 크게 확대할 때 가장자리에서 왜곡이 되는 것을 방지하기 위해서 인데 큰 영상을 확대/축소할 때는 영향이 없다. 이 interpolation은 separable이므로 2차원의 경우 행방향으로 진행한 후 그 결과를 다시 열 방향으로 진행하는 2-pass 방식으로 사용해도  된다.

// windowed sinc(x); window=sinc(x/3);
static double Lanczos3(double x){
    if (x < 0) x = -x; // symmetric;
    x *= PI;
    if (x < 0.01)    return 1. + (- 5./ 27. + 
        (14. / 1215. - 17. * x * x / 45927.) * x * x) * x * x; 
    else if (x < 3.) return 3. * sin(x) * sin(x / 3.) / x / x;
    else     	     return 0;
};
// interpolation in the horizonal direction: single channel or gray image;
// x = pixel position in a destination image;
double Lanczos3_line(BYTE *src_line, int srcWidth, int x, double xscale) {
    double halfWidth;
    if (xscale > 1.) halfWidth = 3.;
    else             halfWidth = 3. / xscale;

    double centerx = double(x) / xscale - 0.5;  //center loc of filter in the src_line;
    int left  = (int)floor(centerx - halfWidth);
    int right = (int)ceil(centerx + halfWidth);
    if (xscale > 1) xscale = 1;
    double s = 0;
    double weightSum = 0;
    for (int ix = left; ix <= right; ix++) {   
        double weight = Lanczos3((centerx - ix) * xscale);
        int xx = ix;         // for ix<0 || ix>=srcWidth: repeat boundary pixels
        xx = min(max(xx, 0), srcWidth - 1));
        s += weight * src_line[xx];
        weightSum += weight;
    }
    return s / weightSum; 
}

 

이차원 이미지의 경우 수평방향의  적용한 후 다시 수직방향으로 적용하면 된다. 

bicubic downsample(1/2): alias 발생

 
 
 
 
 
 
728x90

'Image Recognition > Fundamental' 카테고리의 다른 글

Poisson Image Editing  (0) 2021.08.03
Sampling Theorem  (0) 2021.05.12
Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Posted by helloktk
,

Fant's Algorithm

Image Recognition 2010. 1. 22. 17:38

이미지 resampling 센서에서 받은 신호를 sampling 해서 만든 영상을 확대하거나 줄이는 경우, 또는 기하학적인 변환을 걸쳐서 다시 영상을 생성하는 과정이다. 이 과정에서 출력 영상의 픽셀 좌표 (x, y)에 대응하는 입력 영상의 픽셀 좌표 (u, v)가 요구되는데 일반적으로 이 좌표값은 정수로 주어지지 않는다. 따라서 입력 영상에서 픽셀 값을 얻기 위해서는 (u, v)에 가장 가까운 정수에 해당하는 지점의 픽셀 값을 사용하거나(nearest neighbor interpolation), 아니면 (u, v) 주변의 픽셀 값들을 보간하는 함수를 찾아서 필요한 픽셀 값을 추정할 수 있다: bilinear interpolation, bicubic interpolation,... 

실제 영상 데이터는 하드웨어적으로 스캔라인 방식으로 접근하므로 resampling 알고리즘도 스캔라인 방식에 맞추어서 수평방향으로 처리하고, 그 결과를 다시 수직방향으로 처리하는 separable 알고리즘이 요구된다Fant 알고리즘은 스캔라인 방식이어서 하드웨어적으로 구현이 쉽도록 설계되었다. 

resampling 중에서도 특히 down sampling을 할 때는 입력 픽셀 정보의 손실로 인해서 나타나는 계단 현상이나 모아레 무늬와 같은 alias를 줄이는 알고리즘이 필요하다. Fant 알고리즘에서는 각각의 입력 픽셀들이 출력 픽셀을 만드는데 얼마나 기여하는지를 고려하여 그 기여만큼 가중치를 주어 합산한다. 이렇게 하면 입력 픽셀의 정보의 손실을 줄어 alias 현상을 억제할 수 있다.

Fant 알고리즘은 일반적으로 출력 영상에 해당하는 입력 픽셀의 좌표가 주어지는 역변환(inverse mapping)에 대해서 사용이 된다. (참고: 이곳에서는 입력 영상에 대응하는 출력 영상의 좌표가 주어지는 정변환(forward mapping)에 대한 Fant 알고리즘의 구현이 있다)

참고:

// dst--->src;
// index ++ ;
// pointer += stride ;
/* stride = horizontal_scan=3; vertical_scan=bytes_per_line;*/
void fant_resample24_inverse(const BYTE src[], const int src_len, const int src_stride,
                              BYTE dest[], const int dest_len, const int dest_stride, 
                              const float f[/*dest_len+1*/]) /*dest_pixel=>src_pixel map*/ 
{
    float inseg, outseg ;
    //dest scanline pixel projects into [0, src_len];
    if (f[dest_len] < 0 || f[0] > src_len) return ;
    //advance to;
    for (int x = 0; f[x] < 0; x++) ;
    int xl = x > 0 ? x - 1 :  x ;
    //
    for (x = dest_len; f[x] > src_len; x--) ;
    int xr = (x == dest_len) ? dest_len - 1 : x ;
    //
    float isf = f[xl + 1] - f[xl];  //inverse_size_factor;
    if (f[x] < 0) {
        inseg = 1;
        outseg = isf + f[xl] ;
    } else {
        inseg = 1 - (f[0] - int(f[0]));
        outseg = isf ;
    }
    //src_initial_pos
    int u = f[x] < 0 ? 0 : int(f[xl]) ;
    src += u * src_stride ;
    int bval = src[0];
    int gval = src[1];
    int rval = src[2];
    src += src_stride ;
    u++ ;
    //src_next_pos ;
    int bnext = src[0];
    int gnext = src[1];
    int rnext = src[2];
    src += src_stride ;
    u++;
    //dest_inital_pos;
    dest += xl * dest_stride ;
    float bsum = 0, gsum = 0, rsum = 0;
    for (x = xl ; x <= xr; ) {
        float binten = inseg * bval + (1 - inseg) * bnext;
        float ginten = inseg * gval + (1 - inseg) * gnext;
        float rinten = inseg * rval + (1 - inseg) * rnext;
        if (inseg < outseg) {
            bsum += binten * inseg ;
            gsum += ginten * inseg ;
            rsum += rinten * inseg ;
            //
            outseg -= inseg ;
            inseg = 1;
            //copy pixel values for next use;
            bval = bnext ;
            gval = gnext ;
            rval = rnext ;
            if (u < src_len) {   //dest를 채우기에 부족한 부분은 마지막에서 계속 반복 사용;
                bnext = src[0];  //여기서 끝낼수도 있음;
                gnext = src[1];
                rnext = src[2];
            }
            src += src_stride ;
            u++ ;
        } else {
            bsum += binten * outseg ; bsum /= isf ;
            gsum += ginten * outseg ; gsum /= isf ;
            rsum += rinten * outseg ; rsum /= isf ;
            dest[0] = (BYTE)min(max(bsum, 0), 255);
            dest[1] = (BYTE)min(max(gsum, 0), 255);
            dest[2] = (BYTE)min(max(rsum, 0), 255);
            dest += dest_stride ;
            //
            bsum = gsum = rsum = 0;
            inseg -= outseg ;
            x++ ;
            if (x == dest_len) break ;
            isf = f[x + 1] - f[x] ;
            outseg = isf ;
        }
    }
}

kipl.tistory.com/41

 

Fant's Resampling

// 배열첨자(dj)와 픽셀의 실제위치(srcindex, dstindex)를 따로 분리하여서 // 열방향을 따라서 작업을 하더라도 메모리 복사가 필요없이 처리가 가능하도록 하였음. BOOL resampleRGB(BYTE *src, BYTE* dst, in..

kipl.tistory.com

 

728x90

'Image Recognition' 카테고리의 다른 글

Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Affine Transformation  (0) 2010.01.20
Color Counting  (0) 2010.01.18
Isometric Transformation  (0) 2010.01.11
Posted by helloktk
,

// 배열첨자(dj)와 픽셀의 실제위치(srcindex, dstindex)를 따로 분리하여서
// 열방향을 따라서 작업을 하더라도 메모리 복사가 필요없이 처리가 가능하도록 하였음.
BOOL 
resampleRGB(BYTE *src, BYTE* dst, int length,      // 
               int stride,                         // RGB = 3(행방향), 행의 총 BYTE수(열방향) 
               double outpos[length + 1])          // forward map: src -> dst; 
{
    double *inpos = new double[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.0;            //첫 입력 픽셀은 완전히 사용가능;
    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++ ;
        }
    }
    delete [] inpos;
    return TRUE;
}

728x90
Posted by helloktk
,