c++ GDAL 完成PCA

c++ GDAL 完成PCA

int main()
{
	const char* imagePath = "D:/data/before.img";
	const char* txtPath = "D:/data/result.txt";
	UpdateRasterFile(imagePath, txtPath);
	return 0;
}

通过UpdateRasterFile(imagePath, txtPath);完成PCA的计算和输出

一、获取图像信息

void PCAProcess(const char* imagePath, const char* txtPath)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	
	GDALDataset* Img = (GDALDataset*)GDALOpen(imagePath, GA_Update);
	if (Img == NULL)
	{
		cout << "打开图像" << imagePath << "失败";
		return;
	}
	else
	{
		cout << "打开图像成功！" << endl;
	}

	double GeoTransform[6];
	Img->GetGeoTransform(GeoTransform);
	const char* Project = Img->GetProjectionRef();//读取投影信息
	int nBands = Img->GetRasterCount();
	int XSize = Img->GetRasterXSize();
	int YSize = Img->GetRasterYSize();
	int XOff = 0;//偏移
	int YOff = 0;
	int BufXSize = XSize;//屏幕显示大小
	int BufYSize = YSize;

	Mat Dataset(nBands, XSize * YSize, CV_8U, Scalar::all(0));//创建一个矩阵

	for (int n = 1; n <= nBands; ++n)//逐波段
	{
		GDALRasterBand* xBand = Img->GetRasterBand(n);//获得第i个波段
		unsigned char* pBuffer = (unsigned char*)CPLMalloc(sizeof(unsigned char) * (XSize) * (YSize) * 1);
		xBand->RasterIO(GF_Read, XOff, YOff,XSize, YSize, pBuffer, BufXSize, BufYSize, GDT_Byte, 0, 0);
		for (int i = 0; i < Dataset.cols; i++)
		{
			Dataset.at<uchar>(n - 1, i) = pBuffer[i];//将多个一维数组写成二维数组
		}
	}
	Mat Eigenvector;//建立特征向量矩阵
	Calculate(nBands, Dataset, 1, txtPath, Eigenvector);//第一步：进行主成分分析
	Output(Dataset, txtPath, Eigenvector, Img, YSize, XSize, nBands, YSize * XSize, GeoTransform, Project);	//第二步：输出主成分图像

	
	system("pause");
}

思路：

1. 支持中文路径
2. 注册栅格驱动
3. 获取打开图像的指针
4. 读取相应信息：坐标、投影、图像长宽、波段个数
5. 创建一个二维数组
6. 逐波段，将图像信息线转变成一维数组，再写入二维数组中
7. 建立特征向量的矩阵
8. 计算主成分（另一个函数）
9. 输出计算结果（另一个函数）

Notes:

1. GDAL中有一个函数CPLGetConfigOption( “GDAL_FILENAME_IS_UTF8”, “YES” )判断，通过判断是否是UTF8的编码，而且指定的默认值还是UTF8编码，在含有中文路径的字符串大多数的编码应该是GBK的编码，这样，系统就将GBK的编码当做UTF8的编码来进行转换，结果就是汉字全部是乱码，导致的结果就是找不到文件，所以打不开。
   解决方法：
   不改变GDAL源代码，在自己调用GDALRegisterAll()和OGRAllRegiser()函数后，加上下面一句即可。
   CPLSetConfigOption(“GDAL_FILENAME_IS_UTF8”,“NO”);
2. 注册驱动？？？GDAL所有操作都需要先注册格式
3. 坐标转换：
   影像左上角横坐标：geoTransform[0]
   影像左上角纵坐标：geoTransform[3]
   遥感图像的水平空间分辨率为geoTransform[1]
   遥感图像的垂直空间分辨率为geoTransform[5]
   通常geoTransform[5] 与 geoTransform[1]相等
   如果遥感影像方向没有发生旋转，即上北、下南，则
   geoTransform[2] 与 row *geoTransform[4] 为零。

原文链接：https://blog.csdn.net/schmiloo/article/details/107042643
链接: link.
Xp = padfTransform[0] + PpadfTransform[1] + LpadfTransform[2];
Yp = padfTransform[3] + PpadfTransform[4] + LpadfTransform[5];

原文链接：http://blog.sina.com.cn/s/blog_65b644f90100jmum.html
链接: link.

4. 获得投影信息：
   函数 GetProjectionRef（）通过这个可以获得仿射坐标系统的描述，获得的是一个指向内部投影参考字符串的指针。

