#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;  
} 

效果