Geometry & Recognition

컬러 영상을 처리할 때 가장 흔히 사용하는 컬러 표현은 RGB 컬러이다. 이것은 R, G, B에 각각 8-비트를 할당하여 256-단계를 표현할 수 있게 하여, 전체적으로 256x256x256=16777216가지의 컬러를 표현할 수 있게 하는 것이다. 그러나 인간의 눈은 이렇게 많은 컬러를 다 구별할 수 없으므로 24-비트 RGB 컬러를 사용하는 경우는 대부분의 경우에 메모리의 낭비와 연산에 오버헤드를 가져오는 경우가 많이 생긴다. RGB 컬러 영상을 R, G, B를 각각 한축으로 하는 3차원의 컬러 공간에서의 벡터(점)로 표현이 가능하다. 컬러 영상의 픽셀들이 RGB삼차원의 공간에 골고루 펴져 있는 경우는 매우 드물고, 많은 경우에는 이 컬러 공간에서 군집(groups)을 이루면서 분포하게 된다. 하나의 군(group)은 유사한 RGB 값을 갖는 픽셀들로 구성이 되므로, 이 군에 포함이 되는 픽셀들에게 (군의 중앙에 해당하는) 대표적인 컬러 값을 대체하면 그 군에 포함이 된 픽셀은 이젠 RGB공간에서 한 점으로 표현이 되고, RGB공간상에는 픽셀 수만큼의 점이 있는 것이 아니라, 대표 RGB 값에 해당하는 점만이 존재하게 된다. 따라서 적당한 Lookup테이블(colormap)을 이용하면, 적은 메모리 공간만을 가지고도 원본의 컬러와 유사한 컬러를 구현할 수 있다.

이제 문제는 원래의 컬러 공간을 어떻게 군집화 하는가에 대한 것으로 바뀌었다. 간단한 방법으로는 원래의 컬러 영상이 차지는 하는 RGB공간에서의 영역을 감싸는 최소의 박스를 찾고, 이것을 원하는 최종적으로 원하는 컬러수만큼의 박스로 분할을 하는 것이다. 그러나 박스를 어떨게 분할을 해야만 제대로 컬러를 나누고, 또한 효율적으로 할 수 있는가를 고려해야 한다. 분할된 박스의 크기가 너무 크면 제대로 된 대푯값을 부여하기가 힘들어지고, 너무 작게 만들면 원하는 수에 맞추기가 어렵다.

Median-Cut 양자화(quantization)에서 사용하는 방법은 박스의 가장 긴축을 기준으로 그 축으로  projection 된 컬러 히스토그램의 메디안 값을 기준으로 분할을 하여서 근사적으로 픽셀들을 절반 정도 되게 분리를 한다 (한축에 대한 메디안이므로 정확히 반으로 분리되지 않는다). 두 개의 박스가 이렇게 해서 새로 생기는데, 다시 가장 많은 픽셀을 포함하는 박스를 다시 위의 과정을 통해서 분할을 하게 된다. 이렇게 원하는 수의 박스를 얻을 때까지 위의 과정을 반복적으로 시행을 하게 된다.

여기서 원래의 컬러 값을 모두 이용하게 되면 계산에 필요한 히스토그램을 만들기 위해서 너무 많은 메모리를 사용하게 되고 이것이 또한 연산의 오버헤드로 작용하게 되므로 RGB 컬러 비트에서 적당히 하위 비트를 버리면, 초기의 RGB공간에서의 histogram의 크기를 줄일 수 있게 된다.(보통  하위 3-비트를 버려서, 각각 5-비트만 이용하여, 전체 컬러의 개수를 32x32x32= 32768로 줄인다)

이렇게 RGB공간에서의 컬러 분포가 몇 개의 대표적인 컬러(예:박스의 중앙값)로 줄어들면(양자화 과정:: 공간에 연속적인 분포가 몇 개의 점으로 대체됨), 원본 영상의 각각의 픽셀에서의 대체 컬러 값은 원래의 컬러와 가장 유사한, 즉 RGB 공간에서 Euclidean 거리가 가장 작은 박스의 컬러로 대체하면 된다. 

그러나 너무 적은 수의 컬러로 줄이게 되면 인접 픽셀 간의 컬러 값의 차이가 눈에 띄게 나타나는 현상이 생기게 된다. 이러한 것을 줄이기 위해서는 디더링(dithering) 과정과 같은 후처리가 필요하게 된다.

**네이버 블로그에서 이전;

