matrixcompute.cpp

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

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

extern int allclassobjectmembers[MAXCLASSNUM][MAXOBJECTNUM];
extern int numofoneclass[MAXCLASSNUM];
extern double originfeaturematrix[OBJECTSNUMBER][FEATURELENGTH];
extern double kernelfeaturematrix[OBJECTSNUMBER][OBJECTSNUMBER];
extern int positiveobjects[OBJECTSNUMBER];
extern int negativeobjects[OBJECTSNUMBER];
extern int pandnobjects[OBJECTSNUMBER];
extern int positivenum;
extern int negativenum;
extern int pandnnum;
extern double ** trans;
extern double ** newfeaturematrix;


CMatrixCompute::CMatrixCompute(void)
: Kx(NULL)
, Ky(NULL)
{
	//int i,j,k;
	//double sum;
	//for (i=0;i<OBJECTSNUMBER;i++)
	//{
	//	for (j=0;j<OBJECTSNUMBER;j++)
	//	{
 //            sum = 0;
	//		 for (k=0;k<FEATURELENGTH;k++)
	//		 {
	//			 sum += originfeaturematrix[i][k]*originfeaturematrix[j][k];
	//		 }
	//		 dotfeaturematrix[i][j] = sum;
	//	}
	//}


	double sum;
	double sumtemp;

	int i,object1,object2;

	for (object1=0;object1<OBJECTSNUMBER;object1++)
	{
		for (object2=0;object2<OBJECTSNUMBER;object2++)
		{
			sum = 0;
			sumtemp = 0;
			for (i=0;i<FEATURELENGTH;i++)
			{
				sumtemp += originfeaturematrix[object1][i]*originfeaturematrix[object1][i];
			}
			sum += sumtemp;
			sumtemp = 0;
			for (i=0;i<FEATURELENGTH;i++)
			{
				sumtemp += originfeaturematrix[object2][i]*originfeaturematrix[object2][i];
			}
			sum +=sumtemp;
			sumtemp = 0;
			for (i=0;i<FEATURELENGTH;i++)
			{
				sumtemp += originfeaturematrix[object1][i]*originfeaturematrix[object2][i];
			}
			sum = sum - 2*sumtemp;
            kernelfeaturematrix[object1][object2] = exp(-GAMMA*sum);

		}
	}

	//tempsb = new double[OBJECTSNUMBER*OBJECTSNUMBER];

	//using namespace std;
	//ofstream tempfile("test.txt");
	//int j;
	//for (i=0;i<OBJECTSNUMBER;i++)
	//{
	//	for (j=0;j<OBJECTSNUMBER;j++)
	//	{
	//		tempfile<<kernelfeaturematrix[i][j]<<" ";
	//	}
	//	tempfile<<endl;
	//}
	//tempfile.close();

	//double tempi = exp(-GAMMA*sum);

	//return exp(-GAMMA*sum);
	ep = engOpen(NULL);


}

CMatrixCompute::~CMatrixCompute(void)
{
}

bool CMatrixCompute::computeKx(void)
{
	using namespace std;
	Kx = new double * [pandnnum];
	int i,j;
	for (i=0;i<pandnnum;i++)
	{
        Kx[i] = new double[positivenum];
	}

	for (i=0;i<positivenum;i++)
	{
		for (j=0;j<pandnnum;j++)
		{
			Kx[j][i] = kernelfeaturematrix[pandnobjects[j]][positiveobjects[i]];
		}
	}

 //   ofstream tempfile("test.txt");
 //   for (i=0;i<pandnnum;i++)
 //   {
	//	for (j=0;j<positivenum;j++)
	//	{
	//		tempfile<<Kx[i][j]<<" ";
	//	}
	//	tempfile<<endl;
 //   }
	//tempfile.close();

	return true;
}


double CMatrixCompute::k_function(int object1, int object2)
{
	double sum = 0;
	double sumtemp = 0;

	int i;

	for (i=0;i<FEATURELENGTH;i++)
	{
		sumtemp += allclassobjectmembers[object1][i]*allclassobjectmembers[object1][i];
	}
    sum += sumtemp;
	sumtemp = 0;
	for (i=0;i<FEATURELENGTH;i++)
	{
        sumtemp += allclassobjectmembers[object2][i]*allclassobjectmembers[object2][i];
	}
	sum +=sumtemp;
	sumtemp = 0;
	for (i=0;i<FEATURELENGTH;i++)
	{
		sumtemp += allclassobjectmembers[object1][i]*allclassobjectmembers[object2][i];
	}
    sum = sum - 2*sumtemp;
    
	//double tempi = exp(-GAMMA*sum);

	return exp(-GAMMA*sum);
}

