이미지의 히스토그램을 이용하여 전경과 배경을 분리하는 이진화는 가우시안 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;
for (int k = 0; k < nbins; k++) ntot += hist[k];
for (int i = 0; i < nbins; i++) mean1 += hist[i] * i;
mean1 /= ntot;
for (int 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
728x90
'Image Recognition' 카테고리의 다른 글
| KMeans Algorithm (0) | 2008.07.19 |
|---|---|
| Robust Line Fitting (0) | 2008.07.08 |
| EM Algorithm: Line Fitting (0) | 2008.06.29 |
| Gaussian Mixture Model (2) | 2008.06.07 |
| Rasterizing Voronoi Diagram (0) | 2008.05.26 |
