globalfunc.cpp

这是一个实现Adaboost算法的程序

// 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"

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;
}

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;
}

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((double)(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 BlockDistance(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 += fabs(*pf1 - *pf2);

	dis /= n;

	return dis;
}