'EM Algorithm'에 해당되는 글 3건

  1. 2008.07.01 EM : Binarization
  2. 2008.06.29 EM Algorithm : Line Fitting 예
  3. 2008.06.07 Gaussian Mixture Model (2)

이미지의 히스토그램을 이용하여서 전경과 배경을 분리하는 이진화 과정은 가우시안 mixture model을 이용하여서 EM 알고리즘을 적용하기에 좋은 예다. 전경에 해당하는 픽셀값의 분포와 배경에 해당하는 픽셀값의 분포는 히스토그램상에 섞여서 나타나는데. 이것을 두 가우시안의 혼합으로 생각할 때  em 알고리즘은 각각의 믹싱정도와 각 가우시안의 평균 및 표준편차를 추정한다. mixing parameter πa (a=1, 2,..., nclass)로 표시하는 경우에 특정한 픽셀값에 대한 posterior는


로 쓸수 있다. 여기서 f(x; θ)는 정규분포를 의미한다. 

posterior 정보를 이용하면 mixing parameter πa와 평균, 분산은 다음식으로 갱신이 된다. H[i]는 이미지의 히스토그램을 나타낸다.


   


log-likelihood:

여기서 a는 클래스 인덱스이고, i는 히스토그램 bin 인덱스이다.



// mixing 클래스를 기술하는 클래스;

struct mixclass {
    double prob ;               //mixing parameter;
    double mean ;               // mean
    double var ;                //variance;
} ;
// N(mean, var);

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) {

// posterior (class_prob[i][c]) table 만들기;
void update_class_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// E-step;  pi[c] = mixture parameter for class c;
// posterior를 이용해서 특정클래스의 mixing 정도를 계산;==> next prior;
void update_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// mu[c]; 클래스의 평균;
void update_mean(int nbins, double * hist, int nclass, mixclass* mclass,  double ** class_prob) {
// var[c]; 클래스의 분산;
void update_var(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// M-step; 
void update_parameters(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
// initialization;
void init_em(int nbins, double * hist, int nclass, mixclass* mclass) {
// calculate log-likelihood;
double mixLLK(int nclass, mixclass* mclass) { 
// check termination condition;
bool check_tol(double llk, double llk_p, double  eps) {
// 입력은 이미지의 히스토그램;
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 정도임.

사용자 삽입 이미지

'Image Recognition' 카테고리의 다른 글

Robust Line Fitting  (0) 2008.07.08
Bayesian Spam Filtering  (0) 2008.07.03
EM : Binarization  (0) 2008.07.01
EM Algorithm : Line Fitting 예  (0) 2008.06.29
Shuffling  (0) 2008.06.21
Bayesian Decision Theory  (1) 2008.06.17
Posted by helloktk

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

사용자 삽입 이미지


두 개의 직선을 각각 y=a1 * x + b1 , y = a2 * x + b2,라고 하면, 문제는 (a1, b1), (a2, b2)를 찾는 구하는 문제이다. 만약에 각각의 data 점의 측정에 의한 에러분포가 정규분포를 따른다면(여기서는 두 개의 모델 모두 갖은 표준편차를 갖는다고 가정했다) i-번째의 데이터 점이 각각의 직선모델1 과 2에 속할 확률은 (posterior with equal prior) 베이즈 공식에 의해서 

로 주어진다. 여기서 r1(i)와 r2(i)는 residual error이다:

(*이 값 대신에 직선에서의 거리=로 대체하여도 된다)

이제 각각의 데이터에 대해서 posterior를 구하였으므로(E-step) 이 값을 가중치로 하여서 직선의 방정식을 다시 갱신한다. 즉, 각각의 data 점들에 대한 w1(i)를 가중치로 하여서 다시 직선모델1의 파라미터를 재계산하고, 동일한 과정을 w2(i)를 가지고 직선모델2를 계산하여서 (a1,b1), (a2, b2)를 재계산한다(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=0; i<data.size(); i++){
        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);
        double n2 = SQR(r2) / SQR(sigma);
        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=0; i<data.size(); i++){
        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 ;
}

'Image Recognition' 카테고리의 다른 글

Bayesian Spam Filtering  (0) 2008.07.03
EM : Binarization  (0) 2008.07.01
EM Algorithm : Line Fitting 예  (0) 2008.06.29
Shuffling  (0) 2008.06.21
Bayesian Decision Theory  (1) 2008.06.17
Gaussian Mixture Model  (2) 2008.06.07
Posted by helloktk

N개의 관측값이 주어지고, 이것들이 k개의 가우시안의 mixture로 생각을 할 때, EM알고리즘을 적용하여서 모델을 결정하는 파라미터를 추저하고자 한다. 모델 파라미터는 K개의 가우시안을 기술하는 평균값과 표준편차(Covariance Matrix)와 더블어 mixing정도를 나타내는 파라미터 K개를 추정하여야 한다.

사용자 삽입 이미지


아래는 이 일련의 과정을 수행하는 클래스와 사용예제, 및 결과샘플을 보여준다.

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

    void init(POINT* pts, int np, int sx, int sy, int w) {
        posterior.resize(np);
        wgauss.resize(np);
        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*sdxy)*idet);
    }
    void getParams(POINT *pts, int np, double prior=0)
    {
        double sx=0, sy=0, sxx=0, syy=0, sxy=0;
        for(int j = 0; j < np; j++) {
            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(double *px, int np) {
        weight = 0;
        for(int j = 0; j < np; j++) {
            posterior[j] = wgauss[j] / px[j];   
            weight += posterior[j];
        }
        weight /= np;
    } 

    //P(x|thetal) * prior;
    void setProb(POINT *pts, int np) {
        for(int i=0; i<np; i++){
            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(POINT *pts, int np) {
    GaussKernel2D kernel[NKERNELS];
    int nclass=NKERNELS;
    double weights[20] = {1,1,1,1,1,1};
    std::vector<double> px(np);
    double wsum = 0 ;

    for(int i=0;i<nclass; i++){
        kernel[i].init(pts, np, 400, 400, weights[i]);
        wsum += weights[i];
    };    
#define MAX_ITER 50
    for(int iter = 0; iter < MAX_ITER; ++iter){  
        for(i=0; i<np; i++) px[i]=0;
        for(int k=0; k<nclass; k++){
            GaussKernel2D & gker = kernel[k];
            gker.setProb(pts, np);
            for(int i=0; i<np; i++){
                px[i] += gker.wgauss[i] ;
            }
        }        
        for(k=0; k<nclass; k++) {
            kernel[k].estimate(&px[0], np);
            kernel[k].getParams(pts, np);
        }
        //또는 log-likelihood를 계산하여서 그 변화가 적으면 loop-끝내면 된다..
    }
}
//참고 : 아래의 데이터는 사전에 라벨링이 된것이 아니다. 컬러링은 한번 계산한 후에 분포에 맞게 컬러를 조절하여서 다시 계산한 것이다.

사용자 삽입 이미지

 

f(y|θ);

사용자 삽입 이미지


'Image Recognition' 카테고리의 다른 글

Shuffling  (0) 2008.06.21
Bayesian Decision Theory  (1) 2008.06.17
Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
RANSAC Algorithm  (0) 2008.05.24
Contour Tracing  (0) 2008.05.22
Posted by helloktk