图像模式识别之几何分类器(二)

#include "cv.h"
#include "highgui.h"
#include "utils.h"

#include<stack>

using namespace cv;
using namespace std;

#define PI 3.1415926
#define RADIAN(angle) ((angle)*PI/180.0)

struct pattern
{
	int number; //该类型数字样品个数
	double feature[200][25];//该类型数字样品特征，每个特征为25维
	pattern()
	{
		number=0;
	}
};

struct Position
{
	int bottom;
	int top;
	int left;
	int right;
};

pattern ptns[10]; //一共10个数字
double testsample[25]; //手写数字的特征 25维

//判断手写数字可否保存为cls类别的一个样品 同一类别不能有重复
bool savable(int cls)
{
	for(int i=0;i<ptns[cls].number;i++)
	{
		bool flag=true;
		for(int j=0;j<25;j++)
			flag&=(fabs(ptns[cls].feature[i][j]-testsample[j])<0.0000001 ? true:false);
		if(flag) return false;
	}

	return true;
}

//将手写数字保存为cls类别中的一个样品
void Save(int cls)
{
	for(int i=ptns[cls].number;i>0;i--)
		for(int j=0;j<25;j++)
			ptns[cls].feature[i][j]=ptns[cls].feature[i-1][j];

	for(int j=0;j<25;j++)
		ptns[cls].feature[0][j]=testsample[j];

	ptns[cls].number++;
}

//计算手写数字的在图片中出现的位置
Position GetPosition(Mat& mat)
{
	int M=mat.rows;
	int N=mat.cols;

	bool flag=false;
	int i,j;
	Position ret;

	//top
	for(i=0;i<M-1;i++)
	{
		flag=false;
		for(j=0;j<N;j++)
		{
			if(mat.at<uchar>(i,j)<128)
			{
				flag=true;
				break;
			}
		}
		if(flag) break;
	}
	ret.top=i;

	//bottom
	for(i=M-1;i>0;i--)
	{
		flag=false;
		for(j=0;j<N;j++)
		{
			if(mat.at<uchar>(i,j)<128)
			{
				flag=true;
				break;
			}
		}
		if(flag) break;
	}
	ret.bottom=i;

	//left
	for(j=0;j<N-1;j++)
	{
		flag=false;
		for(i=0;i<M;i++)
		{
			if(mat.at<uchar>(i,j)<128)
			{
				flag=true;
				break;
			}
		}
		if(flag) break;
	}
	ret.left=j;

	for(j=N-1;j>0;j--)
	{
		flag=false;
		for(i=0;i<M;i++)
		{
			if(mat.at<uchar>(i,j)<128)
			{
				flag=true;
				break;
			}
		}
		if(flag) break;
	}
	ret.right=j;

	return ret;
}

double Cal(int row,int col,Position pos,Mat& mat)
{
	int w,h;
	w=(pos.right-pos.left)/5;
	h=(pos.bottom-pos.top)/5;

	int count=0;

	for(int i=pos.top+row*h;i<pos.top+(row+1)*h;i++)
	{
		for(int j=pos.left+col*w;j<pos.left+(col+1)*w;j++)
		{
			if(mat.at<uchar>(i,j)<128)
				count++;
		}
	}
	return ((double)count)/(w*h);
}
//将手写数字特征保存在变量testsample中
void SetFeature(Mat& mat,Position pos)
{
	for(int i=0;i<5;i++)
	{
		for(int j=0;j<5;j++)
		{
			testsample[5*i+j]=Cal(i,j,pos,mat);
		}
	}
}

