wi.cpp

该程序用于在线数字签名,可以将读入的图像与用户开户时录入的图像比较,从而确定身份的合法性.

/* 功能：求出签名每个数据点的速度
/* 参数：s－输入签名样本
/****************************************************************************/
void CWi::SpeedCurve(CWi::Sample* s)
{
	int num = s->effectDataNum;
	int subx,suby;
	double dis;
	for(int i = 0; i < num - 1; i++)
	{
		subx = s->point[i+1].x - s->point[i].x;
		suby = s->point[i+1].y - s->point[i].y;
		dis = SQUARE(subx) + SQUARE(suby);
		dis = sqrt(dis);
		s->speed[i] = dis;
	}
	s->speed[num-1] = s->speed[num-2];
	NormalizeSample(s->speed, s->effectDataNum);
}

/****************************************************************************/
/* 功能：将签名输入数据转换为矢量量化所用的向量
/* 参数：sampleArrayRef－输入签名样本
/****************************************************************************/
void CWi::GenVector(CPtrArray* sampleArrayRef)
{
	int sampleNum;
	int i,count;
 	sampleNum = sampleArrayRef->GetSize();

	for (count = 0; count < sampleNum; count++)
	{
		struct Sample *s;
		s = (Sample *)sampleArrayRef->GetAt (count);
		SpeedCurve(s);
		for(i = 0; i < s->effectDataNum; i++)
		{
			s->dataVec[i].num = 4;
			s->dataVec[i].data[0] = s->pressTotal2[i];
			s->dataVec[i].data[1] = s->xLoc2[i];
			s->dataVec[i].data[2] = s->yLoc2[i];
			s->dataVec[i].data[3] = s->speed[i];
		}
	}
}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::GenVQVector(CPtrArray* sampleArrayRef)
{
	int sampleNum;
	int i,count;
 	sampleNum = sampleArrayRef->GetSize();

	for (count = 0; count < sampleNum; count++)
	{
		struct Sample *s;
		s = (Sample *)sampleArrayRef->GetAt (count);
		SpeedCurve(s);
		for(i = 0; i < s->effectDataNum; i++)
		{
			s->vq_vector[i].Data = (double *)malloc(sizeof(double)*4);
			s->vq_vector[i].nDimension = 4;
			s->vq_vector[i].Data[0] = s->pressTotal2[i];
			s->vq_vector[i].Data[1] = s->xLoc2[i];
			s->vq_vector[i].Data[2] = s->yLoc2[i];
			s->vq_vector[i].Data[3] = s->speed[i];
			s->vq_vector[i].nCluster = 0;     //聚类结果初值 0
		}
	}

}

/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::GenVQVector()
{
	int i;
	for(i = 0; i < templenvec; i++)
	{
		vq_vector[i].Data = (double *)malloc(sizeof(double)*4);
		vq_vector[i].nDimension = 4;
		vq_vector[i].Data[0] = dtwtemplatev[i].data[0];
		vq_vector[i].Data[1] = dtwtemplatev[i].data[1];
		vq_vector[i].Data[2] = dtwtemplatev[i].data[2];
		vq_vector[i].Data[3] = dtwtemplatev[i].data[3];
		vq_vector[i].nCluster = 0;     //聚类结果初值 0
	}

}

/****************************************************************************/
/* 功能：生产VQ所需的模板
/****************************************************************************/
void CWi::GenVQTemplate()
{
	GenTemplateVector();
	GenVQVector();
	VQCluster();
	GenHMM();
	int num = trainSampleArray.GetSize();
	double posi;
	double min_pos = 0.0;
	double max_pos = -10000.0;
	for(int i = 0; i < num; i++)
	{
		CWi::Sample *s = (CWi::Sample *)trainSampleArray.GetAt(i);
		posi = PoHMM(s);
		if((posi < min_pos) && (posi > -10000.0))
			min_pos = posi;
		if(posi > max_pos)
			max_pos = posi;
	}
	threshold = thresholdhmm = min_pos;// - (max_pos - min_pos) / 2;
}

/****************************************************************************/
/* 功能：矢量量化
/****************************************************************************/
void CWi::VQCluster()
{
	int i;
	for(i = 0; i < VQ_M; i++)
	{
		vq_center[i].nDimension = 1;
		vq_center[i].Data = (double *)malloc(sizeof(double) * 1);
		vq_center[i].Num = 0;	   //初始化为 0
	}
	LBGCluster(vq_vector,templenvec , vq_center, VQ_M);
	for(i = 0; i < templenvec; i++)
	{
		vq_vector[i].nCluster++;
		vq_feature[i] = vq_vector[i].nCluster;
	}
}

