Geometry & Recognition

이미지의 히스토그램을 이용하여 전경과 배경을 분리하는 이진화는 가우시안 mixture model과 EM 알고리즘을 적용하기에 좋은 예다. 히스토그램에는 전경에 해당하는 픽셀 분포와 배경에 해당하는 픽셀 분포가 혼합되어 있다. 이를 두 가우시안의 혼합으로 모델링하고 EM 알고리즘을 사용해서 mixing parameter(π_a), 각 클래스의 평균(μ_a)과 표준편차(σ_a)를 추정한다. $N$개의 Gaussian mixture일 때, 

$$\text{Mixture Model:}~~h(x) = \sum_a \pi_a g_a(x)~~~~\Big(\sum _a \pi_a=1\Big)$$

$$\text{Gaussian Mixture:}~~g_a(x) = \phi( x| \theta_a)~~~~~\theta_a =N(\mu_a , \sigma_a^2) $$

Mixing parameter가 π_a $(a=1, 2,..., nclass)$일 때 특정 픽셀 값(=$x_i$)이 클래스 $a$ 소속일 posterior는

$$\text{Posterior:}~~~\gamma_{i,a} \equiv Pr( Z_i =a| x_i, \Theta) = \frac{\pi_a \phi (x_i | \theta_a)}{ \sum_b \pi_b \phi(x_i | \theta_b) }$$

로 쓸 수 있다. posterior 정보를 이용하면 mixing parameter, 평균 그리고 분산은 다음 식으로 주어진다. $H[i] = H_i$는 이미지의 히스토그램을 나타내고, bin 인덱스 $i$는 픽셀 값 $x_i$를 나타낸다. 그러면

$$\pi_a = \frac{\sum_i H_i \gamma_{i, a}}{\sum_k H_k }$$

$$\mu_a = \frac{\sum_i x_i H_i \gamma_{i, a} }{\sum_k H_k \gamma_{k,a}}$$

$$\sigma_a^2 = \frac{\sum_i (x_i - \mu_a) ^2 H_i \gamma_{i,a}}{\sum_k H_k \gamma_{k,a}}$$

log-likelihood:

$$\log (L) = \sum _i \log \Big[\sum_a \pi_a \phi( x_i | \theta_a)  \Big]$$

// mixing 클래스를 기술하는 클래스;
struct mixclass {
    double prob ;               // mixing parameter;
    double mean ;               // mean
    double var ;                // variance;
};

double gauss1d(double x, double mean, double var)

// posterior; Pr(Zi = c | xi, Theta);
// 주어진 관측값 x이 클래스 cid에 속할 posterior;
double classprob(double x, int nclass, mixclass*  mclass, int cid)

// 입력은 이미지의 히스토그램;
double em(int nbins/*=256*/, double hist[/*256*/],
    int nclass/*=2*/, mixclass mclass[/*=2*/], double eps/*=1.e-10*/) {
    double llk = 0, prev_llk = 0;
    // allocate memory buffers for the posterior information;
    double ** class_prob = (double**)malloc(sizeof(double*) * nbins);
    class_prob[0] = (double*)malloc(sizeof(double) * nbins * nclass) ;
    for (int i = 1; i < nbins; i++) class_prob[i] = class_prob[i - 1] + nclass;

    // initialization of algorithm;
    init_em(nbins, hist, nclass, mclass);
    //
    do {
        prev_llk = llk;
        // E-step ;
        update_class_prob(nbins, hist, nclass, mclass, class_prob);
        // M-step;
        update_parameters(nbins, hist, nclass, mclass, class_prob);
        llk = mixLLK(nclass, mclass);
        // TRACE("mean1=%f, mean2=%f\n", mclass[0].mean, mclass[1].mean);
        TRACE("log-likelihood=%e\n", llk);
    } while (!check_tol(llk, prev_llk, eps));
    // clean ;
    free(class_prob[0]);
    free(class_prob) ;
    return llk;
};

