c++中使用GDAL将两张影像相减

void image_process::do_diff(int master_date, int slave_date, string img_path, string out_path)
{

	img_path.append("\\");
	out_path.append("\\");
	string filepath1 = img_path + to_string(master_date) + ".tif";
	string filepath2 = img_path + to_string(slave_date) + ".tif";

	GDALAllRegister(); //注册驱动
	GDALDataset* poSrcDS;//建立数据集读第一张图片
	poSrcDS = (GDALDataset*)GDALOpen(filepath1.c_str(), GA_ReadOnly);//指针指向文件路径
	if (poSrcDS == NULL)
	{
		cout << "open file error!" << endl;
	}

	int	iXSizeDS = poSrcDS->GetRasterXSize();//获取图像的宽
	int iYSizeDS = poSrcDS->GetRasterYSize();//获取图像的高
	int iBandCountDS = poSrcDS->GetRasterCount();//获取图像的波段数

	GDALDataset* poSrcDF;//建立数据集读第二张图片
	poSrcDF = (GDALDataset*)GDALOpen(filepath2.c_str(), GA_ReadOnly);//指针指向文件路径
	if (poSrcDF == NULL)
	{
		cout << "open file error!" << endl;
	}

	int	iXSizeDF = poSrcDF->GetRasterXSize();//获取图像的宽
	int iYSizeDF = poSrcDF->GetRasterYSize();//获取图像的高
	int iBandCountDF = poSrcDF->GetRasterCount();//获取图像的波段数

	GDALDriver* poDriver;
	GDALDataset* poDstDS;
	const char* pszFormat = "GTiff";
	//const char* pszDstFile = "NewPic.tiff";
	string pszDstFile = out_path + to_string(master_date) + "-" + to_string(slave_date) + ".ztd.tif";
	poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
	poDstDS = poDriver->Create(pszDstFile.c_str(), iXSizeDS, iYSizeDS, iBandCountDS, GDT_Float32, NULL);//创建文件用来写入计算后的数据
	double value[6];
	poSrcDF->GetGeoTransform(value);
	poDstDS->SetGeoTransform(value);
	float* pSrcData = new float[iXSizeDS * iYSizeDS];//创建指针存第一张图片的数据
	float* pDstData = new float[iXSizeDS * iYSizeDS];//创建指针存结果图像的数据
	float* pFrcData = new float[iXSizeDS * iYSizeDS];//创建指针存第二张图片的数据


	if (iXSizeDS == iXSizeDF && iYSizeDS == iYSizeDF && iBandCountDS == iBandCountDF)//鲁棒
	{

		//循环波段,按波段来读写图像
		for (int iBand = 1; iBand <= iBandCountDS; iBand++)
		{
			GDALRasterBand* pSrcBandDS = poSrcDS->GetRasterBand(iBand);//获取第一张图片的第iBand个波段
			GDALRasterBand* pSrcBandDF = poSrcDF->GetRasterBand(iBand);//获取第二张图片的第iBand个波段
			GDALRasterBand* pDstBand = poDstDS->GetRasterBand(iBand);
			for (int i = 0; i < iYSizeDS; i++)//循环图像高
			{
				pSrcBandDS->RasterIO(GF_Read, 0, i, iXSizeDS, 1, pSrcData, iXSizeDS, 1, GDT_Float32, 0, 0);
				pSrcBandDF->RasterIO(GF_Read, 0, i, iXSizeDS, 1, pFrcData, iXSizeDS, 1, GDT_Float32, 0, 0);
				for (int j = 0; j < iXSizeDS; j++)//循环图像宽
				{
					pDstData[j] = pSrcData[j] - pFrcData[j];
					pDstBand->RasterIO(GF_Write, 0, i, iXSizeDS, 1, pDstData, iXSizeDS, 1, GDT_Float32, 0, 0);//把数据存入pDstData
				}
			}
		}
	}
	
	GDALClose((GDALDatasetH)poSrcDS);//关闭第一张图片
	GDALClose((GDALDatasetH)poSrcDF);//关闭第二张图片
	GDALClose((GDALDatasetH)poDstDS);//关闭结果图片
	delete[]pSrcData;
	delete[]pFrcData;
	delete[]pDstData;
}