isodata.cpp

用C语言实现了ISODATA算法,包括ISODATA.vcproj

// ISODATA.cpp

#include "stdafx.h"
#include "ISODATA.h"
#include "Sort.h"

using namespace std;

// class ISODATA

extern int N;
extern int dim;

ISODATA::ISODATA()
{
	c = 2;	// 预期的类数；
	Nc = 1;	// 初始聚类中心个数（可以不等于c）；
	theta_n = 2;	// 每一类中允许的最少模式数目（若少于此数，就不能单独成为一类）；
	theta_s = 1;	// 类内各分量分布的标准差上限（大于此数就分裂）；
	theta_D = 4;	// 两类中心间的最小距离下限（若小于此数，这两类应合并）；
	L = 1;	// 在每次迭代中可以合并的类的最多对数；
	I = 4;	// 允许的最多迭代次数；

	_d = 100;			// 总体平均距离
	Ip = 0;				// 迭代次数
	
	double D[MAXNUM][MAXNUM];	// 各类对中心间的距离

	for(int i=0;i<MAXNUM;i++)
		for(int j=0;j<MAXNUM;j++)
			D[i][j] = MAXDOUBLE;
}

ISODATA::~ISODATA()
{
}

// 设置参数
int ISODATA::SetupPattern(Pattern * pattern)
{
	for(int i=0;i<N;i++)
		x[i] = pattern[i];

	return N;
}

// 算法实现步骤
int ISODATA::Process()
{
	bool changed = true;
	// 1.预置

	// 2)将待分类的模式特征矢量x1,x2,...,xn读入；
	// SetupPattern();

	
	// 3)选定初始聚类中心，可从待分类的模式特征矢量集{xi}中任选Nc个模式特征矢量作为初始聚类中心zj（j=1,2,...,Nc).
	InitCenter();

step1:
	// 1)设定聚类分析控制参数：
	SetupParameter();

step2:
	// 2.按最小距离原则将模式集（xi）中每个模式分到某一类中
	changed = false;
	Clustering();
	if(Ip == 0)
		cout << endl << "------------- 选取初始聚类中心 ---------------" << endl;
	else
		cout << endl << "-------------第" << Ip << "次迭代-------------" << endl;
	PrintSort();

step3:
	// 3.依据theta_n判断合并。如果类wj中样本数nj<theta_n，则取消该类的中心zj，Nc = Nc - 1；转至 2.
	if(CombinBytheta_n())
		goto step2;
step4:
	// 计算分类后的参数：各类中心、类内平均距离及总体平均距离。
	CalParameter();

step5:
	// 依据Ip, Nc 判断停止\分裂或合并。
	if(Ip == I)
	{
		theta_D = 0;
		goto step9;
	}
	else if(Nc <= c/2)
		goto step6;
	else if(Nc >= 2*c)
		goto step9;
	else if(Ip%2 == 1)	//Nc> c/2 && Nc < 2*c
		goto step6;
	else
		goto step9;

step6:
	// 计算各类类内距离的标准差矢量
	CalSigma();

step7:
	// 求出每一聚类类内距离标准差矢量sigma中的最大分量sigma_max
	// CalMaxSigma();	// 由于在Sort类中已进行计算，这里不做任何处理

step8:
	// 判断分裂
	if(Split())
	{
		changed = true;
		Ip++;
		goto step2;
	}

step9:
	// 计算各类对中心间的距离
	CalCenterDis();

step10:
	// 依据theta_D判断合并
	if(CombinBytheta_D())
		changed = true;

step11:
	// 判断循环还是退出
	if(Ip >= I)
	{
		cout << endl << endl << "==========经过" << Ip << "次迭代，达到迭代次数===========" << endl;
		goto over;
	}
	else if(changed == false)
	{
		Ip++;
		cout << endl << endl << "==========经过" << Ip << "次迭代，算法收敛===========" << endl;
		goto over;
	}
	else
	{
		Ip++;

		char ch;
		cout << "本次迭代完成，是否需要改变参数（Y/N）?: ";
		cin >> ch;
		if(ch == 'y' || ch == 'Y')
			goto step1;
		else 
			goto step2;
	}

over:	
	PrintSort();
	return Ip;
}

// 3)选定初始聚类中心
void ISODATA::InitCenter()
{
	srand( (unsigned)time( NULL ) );
	int num = N;
	int no[MAXNUM];

	for(int i=0;i<N;i++)
		no[i] = i;

	for(int j=0;j<Nc;j++)
	{		
		int k = rand()%num;
		w[j].z = x[no[k]];
		w[j].z.n = 0;
		for(int l=k;l<num;l++)
			no[l] = no[l+1];
		num--;
	}
}

// step1:	设定聚类分析控制参数
void ISODATA::SetupParameter()
{
	cout << "设定聚类分析控制参数:" << endl;
	cout << "预期的类数 c:";
	cin >> c;
	cout << "初始聚类中心个数（可以不等于c）Nc:";
	cin >> Nc;
	cout << "每一类中允许的最少模式数目（若少于此数，就不能单独成为一类）theta_n:";
	cin >> theta_n;
	cout << "类内各分量分布的标准差上限（大于此数就分裂）theta_s:";
	cin >> theta_s;
	cout << "两类中心间的最小距离下限（若小于此数，这两类应合并）theta_D:";
	cin >> theta_D;
	cout << "在每次迭代中可以合并的类的最多对数 L：";
	cin >> L;
	cout << "允许的最多迭代次数 I：";
	cin >> I;
}

