K-means之C++及OpenCV实现

这个算法的步骤如下：

1.随机选取样本中的K个点作为聚类中心
2.计算所有样本到各个聚类中心的距离，将每个样本规划在最近的聚类中
3.计算每个聚类中所有样本的中心，并将新的中心代替原来的中心
4.检查新老聚类中心的距离，如果距离超过规定的阈值，则重复2-4，直到小于阈值

那么，现在，我实现的程序的步骤也是按照上面一步一步来的，

为了方便，我直接在平面上随机产生n个点，选取前K个点作为聚类中心，

距离就定义为平面上的欧式距离，

然后为了形象化地观察过程和结果，我将过程以图像的方式显示。

代码如下：

首先是主体：

[cpp] view plain copy

 print?

1. int iter_times = 0;//迭代次数  
2.     while(!good_result())//检查是否是需要的聚类中心  
3.     {  
4.         for(int i = 0;i < POINT_NUM;i++)  
5.         {  
6.             double min = 10000;  
7.             int min_k = 0;  
8.             for(int j = 0;j < K;j++)  
9.             {  
10.                 double tmp = DIS(c[j].center,s[i].p);   
11.                 if(tmp < min)  
12.                 {  
13.                     min = tmp;  
14.                     min_k = j;   
15.                 }  
16.             }  
17.             s[i].cluster = min_k;  
18.             c[min_k].count++;  
19.         }  
20.         update_center();//更新聚类中心  
21.         iter_times++;  
22.         show_outcome();  
23.     }  

然后是其他函数的实现：

[cpp] view plain copy

 print?

1. void update_center()  
2. {  
3.     double x[K],y[K];  
4.     memset(x,0,sizeof(x));  
5.     memset(y,0,sizeof(y));  
6.     for(int i = 0;i < POINT_NUM;i++)  
7.     {  
8.         x[s[i].cluster] += s[i].p.x;  
9.         y[s[i].cluster] += s[i].p.y;  
10.     }  
11.     for(int i = 0;i < K;i++)  
12.     {  
13.         c[i].pre_center = c[i].center;  
14.         c[i].center.x = x[i] / c[i].count;  
15.         c[i].center.y = y[i] / c[i].count;  
16.         c[i].count = 0;  
17.     }  
18. }  

判断是否是需要的：

[cpp] view plain copy

 print?

1. bool good_result()  
2. {  
3.     for(int i = 0;i < K;i++)  
4.     {  
5.         if(DIS(c[i].center,c[i].pre_center) > TH)  
6.             return false;  
7.     }  
8.     return true;  
9. }  

显示结果：

[cpp] view plain copy

 print?

1. void show_outcome()  
2. {  
3.     for(int y = 0;y < HEIGHT;y++)//这里将平面中所有的点都标记，就可以看到平面是怎样被划分的了  
4.         for(int x = 0;x < WIDTH;x++)  
5.         {  
6.             double min = 1000;  
7.             int min_k = 0;  
8.             CvPoint p = cvPoint(x,y);  
9.             for(int i = 0;i < K;i++)  
10.             {  
11.                 double tmp = DIS(c[i].center,p);   
12.                 if(tmp < min)  
13.                 {  
14.                     min = tmp;  
15.                     min_k = i;   
16.                 }  
17.             }  
18.             IMG_B(img,x,y) = color[min_k].val[0];  
19.             IMG_G(img,x,y) = color[min_k].val[1];  
20.             IMG_R(img,x,y) = color[min_k].val[2];  
21.             IMG_A(img,x,y) = 200;//4通道图像，就是说可以是透明的，纯试验而已，哪知道直接显示没效果，要保存之后才能看出来。  
22.         }  
23.     CvScalar scalar = cvScalar(255,255,255,255);  
24.     for(int i = 0;i < POINT_NUM;i++)//画每个样本点  
25.     {  
26.         int x = s[i].p.x;  
27.         int y = s[i].p.y;  
28.         cvLine(img,cvPoint(x - 5,y),cvPoint(x + 5,y),scalar,2);  
29.         cvLine(img,cvPoint(x,y - 5),cvPoint(x,y + 5),scalar,2);  
30.     }  
31.     for(int i = 0;i < K;i++)//画聚类中心  
32.     {  
33.         int x = c[i].center.x;  
34.         int y = c[i].center.y;  
35.         cvCircle(img,cvPoint(x,y),20,scalar,2);  
36.     }  
37.     cvShowImage("Image",img);  
38.     cvWaitKey(100);//100毫秒是个差不多的数值，可以完整的看到聚类过程  
39. }  