• 적색 : 히스토그램 
• 청색, 녹색 : posterior(membership); 
• Otsu 알고리즘을 쓰는 경우에 100에서 threshold 값이 결정되고 EM은 110 정도임.
• Otsu Threshold source code: kipl.tistory.com/17

아래 그림을 보면 우리는 데이터를 연결하는 두 개의 직선을 생각할 수 있을 것이다.  그럼 두 개의 직선을 어떻게 얻을 것인가? 물론, ICA(independent component analysis)를 이용하는 것이 한 가지 방법이 될 것이다. 여기서는 EM 알고리즘을 이용하여서 두 개의 직선을 기술하는 기울기와 $y$-절편의 값을 구하는 방법을 알아보자.

두 직선을 각각 $y=a_1 x + b_1$, $y = a_2  x + b_2$라고 하면, $(a_1, b_1)$, $(a_2, b_2)$를 구하는 문제이다. 만약 각각의 data가 에러가 수반되는 측정에 의해서 얻어졌다고 하자. 에러 분포가 정규분포를 따른다면(여기서는 두 개의 모델 모두 갖은 표준편차를 갖는다고 가정했다) $i$-번째의 데이터가 각각의 직선 모델 1과 2에 속할 확률은 (posterior with equal prior) Bayes 공식에 의해서 

$$ w_1[i] = \frac{ e^{ - \frac{ r_1^2[i]}{2\sigma^2} } } {e^{ - \frac{ r_1^2[i]}{2\sigma^2} } + e^{ - \frac{ r_2^2[i]}{2\sigma^2} } }, \quad w_1[i] = \frac{ e^{ - \frac{ r_2^2[i]}{2\sigma^2} } } {e^{ - \frac{ r_1^2[i]}{2\sigma^2} } + e^{ - \frac{ r_2^2[i]}{2\sigma^2} } }, \quad i = 0,1,2,... $$

로 주어진다. 여기서 $r_1(i)$와 $r_2(i)$는 residual error이다:

$$r_1[i] = a_1 x[i] + b_1 - y[i],\quad r_2[i] = a_2 x[i] + b_2 - y[i], \quad i=0,1,2,...$$

(*이 값 대신에 직선까지 거리=$\frac{|a_k x + b_k - y|}{\sqrt{1+ a_k^2}},~ k=1,2)$로 대체해도 된다)

이제 각각의 데이터에 대해서 posterior를 구하였으므로(E-step) 이 값을 가중치로 하여서 직선의 방정식을 다시 갱신한다. 즉, 각각의 data 점들에 대한 $w_1(i)$를 가중치로 하여서 다시 직선 모델 1의 파라미터를 재계산하고, 동일한 과정을 $w_2(i)$를 가지고 직선 모델 2를 계산하여서 $(a_1, b_1)$, $(a_2, b_2)$를 재계산한다(M-step). 이 갱신된 파라미터를 이용하여서 다시 가중치를 구하는 작업과, 직선의 파라미터를 구하는 작업을 일정하게 수렴할 때까지 반복을 하는 과정을 수행한다.

아래의 그림은 3번 만에 원하는 결과를 얻는 것을 보여준다. 직선의 파라미터에 대한 초기값과 residual error의 표준편차 파라미터에 대한 적절한 값의 선택이 중요하다.