/****************************************************************************/
/*功能：训练HMM模型
/****************************************************************************/
void CWi::GenHMM()
{
	int niter;
	double logprobinit, logprobfinal;

	double **alpha = dmatrix(1, templenvec, 1, hmm.N);
	double **beta = dmatrix(1, templenvec, 1, hmm.N);
	double **gamma = dmatrix(1, templenvec, 1, hmm.N);

	BaumWelch(&hmm, templenvec, vq_feature, alpha, beta, gamma, &niter, &logprobinit, &logprobfinal);

	free_dmatrix(alpha, 1, templenvec, 1, hmm.N);
	free_dmatrix(beta, 1, templenvec, 1, hmm.N);
	free_dmatrix(gamma, 1, templenvec, 1, hmm.N);
}

/****************************************************************************/
/*功能：利用HMM模型进行识别
/*参数：distance－保存测试样本与模板的距离
/*      index－测试样本的下标
/****************************************************************************/
int CWi::RecogHMM(double *distance,int index)
{
	CWi::Sample *s = (CWi::Sample *)testSampleArray.GetAt(index);
	double pprob = PoHMM(s);
	*distance = pprob;
	if(pprob > thresholdhmm)
	{
		return 1;
	}
	else
	{
		return -1;
	}
}

/****************************************************************************/
/*功能：计算测试样本的出现概率
/*参数：s－测试样本
/****************************************************************************/
double CWi::PoHMM(CWi::Sample *s)
{
	VQClassify(s);
	double pprob;
	int *q = ivector(1,s->effectDataNum);
	double **delta = dmatrix(1, s->effectDataNum, 1, hmm.N);
	int **psi = imatrix(1, s->effectDataNum, 1, hmm.N);
	Viterbi(&hmm, s->effectDataNum, s->vq_feature, delta, psi,q, &pprob);
	if(pprob != 0)
		pprob = log(pprob)/(s->effectDataNum * log(VQ_M));
	else 
		pprob = -10000.0;
	return pprob;
}

/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::VQClassify(CWi::Sample *s)
{
	int i;
	for(i = 0; i < s->effectDataNum; i++)
	{
		s->vq_vector[i].nCluster = VQClassify(&vq_vector[i], vq_center);
		s->vq_feature[i] = s->vq_vector[i].nCluster;
	}

}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

int CWi::VQClassify(struct VQ_VECTOR *vq, struct VQ_CENTER *vc)
{
	int i,j,min_index = -1;
	double distance = 0.0,temp = 0.0;
	double min_dis = 1000000.0;			//初始值
	for(i = 0; i < VQ_M; i++)
	{	
		distance = 0.0;
		for(j = 0; j < vq->nDimension; j++)
		{
			temp = vc[i].Data[j] - vq->Data[j];
			temp = SQUARE(temp);
			distance += temp;
		}
		distance = sqrt(distance);
		if(distance < min_dis)
		{
			min_dis = distance;
			min_index = i;
		}

	}

	min_index++;
	return min_index;
}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::SamplePress1(double *src,double *dst,int datalen,int inter)
{
	int i;
	for(i = 0;i<datalen;i++)
	{
		dst[i]= src[i];
	}
}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::SamplePress2(double *src,double *dst,int datalen,int inter)
{
	int i,j;
	double temp=0.0;
	for(i = 0;i<datalen-1;i++)
	{
		for(j = 0; j<inter ;j++)
		{
			temp+=src[i*inter+j];
		}
		temp/=inter;
		dst[i] = temp;
		temp = 0.0;

	}
	dst[datalen-1] = src[datalen-1]; 

}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

//提取波峰，波谷的信息 threshold 为阈值
int CWi::WaveCrest(double *src,double *dst,int *pos,int datalen,int threshold)
{
	int i,count=0;
	int num =0;
	double temp1=0.0,temp2=0.0;
	double *diff = (double *)malloc(sizeof(double)*datalen);
 	DiffSample (src, diff, datalen);	//微分
	for(i = 0;i<datalen-1;i++)
	{
		if(diff[i]*diff[i+1]<=0)
		{
			pos[count] = i;
			count++;
		}
	}
	for(i = 0;i<count;i++)
	{
		if((src[(pos[i]-threshold)>0?(pos[i]-threshold):0]-src[pos[i]])
			*(src[(pos[i]+threshold)<datalen?(pos[i]+threshold):datalen]-src[pos[i]])<0)
		{
			pos[i] = -1;
			continue;
		}
		dst[num] = src[pos[i]];
		num++;

	}

	return num;
}

