Kmeans聚类算法

K-means算法是最为经典的基于划分的聚类方法 。K-means算法的基本思想是：以空间中k个点为中心进行聚类，对最靠近他们的对象归类。通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。
　　假设要把样本集分为c个类别，算法描述如下：
　　（1）适当选择c个类的初始中心；
　　（2）在第k次迭代中，对任意一个样本，求其到c个中心的距离，将该样本归到距离最短的中心所在的类；
　　（3）利用均值等方法更新该类的中心值；

　　（4）对于所有的c个聚类中心，如果利用(2)(3)的迭代法更新后，值保持不变，则迭代结束，否则继续迭代。

用opencv实现的代码如下：

#include "stdafx.h"  
   
#include <opencv2/opencv.hpp>  
#include "highgui.h"  
#include <math.h>  
typedef  unsigned long uint32;
typedef  unsigned int  uint16;
typedef  unsigned char uint8;

IplImage *src_gray1, *src_gray2, *src_gray3;    
IplImage* src_img;
uint32 dim=0;
CvMat *Label;
uint16 mt[1000][2];

void AllocateImage(IplImage* I)   //给图像分配大小    
{    
    CvSize sz   = cvGetSize(I);    
    
    src_gray1 = cvCreateImage( sz, IPL_DEPTH_8U, 1);    //原图的三个通道    
    src_gray2 = cvCreateImage( sz, IPL_DEPTH_8U, 1);    
    src_gray3 = cvCreateImage( sz, IPL_DEPTH_8U, 1);       
    
}  

void DeallocateImage()
{
	cvReleaseImage(&src_img);  
    cvReleaseImage(&src_gray1);  
    cvReleaseImage(&src_gray2);  
    cvReleaseImage(&src_gray3);  


}

void AllocateMatrix(uint32 num)                //分配矩阵大小
{
	Label         = cvCreateMat( num, 3, CV_16UC1 );
}

void GenerateMatrix()                         //初始化矩阵，矩阵前两列是非零点的坐标，第三列是标签值，0表示未被分类
{
	int i;
	for(i=0; i<Label->rows; i++)
	{
			cvSetReal2D(Label, i, 0, mt[i][0]);	
			cvSetReal2D(Label, i, 1, mt[i][1]);	
			cvSetReal2D(Label, i, 2, 0);	
	}
}

int myabs(int x)                               //取绝对值
{
	if(x<0)
		x=-x;
	return x;
}

void myKmeans(CvMat* L, CvPoint* pt)            //kmeans算法
{
	uint32 PCondition;
	double px,py;
	double MinDis,Dis;
	int i,j;
	double num[4];
	double XValue[4];
	double YValue[4];
	int FLabel;

	CvPoint Last_Kpt[4]={cvPoint(1,1), cvPoint(1,1), cvPoint(1,1), cvPoint(1,1)};
	PCondition = myabs(pt[0].x-Last_Kpt[0].x) + myabs(pt[0].y-Last_Kpt[0].y)
		+ myabs(pt[1].x-Last_Kpt[1].x) + myabs(pt[1].y-Last_Kpt[1].y)
		+ myabs(pt[2].x-Last_Kpt[2].x) + myabs(pt[2].y-Last_Kpt[2].y)
		+ myabs(pt[3].x-Last_Kpt[3].x) + myabs(pt[3].y-Last_Kpt[3].y);

	while( PCondition>2 )
	{
		Last_Kpt[0] = pt[0];
		Last_Kpt[1] = pt[1];
		Last_Kpt[2] = pt[2];
		Last_Kpt[3] = pt[3];
		
		//对每个点贴上标签
		for(i=0; i<Label->rows; i++)                                  
		{	
			MinDis = 10000;
			px = cvGetReal2D(L,i,0);
			py = cvGetReal2D(L,i,1);
			for(j=0;j<4;j++)
			{
				Dis = sqrt( (px-pt[j].x)*(px-pt[j].x) + (py-pt[j].y)*(py-pt[j].y) );
				if(Dis<MinDis)
				{
					MinDis = Dis;
					cvSetReal2D(L,i,2,j+1);
				}
			}
		}
		//变量清零
		for(i=0;i<4;i++)
		{
			XValue[i] = 0;
			YValue[i] = 0;
			num[i] = 0;
		}
		//重新计算中心
		for(i=0; i<Label->rows; i++)                                  
		{
			FLabel = cvGetReal2D(L, i, 2);
			px = cvGetReal2D(L, i, 0);
			py = cvGetReal2D(L, i, 1);

			switch(FLabel)
			{
			case 1:
				XValue[0] +=  px;
				YValue[0] +=  py;
				num[0]++;
				break;
			case 2:
				XValue[1] +=  px;
				YValue[1] +=  py;
				num[1]++;
				break;
			case 3:
				XValue[2] +=  px;
				YValue[2] +=  py;
				num[2]++;
				break;
			case 4:
				XValue[3] +=  px;
				YValue[3] +=  py;
				num[3]++;
				break;
			default:
				break;

			}
		}
		for(i=0;i<4;i++)
		{
			pt[i].x = XValue[i]/num[i];
			pt[i].y = YValue[i]/num[i];
		}
		
		PCondition = myabs(pt[0].x-Last_Kpt[0].x) + myabs(pt[0].y-Last_Kpt[0].y)    //计算终止条件
		+ myabs(pt[1].x-Last_Kpt[1].x) + myabs(pt[1].y-Last_Kpt[1].y)
		+ myabs(pt[2].x-Last_Kpt[2].x) + myabs(pt[2].y-Last_Kpt[2].y)
		+ myabs(pt[3].x-Last_Kpt[3].x) + myabs(pt[3].y-Last_Kpt[3].y);


	}

}

