由于公司基础软件一直缺乏影像校正能力,遂基于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的测试 程序,代码较乱,请理解