cluster.cpp

该程序通过遗传算法对图像进行聚类分析

		{
			a+=t1[i]*t2[i];
			b1+=t1[i]*t1[i];
			b2+=t2[i]*t2[i];
		}
		delete	[]t1;
		delete	[]t2;
		if ((b2+b1-a)!=0)
		result=(double)(a/(b1+b2-a));
		else 
		{
			return -1;
		}
		return (1-result);
	}
	else
		return -1;
}
/*************************************************************
*函数名称		GetDistance(CCluster::Center mCenter1,CCluster::Center mCenter2,int distype)
*参数			CCluster::Center mCenter1   中心1
*				CCluster::Center mCenter2   中心2
*				const int distype	         距离模式   1：欧氏距离；2：夹角余弦距离；		
*											            3：特征是二值时的夹角余弦距离；
*														4：具有二值特征的Tanimoto测度
*返回值			double
*函数功能		计算两个聚类中心间的距离，距离模式由distype给定
*************************************************************/
double CCluster::GetDistance(CCluster::Center mCenter1,CCluster::Center mCenter2,int distype)
{
	double result;
	result=0;

	if (distype==1)//欧氏距离
	{
		for ( int m=0;m<N*N;m++)
		result+=(mCenter1.feature[m]-mCenter2.feature[m])*(mCenter1.feature[m]-mCenter2.feature[m]);
		return (double)sqrt(result);
	}
	else if (distype==2)//夹角余弦
	{
		double	a,b1,b2;
		a=0;
		b1=0;
		b2=0;
		for (int i=0; i<N*N; i++)
		{
			a+=mCenter1.feature[i]*mCenter2.feature[i];
			b1+=mCenter1.feature[i]*mCenter1.feature[i];
			b2+=mCenter2.feature[i]*mCenter2.feature[i];
		}
		if (b2*b1!=0)
		result=a/sqrt(b1*b2);
		else 
		{
			return -1;
		}

		return (1-result);
	}
	else if (distype==3)//二值夹角余弦
	{
		int *t1,*t2;
		int a,b1,b2;
		 
		a=0;b1=0;b2=0;
		t1=new int [N*N];
		t2=new int [N*N];
		for(int i=0; i<N*N; i++)
		{
			t1[i]=mCenter1.feature[i]>0.2? 1:0;
			t2[i]=mCenter2.feature[i]>0.2? 1:0;
		}
		
		for ( i=0; i<N*N; i++)
		{
			a+=t1[i]*t2[i];
			b1+=t1[i]*t1[i];
			b2+=t2[i]*t2[i];
		}
		delete	[]t1;
		delete	[]t2;
		if (b2*b1!=0)
		result=(double)(a/sqrt(b1*b2));
		else 
		{
			return -1;
		}

		return (1-result);
	}
	else if (distype==4)//Tanimoto
	{
		int *t1,*t2;
		int a,b1,b2;
		 
		a=0;b1=0;b2=0;
		t1=new int [N*N];
		t2=new int [N*N];
		for(int i=0; i<N*N; i++)
		{
			t1[i]=mCenter1.feature[i]>0.2? 1:0;
			t2[i]=mCenter2.feature[i]>0.2? 1:0;
		}
		
		for ( i=0; i<N*N; i++)
		{
			a+=t1[i]*t2[i];
			b1+=t1[i]*t1[i];
			b2+=t2[i]*t2[i];
		}
		delete	[]t1;
		delete	[]t2;
		if ((b2*b1-a)!=0)
		result=(double)(a/(b1+b2-a));
		else 
		{
			return -1;
		}
		return (1-result);
	}
	else
		return -1;	
}

