//#include "stdafx.h"
//#include "ImageObj.h"
#include "Clustering.h"
#include "math.h"

/******************************************************************************
/*	Name:		LBGCluster
/*	Function:	Clustering input vectors using LBG algorithm
/*				Using Euclidean distance
/*	Parameter:	X -- Input vecters
/*				N -- Number of input vectors
/*				Y -- Clustering result
/*				M -- Number of clustering center
/*	Return:		0 -- Correct
/*				1 -- Error
/*
/******************************************************************************/
int LBGCluster(VQ_VECTOR *X, int N, VQ_CENTER *Y, int M)//x分配多大空間，為何不初始化指針？
{
	if(N<M)			 return -1;
	int		L=1000, m=1, nCenter, i, j, k;
	int		nDimension = X[0].nDimension;//X[0].nDimension仍就不確定？
	double  D0, D;
	struct  VQ_CENTERINFO
	{
		double*	  Data;//聚類的形心
		int		  nDimension;//矢量維數
		double*   SumData;//屬于此類的矢量各維之和
		int		  Num;//屬于此類的矢量個數
	};
	VQ_CENTERINFO	*Center = (VQ_CENTERINFO*)malloc(M*sizeof(VQ_CENTERINFO));
	if(Center == NULL)		return -1;
	double			*Distance = (double*)malloc(N*sizeof(double));
	if(Distance == NULL)	return -1;

	for( i=0; i<M; i++)
	{
		Center[i].nDimension = nDimension;
		Center[i].Data = (double*)malloc(sizeof(double)*nDimension);
		Center[i].SumData = (double*)malloc(sizeof(double)*nDimension);
		if( Center[i].Data == NULL || Center[i].SumData == NULL )
		{
			AfxMessageBox( "Memory used up!" );
			return -1;
		}
		for( j=0; j<nDimension; j++ )
		{
			Center[i].Data[j] = X[i*N/M].Data[j];//形成初始碼書
			Center[i].SumData[j] = 0;//？
		}
		Center[i].Num = 0;//？
	}

	D0=1;         D=1e+10;	
	while(m<L && fabs(D0-D)/D0>1e-5)//m=1 l=1000 平均失真？
	{
		for(i=0; i<M; i++)
		{
			for( j=0; j<nDimension; j++ )
				Center[i].SumData[j] = 0;
			Center[i].Num = 0;
		}
		D0 = D;			D = 0;		m++;//m控制迭代次數，不超過1000?
		for(i=0; i<N; i++)
		{
			Distance[i] = 1e+10;
			for(int j=0; j<M; j++)
			{
				double  Dist = 0;
				for( k=0; k<nDimension; k++ )
					Dist += (X[i].Data[k]-Center[j].Data[k])*(X[i].Data[k]-Center[j].Data[k]);
				if( Dist < Distance[i])
				{			
					nCenter = j;//?
					Distance[i] = Dist;
				}
			}
			X[i].nCluster = nCenter;//以上為確定每一個矢量屬于哪個胞腔。
			for( k=0; k<nDimension; k++ )
				Center[nCenter].SumData[k] += X[i].Data[k];
			Center[nCenter].Num++;
			D += Distance[i];//計算平均失真，為何不除以M
		}
		for(i=0; i<M; i++)
		{
			if(Center[i].Num != 0)//確定不是空胞腔
				for( k=0; k<nDimension; k++ )
					Center[i].Data[k] = Center[i].SumData[k]/Center[i].Num;//重新計算形心
			else
			{	
				int MaxNum=0;
				for( k=1; k<M; k++)
					MaxNum = Center[i].Num > Center[MaxNum].Num ? i: MaxNum;
				int   Num = Center[MaxNum].Num/2;
				for( k=0; k<nDimension; k++ )
					Center[MaxNum].SumData[k] = 0;
				Center[MaxNum].Num = 0;
				for(k=0; k<N; k++)
				{	
					if(X[k].nCluster != MaxNum)		continue;
					if(Center[i].Num < Num)
					{   
						X[k].nCluster = i;
						for( m=0; m<nDimension; m++)
							Center[i].SumData[m] += X[k].Data[m];
						Center[i].Num++;
					}
					else
					{
						for( m=0; m<nDimension; m++ )
							Center[MaxNum].SumData[m] += X[k].Data[m];
						Center[MaxNum].Num++;
					}
				}
				for( m=0; m<nDimension; m++ )
					Center[i].Data[m] = Center[i].SumData[m] / Center[i].Num;
				if(MaxNum < i)
					for( m=0; m<nDimension; m++ )
						Center[MaxNum].Data[m] = Center[MaxNum].SumData[m] / Center[MaxNum].Num;
			}
		}
	}
	for(i=0; i<M; i++)
	{
		for( m=0; m<nDimension; m++ )
			Y[i].Data[m] = Center[i].Data[m];
		Y[i].Num = Center[i].Num;
	}
	for( i=0; i<M; i++ )
	{
		free( Center[i].Data );
		free( Center[i].SumData );
	}
	free(Center);
	free(Distance);
	return 0;
}

/******************************************************************************
/*	Name:		DumpClusterData
/*	Function:	Dump clustering result to a text file for debugging
/*	Parameter:	FileName -- Dump text file name
/*				X -- Input vecters
/*				N -- Number of input vectors
/*				Y -- Clustering result
/*				M -- Number of clustering center
/*	Return:		0 -- Correct
/*				1 -- Error
/*
/******************************************************************************/
int DumpClusterData(CString FileName, VQ_VECTOR *X, int N, VQ_CENTER *Y, int M)
{
	int		i, j, k;
	int		nDimension = X[0].nDimension;
	FILE *fp = fopen(FileName, "wt");
	for( i=0; i<M; i++)
	{
		fprintf(fp, "Center%02d: ", i);
		for(  k=0; k<nDimension; k++ )
			fprintf( fp, "%5.1f ", Y[i].Data[k] );
		fprintf( fp, "Num=%03d\n", Y[i].Num);
		for( j=0; j<N; j++)
		{	if(X[j].nCluster != i)		continue;
			double Distance = 0;
			for( k=0; k<nDimension; k++)
				Distance += (X[j].Data[k]- Y[i].Data[k])*(X[j].Data[k]-Y[i].Data[k]);
			Distance = sqrt( Distance );
			for( k=0; k<nDimension; k++ )
				fprintf(fp, "  %03d  ", (int)X[j].Data[k] );
			fprintf( fp, " D=%5.1f\n", Distance);
		}
	}
	fclose(fp);
	return 0;
}