fcm code C代码

 

/**/ ////
// function    : Fuzzy C Means clustering
// parameter 1 : fileName to image file to be processed
// parameter 2 : the number of clustering 


void  FuzzyCMeans( char   * fileName,  int  k)
... {
    //-----------算法本身参数---------------
    double m ; //coefficient  default = 2.0;
    double sigma ;  // to stop iteration 
    
    m = 2.0;
    sigma = 0.000001 ;
    
    //----------打开图象相关----------------
    IplImage *pImg;
    IplImage *pImgK;   //after clustering
    
    int width, height, step;
    int i, j;
    int num;      //total pixels
    
    int *centers;   //clustering centers
    float *u;         //membership matrix  num X k
    char *flag;      //所属类别标志
    
    
    pImg = cvLoadImage(fileName, CV_LOAD_IMAGE_GRAYSCALE);
    
    if(NULL == pImg)
    ...{
        printf("file open error!");
        return ;
    }
    
    pImgK = cvCloneImage(pImg);
    width = pImg->width;
    height = pImg->height;
    step = pImg->widthStep;
    num = width * height;
    
    centers = (int*) malloc(k * sizeof(int));   //
    u = (float*) malloc(num * k * sizeof(float));
    flag = (char*) malloc(num * sizeof(char));
    
    //-----------var initialise------------------
    
    //-----counts & centers-----
    for(i = 0; i < k; i++)
    ...{
        centers[i] = rand() % 256;          //随机初始化中心
    }
    
    
    //循环内部使用到的变量
    float maxU, maxU_old;
    
    
    
    //--------membership matrix initial--------  
    //暂时用随机数来初始化
    maxU = FLT_MIN ;
    
    
    //calculate u[i][j] i for clustering  j is the N points
    //更新U矩阵用到的变量
    maxU = FLT_MIN ;
    float sum1;
    float fenzi;
    int jj;
    
    for(i = 0; i <num; i++)
        for(j = 0; j < k; j++ )
        ...{
            unsigned char data = pImgK->imageData[i];
            fenzi = 1.0 / (pow((data - centers[j]), 2.0) + FLT_MIN);
            
            sum1 = 0.0;
            for(jj = 0; jj <k; jj++)
            ...{
                sum1 +=  1.0 / (pow((data - centers[jj]), 2.0) + FLT_MIN);
                
            }
            
            u[j*num + i] = fenzi / (sum1+ FLT_MIN);
            if(u[j * num + i] > maxU)
            ...{
                maxU = u[j * num +i] ;
            }
            
            
        }
        //随机初始化U矩阵
        //srand((unsigned) time(NULL));   
        /**//*    for(i = 0; i < k; i++)
        for(j = 0; j < num; j++)
        {
        
          u[i * num + j] = (float) (rand() % 11) / 10.0 ;//1.0 / temp;     //(rand() % 10) / 10 ;   //all between 0 and 1
          
} */
        
        float *sum = (float*) malloc(num * sizeof(float)) ;
        //-----------规格化处理----------------
        for(i = 0; i < num; i++)     //使矩阵满足约束条件   c X n  matrix
        ...{
            sum[i] = 0.0;
            for(j = 0; j < k; j++)  //一个数据集的隶属度恒等于 1
            ...{
                sum[i] += u[j * num + i];
            }
        }
        
        for(i = 0; i < num; i++)
        ...{
            for(j = 0; j < k; j++)
            ...{
                u[j * num + i] /= sum[i] ; 
                if(u[j * num + i] > maxU)
                ...{
                    maxU = u[j * num +i] ;
                }
                
                //printf(" u = %f", u[j * num + i]);   //test
            }
            //printf(" ");
        }
        
        
        
        
        //-----------------主循环---------------
        int iteration = 0;
        
        printf("before maxu = %f ", maxU);   //test
        
        do
        ...{
            //保存前后结果
            maxU_old = maxU ;
            
            //计算U矩阵用到的变量
            float temp;
            float tempB;
            
            //calculate  C[i]
            for(i = 0; i < k; i++)
            ...{
                temp = 0.0;
                tempB = 0.0;
                
                for(j = 0; j < num; j++)
                ...{
                    unsigned char data = pImgK->imageData[j];
                    temp +=   pow(u[i * num + j] , m) * data;// Data(i,j) data confine to 0--255
                    tempB +=  pow(u[i * num + j], m);
                }
                
                centers[i] =  temp / tempB;
                printf("center[%d]= %d ", i, centers[i]);
            }
            
            
            //calculate u[i][j] i for clustering  j is the N points
            //更新U矩阵用到的变量
            maxU = FLT_MIN ;
            //float sum1;
            //float fenzi;
            //int jj;
            
            for(i = 0; i <num; i++)
                for(j = 0; j < k; j++ )
                ...{
                    unsigned char data = pImgK->imageData[i];
                    fenzi = 1.0 / (pow((data - centers[j]), 2.0) + FLT_MIN);
                    
                    sum1 = 0.0;
                    for(jj = 0; jj <k; jj++)
                    ...{
                        sum1 +=  1.0 / (pow((data - centers[jj]), 2.0) + FLT_MIN);
                        
                    }
                    
                    u[j*num + i] = fenzi / (sum1+ FLT_MIN);
                    if(u[j * num + i] > maxU)
                    ...{
                        maxU = u[j * num +i] ;
                    }
                    
                    
                }
                
                
                
                //printf("in loop maxu = %f  maxu_old = %f  ", maxU, maxU_old);        //test
                
                
                iteration++ ;
                
                
        }while(fabs(maxU - maxU_old) > 0.0);//sigma);
        
        printf("iteration = %d  ", iteration);
        
        printf("matrix U  ");
        
        
        //---------display u matrix---------- for test
        float temp;
        
        for(i = 0; i < num; i++)
        ...{
            temp = FLT_MIN;
            
            for(j = 0; j < k; j++)
            ...{
                //    printf(" %f ", u[j*num + i]);
                if(u[j*num + i] > temp)
                ...{
                    temp = u[j*num + i];
                    flag[i] = j;
                }
            }
            //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("Fcmeans", CV_WINDOW_AUTOSIZE);
            
            cvShowImage("source", pImg);
            cvShowImage("Fcmeans",pImgK);
            
            cvWaitKey(0);
            
            cvDestroyWindow("source");
            cvDestroyWindow("Fcmeans");
            
            cvReleaseImage(&pImg);
            cvReleaseImage(&pImgK);
            
            free(centers);
            centers = NULL;
            free(u);
            u = NULL;
            
            free(sum);
            sum = NULL;
}