fnmatrix.cpp

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

						multp = b(0,k+strt)*a(j ,indx, n);
						outfilt(j,i,n) += multp;
					}
				}
			}
		}
	} 
	else 
	{
		// make 2D convolution
		for( n = 0; n < numbands; n++)
		{

			for ( j=0; j<ysize; j++) 
			{
				for ( i=0; i<xsize; i++)
				{
					for ( k=-strt1; k<=strt1; k++) 
					{
						jndx = j-k;
						if(jndx<0)
						{
							jndx = 0;
						}
						else if(jndx >= a.ysize)
						{
							jndx = a.ysize - 1;
						}
						for ( l=-strt2; l<=strt2; l++) 
						{
							indx = i-l;
							if(indx<0)
							{
								indx = 0;
							}
							else if(indx >= a.xsize)
							{
								indx = a.xsize - 1;
							}
							multp = b(k+strt1,l+strt2)*a(jndx,indx,n);
							outfilt(j,i,n) += multp;
						}
					}
				}
			}
		}
	}
	//  cerr<<"\nstrt = "<<strt<<"\n";
	return(outfilt);
}

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

matrix operator!(const matrix &in)
{
	ASSERT( in.numbands == 1);
	return(trps(in));
}

//---------------------------
//     Special matrices
//---------------------------
matrix zeros(int k, int l, int n)
{
	matrix res(k, l, n);
	return(res);
}
//----------------------------------------------------------------------------------

matrix eye(int k)
{
	matrix res(k,k);
	int m;
	for(m=0; m<k; m++) 
		res[m][m]=1.0;
	return(res);
}
//----------------------------------------------------------------------------------

matrix	ones(int k, int l, int n)
{
	matrix res(k, l, n);
	register double *t1;
	int i;
	for(i = 0, t1 = res.elem; i < k * l * n; i++) 
		*t1++ = 1.0;
	return(res);
}

//-----------------------------------
// row/column, cut/paste, diag
//----------------------------------- 
matrix	cutr(const matrix& a, int k)   
//cut k'th row and assign to a row vector
{
	ASSERT( a.numbands == 1 );
	matrix res(1,a.xsize);
	register double *t1, *t2;
	int i;
	for(i=0, t1=res.elem, t2=a.rowptr[k]; i<a.xsize; i++) 
		*t1++=*t2++;
	return(res);
}
//----------------------------------------------------------------------------------
matrix	cutc(const matrix& a, int k)
//cut k'th column and assign to a column vector
{
	ASSERT( a.numbands == 1 );
	
	matrix res(a.ysize,1);
	register double *t1, *t2;
	int i;
	for(i=0, t1=res.elem, t2=a.elem+k; i<a.ysize; i++, t2+=a.xsize) 
		*t1++=*t2;
	return(res);
}
//----------------------------------------------------------------------------------

matrix	diag(const matrix& a)
//creates or extracts main diagonal
{
	ASSERT( a.numbands == 1 );
	matrix res;
	if ((a.ysize >1) && (a.xsize>1))  // a is matrix
	{
		int i,len=(a.ysize < a.xsize) ? a.ysize : a.xsize;
		matrix res1(len,1);
		for(i=0;i<len;i++) res1(i)=a.rowptr[i][i];
		res=res1;
	}
	else if ((a.ysize !=0) && (a.xsize!=0))  //a is vector
	{
		int i, len=(a.ysize > a.xsize) ? a.ysize : a.xsize;
		matrix res2(len,len);
		for(i=0;i<len;i++) res2[i][i]=a.elem[i];
		res=res2;
	}
	return(res);
}
//----------------------------------------------------------------------------------
matrix 	sort(const matrix& A)
//sorts a the double elements of a vector in descending order of
//magnitude (abs value of a double number).
{
	ASSERT( A.numbands == 1 );
	
	int 	tmpint, maxk, i, k, l,
		P=(A.ysize > A.xsize) ? A.ysize : A.xsize,
		*indexx=new int [P];
	for (i=0;i<P;i++) indexx[i]=i;
	double comtmp;
	double max;
	matrix B=A;
	for (i=0;i<P;i++)
	{
		max = fabs(B.elem[i]);
		for (k=i+1; k<P;k++)
			if (max < fabs(B.elem[k]))
			{
				maxk=k;
				max=fabs(B.elem[k]);
			}
			tmpint=indexx[maxk];
			indexx[maxk]=indexx[i];
			indexx[i]=tmpint;
			
			//			cout << "i=" << i << "\n";
			for (l=0;l<P;l++) printf("%d", indexx[l]); printf("\n");
			
			if (indexx[i]!=i)
			{
				comtmp=B.elem[i];
				B.elem[i]=B.elem[maxk];
				B.elem[maxk]=comtmp;
			}	
	}
	delete indexx;
	return(B);
}
//------------------------------
// some math & matrix functions
//------------------------------
double sum(const matrix& a)
{
	ASSERT( a.numbands ==1 );
	
	register double *t1;
	double sum = 0.0;
	int i;
	for(i = 0, t1 = a.elem; i < a.ysize * a.xsize; i++) 
		sum += *t1++;
	return(sum);
}
//----------------------------------------------------------------------------------

matrix 	exp(const matrix& a)
{
	ASSERT( a.numbands == 1 );
	
	register double *t1, *t2;
	matrix res(a.ysize,a.xsize);
	t1=a.elem;
	t2=res.elem;
	for(int q=0;q<a.xsize*a.ysize;q++) *t2++=exp(*t1++);
	return(res);
}
//----------------------------------------------------------------------------------

