isodata.cpp

这是ISODATA算法,这个算法具有自组织性,跟大家分享.

// ISODATA.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "math.h"


////模式识别与智能技术 李刚 0720080195

//规定控制参数
#define K 4
#define setaN 10
#define setaS 3
int setaC = 2;
#define L 3
#define I 1000
#define C 5

#define MAX_ROW 150
#define MAX_COL 4

typedef struct Data
{
	float data[MAX_ROW][MAX_COL];
}Data;

typedef struct Cluster
{
	bool data[MAX_ROW];
	float center[MAX_COL];
	float dif[MAX_COL];
	float maxDif;
	int maxIndex;
	int count;
	bool used;
	float dis;
	bool isCom;
}Cluster;

typedef struct Distance
{
	float dis;
	int i;
	int j;
}Distance;

// 从文件中读取原始数据
Data* ReadFile(char* path)
{
	Data* data = new Data();
	FILE* file;
	if((file = fopen(path,"rt")) == NULL)
	{
		printf("Can not open file.\n");
		return NULL;
	}
	int i = 0;
	char *name = "";
	for(; i < MAX_ROW; i++)
	{
		char buffer[150];
		fgets(buffer,150,file);
		sscanf(buffer,"%f,%f,%f,%f",&data->data[i][0],&data->data[i][1],&data->data[i][2],&data->data[i][3]);
	}
	fclose(file);

	return data;
}

// 初始化
void Init(Data* data,Cluster* clusters[], int c)
{
	int f = 0;
	int j = 0;
	for(f = 0; f < C ; f++)
	{
		clusters[f] = new Cluster();
		clusters[f]->data[f] = true;
		clusters[f]->used = true;
		clusters[f]->count = 1;
		clusters[f]->dis = 0;
		clusters[f]->isCom = false;
		clusters[f]->maxDif = 0;
		for(j = 0; j< MAX_COL; j++)
		{
			clusters[f]->center[j] = data->data[f][j];
			clusters[f]->dif[j] = 0;
		}
	}
}

void RefreshCenter(Cluster* clusters[], int length, Data* data);

// 把所有样本分到各个聚类中
void Distribute(Cluster* clusters[], int length, Data* data)
{
	int i = 0;
	int j = 0;
	for(i = 0; i < MAX_ROW; i++)
	{
		float Min = 1000;
		int f = 0;
		int index = -1;
		// 查找最小值
		for(f = 0; f <length; f++)
		{
			if(clusters[f]->used)
			{
				float dis = 0;
				for(j = 0 ; j < MAX_COL; j++)
				{
					dis = dis + (data->data[i][j] - clusters[f]->center[j])*(data->data[i][j] - clusters[f]->center[j]);
				}
				if(dis < Min)
				{
					index = f;
					Min = dis;
				}
			}
		}
		// 把向量i分配到聚类中
		for(f = 0;f < length; f++)
		{
			clusters[f]->data[i] = false;
		}
		clusters[index]->data[i] = true;
		RefreshCenter(clusters,length,data);
	}
}

// 更新聚类的均值向量
void RefreshCenter(Cluster* clusters[], int length, Data* data)
{
	int i = 0;
	int j = 0; 
	int f = 0;
	for(f = 0; f < length; f++)
	{
		if(clusters[f]->used)
		{
			// 针对每一个列计算其均值
			for(j = 0 ; j < MAX_COL; j++)
			{
				float result = 0;
				int count = 0;
				for(i = 0 ; i < MAX_ROW; i++)
				{
					if(clusters[f]->data[i])
					{
						count++;
						result = result + data->data[i][j];
					}
				}
				clusters[f]->center[j] = result / count;
				clusters[f]->count = count;
			}
		}
	}
}