source code: dev.w3.org/Amaya/libjpeg/jquant2.c

int MedCutQuantizer(CRaster& raster, CRaster& out8) {
    const int MAX_CUBES = 256;
    const int HSIZE =32 * 32 * 32; // # of 15-bit colors;
    int hist[HSIZE] = {0};
    CSize sz = raster.GetSize();
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            int b = *p++; int g = *p++;	int r = *p++;
            hist[bgr15(b, g, r)]++;
        }
    }
    std::vector<int> hist_ptr;
    for (int color15 = 0; color15 < HSIZE; color15++)
        if (hist[color15] != 0) { 
            // 0이 아닌 히스토그램의 bin에 해당하는 15비트 컬러를 줌;
            hist_ptr.push_back(color15);
        }
    
    std::vector<Cube> cubeList;
    Cube cube(0, hist_ptr.size()-1, sz.cx * sz.cy, 0);
    shrinkCube(cube, &hist_ptr[0]);
    cubeList.push_back(cube);
    while (cubeList.size() < MAX_CUBES) {
        int level = 255, splitpos = -1; 
        for (int k = cubeList.size(); k-->0;) {
            if (cubeList[k].lower == cubeList[k].upper) ;	// single color; cannot be split
            else if (cubeList[k].level < level) {
                level = cubeList[k].level;
                splitpos = k;
            }
        }
        if (splitpos == -1)	// no more cubes to split
            break;
        // Find longest dimension of this cube
        Cube cube = cubeList[splitpos];
        int dr = cube.rmax - cube.rmin;
        int dg = cube.gmax - cube.gmin;
        int db = cube.bmax - cube.bmin;
        int longdim = 0;
        if (dr >= dg && dr >= db) longdim = 0;
        if (dg >= dr && dg >= db) longdim = 1;
        if (db >= dr && db >= dg) longdim = 2;
        // 가장 넓게 분포한 color를 colorTab의 15비트 컬러의 상위비트 이동;
        reorderColors(&hist_ptr[0], cube.lower, cube.upper, longdim);
        // hist_ptr값을 증가 순서로 바꿈->가장 폭이 넓은 컬러 방향으로 정렬됨;
        quickSort(&hist_ptr[0], cube.lower, cube.upper);
        // hist_ptr을 다시 원컬러로 복구;
        restoreColorOrder(&hist_ptr[0], cube.lower, cube.upper, longdim);
        // Find median
        int count = 0;
        int med = 0;
        for (med = cube.lower; med < cube.upper; med++) {
            if (count >= cube.count/2) break;
            count += hist[hist_ptr[med]];
        }
        // cube에 들어있는 color의 median을 기준으로 두개 cube로 쪼김;
        // 낮은 쪽은 원래 cube를 대체하고 높은 쪽은 새로 list에 추가;
        Cube cubeLow(cube.lower, med - 1, count, cube.level + 1);
        shrinkCube(cubeLow, &hist_ptr[0]);
        cubeList[splitpos] = cubeLow;				// add in old slot
        // 
        Cube cubeHigh(med, cube.upper, cube.count - count, cube.level + 1);
        shrinkCube(cubeHigh, &hist_ptr[0]);
        cubeList.push_back(cubeHigh);	
    }
    // 각 cube에 포함된 컬러를 하나의 대표값으로 설정; 
    // 대표값은 cube내 컬러의 평균이다;
    BYTE rLUT[256], gLUT[256], bLUT[256];
    for (int k = cubeList.size(); k-->0;) {
        int rsum = 0, gsum = 0, bsum = 0;
        Cube &cube = cubeList[k];
        for (int i = cube.lower; i <= cube.upper; i++) {
            int color15 = hist_ptr[i];
            rsum += getRed(color15) * hist[color15];
            gsum += getGreen(color15) * hist[color15];
            bsum += getBlue(color15) * hist[color15];
        }
        rLUT[k] = int(rsum / (double)cube.count + .5);
        gLUT[k] = int(gsum / (double)cube.count + .5);
        bLUT[k] = int(bsum / (double)cube.count + .5);
    }
        
    // 각 cube에 포함된 color를 대표하는 컬러값을 쉽게 찾도록 lookup table
    // 설정함;
    BYTE inverseMap[HSIZE]; 
    for (int k = cubeList.size(); k-->0;) {
        // 각 cube 내의 컬러는 하나의 대표컬러로...;
        Cube& cube = cubeList[k];
        for (int i = cube.lower; i <= cube.upper; i++)
            inverseMap[hist_ptr[i]] = (BYTE)k;
    }	

    // 8-bit indexed color 이미지 생성;
    out8.SetDimensions(sz, 8);
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)out8.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            int b = *p++; int g = *p++;	int r = *p++;
            *q++ = inverseMap[bgr15(b, g, r)];
        }
    }
    // 8비트 이미지 팔레트 설정;
    for (int i = 0; i < 256; i++) {
        out8.m_palette[i].rgbBlue = bLUT[i];
        out8.m_palette[i].rgbGreen = gLUT[i];
        out8.m_palette[i].rgbRed = rLUT[i];
    }
    // SaveRaster(out8, "res_medcut.bmp");
    return 1;
}

