GMM EM C CODE c代码

/**/ ///
// function: calculate u and sigma with histogram
// summary : get u and sigma
// parameters:
//           pImg : the image
//           num  : total pixels
//           u    : 均值  引用传递 返回
//           sigma2: 方差  注意不是平均差
// 引用传递 返回u和sigma^2 

int  GetUAndSigma( const  IplImage  * pImg,  int  num,  float   & u,  float   & sigma2)  
... {
    int i;
    int histo[256];    //个数统计
    float sum;
    
    
    sum = 0.0;
    
    for(i = 0; i < 256; i++)
    ...{
        histo[i] = 0;
    }
    
    //histogram
    for(i = 0; i < num; i++ )
    ...{    
        histo[(unsigned char)pImg->imageData[i]]++;
    }
    
    for(i = 0; i < 256; i++)
    ...{
        sum += i * histo[i];    
    }
    
    u = sum / num;       //平均值
    
    sum = 0.0;
    
    for(i = 0; i < 256; i++)
    ...{
        sum += (u - i) * (u - i) * histo[i] ;
    }
    
    sigma2 = sum / num;  //方差   注意：不是标准差
    
    return true;
    
}

/**/ //
// summary: 利用flag 标志更新各个类的均值和方差
//
//
void  UpdateUandSigma( const  IplImage  * pImg,  int  num,  int  k,  char   * flag,  int   * center,  int   * sigma)
... {
    int i,j;
    int histo[256];    //个数统计
    float sum;
    int n ;
    
    for(j = 0; j < k; j++)
    ...{
        sum = 0.0;
        
        n = 0;
        
        for(i = 0; i < 256; i++)
        ...{
            histo[i] = 0;
        }
        
        //histogram
        for(i = 0; i < num; i++ )
        ...{    
            if(flag[i] == j)
            ...{
                histo[(unsigned char)pImg->imageData[i]]++;
                n++;
            }
        }
        
        for(i = 0; i < 256; i++)
        ...{
            sum += i * histo[i];    
        }
        
        center[j] = sum / n;       //平均值
        
        sum = 0.0;
        
        for(i = 0; i < 256; i++)
        ...{
            sum += (center[j] - i) * (center[j] - i) * histo[i] ;
        }
        
        sigma[j] = sum / n;  //方差   注意：不是标准差
    }
    
}

