GDAL图像重采样

#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include "gdalwarper.h"

/**
* 重采样函数(GDAL)
* @param pszSrcFile        输入文件的路径
* @param pszOutFile        写入的结果图像的路径
* @param fResX             X转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小
* @param fResY             Y转换采样比，默认大小为1.0
* @param nResampleMode     采样模式，有五种，具体参见GDALResampleAlg定义，默认为双线性内插
							GRA_NearestNeighbour=0	最近邻法，算法简单并能保持原光谱信息不变；缺点是几何精度差，灰度不连续，边缘会出现锯齿状	
							GRA_Bilinear=1				双线性法，计算简单，图像灰度具有连续性且采样精度比较精确；缺点是会丧失细节；
							GRA_Cubic=2					三次卷积法，计算量大，图像灰度具有连续性且采样精度高；
							GRA_CubicSpline=3			三次样条法，灰度连续性和采样精度最佳；
							GRA_Lanczos=4				分块兰索斯法，由匈牙利数学家、物理学家兰索斯法创立，实验发现效果和双线性接近；

* @return 0	成功；-1 打开源文件失败；-2	创建新文件失败；-3	处理过程中出错
*/
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX = 1.0, float fResY = 1.0, GDALResampleAlg eResample = GRA_Bilinear)
{
	GDALAllRegister();// 
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");//设置支持中文路径
	GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
	if (pDSrc == NULL)
	{
		return -1;
	}


	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");// 设置文件类型
	if (pDriver == NULL)
	{
		GDALClose((GDALDatasetH)pDSrc);
		return -2;
	}
	int width = pDSrc->GetRasterXSize();// 宽度
	int height = pDSrc->GetRasterYSize();// 高度
	int nBandCount = pDSrc->GetRasterCount();// 波段数
	GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();// 栅格数据类型


	char* pszSrcWKT = NULL;
	pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef()); // 投影信息


	double dGeoTrans[6] = { 0 };
	int nNewWidth = width, nNewHeight = height;
	pDSrc->GetGeoTransform(dGeoTrans); // 仿射变换 

	//adfGeoTransform[0] /* top left x */
	//adfGeoTransform[1] /* w-e pixel resolution */
	//adfGeoTransform[2] /* rotation, 0 if image is "north up" */
	//adfGeoTransform[3] /* top left y */
	//adfGeoTransform[4] /* rotation, 0 if image is "north up" */
	//adfGeoTransform[5] /* n-s pixel resolution */

	//计算采样后影像大小
	dGeoTrans[1] = dGeoTrans[1] / fResX;
	dGeoTrans[5] = dGeoTrans[5] / fResY;
	nNewWidth = static_cast<int>(nNewWidth * fResX + 0.5);
	nNewHeight = static_cast<int>(nNewHeight * fResY + 0.5);


	//创建结果数据集
	GDALDataset* pDDst = pDriver->Create(pszOutFile, nNewWidth, nNewHeight, nBandCount, dataType, NULL);
	if (pDDst == NULL)
	{
		GDALClose((GDALDatasetH)pDSrc);
		return -2;
	}

	pDDst->SetProjection(pszSrcWKT);
	pDDst->SetGeoTransform(dGeoTrans);


	void* hTransformArg = NULL;
	hTransformArg = GDALCreateGenImgProjTransformer2((GDALDatasetH)pDSrc, (GDALDatasetH)pDDst, NULL); //GDALCreateGenImgProjTransformer((GDALDatasetH) pDSrc,pszSrcWKT,(GDALDatasetH) pDDst,pszSrcWKT,FALSE,0.0,1);


	if (hTransformArg == NULL)
	{
		GDALClose((GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)pDDst);
		return -3;
	}

	GDALWarpOptions* psWo = GDALCreateWarpOptions();


	psWo->papszWarpOptions = CSLDuplicate(NULL);
	psWo->eWorkingDataType = dataType;
	psWo->eResampleAlg = eResample;


	psWo->hSrcDS = (GDALDatasetH)pDSrc;
	psWo->hDstDS = (GDALDatasetH)pDDst;


	psWo->pfnTransformer = GDALGenImgProjTransform;
	psWo->pTransformerArg = hTransformArg;


	psWo->nBandCount = nBandCount;// 波段数
	psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));// 输入文件的波段数所占用的内存
	psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));// 输出文件的波段数所占用的内存
	for (int i = 0; i < nBandCount; i++)
	{
		psWo->panSrcBands[i] = i + 1;
		psWo->panDstBands[i] = i + 1;
	}


	GDALWarpOperation oWo;
	if (oWo.Initialize(psWo) != CE_None)
	{
		GDALClose((GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)pDDst);
		return -3;
	}


	oWo.ChunkAndWarpImage(0, 0, nNewWidth, nNewHeight);


	GDALDestroyGenImgProjTransformer(hTransformArg);
	GDALDestroyWarpOptions(psWo);

	GDALFlushCache(pDDst);
	GDALClose((GDALDatasetH)pDSrc);
	GDALClose((GDALDatasetH)pDDst);


	return 0;
}

int main()
{
	ResampleGDAL("数据（测试数据）/label_RGB.tif","test.tif", 0.5, 0.5, GRA_Bilinear);
}