EM: Binarization

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

$$\text{Mixture Model:}h(x)=\sum _a{\pi }_a{g}_a(x)(\sum _a{\pi }_a=1)$$

$$\text{Gaussian Mixture:}{g}_a(x)=\varphi (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,\mathrm{\Theta })=\frac{{\pi }_a\varphi (x_i|{\theta }_a)}{\sum _b{\pi }_b\varphi (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 [\sum _a{\pi }_a\varphi (x_i|{\theta }_a)]$$

// 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