globalfunc.cpp

图像标注程序

// GlobalFunc.cpp: implementation of the CGlobalFunc class.
//
//////////////////////////////////////////////////////////////////////

/*
   This file is created by Shiguang Shan at 01.06.2002 to contain 
   Some common-used global functions
*/

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


#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

#define MAX_GRAY 256


int		gnNormalFaceWidth = 32;
int		gnNormalFaceHeight= 32;
CString	gWorkDir;
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CGlobalFunc::CGlobalFunc()
{

}

CGlobalFunc::~CGlobalFunc()
{

}

/* 

  *******************************************************	
  The following functions are used for Matrix Compuation
  Induced in 2002.01.04 by Shiguang Shan
  *******************************************************
  
*/

////////////////////////////////////////////////////////////////
//对矩阵进行按行置换的LU分解，结果存在该矩阵中
BOOL CGlobalFunc::ludcmp(double *a, int n, int *indx, double *d)
{
	int      i,imax,j,k;
	double   big,dum,sum,temp;
	double   *vv;
	double   swap;
	
	vv = new double[n];
	*d = 1.0;
	for(i=0; i<n; i++)
	{
		big = 0.0;
		for(j=0; j<n; j++)
			if((temp = fabs( *(a+i*n+j))) > big)  big =temp;
			if(big == 0.0) 
			{
				AfxMessageBox("Singular matrix in routine ludcmp");
				delete [] vv;
				return false;
			}
			vv[i] = 1.0/big;
	}
	
	for(j=0; j<n; j++)
	{
		for(i=0; i<j; i++)
		{
			sum = *(a+i*n+j);
			for(k=0; k<i; k++)  sum -= (*(a+i*n+k))*(*(a+k*n+j));
			*(a+i*n+j) = sum;
		}
		
		big = 0.0;
		for(i=j; i<n; i++)
		{
			sum = *(a+i*n+j);
			for(k=0; k<j; k++)
				sum -= (*(a+i*n+k))*(*(a+k*n+j));
			*(a+i*n+j) = sum;
			if((dum=vv[i] * fabs(sum)) >= big)
			{
				big = dum;
				imax = i;
			}
		}
		
		if(j != imax)
		{
			for(k=0; k<n; k++)
			{
				dum = *(a+imax*n+k);
				*(a+imax*n+k) = *(a+j*n+k);
				*(a+j*n+k) = dum;
			}
			*d = -(*d);
			swap = vv[imax];
			vv[imax] = vv[j];
			vv[j] = swap;
		}
		indx[j] = imax;
		if(*(a+j*n+j) == 0.0)  *(a+j*n+j) = 1.0e-20;
		if(j != n -1 )
		{
			dum = 1.0/(*(a+j*n+j));
			for(i=j+1; i<n; i++)  *(a+i*n+j) *= dum;
		}
	}
	delete [] vv;
	
	return  true;
}

//////////////////////////////////////////////////////////////////
//求解n维线性方程A*x=B
void CGlobalFunc::lubksb(double *a, int n, int *indx, double b[])
{
	int      i, ii=-1, ip, j;
	double   sum;
	
	for(i=0; i<n; i++)
	{
		ip = indx[i];
		sum = b[ip];
		b[ip] = b[i];
		if(ii>=0)
			for(j=ii; j<=i-1; j++)  sum -= (*(a+i*n+j))*(b[j]);
			else if  (sum) ii = i;
			b[i] = sum;
	}
	
	for(i=n-1; i>=0; i--)
	{
		sum = b[i];
		for(j=i+1; j<n; j++)  sum -= (*(a+i*n+j))*(b[j]);
		b[i] = sum/(*(a+i*n+i));
	}
}

//////////////////////////////////////////////////////////////////////
// In MatrixComputReverse : a will be distroyed
BOOL CGlobalFunc::MatrixComputReverse(double *a, double *y,int n)
{
	
	double   d, *col;
	int      i, j, *indx;
	
	col = new double[n];
	indx = new int[n];
	
	if (!ludcmp(a, n, indx, &d)) 
	{
		delete [] col;
		delete [] indx;
		return false;
	}
	
	for(j=0; j<n; j++)
	{
		for(i=0; i<n; i++)  col[i] = 0.0;
		
		col[j] = 1.0;
		lubksb(a, n, indx,col);
		for(i=0; i<n; i++)  *(y+i*n+j) = col[i];
	}
	
	delete [] col;
	delete [] indx;
	
	return true;
	
}


