Mean Shift Filter

Mean-shift는 확률 밀도 함수의 극대점(mode)을 연속적으로 찾는 방법이다. 주어진 데이터(특징점)가 샘플링된 확률 밀도 함수의 평균과 분산을 알 수 없을 때 적용할 수 있는 non-parametric 알고리즘이다. Clustering, tracking, 또는 filtering 등에 사용될 수 있다.
Mean-shift filter는 각 픽셀을 중심으로 하는 윈도 내의 모든 픽셀에 대해서 특징점 공간(= 픽셀 좌표 공간 + 픽셀 컬러 공간)에서 밀도가 극대인 특징점(특징점 공간에서 질량중심)을 찾고, 이 특징점에 해당하는 픽셀의 컬러로 윈도 중심 픽셀의 컬러를 대체한다. 

아래 예제는 컬러 공간을 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 함이 극명하게 보인다.