Mean Shift Filter

Mean-shift는 확률 밀도 함수(probability density function)의 극대점(mode)을 반복적으로 찾아가는 알고리즘이다. 데이터가 어떤 확률분포를 따른다고 가정하지 않고, 샘플링된 데이터만으로 밀도 분포를 추정하는 비모수(non-parametric) 방법이므로 데이터의 평균이나 분산과 같은 분포의 매개변수를 미리 알 필요가 없다. 이러한 특성 때문에 Mean-shift는 clustering, tracking, segmentation, filtering 등 다양한 컴퓨터 비전 문제에 널리 사용된다.

Mean-shift의 기본 아이디어는 현재 위치를 중심으로 일정한 크기의 윈도(kernel)를 설정한 뒤, 윈도 안에 포함된 데이터의 centroid를 계산하고 윈도의 중심을 그 위치로 이동시키는 것이다. 이 과정을 질량중심의 위치가 더 이상 변하지 않을 때까지 반복하면, 최종적으로 데이터의 국소적인 밀도 극대점(mode)에 수렴하게 된다. 따라서 서로 다른 초기 위치에서 시작한 점들은 같은 mode로 수렴하며 하나의 클러스터를 형성한다.

Mean-shift filter는 이러한 원리를 영상 처리에 적용한 필터이다. 각 픽셀을 중심으로 하는 윈도를 설정한 뒤, 윈도 안의 모든 픽셀을 특징점 공간(feature space) 에서 고려한다. 특징점 공간은 일반적으로 픽셀의 위치 \((x,y)\) 와 색상(RGB 또는 Lab 등) 을 함께 사용하는 공간으로 정의된다. 이 공간에서 Mean-shift 알고리즘을 적용하여 국소적인 밀도 극대점, 즉 centroid를 찾은 후, 그 특징점에 대응하는 색상값으로 윈도 중심 픽셀의 색상을 대체한다.

이 과정을 영상의 모든 픽셀에 대해 수행하면 공간적으로 가까우면서도 색상이 유사한 픽셀들은 동일한 mode로 수렴하게 된다. 따라서 객체의 경계는 비교적 잘 유지하면서도 동일한 영역 내부의 잡음은 제거되는 효과를 얻을 수 있다. 즉, Mean-shift filter는 에지를 보존(edge-preserving)하면서 영상을 smoothing 하는 대표적인 비선형 필터로 볼 수 있다. 또한 반복적으로 적용하면 비슷한 색상을 갖는 영역들이 하나의 균일한 영역으로 합쳐져 영상 분할의 전처리 과정으로도 자주 사용된다.

아래 예제는 컬러 공간을 RGB에서 YIQ로 바꾸어서 알고리즘을 적용한 후에 그 결과를 RGB로 다시 변환하였다. 윈도를 크게 하면 알고리즘의 수행 속도가 매우 느려진다. 공간 윈도를 원에서 정사각형으로 바꾸어도 된다.

Clustering 알고리즘으로 볼 때 Mean-Shift 알고리즘의 특징은

1. 단순한 알고리즘 구조를 갖지만 연산 비용이 매우 크다. 고차원 특징에 대해서는 부적합
2. 클러스터 수에 대한 제한이 없다.(k-means와 비교)
3. 클러스터 초기 위치 설정이 필요없다.(k-means와 비교)
4. outliers의 영향이 적다.(k-means와 비교)
5. 유일한 매개변수는 윈도 크기다. 원도 크기를 미리 설정해야 한다.

// 구현 code

void RGB2YIQ(BYTE *color, int w, int h, float *ycomp, float *icomp, float *qcomp);

void RGB2YIQ(BYTE *color, int w, int h, float *ycomp, float *icomp, float *qcomp) {
    for (int i = w * h; i-- > 0;) {
        int b = color[3*i+0], g = color[3*i+1], r = color[3*i+2];
        ycomp[i] = 0.299f  *r + 0.587f *g + 0.114f  *b; // Y
        icomp[i] = 0.5957f *r - 0.2744f*g - 0.3212f *b; // I
        qcomp[i] = 0.2114f *r - 0.5226f*g + 0.3111f *b; // Q
    }
}