/////////////////////////////////////////////////////////////////////////////////////
//comput MatrixMultiply
BOOL	CGlobalFunc::MatrixMultiply(double* pScr1, double* pScr2, double* pDst, int m, int n, int l, int p)
{
	int		i, j, k;
	
	ASSERT(n == l);
	
	double	sum = 0.0;
	double* pDstMatrix = pDst;
	
	for(i=0; i<m; i++)
	{
		for(j=0; j<p; j++)
		{
			sum = 0;
			for(k=0; k<n; k++)
			{
				sum += (*(pScr1+i*n+k)) * (*(pScr2+k*p+j));
			}
			*pDstMatrix++ = sum;
		}
	}
	
	return TRUE;
	
}

void CGlobalFunc::ScaleVector(double * src,double * dest, int length, double value)
{
	int i;
	double val = value;
	
	for( i = 0; i < length; i++ )
	{
		dest[i] = src[i] * val;
	}
	
}

void CGlobalFunc::SubVector(double * src1,double * src2, double * dest, int length) 
{
	int i;
	
	for( i = 0; i < length; i++ )
	{
		dest[i] = src1[i] - src2[i];
	}
	
}

/*

  Note: 
  pA: should be a n-row, m-col matrix, vectorized by row-first order
  
*/

void	CGlobalFunc::PseudoInverseMatrix(double *pA, double* pA_PT, long n, long m)
{
	int     i,j;
	
	double 	*p, *p2;
	p = pA;
	
	//compute the reverse matrix of A
	double  *pAT = new double[m*n];
	p = pAT; p2 = pA;
	for(i=0; i<m; i++)
		for(j=0; j<n; j++)
			*(p+i*n+j) = *(p2+j*m+i);
		
	//multiply A' and A
	double   *pATA = new double[m*m];
	MatrixMultiply(pAT, pA, pATA, m, n, n, m);
	
	//compute matrix D = (A'A)~1
	double*	pD = new double[m*m];
	MatrixComputReverse(pATA, pD, m);
	
	//compute matrix the pseudo-reverse of A pA_PT = D * A'
	MatrixMultiply(pD, pAT, pA_PT, m, m, m, n);
	
	delete[] pAT;
	delete[] pATA;
	delete[] pD;

}



float ran1(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum <= 0 || !iy) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if (*idum < 0) *idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = *idum;
	if ((temp=(float)(AM*iy)) > RNMX) return (float)RNMX;
	else return temp;
}