描述信息： const char* GDALDataset.GetDriver().GetDescription()，通常是图像的格式
图像大小： 图像宽度 int GDALDataset.GetRasterXSize()
图像高度 int GDALDataset.GetRasterYSize()
波段数： int GDALDataset.GetRasterCount()
波段少则一个，多则很多个，在遥感影象中波段通常有多个。
投影信息： GDALDataset.GetProjectionRef()
有的图像没有投影信息，如一般的JPG、BMP格式图像。
地理坐标信息：double adfGeoTransform[6]
GDALDataset.GetGeoTransform(adfGeoTransform)
地理坐标信息是一个含6个double型数据的数组，adfGeoTransform[1]和adfGeoTransform[5]表示东西和南北方向 一个像素对应的距离，adfGeoTransform[0]和adfGeoTransform[3]表示左上角的坐标。
波段信息： 数据集中重要的信息，有波段尺寸、数据类型、颜色信息等。
获取波段的方法：
Band band= Dataset.GetRasterBand(i)
波段尺寸： int band.GetXSize()
int band.GetYSize()
数据类型： const char* GDALGetDataTypeName(poBand.GetRasterDataType())
颜色信息： const char* GDALGetColorInterpretationName(poBand.GetColorInterpretation())

原文链接：https://blog.csdn.net/weixin_41340063/article/details/82733527

5. 一般使用的灰度图像元素是8位的uchar型（注意不是无符号整型，而是无符号字符型。元素值0~255），自定义Mat矩阵时，推荐使用8位无符号字符型，即CV_8U、CV_8UC3。
   再说一遍，条件反射式的记住：一般见到的Mat的元素是8位无符号字符型，不是8位无符号整型。
   链接: link.
6. Scalar::all(0)就是给每个通道都赋值0
7. xBand->RasterIO(GF_Read, 0, 0, iWidth, iHeight, pBuf, iWidth, iHeight, GDT_Byte, 0, 0);
   RasterOI函数：
   CPLErr GDALDataset::RasterIO(
   GDALRWFlag eRWFlag, //读写标记
   int nXOff, int nYOff,//相对于图像左上角顶点(从零开始)的行列偏移量
   int nXSize, int nYSize,//要读写的块在x方向的象素个数和y方向的象素列数
   void * pData,//指向目标缓冲区的指针，由用户分配
   int nBufXSize, int nBufYSize,//目标块在x方向上和y方向上的大小
   GDALDataType eBufType, //目标缓冲区的数据类型，原类型会自动转换为目标类型
   int nBandCount, //要处理的波段数
   int * panBandMap,//记录要操作的波段的索引(波段索引从1开始)的数组，若为空，//则数组中存放的是前nBandCount个波段的索引
   int nPixelSpace,//x方向上两个相邻象素之间的字节偏移，默认为0,则列间的实际//字节偏移由目标数据类型eBufType确定
   int nLineSpace, //y方向上相邻两行之间的字节偏移, 默认为0,则行间的实际字节//偏移为eBufType * nBufXSize
   int nBandSpace //相邻两波段之间的字节偏移，默认为0,则意味着波段是顺序结构//的，其间字节偏移为nLineSpace * nBufYSize
   );
   具体参照，写得非常详细！！！
   链接: link.
   8.data.at<uchar>(r - 1, i) = pBuf[i];
   对于单通道，可以直接使用mat.at(row,col)即可获取像素值。

二、计算PCA

int Calculate(int bandnum,Mat data, float percent, String path, Mat& Eigenvector)
{
	
	Mat T = data.t();//矩阵转置
	PCA pca(T, Mat(), CV_PCA_DATA_AS_ROW, percent);
	Mat pca_tMean = pca.mean;//平均值
	cout << "平均值" << pca.mean << endl;

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

	//计算主成分
	PCA pca_final(data_TF, Mat(), CV_PCA_DATA_AS_ROW, percent);
	Mat Eigenvalues = pca_final.eigenvalues.t();//获取特征值并转置

	Eigenvector = pca_final.eigenvectors;
	float sum_Eivactor;
	for (int i = 0; i < Eigenvector.rows; i++)
	{
		sum_Eivactor = 0;
		for (int j = 0; j < Eigenvector.cols; j++)
		{
			sum_Eivactor += 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_Eivactor;
		}
	}
	cout << "特征向量：" << Eigenvector << endl;
	cout << "特征值：" << Eigenvalues << endl;

	//通过设置的信息量确定输出的主成分个数
	double sum_Eivalue = sum(Eigenvalues).val[0];//数字转换成相应的数值型数据并求和
	float Percent = 0;
	int i = 0;
	while (Percent <= percent && i < bandnum) 
	{
		Percent = Percent + Eigenvalues.ptr<float>(0)[i] / sum_Eivalue;	
		i++;
	}
	//存入数据
	FileStorage result(path, FileStorage::WRITE);//把计算结果写入txt中
	result << "num" << i;//输出几个主成分

	//result << "Eigenvector" << Eigenvector;//输出特征向量

	Eigenvalues.release();	//释放内存

	result.release();
	return 0;
}

