matrixcompute.cpp

使用BDA实现三维模型相关反馈的算法 注：1、使用本系统必须在本地计算机上安装matlab 2、算法使用的三维模型特征向量是从PSB模型库中自动提取的DESIRE三维模型特征向量 3、本系统自带

#include "StdAfx.h"
#include ".\matrixcompute.h"

extern int positiveobjects[OBJECTNUMBER];
extern int negativeobjects[OBJECTNUMBER];
extern int pandnobjects[OBJECTNUMBER];
extern int pandnnum;
extern int positivenum;
extern int negativenum;
extern double originfeaturesmatrix[OBJECTNUMBER][FEATHERLENTH];
extern double mx[FEATHERLENTH];
extern double tempfeature[FEATHERLENTH];
extern double tempmatrix[FEATHERLENTH][FEATHERLENTH];
extern double Sb[FEATHERLENTH][FEATHERLENTH];
extern double Sw[FEATHERLENTH][FEATHERLENTH];
extern double tempsb[FEATHERLENTH*FEATHERLENTH];
extern double tempsw[FEATHERLENTH*FEATHERLENTH];
//extern double temptrans[MAXCLASSNUM*FEATHERLENTH];
//extern double trans[MAXCLASSNUM][FEATHERLENTH];
extern double ** trans;
extern double * transcut;


CMatrixCompute::CMatrixCompute(void)
{
	ep = engOpen(NULL);
	arraysb = mxCreateDoubleMatrix(FEATHERLENTH,FEATHERLENTH,mxREAL);
	mxSetClassName(arraysb,"sb");
	arraysw = mxCreateDoubleMatrix(FEATHERLENTH,FEATHERLENTH,mxREAL);
	mxSetClassName(arraysw,"sw");

	temptrans = new double[MAXCLASSNUM*FEATHERLENTH];
}

CMatrixCompute::~CMatrixCompute(void)
{
}

void CMatrixCompute::computemx(void)
{
	int i,j;
	double sumtemp;
	for (j=0;j<FEATHERLENTH;j++)
	{
		sumtemp = 0;
		for (i=0;i<positivenum;i++)
		{
			sumtemp += originfeaturesmatrix[positiveobjects[i]][j];

		}
		sumtemp = sumtemp/positivenum;
		mx[j] = sumtemp;

	}

}

void CMatrixCompute::computesb(void)
{
	int i,j;

	for (i=0;i<FEATHERLENTH;i++)
	{
		for (j=0;j<FEATHERLENTH;j++)
		{
			Sb[i][j] = 0;
		}
	}

	for (i=0;i<negativenum;i++)
	{
		for (j=0;j<FEATHERLENTH;j++)
		{
			tempfeature[j] = originfeaturesmatrix[negativeobjects[i]][j] - mx[j];
		}
		featuremultiself(tempfeature,tempmatrix);
		summatrix(Sb,tempmatrix);
	}

}

void CMatrixCompute::featuremultiself(double feature[FEATHERLENTH], double matrix[FEATHERLENTH][FEATHERLENTH])
{
	int i,j;
	for (i=0;i<FEATHERLENTH;i++)
	{
		for (j=0;j<FEATHERLENTH;j++)
		{
			matrix[i][j] = feature[i]*feature[j];
		}
	}
}

void CMatrixCompute::summatrix(double desmatrix[FEATHERLENTH][FEATHERLENTH], double tempmatrix[FEATHERLENTH][FEATHERLENTH])
{
	int i,j;
	for (i=0;i<FEATHERLENTH;i++)
	{
		for (j=0;j<FEATHERLENTH;j++)
		{
			desmatrix[i][j] += tempmatrix[i][j];
		}
	}
}

void CMatrixCompute::computesw(void)
{
	int i,j;

    for (i=0;i<FEATHERLENTH;i++)
    {
		for (j=0;j<FEATHERLENTH;j++)
		{
			Sw[i][j] = 0;
			
		}
    }

	for (i=0;i<positivenum;i++)
	{
		for (j=0;j<FEATHERLENTH;j++)
		{
			tempfeature[j] = originfeaturesmatrix[positiveobjects[i]][j] - mx[j];
		}
		featuremultiself(tempfeature,tempmatrix);
		summatrix(Sw,tempmatrix);
	}

}