// step2:	按最小距离原则将模式集（xi）中每个模式分到某一类中
void ISODATA::Clustering()
{
	double temp = 0.0, min = MAXDOUBLE;
	int l = 0;

	for(int j=0;j<Nc;j++)
	{
		w[j].n = 0;
		w[j].z.n = 0;
	}

	for(int i=0;i<N;i++)
	{
		min = MAXDOUBLE;
		l=0;

		for(int j=0;j<Nc;j++)
		{
			temp = Pattern::Distance(x[i], w[j].z);
			if(min > temp)
			{
				min = temp;
				l = j;
			}
		}
		w[l].Insert(x[i]);
	}
}

// step3:	依据theta_n判断合并。如果类wj中样本数nj<theta_n，则取消该类的中心zj，Nc = Nc - 1；转至 step2.
bool ISODATA::CombinBytheta_n()
{
	int j=0;
	do
	{
		if(w[j].n < theta_n)
		{
			for(int k=j;k<Nc;k++)
				//w[k] = w[k+1];
				w[k].z = w[k+1].z;
			Nc--;
			// 循环中跳出！？
			return true;
		}
		j++;
	}while(j < Nc);

	return false;
}

// step4:	计算分类后的参数：各类中心、类内平均距离及总体平均距离。
void ISODATA::CalParameter()
{
	_d = 0.0;

	for(int j=0;j<Nc;j++)
	{
		w[j].CalCenter();		
		_d += w[j].n * w[j].Cal_D();
	}
	_d /= N;
}

// step6:	计算各类类内距离的标准差矢量
void ISODATA::CalSigma()
{
	for(int j=0;j<Nc;j++)
		w[j].CalSigma();
}

// step7:	求出每一聚类类内距离标准差矢量sigma中的最大分量sigma_max
void ISODATA::CalMaxSigma()
{
	// 由于在Sort类中已进行计算，这里不做任何处理
}

// step8:	判断分裂
bool ISODATA::Split()
{
	for(int j=0;j<Nc;j++)
	{
		if( ( (w[j]._d>_d) && (w[j].n>2*(theta_n+1)) ) || (Nc <= c/2) )
		{
			int i = w[j].max;
			double sigma = w[j].sigma_max;

			for(int l=Nc;l>j;l--)
				w[l].z = w[l-1].z;

			w[j].z.x[i] -= K*sigma;
			w[j+1].z.x[i] += K*sigma;

			Nc++;
			return true;
		}
	}

	return false;
}

// step9:	计算各类对中心间的距离
void ISODATA::CalCenterDis()
{
	for(int i=0;i<Nc-1;i++)
		for(int j=i+1;j<Nc;j++)
			D[i][j] = Pattern::Distance(w[i].z, w[j].z);
}

// step10:	依据theta_D判断合并
bool ISODATA::CombinBytheta_D()
{
	int i, j, k, l;	// 循环变量
	int num = 0;
	bool b = false;

	// 较小的类对中心
	struct
	{
		double d;
		int i;
		int j;
	}Dmin[MAXNUM];

	for(i=0;i<L;i++)
	{
		Dmin[i].d = 0.0;
		Dmin[i].i = -1;
		Dmin[i].j = -1;
	}

	// 将D[i][j]与theta_D比较，并将小于theta_D的那些D[i][j]按递增次序排列，取前L个，Dmin[0]<Dmin[1]<Dmin[2]...
	for(i=0;i<Nc-1;i++)
		for(j=i+1;j<Nc;j++)
			if(D[i][j] < theta_D)
				for(k=0;k<L;k++)
					if(D[i][j] < Dmin[k].d)
					{
						for(l=L-1;l>k;l--)
							Dmin[l] = Dmin[l-1];
						Dmin[k].d = D[i][j];
						Dmin[k].i = i;
						Dmin[k].j = j;
						break;
					}

	// 从最小的Dmin开始，将相应的两类合并。
	for(i=0;i<L;i++)
		if(Dmin[i].i > -1 && Dmin[i].j > -1)
			// 合并两类
			if( w[Dmin[i].i].Combin(w[Dmin[i].j]) )
				b = true;

	// 去掉已取消的类心
	for(j=0;j<Nc;j++)
		if(w[j].z.n == -2)
		{
			for(k=j;k<Nc;k++)
				w[k].z = w[k+1].z;
			Nc--;
		}

	return b;
}

// step11:	判断循环还是退出

// 输出当前模式分类情况
void ISODATA::PrintSort()
{	
	cout << "共分为分为" << Nc << "类。" << endl;
	for(int i=0;i<Nc;i++)
	{
		cout << endl << "第" << i+1 << "类类心为:\t(";
		for(int j=0;j<dim-1;j++)
			cout << setw(5) << w[i].z.x[j] << ", ";
		cout << setw(5) << w[i].z.x[dim -1] << ")" << endl;

		cout << "包含的模式为:" << endl;
		for(int k=0;k<w[i].n;k++)
		{
			// Xi(12, 12, 12); 
			cout << "\tX" << w[i].x[k].n << "(";
			for(int j=0;j<dim-1;j++)
				cout << setw(5) << w[i].x[k].x[j] << ", ";
			cout << setw(5) << w[i].x[k].x[dim-1] << ");" << endl;
		}
		cout << endl;
	}
}