// YIQ2RGB color conversion;

void YIQ2RGB(float Yc, float Ic, float Qc, int *b, int *g, int *r) ;

void YIQ2RGB(float Yc, float Ic, float Qc, int *b, int *g, int *r) {
     *r = int(Yc + 0.9563f*Ic + 0.6210f*Qc);
     *g = int(Yc - 0.2721f*Ic - 0.6473f*Qc);
     *b = int(Yc - 1.1070f*Ic + 1.7046f*Qc);
}

#define SQR(x) ((x) * (x))
void MeanShiftFilterRGB(BYTE *image/*color*/, int width, int height,
                    int rad/*=spatial*/, int radcol/*=color space(YIQ)*/,
                    BYTE *out) {
    std::vector<float> ycomp(width*height) ;
    std::vector<float> icomp(width*height) ;
    std::vector<float> qcomp(width*height) ;
    int rad2    = rad * rad ;
    int radcol2 = radcol*radcol;
    // RGB2YIQ;
    RGB2YIQ(image, width, height, &ycomp[0], &icomp[0], &qcomp[0]);
    CRect rect(0, 0, width, height);
    rect.DeflateRect(0, 0, 1, 1);
    for (int y=0, pos=0; y<height; y++) {
        for (int x=0; x<width; x++, pos++) {
            int xc = x, yc = y;
            float Yc = ycomp[pos], Ic = icomp[pos], Qc = qcomp[pos];
            int iters = 0;
            while (1) {
                int oldxc = xc, oldyc = yc;
                float YcOld = Yc, IcOld = Ic, QcOld = Qc;
                float xsum = 0, ysum = 0; 
                float sumY = 0, sumI = 0, sumQ = 0;
                int count = 0;
                CRect rc = CRect(xc-rad, yc-rad, xc+rad, yc+rad) & rect;
                for (int y2 = rc.top; y2 <= rc.bottom; y2++) {
                    int ry = y2 - yc;
                    for (int x2 = rc.left; x2 <= rc.right; x2++) {
                        int rx = x2 - xc;
                        if ((SQR(ry) + SQR(rx)) <= rad2) {    
                            int curpos = y2*width + x2;
                            float Y2 = ycomp[curpos];
                            float I2 = icomp[curpos];
                            float Q2 = qcomp[curpos];
                            float dY = Y2 - Yc, dI = I2 - Ic, dQ = Q2 - Qc;
                            if ((SQR(dY) + SQR(dI) + SQR(dQ)) <= radcol2) {
                                xsum += x2; ysum += y2;
                                sumY += Y2; sumI += I2; sumQ += Q2;
                                count++;
                            }
                        }
                    }
                }
                if (count==0) count = 1;
                //average(YIQ);
                Yc = sumY / count; Ic = sumI / count; Qc = sumQ / count;
                // average(x,y);
                xc = int(xsum / count + 0.5); 
                yc = int(ysum / count + 0.5);
                int dx = xc - oldxc, dy = yc - oldyc;
                float dY = Yc - YcOld, dI = Ic - IcOld, dQ = Qc - QcOld;
                // shift;
                float shift = SQR(dx) + SQR(dy) + SQR(dY) + SQR(dI) + SQR(dQ); 
                if (shift <= 3 || ++iters >= 100) break;
            }
            int r, g, b;
            YIQ2RGB(Yc, Ic, Qc, &b, &g, &r);
            out[3*pos+0] = b > 255 ? 255: b < 0 ? 0: b;
            out[3*pos+1] = g > 255 ? 255: g < 0 ? 0: g;
            out[3*pos+2] = r > 255 ? 255: r < 0 ? 0: r;
        }       
    }
};

(rad=6, radcol=50)인 예.

YIQ space에서 컬러분포:

(rad=6, radcol=20) 인 예.

아래 이미지의 검정색 부분은 unstable point이다. 즉, 이 점들은 항상 다른 점으로 이동한다. 에지점들은 unstable 함이 극명하게 보인다.