// 计算各个样本离开其中心的平均距离
void GetDis(Cluster* clusters[], int length, Data* data)
{
	int i = 0;
	int j = 0;
	int f = 0;

	for(f = 0 ; f < length; f++)
	{
		if(clusters[f]->used)
		{
			int count = 0;
			float result = 0;
			for(i = 0; i < MAX_ROW; i++)
			{
				if(clusters[f]->data[i])
				{
					count++;
					for(j = 0; j < MAX_COL; j++)
					{
						result = result + (clusters[f]->center[j] - data->data[i][j])*(clusters[f]->center[j] - data->data[i][j]);
					}
				}
			}
			result = result / count;
			clusters[f]->dis = result;
		}
	}
}

// 计算所有样本离开聚类中心的平均值
float GetAllDis(Cluster* clusters[], int length)
{
	int i = 0;
	int count =0;
	float result = 0;
	for(i  = 0; i < length; i++)
	{
		if(clusters[i]->used)
		{
			count ++;
			result += clusters[i]->dis;
		}
	}
	return result/count;
}

// 返回当前聚类的个数
int GetCountOfClusters(Cluster* clusters[], int length)
{
	int i = 0;
	int count = 0;
	for(i = 0; i < length; i++)
	{
		if(clusters[i]->used)
		{
			count++;
		}
	}
	return count;
}

//////////////////////////////////////////////分裂//////////////////

// 计算某一个聚类的标准差向量
void GetDif(Cluster* cluster, Data* data)
{
	int i = 0;
	int j = 0;
	int count = 0;
	float result = 0;

	for(j = 0; j < MAX_COL; j++)
	{
		count = 0;
		result = 0;
		for(i = 0 ; i < MAX_ROW; i++)
		{
			if(cluster->data[i])
			{
				count++;
				result = result + (data->data[i][j] - cluster->center[j]) * (data->data[i][j] - cluster->center[j]);
			}
		}
		cluster->dif[j] = sqrt(result/count);
	}
}

// 计算全部聚类的标准差向量
void GetDifOfAllCluster(Cluster* clusters[], int length, Data* data)
{
	int i = 0;
	for(i = 0; i < length; i++)
	{
		if(clusters[i]->used)
		{
			GetDif(clusters[i],data);
		}
	}
}

// 计算某一个聚类的最大标准偏差
void GetMaxDif(Cluster* cluster)
{
	int j = 0;
	int index = 0;
	float max = -100;
	for(j = 0; j < MAX_COL; j++)
	{
		if(cluster->dif[j] > max)
		{
			index = j;
			max = cluster->dif[j];
		}
	}
	cluster->maxDif = max;
	cluster->maxIndex = index;
}

// 计算全部聚类的标准差
void GetAllMaxDif(Cluster* clusters[], int length, Data* data)
{
	int f = 0;

	for(f = 0; f <length; f++)
	{
		GetDif(clusters[f],data);
	}
	for(f = 0; f < length; f++)
	{
		if(clusters[f]->used)
		{
			GetMaxDif(clusters[f]);
		}
	}
}

// 分裂某一聚类
void DivideCluster(Cluster* clusters[], int length, Cluster* cluster, float allDis, int k)
{
	if(cluster->maxDif > setaS &&
		(cluster->dis >= allDis && cluster->count > 2*(setaN + 1)
		|| GetCountOfClusters(clusters,C) <= K/2))
	{
		int i = 0; 
		for(i = 0; i < length; i++)
		{
			if(clusters[i] != cluster && !clusters[i]->used)
			{
				clusters[i]->used = true;
				clusters[i]->center[cluster->maxIndex] =
					cluster->center[cluster->maxIndex] - k* cluster->maxDif;
				cluster->center[cluster->maxIndex] = cluster->center[cluster->maxIndex] + k* cluster->maxDif;
				break;
			}
		}
	}
}

// 分裂阶段
void DivideAllClusters(Cluster* clusters[], int length, Data* data)
{
	int f = 0;
	float allDis = GetAllDis(clusters,length);
	GetAllMaxDif(clusters,length,data);
	for(f = 0; f < length; f++)
	{
		if(clusters[f]->used)
		{
			DivideCluster(clusters,length,clusters[f],allDis,0.3);
		}
	}
}

