在编写重采样图像时,可以使用GDAL来读写图像,然后自己编写重采样算法(最邻近像元法,双线性内插法,三次立方卷积法等)【关于这采样算法有时间我会单独写一篇文章来说明原理的】将计算的结果写入图像中来实现。
在GDAL的算法中,已经提供了五种重采样算法,其定义如下(位置gdalwarper.h 的46行):
/*! Warp Resampling Algorithm */
typedef enum {
/*! Nearest neighbour (select on one input pixel) */ GRA_NearestNeighbour=0,
/*! Bilinear (2x2 kernel) */ GRA_Bilinear=1,
/*! Cubic Convolution Approximation (4x4 kernel) */ GRA_Cubic=2,
/*! Cubic B-Spline Approximation (4x4 kernel) */ GRA_CubicSpline=3,
/*! Lanczos windowed sinc interpolation (6x6 kernel) */ GRA_Lanczos=4
} GDALResampleAlg;
在查看Gdalwarp的源代码发现,warp的功能非常强大,可以用来做投影转换,重投影,投影定义,重采样,镶嵌,几何精校正和影像配准等。一句话,很好很强大。下面就看看其中的一点点皮毛,使用warp来编写一个重采样的接口,代码如下:
/**
* 重采样函数(GDAL)
* @param pszSrcFile 输入文件的路径
* @param pszOutFile 写入的结果图像的路径
* @param fResX X转换采样比,默认大小为1.0,大于1图像变大,小于1表示图像缩小
* @param fResY Y转换采样比,默认大小为1.0
* @param nResampleMode 采样模式,有五种,具体参见GDALResampleAlg定义,默认为双线性内插
* @param pExtent 采样范围,为NULL表示计算全图
* @param pBandIndex 指定的采样波段序号,为NULL表示采样全部波段
* @param pBandCount 采样的波段个数,同pBandIndex一同使用,表示采样波段的个数
* @param pszFormat 写入的结果图像的格式
* @param pProgress 进度条指针
* @return 成功返回0,否则为其他值
*/
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX , float fResY, LT_ResampleMode nResampleMode,
LT_Envelope* pExtent, int* pBandIndex, int *pBandCount, const char* pszFormat, LT_Progress *pProgress)
{
if(pProgress != NULL)
{
pProgress->SetProgressCaption("重采样");
pProgress->SetProgressTip("正在执行重采样...");
}
GDALAllRegister();
GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile, GA_ReadOnly);
if (pDSrc == NULL)
{
if(pProgress != NULL)
pProgress->SetProgressTip("指定的文件不存在,或者该格式不被支持!");
return RE_NOFILE;
}
GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if (pDriver == NULL)
{
if(pProgress != NULL)
pProgress->SetProgressTip("不能创建该格式的文件!");
GDALClose((GDALDatasetH) pDSrc);
return RE_CREATEFILE;
}
int iBandCount = pDSrc->GetRasterCount();
string strWkt = pDSrc->GetProjectionRef();//返回坐标系统
GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
double dGeoTrans[6] = {0};
pDSrc->GetGeoTransform(dGeoTrans); //六个参数其实是图像行列号坐标和地理坐标转换系数
int iNewBandCount = iBandCount;
if (pBandIndex != NULL && pBandCount != NULL)
{
int iMaxBandIndex = pBandIndex[0]; //找出最大的波段索引序号
for (int i=1; i<*pBandCount; i++)
{
if (iMaxBandIndex < pBandIndex[i])
iMaxBandIndex = pBandIndex[i];
}
if(iMaxBandIndex > iBandCount)
{
if(pProgress != NULL)
pProgress->SetProgressTip("指定的波段序号超过图像的波段数,请检查输入参数!");
GDALClose((GDALDatasetH) pDSrc);
return RE_PARAMERROR;
}
iNewBandCount = *pBandCount;
}
LT_Envelope enExtent;
enExtent.setToNull();
if (pExtent == NULL) //全图计算
{
double dPrj[4] = {0}; //x1,x2,y1,y2
ImageRowCol2Projection(dGeoTrans, 0, 0, dPrj[0], dPrj[2]); //将图像左上角行列号(0,0)坐标转为地理坐标
ImageRowCol2Projection(dGeoTrans, pDSrc->GetRasterXSize(), pDSrc->GetRasterYSize(), dPrj[1], dPrj[3]);//将图像右下角的坐标转为地理坐标
enExtent.init(dPrj[0], dPrj[1], dPrj[2], dPrj[3]);
pExtent = &enExtent;
}
dGeoTrans[0] = pExtent->getMinX();
dGeoTrans[3] = pExtent->getMaxY(); //dGeoTrans[0],dGeoTrans[3]表示的是左上角的地理坐标
dGeoTrans[1] = dGeoTrans[1] / fResX;
dGeoTrans[5] = dGeoTrans[5] / fResY; //dGeoTrans[1],dGeoTrans[5]表示的是图像横向和纵向的分辨率
int iNewWidth = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[1]) + 0.5) );
int iNewHeight = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[5]) + 0.5) );
GDALDataset *pDDst = pDriver->Create(pszOutFile, iNewWidth, iNewHeight, iNewBandCount, dataType, NULL);
if (pDDst == NULL)
{
if(pProgress != NULL)
pProgress->SetProgressTip("创建输出文件失败!");
GDALClose((GDALDatasetH) pDSrc);
return RE_CREATEFILE;
}
pDDst->SetProjection(strWkt.c_str());
pDDst->SetGeoTransform(dGeoTrans);
GDALResampleAlg eResample = (GDALResampleAlg) nResampleMode;
if(pProgress != NULL)
{
pProgress->SetProgressTip("正在执行重采样...");
pProgress->SetProgressTotalStep(iNewBandCount*iNewHeight);
}
int *pSrcBand = NULL;
int *pDstBand = NULL;
int iBandSize = 0;
if (pBandIndex != NULL && pBandCount != NULL)
{
iBandSize = *pBandCount;
pSrcBand = new int[iBandSize];
pDstBand = new int[iBandSize];
for (int i=0; i<iBandSize; i++)
{
pSrcBand[i] = pBandIndex[i];
pDstBand[i] = i+1;
}
}
else
{
iBandSize = iBandCount;
pSrcBand = new int[iBandSize];
pDstBand = new int[iBandSize];
for (int i=0; i<iBandSize; i++)
{
pSrcBand[i] = i+1;
pDstBand[i] = i+1;
}
}
void *hTransformArg = NULL, *hGenImgPrjArg = NULL;
hTransformArg = hGenImgPrjArg = GDALCreateGenImgProjTransformer2((GDALDatasetH) pDSrc, (GDALDatasetH) pDDst, NULL);
if (hTransformArg == NULL)
{
if(pProgress != NULL)
pProgress->SetProgressTip("转换参数错误!");
GDALClose((GDALDatasetH) pDSrc);
GDALClose((GDALDatasetH) pDDst);
return RE_PARAMERROR;
}
GDALTransformerFunc pFnTransformer = GDALGenImgProjTransform;
GDALWarpOptions *psWo = GDALCreateWarpOptions();
psWo->papszWarpOptions = CSLDuplicate(NULL);
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH) pDSrc;
psWo->hDstDS = (GDALDatasetH) pDDst;
psWo->pfnTransformer = pFnTransformer;
psWo->pTransformerArg = hTransformArg;
psWo->pfnProgress = GDALProgress;
psWo->pProgressArg = pProgress;
psWo->nBandCount = iNewBandCount;
psWo->panSrcBands = (int *) CPLMalloc(iNewBandCount*sizeof(int));
psWo->panDstBands = (int *) CPLMalloc(iNewBandCount*sizeof(int));
for (int i=0; i<iNewBandCount; i++)
{
psWo->panSrcBands[i] = pSrcBand[i];
psWo->panDstBands[i] = pDstBand[i];
}
RELEASE(pSrcBand);
RELEASE(pDstBand);
GDALWarpOperation oWo;
if (oWo.Initialize(psWo) != CE_None)
{
if(pProgress != NULL)
pProgress->SetProgressTip("转换参数错误!");
GDALClose((GDALDatasetH) pDSrc);
GDALClose((GDALDatasetH) pDDst);
return RE_PARAMERROR;
}
oWo.ChunkAndWarpImage(0, 0, iNewWidth, iNewHeight);
GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
GDALDestroyWarpOptions( psWo );
GDALClose((GDALDatasetH) pDSrc);
GDALClose((GDALDatasetH) pDDst);
if(pProgress != NULL)
pProgress->SetProgressTip("重采样完成!");
return RE_SUCCESS;
}