//利用fisher算法进行2类分类
int Fisher2Classes(int class0,int class1)
{
	double xMeans[2][25]; //均值
	double S[2][25][25]; //类内离散度矩阵
	double SW[25][25];//类间离散度矩阵
	double SW_[25][25]; //SW的逆矩阵
	double w[25]; //解向量
	double diffXMeans[25]; //均值差
	double X[25]; //测试样品
	
	//计算两类样本的均值
	for(int i=0;i<2;i++)
		for(int j=0;j<25;j++)
			xMeans[i][j]=0;
	for(int i=0;i<ptns[class0].number;i++)
		for(int j=0;j<25;j++)
			xMeans[0][j]+=ptns[class0 ].feature[i][j];
	for(int i=0;i<ptns[class1].number;i++)
		for(int j=0;j<25;j++)
			xMeans[1][j]+=ptns[class1].feature[i][j];
	for(int i=0;i<25;i++)
	{
		xMeans[0][i]/=(double)ptns[class0].number;
		xMeans[1][i]/=(double)ptns[class1].number;
	}

	//计算类内离散度矩阵
	for(int i=0;i<25;i++)
		for(int j=0;j<25;j++)
		{
			double s0=0,s1=0;
			for(int k=0;k<ptns[class0].number;k++)
				s0+=(ptns[class0].feature[k][i]-xMeans[0][i])*(ptns[class0].feature[k][j]-xMeans[0][j]);
			s0/=(double)(ptns[class0].number-1);
			S[0][i][j]=s0;

			for(int k=0;k<ptns[class1].number;k++)
				s1+=(ptns[class1].feature[k][i]-xMeans[1][i])*(ptns[class1].feature[k][j]-xMeans[1][j]);
			s1/=(double)(ptns[class1].number-1);
			S[1][i][j]=s1;
		}

	//类间离散度矩阵
	for(int i=0;i<25;i++)
		for(int j=0;j<25;j++)
		{
			SW[i][j]=S[0][i][j]+S[1][i][j];
			SW_[i][j]=SW[i][j];
		}

	//SW的逆矩阵
	double (*p)[25]=SW_;
	brinv(*p,25);

	//计算W
	for(int i=0;i<25;i++)
		diffXMeans[i]=xMeans[0][i]-xMeans[1][i];
	for(int i=0;i<25;i++)
		w[i]=0.0;
	brmul(SW_,diffXMeans,25,w);

	//映射后的样品均值
	double m0=0,m1=0;
	for(int i=0;i<ptns[class0].number;i++)
		m0+=brmul(w,ptns[class0].feature[i],25);
	for(int i=0;i<ptns[class1].number;i++)
		m1+=brmul(w,ptns[class1].feature[i],25);

	m0/=(double)ptns[class0].number;
	m1/=(double)ptns[class1].number;
	double y0=(ptns[class0].number*m0+ptns[class1].number*m1)/(ptns[class0].number+ptns[class1].number);

	for(int i=0;i<25;i++)
		X[i]=testsample[i];
	double y=brmul(w,X,25);
	if(y>=y0) return class0;
	else return class1;
}

int Fisher()
{
	int num[10];
	for(int i=0;i<10;i++)
		num[i]=0;
	//每两类进行比较
	for(int i=0;i<10;i++)
		for(int j=0;j<10;j++)
			num[Fisher2Classes(i,j)]++;
	//求最大值
	int maxVal=num[0];
	int ret=0;
	for(int i=1;i<10;i++)
		if(maxVal<num[i])
		{
			maxVal=num[i];
			ret=i;
		}

	return ret;
}

//势函数法
double K(double X[],double Xk[])
{
	double d=0;
	for(int i=0;i<25;i++)
		d+=(X[i]-Xk[i])*(X[i]-Xk[i]);
	return exp(-d);
}

int main()
{
	//读取图像
	Mat image;

	//样本库特征提取
	int cls;
	char filename[100]="a.txt";
	char filePath[100];
	FILE* fp=fopen(filename,"r");
	while(fscanf(fp,"%d%s",&cls,filePath))
	{
		if(cls==-1)
			break;
		Mat image=imread(filePath,CV_LOAD_IMAGE_GRAYSCALE);
		Position pos=GetPosition(image);
		SetFeature(image,pos);
		if(savable(cls))
			Save(cls);
	}
	fclose(fp);

	image=imread("testsample.png",CV_LOAD_IMAGE_GRAYSCALE);
	//计算手写数字的位置
	Position pos=GetPosition(image);
	//计算手写数字的特征
	SetFeature(image,pos);

	printf("%d\n",Fisher());
	//cvNamedWindow("test",CV_WINDOW_AUTOSIZE);
	//imshow("test",image);

	waitKey(0);
	return 0;
}