// ISODATA.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#define  N 10
#define  eps   0.00001

struct Pointf
{
	int sequence;
	float x;
	float y;
};

struct PointZ
{
	float x;
	float y;
};

float CalDistancef(Pointf x1,Pointf x2)
{
	return sqrtf((x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y));
}

float CalDistanceZ(PointZ x1,PointZ x2)
{
	return sqrtf((x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y));
}

float CalDistancefZ(Pointf x1,PointZ x2)
{
	return sqrtf((x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y));
}




int main(int argc, char* argv[])
{
	Pointf pts[N]={
		{0,0.0,0.0},{1,3.0,8.0},{2,2.0,2.0},{3,1.0,1.0},{4,5.0,3.0},
		{5,4.0,8.0},{6,6.0,3.0},{7,5.0,4.0},{8,6.0,4.0},{9,7.0,5.0},
	};
	int i,j,m;

	printf("樣本集為:\n");
	
	for(i=0;i<N;i++)
	{
		printf("X%d(%.1f,%.1f)  ",i,pts[i].x,pts[i].y);
		
		if((i+1)%5==0)
		{
			printf("\n");
		}
		
	}
	printf("\n");
	printf("\n");
	
	int Nc=0;
	printf("please input Nc(0-10): ");
	scanf("%d",&Nc);

	int Z[N];

	for(i=0;i<Nc;i++)
	{
		printf("輸入初始第%d聚類中心的序號(0-9):",i);
		scanf("%d",&Z[i]);
	}
	
    int    Nj[N]; //記錄每個類中元素的個數
	PointZ ZArray[N];
	Pointf SAArray[N][N];
	float  DjAv[N];
	float  Deltaj[N][2];
	float  Deltajmax[N];
	int  DeltajmaxCor[N];
	float  DAv;
	int    Nreal=N;
	int    count=0;

	float  Dij[N*N/2];
	int    Diji[N];
	int    Dijj[N];
	int    q=0;
	int    p=0;
    float ft;
	int   it;
	int   jt;
	int   flag;
	int ss=0;
	PointZ Ztp;
	PointZ ZArraytp[N];
	int    Nctp;

	char ch;
	
	int cur=0;

	for(i=0;i<N;i++)
	{
		Nj[i]=0; 
	}

	//聚類中心的特征值
	for(i=0;i<Nc;i++)
	{
		int ihere=Z[i];
		ZArray[i].x=pts[ihere].x;
		ZArray[i].y=pts[ihere].y;
	}

	int K,ThetaN;
	float ThetaS,ThetaC;
	int L,I;

Step1:
	printf("輸入預期聚類中心數目 K :");
	scanf("%d",&K);
	printf("輸入每個聚類域中最少的樣本數ThetaN: ");
	scanf("%d",&ThetaN);

	

	printf("輸入同一聚類域中樣本標準差的最大值: ");
	scanf("%f",&ThetaS);

	printf("輸入不同聚類域距離最小值: ");
	scanf("%f",&ThetaC);

	printf("輸入一次可以合并的聚類中心的最多對數: ");
	scanf("%d",&L);

	printf("輸入最大迭代次數: ");
	scanf("%d",&I);
Step2:
	for(i=0;i<Nc;i++)
	{
		Nj[i]=0;
	}

		printf("\n");
		printf("這是第%d次歸類\n",count+1);
		
	for(i=0;i<N;i++)                          //將模式樣本歸類
	{
		if(pts[i].sequence==-1)continue;      //若該點的序號為－1則說明它是被剔除的
		float dis=1.0e+10;
		int xx=0;
		float ftemp;

		for(j=0;j<Nc;j++)                     
		{
			ftemp=CalDistancefZ(pts[i],ZArray[j]);
			if(ftemp<dis||fabs(dis-ftemp)<eps)
			{
				xx=j;
				dis=ftemp;
			}
		}
				
		SAArray[xx][Nj[xx]].x=pts[i].x;
		SAArray[xx][Nj[xx]].y=pts[i].y;
		SAArray[xx][Nj[xx]].sequence=pts[i].sequence;	
		Nj[xx]=Nj[xx]+1;	
	}

	   for(i=0;i<Nc;i++)
	   {
		   printf("第%d個聚類中心是:(%.2f,%.2f)   ",i,ZArray[i].x,ZArray[i].y);
		   printf("包含的元素有:",i);
		   for(j=0;j<Nj[i];j++)
		   {
			   printf(" X%d ",SAArray[i][j].sequence);
		   }
		   printf("\n");
		  // printf("Nj(%d) is %d\n",i,Nj[i]);
	   }

	count++;

Step3:

	for(j=0;j<Nc;j++)               //是否可以去掉一些數據
	{
		if(Nj[j]<ThetaN)
		{
			for(i=0;i<Nj[j];i++)
			{
				pts[SAArray[j][i].sequence].sequence=-1;
			}
			i=j;
			int tr=j;
			Nreal-=Nj[j];
			while(j<Nc-1)
			{
				
				for(m=0;m<Nj[j+1];m++)
				{
					SAArray[j][m].x=SAArray[j+1][m].x;
					SAArray[j][m].y=SAArray[j+1][m].y;
					SAArray[j][m].sequence=SAArray[j+1][m].sequence;
				} 
				j++;			
			}

			while(i<Nc-1)
			{
				Nj[i]=Nj[i+1];
				i++;
			}
			Nc--;
			j=tr;
		}		
	}

Step4:              //修正各聚類中心
	for(j=0;j<Nc;j++)
	{
		float temx=0,temy=0;
		for(i=0;i<Nj[j];i++)
		{
			temx+=SAArray[j][i].x;
			temy+=SAArray[j][i].y;		
		}
		ZArray[j].x=temx/Nj[j];
		ZArray[j].y=temy/Nj[j];
	}

Step5://計算各聚類域中諸樣本與聚類中心的平均距離
	float temp=0.0;
	for(j=0;j<Nc;j++)
	{
		for(i=0;i<Nj[j];i++)
		{
			temp+=CalDistancefZ(SAArray[j][i],ZArray[j]);
		}

		DjAv[j]=temp/Nj[j];
		
		temp=0.0;
	}
Step6://計算全部模式樣本對應聚類中心的總平均距離
	DAv=0;
	for(j=0;j<Nc;j++)
	{
		DAv+=Nj[j]*DjAv[j];
	}
	DAv/=Nreal;

Step7:

	if(count>=I) goto Step14;

	if(Nc<=K/2)goto Step8;

	if((count%2==0)||Nc>=2*K) goto Step11;
Step8://計算各聚類中樣本距離標準差
	
	for(j=0;j<Nc;j++)
	{
		float temx=0.0,temy=0.0;
		for(i=0;i<Nj[j];i++)
		{
           temx+=(SAArray[j][i].x-ZArray[j].x)*(SAArray[j][i].x-ZArray[j].x);
		   temy+=(SAArray[j][i].y-ZArray[j].y)*(SAArray[j][i].y-ZArray[j].y);
		}
        Deltaj[j][0]=sqrtf(temx/Nj[j]);
		Deltaj[j][1]=sqrtf(temy/Nj[j]);	
		temx=0.0,temy=0.0;	
	}	
Step9://求每個標準差向量中的最大分量 
	for(j=0;j<Nc;j++)
	{
		Deltajmax[j]=Deltaj[j][0]>Deltaj[j][1]?Deltaj[j][0]:Deltaj[j][1];
		DeltajmaxCor[j]=Deltaj[j][0]>Deltaj[j][1]?0:1;
	}
Step10://分裂判斷和計算
	for(j=0;j<Nc;j++)
	{
		if(Deltajmax[j]>ThetaS)
		{
			if((DjAv[j]>DAv&&Nj[j]>2*(ThetaN+1))||Nc<=K/2)
			{
				float Garma=0.5;
				PointZ Zj1,Zj2;

				if(DeltajmaxCor[j]==0)
				{
					Zj1.x=ZArray[j].x+Deltajmax[j]*Garma;
					Zj1.y=ZArray[j].y;
					Zj2.x=ZArray[j].x-Deltajmax[j]*Garma;
					Zj2.y=ZArray[j].y;
				}
				else if(DeltajmaxCor[j]==1)
				{
					Zj1.x=ZArray[j].x;
					Zj1.y=ZArray[j].y+Deltajmax[j]*Garma;
					Zj2.x=ZArray[j].x;
					Zj2.y=ZArray[j].y-Deltajmax[j]*Garma;
				}
				
                ZArray[j].x=Zj1.x;
				ZArray[j].y=Zj1.y;

				ZArray[Nc].x=Zj2.x;
				ZArray[Nc].y=Zj2.y;
				Nc++;
				goto Step2;
			}
		}
	}

Step11://計算全部聚類中心的距離
	ss=0;
	for(i=0;i<Nc-1;i++)
	{
		for(j=i+1;j<Nc;j++)
		{
			Dij[ss]=CalDistanceZ(ZArray[i],ZArray[j]);
			Diji[ss]=i;
			Dijj[ss]=j;
			ss++;		
		}
	}
		
Step12: //簡單起見，只考慮一次只合并一對聚類中心的情況
	    //找出類間距離最小的
	ft=Dij[0];
	it=Diji[0];
	jt=Dijj[0];
	for(i=1;i<ss;i++)
	{
		if(Dij[i]<ft)
		{
			ft=Dij[i];
			it=Diji[i];
			jt=Dijj[i];
		}
	}
Step13:
	if(ft<ThetaC)
	{
		Ztp.x=(Nj[it]*ZArray[it].x+Nj[jt]*ZArray[jt].x)/(Nj[it]+Nj[jt]);
		Ztp.y=(Nj[it]*ZArray[it].y+Nj[jt]*ZArray[jt].y)/(Nj[it]+Nj[jt]);
		ZArray[it].x=Ztp.x;
		ZArray[jt].y=Ztp.y;

		j=jt;
		while(j<Nc-1)
		{
			ZArray[j].x=ZArray[j+1].x;
			ZArray[j].y=ZArray[j+1].y;
			j++;
		}
		j=jt;
		
		while(j<Nc-1)
		{
			Nj[j]=Nj[j+1];
			j++;
		}
		Nc--;
	}
Step14:
	if(count>=I) 
	{
		count=0;	
		printf("\n");
		printf("共分為%d類\n",Nc);
			return 0;
		
	}
	else
	{
		goto Step2;	
	}

}