主成分分析（C++ & Eigen/OpenCV库）

一、原理

KL（Karhunen-Loeve）变换，也称为Hotelling变换、特征向量变换(Eigenvector-Based Transform)、主成分分析（PCA Principle Component Analysis）等。

KL变换方法是应用最广泛的一种特征提取方法之一，它是一种利用图像的统计性质的变换。在信号处理、模式识别、数字图像处理等领域已经得到了广泛的应用。其基本思想是提取出空间原始数据中的主要特征，减少数据冗余，使得数据在一个低维的特征空间被处理，同时保持原始数据的绝大部分的信息，从而解决数据空间维数过高的瓶颈问题。常用在数据压缩、特征提取等方面。

图中水平轴和垂直轴表示数据集的自然坐标轴。标号为1和2的旋转坐标轴是应用这个数据集的主变量分析产生的结果。从图可以看出数据集投影到1号轴上抓住了数据的主要特征，即具有双峰(即在它的结构上有两个聚类)的特点。

二、实现流程

输入：多波段的遥感影像，输出：变换后的各成分

三、代码实现

3.1 方法一：基于Eigen库与GDAL库

//eigen实现主成分分析
void featurenormalize(MatrixXd &X)
{
	//计算每一维度均值
	MatrixXd meanval = X.colwise().mean();
	RowVectorXd meanvecRow = meanval;
	//样本均值化为0
	X.rowwise() -= meanvecRow;
}
void computeCov(MatrixXd &X, MatrixXd &C)
{
	//计算协方差矩阵C = XTX / n-1;
	C = X.adjoint() * X;
	C = C.array() / (X.rows() - 1);
}
void computeEig(MatrixXd &C, MatrixXd &vec, MatrixXd &val)
{
	//计算特征值和特征向量，使用selfadjont按照对阵矩阵的算法去计算，可以让产生的vec和val按照有序排列（默认从大到小）
	SelfAdjointEigenSolver<MatrixXd> eig(C);
	vec = eig.eigenvectors();
	val = eig.eigenvalues();
}
int computeDim(MatrixXd &val)
{
	//输出信息量达到95%的前n主成分
	/*int dim;
	double sum = 0;
	for (int i = val.rows() - 1; i >= 0; --i)
	{
	sum += val(i, 0);
	dim = i;

	if (sum / val.sum() >= 0.95)
	break;
	}
	return val.rows() - dim;*/
	return 7;//这里设置输出7个主成分
}

void writePcaImg(const char* path, int width, int height, double *pBuff, double *adfGeo, const char *prj, int bandNum, int imageSize, int pcaInd)
{
	GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); //图像驱动
	char** ppszOptions = NULL;
	int depth = 8;//图像位深
	int dim = 1;//每个图像波段数，这里将每个主成分存储到一个单波段图像
	GDALDataset* dst = pDriver->Create(path, width, height, dim, GDT_Float64, ppszOptions);//创建图像

	if (dst == nullptr)
		printf("Can't Write Image!");
	dst->SetGeoTransform(adfGeo);//设置坐标
	dst->SetProjection(prj);//设置投影
	dst->RasterIO(GF_Write, 0, 0, width, height, &pBuff[(bandNum - pcaInd)*imageSize], width, height,
		GDT_Float64, dim, nullptr, dim*depth, width * dim * depth, depth);//写入图像
	GDALClose(dst);
}
int main(int argc, char *argv[])
{	//读取影像
	char* pszFilename = "G:/before.img";
	char *outPath = "G:/result/Eigen_PC";
	GDALDataset  *poDataset;
	GDALAllRegister();
	poDataset = (GDALDataset *)GDALOpen(pszFilename, GA_ReadOnly);
	if (poDataset == NULL)
	{
		printf_s("read failed!\n");
	}
	else
	{
		printf_s("read successful!\n");
	}
	double adfGeoTransform[6];
	if (poDataset->GetGeoTransform(adfGeoTransform) == CE_Failure)//读取坐标信息
	{
		printf("获取参数失败");
	}
	const char *prj = poDataset->GetProjectionRef();//读取投影信息
	int iWidth = poDataset->GetRasterXSize();//图像宽度
	int iHeight = poDataset->GetRasterYSize();//图像高度
	int iBandCount = poDataset->GetRasterCount();//波段数
	int iImageSize = iWidth * iHeight;//图像像元数

	double *pBuff1 = new double[iImageSize*iBandCount];//开辟空间存储原始图像
	poDataset->RasterIO(GF_Read, 0, 0, iWidth, iHeight, pBuff1,
		iWidth, iHeight, GDT_Float64, iBandCount, 0, 0, 0, 0);//读取原始图像

	MatrixXd staMat = Map<MatrixXd>(pBuff1, iImageSize, iBandCount);//将图像读入eigen矩阵
	MatrixXd X(iImageSize, iBandCount), C(iBandCount, iBandCount);//按波段存储至X矩阵，构建协方差矩阵C
	MatrixXd vec, val;//构建特征向量、特征值矩阵vec、val
	X = MatrixXd(staMat);
	MatrixXd Y = X;

	//零均值化
	featurenormalize(X);
	
	//计算协方差
	computeCov(X, C);
	
	//计算特征值和特征向量
	computeEig(C, vec, val);
	cout << "特征值为：" << val << endl;
	cout << "特征向量为：" << vec << endl;

	//计算损失率，确定降低维数
	int dim = computeDim(val);

	//计算结果
	MatrixXd res = Y * vec.rightCols(dim);

	//验证结果
	featurenormalize(res);
	MatrixXd CY;
	computeCov(res, CY);
	MatrixXd vecY, valY;
	computeEig(CY, vecY, valY);
	
	/*cout << "主成分分析后的特征向量："<<endl << vecY << endl;
	cout << "主成分分析后的特征值：" << endl<< valY << endl;*/

	//将主成分分量存储至pBuff2
	double *pBuff2 = new double[iImageSize*iBandCount];
	for (int i = 0; i < dim; ++i)
	{
		for (int j = 0; j < iImageSize; ++j)
		{
			pBuff2[i*iImageSize + j] = res(j, i);
		}
	}

	//各个主成分写入图像（包含坐标及投影信息）
	for (int i = 0; i < iBandCount; i++)
	{
		char x[] = " ";
		strcpy(x, outPath);
		char dstPath[10] = {};
		sprintf(dstPath, "%d.tif", i + 1);
		strcat(x, dstPath);
		writePcaImg(x, iWidth, iHeight, pBuff2, adfGeoTransform, prj, 7, iImageSize, i + 1);
		cout << "pca " << i + 1 << " complete" << endl;
	}
	cout << "pca complete!" << endl;
	cin.get();
	return 0;
}

