Geometry & Recognition

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 함이 극명하게 보인다.

Histogram equalization(HE)은 주어진 이미지의 픽셀 분포가 모든 픽셀 값에서 같은 확률로 나타나도록 픽셀 값을 변환하여 이미지를 보다 잘 인식되게 만드는 영상 처리 과정이다. 이러한 픽셀 값의 균일한 분포는 엔트로피의 관점에서 보면 최대 엔트로피를 주는 분포이기도 하다. 그러나, HE는 원본 이미지의 밝기를 유지하지 않는다. 따라서, 원본 이미지의 밝기를 유지하면서 엔트로피를 최대화시키는 히스토그램 분포를 찾은 후 그것으로 변환(histogram specification)을 시도하는 방법을 고려하는 것이 보다 현실적일 수 있다. 그러한 목표 히스토그램을 $f(s)$ (연속적인 확률 밀도 함수(pdf)로 취급)라고 하면,

1. 정규화된  조건 (픽셀 값을 $[0,255] \rightarrow [0,1]$ 로 rescale 함)
   $$f(s) \ge 0, \quad \int_0^1 f(s)ds = 1.$$
2. 밝기 보존: 
   $$\int_0^1 sf(s) ds = \mu=\text{pixel mean of input image}.$$

이러한 제약 조건에서 아래의 최대 엔트로피 조건을 만족시켜야 한다.

$$\underset{f}{\text{argmax}}\int_0^1 \Big[ -f(s) \ln f(s)\Big] ds.$$

위 조건를 만족시키는 $f(s)$는 Lagrange-multiplier를 이용하고, 변분 방법을 쓰면 쉽게 구해진다: 다음의 범함수 $J [f]$

$$J=-\int_0^1  f(s) \ln f(s) ds + \lambda_1  \left( \int_0^1 f(s) ds-1\right) + \lambda \left(\int_0^1 s f(s) ds -\mu\right)$$

의 극값을 구하면 된다.

$$\frac{\delta J}{\delta f}=0$$

에서

$$ -\log f - 1 + \lambda_1 f +\lambda  s = 0$$

이를 정리하면,

$$f(s)=e^{\lambda s } e^{\lambda_1-1}$$

을 얻고, 제약 조건을 적용하면 ($ e^{\lambda_1 -1}=\lambda/(e^{\lambda}-1))$

$$f(s) = \left\{\begin{array}{ll} 1,& \mu = 0.5 \\ \frac{\lambda e^{\lambda s}}{e^\lambda -1},& \mu \ne 0.5 \end{array} \right.$$

그리고 $λ$는 

$$ \mu = \frac{\lambda e^\lambda - e^\lambda +1}{\lambda (e^\lambda -1)},\quad(\mu \ne 0.5)$$

의 근이다. 오른쪽은 $λ$에 대해서 단조함수이므로 근이 하나만 존재한다.  원본 이미지에서 픽셀 평균을 계산하여 위 방정식에 대입하면 $λ$를 구할 수 있고, 목표 히스토그램을 얻을 수 있다. 목표 히스토그램의 cdf와 ($\text {chist}(x)=\int_0^x sf(s) ds$), 원본 이미지의 cdf를 구해서, 그 차이를 최소로 하는 픽셀 값의 대응관계를 찾으면 된다 (histogram specification):

$$x \rightarrow y = F(x) = \int_0^x f(s) ds$$

참고: Bright Preserving Histogram Equalization with Maximum Entropy: A Variational Perspective.
        C. Wang and Z.Ye, IEEE Trans. Consumer Electronics. V51. No4. (2005);
// 알고리즘 적용 단계에서 픽셀 값의 범위를 [0,255] -> [0,1]로 변환해야 한다.
// BPHEME의 cumulative histogram(continuous version)::integral of f(s) over s;

double cdf(double s, double mu, double lambda) {
    if (mu == 0.5) return s ;
    else return (exp(lambda * s) - 1) / (exp(lambda) - 1);
}

std::vector<int> hist_spec(const std::vector<double> &chist1,
                           const std::vector<double> &chist2) 
{
    std::vector<int> lut(chist1.size());                           
    for (int i = 0; i < chist1.size(); i++) {
        int j = chist1.size() - 1;
        while (chist2[j] > chist1[i]) j--;
        lut[i] = j < 0 ? 0 : j
    }
    return lut;
}

std::vector<BYTE> BPHEME(const std::vector<BYTE> &src) {
    std::vector<double> hist = make_hist(src);
    normalize_hist(hist);
    // cumulative histogram;
    std::vector<double> chist = make_cumulative_hist0(hist);
    // gray-mean; note, pixel range should be changed [0,255] -> [0,1]
    double mu = hist_mean(hist);
    // determine lambda;
    double lambda, th = 1.e-15;
    FindRoot(mu, th, lambda);
    // dst-cdf;
    std::vector<double> chist2(256);
    for (int i = 0; i < 256; i++) 
        chist2[i] = cdf(double(i) / 255., mu, lambda);
    // histogram-specification;
    std::vector<int> lut = hist_spec(chist, chist2);
    // make dst image;
    std::vector<BYTE> dst(src.size());
    for (int i = src.size(); i-->0;) dst[i] = lut[src[i]];
    return dst;
};

// Root finding procedure ::

double F(double s, double mu) {
    if (s == 0.) return (mu - 0.5);
    double ev = exp(s);
    return mu - (s * ev - ev + 1.)/(s * (ev - 1.)) ;
};
// derivative of F(s, mu_src) w.r.t. s;
double DF(double s, double mu) {   
    if (s == 0.) return -1.0 / 12.0;
    double ev = exp(s), ss = s * s;
    return -(ev * (ev - ss - 2.) + 1.) / (ss * (ev - 1.) * (ev - 1.));
}
// Find root of F(s,mu_src) based on Newton-Raphson method.
int FindRoot(double mu_src, double threshold, double& s) {
    int iter = 0, max_iter = 1000;
    s = 0. ; //initial guess; problem specific!
    for (iter = 0; iter < max_iter; ++iter) {
        double sold = s;
        // ASSERT(DF)
        s = s - F(s, mu_src ) / DF(s, mu_src) ;
        if (fabs(s - sold) < threshold)
            break;
    } ;
    TRACE("terminating iters=%d\n", iter);
    return (iter < max_iter) ;
};

histogram equaliztion 결과;

BPHEME  결과: