GDAL 做影像校正, 支持RPC和GCP

由于公司基础软件一直缺乏影像校正能力,遂基于GDAL写了一个小程序,可以做小数据的影像校正, 网上搜到的都是骤点校正,速度太慢,本程序使用全内存模式, 校正个几个G的数据 的效率比原生GDAL提供的gdalWarp.exe还要快,主要代码如下:

主要函数声明:

	int ImageWarpRCPorGCP2(GDALDatasetH hSrcDS,
		const char * pszDstFile,
		RasterWarpInfo & info,
		int blockw = 256,
		int blockh = 256,
		int transtype = 0, /*0 RPC, 1 GCP*/
		int iOrder = 0,
		double dResX = 0.0,
		double dResY = 0.0,
		GDALRIOResampleAlg eResampleMethod = GRIORA_NearestNeighbour,
		const char * pszFormat = "GTiff");

结构体声明:

///\brief 栅格类基本信息
struct  RasterWarpInfo
{

	///\brief src X，Y维度范围
	double			XYDomain[4] = {0};
	///\brief dst X，Y维度范围
	double			dstXYDomain[4] = { 0 };

	///\brief src X，Y维度范围
	double			Transform[6] = { 0 };
	///\brief dst X，Y维度范围
	double			dstTransform[6] = { 0 };

	void*	pGDALRPCorGCPTransform = 0;
	///\brief 像素宽度
	int				Width =0;
	///\brief 像素高度
	int				Height = 0;

	///\brief 像素宽度
	int				dstWidth = 0;
	///\brief 像素高度
	int				dstHeight = 0;

	///\brief 数据存储块宽度大小
	int				BlockWidth = 0;
	///\brief 数据存储块高度大小
	int				BlockHeight =0;

	///\brief 波段数据类型
	GDALDataType DataType = GDT_Byte;

	///\brief 波段的名称
	std::vector<GDALColorInterp> BandTypes = {};;
	//空间参考
	char * srcRef = 0;;
	char * dstRef = 0;;

	int nPixelSize = 0;

	std::vector<int>	m_vBand = {};
};

 

主要实现:

/**
* @brief 影像几何校正
* @param GDALDatasetH hSrcDS,			输入数据集
* @param pszDstFile			输出文件路径
* @param RasterWarpInfo		校正信息
* @param	int blockw ,分块大小
* @param	int blockh ,分块大小
* @param	int transtype,RPC或者GCP
* @param	int iOrder,
* @param	double dResX, 分辨率  目前为原始分辨率
* @param	double dResY 分辨率 目前为原始分辨率
* @param
* @param eResampleMethod	重採样方式
* @param @param pszFormat			输出文件格式
* * @return  返回值，返回true或者false
* */
int CImageCorrectionDlg::ImageWarpRCPorGCP2(GDALDatasetH hSrcDS,
	const char * pszDstFile,
	RasterWarpInfo & info,
	int blockw ,
	int blockh ,
	int transtype,
	int iOrder,
	double dResX,
	double dResY,
	GDALRIOResampleAlg eResampleMethod,
	const char * pszFormat)
{

	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	GDALDataType eDataType = info.DataType;
	int nBandCount = info.BandTypes.size();

	// 创建几何多项式坐标转换关系
	void *hTransform = info.pGDALRPCorGCPTransform;
	if (NULL == hTransform)
	{
		GDALClose(hSrcDS);
		return 0;
	}

	// 计算输出图像四至范围、大小、仿射变换六參数等信息
	double* adfGeoTransform = info.dstTransform;
	double* adfExtent = info.dstXYDomain;
	int    nPixels = info.dstWidth, nLines = info.dstHeight;

	// 以下開始依据用户指定的分辨率来反算输出图像的大小和六參数等信息
	//double dResXSize = dResX;
	//double dResYSize = dResY;

	//dResXSize = (info.dstTransform[1]);
	//dResYSize = (info.dstTransform[5]);


	// 计算输出图像的范围
	//double minX = adfGeoTransform[0];
	//double maxX = adfGeoTransform[0] + adfGeoTransform[1] * nPixels;
	//double maxY = adfGeoTransform[3];
	//double minY = adfGeoTransform[3] + adfGeoTransform[5] * nLines;

	//nPixels = ceil((maxX - minX) / dResXSize);
	//nLines = abs(ceil((minY - maxY) / dResYSize));
	//adfGeoTransform[0] = minX;
	//adfGeoTransform[3] = maxY;
	//adfGeoTransform[1] = dResXSize;
	//adfGeoTransform[5] = dResYSize;

	// 创建输出图像
	GDALDriverH hDriver = GDALGetDriverByName(pszFormat);
	if (NULL == hDriver)
	{
		::MessageBox(this->m_hWnd, _T("获取驱动失败"), _T("错误"),
			MB_ICONERROR);
		return  0;
	}
	GDALDatasetH hDstDS = GDALCreate(hDriver, pszDstFile, nPixels, nLines, nBandCount, eDataType, NULL);
	if (NULL == hDstDS)
	{
		::MessageBox(this->m_hWnd, _T("创建目标文件失败"), _T("错误"),
			MB_ICONERROR);
		return  0;
	}
	GDALDataset* gsDS = (GDALDataset*)hSrcDS;
	char* pszWkt1 = (char*)gsDS->GetProjectionRef();
	if (0 == strlen(pszWkt1))
	{
		pszWkt1 = (char*)gsDS->GetGCPProjection();
		if (0 == strlen(pszWkt1))
		{
			pszWkt1 = "GEOGCS[\"GCS_WGS_1984\", DATUM[\"D_WGS_1984\", SPHEROID[\"WGS_1984\", 6378137, 298.257223563]], PRIMEM[\"Greenwich\", 0], UNIT[\"degree\", 0.0174532925199433]]";
		}
	}
	GDALSetProjection(hDstDS, pszWkt1);
	GDALSetGeoTransform(hDstDS, adfGeoTransform);// info.dstTransform);



	//获得原始图像的行数和列数
	int nXsize = GDALGetRasterXSize(hSrcDS);
	int nYsize = GDALGetRasterYSize(hSrcDS);

	//然后是图像重採样
	int nFlag = 0;
	float dfValue = 0;
	CPLErr err = CE_Failure;
	//进度
	INIT_RASTERIO_EXTRA_ARG(m_Arg);
	m_Arg.eResampleAlg = GRIORA_NearestNeighbour;
	m_Arg.pfnProgress = RasterIOProgress;
	//m_Arg.pProgressData = pClass;

	long long srcLen = (long long)info.nPixelSize* info.Width * info.Height;
	unsigned char* srcbuff = (unsigned char*)malloc(srcLen);

	long long targetLen = (long long)info.nPixelSize * info.dstWidth* info.dstHeight;
	unsigned char* tgbuff = (unsigned char*)malloc(targetLen);
	memset(tgbuff, 0, sizeof(unsigned char));
	GDALDataset* srcDS = (GDALDataset*)hSrcDS;
	GDALDataset* tgDS = (GDALDataset*)hDstDS;
	for (int i = 0; i < tgDS->GetRasterCount(); i++)
	{
		tgDS->GetRasterBand(i + 1)->SetNoDataValue(0.);
	}
	//直接读取所有原始数据,  
	 err = srcDS->RasterIO(GF_Read, 0, 0, info.Width, info.Height, srcbuff,
		info.Width, info.Height, info.DataType, info.BandTypes.size(), &info.m_vBand[0],
		info.nPixelSize, info.nPixelSize * info.Width,
		info.nPixelSize / info.m_vBand.size(), &m_Arg);

	//unsigned char* emptyvalue = (unsigned char*)malloc(info.nPixelSize*8);
	//memset(emptyvalue, 0, info.nPixelSize * 8);
	//遍历目标像素值,给目标像素buff赋值

	 nPixels = info.dstWidth; 
	 nLines = info.dstHeight;
#pragma omp parallel for
	for (int nRow = 0; nRow < nLines; nRow++)
	{
		for (int nCol = 0; nCol < nPixels; nCol++)
		{
			double dbX = adfGeoTransform[0] + nCol * adfGeoTransform[1]
				+ nRow * adfGeoTransform[2];
			double dbY = adfGeoTransform[3] + nCol * adfGeoTransform[4]
				+ nRow * adfGeoTransform[5];

			//由输出的图像地理坐标系变换到原始的像素坐标系

			if (0 == transtype)
				GDALRPCTransform(hTransform, TRUE, 1, &dbX, &dbY, NULL, &nFlag);
			else if (1 == transtype)
				GDALGCPTransform(hTransform, TRUE, 1, &dbX, &dbY, NULL, &nFlag);
			int nXCol = (int)(dbX );
			int nYRow = (int)(dbY);

			//超出范围的用0填充
			if (nXCol < 0 || nXCol >= nXsize || nYRow < 0 || nYRow >= nYsize)
			{
				int tmp = 0;
			
				tgbuff[nRow * nPixels * info.nPixelSize + nCol* info.nPixelSize]=0;
			}
			else
			{
				 memcpy(&tgbuff[nRow * nPixels * info.nPixelSize + nCol * info.nPixelSize],
					 &srcbuff[nYRow   * nXsize * info.nPixelSize + nXCol * info.nPixelSize],
					 info.nPixelSize);
	
				 //memcpy(&tgbuff[nCol*nLines + nRow], &srcbuff[nXCol  * nXsize + nYRow], sizeof(short)*4);

			}
		}
	}
	//free(emptyvalue);
	//写好的buff直接IO进文件
	tgDS->RasterIO(GF_Write,0,0,info.dstWidth,info.dstHeight, tgbuff,
		info.dstWidth, info.dstHeight, info.DataType, info.BandTypes.size(), &info.m_vBand[0],
		info.nPixelSize, info.nPixelSize * info.dstWidth,
		info.nPixelSize / info.m_vBand.size(), &m_Arg);

	//err = tgDS->RasterIO(GF_Write, 0, 0, info.Width, info.Height, tgbuff,
	//	info.Width, info.Height, info.DataType, info.BandTypes.size(), &info.m_vBand[0],
	//	0,0,0, &m_Arg);

	if (hTransform != NULL)
	{
		GDALDestroyGCPTransformer(hTransform);
		hTransform = NULL;
	}

	GDALClose(hSrcDS);
	GDALClose(hDstDS);


	return 1;
}

 

 

