isodata.cpp

利用VC++6.0实现ISODATA算法进行聚类分析

// Isodata.cpp : implementation file
//

#include "stdafx.h"
#include "iso.h"
#include "Isodata.h"
#include "math.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CIsodata dialog


CIsodata::CIsodata(CWnd* pParent /*=NULL*/)
	: CDialog(CIsodata::IDD, pParent)
{
	//{{AFX_DATA_INIT(CIsodata)
	m_K = 3;
	m_N = 1;
	m_S = 1.0;
	m_C = 4.0;
	m_L = 0;
	m_I = 20;
	//}}AFX_DATA_INIT
}


void CIsodata::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CIsodata)
	DDX_Text(pDX, IDC_K, m_K);
	DDX_Text(pDX, IDC_N, m_N);
	DDX_Text(pDX, IDC_S, m_S);
	DDX_Text(pDX, IDC_C, m_C);
	DDX_Text(pDX, IDC_L, m_L);
	DDX_Text(pDX, IDC_I, m_I);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CIsodata, CDialog)
	//{{AFX_MSG_MAP(CIsodata)
	ON_BN_CLICKED(IDC_JS, OnJs)
	ON_BN_CLICKED(IDC_LIULAN, OnLiulan)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CIsodata message handlers

void CIsodata::OnJs() 
{
	// TODO: Add your control notification handler code here
	 HWND hwnd;
	 float x[14],y[14];
     int i=0,j=0,m=0;
     FILE *pFile;
     pFile=fopen( strpath, "r" );
     for( i=0;i<14;i++)
	 {
		 fscanf(pFile,"%f\n",&x[i]);
	 }
	 for( i=0;i<14;i++)
	 {
		 fscanf(pFile,"%f\n",&y[i]);
	 }
     
	 
	 float smin[14];   //第j个样本到各个聚类中心的最小距离
	 float rx[5][14];  //第i个聚类中心的第j个样本
     float ry[5][14];
	 int t[5];        //第n个聚类中心的样本个数
	 float zx[5];     //某聚类中样本坐标值的和
     float zy[5];
	 float avx[5];    //某聚类中心值
	 float avy[5];
	 float u[5];     //某聚类中样本到中心的距离之和
	 float avu[5];   //某聚类中样本到中心的距离的平均值
	 float ax[5];    //每一类别中样本与聚类中心距离的标准差向量
     float ay[5];
	 float amax[5];  //每一类中标准差的最大分量
	 float avs;      //所有样本到其相应聚类中心的距离的平均值
	 float k=0.5;
	 float bt[5][5];    //每两个聚类中心之间的距离
	   
     float h=0,q=0;
	 int pp=0;
	 int P=0;
	 int time=0;
	 
	 
	 float s[5][14];    //第j个样本到第i个聚类中心的距离
	 int c = 2;
	  avx[0]=0.0; avy[0]=0.0;
      avx[1]=0.0; avy[1]=1.0;
     
     int K=m_K;
	 int N=m_N;
     float S=m_S;
	 float C=m_C;
	 int L=m_L;
	 int I=m_I;

	 

	 while(time++<I)
	 {
      for( j=0;j<14;j++)            //第三步
	  {

		  smin[j]=20.0;
		  for( i=0;i<c;i++)
		 {
			 
			 s[i][j]=sqrt((x[j]-avx[i])*(x[j]-avx[i])+(y[j]-avy[i])*(y[j]-avy[i]));
             smin[j]=min(smin[j],s[i][j]);
             
		  }
	  }
         for( i=0;i<c;i++)
		 {
             t[i]=0;
			 for(j=0;j<14;j++)
			 {
				
				 if(smin[j]==s[i][j])
				 {
					 rx[i][t[i]]=x[j];
				     ry[i][t[i]]=y[j];
					 t[i]++;
				 }
			 }
			 
		 }
	  
		 BOOL klq1=false;
		 for( i=0;i<c;i++)          //第四步
		 {
			 if(t[i]<N)
			 {
				 t[i]=0;
			     rx[i][t[i]]=rx[i+1][t[i+1]];
				 c--;
				 klq1 = true;
				 
			 }
            
		
		 }
			if(klq1==true)
			{
			 continue;
         	}

		 for( i=0;i<c;i++)        //第五步
		 {
                 zx[i]=0;
                 zy[i]=0;
			 for( m=0;m<t[i];m++)
			 {
                 
				 zx[i]+=rx[i][m];
                 zy[i]+=ry[i][m];
			 }
			 avx[i]=zx[i]/t[i];
             avy[i]=zy[i]/t[i];
		 }
         for( i=0;i<c;i++)      //第六步
		 {
              u[i]=0;
              for( m=0;m<t[i];m++)
			  {
				  u[i]+=sqrt((rx[i][m]-avx[i])*(rx[i][m]-avx[i])+(ry[i][m]-avy[i])*(ry[i][m]-avy[i]));
			  }
			  avu[i]=u[i]/t[i];
		 }
		float he=0;   //第7步
		 int g=0;
		 for( i=0;i<c;i++)
		 {
			  he+=avu[i]*t[i];
			  g+=t[i];
		 }
			  avs=he/g;
		 
			pp = 0;
           for( i=0;i<c;i++)
		   {
			   h = 0.0;
			   q = 0.0;
			   for( m=0;m<t[i];m++)
			   {
				   
				   h+=(rx[i][m]-avx[i])*(rx[i][m]-avx[i]);
				   q+=(ry[i][m]-avy[i])*(ry[i][m]-avy[i]);
			   }
               if(c<=K/2||(c>K/2&&c<2*K&&time%2==1))   //第8步 第2种情况 第3种情况1
			   {
                  ax[i]=sqrt(h/t[i]);       //第9步
			      ay[i]=sqrt(q/t[i]);
			      amax[i]=max(ax[i],ay[i]);  //第10步
			   }

              if(amax[i]>S&&((avu[i]>avs&&t[i]>(2*(N+1)))||c<=K/2)) 
			  {
					pp++;
					P=i;
			  }
				  
		   }
		    
		 if(pp>0) //第11步
		 {
			  
				  c++;
			      if(ax[P]==amax[P])
					   {
						   avx[P]+=k*amax[P];
					       avx[c-1]=avx[P]-amax[P];
						   avy[c-1]=avy[P];
					   }
                       else if(amax[P]==ay[P])
					   {
						   avy[P]+=k*amax[P];
					       avy[c-1]=avy[P]-amax[P];
                           avx[c-1]=avx[P];
					   }
					   continue;
		 }
		   
			   else if(pp==0)              //第12步
			   {
				   float dd[7];
				   int ndd1[7],ndd2[7],countdd=0;
				   for(i=0;i<c-1;i++)
					   for(j=i+1;j<c;j++)
					   {
				          bt[i][j]=sqrt((avx[i]-avx[j])*(avx[i]-avx[j])+(avy[i]-avy[j])*(avy[i]-avy[j]));
                          if(bt[i][j]<C)      //第13步
							   
						  {
								dd[countdd] = bt[i][j];
								ndd1[countdd] = i;
								ndd2[countdd++] = j;
						  }
					   }
					   for (i=1;i<countdd;i++)	//排序
					   {
						   float temp = dd[i];
						   int temp1 = ndd1[i];
						   int temp2 = ndd2[i];
						   j = i-1;
						   while (dd[j]>temp)
						   {
							   dd[j+1] = dd[j];
							   ndd1[j+1] = ndd1[j];
							   ndd2[j+1] = ndd2[j];
							   j--;
						   }
					   }
					   int used[14];
					   int num_used = 0;

					   for(i=0;i<countdd;i++)//合并
					   {
							bool if_used = false;
							for(j=0;j<num_used;j++)
							{
								if(ndd1[i] == used[j] || ndd2[i] == used[j])
									if_used = true;
							}
							if(!if_used)
							{
								avx[ndd1[i]]=(zx[ndd1[i]]+zx[ndd2[i]])/(t[ndd1[i]]+t[ndd2[i]]);   //第14步
								avy[ndd1[i]]=(zy[ndd1[i]]+zy[ndd2[i]])/(t[ndd1[i]]+t[ndd2[i]]);
								avx[ndd2[i]]=avx[c-1];
								avy[ndd2[i]]=avy[c-1];
								t[c-1] = 0;
								used[num_used++] = ndd1[i];
								used[num_used++] = ndd2[i];
								c--;
							}
					   }

					   for( j=0;j<14;j++)            //第三步
	  {

		  smin[j]=20.0;
		  for( i=0;i<c;i++)
		 {
			 
			 s[i][j]=sqrt((x[j]-avx[i])*(x[j]-avx[i])+(y[j]-avy[i])*(y[j]-avy[i]));
             smin[j]=min(smin[j],s[i][j]);
             
		  }
	  }
         for( i=0;i<c;i++)
		 {
             t[i]=0;
			 for(j=0;j<14;j++)
			 {
				
				 if(smin[j]==s[i][j])
				 {
					 rx[i][t[i]]=x[j];
				     ry[i][t[i]]=y[j];
					 t[i]++;
				 }
			 }
			 
		 }
	  
		 BOOL klq1=false;
		 for( i=0;i<c;i++)          //第四步
		 {
			 if(t[i]<N)
			 {
				 t[i]=0;
			     rx[i][t[i]]=rx[i+1][t[i+1]];
				 c--;
				 klq1 = true;
				 
			 }
            
		
		 }
			if(klq1==true)
			{
			 continue;
         	}

		 for( i=0;i<c;i++)        //第五步
		 {
                 zx[i]=0;
                 zy[i]=0;
			 for( m=0;m<t[i];m++)
			 {
                 
				 zx[i]+=rx[i][m];
                 zy[i]+=ry[i][m];
			 }
			 avx[i]=zx[i]/t[i];
             avy[i]=zy[i]/t[i];
		 }
         for( i=0;i<c;i++)      //第六步
		 {
              u[i]=0;
              for( m=0;m<t[i];m++)
			  {
				  u[i]+=sqrt((rx[i][m]-avx[i])*(rx[i][m]-avx[i])+(ry[i][m]-avy[i])*(ry[i][m]-avy[i]));
			  }
			  avu[i]=u[i]/t[i];
		 }
					   
			   }		
	}
	 
	if(time>=I)    //第8步第一种情况
	{
		C=0;
        CString str="ISODATA算法的分类结果.txt";
		FILE *fp;
		fp=fopen("ISODATA算法的分类结果.txt","w+");
		
		for(i=0;i<c;i++)
		{
			fprintf(fp,"聚类中心的坐标为：\n");
			fprintf(fp,"%f\t",avx[i]);
			fprintf(fp,"%f\n",avy[i]);
			fprintf(fp,"第%d类的点的坐标为：\n",i+1);
		   for(j=0;j<t[i];j++)
		   {
			   fprintf(fp,"%f\t",rx[i][j]);
			   fprintf(fp,"%f\n",ry[i][j]);
		   }
           fprintf(fp,"\n");
		}		
		fclose(fp);
		ShellExecute(hwnd,NULL,"ISODATA算法的分类结果.txt",NULL,NULL,SW_SHOWNORMAL);
	}
		 
		 
		
}





 
			 

	


void CIsodata::OnLiulan() 
{
	// TODO: Add your control notification handler code here
    UpdateData(true);
    CFileDialog dlg(true,"txt",".txt");	
	if(dlg.DoModal()==IDOK)
	{
	  strpath= dlg.GetPathName();//获取任意文件路径
	}
	
	UpdateData(false); 
}