思路：

1. 矩阵转置
2. 计算PCA后的平均值
3. 进行零均值化
4. 再次进行PCA
5. 获取计算结果的特征值、特征向量
6. 计算新的特征向量
7. 根据需要达到的百分比计算成分个数
8. 将计算结果存入txt文档中

Notes:

1. 参数：Mat data是建立的二维数组、T是？？？、path是文本文档的路径、Mat& Eigenvector是特征向量
2. .t()矩阵转置
3. 关于PCA函数，详见链接: link.

cv::PCA::PCA	(	InputArray 	data,
InputArray 	mean,
int 	flags,
int 	maxComponents = 0 
)	

parameters:
data 输入存储为矩阵行或矩阵列的样本。
mean 可选平均值；如果矩阵为空（noArray（）），则根据数据计算平均值。
flags 操作标志；当前该参数仅用于指定数据布局（PCA:：Flags）
maxComponents PCA 应保留的最大组件数量；默认情况下，将保留所有零部件。

retainedVariance 主成分分析应保留的方差百分比。使用此参数将使PCA决定保留多少组件，但它将始终至少保留2个组件。???

1. convertTo函数
   链接: link.

void cv::cuda::GpuMat::convertTo(OutputArray dst,int rtype 
)const

将Mat转换为其他数据类型（阻塞调用）
CV_32FC1，那么后面对该矩阵的输入输出的数据指针类型都应该是float，这在32位编译器上是32位浮点数，也就是单精度。
5. 零均值化
将每个像素的值减去训练集上所有像素值的平均值，比如已计算得所有像素点的平均值为128，所以减去128后，现在的像素值域即为[-128,127]，即满足均值为零。

data_tFloat.at<float>(j, i) -= pca_tMean.at<float>(0, i);

每个像素上的值减去该列的平均值

eigenvalues
Mat cv::PCA::eigenvalues
eigenvalues of the covariation matrix//协方差矩阵的特征值

eigenvectors
Mat cv::PCA::eigenvectors
eigenvectors of the covariation matrix//协方差矩阵的特征向量

7. double Eigenvalues_sum = sum(Eigenvalues).val[0];
   求特征值的和
   .val[0]是用来获得向量中的数字，将一组字符型数据的数字部分转换成相应的数值型数据
8. Mat像素访问指针ptr<>与at<>的用法区别
   具体的区别和用法参见链接: link.
   最快区别为：
   使用ptr指针访问像素时访问任意一行像素的首地址，特别方便图像的一行一行的横向访问，如果需要一列一列的纵向访问图像，就稍微麻烦一点，<>中一般是uchar。
   使用at访问像素时，使用调用方法image.atcv::Vec3b(j,i)，at方法方便，直接给i，j赋值就可以随意访问图像中任何一个像素，其中j表示第j行，i表示该行第i个像素。<>中一般是Vec3b
   原文链接：link.

三、输出图像

void Output(Mat data, String path, Mat Eigenvector, GDALDataset* Img, int iHeight, int iWidth, int nBands, int imageSize, double* adfGeo, const char* prj)
{
	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); //图像驱动
	FileStorage Save(path, FileStorage::READ);//存储读取的信息
	int num;
	

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

	Mat data_float;
	data.convertTo(data_float, CV_32FC1);

	Mat Y = Eigenvector * data_float;//转回去

	for (int band = 1; band <= num; ++band)
	{
		GDALRasterBand* nBand = Img->GetRasterBand(band);
		Mat Show = Y.row(band - 1).reshape(data.channels(), iHeight); // 从行向量重塑为图像形状

		//灰度化
		InputArray _src = Show;
		Mat src = _src.getMat();
		if (src.channels() != 1)
		{
			CV_Error(CV_StsBadArg, "只支持单通道");//处理异常的机制
		}
		cv::normalize(_src, Show, 0, 255, NORM_MINMAX, CV_8UC1);


		if (band == 1)
		{
			imwrite("D: / data2 / PC1.tif", Show);
			cout<<"主成分1完成"<<endl;
			/*Mat img1 = imread("D: / data2 / PC1.tif");
			imshow("PC1", img1);
			waitKey(0);*/
			
		}
		else if (band == 2)
		{
			imwrite("D: / data2 / PC2.tif", Show);
			cout << "主成分2完成" << endl;
		}
		else if (band == 3)
		{
			imwrite("D: / data2 / PC3.tif", Show);
			cout << "主成分3完成" << endl;
		}
		else if (band == 4)
		{
			imwrite("D: / data2 / PC4.tif", Show);
			cout << "主成分4完成" << endl;
		}
		else if (band == 5)
		{
			imwrite("D: / data2 / PC5.tif", Show);
			cout << "主成分5完成" << endl;
		}
		else if (band == 6)
		{
			imwrite("D: / data2 / PC6.tif", Show);
			cout << "主成分6完成" << endl;
		}
		else
		{
			imwrite("D: / data2 / PC7.tif", Show);
			cout << "主成分7完成" << endl;
		}

		
	}
}

思路：

• Output_PCA目的：用于将二维数组和txt文档中的内容转换成图像输出

1. 注册图像注册
2. 存储要读取的信息
3. 读取文档中的信息
4. 建立特征向量矩阵
5. 将图像以行表示出来
6. 行向量重组成图像

Notes
1.FileStorage类
FileStorage类将各种OpenCV数据结构的数据存储为XML 或 YAML格式。同时，也可以将其他类型的数值数据存储为这两种格式。
构造函数
FileStorage类的构造函数为：
cv::FileStorage(const string& source, int flags， const string& encoding=string())

source –存储或读取数据的文件名（字符串），其扩展名(.xml 或 .yml/.yaml)决定文件格式。
flags – 操作模式，包括：
FileStorage::READ 打开文件进行读操作
FileStorage::WRITE 打开文件进行写操作
FileStorage::APPEND打开文件进行附加操作
FileStorage::MEMORY 从source读数据，或向内部缓存写入数据（由FileStorage::release返回）
encoding – 文件编码方式。目前不支持UTF-16 XML 编码，应使用 8-bit 编码。

2.Resharpe函数
在opencv中，reshape函数比较有意思，它既可以改变矩阵的通道数，又可以对矩阵元素进行序列化，非常有用的一个函数。
函数原型：

C++: Mat Mat::reshape(int cn, int rows=0) const

参数比较少，但设置的时候却要千万小心。
cn: 表示通道数(channels), 如果设为0，则表示保持通道数不变，否则则变为设置的通道数。

rows: 表示矩阵行数。 如果设为0，则表示保持原有的行数不变，否则则变为设置的行数。

首先设置一个初始矩阵：一个20行30列1通道的一个矩阵

3.图像的通道（channels）问题

整理一下OpenCV中文论坛里关于图像通道的问题，如下：
（1）图像的通道指的是什么？是不是灰度图的通道数为1，彩色图的通道为3？（zhuker）
正确！
基本上，描述一个像素点，如果是灰度，那么只需要一个数值来描述它，就是单通道。
如果一个像素点，有RGB三种颜色来描述它，就是三通道。（ollydbg23）
（2）对于2通道和4通道如何看待？哪位帮忙解释一下？（feixue）
我见过四通道的，两通道暂时没见过，估计只是编程的方便吧。windows的bmp有时候是一个四通道图像，R、G、B加上一个A通道，表示透明度。（eralvc）
是的，最后这个，一般叫做alpha通道，表示透明度的。（ollydbg23）
4通道通常为RGBA，在某些处理中可能会用到。
2通道图像不常见，通常在程序处理中会用到，如傅里叶变换，可能会用到，一个通道为实数，一个通道为虚数，主要是编程方便。
还有一种情况就是16位图像，本来是3通道，但是为了减少数据量，压缩为16位，刚好两个通道，常见格式有RGB555或RGB565，也就是说R占5位，G占5或6位，B占5位，也有RGBA5551格式。古老的格式，不用也罢。（Loren）
主要是有些摄像头常采用一些比较“古怪”的格式，没办法。
补充一种情况，目前常见的一些摄像头喜欢采用YUV2等格式，格式如下YUYV，在处理的时候可以用4通道或者2通道来处理。
如原格式为：Y1UY2V，插值成为Y1UV，Y2UV 就成两个彩色点了。
YCrCb也有类似压缩情况。（Loren）
详见：链接: link. 链接: link.

4 normalize()函数详解
opencv中的normalize()函数可以对向量进行范数归一化或者范围（最大最小）归一化。

void cv::normalize	(	InputArray 	src,
 	 	InputOutputArray 	dst,
 	 	double 	alpha = 1,
 	 	double 	beta = 0,
 	 	int 	norm_type = NORM_L2,
 	 	int 	dtype = -1,
 	 	InputArray 	mask = noArray() 
 	)

Parameters
src 原始数组（向量）或者矩阵（将其展开当作数组）.
dst 归一化后的结果
alpha 要归一化的范数值 或者为范围归一化的下边界
beta 范围归一化的的下边界，在范数归一化时不会使用

注意： 对于多通道数据，normalize()函数直接将其按内存中的顺序展开为数组，及当作一个向量进行处理。
链接: link.