isodataclassifier.cpp

这是我利用STL写的ISODATA分类算法

// ISODATAClassifier.cpp: implementation of the CISODATAClassifier class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ISODATAClassifier.h"
#include "math.h"
#include "algorithm"
#include "resource.h"
#include "dlgisodataopt.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CISODATAClassifier::CISODATAClassifier()
{
	m_pdSwatchData = NULL;
	m_iPattDimen = 0;
	m_iSwatchNum = 0;
}

CISODATAClassifier::~CISODATAClassifier()
{
	delete []m_pdSwatchData;
}

int CISODATAClassifier::InitSwatchData(double* pSwatchData, int iSwatchNum, int iPattDimen)
{
	m_iSwatchNum = iSwatchNum;
	m_iPattDimen = iPattDimen;
	if(m_pdSwatchData != NULL)
		delete []m_pdSwatchData;
	m_pdSwatchData = new double[iPattDimen*iSwatchNum];
	memcpy(m_pdSwatchData, pSwatchData, sizeof(double)*iPattDimen*iSwatchNum);
	return 1;
}

int CISODATAClassifier::ClassifyProcess()
{
	if(NULL == m_pdSwatchData)
		return 0;
	CTypeClass FirClass(2);
	//初始聚类中心
	for(int i=0; i<m_iPattDimen; i++)
		FirClass.m_pdCenter[i] = m_pdSwatchData[i];

	m_VETypeClass.push_back(FirClass);
	double dMinDis;
	int iClassType;
	int n = 0;
	
	do 
	{
		n++;
		CDlgISODATAOpt DlgISODATAOpt;
		DlgISODATAOpt.SetISODATAOpt(m_ISODATAOpt, n);
		if(DlgISODATAOpt.DoModal() == IDOK)
		{
			DlgISODATAOpt.GetISODATAOpt(m_ISODATAOpt);
		}

		//第三步：最小距离分类
		for(int i=0; i<m_VETypeClass.size(); i++)
			m_VETypeClass[i].m_VESwatch.clear();
		for(i=0; i<m_iSwatchNum; i++)		
		{
			dMinDis = 0;

			for(int j=0; j<m_VETypeClass.size(); j++)
			{
				double dTmpDis = 0;
				for(int k=0; k<m_iPattDimen; k++)
					dTmpDis += (m_pdSwatchData[i*m_iPattDimen+k]-m_VETypeClass[j].m_pdCenter[k])*
							   (m_pdSwatchData[i*m_iPattDimen+k]-m_VETypeClass[j].m_pdCenter[k]);
				if(j == 0)
				{
					dMinDis = dTmpDis;
					iClassType = j;
					continue;
				}
				if(dTmpDis < dMinDis)
				{
					dMinDis = dTmpDis;
					iClassType = j;
				}
			}
			m_VETypeClass[iClassType].m_VESwatch.push_back(i);
		}
		//第四步：取消样本数目不足的类
		for(int j=0; j<m_VETypeClass.size(); j++)
			if(m_VETypeClass[j].m_VESwatch.size() < m_ISODATAOpt.iClassMinNum)
				m_VETypeClass.erase(m_VETypeClass.begin() + j);
		//第五步：修正各聚类中心值
		for(i=0; i<m_VETypeClass.size(); i++)
		{
			for(int k=0; k<m_iPattDimen; k++)
			{	
				double dMean = 0;
				for(int j=0; j<m_VETypeClass[i].m_VESwatch.size(); j++)
					dMean += m_pdSwatchData[m_VETypeClass[i].m_VESwatch[j]*m_iPattDimen+k];
				m_VETypeClass[i].m_pdCenter[k] = dMean/m_VETypeClass[i].m_VESwatch.size();
			}
		}	
		//第六步：所有样本到中心的距离平均值
		for(i=0; i<m_VETypeClass.size(); i++)
		{
			m_VETypeClass[i].m_dSCDist = 0;

			double dSCDis = 0;
			for(int j=0; j<m_VETypeClass[i].m_VESwatch.size(); j++)
			{	  
				for(int k=0; k<m_iPattDimen; k++)
				{	
					dSCDis += (m_pdSwatchData[m_VETypeClass[i].m_VESwatch[j]*m_iPattDimen+k]-m_VETypeClass[i].m_pdCenter[k])*
							  (m_pdSwatchData[m_VETypeClass[i].m_VESwatch[j]*m_iPattDimen+k]-m_VETypeClass[i].m_pdCenter[k]);

					
				//	m_VETypeClass[i].m_pdCenter[k] = dSCDis/m_VETypeClass[i].m_VESwatch.size();
				}
				dSCDis = sqrt(dSCDis);
				m_VETypeClass[i].m_dSCDist += dSCDis;
			}
			m_VETypeClass[i].m_dSCDist = m_VETypeClass[i].m_dSCDist/m_VETypeClass[i].m_VESwatch.size();
		}
		//第七步：所有样本到中心的平均值
		double dMeanDis = 0;
		for(i=0; i<m_VETypeClass.size(); i++)
			dMeanDis += m_VETypeClass[i].m_dSCDist*m_VETypeClass[i].m_VESwatch.size();
		dMeanDis = dMeanDis/m_iSwatchNum;
		//第八步：选择
	//	if(m_VETypeClass.size() <= m_ISODATAOpt.iDemaK/2)
		int iChoice = 0;
		if(m_VETypeClass.size() <= m_ISODATAOpt.iDemaK/2)
			iChoice = 1;
		else if(n%2 == 0 || m_VETypeClass.size() >= m_ISODATAOpt.iDemaK*2 || n == m_ISODATAOpt.iIterMaxNum)
			iChoice = 2;
		else
			iChoice = 1;
		if(1 == iChoice)
		{	
	    //第九步：标准差计算
			for(int i=0; i<m_VETypeClass.size(); i++)
			{
				m_VETypeClass[i].m_iMaxVari = 0;
				for(int j=0; j<m_iPattDimen; j++)
				{
					for(int k=0; k<m_VETypeClass[i].m_VESwatch.size(); k++)
					{
						m_VETypeClass[i].m_pdVariance[j] += 
							(m_pdSwatchData[m_VETypeClass[i].m_VESwatch[k]*m_iPattDimen+j]-m_VETypeClass[i].m_pdCenter[j])*
							(m_pdSwatchData[m_VETypeClass[i].m_VESwatch[k]*m_iPattDimen+j]-m_VETypeClass[i].m_pdCenter[j]);
						
					}
					m_VETypeClass[i].m_pdVariance[j] = sqrt(m_VETypeClass[i].m_pdVariance[j]/m_VETypeClass[i].m_VESwatch.size() );
					//第十步：最大标准差分量
					if(m_VETypeClass[i].m_pdVariance[j] > m_VETypeClass[i].m_iMaxVari)
						m_VETypeClass[i].m_iMaxVari = j;
				}
			}
		//第十步：		
			int iOriNum = m_VETypeClass.size();
			int* pFlag = new int[iOriNum];
			memset(pFlag, 0, sizeof(int));
			for(i=0; i<iOriNum; i++)
			{
				if(m_VETypeClass[i].m_pdVariance[m_VETypeClass[i].m_iMaxVari] > m_ISODATAOpt.dSigThres)
				{
					if( (m_VETypeClass[i].m_dSCDist > dMeanDis && m_VETypeClass[i].m_VESwatch.size() > 2*(m_ISODATAOpt.iClassMinNum+1)) ||
					   m_VETypeClass.size() <= m_ISODATAOpt.iDemaK/2)
					{
						CTypeClass ClassZAdd(2), ClassZMinu(2);
						ClassZAdd = m_VETypeClass[i];
						ClassZMinu = m_VETypeClass[i];
						ClassZAdd.m_pdCenter[m_VETypeClass[i].m_iMaxVari] += m_VETypeClass[i].m_pdVariance[m_VETypeClass[i].m_iMaxVari]*m_ISODATAOpt.dSigInter;
						ClassZMinu.m_pdCenter[m_VETypeClass[i].m_iMaxVari] -= m_VETypeClass[i].m_pdVariance[m_VETypeClass[i].m_iMaxVari]*m_ISODATAOpt.dSigInter;
						//m_VETypeClass.erase(m_VETypeClass.begin()+i);
						pFlag[i] = 1;
						m_VETypeClass.push_back(ClassZAdd);
						m_VETypeClass.push_back(ClassZMinu);
					}
	
				}
			}
			for(i=0; i<iOriNum; i++)
			{
				if(pFlag[i] == 1)
					m_VETypeClass.erase(m_VETypeClass.begin() + i);
			}	
			delete[] pFlag;
		}	
		//第十三步
		else if(2 == iChoice)
		{	
			if(n != m_ISODATAOpt.iIterMaxNum)
			{
				vector<MERGE> VEMerge;
				int iOriNum = m_VETypeClass.size();
				int* pFlag = new int[iOriNum];
				memset(pFlag, 0, sizeof(double));
				//计算类间距离	
				for(int i=0; i<m_VETypeClass.size(); i++)
				{
					for(int j=i+1; j<m_VETypeClass.size(); j++)
					{
						double pDis = 0;
						for(int k=0; k<m_iPattDimen; k++)
						{
							pDis += (m_VETypeClass[i].m_pdCenter[k]-m_VETypeClass[j].m_pdCenter[k])*
									(m_VETypeClass[i].m_pdCenter[k]-m_VETypeClass[j].m_pdCenter[k]);
						}
						if(sqrt(pDis) > m_ISODATAOpt.dCenThres)
						{
							MERGE Merge;
							Merge.iFir = i;
							Merge.iSec = j;
							Merge.dDis = pDis;
							VEMerge.push_back(Merge);
						}
						
					}
				}
				//类间距离排序
				sort(VEMerge.begin(), VEMerge.end(), MergeCompare);
				int iEnd = VEMerge.size();
				if(iEnd > m_ISODATAOpt.iMergMaxNum)
					iEnd = m_ISODATAOpt.iMergMaxNum;
				for(i=0; i<iEnd; i++)
				{
					if(pFlag[VEMerge[i].iSec] != 1 && pFlag[VEMerge[i].iFir] != 1)
					{
						CTypeClass MerClass(2);
						for(int j=0; j<m_iPattDimen; j++)
						{
							MerClass.m_pdCenter[j] =(m_VETypeClass[VEMerge[i].iFir].m_VESwatch.size()*m_VETypeClass[VEMerge[i].iFir].m_pdCenter[j] +
												     m_VETypeClass[VEMerge[i].iFir].m_VESwatch.size()*m_VETypeClass[VEMerge[i].iSec].m_pdCenter[j])/
													 (m_VETypeClass[VEMerge[i].iFir].m_VESwatch.size()+m_VETypeClass[VEMerge[i].iSec].m_VESwatch.size());
						}
						m_VETypeClass.push_back(MerClass);
					}
				}
				for(i=0; i<iOriNum; i++)
				{
					if(pFlag[i] == 1)
						m_VETypeClass.erase(m_VETypeClass.begin() + i);
				}

				delete[] pFlag;
			}
		}
		if(m_ISODATAOpt.bNoExamMidFile == FALSE)
		{
			char strPath[100];
			::GetModuleFileName(NULL, strPath, 100);
			CString strOutFile;
			strOutFile.Format("%s", strPath);
			strOutFile.Delete(strOutFile.GetLength()-17, 17);
			CString strDataPath("Data\\中间状态文件.txt");
			strOutFile += strDataPath;
			char* strFile = strOutFile.LockBuffer();
			strOutFile.UnlockBuffer();

			ShowRecentCondition(m_VETypeClass, strFile);
			::ShellExecute(NULL, "open" ,strFile, "", "", SW_SHOW);
		}
	}while(n != m_ISODATAOpt.iIterMaxNum);

	return 1;
}

