/************************************************************************************************/
/*功能：kmeans分割 */
/*參數說明：lpDIBBits－指向DIB位圖數據的指針, lWidth－BMP圖像寬度,lHeight－BMP圖像高度,k－聚類個數*/
void KMeans(BYTE *lpDIBBits, LONG lWidth, LONG lHeight,int k)
{
	unsigned char * lpSrc;
	LONG lLineBytes;	//圖像每行的字節數
	lLineBytes=WIDTHBYTES(lWidth*8);       //計算圖像每行的字節數,假設為256色
	
	int* mark=(int *)HeapAlloc(GetProcessHeap(),0,lWidth*lHeight*sizeof(int));
	float* formalCluster=(float*)HeapAlloc(GetProcessHeap(),0,k);//前次的聚類中心
	float* lastCluster=(float*)HeapAlloc(GetProcessHeap(),0,k);//本次的聚類中心
	const float e=(float)0.1;//誤差
	bool flag=true;//迭代結束標志
	LONG i,j,m; //循環變量

	float sum;
	LONG totalPixel;
	BYTE inteval=255/(k+1);
	//初始的聚類中心
	formalCluster[0]=inteval;
	for(i=1;i<k;i++)
	{
		formalCluster[i]=formalCluster[i-1]+inteval;
	}
	
	//迭代計算聚類中心
	while (flag) 
	{
		//根據前次的聚類中心，把各個象素點分類
		for(m=0;m<k;m++)//每個聚類
		{
			for(i=0;i<lHeight;i++)//每個象素點
				for(j=0;j<lWidth;j++)
				{
					lpSrc=pixelValue(i,j);
					mark[i*lWidth+j]=mostSim(*lpSrc,formalCluster,k);
				}
		}
		//根據分類的象素點，重新計算聚類中心
		
		for(m=0;m<k;m++)//每個聚類
		{
			sum=0;
			totalPixel=0;
			for(i=0;i<lHeight;i++)//每個象素點
				for(j=0;j<lWidth;j++)
				{
					lpSrc=pixelValue(i,j);
					if(mark[i*lWidth+j]==m)//該點屬于這個聚類
					{	
						sum=sum + *lpSrc;
						totalPixel++;
					}

				}
			lastCluster[m]=sum/totalPixel;//用所有屬于這個聚類的點的灰度均值作為新的聚類中心
		}
		//比較原聚類中心和新的聚類中心的誤差
		bool unchanged=true;
		for(m=0;m<k;m++)
		{
			if((float)fabs(formalCluster[m]-lastCluster[m])>e)
			{
				unchanged=false;
				break;
			}
		}
		flag=!unchanged;
		
		//用新的聚類中心代替原聚類中心，再開始下次迭代
		if(flag)
		{
			memcpy(formalCluster,lastCluster,k*sizeof(float));
		}
	}
	//根據劃分好的類別，輸出圖像
		for(i=0;i<lHeight;i++)//每個象素點
			for(j=0;j<lWidth;j++)
			{
				lpSrc=pixelValue(i,j);
				*lpSrc=(unsigned char)lastCluster[mark[i*lWidth+j]];//mark[i][j]就是ij點的類別編號
			}
	//釋放空間
	HeapFree(GetProcessHeap(),0,formalCluster);
	HeapFree(GetProcessHeap(),0,lastCluster);
	HeapFree(GetProcessHeap(),0,mark);
}
/************************************************************************************************/
/*函數功能：找到最相近的類別編號*/
/*參數含義：pixel－象素值, clusterCen－聚類中心, k聚類個數*/
/*函數返回值：返回與象素值pixel最相近的聚類編號*/
int mostSim(BYTE pixel, float *clusterCen, int k)
{
	int i,result=0;
	float* sim=(float *)HeapAlloc(GetProcessHeap(),0,k*sizeof(float));
	//cal the similarities to those clusters
	for(i=0;i<k;i++)
	{
		sim[i]=(float)fabs(clusterCen[i]-pixel*1.0);		
	}
	//find the most sim cluster
	for(i=0;i<k;i++)
	{
		if(sim[i]<sim[result])
			result=i;
	}
	HeapFree(GetProcessHeap(),0,sim);
	return result;
}