/************************************************************
*函数名称		GetDistance(double *dis, int i, int j)
*参数			double *dis    指向模糊距离数组的指针大小为[patternnum*patternnum]
*				int i			数组的第i行
*				int j			数组的第j列
*返回值		    double
*函数功能		返回模糊距离数组中第i行第j列的模糊积
*************************************************************/
double CCluster::GetDistance(double *dis, int i, int j)
{
	double result;
	result=0;

	for(int t=0;t<patternnum;t++)
	{
		double td=dis[i*patternnum+t]<dis[t*patternnum+j]? dis[i*patternnum+t]:dis[t*patternnum+j];
		if (result<td)
			result=td;
	}
	return result;
}
/************************************************************
*函数名称		GetFuzzyDistance(CCluster::Pattern pattern1,CCluster::Pattern pattenr2, int distype)
*参数			CCluster::Pattern pattern1    样品1
*				CCluster::Pattern pattern2    样品2
*				int distype					  距离计算形式：1，欧氏距离；2，数量积；
*															3，相关系数
*															4，最大最小法；5，算数平均法; 
*															6,几何平均最小法
*返回值         double
*功能			计算两个样品的模糊距离，距离模式由参数distype给出
************************************************************/
double CCluster::GetFuzzyDistance(CCluster::Pattern pattern1,CCluster::Pattern pattern2, int distype)
{
	switch	(distype)
	{
	case 1://欧氏距离
		{
				double max=0;
				for (int i=0; i<patternnum-1; i++)
					for (int j=i+1; j<patternnum; j++)
						if (max<GetDistance(m_pattern[i],m_pattern[j],1))
							max=GetDistance(m_pattern[i],m_pattern[j],1);
				return (max-GetDistance(pattern1,pattern2,1))/max;
		}
	case 2://数量积
		{
			double temp,max;
			max=0;
			for (int i=0; i<patternnum-1; i++)
				for (int j=0; j<patternnum; j++)
				{
					temp=0;
					for (int k=0; k<N*N; k++)
					{
						temp+=m_pattern[i].feature[k]*m_pattern[j].feature[k];
					}
					if (max<temp)
						max=temp;
				}
			temp=0;
			for ( i=0; i<N*N; i++)
			{
				temp+=pattern1.feature[i]*pattern2.feature[i];
			}
			return (temp/max);
		}
	case 3://相关系数
		{
			double ap1,ap2;
			ap1=0;ap2=0;
			for (int i=0; i<N*N; i++)
			{
				ap1+=pattern1.feature[i];
				ap2+=pattern2.feature[i];
			}
			ap1/=N; 
			ap2/=N;

			double a,b1,b2;
			a=0;b1=0;b2=0;

			for (i=0; i<N*N; i++)
			{
				a+=(pattern1.feature[i]-ap1)*(pattern2.feature[i]-ap2);
				b1+=(pattern1.feature[i]-ap1)*(pattern1.feature[i]-ap1);
				b2+=(pattern2.feature[i]-ap2)*(pattern2.feature[i]-ap2);
			}
			if (b2*b1!=0)
				return (a/sqrt(b1*b2));
		}
	case 4://最大最小法
		{
			double min ,max;
			min=0; max=0;
			for (int i=0; i<N*N; i++)
			{
				min+=pattern1.feature[i]<pattern2.feature[i]? pattern1.feature[i]:pattern2.feature[i];
				max+=pattern1.feature[i]<pattern2.feature[i]? pattern2.feature[i]:pattern1.feature[i];
			}
			if (max!=0)
			return (min/max);
		}
	case 5://算数平均法
		{
			double min ,max;
			min=0; max=0;
			for (int i=0; i<N*N; i++)
			{
				min+=pattern1.feature[i]<pattern2.feature[i]? pattern1.feature[i]:pattern2.feature[i];
				max+=pattern1.feature[i]+pattern2.feature[i];
			}
			if (max!=0)
			return (2*min/max);	
		}
	case 6://几何平均最小法
		{
			double min ,max;
			min=0; max=0;
			for (int i=0; i<N*N; i++)
			{
				min+=pattern1.feature[i]<pattern2.feature[i]? pattern1.feature[i]:pattern2.feature[i];
				max+=sqrt(pattern1.feature[i]*pattern2.feature[i]);
			}
			if (max!=0)
			return (min/max);			
		}
	default:
		return -1;
	}
}


