728x90

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

Mixing parameter가 πa (a=1, 2,..., nclass)일 때 특정 픽셀 (=xi)이 클래스 a 소속일 posterior 확률은

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

   

 

 

log-likelihood:

.


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

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

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

더보기

    double a = 1 / sqrt(2*M_PI * var);
    double b = 0.5*(x-mean)*(x-mean)/var;
    return a * exp(-b);
};

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

더보기
    double marginal = 0;
    for (int c = 0; c < nclass; c++) {
        marginal += mclass[c].prob * gauss1d(x, mclass[c].mean, mclass[c].var) ;
    };
    // Bayes 공식 = prior * PDF;
    return mclass[cid].prob * gauss1d(x, mclass[cid].mean, mclass[cid].var) / marginal;
}
// posterior (class_prob[i][c]) table 만들기;
void update_class_prob(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
더보기
        for (int i = 0; i < nbins; i++) {
            for (int c = 0; c < nclass; c++) {
                class_prob[i][c] = classprob(double(i), nclass, mclass, c);
            }
        }
};
// 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) {
더보기
        double ntot = 0;
        for (int i = 0; i < nbins; i++) ntot += hist[i];
        for (int c = 0; c < nclass; c++) {
            double s = 0;
            for (int i = 0; i < nbins; i++) s += hist[i] * class_prob[i][c];
            mclass[c].prob = s / ntot;
        }
};
// mu[c]; 클래스의 평균;
void update_mean(int nbins, double * hist, int nclass, mixclass* mclass,  double ** class_prob) {
더보기
        double ntot = 0;
        for(int i=0; i<nbins; i++) ntot += hist[i];
        for (int c = 0; c < nclass; c++) {
            double sx = 0.0;
            for (int i = 0; i < nbins; i++) sx += hist[i] * i * class_prob[i][c];
            mclass[c].mean = sx / (ntot * mclass[c].prob);
        }
};
// var[c]; 클래스의 분산;
void update_var(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
더보기
    double ntot = 0;
    for (int i = 0; i < nbins; i++) ntot += hist[i];
    for (int c = 0; c < nclass; c++) {
        double m= mclass[c].mean ;
        double sxx = 0;
        for (int i = 0; i < nbins; i++) sxx += hist[i] * SQR(i - m) * class_prob[i][c];
        mclass[c].var = sxx / (ntot * mclass[c].prob);
    }
};
// M-step; 
void update_parameters(int nbins, double * hist, int nclass, mixclass* mclass, double ** class_prob) {
더보기
    // mixture 파라미터를 갱신;
    update_prob(nbins, hist, nclass, mclass, class_prob);
    // 각 클래스의 평균을 갱신;
    update_mean(nbins, hist, nclass, mclass, class_prob);
    // 각 클래스의 분산을 갱신;
    update_var(nbins, hist, nclass, mclass, class_prob);
};
// initialization;
void init_em(int nbins, double * hist, int nclass, mixclass* mclass) {
더보기
        srand(unsigned(time(0)));
        double mean1 = 0, var1 = 0, ntot = 0;
        int i, k;
        for (k = 0; k < nbins; k++) ntot += hist[k];
        for (i = 0; i < nbins; i++) mean1 += hist[i] * i;
        mean1 /= ntot;
        for (i = 0; i < nbins; i++) var1 += hist[i] * SQR(i - mean1);
        var1 /= ntot;
        for (int c = 0; c < nclass; c++) {
            mclass[c].prob = 1.0 / nclass;          //same mixing parameter;
            mclass[c].mean = rand() % nbins; // random mean;
            mclass[c].var = var1;                     // same standard deviation;
        }
};
// calculate log-likelihood;
double mixLLK(int nclass, mixclass* mclass) { 
더보기
    double llk = 0;
    for (int i = 0; i < nbins; i++) {
        double s = 0 ;
        for (int c = 0; c < nclass; c++) 
            s += mclass[c].prob * gauss1d(double(i), mclass[c].mean, mclass[c].var);
        llk+= log(s);
    }
    return llk;
};
// check termination condition;
bool check_tol(double llk, double llk_p, double  eps) {
더보기
    return (fabs(llk - llk_p) / fabs(llk)) > 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 정도임.

  • Otsu Threshold source code: kipl.tistory.com/17

 

사용자 삽입 이미지

 

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

KMeans Algorithm  (0) 2008.07.19
Robust Line Fitting  (0) 2008.07.08
EM: Binarization  (0) 2008.07.01
EM Algorithm: Line Fitting 예  (0) 2008.06.29
Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
Posted by helloktk

댓글을 달아 주세요