BOOL CISODATAClassifier::MergeCompare(const MERGE& MergeL, const MERGE& MergeR)
{
	if(MergeL.dDis > MergeR.dDis)
		return FALSE;
	return TRUE;
}

int CISODATAClassifier::ResultFile(char* strFile)
{
	FILE* fp = fopen(strFile, "wt");
	fprintf(fp, "			             结果输出文件\n");
	fprintf(fp, "//////////////////////////////////////////////////////////////////////////\n");
	fprintf(fp, "//                  作者:          学号                                 //\n");
	fprintf(fp, "//                  林立文     200432590138                             //\n");			
	fprintf(fp, "//                  研发单位：武汉大学遥感学院04023                     //\n");
	fprintf(fp, "//                  版权所有    翻版必究                                //\n");
	fprintf(fp, "//////////////////////////////////////////////////////////////////////////\n\n");
	fprintf(fp, "                         文件基本信息\n\n");
	fprintf(fp, "总样本数:%d              总分类数:%d\n", m_iSwatchNum, m_VETypeClass.size());

	fprintf(fp, "///////////////////////////////////////////////////////\n");

	for(int i=0; i<m_VETypeClass.size(); i++)
	{
		fprintf(fp, "第%d类样本信息：\n", i+1);
		fprintf(fp, "样本中心：");
		for(int k=0; k<m_iPattDimen; k++)
		{
			fprintf(fp, "%4.2lf ", m_VETypeClass[i].m_pdCenter[k]);
		}

		fprintf(fp, "\n");
		fprintf(fp, "样本数：%d\n", m_VETypeClass[i].m_VESwatch.size());
		fprintf(fp, "样本号:样本特征向量  ");
		for(int j=0; j<m_VETypeClass[i].m_VESwatch.size(); j++)
		{
			fprintf(fp, " %d: ", m_VETypeClass[i].m_VESwatch[j]+1);
			for(int k=0; k<m_iPattDimen; k++)
			{
				fprintf(fp, "%4.2lf ", m_pdSwatchData[m_VETypeClass[i].m_VESwatch[j]*m_iPattDimen+k]);
			}
			fprintf(fp, "\n                     ");
		}
		fprintf(fp, "\n");
		fprintf(fp, "///////////////////////////////////////////////////////\n");
	}
	fclose(fp);
	return 1;
}

