只遍历一次图像如何求图像个波段之间的协方差矩阵

根据图像的协方差矩阵的求解公式，采用传统的方法至少要遍历2次图像才能求得图像的协方差矩阵：

第一次是求均值，第二次是求每个样本与均值样本之差的积。

如果图像比较大，存储在磁盘上，频繁的读取磁盘数据是很费时的一项工作。

有没有只遍历一次图像就能求得图像的协方差矩阵的方法呢？

经过研究还真发现了这样的一种方法，废话少说直接上代码：

        GDALAllRegister();                                  //利用GDAL读取图片，先要进行注册
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");   //设置支持中文
	//准备读取图片
	GDALDataset *ReadDataSet=(GDALDataset*)GDALOpen(inputFileName,GA_ReadOnly);
	if(ReadDataSet==NULL) return -1;

	int x=0;                                     //图像的起始位置
	int y=0;
	int width=ReadDataSet->GetRasterXSize();     //图像的宽度
	int height=ReadDataSet->GetRasterYSize();    //图像的高度
	int bandsCount=ReadDataSet->GetRasterCount();//图像的波段数

	if(bandsCount<2)
	{
		delete ReadDataSet;ReadDataSet=NULL;
		return -3;       //波段个数必须大于等于2个
	}
	//缓存波段列表
	int *bandsList=new int[bandsCount];
	if (dims!=NULL && dims->bandsIndex!=NULL)
	{
		for(int i=0;i<bandsCount;i++)
		{
			bandsList[i]=dims->bandsIndex[i];
		}
	}
	else
	{
		for(int i=0;i<bandsCount;i++)
		{
			bandsList[i]=i+1;
		}
	}//缓存波段列表完毕

	int ppt=0;
	int progress=0;
	//协方差矩阵
	long pixelsCount=width*height;
	//有效像元的个数，不包括背景像元
	long validPixelsCount=0;
	//分配协方差矩阵内存
	double *covMat=new double[bandsCount*bandsCount]();	
	//分配波段数据平均值的内存
	double *bandsAverage=new double[bandsCount];
	{//协方差矩阵的计算
		//用于保存每行像元值的内存
		float *rowData=new float[width*bandsCount];
		//用于保存每个像元值的内存
		double *pixelData=new double[bandsCount];
		//分配波段数据之和的内存
		double *bandsSum=new double[bandsCount]();		
		//开始循环
		for(int r=0;r<height;++r)
		{
			// 获取一行数据
			if(ReadDataSet->RasterIO(GF_Read,x,y+r,width,1,rowData,width,1,
				GDT_Float32,bandsCount,bandsList,0,0,0)==CE_Failure )
			{
				delete ReadDataSet;
				delete[] covMat;
				delete[] bandsAverage;
				delete[] rowData;
				delete[] pixelData;
				delete[] bandsSum;
				delete[] bandsList;
				return -1;
			}	  
			//循环一行中的每一个像元
			for (int w=0;w<width;w++)
			{
				//判断该像元是否是背景像元
				if(fabs(rowData[w]-backgroundValue)>FLT_EPSILON)
				{
					//有效像元个数加1
					validPixelsCount++;
					//获取一个像元的数据
					for(int n=0;n<bandsCount;n++)
					{
						pixelData[n]=rowData[n*width+w];
					}
					//计算协方差矩阵
					for (int h=0;h<bandsCount;++h)
					{
						for (int l=h;l<bandsCount;++l)
						{
							covMat[h*bandsCount+l]+=pixelData[h]*pixelData[l];
						}
						bandsSum[h]+=pixelData[h];
					}
				}
			}
		}
		// 计算均值
		for(int n=0;n<bandsCount;n++)
		{
			bandsAverage[n]=bandsSum[n]/validPixelsCount;
		}
		//计算协方差矩阵
		for (int h=0;h<bandsCount;++h)
		{
			for (int l=h;l<bandsCount;++l)
			{
				covMat[h*bandsCount+l]=covMat[h*bandsCount+l]-bandsSum[h]*bandsAverage[l]
				-bandsSum[l]*bandsAverage[h]+bandsAverage[h]*bandsAverage[l]*validPixelsCount;
			}
		}
		// 将协方差矩阵的和除以像元个数，并将上三角复制成下三角
		for (int h=0;h<bandsCount;++h)
		{
			for (int l=h;l<bandsCount;++l)
			{
				covMat[h*bandsCount+l]/=(validPixelsCount-1);
				covMat[l*bandsCount+h]=covMat[h*bandsCount+l];
			}
		}
		//删除临时内存
		delete[] pixelData;
		delete[] rowData;
		delete[] bandsSum;
	}//协方差矩阵计算完毕
	delete ReadDataSet;
	delete[] covMat;
	delete[] bandsAverage;
	delete[] bandsList;
	return 0;

代码是从项目中摘抄出来的，要用的话可能的自己修改一下，用到了GDAL读取图像数据。有问题欢迎交流！