// 코드의 일부...
std::vector<CPoint> data ;                             // data,
std::vector<double> w1(data.size()), w2(data.size());  // weights.
double a1, b1, a2, b2 ;                                // line params;
double sigma ;
// E-step;
void calcWeights() {
    for (int i = data.size(); i-- > 0;) {
        double  x = data[i].x, y = data[i].y ;
        double r1 = a1 * x + b1 - y ;
        double r2 = a2 * x + b2 - y ;
        double n1 = SQR(r1) / SQR(sigma) / 2;
        double n2 = SQR(r2) / SQR(sigma) / 2;
        double p1 = exp( - n1);
        double p2 = exp( - n2);
        w1[i] = p1 / (p1 + p2);
        w2[i] = p2 / (p1 + p2);
    }
};
//  M-step
void estimModels() {
    double s1xx = 0, s1x = 0, s1 = 0, s1y = 0, s1xy = 0;
    double s2xx = 0, s2x = 0, s2 = 0, s2y = 0, s2xy = 0;
    for (int i = data.size(); i-- > 0;) {
        double  x = data[i].x,
                y = data[i].y;
            s1xx += w1[i] * SQR(x);
            s1xy += w1[i] * x * y;
            s1x  += w1[i] * x;
            s1y  += w1[i] * y;
            s1   += w1[i];
            //
            s2xx += w2[i] * SQR(x);
            s2xy += w2[i] * x * y;
            s2x  += w2[i] * x;
            s2y  += w2[i] * y;
            s2   += w2[i];
    };
    double det1 = s1xx * s1 - SQR(s1x);
    double det2 = s2xx * s2 - SQR(s2x);
    a1 = (s1 * s1xy  - s1x * s1y ) / det1;
    b1 = (s1xx * s1y - s1x * s1xy) / det1;
    a2 = (s2 * s2xy  - s2x * s2y ) / det2;
    b2 = (s2xx * s2y - s2x * s2xy) / det2;
}

Gaussian Mixture Model(GMM)은 주어진 데이터가 하나의 가우시안 분포가 아니라 여러 개의 가우시안 분포가 혼합되어 생성되었다고 가정하는 확률 모델이다. 각 가우시안은 하나의 데이터 군집(cluster)을 나타내며, 전체 확률밀도함수는 여러 개의 가우시안 분포 $\mathcal N(y\, ;\mu,\Sigma)$의 가중합으로 표현된다.

$$p(y)=\sum_{j=1}^{k}\alpha_j  \mathcal N(y\, ;\mu_j,\Sigma_j)$$

여기서 \(k\)는 가우시안의 개수, \(\alpha_j\)는 각 가우시안의 혼합비(mixture weight), \(\mu_j\)는 평균, \(\Sigma_j\)는 공분산 행렬을 나타낸다. 혼합비는 \(\sum_j\alpha_j=1\)을 만족한다.

GMM은 각 데이터가 하나의 군집에 완전히 속한다고 가정하는 K-means와 달리, 각 데이터가 여러 가우시안에 속할 확률을 계산하는 확률적(soft) 군집화 방법이다. 따라서 서로 겹치는 군집도 자연스럽게 표현할 수 있으며, 군집의 크기나 방향, 분산이 서로 다른 경우에도 효과적으로 모델링할 수 있다.

GMM의 모수(평균, 공분산, 혼합비)는 일반적으로 EM 알고리즘을 이용하여 추정한다. E-step에서는 각 데이터가 각 가우시안에 속할 확률을 계산하고, M-step에서는 이 확률을 이용하여 평균, 공분산, 혼합비를 다시 추정한다. 이 두 과정을 반복하면 likelihood가 증가하는 방향으로 모수가 갱신되어 최종 모델이 얻어진다.

영상 처리에서는 컬러 분할(color segmentation), 배경 제거(background subtraction), 객체 추적, 영상 압축 등 다양한 분야에서 GMM이 널리 사용된다. 특히 배경 제거에서는 각 픽셀의 시간에 따른 밝기 변화를 여러 개의 가우시안으로 모델링하여, 반복적으로 나타나는 배경은 높은 확률의 가우시안으로 표현하고, 새롭게 등장한 물체는 전경으로 검출하는 데 많이 활용된다.

다음은 이 과정을 수행하는 C++ class와 사용 예제, 및 결과를 보여준다.

//가우시안 kernel Class;
class GaussKernel2D {
    double mx, my;
    double sdx, sdy, sdxy;
    double weight;
public:
    std::vector<double> posterior;
    std::vector<double> wgauss;