bool CMatrixCompute::computeKy(void)
{
	using namespace std;
	Ky = new double * [pandnnum];
	int i,j;
	for (i=0;i<pandnnum;i++)
	{
		Ky[i] = new double[negativenum];
	}

	for (i=0;i<negativenum;i++)
	{
		for (j=0;j<pandnnum;j++)
		{
			Ky[j][i] = kernelfeaturematrix[pandnobjects[j]][negativeobjects[i]];
		}
	}

	//  ofstream tempfile("test.txt");
	//  for (i=0;i<pandnnum;i++)
	//  {
	//for (j=0;j<negativenum;j++)
	//{
	//	tempfile<<Ky[i][j]<<" ";
	//}
	//tempfile<<endl;
	//  }
	//  tempfile.close();


	return true;
}

bool CMatrixCompute::computeSbSw(void)
{
	//formerSb = new double*[pandnnum];
	//int i,j;
	//for (i=0;i<pandnnum;i++)
	//{
	//	formerSb[i] = new double[negativenum];
	//}
	using namespace std;
	int i,j,k;
	//计算Kx×Inx
	double** tempTNy = new double*[pandnnum];
    for (i=0;i<pandnnum;i++)
    {
		tempTNy[i]= new double[negativenum];
    }
	double tempsum;
	for (i=0;i<pandnnum;i++)
	{
        tempsum=0;
		for (k=0;k<positivenum;k++)
		{
			tempsum += Kx[i][k];
		}
		tempsum = tempsum * (1/double(positivenum));
		for (j=0;j<negativenum;j++)
		{
			tempTNy[i][j] = tempsum;
		}
	}
    //计算Ky－Kx×Inx
	for (i=0;i<pandnnum;i++)
	{
		for (j=0;j<negativenum;j++)
		{
			tempTNy[i][j] = Ky[i][j] - tempTNy[i][j];
		}
	}

	//为sb和sw分配空间
	Sb = new double*[pandnnum];
	for (i=0;i<pandnnum;i++)
	{
		Sb[i] = new double[pandnnum];
	}
	Sw =new double*[pandnnum];
	for (i=0;i<pandnnum;i++)
	{
		Sw[i] = new double[pandnnum];
	}

	for (i=0;i<pandnnum;i++)
	{
		for (j=0;j<pandnnum;j++)
		{
			tempsum = 0;
			for (k=0;k<negativenum;k++)
			{
				tempsum += tempTNy[i][k] * tempTNy[j][k];
			}
			Sb[i][j] = tempsum;
		}
	}

	//tempty没用了
	for (i=0;i<pandnnum;i++)
	{
		delete[] tempTNy[i];
	}
	delete[] tempTNy;


	double** tempNxNx = new double*[positivenum];
	for (i=0;i<positivenum;i++)
	{
		tempNxNx[i] = new double[positivenum];
	}
	double** tempTNx = new double*[pandnnum];
	for (i=0;i<pandnnum;i++)
	{
		tempTNx[i]= new double[positivenum];
	}

    for (i=0;i<positivenum;i++)
    {
		for (j=0;j<positivenum;j++)
		{
			if(i==j)
				tempNxNx[i][j] = 1- (1/double(positivenum));
			else
				tempNxNx[i][j] = -1/double(positivenum);
		}
    }

	for (i=0;i<pandnnum;i++)
	{
		for (j=0;j<positivenum;j++)
		{
			tempsum = 0;
			for (k=0;k<positivenum;k++)
			{
				tempsum += Kx[i][k]*tempNxNx[k][j];
			}
			tempTNx[i][j] = tempsum;
		}
	}
	//计算sw
	for (i=0;i<pandnnum;i++)
	{
		for (j=0;j<pandnnum;j++)
		{
			tempsum = 0;
			for (k=0;k<positivenum;k++)
			{
				tempsum += tempTNx[i][k]*Kx[j][k];
			}
			Sw[i][j] = tempsum;
		}
	}

	//清除变量
	for (i=0;i<pandnnum;i++)
	{
		delete[] tempTNx[i];
	}
	delete[] tempTNx;
	for (i=0;i<positivenum;i++)
	{
		delete[] tempNxNx[i];
	}
	delete[] tempNxNx;
	for (i=0;i<pandnnum;i++)
	{
		delete[] Kx[i];
	}
	delete[] Kx;
	for (i=0;i<pandnnum;i++)
	{
		delete[] Ky[i];
	}
	delete[] Ky;

	//ofstream tempfile("test.txt");
	//for (i=0;i<pandnnum;i++)
	//{
	//	for (j=0;j<pandnnum;j++)
	//	{
	//		tempfile<<Sw[i][j]<<" ";
	//	}
	//	tempfile<<endl;
	//}

	//double** tempiiii=NULL;

	//delete[] tempiiii;
	//compute formerSb
	return true;
}