Union-Find 알고리즘을 이용하여 영역 분할 과정이다. 각각의 픽셀에서 4방향으로 연결된 픽셀이 속한 영역 merge 할 수 있는지를 시도한다. merge 조건은 현재 픽셀의 그레이 값과 인접 영역의 평균 그레이 값의 차이가 주어진 임계값보다 작은 경우다.

$$\tt \text{merge 조건: }\quad | 그레이 값- \text{인접 영역 평균값}| < threshold$$

컬러 영상의 경우에는 RGB 채널에 조건을 부여하거나 gray level만 이용해서 판단할 수 있다. root node 수만큼의 분할 영역이 만들어지고, 임계값이 클수록 분할된 각 영역의 사이즈가 커진다.  

gray = 0.2989 * r + 0.5870 * g + 0.1140 * b ;

같은 평균 픽셀값을 가지고 있더라도 4-방향으로 서로 연결된 영역이 아니면 합병되지 못하고 서로 다른 영역으로 남는다. 이를 하나의 영역으로 만들려면 분할된 영역을 다시 검사를 하는 과정이 필요하다. 자세한 과정은 Statistical Region Merge(SRM) 알고리즘(kipl.tistory.com/98)에 잘 구현이 되어 있다.

struct Universe {
    int n;
    std::vector<int> sizes;
    std::vector<int> sums;
    std::vector<int> parent;
    Universe(int n, BYTE *image) {
        this->n = n;
        parent.resize(n, -1);  // all nodes are root;
        sizes.resize(n, 1); 
        sums.resize(n);
        for (int i = 0; i < n; i++)
            sums[i] = image[i];
    }
    int FindCompress (int node) { // find root node;
        if ( parent[node] < 0 ) return node;	 // root-node;
        else  return parent[node] = FindCompress (parent[node]);
        // 찾음과 동시에 depth을 줄인다;   
    }
    int Find(int node) { // find root;
        while (parent[node] >= 0) node = parent[node];
        return node;
    }
    bool Predicate(int root1, int root2, int thresh) {
        double diff = double(sums[root2]) / sizes[root2] 
                      - double(sums[root1]) / sizes[root1];
        return fabs(diff) < thresh;
    }
    void Merge(int root1, int root2) {
        sums[root2]  += sums[root1];
        sizes[root2] += sizes[root1];
        parent[root1] = root2;			 // r2 트리로 합병함;
    }
};
// root 노드는 분할영역의 픽셀 갯수와 픽셀 값의 합을 저장한다.
// 처음 각 픽셀이 root 이고, 픽셀 갯수는 1, 픽셀값은 자신의 픽셀값을 저장;
BOOL UnionFindAverage(BYTE *image, int width, int height, int thresh) {
    Universe *UF = new Universe(width * height, image);
    // 4-connectivity:
    // 첫행; LEFT 픽셀이 속한 영역의 평균값과 차이가 thresh 내이면 left로 합병;
    for ( int x = 1; x < width; x++ ) {
        int left = UF->FindCompress (x - 1);
        int curr = UF->FindCompress (x);
        if (UF->Predicate(curr, left, thresh))
            UF->Merge(curr, left);
    }
    //Flatten(y=0);
    for (int x = 0; x < width; x++) UF->FindCompress(x);
    for ( int y = 1, pos = y * width; y < height; y++ ) {
        // x = 0; TOP 픽셀이 속학 영역과 합병 체크;
        int up   = UF->FindCompress (pos - width);
        int curr = UF->FindCompress (pos);
        if (UF->Predicate(curr, up, thresh))
            UF->Merge(curr, up);
        pos++;
        // TOP-LEFT 픽셀 영역과 합병 체크;
        for ( int x = 1; x < width; x++, pos++ ) {
            int left = UF->FindCompress(pos - 1);
            int curr = UF->FindCompress(pos);
            // left와 합병이 되면 left가 up과 합병이 되는지 확인;
            if (UF->Predicate(curr, left, thresh)) {
                UF->Merge(curr, left);
                curr = left;
            }
            int up = UF->FindCompress (pos - width);
            if ((curr != up) && (UF->Predicate(curr, up, thresh))) 
                UF->Merge(curr, up);
        }
        // Flatten(y>0)
        for (int x = 0; x < width; x++) UF->FindCompress(x + y * width);
    }

    // 평균 이미지 생성;
    for ( int k = width*height; k-->0;) {
        int root = UF->Find(k);
        int avg = int(double(UF->sums[root]) / UF->sizes[root] + 0.5);
        image[k] = avg > 255 ? 255 : avg;
    }
    delete UF;
    return TRUE;
};

Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.

