fnmatrix.cpp

在人脸检测的基础之上,对嘴部的运动表情进行分析,进行语音模拟.

		rowPtrNew = a[i];
		rowPtrMat = mat[i];
		for (j = 0; j < a.xsize; j++) 
		{
			rowPtrNew[j] = rowPtrMat[j*factor];
		}
	}
	

	return a;
}
//----------------------------------------------------------------------------------
matrix UpsampleColumn(matrix& mat, int factor)
// x x --->  x 0 x 0   if factor = 2 
// x x       x 0 x 0 
{
	ASSERT( factor > 0);
	
	int i, j;
	int newxsize;
	
	
	newxsize = (int) mat.xsize*factor;
	
	matrix a(mat.ysize, newxsize, mat.numbands);
	
	double* rowPtrNew;
	double* rowPtrMat;

	for (i = 0; i < a.ysize*a.numbands; i++)
	{
		rowPtrNew = a[i];
		rowPtrMat = mat[i];
		for (j = 0; j < mat.xsize; j++) 
		{
			rowPtrNew[j*factor] = rowPtrMat[j];
		}
	}
	

	return a;
}

//----------------------------------------------------------------------------------
matrix reduce( matrix& image)
// reduce the image size by two in both directions after smoothing 
// with filterReduce()
{
	// matrix filter = filterReduce();
	
	matrix filter = filterGaussSmooth(7.0);
	matrix temp = image;
	temp = temp | filter; // filter columns
	temp = DownsampleRow( temp, 2); // downsample rows
	temp = temp | (!filter); // filter rows
	temp = DownsampleColumn( temp, 2);  // downsample columns 

	return temp; 

}

//----------------------------------------------------------------------------------
matrix expand( matrix& image)
{
	//matrix filter = filterExpand();

	matrix filter =  2*filterGaussSmooth(7.0);

	matrix temp = image;
	temp = UpsampleRow( temp, 2);
	temp = temp | filter;
	temp = UpsampleColumn( temp, 2);
	temp = temp | (!filter); 

	return temp; 
}
//-----------------------------------------------------------------------------------
matrix calcSmoothBoth( matrix& image, double spread)
{
	matrix filter = filterGaussSmooth( spread);
	matrix temp = image;

	temp = temp | filter;
	temp = temp | (!filter);

	return temp;
}

//-----------------------------------------------------------------------------------
matrix calcDerivHor( matrix& image, double spread)
{
	matrix temp = image;

	temp = temp | filterGaussSmooth( spread); //vertical smoothing
	temp = temp | (! ( -1 * filterGaussDeriv( spread) ) ); // horizontal derivative 

	return temp; 
}
//-----------------------------------------------------------------------------------
matrix calcDerivVert( matrix& image, double spread)
{
	matrix temp = image;

	temp = temp | ( -1 * filterGaussDeriv( spread) ); // vertical derivative
	temp = temp | ( ! filterGaussSmooth( spread) ); // horizontal smoothing 
	
	return temp; 
}
//-----------------------------------------------------------------------------------