看下效果：

而通过几次运行和观察，阈值不必取的很小，首先是迭代次数越来越多，时间越来越长，但结果差别却是越来越小，即，到几次迭代之后就能取得好的效果了，再迭代下去取的结果跟原来相差不大。

const int nClusters = 4;//这是Kmeans算法的一个缺点，在聚类之前需要指定类别个数
int _tmain(int argc, _TCHAR* argv[])
{
Mat src; 
src = imread("E:\\bad\\belt  (1).jpeg"); 
imshow("original", src);

blur(src, src, Size(11,11));//使用blur对图像进行平滑处理，这种方法就是最简单的求平均数
                           //size:定义滤波器的大小
    imshow("blurred", src);

//p是特征矩阵，每行表示一个特征，每个特征对应src中每个像素点的(x,y,r,g,b共5维)
    Mat p = Mat::zeros(src.cols*src.rows, 5, CV_32F);    //初始化全0矩阵(行和列还有类型)
                    //所有元素行，5列的矩阵
Mat bestLabels, centers, clustered;
vector<Mat> bgr;   //定义一个Mat向量容器保存拆分后的数据 
    split(src, bgr);    //分隔出src的三个通道


for(int i=0; i<src.cols*src.rows; i++) 
{
p.at<float>(i,0) = (i/src.cols) / src.rows;   // p.at<uchar>(y,x) 相当于 p->Imagedata[y *p->widthstep + x], p是8位uchar
p.at<float>(i,1) = (i%src.cols) / src.cols;   // p.at<float>(y,x) 相当于 p->Imagedata[y *p->widthstep + x], p是32位float
p.at<float>(i,2) = bgr[0].data[i] / 255.0;
p.at<float>(i,3) = bgr[1].data[i] / 255.0;
p.at<float>(i,4) = bgr[2].data[i] / 255.0;
}
    //计算时间
    double t = (double)cvGetTickCount();//GetTickcount函数：它返回从操作系统启动到当前所经的计时周期数

//kmeans聚类，每个样本的标签保存在bestLabels中
kmeans(p, nClusters, bestLabels,
        TermCriteria( CV_TERMCRIT_EPS+CV_TERMCRIT_ITER, 10, 1.0),
        3, KMEANS_PP_CENTERS, centers);


    t = (double)cvGetTickCount() - t;
    float timecost = t/(cvGetTickFrequency()*1000);//getTickFrequency函数：返回每秒的计时周期数


    //给每个类别赋颜色，其值等于每个类第一个元素的值
    Vec3b    colors[nClusters];//Vec3b 是定义一个uchar类型的数组长度为3
    bool    colormask[nClusters]; 
memset(colormask, 0, nClusters*sizeof(bool));//将colormask中前nClusters*sizeof(bool)个字节用0替换并返回colormask
//memset：作用是在一段内存块中填充某个给定的值，
//它是对较大的结构体或数组进行清零操作的一种最快方法




    int        count = 0;
    for(int i=0; i<src.cols*src.rows; i++) 
    {
        int clusterindex = bestLabels.at<int>(i,0);
          for (int j=0; j<nClusters; j++)
        {
            if(j == clusterindex && colormask[j] == 0)
            {
                int y = i/src.cols;
                int x = i%src.cols;
                 colors[j] = src.at<Vec3b>(y,x);
               colormask[j] = 1;
                count++;
                 break;
            }
        }
         if(nClusters == count)break;
    }


    //显示聚类结果
    clustered = Mat(src.rows, src.cols, CV_8UC3);
    for(int i=0; i<src.cols*src.rows; i++)
{
        int y = i/src.cols;
        int x = i%src.cols;
        int clusterindex = bestLabels.at<int>(i,0);
        clustered.at<Vec3b>(y, x) = colors[clusterindex];
    }


    imshow("clustered", clustered);

cvWaitKey(0);