이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram($\tt ch$)과 평균을 계산할 수 있는 1차 culmuative histogram($\tt cxh$)을 입력으로 사용한다. 

$$ {\tt thresholds}= \text {argmax} \left( \sigma^2_B = \sum_{j=0}^{nclass-1} \omega_j m_j^2 \right)$$

* Otsu 알고리즘을 이용한 이미지 이진화 코드: kipl.tistory.com/17

* 좋은 결과를 얻으려면 히스토그램에 적당한 필터를 적용해서 smooth하게 만드는 과정이 필요하다.

// 0 <= start < n;
double histo_partition(int nclass, double cxh[], int ch[], int n, int start, int th[]) {
    if (nclass < 1) return 0;
    if (nclass == 1) {
        int ws; double ms;
        if (start == 0) {
            ws = ch[n - 1];
            ms = cxh[n - 1] / ws;
        } else {
            ws = ch[n - 1] - ch[start - 1];             // start부터 끝까지 돗수;
            ms = (cxh[n - 1] - cxh[start - 1]) / ws;    // start부터 끝까지 평균값;
        }
        th[0] = n - 1;
        return ws * ms * ms;                            // weighted mean;
    }

    double gain_max = -1;
    int *tt = new int [nclass - 1];
    for (int j = start; j < n; j++) {
        int wj; double mj;
        if (start == 0) {
            wj = ch[j]; 
            mj = cxh[j];                    //mj = cxh[j] / wj;
        }
        else {
            wj = ch[j] - ch[start - 1];     //start부터 j까지 돗수;
            mj = (cxh[j] - cxh[start - 1]); //mj = (cxh[j] - cxh[start - 1]) / wj;
        }
        if (wj == 0) continue;
        mj /= wj;                           //start부터 j까지 평균;
        double gain = wj * mj * mj + histo_partition(nclass - 1, cxh, ch, n, j + 1, tt);
        if (gain > gain_max) {
            th[0] = j;
            for (int k = nclass - 1; k > 0; k--) th[k] = tt[k - 1];
            gain_max = gain;
        }
    }
    delete [] tt;
    return gain_max;
};

trimodal 분포의 분리;

class0: 0~th[0]

class1: (th[0]+1)~th[1],

class2: (th[1]+1)~th[2]=255

// recursive histo-partition 테스트;
// 0--t[0];--t[1];...;--t[nclass-2];t[nclass-1]=255=n-1;
// nclass일 때 threshold 값은 0---(nclss-2)까지;
double GetThreshValues(int hist[], int n, int nclass, int th[]) {
    if (nclass < 1) nclass = 1;
    // preparing for 0-th and 1-th cumulative histograms;
    int *ch = new int [n];          // cdf;
    double *cxh = new double [n];   //1-th cdf;
    ch[0] = hist[0];
    cxh[0] = 0; // = 0 * hist[0]
    for (int i = 1; i < n; i++) {
        ch[i] = ch[i - 1] + hist[i] ;
        cxh[i] = cxh[i - 1] + i * hist[i];
    }
    // nclass=1인 경우도 histo_partition()내에서 처리할 수 있게 만들었다.
    double var_b = histo_partition(nclass, cxh, ch, n, 0, th);
    delete [] ch;
    delete [] cxh;
    return var_b;
}