/****************************************************************************/
/* 功能：旋转
/* 参数：posX－X坐标序列
/*       posY－Y坐标序列
/*       datalen－序列的长度
/*       rotAngle－序列的长度
/****************************************************************************/
void CWi::Rotation(double *posX, double *posY,int dataLen,double rotAngle)
{
	int i;
	double tempX,tempY;
	double cosa = cos(rotAngle);
	double sina = sin(rotAngle);
	for(i = 0; i<dataLen;i++)
	{
		tempX = posX[i]*cosa+posY[i]*sina;
		tempY = posY[i]*cosa-posX[i]*sina;
		posX[i] = tempX;
		posY[i] = tempY;
	}
}

/****************************************************************************/
/*功能：中值滤波
/*参数：data－输入数组
/*      datalen－数组长度
/****************************************************************************/
void CWi::MedianFilter(double *data,int dataLen)
{

	int i;
	double temp;
	for(i = 2; i<dataLen-2;i++)
	{
		temp = (data[i-2]+data[i-1]+data[i+1]+data[i+2])/4;
		if(ABS(data[i]-temp)>1.0)
			data[i]=temp;
	}
}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::ProjectionFeature(CWi::Sample* s, int num)
{
	int i,j,fea,count = 0;
	double hor[IMAGEWIDTH],ver[IMAGEHEIGHT],total=0.0,piece=0.0;
	s->imgPixel[0][0]=0.0;
	for(i = 0;i<IMAGEWIDTH;i++)
	{
		hor[i] = 0.0;
	}
	for(i = 0;i<IMAGEHEIGHT;i++)
	{
		ver[i] = 0.0;
	}
	for(i = 0;i<IMAGEWIDTH;i++)
	{
		for(j = 0;j<IMAGEHEIGHT;j++)
		{
			hor[i]+=(255.0-s->imgPixel[i][j]);
			ver[j]+=(255.0-s->imgPixel[i][j]);
			total+=(255.0-s->imgPixel[i][j]);
		}
	}
	total/=num;
	fea = 0;
	for(i = 0;i<IMAGEWIDTH;i++)
	{	
		count++;
		piece+=hor[i];
		if(piece>=total)
		{
			s->feature[fea]=count;
			fea++;
			count = 0;
			piece=0;
		}
	}
	if(piece<total)
		s->feature[fea]=count;
	fea++;
	count = 0;
	piece = 0;
	for(i = 0;i<IMAGEHEIGHT;i++)
	{	
		count++;
		piece+=ver[i];
		if(piece>=total)
		{
			s->feature[fea]=count;
			fea++;
			count = 0;
			piece = 0;
		}
	}
	if(piece<total)
		s->feature[fea]=count;

}

/****************************************************************************/
/*功能：计算方向分布
/*参数：angle－保存方向计算的结果
/*      xLog－X坐标序列
/*      yLog－Y坐标序列
/*      datalen－序列长度
/****************************************************************************/
void CWi::CalDirection(double *angle,double *xLoc,double *yLoc,int dataLen)
{	
	int i;
	double *diff0 = (double *)malloc(sizeof(double)*dataLen);
	double *diff1 = (double *)malloc(sizeof(double)*dataLen);
	DiffSample (xLoc, diff0, dataLen);	//微分
	DiffSample (yLoc, diff1, dataLen);
	for(i = 0;i<dataLen;i++)
	{
		if(diff0[i] == 0)
		{
			if(diff1[i]>0)
				angle[i] = 90+90;
			if(diff1[i]<0)
				angle[i] = -90+90;
		}
		else
			angle[i] = 180*(atan(diff1[i]/diff0[i])/PI)+90;
	}
}
/****************************************************************************/
/*
/*	
/*  
/****************************************************************************/

void CWi::Direction(int sample)
{

	CWi::Sample *s0 =(CWi::Sample*)trainSampleArray.GetAt(sample);

	CalDirection(s0->direction,s0->xLoc2,s0->yLoc2,s0->effectDataNum);
}

/****************************************************************************/
/*功能：低通滤波
/*参数：data－输入数据
/*      threshold－阈值
/*      datalen－数组长度
/****************************************************************************/
void CWi:: lowpassfilter(double *data,double threshold,unsigned long dataLen)
{
	int r=0;
	unsigned long i,mask=0xffffffff;
	while(mask&dataLen)
	{
		mask<<=1;
		r++;
	}
	unsigned long count = 1<<r;
	complex *com1,*com2;
	com1 = (complex *)malloc(sizeof(complex)*count);
	com2 = (complex *)malloc(sizeof(complex)*count);
	//Initiallize
	for (i = 0; i < count; i++)
	{
		com1[i].r = i>=dataLen?0:data[i];
		com1[i].i = 0.0;
	}

	fft(com1,com2, r);
	
	for(i = (int)((threshold/2)*count); i< count; i++)
	{
		com2[i].r = 0.0;
	}

	ifft (com2,com1, r);