그레이 영상의 히스토그램 $h(x)$를 두 개의 가우시안 분포($g_1(x)$, $g_2(x)$)의 혼합으로 모델링하여 분리하려고 할 때 기준인 decision boundary 값 (threshold value)을 expectation maximization(EM) 알고리즘을 적용하여 구한다. 

 

E-step: compute responsibility of class 2; (for class 1, 1-γ_i)

 

 

M-step: compute the weighted means (μ1, μ2), variances (σ1, σ2) and mixing probability (π)

 

 

log-likelihood: 

$$\log L = \sum _{i} \log \left[ (1- \pi) \phi_{\theta_1 } (x_i) + \pi \phi_{\theta_2 }(x_i) \right] $$

decision boundary 값은 responsibility = 0.5인  bin 인덱스를 선택하면 된다.

아래 그림의 왼쪽은 히스토그램, 오른쪽은 최대우도 gaussian fitting 결과와 왼쪽 분포의 responsibility($1-\gamma_i$)를 그린 것이다.

void estimGaussParams(std::vector<double>& data, int start, int end, double *mean, double *var) ;

더보기
void estimGaussParams(std::vector<double>& data, int start, int end, double *mean, double *var) {
    double s = 0, sx = 0, sxx = 0;
    for (int i = start; i <= end; i++) {
        s += data[i];
        sx += data[i] * i;
        sxx += data[i] * i * i;
    }
    *mean = sx / s;
    *var = (sxx - sx * sx / s) / s;
};

void initGuess(std::vector<double>& data, double mean[], double var[], double *mixprob);

더보기
void initGuess(std::vector<double>& data, double mean[], double var[], double *mixprob) {
    int start = -1, end = data.size(); 
    // trim null data;
    while (data[++start] <= 0) ;
    while (data[--end] <= 0) ;
    // split given data into two equal size sets;
    int mid = (end + start) / 2;
    // simple mean and variance;
    estimGaussParams(data, start, mid, &mean[0], &var[0]);
    estimGaussParams(data, mid + 1, end, &mean[1], &var[1]);
    // initial guess for mixing probability;
    *mixprob = 0.5; 
};

#define PI (4.0 * atan(1.))

double gaussDist(double x, double mean, double var) ;  

더보기
double gaussDist(double x, double mean, double var) { 
    // N(mean, var);
    double arg = 0.5 * (x - mean) * (x - mean) / var;
    double factor = 1 / sqrt(2.* PI * var);
    return factor * exp(-arg); 
}

double responsibility2(double x, double mean[], double var[], double mixprob) ;   

더보기
double responsibility2(double x, double mean[], double var[], double mixprob) {   
    double a = (1 - mixprob) * gaussDist(x, mean[0], var[0]);
    double b = mixprob * gaussDist(x, mean[1], var[1]);  
    return b / (a + b); 
}

double weightedMeanVar(std::vector<double>& data, std::vector<double> & gamma, double mean[], double var[]) ;

더보기
double weightedMeanVar(std::vector<double>& data, std::vector<double>& gamma, double mean[], double var[]) { 
	// estimate new means;
    double s = 0, sx0 = 0, sx1 = 0, sg = 0;
    for (int i = data.size(); i-- > 0; ) {
        s   += data[i];
        sg  += data[i] * gamma[i]; 
        sx0 += data[i] * i * (1 - gamma[i]);
        sx1 += data[i] * i * gamma[i];
    }
    mean[0] = sx0 / (s - sg);
    mean[1] = sx1 / sg;
    // variances with new mean;
    double sv0 = 0, sv1 = 0;
    for (i = data.size(); i-- > 0; ) {
        sv0 += data[i] * (i - mean[0]) * (i - mean[0]) * (1 - gamma[i]);
        sv1 += data[i] * (i - mean[1]) * (i - mean[1]) * gamma[i];
    }
    var[0] = sv0 / (s - sg);
    var[1] = sv1 / sg;
    // return mixing probability = mixing ratio for class 2;
    return (sg / s);
};
#define EPSILON  1e-6
// Expectation Maximization algorithm applied to Two component Gaussian Mixture Model;
double emTwoCompGMM(std::vector<double>& data) {
    double mean[2], var[2], mixprob;
    std::vector<double> gamma(data.size());     // responsibilities for class 2;
    initGuess(data, mean, var, &mixprob);
    // begin algorithm;
    while (1) {
        // E-step;
        for (int i = data.size(); i-- > 0; ) 
            gamma[i] = responsibility2(i, mean, var, mixprob);
        double old_mixprob = mixprob;
        // M-step;
        mixprob = weightedMeanVar(data, gamma, mean, var);
        TRACE("mixing probability= %f\n", mixprob);
        // check convergence(usually loglikelihood is tested);
        if (fabs(mixprob - old_mixprob) < EPSILON)
            break;
    }
    // estimate decision boundary;
    int k = data.size();
    while (gamma[--k] >= 0.5) ;
    return (2 * k + 1) / 2.; // = average of ;
};

