이미지의 히스토그램을 이용하여서 전경과 배경을 분리하는 이진화 과정은 가우시안 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) {
더보기
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 정도임.
