K-means之C++及OpenCV实现

这个算法的步骤如下：

1.随机选取样本中的K个点作为聚类中心

2.计算所有样本到各个聚类中心的距离，将每个样本规划在最近的聚类中

3.计算每个聚类中所有样本的中心，并将新的中心代替原来的中心

4.检查新老聚类中心的距离，如果距离超过规定的阈值，则重复2-4，直到小于阈值

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

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

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

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

代码如下：

首先是主体：

[cpp] view
plain copy

 print?

int iter_times = 0;//迭代次数  

    while(!good_result())//检查是否是需要的聚类中心  

    {  

        for(int i = 0;i < POINT_NUM;i++)  

        {  

            double min = 10000;  

            int min_k = 0;  

            for(int j = 0;j < K;j++)  

            {  

                double tmp = DIS(c[j].center,s[i].p);   

                if(tmp < min)  

                {  

                    min = tmp;  

                    min_k = j;   

                }  

            }  

            s[i].cluster = min_k;  

            c[min_k].count++;  

        }  

        update_center();//更新聚类中心  

        iter_times++;  

        show_outcome();  

    }  

然后是其他函数的实现：

[cpp] view
plain copy

 print?

void update_center()  

{  

    double x[K],y[K];  

    memset(x,0,sizeof(x));  

    memset(y,0,sizeof(y));  

    for(int i = 0;i < POINT_NUM;i++)  

    {  

        x[s[i].cluster] += s[i].p.x;  

        y[s[i].cluster] += s[i].p.y;  

    }  

    for(int i = 0;i < K;i++)  

    {  

        c[i].pre_center = c[i].center;  

        c[i].center.x = x[i] / c[i].count;  

        c[i].center.y = y[i] / c[i].count;  

        c[i].count = 0;  

    }  

}  

判断是否是需要的：

[cpp] view
plain copy

 print?

bool good_result()  

{  

    for(int i = 0;i < K;i++)  

    {  

        if(DIS(c[i].center,c[i].pre_center) > TH)  

            return false;  

    }  

    return true;  

}  

显示结果：

[cpp] view
plain copy

 print?

void show_outcome()  

{  

    for(int y = 0;y < HEIGHT;y++)//这里将平面中所有的点都标记，就可以看到平面是怎样被划分的了  

        for(int x = 0;x < WIDTH;x++)  

        {  

            double min = 1000;  

            int min_k = 0;  

            CvPoint p = cvPoint(x,y);  

            for(int i = 0;i < K;i++)  

            {  

                double tmp = DIS(c[i].center,p);   

                if(tmp < min)  

                {  

                    min = tmp;  

                    min_k = i;   

                }  

            }  

            IMG_B(img,x,y) = color[min_k].val[0];  

            IMG_G(img,x,y) = color[min_k].val[1];  

            IMG_R(img,x,y) = color[min_k].val[2];  

            IMG_A(img,x,y) = 200;//4通道图像，就是说可以是透明的，纯试验而已，哪知道直接显示没效果，要保存之后才能看出来。  

        }  

    CvScalar scalar = cvScalar(255,255,255,255);  

    for(int i = 0;i < POINT_NUM;i++)//画每个样本点  

    {  

        int x = s[i].p.x;  

        int y = s[i].p.y;  

        cvLine(img,cvPoint(x - 5,y),cvPoint(x + 5,y),scalar,2);  

        cvLine(img,cvPoint(x,y - 5),cvPoint(x,y + 5),scalar,2);  

    }  

    for(int i = 0;i < K;i++)//画聚类中心  

    {  

        int x = c[i].center.x;  

        int y = c[i].center.y;  

        cvCircle(img,cvPoint(x,y),20,scalar,2);  

    }  

    cvShowImage("Image",img);  

    cvWaitKey(100);//100毫秒是个差不多的数值，可以完整的看到聚类过程  

}  

看下效果：



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

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);