//=======================================================================
// SVD from Numerical Recipes
//
// This routine computes the singular value decomposition 
// M =  U W V~ of the matrix M. The diagonal matrix W
// of singular values is output as vector w, 
// and the square matrix V (not the transpose V~) is output as matrix v.
//=======================================================================
int svdcmp(matrix& inmatrix, matrix& U_mat, matrix& diagonal, matrix& square)
{

	U_mat.resize(inmatrix.ysize, inmatrix.xsize);
	diagonal.resize(inmatrix.xsize,1);
	square.resize(inmatrix.xsize, inmatrix.xsize);

	// initialization of the parameters
	// used by the original KRC routine
	// Note that KRC uses arrays with index
	// starting from 1 instead of 0.
	int m = inmatrix.ysize; // number of rows
	int n = inmatrix.xsize; // number of columns

	// allocate memory for 1-indexed arrays
	double** a = new double*[ m + 1];

	for( int M = 0; M <= m; M++) 
		a[M] = new double[ n + 1];

	double** v = new double*[ n + 1];

	for( int N = 0; N <= n; N++) 
		v[N] = new double[ n + 1];

	double* w = new double[ n + 1];

	// copy the entries of the input matrix
	// into the 1-indexed array
	for( M = 1; M <= m; M++)
	{
		for( N = 1; N <= n; N++)
		{
			a[M][N] = inmatrix(M-1,N-1);
		}
	}
	
	// ORIGINAL ROUTINE STARTS HERE
	int flag,i,its,j,jj,k,l,nm;
	double anorm,c,f,g,h,s,scale,x,y,z,*rv1;

	//rv1=dvector(1,n);
	rv1 = new double[ n + 1];
	g=scale=anorm=0.0;
	for (i=1;i<=n;i++) {
		l=i+1;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i <= m) {
			for (k=i;k<=m;k++) scale += fabs(a[k][i]);
			if (scale) {
				for (k=i;k<=m;k++) {
					a[k][i] /= scale;
					s += a[k][i]*a[k][i];
				}
				f=a[i][i];
				g = -fMath.Sign(sqrt(s),f);
				h=f*g-s;
				a[i][i]=f-g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
					f=s/h;
					for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
				}
				for (k=i;k<=m;k++) a[k][i] *= scale;
			}
		}
		w[i]=scale *g;
		g=s=scale=0.0;
		if (i <= m && i != n) {
			for (k=l;k<=n;k++) scale += fabs(a[i][k]);
			if (scale) {
				for (k=l;k<=n;k++) {
					a[i][k] /= scale;
					s += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g = -fMath.Sign(sqrt(s),f);
				h=f*g-s;
				a[i][l]=f-g;
				for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
				for (j=l;j<=m;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
					for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
				}
				for (k=l;k<=n;k++) a[i][k] *= scale;
			}
		}
		anorm=fMath.Max(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n;i>=1;i--) {
		if (i < n) {
			if (g) {
				for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
					for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
				}
			}
			for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
		}
		v[i][i]=1.0;
		g=rv1[i];
		l=i;
	}
	for (i=fMath.Min(m,n);i>=1;i--) {
		l=i+1;
		g=w[i];
		for (j=l;j<=n;j++) a[i][j]=0.0;
		if (g) {
			g=1.0/g;
			for (j=l;j<=n;j++) {
				for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
				f=(s/a[i][i])*g;
				for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
			}
			for (j=i;j<=m;j++) a[j][i] *= g;
		} else for (j=i;j<=m;j++) a[j][i]=0.0;
		++a[i][i];
	}
	for (k=n;k>=1;k--) {
		for (its=1;its<=30;its++) {
			flag=1;
			for (l=k;l>=1;l--) {
				nm=l-1;
				if ((double)(fabs(rv1[l])+anorm) == anorm) {
					flag=0;
					break;
				}
				if ((double)(fabs(w[nm])+anorm) == anorm) break;
			}
			if (flag) {
				c=0.0;
				s=1.0;
				for (i=l;i<=k;i++) {
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if ((double)(fabs(f)+anorm) == anorm) break;
					g=w[i];
					h=fMath.dpythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=1;j<=m;j++) {
						y=a[j][nm];
						z=a[j][i];
						a[j][nm]=y*c+z*s;
						a[j][i]=z*c-y*s;
					}
				}
			}
			z=w[k];
			if (l == k) {
				if (z < 0.0) {
					w[k] = -z;
					for (j=1;j<=n;j++) v[j][k] = -v[j][k];
				}
				break;
			}
			//nrerror("no convergence in 30 dsvdcmp iterations");
			if(its == 30)
			{
				return 0;
			}
			x=w[l];
			nm=k-1;
			y=w[nm];
			g=rv1[nm];
			h=rv1[k];
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
			g=fMath.dpythag(f,1.0);
			f=((x-z)*(x+z)+h*((y/(f+fMath.Sign(g,f)))-h))/x;
			c=s=1.0;
			for (j=l;j<=nm;j++) {
				i=j+1;
				g=rv1[i];
				y=w[i];
				h=s*g;
				g=c*g;
				z=fMath.dpythag(f,h);
				rv1[j]=z;
				c=f/z;
				s=h/z;
				f=x*c+g*s;
				g = g*c-x*s;
				h=y*s;
				y *= c;
				for (jj=1;jj<=n;jj++) {
					x=v[jj][j];
					z=v[jj][i];
					v[jj][j]=x*c+z*s;
					v[jj][i]=z*c-x*s;
				}
				z=fMath.dpythag(f,h);
				w[j]=z;
				if (z) {
					z=1.0/z;
					c=f*z;
					s=h*z;
				}
				f=c*g+s*y;
				x=c*y-s*g;
				for (jj=1;jj<=m;jj++) {
					y=a[jj][j];
					z=a[jj][i];
					a[jj][j]=y*c+z*s;
					a[jj][i]=z*c-y*s;
				}
			}
			rv1[l]=0.0;
			rv1[k]=f;
			w[k]=x;
		}
	}
	//free_dvector(rv1,1,n);
	delete[] rv1;

	//ORIGINAL KRC ROUTINE ENDS HERE

	// copy back from the 1-indexed arrays
	// and deallocate memory
	for( M = 1; M <= m; M++)
	{
		for( N = 1; N <= n; N++)
		{
			U_mat(M-1, N-1) = a[M][N];
		}
		delete[] a[M];
	}
	delete[] a[0];
	delete[] a;

	for( int NN = 1; NN <= n; NN++)
	{
		diagonal(NN-1) = w[NN];
		for( N = 1; N <= n; N++)
		{
			square(NN-1, N-1) = v[NN][N];
		}
		delete[] v[NN];
	}
	delete[] v[0];
	delete[] v;
	delete[] w;

	return 1;
}
//---------------------------------------------------------------------
matrix inverse(matrix& inMatrix)
{
	// calculates the inverse of a matrix using the singular 
	// value decomposition of the input matrix inMatrix 

	ASSERT( inMatrix.ysize == inMatrix.xsize && inMatrix.xsize > 0); 

	matrix U, S, V, T;
	register int i, j;
	matrix temp; 

	if ( svdcmp(inMatrix, U, S, V))
	{
		ASSERT ( U.xsize == S.ysize && S.ysize == V.xsize && S.ysize == V.ysize);

		for ( i = 0; i < S.ysize; i++)
			if ( S(i) == 0.0)
			{
				CWnd mess;
				mess.MessageBox("Singular matrix! Cannot invert.");
				return temp; 
			}

		T = zeros( U.xsize, U.ysize);
		for ( j = 0; j < T.ysize; j++)
		{
			for ( i = 0; i < T.xsize; i++)
			{
				T(j, i) = (1/S(j)) * U(i, j);
			}
		}
		return ( V*T);
	}
	else
	{
		CWnd mess;
		mess.MessageBox("Cannot invert matrix."); 	
		return ( temp); 
	}
}
//---------------------------------------------------------------------
matrix OuterProduct(matrix& a, matrix& b)
{
	ASSERT( a.xsize == 1 && b.ysize == 1);

	matrix result(a.ysize, b.xsize);
	
	for ( int j = 0; j < result.ysize; j++)
	{
		for ( int i = 0; i < result.xsize; i++)
		{
			result(j, i) = a(j) * b(i);
		}
	}

	return result; 
}