void CMatrixCompute::computetransformmatrix(void)
{
	int i,j;
	//转换sb和sw的输入矩阵
	for (i=0;i<FEATHERLENTH;i++)//列标号
	{
		for (j=0;j<FEATHERLENTH;j++)//行标号
		{
			tempsb[i*FEATHERLENTH + j]  = Sb[j][i];
		}
	}
	for (i=0;i<FEATHERLENTH;i++)
	{
		for (j=0;j<FEATHERLENTH;j++)
		{
			tempsw[i*FEATHERLENTH + j] = Sw[j][i];
		}
	}

	memcpy((char*)mxGetPr(arraysb),(char*)tempsb,FEATHERLENTH*FEATHERLENTH*sizeof(double));
	memcpy((char*)mxGetPr(arraysw),(char*)tempsw,FEATHERLENTH*FEATHERLENTH*sizeof(double));
	engPutVariable(ep,"sb",arraysb);
	engPutVariable(ep,"sw",arraysw);

	//计算特征向量
	engEvalString(ep,"[v,d]=eig(sb,sw);");


	//筛选特征向量
	//engEvalString(ep,"[l c]=size(d);for(i=1:90);temp=-inf;pos=1;for(j=1:c);if(d(j,j)~=nan);if(d(j,j)>=temp);temp=d(j,j);pos=j;end;end;end;if(d(pos,pos)<=0);trans(1:c,i)=0;elseif(d(pos,pos)>1000000);trans(1:c,i)=1000*v(1:c,pos);else;trans(1:c,i)=(d(pos,pos)^0.5)*v(1:c,pos);end;d(pos,pos)=nan;end;");
	engEvalString(ep,"eigensum=0;[l c]=size(d);for(i=1:c);temp=-inf;pos=1;for(j=1:c);if(d(j,j)~=nan);if(d(j,j)>=temp);temp=d(j,j);pos=j;end;end;end;if(d(pos,pos)<=0);trans(1:c,i)=0;eigenv(i,i)=d(pos,pos);elseif(d(pos,pos)>1000000);trans(1:c,i)=1000*v(1:c,pos);eigensum=eigensum+1000;eigenv(i,i)=1000;else;trans(1:c,i)=(d(pos,pos)^0.5)*v(1:c,pos);eigensum=eigensum+d(pos,pos);eigenv(i,i)=d(pos,pos);end;d(pos,pos)=nan;end;");


	engEvalString(ep,"eigensum=0.9*eigensum;sumtemp=0;for(eigencut=1:c);sumtemp=sumtemp+eigenv(eigencut,eigencut);if(sumtemp>=eigensum);break;end;end;");
	////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////
	////POWER
	//engEvalString(ep,"[v,d]=eig(sw);");

	//engEvalString(ep,"[l c]=size(d);for(i=1:90);temp=-inf;pos=1;for(j=1:c);if(d(j,j)~=nan);if(d(j,j)>=temp);temp=d(j,j);pos=j;end;end;end;eigenv(i,i)=d(pos,pos);if(d(pos,pos)<=0);trans(1:c,i)=0;elseif(d(pos,pos)>1000000);trans(1:c,i)=1000*v(1:c,pos);else;trans(1:c,i)=(d(pos,pos)^0.5)*v(1:c,pos);end;d(pos,pos)=nan;end;");

	////比较eigenvalue的power的代码
	////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////
	//engEvalString(ep,"for(i=1:90);vv(1,i)=trans(1:c,i)'*sb*trans(1:c,i)/eigenv(i,i);end;");
	////如果真的根据power来选值，会怎么样？效果大幅度提高？
	////////////////////////////////////////////////////////////////////////////
	////////////////////////////////////////////////////////////////////////////

	//矩阵转置
	engEvalString(ep,"trans = trans';");//有问题

	arraytrans = engGetVariable(ep,"trans");

	arraytranscut = engGetVariable(ep,"eigencut");
	//double * temppoint = &temptrans[0];
	temptrans = mxGetPr(arraytrans);
	transcut = mxGetPr(arraytranscut);


	trans = new double*[(int)(*transcut)];
	for (i=0;i<(int)(*transcut);i++)
	{
		trans[i]= new double[FEATHERLENTH];
	}

	for (i=0;i<FEATHERLENTH;i++)
	{
		for (j=0;j<(int)(*transcut);j++)
		{
			trans[j][i] = temptrans[i*FEATHERLENTH + j];
		}
	}

	temptrans = NULL;
	mxDestroyArray(arraytrans);
	//using namespace std;
	//ofstream tempfile("test.txt");
	//for (i=0;i<MAXCLASSNUM;i++)
	//{
	//	for (j=0;j<FEATHERLENTH;j++)
	//	{
	//		tempfile<<trans[i][j]<<" ";

	//	}
	//	tempfile<<endl;
	//}

	//int tempi=0;
}