    void init(std::vector<POINT> &pts, int sx, int sy, int w) {
        posterior.resize(pts.size());
        wgauss.resize(pts.size());
        weight = w ;
        // initialize model parameters(random하게 선택함-->try another method!!::
        // 주어진 데이터로 전체 분포의 범위와 중심을 알 수 있고, 이를 주어진 클래스 수만큼 임의로 
        // 분배하는 방식으로 초기조건을 설정하면 보다 안정적으로 동작한다)
        mx = sx * double(rand()) / RAND_MAX;
        my = sy * double(rand()) / RAND_MAX;
        sdx = sx / 4 + 1;
        sdy = sy / 4 + 1;
        mx += rand() % 100 ;
        my += rand() % 100 ;
        sdxy = 0;
    }
    double gauss2d(double x, double y) { 
        double varx = sdx * sdx;
        double vary = sdy * sdy;
        double det = varx * vary - sdxy * sdxy;
        double idet = 1.0 / det ;
        double dxx = (x - mx) * (x - mx);
        double dyy = (y - my) * (y - my);
        double dxy = (x - mx) * (y - my);
        return (1.0 / sqrt(det) / 6.28319) * exp(-0.5 * (dxx * vary \
                               + dyy * varx - 2. * dxy * dxy) * idet); 
    }
    void getParams(std::vector<POINT> &pts, double prior = 0) {
        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
        for (int j = pts.size(); j-->0;) {
            double x = pts[j].x, y = pts[j].y;
            sx  += posterior[j] * x;
            sy  += posterior[j] * y;
            sxx += posterior[j] * x * x;
            syy += posterior[j] * y * y;
            sxy += posterior[j] * x * y;
        }
        double denorm = np * weight;
        mx = sx/denorm; my= sy / denorm;
        double devx = sxx / denorm - mx * mx ;
        if (devx <= 0) devx = 0.001;
        sdx = sqrt(devx);
        double devy=syy / denorm - my * my;
        if (devy <= 0) devy = 0.001;
        sdy = sqrt(devy);
        sdxy = sxy / denorm - mx * my;
        // if prior = non-zero -> weight = weight*(1-alpha)+alpha*prior; alpha=0.1?
    };
    // weight; // posterior;
    void estimate(std::vector<double> &px) {
        weight = 0;
        for (int j = px.size(); j-->0;) {
            posterior[j] = wgauss[j] / px[j];    
            weight += posterior[j];
        }
        weight /= px.size(); 
    } 
    //P(x|thetal) * prior;
    void setProb(std::vector<POINT> &pts) {
        for (int i = pts.size(); i-->0;)
            wgauss[i] = weight * gauss2d(pts[i].x, pts[i].y);
    }
    void Draw(CDC* pDC, DWORD color = RGB(0xFF, 0, 0)) {
        CPen pen0(PS_SOLID, 1, color);
        CPen* pOld = pDC->SelectObject(&pen0);
        drawCon(pDC, mx, my, sdx, sdy, sdxy);        // draw ellipses;
        pDC->Ellipse(mx - 2, my - 2, mx + 2, my + 2);// draw center of ellipse;
        pDC->SelectObject(pOld);            
    }
};

void em_main(std::vector<POINT> &pts) {
    GaussKernel2D kernel[NKERNELS];
    int nclass = NKERNELS;
    double weights[20] = {1};
    std::vector<double> px(pts.size());
    double wsum = 0 ;

    for (int i = 0; i < nclass; i++) {
        kernel[i].init(pts, 400, 400, weights[i]);
        wsum += weights[i];
    };    
#define MAX_ITER 50
    for (int iter = 0; iter < MAX_ITER; ++iter) {   
        for (i = px.size(); i-->0;) px[i] = 0;
        for (int k = 0; k < nclass; k++){
            GaussKernel2D &gker = kernel[k];
            gker.setProb(pts); 
            for (int i = px.size(); i-->0;)
                px[i] += gker.wgauss[i] ;
        }        
        for (k = 0; k < nclass; k++) {
            kernel[k].estimate(px);
            kernel[k].getParams(pts);
        }
        //또는 log-likelihood를 계산하여서 그 변화가 적으면 loop-끝내면 된다..
    }
}

//참고 : 아래의 데이터는 사전에 라벨링이 된 것이 아니다. 컬러링은 한번 계산한 후에 분포에 맞게 컬러를 조절하여서 다시 계산한 것이다.

 

f(y|θ);