//---------------------------------------------------------------------

void ConvertMatrixTo2DArray( matrix& inImage, double** outImage)
{
	// Converts the image stored in the matrix inImage to a 
	// 2D array to be stored in outImage

	register int i, j;

	for ( j = 0; j < inImage.ysize; j++)
	{
		for ( i = 0; i< inImage.xsize; i++)
		{
			outImage[j][i] = inImage(j, i);
		} // for i
	} // for j
}
//---------------------------------------------------------------------
matrix Convert2DArrayToMatrix( double** inImage, int ysi, int xsi)
{
	// Converts the image stored in the 2D array  inImage to a 
	// matrix to be stored in outImage, and returns it. 

	register int i, j;

	matrix outImage(ysi, xsi);

	for ( j = 0; j < ysi; j++)
	{
		for ( i = 0; i < xsi; i++)
		{
			outImage(j, i) = inImage[j][i];
		} // for i
	} // for j
	return outImage;
}
//---------------------------------------------------------------------
double determinant(matrix& a)
{
	ASSERT(a.numbands == 1);

	double det = 0;

	if (a.xsize == 3 && a.ysize ==3)
	{
		det += a(0, 0) * ( a(1, 1) * a(2,2) - a(1,2) * a(2,1) );
		det += -a(0,1) * ( a(1, 0) * a(2,2) - a(1,2) * a(2,0) );
		det += a(0, 2) * ( a(1, 0) * a(2,1) - a(1,1) * a(2,0) );
	}
	return det;
}
//---------------------------------------------------------------------