/**/ //
//
//
//
void  GMM( char   * fileName,  int  k)
... {
    //-------opencv 读取图象相关数据----------
    IplImage *pImg;
    IplImage *pImgK;
    
    int width, height, step;
    int num;          //total pixels  width * height
    
    pImg = cvLoadImage(fileName, CV_LOAD_IMAGE_GRAYSCALE);
    
    if(NULL == pImg)
    ...{
        printf("image open fail ");
        return ;
    }
    
    pImgK = cvCloneImage(pImg);
    width = pImg->width;
    height = pImg->height;
    step = pImg->widthStep;
    num = width * height;
    
    //--------------------------------------------------
    
    int i, j, ii;  //loop variant 
    
    int *center;   //clustering center  均值
    int *sigma;    //方差 每个类一个方差 注意：是方差不是标准差
    float *a;        //ai i = 0..k-1 sum(ai) == 1
    float *z;        // k * num matrix  p(xk|wj)
    float *p;        // k*num matrix  p(wj|xk)
    char *flag;      //indicate which class the points are in
    unsigned char data; //pixel value
    
    float u, sigma2;  //
    int iteration;
    
    center = (int*) malloc(k * sizeof(int));  //for memory reason , substitute int for char
    sigma = (int*) malloc(k * sizeof(int));  //方差数组  注意：是方差不是标准差
    a = (float*) malloc(k*sizeof(float));        //
    z = (float*) malloc(k * num * sizeof(float));  //
    p = (float*) malloc(k * num * sizeof(float));
    flag = (char*) malloc(num * sizeof(char));
    
    //------函数初始化------
    
    GetUAndSigma(pImgK, num, u, sigma2);
    
    printf("the whole image's u = %f  sigma2 = %f  ", u, sigma2);   //test
    
    for(i = 0; i < k; i++)
    ...{
        a[i] = 1.0 / k;     //初始化每部分概率为1/k 各部分概率相同 for p
    }
    
    for(i = 0; i < (k*num); i++)
    ...{
        z[i] = 0.0;   
        p[i] = 0.0;
    }
    
    //srand((unsigned) time(NULL));
    for(i = 0; i < num; i++)
    ...{
        flag[i] = rand() % k;                //随机初始化各个类 
    }
    //---------------------------------------
    
    UpdateUandSigma(pImg, num, k, flag, center, sigma);
    
    
    //主循环
    iteration = 0;
    float maxLike_old; 
    float    maxLike;
    
    maxLike = FLT_MAX;
    
    
    do
    ...{
        
        printf("center[0] = %d, center[1] = %d  ", center[0], center[1]); //test
        printf("sigma[0] = %d, sigma[1] = %d  ",sigma[0], sigma[1]);
        printf("a[0] = %f, a[1] = %f ", a[0], a[1]);
        
        maxLike_old = maxLike;
        
        //E-step
        //---------------p(xn|wk)--------------------
        for(i = 0; i < num; i++)
        ...{
            data  = pImg->imageData[i] ;
            
            for(j = 0; j < k; j++)
            ...{    
                z[j*num + i] = exp(-pow((data - center[j]), 2.0) / (2* sigma[j] + FLT_MIN)) / (sqrt(2*PI*sigma[j]) + FLT_MIN) ;
            }
        }    //计算每个点在对应类的正态分布值
        
        //e-step
        float maxProb = FLT_MIN;
        char currentClass = -1;
        float temp;
        
        for(i = 0; i < num; i++)
        ...{
            currentClass = -1;  //进入前清空
            maxProb = FLT_MIN;
            temp = 0.0;
            
            for(ii = 0; ii < k; ii++)
            ...{
                temp += z[ii*num + i] * a[ii];    //分母啊
            }
            
            for(j = 0; j < k; j++)
            ...{
                p[j*num + i] = z[j*num + i]* a[j] / temp ;      //p(wj|xk)
                if(p[j*num + i] > maxProb)
                ...{
                    maxProb = p[j*num + i];
                    currentClass = j ;
                }
            }    //bayes公式 点在哪个类的概率大 划到哪个类去
            
            flag[i] = currentClass;
            
        }
        
        //        for(i = 0; i < num; i++)
        //        {
        //            printf("flag = %d ", flag[i]);
        //        }
        // UpdateUandSigma(pImg, num, k, flag, center, sigma);
        //------------------------------------------------------------------
        
        
        //M-step
        float denominator; //分母
        float numerator ;  //分子
        //m-step
        
        for(j = 0; j < k; j++)
        ...{
            numerator = 0.0;
            denominator = 0.0;
            
            for(i = 0; i < num; i++)
            ...{
                data = pImg->imageData[i];
                numerator += p[j*num +i] * data;
                denominator += p[j*num +i];
            }
            
            center[j] = numerator / denominator ;
            a[j] = denominator / num ;
        }
        
        //UpdateUandSigma(pImg, num, k, flag, center, sigma);
        
        
        // 计算max likehood for GMM
        maxLike = 0.0;
        
        
        for(i = 0; i < num; i++)
        ...{
            temp = 0.0;
            for(j = 0; j < k; j++)
            ...{
                temp += a[j]*z[j*num + i] ;
            }
            maxLike += log(temp);
        }
        
        iteration++;
        
        printf("maxLike = %f ", maxLike);
        printf("maxLike_old = %f  ", maxLike_old);
        
    }while(fabs(maxLike - maxLike_old) > 0.000001);   //fabs(maxLike - maxLike_old) > 0
    
    printf("center[0] = %d, center[1] = %d  ", center[0], center[1]); //test
    printf(" iteration = %d  ", iteration);
    
    
    
    printf(" z[i][j]  ");
    /**//*    for(i = 0; i < num; i++)
    {
    
      for(j = 0; j < k; j++)
      {
      printf(" %f ", z[j*num +i]);
      }
      printf("  ");
      }
      
        printf(" p[i][j]  ");
        for(i = 0; i < num; i++)
        {
        
          for(j = 0; j < k; j++)
          {
          printf(" %f ", p[j*num +i]);
          }
          printf("  ");
    } //*/
    
    //-------------display clustering results-----------------
    int color = 255 / k;
    for(i = 0; i < height; i++)
        for(j = 0; j < width; j++)
        ...{
            pImgK->imageData[i*step +j] = (unsigned char) flag[i*step + j] *  color ;
            
        }
        
        
        //--------display image-----------------------------------------
        cvNamedWindow("source", CV_WINDOW_AUTOSIZE);
        cvNamedWindow("GMM", CV_WINDOW_AUTOSIZE);
        
        cvShowImage("source", pImg);
        cvShowImage("GMM",pImgK);
        
        cvWaitKey(0);
        
        cvDestroyWindow("source");
        cvDestroyWindow("GMM");
        
        cvReleaseImage(&pImg);
        cvReleaseImage(&pImgK);
        
        
        free(center);    //release memory
        center = NULL;
        free(a);
        a = NULL;
        free(z);
        z = NULL;
        free(flag);
        flag = NULL;
        
}