int CISODATAClassifier::ShowRecentCondition(vector<CTypeClass> VETypeClass, char* strFile)
{
	FILE* fp = fopen(strFile, "wt");

	fprintf(fp, "总样本数:%d             当前聚类数:%d\n", m_iSwatchNum, m_VETypeClass.size());

	fprintf(fp, "///////////////////////////////////////////////////////\n");

	for(int i=0; i<m_VETypeClass.size(); i++)
	{
		fprintf(fp, "第%d类样本信息：\n", i+1);
		fprintf(fp, "类中心：");
		for(int k=0; k<m_iPattDimen; k++)
		{
			fprintf(fp, "%4.2lf ", m_VETypeClass[i].m_pdCenter[k]);
		}
		fprintf(fp, "\n");

		fprintf(fp, "类标准方差：");
		for(k=0; k<m_iPattDimen; k++)
		{
			fprintf(fp, "%4.2lf ", m_VETypeClass[i].m_pdVariance[k]);
		}	
		fprintf(fp, "\n");

		fprintf(fp, "样本数：%d\n", m_VETypeClass[i].m_VESwatch.size());
		fprintf(fp, "样本号:样本特征向量  ");
		for(int j=0; j<m_VETypeClass[i].m_VESwatch.size(); j++)
		{
			fprintf(fp, " %d: ", m_VETypeClass[i].m_VESwatch[j]+1);
			for(int k=0; k<m_iPattDimen; k++)
			{
				fprintf(fp, "%4.2lf ", m_pdSwatchData[m_VETypeClass[i].m_VESwatch[j]*m_iPattDimen+k]);
			}
			fprintf(fp, "\n                     ");
		}
		fprintf(fp, "\n");
		fprintf(fp, "///////////////////////////////////////////////////////\n");
	}
	fclose(fp);
	return 1;
}