matrix 	trps(const matrix &a)
{
	ASSERT( a.numbands == 1 );
	
	register double *t1,*t2;
	int i,k;
	
	matrix temp(a.xsize,a.ysize);
	for(i=0;i<a.xsize;i++)
		for(k=0,t1=temp.rowptr[i],t2=a.elem+i;k<a.ysize;k++,t2+=a.xsize)
			*t1++ =*t2;
		return(temp);
}
//----------------------------------------------------------------------------------
matrix power(const matrix &a, double p)
{
	
	register double *t1, *t2;
	matrix res(a.ysize, a.xsize, a.numbands);
	
	t1 = a.elem;
	t2 = res.elem;
	for(int q = 0; q < a.xsize * a.ysize * a.numbands; q++, t2++, t1++) 
		*t2 = pow(*t1, p);
	return(res);
}
//----------------------------------------------------------------------------------
matrix abs(const matrix& a)
{
	register double *t1, *t2;
	matrix res(a.ysize, a.xsize, a.numbands);
	
	t1 = a.elem;
	t2 = res.elem;
	for(int q=0; q < a.xsize * a.ysize * a.numbands; q++, t2++, t1++) 
		*t2 = fabs(*t1);
	return(res);
}

//----------------------------------------------------------------------------------
double 	norm2(const matrix& a)
{
	ASSERT( a.numbands == 1 );
	
	ASSERT( (a.xsize == 1) ||  (a.ysize == 1) );	
	int i, s = (a.ysize > a.xsize) ? a.ysize : a.xsize;
	register double *t;
	double sum=0.0;
	for(t = a.elem, i=0; i < s; i++) 
		sum += fabs(*t) * fabs(*t++);
	return sqrt(sum);
}

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


matrix 	inv(const matrix& A)
{//Faddeev-Leverier method
	
	ASSERT( A.numbands == 1 );
	
	ASSERT(A.xsize == A.ysize); 
	
	int N = A.xsize;
	matrix B(N,N), invA(N,N);
	
	if (N==1) 
	{
		invA.rowptr[0][0]=1.0/A.rowptr[0][0];
		return(invA);
	}
	
	B=A;
	double p=0;      //trace
	for (int i=0; i<N; i++) 
		p+=B[i][i];   
	for (int k=N-2; k>=0; k--)   //loop
		  {
		for (int l=0; l<N; l++) 
			B[l][l] -= p;
		if (k==0) 
			invA=B;
		B=A*B;
		p=0;
		for (int i=0; i<N; i++) 
			p+=B[i][i];   //trace
		p /= N-k;
		  }
	return(invA / p);
}
//----------------------------------------------------------------------------------

matrix	cov(const matrix& x)
// x is m-by-n, each row is an observation, each column is a variable
{
	ASSERT( x.numbands == 1 );
	
	register double *t1;
	int i,k, m=x.ysize, n=x.xsize;
	matrix y;
	double mean;
	//if (m==1) {printf("not an appropriate matrix (cov)\n");exit(1);}
	
	ASSERT( m > 1);
	
	y=x;
	for(k=0;k<n;k++)
		  {
		mean=0.0;
		for(i=0, t1=y.elem+k; i<m; i++, t1+=y.xsize) 
			mean+=*t1;
		mean /=m;
		for(i=0, t1=y.elem+k; i<m; i++, t1+=y.xsize) 
			*t1 -=mean;
		  }
	return((trps(y)*y)/(m-1));
}
//-------------------------------------------------
double FindMax(matrix a)
{
	ASSERT( (a.xsize > 0) && (a.ysize > 0) && (a.numbands == 1));

	double maximum = a(0,0);

	for (int j = 0; j < a.ysize; j++)
	{
		for (int i = 0; i < a.xsize; i++)
		{
			if ( a(j,i) > maximum)
				maximum = a(j,i);
		}
	}
	return maximum;
}
//-------------------------------------------------
double MaxFind(matrix &a,int bandno)
{
	ASSERT( (a.xsize > 0) && (a.ysize > 0));

	double maximum = a(0,0,bandno);

	for (int j = 0; j < a.ysize; j++)
	{
		for (int i = 0; i < a.xsize; i++)
		{
			if ( abs(a(j,i,bandno)) > maximum)
				maximum = abs(a(j,i,bandno));
		}
	}
	return maximum;
}
//-----------------------------------------------------------
double FindMin(matrix a)
{
	ASSERT( (a.xsize > 0) && (a.ysize > 0) && (a.numbands == 1));

	double minimum = a(0,0);

	for (int j = 0; j < a.ysize; j++)
	{
		for (int i = 0; i < a.xsize; i++)
		{
			if ( a(j,i) < minimum)
				minimum = a(j,i);
		}
	}
	return minimum;
}
//----------------------------------------------------------------------------------
double MinFind(matrix &a,int bandno)
{
	ASSERT( (a.xsize > 0) && (a.ysize > 0));

	double minimum = a(0,0,bandno);

	for (int j = 0; j < a.ysize; j++)
	{
		for (int i = 0; i < a.xsize; i++)
		{
			if ( a(j,i,bandno) < minimum)
				minimum = a(j,i,bandno);
		}
	}
	return minimum;
}
//-------------------------------------------------------------------
matrix DownsampleRow(matrix& mat, int factor)
// x x  ---> x x   if factor = 2
// y y       z z 
// z z          
// u u          
{
	ASSERT( factor > 0);
	
	int i, j;
	int newysize;
	
	
	newysize = (int) mat.ysize/factor;
	
	matrix a(newysize, mat.xsize, mat.numbands);
	
	double* rowPtrNew;
	double* rowPtrMat;

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

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

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

	return a;
}
//----------------------------------------------------------------------------------
matrix DownsampleColumn(matrix& mat, int factor)
// a b c d  --->  a c  if factor = 2
// w x y z        w y
{
	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++)
	{