用C++实现的ISODATA算法

// 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; } }