k-means算法实现图像颜色聚类

k-means算法原理比较简单，有关介绍参考http://blog.youkuaiyun.com/cai0538/article/details/7061922。原始算法最大的缺陷就是初始聚类中心的选取对效果的影响很大，对此有很多研究者提出了解决方法，网上一搜有大量论文。其中k-means++算是比较简单效果较好的了，参考http://blog.youkuaiyun.com/loadstar_kun/article/details/39450615，http://blog.youkuaiyun.com/ihsin/article/details/41809633。
这里实现了对彩色图像的聚类，效果和opencv的KMeans2方法没差别，只是速度很慢。。

#include<stdio.h>  
#include <cstdio>  
#include<string>  
#include<math.h>  
#include<stdlib.h>   
#include<vector>  
#include<string.h>  
#include<iostream>  

#include <cv.h>
#include <highgui.h>
using namespace cv;
using namespace std;

struct Tuple{
    int x;
    int y;
    unsigned char b;
    unsigned char g;
    unsigned char r;
};


const int K=4;//簇的数目  

vector<Tuple> node;//存储原始数据  
Tuple means[K];//存储均值中心  
vector<Tuple> clusters[K];//存储k个聚类  
int number;//输入数据的数目  
int dimension;//输入数据的维数  

//计算欧式距离
 double getDistXY(Tuple t1,Tuple t2)
 {
     return sqrt((float)((t1.b-t2.b)*(t1.b-t2.b)+(t1.g-t2.g)*(t1.g-t2.g)+(t1.r-t2.r)*(t1.r-t2.r)));
 }

 //根据质心,决定当前元组属于哪个簇
 int clusterOfTuple(Tuple means[],Tuple tuple)
 {
     float dist=getDistXY(means[0],tuple);
     float tmp;
     int label=0;//标示属于哪一簇
     for(int i=1;i<K;i++)
     {
         tmp=getDistXY(means[i],tuple);
         if(tmp<dist)
         {
             dist=tmp;
             label=i;
         }
     }
     return label;
 }

 //获得给定簇集的平方误差
 float getVar(vector<Tuple> clusters[],Tuple means[])
 {
     float var=0;
     for(int i=0;i<K;i++)
     {
         vector<Tuple> m=clusters[i];
         for(int j=0;j<m.size();j++)
         {
             var+=getDistXY(m[j],means[i]);
         }
     }
     return var;
 }

 //获得当前簇的均值（质心）
 Tuple getMeans(vector<Tuple> cluster)
 {
     int num=cluster.size();
     double meansX=0,meansY=0,meansZ=0;
     Tuple t;
     for(int i=0;i<num;i++)
     {
         meansX+=cluster[i].b;
         meansY+=cluster[i].g;
         meansZ+=cluster[i].r;
     }
     t.b=meansX/num;
     t.g=meansY/num;
     t.r=meansZ/num;
     return t;
 }
 void ChooseSeeds()//选择k个点作为初始的聚类中心，此处用k-means++算法优化，不是随机选择  
{  
    number=node.size();
    srand((unsigned int) time(NULL));  
    int idx=rand()%number;  
    int cnt=0;  
    means[cnt++]=node[idx];//记录选择的均值中心  
    double* dis=(double*)malloc(number*sizeof(double));//记录每个点距离它最近的均值中心的距离  
    memset(dis,0x3f,sizeof(dis));  
    while(cnt<K)//求出每个点与距离它最近的均值中心的距离  
    {  
        double sum=0;  
       for(int i=0;i<number;i++)  
        {  
            for(int j=0;j<cnt;j++)  
            {  
                if(i==j) continue;  
                dis[i]=min(dis[i],getDistXY(node[i],means[j]));  
            }  
            sum+=dis[i];  
        }  
        for(int i=0;i<number;i++)//归一化，其后可以对应到概率  
        {  
            dis[i]/=sum;  
        }  
        double* cumprobs=(double*)malloc(number*sizeof(double)); 
        cumprobs[0]=dis[0];  
        for(int i=1;i<number;i++)//求出概率的前缀和  
        {  
            cumprobs[i]=dis[i]+cumprobs[i-1];  
        }  
        bool* used=(bool*)malloc(number*sizeof(bool)); 
        memset(used,true,sizeof(used));  
        used[idx]=false;  
        next:  
        double r= (rand()%1000)*0.001;  
        bool flg=true;  
        for(int i=0;i<number;i++)  
        {  
            if(r<cumprobs[i]&&used[i])//选择满足概率的点作为簇中心  
            {  
                idx=i;  
                flg=false;  
                used[i]=false;  
                break;  
            }  
        }  
        if(flg) goto next; //如果没有找到，重新产生随机数r继续找  
        means[cnt++]=node[idx];  
        free(cumprobs);
        free(used);
    }  

    free(dis);

}  

 void KMeans(vector<Tuple> tuples)
 {

     int i=0;
     ChooseSeeds();
     int label=0;
     //根据默认的质心给簇赋值
     for(i=0;i!=tuples.size();++i)
     {
         label=clusterOfTuple(means,tuples[i]);
         clusters[label].push_back(tuples[i]);
     }

     float oldVar=-1;
     float newVar=getVar(clusters,means);
     int times=0;
     while(abs(newVar-oldVar)>=1&&times<100)//当新旧函数值相差不到1即准则函数不发生明显变化时，算法终止
     {
         for(i=0;i<K;i++)//更新每个簇的中心点
         {
             means[i]=getMeans(clusters[i]);

         }
         oldVar=newVar;
         newVar=getVar(clusters,means);//计算新的准则函数值
         for(i=0;i<K;i++)
         {
             clusters[i].clear();
         }
         //根据新的质心获得新的簇
         for(i=0;i!=tuples.size();++i)    
         {
             label=clusterOfTuple(means,tuples[i]);
             clusters[label].push_back(tuples[i]);
         }
         times++;

     }

 }



 int main()  
{  

    IplImage *img=cvLoadImage("horses2.jpg",-1);
    int img_w=img->width;
    int img_h=img->height;

    number=img_w*img_h;
    dimension=3;

    int i,j;
    Tuple t;
    for ( i = 0; i < img_h; i++)
    {
        for ( j = 0; j < img_w; j++)
        {
            t.x=j;
            t.y=i;
            t.b=(unsigned char)*(img->imageData + i*img->widthStep+3*j);
            t.g=(unsigned char)*(img->imageData + i*img->widthStep+3*j+1);
            t.r=(unsigned char)*(img->imageData + i*img->widthStep+3*j+2);
            node.push_back(t); 
        }
    }

    KMeans(node);  

    IplImage *bin=cvCreateImage(cvSize(img->width,img->height),IPL_DEPTH_8U,1);//创建用于显示的图像

    int val=0;
    float step=255/(K-1);
    for(i=0;i<K;i++)
    {
        vector<Tuple> m=clusters[i];
        val=255-i*step;
        for(j=0;j<m.size();j++)
        {
            *(bin->imageData+m[j].y*bin->widthStep+m[j].x)=val;
        }
    }
    cvNamedWindow( "原始图像", 1 ); 
    cvShowImage( "原始图像", img  ); 

     cvNamedWindow( "聚类图像", 1 ); 
     cvShowImage( "聚类图像", bin  ); 
     cvSaveImage("horses2_k4_a.jpg",bin);
     cvWaitKey(0);

     cvDestroyWindow( "原始图像" );
     cvReleaseImage( &img ); 
     cvDestroyWindow( "聚类图像" );
     cvReleaseImage( &bin ); 

     return 0;  
} 

效果