float gasdev(long *idum)
{
	
	ran1(idum);
	static int iset=0;
	static float gset;
	float fac,rsq,v1,v2;

	if  (iset == 0) {
		do {
			v1=(float)(2.0*ran1(idum)-1.0);
			v2=(float)(2.0*ran1(idum)-1.0);
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=(float)sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}

int		Convert2Stride(int	nWidth)
{
	 return ((nWidth%4)?(nWidth+4-nWidth%4):(nWidth));
}

/* To normalize the vector to be length = 1 */
double NormalVector2UnitLength(double* vector, long N)
{
	long j;
	double sum = 0.0;
    double  *pv = vector;
	
	for(j=0; j <N; j++, pv++)
		sum += (*pv) * (*pv);
	sum = sqrt(sum); // it is ||vector|| 
	
	pv  = vector;
	for(j=0; j<N; j++)
		*pv++ /= sum;  

	return sum;
}

/* 
   First, To Normalize the vector to be 0-mean and 1-variance;
   Then,  Normalize the vector to length 1
*/
double NormalVector2ZeroMeanUnitVar(double* vector, long N)
{
	long j;
	double std=0.0, mean = 0.0;
    double  *pv;

	pv = vector;
	for(j=0; j<N; j++)
		mean += *pv++;
	mean /= N;
	
	pv  = vector;
	for(j=0; j<N; j++, pv++)
		std += (*pv-mean)*(*pv-mean);
	std = sqrt(std/(N-1));

	pv  = vector;
	for(j=0; j<N; j++, pv++)
		*pv = (*pv-mean)/std;

//	pv  = vector;
//	NormalVector2UnitLength(pv, N);//Still need to normalize to length '1'

	return mean;
}

int MyStringToInt(CString strNum)
{
	int n = strNum.GetLength();
	int begin = 0;

	for (int i = 0; i < n; i++) {
		if ('0' != strNum[i]) {
			begin = i;
			break;
		}
	}

	if (n == i) return 0;

	return atoi(strNum.Mid(begin, n - begin));
}
/*
double InnerProduct(double* pv1,double* pv2, long N) 
{
	long   i;
	double  sum=(double)0.0;

	double* q1=pv1;
	double* q2=pv2;
	
	for(i=0;i<N;i++)
		sum += (*q1++) * (*q2++);
	
	return sum;
}
*/

double PointDistance(CPoint p1, CPoint p2)
{
	return (double) sqrt((p2.y-p1.y) * (p2.y-p1.y) + (p2.x-p1.x) * (p2.x-p1.x));
}

double EuclidDistance(double* pv1, double* pv2, long N)
{
	long	i;
	double	dis = 0.0;

	double *pf1 = pv1;
	double *pf2 = pv2;

	for(i=0; i<N; i++, pf1++, pf2++)
		dis += (*pf1 - *pf2) * (*pf1 - *pf2);

	dis = (double)sqrt(dis);
	return dis;
}

double InnerProduct(const double* pv1,const double* pv2, long N) 
{
	long   i;
	double  sum=(double)0.0;

	const double* q1=pv1;
	const double* q2=pv2;
	
	for(i=0;i<N;i++)
		sum += (*q1++) * (*q2++);
	
	return sum;
}

float InnerProduct(const float* pv1,const float* pv2, long N) 
{
	long   i;
	float  sum=0.0;

	const float* q1=pv1;
	const float* q2=pv2;
	
	for(i=0;i<N;i++)
		sum += (*q1++) * (*q2++);
	
	return sum;
}
double Chi_Square(double *pInVec1, double *pInVec2, int nDim)
{

	// choose the min as the fit
/*
	double d = 0.0;
	double dTmpSum, dTmpSub, dTmp;
	double *p1;
	double *p2;

	p1 = pInVec1;
	p2 = pInVec2;

//	NormalizeHistogramLbp(p1, nDim);
//	NormalizeHistogramLbp(p2, nDim);

	for (int i=0; i<nDim; i++) {
		dTmpSum = *p1 + *p2;

		if (dTmpSum > 0.0000001) {
			dTmpSub = *p1 - *p2;
			dTmp = dTmpSub*dTmpSub;
			d += dTmp/dTmpSum;
		}
		p1++; p2++;
	}
*/

	double d = 0.0;
	double dTmpSum, dTmpSub, dTmp;

//	NormalizeHistogramLbp(p1, nDim);
//	NormalizeHistogramLbp(p2, nDim);

	for (int i=0; i<nDim; i++) {
		dTmpSum = pInVec1[i] + pInVec2[i];

		if (dTmpSum > 0.0000001) {
			dTmpSub = pInVec1[i] - pInVec2[i];
			dTmp = dTmpSub*dTmpSub;
			d += dTmp/dTmpSum;
		}
	}

	return d;
}
double Histogram_Intersection(double *pInVec1, double *pInVec2, int nDim)
{
	double d = 0.0;
	double *p1;
	double *p2;

	p1 = pInVec1;
	p2 = pInVec2;

	int i;
	for (i=0; i<nDim; i++) {
		if (p1[i] > p2[i]) {
			d += p2[i];
		}
		else {
			d += p1[i];
		}
	}

	return d;
}
double Histogram_Intersection(BYTE *pInVec1, BYTE *pInVec2, int nDim)
{
	double d = 0.0;
	BYTE *p1;
	BYTE *p2;

	p1 = pInVec1;
	p2 = pInVec2;

	int i;
	for (i=0; i<nDim; i++) {
		if (p1[i] > p2[i]) {
			d += p2[i];
		}
		else {
			d += p1[i];
		}
	}

	return d;
}
double Histogram_Intersection(int *pInVec1, int *pInVec2, int nDim)
{
	double d = 0.0;
	int *p1;
	int *p2;

	p1 = pInVec1;
	p2 = pInVec2;

	int i;
	for (i=0; i<nDim; i++) {
		if (p1[i] > p2[i]) {
			d += p2[i];
		}
		else {
			d += p1[i];
		}