RasterWrapInfo 的获取如下:

int CImageCorrectionDlg::ImageCorrection(const char * pszSrcFile, const char * pszDstFile)
{	//数据集。
	GDALDataset* m_pDS = (GDALDataset *)GDALOpen(pszSrcFile, GA_Update);
	if (m_pDS == NULL)
	{
		::MessageBox(this->m_hWnd, _T("打开文件失败"), _T("错误"),
			MB_ICONERROR);
		return  0;
	}
	if (m_pDS == NULL)
		return 0; 
	//获取原始信息
	m_pDS->GetGeoTransform(m_GeoTransform);
	RasterWarpInfo srcInfo;
	//memset(&srcInfo, 1, sizeof(srcInfo));
	/// \brief 像素宽度
	srcInfo.Width = m_pDS->GetRasterXSize();
	/// \brief 像素高度
	srcInfo.Height = m_pDS->GetRasterYSize();
	for (int i = 0; i < m_pDS->GetRasterCount(); i++)
	{
		srcInfo.m_vBand.emplace_back(i + 1);
		srcInfo.DataType = m_pDS->GetRasterBand(i + 1)->GetRasterDataType();
		srcInfo.BandTypes.emplace_back(m_pDS->GetRasterBand(i + 1)->GetColorInterpretation());
		m_pDS->GetRasterBand(i + 1)->GetBlockSize(&srcInfo.BlockWidth, &srcInfo.BlockHeight);
		srcInfo.nPixelSize += GDALGetDataTypeSize(srcInfo.DataType);
	}
	
	srcInfo.nPixelSize /= 8;

	char ** datalist = m_pDS->GetMetadataDomainList();

	char** papszRPC = m_pDS->GetMetadata("RPC");
	if (papszRPC)
	{
		GDALRPCInfo oInfo;
		GDALExtractRPCInfo(papszRPC, &oInfo);
		//设置RPC模型中所需的DEM路径
		char** papszTransOption = NULL;
		double adfGeoTransform[6] = { 0 };

		int    nPixels = 0, nLines = 0;
		//使用RPC信息，DEM等构造RPC转换参数
		m_pGDALRPCorGCPTransform = GDALCreateRPCTransformer(&oInfo, FALSE, 0, papszTransOption);
		if (m_pGDALRPCorGCPTransform != NULL)
		{
			if (GDALSuggestedWarpOutput2(m_pDS, GDALRPCTransform/*GDALRPCTransform*/, m_pGDALRPCorGCPTransform,
				m_RPCorGCPTransform, &m_nPixels, &m_nLines, m_adfExtent, 0) != CE_None)
			{

				GDALClose(m_pGDALRPCorGCPTransform);
				return 0;
			}
			srcInfo.dstHeight = m_nLines;
			srcInfo.dstWidth = m_nPixels;
			memcpy(srcInfo.dstXYDomain, m_adfExtent, sizeof(double) * 4);
			memcpy(srcInfo.dstTransform, m_RPCorGCPTransform, sizeof(double) * 6);
			srcInfo.pGDALRPCorGCPTransform = m_pGDALRPCorGCPTransform;
		}

		int nGCPCount = m_pDS->GetGCPCount();

		const GDAL_GCP *pGCPList = m_pDS->GetGCPs();
		if (nGCPCount > 0)
		{
			m_pGDALRPCorGCPTransform = GDALCreateGCPTransformer(nGCPCount, pGCPList, 0, FALSE);

			if (GDALSuggestedWarpOutput2(m_pDS, GDALGCPTransform, m_pGDALRPCorGCPTransform,
				m_RPCorGCPTransform, &m_nPixels, &m_nLines, m_adfExtent, 0) != CE_None)
			{
				GDALClose(m_pDS);
				return 0;
			}
			srcInfo.dstHeight = m_nLines;
			srcInfo.dstWidth = m_nLines;
			memcpy(srcInfo.dstXYDomain, m_adfExtent, sizeof(double) * 4);
			memcpy(srcInfo.dstTransform, m_RPCorGCPTransform, sizeof(double) * 6);
			srcInfo.pGDALRPCorGCPTransform = m_pGDALRPCorGCPTransform;
		}
	}

	//如果没有空间参考则读取空间参考
	const char*pSrsRef = m_pDS->GetProjectionRef();
	char* pszWkt1 = (char*)m_pDS->GetGCPProjection();
	

	//ImageWarpRCPorGCP(m_pDS, pszDstFile, srcInfo);

	return ImageWarpRCPorGCP2(m_pDS, pszDstFile, srcInfo);
}

 

 

此程序未做分块处理,分块和这整块也无太大区别, 因为是测试程序,所以不做过多代码优化,源码下载请到我的资源下载,只不过是个mfc的测试 程序,代码较乱,请理解