728x90

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

Kuwahara Filter  (2) 2020.12.28
Moving Average을 이용한 Thresholding  (0) 2020.11.26
Union-Find Connected Component Labeling  (0) 2012.11.01
RANSAC: Ellipse Fitting  (1) 2012.10.07
Autofocus Algorithm  (0) 2012.06.03
Posted by helloktk
,

각각 200개의 점들로 이루어진 8개의 2차원 가우시안 군집을 무작위로 만들고, 이를 kmeans 알고리즘을 써서 8개로 분할하였다. 아래의 시뮬레이션은 이 정보를 초기 조건으로 하여서 Gaussian Mixture Model (GMM)에 적용한 결과이다. 두 개의 군집에 대해서 kmeans 결과와 GMM의 결과가 서로 많이 차이가 남을 보여준다.

 


코드 추가: 2010.02.23

struct Data2d {
    double x, y ;
    int id ;
    Data2d() { };
    Data2d(double x, double y) : x(x), y(y), id(-1) { }
};
struct Gauss2d {    
    double cov[4];
    double mx, my ; //mean ;
    double mix;
    //
    double nfactor; //1/(2 * Pi * sqrt(det));
    double det;     //det(cov)
    double icov[4];
    void   prepare();
    //
    double pdf(double x, double y);
} ;
void Gauss2d::prepare() {
    // det(cov);
    det = cov[0] * cov[3] - cov[1] * cov[2];
    if (det < 1e-10) {
        AfxMessageBox("not converging");
        return ;
    };
    nfactor = 1. / (2. * MPI * sqrt(det)) ;
    //inv(cov);
    icov[0] =  cov[3] / det;
    icov[1] = -cov[1] / det;
    icov[2] = -cov[2] / det;
    icov[3] =  cov[0] / det;
}
double Gauss2d::pdf(double x, double y) {
    x -= mx ;
    y -= my ;
    double a = x * (icov[0] * x + icov[1] * y) +
               y * (icov[2] * x + icov[3] * y);
    return (nfactor * exp(-0.5 * a));
};
void init_classes(std::vector<Data2d>& data, std::vector<Gauss2d>& classes) {
    /*
    for (int i = 0; i < classes.size(); i++) {
        Gauss2d& cls = classes[i] ;
        cls.cov[0] = 10 + 50 * rand() / double(RAND_MAX);
        cls.cov[1] = 0;
        cls.cov[2] = 0;
        cls.cov[3] = 10 + 50 * rand() / double(RAND_MAX);
        cls.mx = 100 + 300 * rand() / double(RAND_MAX);   
        cls.my = 100 + 300 * rand() / double(RAND_MAX);   
        cls.mix = 1;
    }
    */
    KMeans(data, classes);
    //use kmeans to locate initial positions;
}
void test_step(std::vector<Data2d>& data,
               std::vector<Gauss2d>& classes,
               std::vector<std::vector<double> >& prob_cls) 
{
    //E-step ;
    for (int k = 0; k < classes.size(); k++) {
        Gauss2d& cls = classes[k];
        cls.prepare();
        //
        for (int i = 0; i < data.size(); i++) {
            prob_cls[i][k] = cls.mix * cls.pdf(data[i].x, data[i].y);
        };
    }
    // normalize-->임의의 데이터는 각 어떤 클레스에 속할 활률의 합=1;
    for (int i = 0; i < data.size(); i++) {
        double s = 0; 
        int bc = 0; double bp = 0;  // to determine membership(debug);
        for (int k = 0; k < classes.size(); ++k) {
            s += prob_cls[i][k];
            // find maximum posterior for each data;
            if (bp < prob_cls[i][k]) {
                bp = prob_cls[i][k] ;
                bc = k ;
            };
        }
        data[i].id = bc;
        // normalize to 1;
        for (int k = 0; k < classes.size(); ++k)
            prob_cls[i][k] /= s;
    }
    //M-step;
    for (int k = 0; k < classes.size(); k++) {
        Gauss2d & cls = classes[k];
        //get mean;
        double meanx    = 0;
        double meany    = 0;
        double marginal = 0; 
        for (int i = 0; i < data.size(); i++) {
            meanx    += prob_cls[i][k] * data[i].x ;
            meany    += prob_cls[i][k] * data[i].y ;
            marginal += prob_cls[i][k];
        };
        cls.mx = meanx = meanx / marginal ; 
        cls.my = meany = meany / marginal ;
        // get mixing;
        cls.mix = marginal / classes.size();
        // get stdev;
        double sxx = 0, syy = 0, sxy = 0;
        for (int i = 0; i < data.size(); i++) {
            double dx = data[i].x - meanx ;
            double dy = data[i].y - meany ;
            sxx += prob_cls[i][k] * dx * dx ;
            syy += prob_cls[i][k] * dy * dy ;
            sxy += prob_cls[i][k] * dx * dy ;
        };
        //set covariance;
        cls.cov[0] = sxx / marginal;
        cls.cov[1] = sxy / marginal;
        cls.cov[3] = syy / marginal;
        cls.cov[2] = cls.cov[1]; //symmetric;
    }   
}
void test() {
    int max_iter = 100;
    int nclass = 8;
    int ndata = 500;
    std::vector<Gauss2d> classes(nclass);
    std::vector<Data2d> data(ndata);
    // prepare posterior space;
    std::vector<std::vector<double> > prob_cls;
    for (int i = 0; i < data.size(); ++i) {
        prob_cls.push_back(std::vector<double>(classes.size()));
    } ;
   // generate data...
   ..................................
    //init_classes
    init_classes(data, classes) ;

    int iter = 0;
    do {
        iter++;
        test_step(data, classes, prob_cls);
    } while (iter < max_iter) ;
};
728x90

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

