GDAL原始影像列均值和方差输出到txt中

在进行图像处理时，经常需要分析图像的每一行，每一列的规律，而这些规律通常是通过每一行或每一列的均值和方差来总结的，在本人的项目中也遇到了这样的事情，所以写了一个exe。exe核心代码如下：（其中的一个遍历文件夹下所有tif的函数FindFileExt 可以在我的博客中找到。

#include "stdafx.h"
#include <io.h> 
#include <iostream>
#include <fstream> 
#include <istream>
#include <ostream>
#include "FindFileExt.h"
#include "gdal_priv.h"
#include "windows.h"


using namespace std;
int _tmain(int argc, _TCHAR* argv[])
{
	return 0;
}

#define  PATH_SIZE 256
#define  RAM_USE   200
//Return normally, no stripe droped frame
#define RET_ERR_NONE 0
//Return normally,have stripe droped frame
#define RET_ERR_1 1
//Return abnormally,can't open the image
#define RET_ERR_2 2
//Return abnormally, have other error 
#define RET_ERR_3 3

double *fMean_New;
double *fSd_New;


BOOL Usage()
{

	printf("\nParameters Wrong!\n\n"
		"Img_Dect.exe [*.tiff | path]\n\n"
		"\nby huxudong 2018/1/30\n");

	return RET_ERR_3;

}
int main(int argc, char* argv[])
{

	//************************************************************************************************************
	//*******Find all the TIFF files under the specified path
	//************************************************************************************************************

	char strExe[PATH_SIZE];
	char strImg2[PATH_SIZE];
	if (argc == 2)
	{

		strcpy(strExe, argv[0]);
		strcpy(strImg2, argv[1]);

	}
	else
	{

		Usage();
		return RET_ERR_3;

	}
	char *srcimg = strImg2;
	string filePath;
	filePath = srcimg;
	vector<string> files;
	const char * distAll = "AllFiles.txt";
	string format = ".tif";
	string format1 = ".tiff";
	FindFileExt(filePath, files, format);
	FindFileExt(filePath, files, format1);
	ofstream ofn(distAll);
	int nFileNum = files.size();
	cout << "The numble of the total images is" << " "<<nFileNum << endl;
	for (int i = 0; i < nFileNum; i++)
	{

		ofn << files[i] << endl; 

	}
	FILE *fp1;
	if ((fp1 = fopen(distAll, "r")) == NULL)
	{

		printf("Open Tif File list failed:\n%s\n", distAll);

		return RET_ERR_3;

	}

	//************************************************************************************************************
	//*******Process Begining
	//************************************************************************************************************

	cout << "******************************Process Begining******************************" << endl;
	
	for (int n = 0; n < nFileNum; n++)
	{

		//**********************Get the basic information of the iamge**********************

		const char*strImg;
		strImg = files[n].c_str();
		GDALAllRegister();
		GDALDataset *ImgBef = (GDALDataset*)GDALOpen(strImg, GA_ReadOnly);
		if (ImgBef == NULL)
		{

			printf("Open Img Failed:\n%s\n", strImg);

			return RET_ERR_2;

		}

		int nCols = ImgBef->GetRasterXSize(); 	        
		int nRows = ImgBef->GetRasterYSize();
		int nBands = ImgBef->GetRasterCount();
		GDALDataType gBand = ImgBef->GetRasterBand(1)->GetRasterDataType();
		int nBits = GDALGetDataTypeSize(gBand);

		double GeoTransform[6];                     	
		ImgBef->GetGeoTransform(GeoTransform);
		const char *sProRef = ImgBef->GetProjectionRef();	

		int nStepSize = (RAM_USE * 1024 * 1024) / (nCols*nBands);
		int nStepNum = nRows / nStepSize; if (nRows%nStepSize) nStepNum++;

		//**********************If the image is 16bits,process the txt**********************

		//**********************Output the rate of the process**********************

		int perc = 0, percc = 0;
		percc = int(double(n * 100) / double(nFileNum));
		if (percc >= perc)
		{

			printf("Generating 4 TXT has done %d%%\n", percc);
			perc = percc;

		}

		//**********************If the image is 16bis,process the txt**********************

		if (nBits == 16)
		{

			int isize = GDALGetDataTypeSize(GDT_UInt16) / 8;
			for (int x = 1; x <= nBands; x++)
			{

				double *fMeanTmp = new double[nCols];
				memset(fMeanTmp, 0, nCols*sizeof(double));
				double *fStdTmp = new double[nCols];
				memset(fStdTmp, 0, nCols*sizeof(double));

				//**********************Compute the mean value of 4 bands**********************

				for (int k = 0; k < nStepNum; k++)
				{

					int ybeg = max(0, min(nStepSize*k, nRows - 1));
					int yend = max(0, min(nStepSize*(k + 1), nRows));
					GDALRasterBand*pBand = ImgBef->GetRasterBand(x);
					WORD *pbuf = new WORD[(yend - ybeg)*nCols];
					memset(pbuf, 0, ((yend - ybeg)*nCols)*sizeof(WORD));
					pBand->RasterIO(GF_Read, 0, ybeg, nCols, (yend - ybeg), pbuf, nCols, (yend - ybeg), GDT_UInt16, isize, isize*nCols);
					for (int c = 0; c < nCols; c++)
					{

						for (int b = 0; b < (yend - ybeg); b++)
						{

							fMeanTmp[c] = fMeanTmp[c] + pbuf[c + b*nCols];

						}

					}

					delete[]pbuf; pbuf = NULL;

				}
				for (int c = 0; c < nCols; c++)
				{

					fMeanTmp[c] = double(fMeanTmp[c] / nRows);

				}

				//**********************Compute the variance value of 4 bands**********************

				for (int k = 0; k < nStepNum; k++)
				{

					int ybeg = max(0, min(nStepSize*k, nRows - 1));
					int yend = max(0, min(nStepSize*(k + 1), nRows));
					GDALRasterBand*pBand = ImgBef->GetRasterBand(x);
					WORD *pbuf = new WORD[(yend - ybeg)*nCols];
					memset(pbuf, 0, ((yend - ybeg)*nCols)*sizeof(WORD));
					pBand->RasterIO(GF_Read, 0, ybeg, nCols, (yend - ybeg), pbuf, nCols, (yend - ybeg), GDT_UInt16, isize, isize*nCols);
					for (int c = 0; c < nCols; c++)
					{

						for (int b = 0; b < (yend - ybeg); b++)
						{

							fStdTmp[c] = fStdTmp[c] + (pbuf[c + b*nCols] - fMeanTmp[c])*(pbuf[c + b*nCols] - fMeanTmp[c]);

						}

					}	

					delete[]pbuf; pbuf = NULL;

				}

				for (int c = 0; c < nCols; c++)
				{

					fStdTmp[c] = double(sqrt(fStdTmp[c] / nRows));

				}

				//*****************Write the mean value and  variance to the 4 txts*****************

				string strTxt = strImg;
				int pos1 = strTxt.find_last_of('.');
				strTxt.erase(pos1);
				string bandTh = std::to_string(x);
				strTxt += bandTh;
				strTxt += "_band.txt";
				ofstream outTxtFile(strTxt);
				for (int r = 0; r < nCols; r++)
				{

					outTxtFile << r << " "<< fMeanTmp[r] <<" "<< fStdTmp[r] << endl;

				}

				outTxtFile.close();

				//*********************************release the RAM**********************************

				delete[]fMeanTmp; fMeanTmp = nullptr;
				delete[]fStdTmp; fStdTmp = nullptr;

			}
		}

		//*****************If the image is not 16bis,jump the outside of this circle*****************

		if(nBits != 16)
		{

			continue;

		}

		GDALClose(ImgBef);

	}

	//************************************************************************************************************
	//*******Process ending
	//************************************************************************************************************

	cout << "Generating 4 TXT has done 100%\n\n" << endl;
	cout << "******************************Process ending******************************\n\n" << endl;
	fclose(fp1);

}