uint32 CalFeatureNum(IplImage* I)           //计算特征点数目，保存坐标
{
	int i,j;
	uint32 num = 0;
	for ( i=0; i<I->height; i++ )
	{
		uint8* ptr = (uint8*)( I->imageData + i*I->widthStep );
		for ( j=0; j<I->width; j++ )
		{
			uint8 Pixel = ptr[j];
			if( Pixel == 255 )
			{
				mt[num][0] = j;
				mt[num][1] = i;
				num++;
			}
			
		}
	}
	return num;
}

void CalRadius(CvMat* L, CvPoint* pt, double* MaxDis)    //计算点群最大半径
{
	int i,FLabel;
	double px,py;
	double dis;
	
	for(i=0; i<Label->rows; i++)                                  
		{
			FLabel = cvGetReal2D(L, i, 2);
			px = cvGetReal2D(L, i, 0);
			py = cvGetReal2D(L, i, 1);

			switch(FLabel)
			{
			case 1:
				dis = sqrt( (px-pt[0].x)*(px-pt[0].x) + (py-pt[0].y)*(py-pt[0].y) );
				if(dis>MaxDis[0])  MaxDis[0] =dis;

				break;
			case 2:
				dis = sqrt( (px-pt[1].x)*(px-pt[1].x) + (py-pt[1].y)*(py-pt[1].y) );
				if(dis>MaxDis[1])  MaxDis[1] =dis;

				break;
			case 3:
				dis = sqrt( (px-pt[2].x)*(px-pt[2].x) + (py-pt[2].y)*(py-pt[2].y) );
				if(dis>MaxDis[2])  MaxDis[2] =dis;

				break;
			case 4:
				dis = sqrt( (px-pt[3].x)*(px-pt[3].x) + (py-pt[3].y)*(py-pt[3].y) );
				if(dis>MaxDis[3])  MaxDis[3] =dis;

				break;
			default:
				break;

			}
		}

}

int _tmain(int argc, _TCHAR* argv[])
{
	CvPoint Kpt[4]={cvPoint(170,96), cvPoint(511,96), cvPoint(170,288), cvPoint(511,288)};      //初始种子点
	double Kmeans_Radius[4] = {0,0,0,0};                                                        //初始半径
	src_img  = cvLoadImage("kmeans_24bit.bmp");
	AllocateImage(src_img);
	cvSplit( src_img, src_gray1, src_gray2, src_gray3, 0);        //将两幅三通道图像分解为3幅单通道图像
	
	dim = CalFeatureNum(src_gray1);                               //计算特征点的个数，并保存特征点的坐标
	AllocateMatrix(dim);										//给矩阵分配大小
	GenerateMatrix();                                            //生成用于Kmeans计算的矩阵
	myKmeans(Label, Kpt);										//Kmeans算法
	CalRadius(Label, Kpt, Kmeans_Radius);                       //计算每个点群的最大半径
	
	cvCircle(src_img,Kpt[0],Kmeans_Radius[0],cvScalar(0,0,255),1,8,0);  //将每个点群标出来
	cvCircle(src_img,Kpt[1],Kmeans_Radius[1],cvScalar(0,0,255),1,8,0);
	cvCircle(src_img,Kpt[2],Kmeans_Radius[2],cvScalar(0,0,255),1,8,0);
	cvCircle(src_img,Kpt[3],Kmeans_Radius[3],cvScalar(0,0,255),1,8,0);

	cvNamedWindow("my picture",CV_WINDOW_AUTOSIZE);  
    
	

	cvShowImage("my picture",src_img);             
   
	

	cvWaitKey(0);  
    DeallocateImage();
	cvDestroyWindow("my picture");  
    
	return 0;
}