Savitzky-Golay Smoothing Filter  (2) 2010.03.24
Retinex 알고리즘  (11) 2010.02.03
Image Morphing  (0) 2010.01.24
Fant's Algorithm  (0) 2010.01.22
Affine Transformation  (0) 2010.01.20
Posted by helloktk
,

 

이미지의 히스토그램을 이용하여 전경과 배경을 분리하는 이진화는 가우시안 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
Posted by helloktk
,

주어진 N개의 관측값을 k개의 가우시안 mixture로 모델링한 후 EM 알고리즘을 적용하여 모델을 결정하는 파라미터를 추정하고자 한다. 추정해야 할 모델 파라미터는 k개의 가우시안을 기술하는 평균값(μ)과 공분산 행렬(Covariance Matrix: Σ)과 이들 가우시안 분포의 mixing(α) 정도를 나타내는 파라미터이다.

사용자 삽입 이미지

다음은 이 과정을 수행하는 C++ class와 사용 예제, 및 결과를 보여준다.

//가우시안 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 * dxy) * 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};
    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|θ);

사용자 삽입 이미지

 

728x90

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

EM: Binarization  (0) 2008.07.01
EM Algorithm: Line Fitting  (0) 2008.06.29
Rasterizing Voronoi Diagram  (0) 2008.05.26
RANSAC Algorithm  (0) 2008.05.24
Contour Tracing  (0) 2008.05.22
Posted by helloktk
,