/////////////////////////////////合并处理///////////////////////////////////////

// 合并阶段
void Combine(Cluster* clusters[], int length)
{
	Distance* distances[500];
	
	// 计算两两聚类之间的距离
	int i = 0;
	int j = 0;
	int f = 0;
	for(i = 0; i < length; i++)
	{
		if(clusters[i]->used)
			for(j = i+1; j < length; j++)
			{
				if(clusters[j]->used)
				{
					int col =0;
					float result = 0;
					for(col = 0; col < MAX_COL; col++)
					{
						result = result + abs((clusters[i]->center[col] - clusters[j]->center[col]));
					}
					Distance* dis = new Distance();
					dis->dis = result;
					dis->i = i;
					dis->j = j;
					distances[f] = dis;
					f++;
				}
			}
	}

	//对小于阈值的距离进行排序.
	int count = f;
	for(f = 0; f < count; f++)
	{
		int a = 0;
		float maxValue = -100;
		int c = -1;
		for(a = f; a < count; a++)
		{
			if(distances[a]->dis > maxValue)
			{
				c = a;
				maxValue = distances[a]->dis;
			}
		}
		Distance* temp = distances[f];
		distances[f] = distances[c];
		distances[c] = temp;
	}

	//将两个中心合并
	for(f = 0; f < count; f++)
	{
		if(distances[f]->dis < setaC &&
			clusters[distances[f]->i]->used && clusters[distances[f]->j]->used &&
			!clusters[distances[f]->i]->isCom && !clusters[distances[f]->j]->isCom)
		{
			clusters[distances[f]->i]->isCom = true;
			clusters[distances[f]->j]->isCom = true;
			clusters[distances[f]->j]->used = false;
			int col = 0;
			for(col = 0; col < MAX_COL; col++)
			{
				clusters[distances[f]->i]->center[col] = clusters[distances[f]->i]->center[col] + clusters[distances[f]->j]->center[col];
			}
			clusters[distances[f]->j]->used = false;
		}
	}
}

int _tmain(int argc, _TCHAR* argv[])
{
	// 初始化数据结构，从文件中读取原始数据。
	Data* data = ReadFile("iris.txt");
	if(data == NULL)
	{
		return 0;
	}
	Cluster* clusters[C];
	Init(data,clusters,C);

	int num = I;
	while(num--)
	{
		// 将所有的样本分到各个聚类中
		Distribute(clusters,C,data);

		//如有任何一个聚类的样本数少于setaN，则舍去该聚类
		int i = 0; 
		for(i = 0; i < C; i++)
		{
			if(clusters[i]->used && clusters[i]->count < setaN)
			{
				clusters[i]->used = false;
			}
		}

		// 更新聚类的均值
		RefreshCenter(clusters,C,data);

		//计算各个样本离开聚类中心的平均值
		GetDis(clusters,C,data);

		// 若迭代运算次数已达到I次，即最后一次迭代，则置θc =0
		if(num !=0)
		{
			// 若聚类中心的数目小于或等于规定值的一半则进行分裂处理
			if(GetCountOfClusters(clusters,C) <= K/2)
			{
				// 分裂处理
				DivideAllClusters(clusters,C,data);
				continue;
			}
			// 若迭代预算次数是偶数次,或者现有聚类数目大于期望的两倍,则不进行分裂,否则进行分裂处理.
			else if(num%2 == 0 && GetCountOfClusters(clusters,C) > 2*K)
			{
				// 分裂处理
				DivideAllClusters(clusters,C,data);
				continue;
			}
			// 上述条件不满足其一,都进行合并处理
		}
		if(num ==0)
		{
			// 若为最后一次则进行合并处理
			setaC = 0;
		}
		//合并处理阶段
		Combine(clusters,C);
		if(num == 0)
		{
			break;
		}
		
	}
	return 0;
}