isodata.cpp

文件中包含多种模式识别常用的算法

#include<stdio.h>
#include<math.h>

typedef struct sample
{
	float x;
	float y;
} SAMPLE;

float distance(SAMPLE a,SAMPLE b)
{
	float d;
	float a1=a.x, a2=a.y;
	float b1=b.x, b2=b.y;
	d=sqrt((a1-b1)*(a1-b1)+(a2-b2)*(a2-b2));
	return d;
}

int MIN(float D[],int n)
{
	int i;
	int p=0;
	float fmin=D[0];
	for(i=1;i<n;i++)
	{
		if(D[i]<fmin) { fmin=D[i]; p=i; }
	}
	return p;
}


void main()
{
	SAMPLE point[10]={{0,0},{3,8},{2,2},{1,1},{5,3},{4,8},{6,3},{5,4},{6,4},{7,5}}; //模式样本
	SAMPLE center[3]; //聚类中心
	SAMPLE sd[3]; //standard deviation标准差向量
	float sdm[3]; //标准差向量中的最大分量
	int flag[10]; //每个样本所属类别
	int count[3]; //各类别样本数
	float D[3]; //样本到每个聚类中心的距离
	float Da[3]; //各聚类域中诸样本与聚类中心间的平均距离
	float Dl; //全部模式样本对其相应聚类中心的总平均距离
	float Dc[3]; //全部聚类中心的距离
	float a=0.5;
	int t,t1,t2;
	int i,j;
	int n=0; //迭代次数
	int K=3,Nl=1,S=1,C=4,L=1,I=6; //预设参数
	int N=10;
	int Nc=1;
	center[0]=point[0];
	for(i=1;i<3;i++)
	{
		center[i].x=-1;
		center[i].y=-1;
	}

	while(1)
	{
		n++;
		//按近邻原则聚类
		for(i=0;i<N;i++)
		{
			for(j=0;j<Nc;j++)
				D[j]=distance(point[i],center[j]);
			flag[i]=MIN(D,Nc);
		}

		//赋初值
		for(j=0;j<Nc;j++)
		{
			count[j]=0;
			center[j].x=0;
			center[j].y=0;
			Da[j]=0;
			Dl=0;
		}
		//修正聚类中心
		for(i=0;i<N;i++)
		{
			j=flag[i];
			center[j].x=center[j].x+point[i].x;
			center[j].y=center[j].y+point[i].y;
			count[j]++;
		}
		for(j=0;j<Nc;j++)
		{
			center[j].x=center[j].x/count[j];
			center[j].y=center[j].y/count[j];
		}
		//计算各聚类域中诸样本与聚类中心间的平均距离
		//和全部模式样本对其相应聚类中心的总平均距离
		for(i=0;i<N;i++)
		{
			j=flag[i];
			Da[j]+=distance(point[i],center[j]);
		}
		for(j=0;j<Nc;j++)
		{
			Dl+=Da[j];
			Da[j]=Da[j]/count[j];
		}
		Dl=Dl/N;
		
		//进行分裂处理
		if(n!=I&&(Nc<=K/2||(Nc>K/2&&((fmod(n,2)!=0)&&Nc<2*K))))
		{
			for(j=0;j<Nc;j++)
			{
				sd[j].x=0;
				sd[j].y=0;
			}
			for(i=0;i<N;i++)
			{
				j=flag[i];
				sd[j].x+=(point[i].x-center[j].x)*(point[i].x-center[j].x);
				sd[j].y+=(point[i].y-center[j].y)*(point[i].y-center[j].y);
			}
			t=Nc;
			for(j=0;j<t;j++)
			{
				sd[j].x=sqrt(sd[j].x/count[j]);
				sd[j].y=sqrt(sd[j].y/count[j]);
				if(sd[j].x>=sd[j].y)
				{
					sdm[j]=sd[j].x;
					if(sdm[j]>S&&((Da[j]>Dl&&count[j]>2*(Nl+1))||Nc<=K/2))
					{
						center[Nc].x=center[j].x-a*sdm[j];
						center[Nc].y=center[j].y;
						center[j].x=center[j].x+a*sdm[j];
						Nc++;
					}
				}
				else
				{
					sdm[j]=sd[j].y;
					if(sdm[j]>S&&((Da[j]>Dl&&count[j]>2*(Nl+1))||Nc<=K/2))
					{
						center[Nc].x=center[j].x;
						center[Nc].y=center[j].y-a*sdm[j];
						center[j].y=center[j].y+a*sdm[j];
						Nc++;
					}
				}
			}
			if(Nc!=t) //进行了分裂处理
				continue;
		}

		if(n==I) C=0; //最后一次迭代

		//进行合并处理
		t=0;
		t1=-1;
		t2=-1;
		for(i=1;i<Nc;i++)
		{
			for(j=0;j<i;j++)
			{
				Dc[t]=distance(center[i],center[j]);
				if(t>0&&Dc[t]<C&&Dc[t-1]>Dc[t])
					Dc[t]=Dc[t-1];
				else if(Dc[t]<C) { t1=j; t2=i; }
				t++;
			}
		}
		if(t1!=-1&&t2!=-1)
		{
			center[t1].x=(count[t1]*center[t1].x+count[t2]*center[t2].x)/(count[t1]+count[t2]);
			center[t1].y=(count[t1]*center[t1].y+count[t2]*center[t2].y)/(count[t1]+count[t2]);
			center[t2].x=-1;
			center[t2].y=-1;
			if(center[Nc-1].x!=-1&&center[Nc-1].y!=-1)
			{
				center[t2].x=center[Nc-1].x;
				center[t2].y=center[Nc-1].y;
			}
			Nc--;
		}

		if(n==I) break;
	}

	printf("the first class:\n");
	for(i=0;i<10;i++)
	{
		if(flag[i]==0)
			printf("point[%d]=(%f,%f)\n",i,point[i].x,point[i].y);
	}
	printf("the class center is (%f,%f)\n",center[0].x,center[0].y);
	printf("\n");
	printf("the secend class:\n");
	for(i=0;i<10;i++)
	{
		if(flag[i]==1)
			printf("point[%d]=(%f,%f)\n",i,point[i].x,point[i].y);
	}
	printf("the class center is (%f,%f)\n",center[1].x,center[1].y);
	printf("\n");
	printf("the third class:\n");
	for(i=0;i<10;i++)
	{
		if(flag[i]==2)
			printf("point[%d]=(%f,%f)\n",i,point[i].x,point[i].y);
	}
	printf("the class center is (%f,%f)\n",center[2].x,center[2].y);
}