void CCluster::FuzzyCluster()
{
	GetFeature();//获得所有样品特征

	double	*dis, *tempdis;
	double	dismax;
	int		i,j;
	int		distype;
	dismax=0;

	dis=new double [patternnum*patternnum];//模糊系数矩阵
	tempdis=new double [patternnum*patternnum];
	DlgFuzzyDistance  dfd;
	if (dfd.DoModal()==IDOK)
	{
		distype=dfd.mfuzzydistance;
	}
	else 
		return;
    //得到初始模糊距离
	for ( i=0; i<patternnum; i++)
		for (j=0; j<patternnum; j++)
	{
		dis[i*patternnum+j]=GetFuzzyDistance(m_pattern[i],m_pattern[j],distype);//distype
	}

   //构造等价类
	bool flag;
	flag=true;
	while (flag)
	{
		flag=false;
		for ( i=0; i<patternnum; i++)
			for ( j=0; j<patternnum; j++)
				if (i==j)//对角线为1
					tempdis[i*patternnum+j]=1;
				else
					tempdis[i*patternnum+j]=GetDistance(dis,i,j);
		for ( i=0; i<patternnum; i++)
		{
			for ( j=0; j<patternnum; j++)
				if ((tempdis[i*patternnum+j]-dis[i*patternnum+j])*(tempdis[i*patternnum+j]-dis[i*patternnum+j])>0.000001)//(tdis[i][j]!=dis[i][j])
				{
					flag=true;
					break;
				}
			if (flag)
				break;
		}
	
		for ( i=0; i<patternnum*patternnum; i++)
				dis[i]=tempdis[i];
	}
	//输出模糊矩阵
	CString ts,sout;
	//第一行标号1，2，。。。，patternnum。
	sout.Insert(sout.GetLength(),"   \t");
	for (i=1; i<=patternnum; i++)
	{
		ts.Format("%d   \t",i);
		sout.Insert(sout.GetLength(),ts);
	}
	sout.Insert(sout.GetLength(),"\n\n");

	for (i=0; i<patternnum; i++)
	{
		//每列开头标号
		ts.Format("%d   \t",i+1);
		sout.Insert(sout.GetLength(),ts);
		//输出模糊系数矩阵
		for(j=0; j<patternnum; j++)
		{
			ts.Format("%0.3f\t",dis[i*patternnum+j]);
			sout.Insert(sout.GetLength(),ts);
			if ((j+1)%patternnum==0)
				sout.Insert(sout.GetLength(),"\n\n\n");

		}
	}
	MessageBox(NULL,sout,"模糊矩阵",MB_OK);

	delete []tempdis;

	double *xishu=new double [patternnum*patternnum];
	for(i=0; i<patternnum*patternnum; i++)
		xishu[i]=-1;

	int pointer=0;
	//记录模糊系数矩阵中不同的系数
	for ( i=0; i<patternnum; i++)
		for ( j=i; j<patternnum; j++)
		{
			bool done=false;
			for (int k=0; k<pointer; k++)
			{
				if ((xishu[k]-dis[i*patternnum+j])*(xishu[k]-dis[i*patternnum+j])<0.000001)
				{
					done=true;
					break;
				}
			}
			if (!done)
			{
				xishu[pointer]=dis[i*patternnum+j];
				pointer++;
			}
		}
	
	for(i=0; i<pointer-1; i++)//对阈值由小到大排序
		for (j=0; j<pointer-i-1; j++)
		{
			if (xishu[j]>xishu[j+1])
			{
				double	temp=xishu[j];
				xishu[j]=xishu[j+1];
				xishu[j+1]=temp;
			}
		}
	CString s,str;
	for (i=0; i<pointer; i++)
	{
		s.Format("%0.4f  ",xishu[i]);
		str.Insert(str.GetLength(),s);
	}
	delete   []xishu;
	//用户输入聚类阈值
	double    yz;


	
	int *result;
	//根据阈值输出聚类结果
	result=new int [patternnum*patternnum];
	for (i=0; i<patternnum*patternnum; i++)
		if (dis[i]>=yz)
			result[i]=1;
		else 
			result[i]=0;

		//分类后输出

	//第一行标号
	sout="";
	sout.Insert(-1,"每行中\"1\"对应的列为同一类");
	sout.Insert(sout.GetLength(),"\n\n");
	sout.Insert(sout.GetLength(),"   \t");
	for (i=1; i<=patternnum; i++)
	{
		ts.Format("%d   \t",i);
		sout.Insert(sout.GetLength(),ts);
	}
	sout.Insert(sout.GetLength(),"\n\n");
		

	for (i=0; i<patternnum; i++)
	{
		//每列开头标号
		ts.Format("%d   \t",i+1);
		sout.Insert(sout.GetLength(),ts);
		for(j=0; j<patternnum; j++)//首行为原模糊系数
		{
			ts.Format("%0.3f\t",dis[i*patternnum+j]);
			sout.Insert(sout.GetLength(),ts);
			if ((j+1)%patternnum==0)
				sout.Insert(sout.GetLength(),"\n");
		}
		sout.Insert(sout.GetLength(),"   \t");
		for(j=0; j<patternnum; j++)//次行为根据阈值修改后的系数“1”或“0”
		{
			ts.Format("%d\t",result[i*patternnum+j]);
			sout.Insert(sout.GetLength(),ts);
			if ((j+1)%patternnum==0)
				sout.Insert(sout.GetLength(),"\n\n");
		}
	}
	MessageBox(NULL,sout,"分类前后的矩阵对照",MB_OK);

	centernum=0;

	for (i=0; i<patternnum; i++)//按照阈值分类
		for (j=i; j<patternnum; j++)
		{
			if (result[i*patternnum+j]==1)
			{
				if (m_pattern[i].category!=0)
					m_pattern[j].category=m_pattern[i].category;
				else if (m_pattern[j].category!=0)
					m_pattern[i].category=m_pattern[j].category;
				else
				{
					centernum++;
					m_pattern[j].category=centernum;
					m_pattern[i].category=centernum;
				}
			}
		}
	delete []dis;
	delete []result;
}