사진 속에 포함된 물체를 판별하기 위해서는 보통 이미지를 이진화시켜 binary 이미지를 만들고, conneted component labeling 등을 이용하여 영상에서 겹치지 않는 다수의 물체들을 분리해 내는 과정을 거친다. 이렇게 분리된 각 물체의 외각을 감싸는 convex hull을 찾으려면 먼저 물체의 경계을 추적해서 폴리곤을 만들고, 그 다음 이 폴리곤의 convex hull을 찾는 알고리즘을 적용하면 된다. 이들 폴리곤은 일정한 규칙에 의해서 정렬이 되어 있으므로 convex hull 알고리즘들 중에 chain-hull 이나 Melkman 알고리즘 이용하면 빠르게 결과를 얻을 수 있다.

Raster-scan 방식은 물체의 경계점들을 일정하게 정렬된 순서로 읽게 되므로 on-site convex hull 찾기가 가능하게 된다 (chain-hull에서도 동일한 순서가 들어오지만, 전체 점들을 다 얻은 다음에 다시 분류작업을 해야 하므로, 읽은 즉시 바로 convex  hull을 구성하는 알고리즘은 제공하지 못한다). 아래의  코드는 binary raster image를 스캐닝하면서 경계의 convex hull을 찾는 과정을 보여준다.

#define bg     0x00             /* foreground color (black) */ 
#define fg     0xFF             /* background color (white) */
//cross(AB, AC) > 0 if C is lying on the left side of edge(AB);
#define CCW(A,B,C) ((B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x))

/* find a convex hull of foreground object in a binary image; */ 
/* in some case L[0] = R[0], or L[ll] = R[rr] if first line or last line of object 
** is composed of a single point (수정본:: 배경 체크 부분);
** convex hull이 중복점을 갖는 경우를 허용하지 않을 때는 comment에 따라서 처리해야 한다.
*/ 
int RasterHull(BYTE *image, int w, int h, CPoint H[/* 2xh */]) { 
    CPoint *L = &H[0]; //stack of left-side hull; 
    CPoint *R = &H[h]; //stack of right side hull; 
    int rr = -1, ll = -1, x; 
    CPoint q; 
    for (int y = 0 ; y < h; y++) {
        BYTE *scanline = &image[y*w]
        for (x = 0; x < w; x++)
            if (scanline[x] != bg) break; 
        if (x == w) continue; 
        q = CPoint(x, y);
        while (ll > 0) {
            // q가 이전 edge의 오른쪽이면 convex hull의 꼭지점이 됨;
            if (CCW(L[ll - 1], L[ll], q) < 0) break; 
            else ll--; 
        } 
        L[++ll] = q; 
        for (x = w; x-- > 0;) //x=-1 never occurs; 
            if (scanline[x] != bg) break; 
        // 행에 싱글 픽셀 있는 경우 왼쪽과 오른쪽이 점이 중복이 되는 경우가
        // 생긴다. 이 중복을 없애고 싶으면
        if (x <= q.x) continue ;
        q = CPoint(x, y); 
        while (rr > 0) { 
            // q가 이전 edge의 왼편에 있으면 convex hull의 꼭지점이 됨;
            if (CCW(R[rr-1], R[rr], q) > 0) break; 
            else rr--; 
        } 
        R[++rr] = q; 
    } 
    // 왼쪽 처음에 점과 오른쪽 처음 점유 중복될 가능성이 있다(싱글 픽셀)
    // 왼쪽 마지막 점과 오른쪽 마지막 점이 중복이 되는 경우도 있다(싱글 픽셀) 
    // convex hull에서 중복점을 완전히 제외하고 싶으면 이를 체크하는 부분이 필요;
    // copy right hull to H;
    int n = ll + 1; // # of pts within left hull;
    // remove bottom degeneracy; 
    if (H[ll] == R[rr]) rr--;
    for (int j = rr; j >= 0; j--) H[n++] = R[j];//right hull; 
    //remove top degeneracy;
    if (n > 1 && H[0] == H[n - 1]) n--;
    return n; 
}

위 함수는 픽셀 값이 bg가 아닌 것은 모두 전경으로 처리하도록 작성하였다. 만약 connected component를 구한 상태에서 특정 component의 convex hull을 구하고자 하는 경우에는 if-문의 조건을 image[y * w + x] == index로 바꾸어야 한다 (여기서 index는 connected component label이다. 그리고 labeling 이미지는 일반적으로 정수 이미지이므로 그에 맞게 수정을 해야 한다.) 

data matrix 코드를 포함하는 이미지에서 코드 영역을 찾기 위해서는 우선 이진화를 시켜야 한다. 다음으로 이진 이미지에서 connected component labeling을 해서 사이즈가 작은 blob을 제거하면 거의 코드 영역만 남길 수 있다. 그림은 코드 영역(노란색)을 감싸는 convex hull(붉은색 라인)을 찾은 것을 보여준다.

