각각 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
,

KMeans Algorithm

Image Recognition 2008. 7. 19. 21:03

KMeans 알고리즘은 주어진 데이터들을 사전에 정해진 수만큼의 클러스터로 자동 분류하는 간단한 클러스터링 알고리즘의 일종이다. 이 알고리즘의 주된 아이디어는 각 클러스터의 중심을 어떻게 하면 잘 잡는가이다. 좋은 선택은 각각의 중심이 가능한 서로 멀리 떨어지게 잡아야 한다. 일단 클러스터 중심이 설정되면 각각의 데이터에 대해서 가장 가까이에 있는 클러스터에 할당시키면 된다. 다음으로 클러스터의 중심을 그 클러스터에 속한 데이터의 무게중심으로 재설정하고, 다시 각각의 데이터를 클러스터에 할당하는 작업을 반복한다. 이 과정은 클러스터의 중심이 더 이상 이동하지 않을 때까지 반복 시행하면 된다. 결과적으로 이 알고리즘은 클러스터의 중심에서 그 클러스터에 할당된 데이터와의 거리의 제곱을 그 클러스터의 모든 점과 전체 클러스터에 대해서 합한 값을 최소화시키는 클러스터 중심과 데이터의 분할을 찾는 것이 된다.

1. k개의 클러스터의 중심을 설정한다.
2. 설정된 중심을 가지고, 각 데이터를 분할한다.
3. 분할된 데이터를 이용하여서 클러스터의 중심을 재설정한다.
4. 중심이 변화가 없을 때까지 2,3 과정을 반복한다.

비용 함수의 관점에서 보면 KMeans 클러스터링은 다음 함수를

$$ \text{cost}=\sum_{k \in classes} \sum_{i \in class_k} \big| \text{data}_k(i) - \text{mean}_k \big|^2 $$

을 최소화시키는 분할을 하는 것이다. 그러나,  KMeans는 항상 최적의 비용 함숫값을 보장하지는 않는다. 그리고 알고리즘의 결과는 초기값의 설정에 의존하여서 다른 결과를 낼 수가 있다. 보다 큰 단점은 클러스터의 수를 사전에 정해주어야 한다는 사실인데, 잘못된 클러스터 개수의 설정은 성능에 크게 영향을 미친다. 자동의 클러스터의 수를 설정하면서 클러스터링을 할 수 있는 간단한 것으로는 isodata 알고리즘이 있다.

typedef double * vector_t ;
double kmeans_label(
             const vector_t *obs,             // observations;
             const int n_obs,                   /* in # of vectors */            
             const vector_t *mean,           // mean vectors;
             const int n_mean,                /* # of mean vectors */
             int* label,                           // cluster labeling;
             const int veclen)                 // dimension of feature;
{
    double sqerr=0;
    for (int i = 0; i < n_obs; i++) {
        vector_t c = obs[i];
        /* Get an estimate of best distance (min_dist) and (min_idx) */
        int min_idx = label[i];
        vector_t m = mean[min_idx];
        double min_sqdist = 0 ;
        for (int l = 0; l < veclen; l++) {
            double t = m[l] - c[l]; min_sqdist += t * t;
        }
        for (int j = 0; j < n_mean; j++) {
            vector_t m = mean[j];
            // save time by checking:: dist < min_sqdist=previous min_val;
            double sqdist = 0;
            for (l = 0; (l < veclen) && (sqdist < min_sqdist); l++) {
                double t = m[l] - c[l]; sqdist += t * t;
            }
            if (sqdist < min_sqdist) {
                min_sqdist = sqdist; min_idx = j;
            }
        }
        label[i] = min_idx;
        sqerr += min_sqdist;
    }
    return sqerr ;
};
int kmeans_update(
              const vector_t *obs, //observations;
              const int n_obs,       //# of observations;     
              vector_t *mean,       // mean vectors;
              const int n_mean,    // # of mean vectors;
              int *label,               // cluster labelling;
              const int veclen      // dimension of features;
              ) {
    std::vector<int> cnt(n_mean);
    for (int i = 0; i < n_mean; i++)
        for (int l = 0; l < veclen; l++) mean[i][l] = 0.0;
        
    for (i = 0; i < n_obs; i++) {
        ASSERT((0 <= label[i]) && (label[i] < n_mean));  
        vector_t m = mean[label[i]];
        cnt[label[i]]++;  
        vector_t c = obs[i];
        if (c == NULL)
            TRACE("No observations for %u, but expected up through %u", i, n_obs-1);
        for (int l = 0; l < veclen; l++) m[l] += c[l];
    }
    for (i = 0; i < n_mean; i++) {
        int j = cnt[i];
        if (j != 0) for (int l = 0; l < veclen; l++) mean[i][l] /= (double) j;
        else {
            TRACE("Empty cluster %u\n", i);
            return FALSE;
        }
    }
    return TRUE;
}
double kmeans(
       const vector_t *obs,     // observations;
       const int n_obs,          /* # of observations */
       vector_t *mean,           /* initial set of means */
       const int n_mean,        /* # of means (should be k_mean?) */
       int *label,                   /* cluster labeling*/
       const int veclen,         /* vector length of means and corpus */    
       const double min_conv_ratio,
       const int max_iter
    ) {
    double p_sqerr = DBL_MAX;
    int ret;
    double sqerr = kmeans_label( obs, n_obs,mean, n_mean, &label[0], veclen);
    double conv_ratio = (p_sqerr - sqerr) / p_sqerr;
    for (int i = 0; (i < max_iter) && (conv_ratio > min_conv_ratio); i++) {
        TRACE("kmeans iter [%u] %e ...\n", i, conv_ratio);
        ret = kmeans_update(obs, n_obs, mean, n_mean, &label[0], veclen );
        if (ret != TRUE)
            return (double)ret;
        p_sqerr = sqerr;
        sqerr = kmeans_label(  obs, n_obs,mean, n_mean, &label[0], veclen);
        conv_ratio = (p_sqerr - sqerr) / p_sqerr;
    }
    TRACE("kmeans n_iter %u sqerr %e conv_ratio %e\n", i, sqerr, conv_ratio);   
    return sqerr;
}



사용자 삽입 이미지

728x90

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

Retinex Algorithm  (2) 2008.07.26
RANSAC: Circle Fit  (0) 2008.07.21
Robust Line Fitting  (0) 2008.07.08
EM: Binarization  (0) 2008.07.01
EM Algorithm: Line Fitting  (0) 2008.06.29
Posted by helloktk
,