bool CMatrixCompute::computetransformmatrix(void)
{
	int i,j;
	double * tempsb ;
	tempsb = new double[pandnnum*pandnnum];
	double * tempsw ;
	tempsw = new double[pandnnum*pandnnum];
	mxArray * arraysb;
	mxArray * arraysw;
	arraysb = mxCreateDoubleMatrix(pandnnum,pandnnum,mxREAL);
	mxSetClassName(arraysb,"sb");
	arraysw = mxCreateDoubleMatrix(pandnnum,pandnnum,mxREAL);
	mxSetClassName(arraysw,"sw");

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

	memcpy((char*)mxGetPr(arraysb),(char*)tempsb,pandnnum*pandnnum*sizeof(double));
	memcpy((char*)mxGetPr(arraysw),(char*)tempsw,pandnnum*pandnnum*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: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;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;");


	arraysb = engGetVariable(ep,"trans");
	//arraysw = engGetVariable(ep,"d");

	double* temptrans = mxGetPr(arraysb);
	//tempsw = mxGetPr(arraysw);
	//这句话影响了指针的释放，这是为什么呢？查查
	//这个不是长度的问题

    trans = new double*[pandnnum];
	for (i=0;i<pandnnum;i++)
	{
		trans[i] = new double[pandnnum];
	}
	for (i=0;i<pandnnum;i++)
	{
		for (j=0;j<pandnnum;j++)
		{
			trans[j][i] = temptrans[i*pandnnum + j];
		}
	}
	//for (i=0;i<pandnnum;i++)
	//{
	//	for (j=0;j<pandnnum;j++)
	//	{
	//		Sw[j][i] = tempsw[i*pandnnum + j];
	//	}
	//}

	//using namespace std;  
	//ofstream tempfile("test.txt");
	//for (i=0;i<pandnnum;i++)
	//{
	//	for (j=0;j<pandnnum;j++)
	//	{
	//		tempfile<<Sw[i][j]<<" ";
	//	}
	//	tempfile<<endl;
	//}


	for (i=0;i<pandnnum;i++)
	{
		delete[] Sb[i];
	}
	delete[] Sb;
	for (i=0;i<pandnnum;i++)
	{
		delete[] Sw[i];
	}
	delete[] Sw;
	delete[] tempsb;
	delete[] tempsw;
	mxDestroyArray(arraysb);
	mxDestroyArray(arraysw);

	return true;
}

bool CMatrixCompute::computenewfeaturematrix(void)
{
    newfeaturematrix = new double*[OBJECTSNUMBER];
	int i,j,k;
	double tempsum;
	for (i=0;i<OBJECTSNUMBER;i++)
	{
		newfeaturematrix[i] = new double[pandnnum];
	}

	for (i=0;i<OBJECTSNUMBER;i++)
	{
		for (j=0;j<pandnnum;j++)
		{
			tempsum = 0;
			for (k=0;k<pandnnum;k++)
			{
				tempsum += kernelfeaturematrix[pandnobjects[k]][i]*trans[j][k];
			}
			newfeaturematrix[i][j] = tempsum;

		}
	}

	for (i=0;i<pandnnum;i++)
	{
		delete[] trans[i];
	}
	delete[] trans;
	return true;
}