C++版本的KL变换

对KL的理解

关键词大概可以总结为去相关。
怎么讲呢？在阅读分析了别人的Matlab代码之后，发现有协方差和去均值操作，之后再去做类似卷积。
题主习惯使用帧间差异，看到去均值立马反应去相关。

实现流程

暂时做两张图的KL变换，多张图也就for循环的max大一点，后边卷积核多一点。
大致流程为：
1、每张图像按照一定规则reshape，并看作一维数据，此后，图像暂时不需要行列的概念；
2、求取每张图的灰度均值；
3、生成协方差矩阵，这个有公式，两张图就是2乘2格式；
4、协方差矩阵求特征矩阵和特征向量；
5、两张图分别去均值处理，其实就是类似帧间差异；
6、差异后的二维矩阵某一维置零；
7、置零后二维矩阵和特征矩阵卷积；
8、提取特征较大的一维作为差异，reshape后得到二维图像。

部分代码

分享几段实用代码：

1、求实对称矩阵的特征值及特征向量的雅克比法（转自其他地方）

/**
* @brief 求实对称矩阵的特征值及特征向量的雅克比法
* 利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量
* @param pMatrix				长度为n*n的数组，存放实对称矩阵
* @param nDim					矩阵的阶数
* @param pdblVects				长度为n*n的数组，返回特征向量(按列存储)
* @param pdbEigenValues			特征值数组
* @param dbEps					精度要求
* @param nJt					整型变量，控制最大迭代次数
* @return
*/	
bool CKL::JacbiCor(double * pMatrix, int nDim, double *pdblVects, double *pdbEigenValues, double dbEps, int nJt)

{
	for (int i = 0; i < nDim; i++)
	{
		pdblVects[i*nDim + i] = 1.0f;
		for (int j = 0; j < nDim; j++)
		{
			if (i != j)
				pdblVects[i*nDim + j] = 0.0f;
		}
	}
	int nCount = 0;		//迭代次数
	while (1)
	{
		//在pMatrix的非对角线上找到最大元素
		double dbMax = pMatrix[1];
		int nRow = 0;
		int nCol = 1;
		for (int i = 0; i < nDim; i++)			//行
		{
			for (int j = 0; j < nDim; j++)		//列
			{
				double d = fabs(pMatrix[i*nDim + j]);

				if ((i != j) && (d> dbMax))
				{
					dbMax = d;

					nRow = i;

					nCol = j;

				}
			}
		}

		if (dbMax < dbEps)     //精度符合要求 
			break;

		if (nCount > nJt)       //迭代次数超过限制

			break;

		nCount++;

		double dbApp = pMatrix[nRow*nDim + nRow];
		double dbApq = pMatrix[nRow*nDim + nCol];
		double dbAqq = pMatrix[nCol*nDim + nCol];

		//计算旋转角度

		double dbAngle = 0.5*atan2(-2 * dbApq, dbAqq - dbApp);

		double dbSinTheta = sin(dbAngle);

		double dbCosTheta = cos(dbAngle);

		double dbSin2Theta = sin(2 * dbAngle);

		double dbCos2Theta = cos(2 * dbAngle);



		pMatrix[nRow*nDim + nRow] = dbApp*dbCosTheta*dbCosTheta +

			dbAqq*dbSinTheta*dbSinTheta + 2 * dbApq*dbCosTheta*dbSinTheta;

		pMatrix[nCol*nDim + nCol] = dbApp*dbSinTheta*dbSinTheta +

			dbAqq*dbCosTheta*dbCosTheta - 2 * dbApq*dbCosTheta*dbSinTheta;

		pMatrix[nRow*nDim + nCol] = 0.5*(dbAqq - dbApp)*dbSin2Theta + dbApq*dbCos2Theta;
		pMatrix[nCol*nDim + nRow] = pMatrix[nRow*nDim + nCol];


		for (int i = 0; i < nDim; i++)
		{
			if ((i != nCol) && (i != nRow))
			{
				int u = i*nDim + nRow;	//p  
				int w = i*nDim + nCol;	//q
				dbMax = pMatrix[u];
				pMatrix[u] = pMatrix[w] * dbSinTheta + dbMax*dbCosTheta;
				pMatrix[w] = pMatrix[w] * dbCosTheta - dbMax*dbSinTheta;
			}
		}
		for (int j = 0; j < nDim; j++)
		{

			if ((j != nCol) && (j != nRow))
			{
				int u = nRow*nDim + j;	//p
				int w = nCol*nDim + j;	//q
				dbMax = pMatrix[u];
				pMatrix[u] = pMatrix[w] * dbSinTheta + dbMax*dbCosTheta;
				pMatrix[w] = pMatrix[w] * dbCosTheta - dbMax*dbSinTheta;
			}
		}

		//计算特征向量

		for (int i = 0; i < nDim; i++)
		{
			int u = i*nDim + nRow;		//p   
			int w = i*nDim + nCol;		//q
			dbMax = pdblVects[u];
			pdblVects[u] = pdblVects[w] * dbSinTheta + dbMax*dbCosTheta;
			pdblVects[w] = pdblVects[w] * dbCosTheta - dbMax*dbSinTheta;
		}
	}

	//对特征值进行排序以及重新排列特征向量,特征值即pMatrix主对角线上的元素
	std::map<double, int> mapEigen;
	for (int i = 0; i < nDim; i++)
	{
		pdbEigenValues[i] = pMatrix[i*nDim + i];
		mapEigen.insert(make_pair(pdbEigenValues[i], i));
	}

	double *pdbTmpVec = new double[nDim*nDim];
	std::map<double, int>::reverse_iterator iter = mapEigen.rbegin();

	for (int j = 0; iter != mapEigen.rend(), j < nDim; ++iter, ++j)
	{
		for (int i = 0; i < nDim; i++)
		{
			pdbTmpVec[i*nDim + j] = pdblVects[i*nDim + iter->second];
		}
		//特征值重新排列
		pdbEigenValues[j] = iter->first;
	}

	//设定正负号
	for (int i = 0; i < nDim; i++)
	{
		double dSumVec = 0;
		for (int j = 0; j < nDim; j++)
			dSumVec += pdbTmpVec[j * nDim + i];
		if (dSumVec < 0)
		{
			for (int j = 0; j < nDim; j++)
				pdbTmpVec[j * nDim + i] *= -1;
		}
	}
	memcpy(pdblVects, pdbTmpVec, sizeof(double)*nDim*nDim);
	delete[]pdbTmpVec;
	return 1;
}

这段代码输出后，特征值已经按照特征向量大小顺序排序。

2、矩阵转置

/**
* @brief 矩阵转置
* @param m[2][2]			长度为2*2的方阵
* @return
*/
void CKL::Matrix_Transpose(double m[2][2])
{
	int i, j;
	for (i = 1; i < 2; i++)
	{
		for (j = 0; j < i; j++)
		{
			swap(m[i][j], m[j][i]);
		}
	}
}

3、协方差矩阵

	//协方差矩阵：截取部分
	for (int i = 0; i < 2; i++)
	{
		for (int j = 0; j < 2; j++)
		{

			double tempSum = 0.0;

			int length = pixels[i].size();
			for (int k = 0; k < length; k++)
			{
				tempSum = tempSum + pixels[i][k] * pixels[j][k];
			}
			cov[i][j] = tempSum / (row*col) - meanV[i] * meanV[j];
			tempSum = 0.0;
		}
	}

最近在规范化代码，于是注释部分看起来很正式。