3.2 方法二：基于openCV库与GDAL库

//图像灰度化函数
static Mat toGrayscale(InputArray _src) {
	Mat src = _src.getMat();
	// only allow one channel
	if (src.channels() != 1) {
		CV_Error(CV_StsBadArg, "Only Matrices with one channel are supported");
	}
	// create and return normalized image
	Mat dst;
	cv::normalize(_src, dst, 0, 255, NORM_MINMAX, CV_8UC1);
	return dst;
}

int Cal_PCA(Mat data, float T, String path, Mat &Eigenvector)
{
	Mat data_t = data.t();	
	PCA pca(data.t(), Mat(), CV_PCA_DATA_AS_ROW,T);		//这个pca就已经算出来了
	Mat pca_tMean = pca.mean;

	//零均值化
	Mat data_tFloat;
	data_t.convertTo(data_tFloat, CV_32FC1);	//类型转换
	for (int i = 0; i < data_t.cols; i++)
	{
		for (int j = 0; j < data_t.rows;j++)
		{
			data_tFloat.at<float>(j, i) -= pca_tMean.at<float>(0, i);
		}
	}

	//零均值化后的数组计算主成分
	PCA pca_tnew(data_tFloat, Mat(), CV_PCA_DATA_AS_ROW,T);
	Mat Eigenvalues = pca_tnew.eigenvalues.t();		//特征值
	
	Eigenvector = pca_tnew.eigenvectors;		//特征向量
	float sum_Ei;
	for (int i = 0; i < Eigenvector.rows; i++)
	{
		sum_Ei = 0;
		for (int j = 0; j < Eigenvector.cols; j++)
		{
			sum_Ei += abs(Eigenvector.at<float>(i,j));
		}

		for (int j = 0; j < Eigenvector.cols; j++)
		{
			Eigenvector.at<float>(i,j) = Eigenvector.at<float>(i,j) / sum_Ei;
		}
	}
	cout <<"new_Eigenvector="<< Eigenvector << endl;
	cout << "Eigenvalues=" << Eigenvalues << endl;
	//根据需要达到的百分比设置成分个数
	double Eigenvalues_sum = sum(Eigenvalues).val[0];	//特征值的和，val是干嘛？回答：将一组字符型数据的数字部分转换成相应的数值型数据
	int Eigenvalues_len = Eigenvalues.cols;
	float flag_sum = 0;
	int i = 0;
	while (flag_sum <= T && i < Eigenvalues_len)                    // 计算达到要求的主成分数
	{
		flag_sum = flag_sum + Eigenvalues.ptr<float>(0)[i] / Eigenvalues_sum;	//从第几个开始就能达到所需要的比重
		i++;
	}
	//主要参数存入
	FileStorage f_pca(path, FileStorage::WRITE);         
	f_pca << "num" << i;
	//f_pca << "coff" << coeff;
	f_pca << "Eigenvector" << Eigenvector;

	Eigenvalues.release();	//释放内存
	//coeff.release();
	f_pca.release();
	return 0;
}
/*	
	data：输入数据
	path：数据保存路径
	Eigenvector:特征向量
	GDALDataset：数据集
	iHeight：图像高
 */