Convex hull 알고리즘은 아래에서 찾을 수 있다.

• Chain hull: kipl.tistory.com/109 
• Quick hull: kipl.tistory.com/108  
• Brute Force Algorithm: kipl.tistory.com/106

Kuwahara filter는 주어진 픽셀을 한 코너로 하는 4개 block 중에서 분산이 가장 작은 블록의 평균으로 중앙 픽셀 값을 대체하는 비선형 filter다. 이 필터를 적용하면  균일한 영역에서  spike noise 픽셀은 주변의 영역의 픽셀 값으로 대체되므로 제거될 수 있고, edge 라인을 기준으로 수직방향으로는 거의 균일한 영역이므로 필터링을 하더라도 edge가 뭉개지지 않음을 알 수 있다. Kuwahara filter는 이미지의 노이즈를 제거하고 edge를(위치) 보존하는 특성을 가진다. 윈도의 크기가 5x5인 경우 4개의 3x3영역에서 평균과 분산을 계산해서 중앙 픽셀의 값을 결정한다. integral 이미지를 사용하면 block의 평균과 분산을 매번 새롭게 계산할 필요가 없으므로 연산 비용이 많이 감소한다.

Kuwahara filter보다 더 향상된 filter로는 symmetric nearest neighbour (SNN) filter가 있다. 이 filter는 기준 픽셀에서 대칭 지점에 있는 두 지점의 픽셀 값 중에서 중앙 픽셀 값에 가까운 것을 취해서 평균값을 취한다.

아래 코드는 gray 이미지에 대해 Kuwahara filter을 구현한 것으로 컬러 이미지의 경우는 RGB 채널별로 따로 적용하면 된다. 그리고 경계 영역에서는 처리를 하지 않도록 만들었다.

#define GET_SUM(x,y)  sum[ (y) * width + (x)]
#define GET_SUM2(x,y) sum2[(y) * width + (x)]
#define BLK_SUM(x,y) ((GET_SUM(x,y) + GET_SUM(x-hs1,y-hs1) \
                        - GET_SUM(x-hs1,y) - GET_SUM(x,y-hs1)))
#define BLK_SUM2(x,y) ((GET_SUM2(x,y) + GET_SUM2(x-hs1,y-hs1) \
                        - GET_SUM2(x-hs1,y) - GET_SUM2(x,y-hs1)))
void integral_image(BYTE *image, int width, int height, int *sum);
void integral_image_sq(BYTE *image, int width, int height, int *sum2);
void kuwahara_filter ( BYTE *image, int width, int height, int size, BYTE *out ) {
    std::vector<int> sum(width * height, 0);
    stt::vector<int> sum2(width * height, 0);
    // integral image; for window mean;
    integral_image(image, width, height, &sum[0]);
    // integral image of square; for window variance;
    integral_image_sq(image, width, height, &sum2[0]);
    
    // do filtering;
    size = ((size >> 1) << 1) + 1; // make window size be an odd number;
    int offset = ( size - 1 ) >> 1;
    int hs1 = ( size + 1 ) >> 1;
    int nn = hs1 * hs1;    // # of pixels within the  window;
    int var[4], mean[4];
    for ( int y = height - hs1; y-- > hs1;) {
        for ( int x = width - hs1; x-- > hs1; ) {
            mean[0] = BLK_SUM ( x, y ) / nn;
            mean[1] = BLK_SUM ( x + offset, y ) / nn;
            mean[2] = BLK_SUM ( x, y + offset ) / nn;
            mean[3] = BLK_SUM ( x + offset, y + offset ) / nn;
            var[0] = BLK_SUM2 ( x, y ) / nn                   - mean[0] * mean[0];
            var[1] = BLK_SUM2 ( x + offset, y ) / nn          - mean[1] * mean[1] ;
            var[2] = BLK_SUM2 ( x, y + offset ) / nn          - mean[2] * mean[2] ;
            var[3] = BLK_SUM2 ( x + offset, y + offset ) / nn - mean[3] * mean[3] ;
            int vmin = var[0], winner = 0;
            for ( int i = 1; i < 4; ++i )
                if ( var[i] < vmin )
                    vmin = var[i], winner = i;
            int a = mean[winner];
            out[y * width + x] = a > 255 ? 255 : a < 0 ? 0 : a;
        }
    }
};