void Output_PCA(Mat data, String path, Mat Eigenvector, GDALDataset *poDS, int iHeight, int iWidth, int nBands, int imageSize, double *adfGeo, const char *prj)
{
	GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); //图像驱动
	FileStorage f_pca(path, FileStorage::READ);
	int num;
	const char* sPath = "G:/result/a";

	//把数据读取出
	f_pca["num"] >> num;
	f_pca.release();

	cout << "Eigenvector=" << Eigenvector << endl;
	Mat data_float;
	data.convertTo(data_float, CV_32FC1);
	/*Mat Eigenvector_8U;
	Eigenvector.convertTo(Eigenvector_8U, CV_8U);*/
	Mat Y = Eigenvector * data_float;

	//进行波段循环，将矩阵形式表示的图像数据以行向量形式表示，如92*112的图像，表示为1*10304,
	for (int r = 1; r <= num; ++r)
	{
		GDALRasterBand *xBand = poDS->GetRasterBand(r);
		Mat ReCons;
		ReCons = Y.row(r - 1).reshape(data.channels(), iHeight); // 从行向量重塑为图像形状
		ReCons = toGrayscale(ReCons);
		
		switch (r)
		{
		case 1:
			imwrite("G:/result/Opencv_PC1.tif", ReCons);break;
		case 2:
			imwrite("G:/result/Opencv_PC2.tif", ReCons);break;
		case 3:
			imwrite("G:/result/Opencv_PC3.tif", ReCons);break;
		case 4:
			imwrite("G:/result/Opencv_PC4.tif", ReCons);break;
		case 5:
			imwrite("G:/result/Opencv_PC5.tif", ReCons);break;
		case 6:
			imwrite("G:/result/Opencv_PC6.tif", ReCons);break;
		case 7:
			imwrite("G:/result/Opencv_PC7.tif", ReCons);break;
		}
		cout << "PC" << r << "complete" << endl;
	}
}

void UpdateRasterFile(const char* pszRasterFile ,const char* txtPath)
{
	//使之支持中文路径
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

	//注册栅格驱动
	GDALAllRegister();

	//打开要更新的数据，第二个参数使用介个
	GDALDataset *poDS = (GDALDataset*)GDALOpen(pszRasterFile, GA_Update);
	if (poDS == NULL)
	{
		cout << "打开图像" << pszRasterFile << "失败";
		return;
	}

	double adfGeoTransform[6];
	if (poDS->GetGeoTransform(adfGeoTransform) == CE_Failure)//读取坐标信息
		{
			printf("获取参数失败");
		}
	const char *prj = poDS->GetProjectionRef();//读取投影信息

	//获取图像大小
	int iWidth = poDS->GetRasterXSize();
	int iHeight = poDS->GetRasterYSize();
	//获取波段个数
	int nBands = poDS->GetRasterCount();

	//创建数组
	Mat data(nBands, iWidth*iHeight, CV_8U, Scalar::all(0));
	//进行波段循环，将矩阵形式表示的图像数据以行向量形式表示
	for (int r = 1; r <= nBands; ++r)
	{
		GDALRasterBand *xBand = poDS->GetRasterBand(r);
		unsigned char *pBuf = new unsigned char[iWidth * iHeight]; //动态指配数组

		xBand->RasterIO(GF_Read, 0, 0, iWidth, iHeight, pBuf, iWidth, iHeight, GDT_Byte, 0, 0);
		for (int i = 0; i < data.cols; i++)
		{
			data.at<uchar>(r - 1, i) = pBuf[i];
		}
	}
	Mat Eigenvector;
	Cal_PCA(data, 1, txtPath, Eigenvector);	//第一步：进行主成分分析
	Output_PCA(data, txtPath, Eigenvector, poDS, iHeight, iWidth, nBands, iHeight*iWidth, adfGeoTransform, prj);	//第二步：输出主成分图像
	system("pause");
}
int main()
{
	const char *imagePath = "G:/before.img";
	const char *txtPath = "G:/b.txt";
	UpdateRasterFile(imagePath,txtPath);
	return 0;
}

四、处理结果

原始图像

结果图像

用Eigen和GDAL库处理结果（PC1 - PC6）：

用openCV和GDAL库处理结果（PC1 - PC6）：

五、结果分析

由生成的主成分分析后的图像可以看出第一二主成分所占的信息最多，依次递减，第七主成分所包含的信息最少，基本只剩下噪声。基本符合主成分分析的结果。
不同的库其调用方式不同，在调用openCV库时应该先将GDAL导入的图像矩阵转置，生成的才是正确的特征值，再将其特征向量单位化，才能得到与eigen库相同的结果，从这个方面来看，openCV库主成分分析函数具有更差的适用性。