1 前言
之前在写影像融合算法的时候，免不了要实现将多光谱影像重采样到全色大小。当时为了不影响融合算法整体开发进度，其中重采样功能用的是GDAL开源库中的Warp接口实现的。
后来发现GDAL Warp接口实现的多光谱到全色影像的重采样主要存在两个问题：1 与原有平台的已有功能不兼容，产生冲突；2 效率较低。因此，决定重新设计和开发一个这样的功能，方便后期软件系统的维护等。
2 图像重采样
图像处理从形式上来说主要包括两个方面：1 单像素或者邻域像素的处理，比如影像的相加或者滤波运算等；2 图像几何空间变换，如图像的重采样，配准等。
影像重采样的几何空间变换公式如下：
其中 为变换系数，常用的重采样算法主要包括以下三种：1 最近邻；2 双线性；3 三次卷积。
2.1 最近邻采样
最近邻采样的原理概况起来就是用采样点位置最近的一个像素值替代采样点位置的像素值。在这里插入一点：
通常图像空间变换有两种方法，直接法或者间接法。以图像重采样为例说明如下：直接法：从原始的图像行列初始值开始，根据变换公式，计算采样后的像素位置，并对位置赋值，但是这种方法会出现，原始图像的多个像素点对应到同一采样后的像素点，从而还要增加额外方法进行处理；间接法：是从重采样后图像的行列初始值开始，计算得到其在原始影像中的位置，并根据一定的算法进行计算，得到采样后的值。这种方法简单直接，本文就是采用这样的方法。
最近邻采样的实现算法如下：
首先遍历采样点，根据公式计算采样点在原始图像中的位置，假设位置为 。然后计算与 最近的点，并将其像素值赋为采样点的像素值。其公式如下：
2.2 双线性
双线性采样和最近邻赋值不同，是找到 最近的四个像素点，然后将距离作为权重分别插值 和 方向，从而插值到采样点位置。具体公式见代码部分。
2.3 三次卷积
三次卷积是利用 最近的16个像素点进行插值计算得到。同样也是分别插值 和 方向。具体公式见代码部分。
3 实验结果
主要实现了两个版本的差值结果：1 CPU版本；2 GPU版本（初学）。考虑到多光谱和全色影像范围大小不一致的事实，算法首先计算全色和多光谱的重叠区域。话不多说，先看看详细的代码实现过程。
CPU版本：
#ifndef RESAMPLECPU_H
#define RESAMPLECPU_H #include <gdal_alg_priv.h>
#include <gdal_priv.h>
#include <assert.h> template<typename T>
void ReSampleCPUKernel(const float *mssData,
T *resampleData,
int mssWidth,
int mssHeight,
int mssBandCount,
int mssOffsetX,
int mssOffsetY,
int panWidth,
int panHeight,
float radioX,
float radioY,
double dfDstNoDataValue,
int MethodType)
{
assert(mssData != NULL);
float eps = 0.00001f;
for(int idx = ;idx < panHeight;idx++){
for(int idy = ;idy<panWidth;idy++){
// 找到对应的MSS像素位置
float curX = (float)idx * radioX;
float curY = (float)idy * radioY;
if(mssData[int(curX)*mssWidth*mssBandCount + int(curY)] == dfDstNoDataValue)
{
int i = ;
while(i < mssBandCount){
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(dfDstNoDataValue);
i++;
}
continue;
}
if(MethodType == ){ // 最近邻
int nearX = (int)(curX + 0.5)>(int)curX?(int)(curX + ):(int)curX;
int nearY = (int)(curY + 0.5)>(int)curY?(int)(curY + ):(int)curY;
if(nearX >= mssHeight - ){
nearX = mssHeight - ;
}
if(nearY >= mssWidth - ){
nearY = mssWidth - ;
}
if(nearX < mssHeight && nearY < mssWidth){
int i = ;
while(i < mssBandCount){
float tmp = mssData[nearX*mssWidth*mssBandCount + i*mssWidth + nearY];
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(tmp);
i++;
}
}
} if(MethodType == ){ // 双线性
float dataX = curX - (int)curX;
float dataY = curY - (int)curY;
int preX = (int)curX;
int preY = (int)curY;
int postX = (int)curX + ;
int postY = (int)curY + ;
if(postX >= mssHeight - ){
postX = mssHeight - ;
}
if(postY >= mssWidth - ){
postY = mssWidth - ;
} float Wx1 = - dataX;
float Wx2 = dataX;
float Wy1 = - dataY;
float Wy2 = dataY;
// 双线性差值核心代码
int i = ;
while(i < mssBandCount){
float pMssValue[] = {,,,};
pMssValue[] = mssData[preX*mssWidth*mssBandCount + i*mssWidth + preY];
pMssValue[] = mssData[preX*mssWidth*mssBandCount + i*mssWidth + postY];
pMssValue[] = mssData[postX*mssWidth*mssBandCount + i*mssWidth + preY];
pMssValue[] = mssData[postX*mssWidth*mssBandCount + i*mssWidth + postY];
float tmp = Wy1*(Wx1*pMssValue[] + Wx2*pMssValue[]) + Wy2*(Wx1*pMssValue[] + Wx2*pMssValue[]);
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(tmp);
i++;
}
} if(MethodType == ){ // 双三次卷积
float dataX = curX - (int)curX;
float dataY = curY - (int)curY;
int preX1 = (int)curX - ;
int preX2 = (int)curX;
int postX1 = (int)curX + ;
int postX2 = (int)curX + ;
int preY1 = (int)curY - ;
int preY2 = (int)curY;
int postY1 = (int)curY + ;
int postY2 = (int)curY + ;
if(preX1 < ) preX1 = ;
if(preY1 < ) preY1 = ;
if(postX1 > mssHeight - ) postX1 = mssHeight - ;
if(postX2 > mssHeight - ) postX2 = mssHeight - ;
if(postY1 > mssWidth - ) postY1 = mssWidth - ;
if(postY2 > mssWidth - ) postY2 = mssWidth - ; float Wx1 = -1.0f*dataX + *dataX*dataX - dataX*dataX*dataX;
float Wx2 = - *dataX*dataX + dataX*dataX*dataX;
float Wx3 = dataX + dataX*dataX - dataX*dataX*dataX;
float Wx4 = -1.0f*dataX*dataX + dataX*dataX*dataX;
float Wy1 = -1.0f*dataY + *dataY*dataY - dataY*dataY*dataY;
float Wy2 = - *dataY*dataY + dataY*dataY*dataY;
float Wy3 = dataY + dataY*dataY - dataY*dataY*dataY;
float Wy4 = -1.0f*dataY*dataY + dataY*dataY*dataY; int i = ;
while(i < mssBandCount){
float *pMssValue = new float[];
memset(pMssValue,,sizeof(float)*);
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + postY2]; pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + postY2]; pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + postY2]; pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + postY2]; float tmp = Wy1*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[])+
Wy2*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[])+
Wy3*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[])+
Wy4*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[]);
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = T(tmp);
delete []pMssValue;pMssValue = NULL;
i++;
}
}
}
}
} int ReSampleCPUApp(const char *mssfileName,
const char *panfileName,
const char *resamplefileName,
int MethodType = )
{
GDALAllRegister();
GDALDataset *poPANDS = (GDALDataset*)GDALOpen(panfileName,GA_ReadOnly);
GDALDataset *poMSSDS = (GDALDataset*)GDALOpen(mssfileName,GA_ReadOnly);
if(!poPANDS || !poMSSDS)
return -; //MSS info
int mssBandCount = poMSSDS->GetRasterCount();
int mssWidth = poMSSDS->GetRasterXSize();
int mssHeight = poMSSDS->GetRasterYSize();
double adfMssGeoTransform[] = {};
poMSSDS->GetGeoTransform(adfMssGeoTransform);
GDALDataType mssDT = poMSSDS->GetRasterBand()->GetRasterDataType(); int bSrcHasNoData;
double dfSrcNoDataValue = ;
dfSrcNoDataValue = GDALGetRasterNoDataValue(poMSSDS->GetRasterBand(),&bSrcHasNoData);
if(!bSrcHasNoData) dfSrcNoDataValue = 0.0;
//DT = mssDT; int *pBandMap= new int[mssBandCount];
for(int i = ;i<mssBandCount;i++){
pBandMap[i] = i+;
} // PAN Info
int panBandCount = poPANDS->GetRasterCount();
int panWidth = poPANDS->GetRasterXSize();
int panHeidht = poPANDS->GetRasterYSize();
double adfPanGeoTransform[] = {};
poPANDS->GetGeoTransform(adfPanGeoTransform);
GDALDataType panDT = poPANDS->GetRasterBand()->GetRasterDataType(); // 创建新数据集=======投影信息
double adfResampleGeoTransform[] = {};
adfResampleGeoTransform[] = adfPanGeoTransform[];
adfResampleGeoTransform[] = adfPanGeoTransform[];
adfResampleGeoTransform[] = adfPanGeoTransform[];
adfResampleGeoTransform[] = adfPanGeoTransform[];
if(adfMssGeoTransform[] >= adfPanGeoTransform[]){
adfResampleGeoTransform[] = adfMssGeoTransform[];
}else{
adfResampleGeoTransform[] = adfPanGeoTransform[];
}
if(adfMssGeoTransform[] > adfPanGeoTransform[]){
adfResampleGeoTransform[] = adfPanGeoTransform[];
}else{
adfResampleGeoTransform[] = adfMssGeoTransform[];
} // 创建新数据集=======影像大小
double panEndX = adfPanGeoTransform[] + panWidth*adfPanGeoTransform[] +
panHeidht*adfPanGeoTransform[];
double panEndY = adfPanGeoTransform[] + panHeidht*adfPanGeoTransform[] +
panHeidht*adfPanGeoTransform[]; double mssEndX = adfMssGeoTransform[] +mssWidth*adfMssGeoTransform[] +
mssHeight*adfMssGeoTransform[];
double mssEndY = adfMssGeoTransform[] + mssWidth*adfMssGeoTransform[] +
mssHeight*adfMssGeoTransform[];
double resampleEndXY[] = {};
if(panEndX > mssEndX)
resampleEndXY[] = mssEndX;
else
resampleEndXY[] = panEndX;
if(panEndY >= mssEndY)
resampleEndXY[] = panEndY;
else
resampleEndXY[] = mssEndY; // 创建新数据集=======MSS AND PAN 有效长宽
int resampleWidth = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfResampleGeoTransform[] + 0.5);
int resampleHeight = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfResampleGeoTransform[] + 0.5);
int mssEffectiveWidth = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfMssGeoTransform[] + 0.5);
int mssEffectiveHeight = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfMssGeoTransform[] + 0.5);
int panEffectiveWidth = resampleWidth;
int panEffectiveHeight = resampleHeight; // 创建新数据集=======位置增益大小
int mssGainX = static_cast<int>((adfResampleGeoTransform[] - adfMssGeoTransform[])/adfMssGeoTransform[] + 0.5);
int mssGainY = static_cast<int>((adfResampleGeoTransform[] - adfMssGeoTransform[])/adfMssGeoTransform[] + 0.5);
int panGainX = static_cast<int>((adfResampleGeoTransform[] - adfPanGeoTransform[])/adfPanGeoTransform[] + 0.5);
int panGainY = static_cast<int>((adfResampleGeoTransform[] - adfPanGeoTransform[])/adfPanGeoTransform[] + 0.5); // 创建新数据集=======创建文件
GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF");
if(!poOutDriver){
return -;
}
GDALDataset *poOutDS = poOutDriver->Create(resamplefileName,resampleWidth,
resampleHeight,mssBandCount,mssDT,NULL);
poOutDS->SetGeoTransform(adfResampleGeoTransform);
poOutDS->SetProjection(poPANDS->GetProjectionRef()); // 重采样核心代码============图像分块
int iNumRow = ;
if(iNumRow > resampleHeight){
iNumRow = ;
}
int loopNum = (resampleHeight + iNumRow - )/iNumRow; //分块数
int nLineSpace,nPixSpace,nBandSpace;
nLineSpace = sizeof(float)*mssEffectiveWidth*mssBandCount;
nPixSpace = ;
nBandSpace = sizeof(float)*mssEffectiveWidth; // 重采样采样比例
float radioX = float(adfPanGeoTransform[]/adfMssGeoTransform[]);
float radioY = float(adfPanGeoTransform[]/adfMssGeoTransform[]); int mssCurPosX = mssGainX;
int mssCurPosY = mssGainY;
int mssCurWidth = ;
int mssCurHeight = ; // 重采样核心代码============
for(int i = ;i<loopNum;i++){
int tmpRowNum = iNumRow;
int startR = i*iNumRow;
int endR = startR + iNumRow - ;
if(endR>resampleHeight -){
tmpRowNum = resampleHeight - startR;
//endR = startR + tmpRowNum - 1;
}
//计算读取的MSS影像区域大小
int mssCurWidth = mssEffectiveWidth;
int mssCurHeight = ;
if(MethodType == )
mssCurHeight = int(tmpRowNum*radioY);
else if(MethodType == )
mssCurHeight = int(tmpRowNum*radioY)+;
else if(MethodType == )
mssCurHeight = int(tmpRowNum*radioY)+; if(mssCurHeight + mssCurPosY > mssHeight - ){
mssCurHeight = mssHeight - mssCurPosY;
} //创建数据
/*float *resampleBuf = (float *)malloc(sizeof(cl_float)*tmpRowNum*resampleWidth*trueBandCount);*/
float *mssBuf = (float *)malloc(sizeof(cl_float)*mssCurHeight*mssCurWidth*mssBandCount);
//memset(resampleBuf,0,sizeof(float)*tmpRowNum*resampleWidth*trueBandCount);
memset(mssBuf,,sizeof(float)*mssCurHeight*mssCurWidth*mssBandCount); // 读取数据
poMSSDS->RasterIO(GF_Read,mssCurPosX,mssCurPosY,mssCurWidth,mssCurHeight,
mssBuf,mssCurWidth,mssCurHeight,GDT_Float32,mssBandCount,NULL,nPixSpace,
nLineSpace,nBandSpace); if(MethodType == )
mssCurPosY += mssCurHeight;
else if(MethodType == )
mssCurPosY += mssCurHeight - ;
else if(MethodType == )
mssCurPosY += mssCurHeight - ; // 数据格式转换
long sz = tmpRowNum*resampleWidth*mssBandCount;
void *resampleBuf = NULL;
switch(mssDT){
case GDT_Byte:
resampleBuf = new unsigned char[sz];
ReSampleCPUKernel<unsigned char>(mssBuf,(unsigned char*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(unsigned char),
resampleWidth*sizeof(unsigned char));
break;
case GDT_UInt16:
resampleBuf = new unsigned short int[sz];
ReSampleCPUKernel<unsigned short int>(mssBuf,(unsigned short int*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(unsigned short int),
resampleWidth*sizeof(unsigned short int));
break;
case GDT_Int16:
resampleBuf = new short int[sz];
ReSampleCPUKernel<short int>(mssBuf,(short int*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(short int),
resampleWidth*sizeof(short int));
break;
case GDT_UInt32:
resampleBuf = new unsigned int[sz];
ReSampleCPUKernel<unsigned int>(mssBuf,(unsigned int*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(unsigned int),
resampleWidth*sizeof(unsigned int));
break;
case GDT_Int32:
resampleBuf = new int[sz];
ReSampleCPUKernel<int>(mssBuf,(int*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(int),
resampleWidth*sizeof(int));
break;
case GDT_Float32:
resampleBuf = new float[sz];
ReSampleCPUKernel<float>(mssBuf,(float*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(float),
resampleWidth*sizeof(float));
break;
case GDT_Float64:
resampleBuf = new double[sz];
ReSampleCPUKernel<double>(mssBuf,(double*)resampleBuf,mssCurWidth,mssCurHeight,mssBandCount,mssGainX,mssGainY,
resampleWidth,tmpRowNum,radioX,radioY,dfSrcNoDataValue,MethodType);
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,resampleBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(double),
resampleWidth*sizeof(double));
break;
}
delete []mssBuf;
delete []resampleBuf;
std::cout<<i<<std::endl;
}
delete []pBandMap;pBandMap = NULL;
GDALClose((GDALDatasetH)poPANDS);
GDALClose((GDALDatasetH)poMSSDS);
GDALClose((GDALDatasetH)poOutDS);
return ;
} #endif
GPU版本：
#ifndef RESAMPLEOPENCL_H
#define RESAMPLEOPENCL_H #include <CL/cl.h>
#include <gdal_alg_priv.h>
#include <gdal_priv.h> #pragma comment(lib,"OpenCL.lib") /*
@ 功能描述
读取源程序，将文本源程序读到内核中
*/
char* LoadProgSource(const char* cFilename, const char* cPreamble, size_t* szFinalLength)
{
FILE* pFileStream = NULL;
size_t szSourceLength; // open the OpenCL source code file
pFileStream = fopen(cFilename, "rb");
if(pFileStream == )
{
return NULL;
} size_t szPreambleLength = strlen(cPreamble); // get the length of the source code
fseek(pFileStream, , SEEK_END);
szSourceLength = ftell(pFileStream);
fseek(pFileStream, , SEEK_SET); // allocate a buffer for the source code string and read it in
char* cSourceString = (char *)malloc(szSourceLength + szPreambleLength + );
memcpy(cSourceString, cPreamble, szPreambleLength);
if (fread((cSourceString) + szPreambleLength, szSourceLength, , pFileStream) != )
{
fclose(pFileStream);
free(cSourceString);
return ;
} // close the file and return the total length of the combined (preamble + source) string
fclose(pFileStream);
if(szFinalLength != )
{
*szFinalLength = szSourceLength + szPreambleLength;
}
cSourceString[szSourceLength + szPreambleLength] = '\0'; return cSourceString;
} template<typename T>
bool DataTypeTrans(const float *pSrcBuf,T *pDesBuf,long size)
{
if(pSrcBuf == NULL){
return false;
}
while(size--){
pDesBuf[size] = T(pSrcBuf[size]);
}
return true;
} int ReSampleOpenCLApp(const char *mssfileName,
const char *panfileName,
const char *resamplefileName,
int MethodType = )
{
GDALAllRegister();
GDALDataset *poPANDS = (GDALDataset*)GDALOpen(panfileName,GA_ReadOnly);
GDALDataset *poMSSDS = (GDALDataset*)GDALOpen(mssfileName,GA_ReadOnly);
if(!poPANDS || !poMSSDS)
return -; //MSS info
int mssBandCount = poMSSDS->GetRasterCount();
int mssWidth = poMSSDS->GetRasterXSize();
int mssHeight = poMSSDS->GetRasterYSize();
double adfMssGeoTransform[] = {};
poMSSDS->GetGeoTransform(adfMssGeoTransform);
GDALDataType mssDT = poMSSDS->GetRasterBand()->GetRasterDataType(); int bSrcHasNoData;
float dfSrcNoDataValue = ;
dfSrcNoDataValue = (float)GDALGetRasterNoDataValue(poMSSDS->GetRasterBand(),&bSrcHasNoData);
if(!bSrcHasNoData) dfSrcNoDataValue = 0.0; // PAN Info
int panBandCount = poPANDS->GetRasterCount();
int panWidth = poPANDS->GetRasterXSize();
int panHeidht = poPANDS->GetRasterYSize();
double adfPanGeoTransform[] = {};
poPANDS->GetGeoTransform(adfPanGeoTransform);
GDALDataType panDT = poPANDS->GetRasterBand()->GetRasterDataType(); // 创建新数据集=======投影信息
double adfResampleGeoTransform[] = {};
adfResampleGeoTransform[] = adfPanGeoTransform[];
adfResampleGeoTransform[] = adfPanGeoTransform[];
adfResampleGeoTransform[] = adfPanGeoTransform[];
adfResampleGeoTransform[] = adfPanGeoTransform[];
if(adfMssGeoTransform[] >= adfPanGeoTransform[]){
adfResampleGeoTransform[] = adfMssGeoTransform[];
}else{
adfResampleGeoTransform[] = adfPanGeoTransform[];
}
if(adfMssGeoTransform[] > adfPanGeoTransform[]){
adfResampleGeoTransform[] = adfPanGeoTransform[];
}else{
adfResampleGeoTransform[] = adfMssGeoTransform[];
} // 创建新数据集=======影像大小
double panEndX = adfPanGeoTransform[] + panWidth*adfPanGeoTransform[] +
panHeidht*adfPanGeoTransform[];
double panEndY = adfPanGeoTransform[] + panHeidht*adfPanGeoTransform[] +
panHeidht*adfPanGeoTransform[]; double mssEndX = adfMssGeoTransform[] +mssWidth*adfMssGeoTransform[] +
mssHeight*adfMssGeoTransform[];
double mssEndY = adfMssGeoTransform[] + mssWidth*adfMssGeoTransform[] +
mssHeight*adfMssGeoTransform[];
double resampleEndXY[] = {};
if(panEndX > mssEndX)
resampleEndXY[] = mssEndX;
else
resampleEndXY[] = panEndX;
if(panEndY >= mssEndY)
resampleEndXY[] = panEndY;
else
resampleEndXY[] = mssEndY; // 创建新数据集=======MSS AND PAN 有效长宽
int resampleWidth = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfResampleGeoTransform[] + 0.5);
int resampleHeight = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfResampleGeoTransform[] + 0.5);
int mssEffectiveWidth = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfMssGeoTransform[] + 0.5);
int mssEffectiveHeight = static_cast<int>((resampleEndXY[] - adfResampleGeoTransform[])/adfMssGeoTransform[] + 0.5);
int panEffectiveWidth = resampleWidth;
int panEffectiveHeight = resampleHeight; // 创建新数据集=======位置增益大小
int mssGainX = static_cast<int>((adfResampleGeoTransform[] - adfMssGeoTransform[])/adfMssGeoTransform[] + 0.5);
int mssGainY = static_cast<int>((adfResampleGeoTransform[] - adfMssGeoTransform[])/adfMssGeoTransform[] + 0.5);
int panGainX = static_cast<int>((adfResampleGeoTransform[] - adfPanGeoTransform[])/adfPanGeoTransform[] + 0.5);
int panGainY = static_cast<int>((adfResampleGeoTransform[] - adfPanGeoTransform[])/adfPanGeoTransform[] + 0.5); // 创建新数据集=======创建文件
GDALDriver *poOutDriver = (GDALDriver*)GDALGetDriverByName("GTIFF");
if(!poOutDriver){
return -;
}
GDALDataset *poOutDS = poOutDriver->Create(resamplefileName,resampleWidth,
resampleHeight,mssBandCount,mssDT,NULL);
//GDALDataset *poOutDS = poOutDriver->Create(resamplefileName,resampleWidth,
// resampleHeight,mssBandCount,GDT_Float32,NULL);
poOutDS->SetGeoTransform(adfResampleGeoTransform);
poOutDS->SetProjection(poPANDS->GetProjectionRef()); int pBandMap[] = {,,,};
// 重采样核心代码============图像分块
int iNumRow = ;
if(iNumRow > resampleHeight){
iNumRow = ;
}
int loopNum = (resampleHeight + iNumRow - )/iNumRow; //分块数
int nLineSpace,nPixSpace,nBandSpace;
nLineSpace = sizeof(float)*mssEffectiveWidth*mssBandCount;
nPixSpace = ;
nBandSpace = sizeof(float)*mssEffectiveWidth; // 重采样采样比例
float radioX = adfPanGeoTransform[]/adfMssGeoTransform[];
float radioY = adfPanGeoTransform[]/adfMssGeoTransform[]; int mssCurPosX = mssGainX;
int mssCurPosY = mssGainY;
int mssCurWidth = ;
int mssCurHeight = ; // 重采样核心代码============
// OpenCL部分 =============== 1 创建平台
cl_uint num_platforms;
cl_int ret = clGetPlatformIDs(,NULL,&num_platforms);
if(ret != CL_SUCCESS || num_platforms < ){
printf("clGetPlatformIDs Error\n");
return -;
}
cl_platform_id platform_id = NULL;
ret = clGetPlatformIDs(,&platform_id,NULL);
if(ret != CL_SUCCESS){
printf("clGetPlatformIDs Error2\n");
return -;
} // OpenCL部分 =============== 2 获得设备
cl_uint num_devices;
ret = clGetDeviceIDs(platform_id,CL_DEVICE_TYPE_GPU,,NULL,
&num_devices);
if(ret != CL_SUCCESS || num_devices < ){
printf("clGetDeviceIDs Error\n");
return -;
}
cl_device_id device_id;
ret = clGetDeviceIDs(platform_id,CL_DEVICE_TYPE_GPU,,&device_id,NULL);
if(ret != CL_SUCCESS){
printf("clGetDeviceIDs Error2\n");
return -;
} // OpenCL部分 =============== 3 创建Context
cl_context_properties props[] = {CL_CONTEXT_PLATFORM,
(cl_context_properties)platform_id,};
cl_context context = NULL;
context = clCreateContext(props,,&device_id,NULL,NULL,&ret);
if(ret != CL_SUCCESS || context == NULL){
printf("clCreateContext Error\n");
return -;
} // OpenCL部分 =============== 4 创建Command Queue
cl_command_queue command_queue = NULL;
command_queue = clCreateCommandQueue(context,device_id,,&ret);
if(ret != CL_SUCCESS || command_queue == NULL){
printf("clCreateCommandQueue Error\n");
return -;
} // OpenCL部分 =============== 6 创建编译Program
const char *strfile = "D:\\PIE3\\src\\Test\\TextOpecCLResample\\TextOpecCLResample\\ReSampleKernel.txt";
size_t lenSource = ;
char *kernelSource = LoadProgSource(strfile,"",&lenSource);
cl_program *programs = (cl_program *)malloc(loopNum*sizeof(cl_program));
memset(programs,,sizeof(cl_program)*loopNum); cl_kernel *kernels = (cl_kernel*)malloc(loopNum*sizeof(cl_kernel));
memset(kernels,,sizeof(cl_kernel)*loopNum); for(int i = ;i<loopNum;i++){
int tmpRowNum = iNumRow;
int startR = i*iNumRow;
int endR = startR + iNumRow - ;
if(endR>resampleHeight -){
tmpRowNum = resampleHeight - startR;
//endR = startR + tmpRowNum - 1;
}
//计算读取的MSS影像区域大小
int mssCurWidth = mssEffectiveWidth;
int mssCurHeight = ;
if(MethodType == )
mssCurHeight = int(tmpRowNum*radioY);
else if(MethodType == )
mssCurHeight = int(tmpRowNum*radioY)+;
else if(MethodType == )
mssCurHeight = int(tmpRowNum*radioY)+; if(mssCurHeight + mssCurPosY > mssHeight - ){
mssCurHeight = mssHeight - mssCurPosY;
} //创建数据
float *resampleBuf = (float *)malloc(sizeof(cl_float)*tmpRowNum*resampleWidth*mssBandCount);
float *mssBuf = (float *)malloc(sizeof(cl_float)*mssCurHeight*mssCurWidth*mssBandCount);
memset(resampleBuf,,sizeof(cl_float)*tmpRowNum*resampleWidth*mssBandCount);
memset(mssBuf,,sizeof(cl_float)*mssCurHeight*mssCurWidth*mssBandCount); // 读取数据
poMSSDS->RasterIO(GF_Read,mssCurPosX,mssCurPosY,mssCurWidth,mssCurHeight,
mssBuf,mssCurWidth,mssCurHeight,GDT_Float32,mssBandCount,pBandMap,nPixSpace,
nLineSpace,nBandSpace); if(MethodType == )
mssCurPosY += mssCurHeight;
else if(MethodType == )
mssCurPosY += mssCurHeight - ;
else if(MethodType == )
mssCurPosY += mssCurHeight - ; // OpenCL部分 =============== 5 创建Memory Object
cl_mem mem_mss = NULL;
mem_mss = clCreateBuffer(context,CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR,
sizeof(cl_float)*mssCurHeight*mssCurWidth*mssBandCount,mssBuf,&ret);
if(ret != CL_SUCCESS || NULL == mem_mss){
printf("clCreateBuffer Error\n");
return -;
} cl_mem mem_resample = NULL;
mem_resample = clCreateBuffer(context,CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR,
sizeof(cl_float)*resampleWidth*tmpRowNum*mssBandCount,resampleBuf,&ret);
if(ret != CL_SUCCESS || NULL == mem_resample){
printf("clCreateBuffer Error\n");
return -;
} // OpenCL部分 =============== 6 创建编译Program
//const char *strfile = "D:\\PIE3\\src\\Test\\TextOpecCLResample\\TextOpecCLResample\\ReSampleKernel.txt";
//size_t lenSource = 0;
//char *kernelSource = LoadProgSource(strfile,"",&lenSource);
//cl_program program = NULL;
programs[i] = clCreateProgramWithSource(context,,(const char**)&kernelSource,
NULL,&ret);
if(ret != CL_SUCCESS || NULL == programs[i]){
printf("clCreateProgramWithSource Error\n");
return -;
}
ret = clBuildProgram(programs[i],,&device_id,NULL,NULL,NULL);
if(ret != CL_SUCCESS){
char* build_log;
size_t log_size;
//查询日志的大小
clGetProgramBuildInfo(programs[i], device_id, CL_PROGRAM_BUILD_LOG, , NULL, &log_size);
build_log = new char[log_size+];
//获得编译日志信息
ret = clGetProgramBuildInfo(programs[i], device_id, CL_PROGRAM_BUILD_LOG, log_size, build_log, NULL);
build_log[log_size] = '\0';
printf("%s\n",build_log);
printf("编译失败!");
delete []build_log;
return -;
} // OpenCL部分 =============== 7 创建Kernel
//cl_kernel kernel = NULL;
kernels[i] = clCreateKernel(programs[i],"ReSampleKernel",&ret);
if(ret != CL_SUCCESS || NULL == kernels[i]){
printf("clCreateProgramWithSource Error\n");
return -;
} // OpenCL部分 =============== 8 设置Kernel参数
ret = clSetKernelArg(kernels[i],,sizeof(cl_mem),&mem_mss);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_mem),&mem_resample);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&mssCurWidth);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&mssCurHeight);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&mssBandCount);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&mssGainX);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&mssGainY);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&resampleWidth);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&tmpRowNum);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_float),&radioX);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_float),&radioY);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_float),&dfSrcNoDataValue);
ret |= clSetKernelArg(kernels[i],,sizeof(cl_int),&MethodType);
if(ret != CL_SUCCESS){
printf("clSetKernelArg Error\n");
return -;
} // OpenCL部分 =============== 9 设置Group Size
cl_uint work_dim = ;
size_t global_work_size[] = {resampleWidth,tmpRowNum};
size_t *local_work_size = NULL; // OpenCL部分 =============== 10 执行内核
ret = clEnqueueNDRangeKernel(command_queue,kernels[i],work_dim,NULL,global_work_size,
local_work_size,,NULL,NULL);
ret |= clFinish(command_queue);
if(ret != CL_SUCCESS){
printf("clEnqueueNDRangeKernel Error\n");
return -;
} // OpenCL部分 =============== 11 读取结果 resampleBuf = (float*)clEnqueueMapBuffer(command_queue,mem_resample,CL_TRUE,CL_MAP_READ | CL_MAP_WRITE,
,sizeof(cl_float)*tmpRowNum*resampleWidth*mssBandCount,,NULL,NULL,&ret);
//ret = clEnqueueReadBuffer(command_queue,mem_resample,CL_TRUE,0,
// sizeof(cl_float)*tmpRowNum*resampleWidth*mssBandCount,(void*)resampleBuf,0,NULL,NULL);
if(ret != CL_SUCCESS){
printf("clEnqueueMapBuffer Error\n");
return -;
} // 数据格式转换
long sz = tmpRowNum*resampleWidth*mssBandCount;
void *pBuf = NULL;
CPLErr err;
switch(mssDT){
case GDT_Byte:
pBuf = new unsigned char[sz];
if(!DataTypeTrans<unsigned char>(resampleBuf,(unsigned char*)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(unsigned char),
resampleWidth*sizeof(unsigned char));
break;
case GDT_UInt16:
pBuf = new unsigned short int[sz];
if(!DataTypeTrans<unsigned short int>(resampleBuf,(unsigned short int*)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
err = poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(unsigned short int),
resampleWidth*sizeof(unsigned short int));
break;
case GDT_Int16:
pBuf = new short int[sz];
if(!DataTypeTrans<short int>(resampleBuf,(short int*)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(short int),
resampleWidth*sizeof(short int));
break;
case GDT_UInt32:
pBuf = new unsigned int[sz];
if(!DataTypeTrans<unsigned int>(resampleBuf,(unsigned int*)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(unsigned int),
resampleWidth*sizeof(unsigned int));
break;
case GDT_Int32:
pBuf = new int[sz];
if(!DataTypeTrans<int>(resampleBuf,(int*)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(int),
resampleWidth*sizeof(int));
break;
case GDT_Float32:
pBuf = new float[sz];
if(!DataTypeTrans<float>(resampleBuf,(float *)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(float),
resampleWidth*sizeof(float));
break;
case GDT_Float64:
pBuf = new double[sz];
if(!DataTypeTrans<double>(resampleBuf,(double *)pBuf,sz))
{
printf("DataTypeTrans Error\n");
return -;
}
poOutDS->RasterIO(GF_Write,,startR,resampleWidth,tmpRowNum,pBuf,
resampleWidth,tmpRowNum,mssDT,mssBandCount,NULL,nPixSpace,resampleWidth*mssBandCount*sizeof(double),
resampleWidth*sizeof(double));
break;
}
delete []pBuf;pBuf = NULL;
free(mssBuf);
free(resampleBuf); // OpenCL部分 =============== 12 释放资源
if(NULL != mem_mss) clReleaseMemObject(mem_mss);
if(NULL != mem_resample) clReleaseMemObject(mem_resample);
std::cout<<i<<std::endl;
}
// OpenCL部分 =============== 12 释放资源
int i = ;
while(i < loopNum){
if(NULL != kernels[i]) clReleaseKernel(kernels[i]);
if(NULL != programs[i]) clReleaseProgram(programs[i]);
i++;
} if(NULL != command_queue) clReleaseCommandQueue(command_queue);
if(NULL != context) clReleaseContext(context);
GDALClose((GDALDatasetH)poPANDS);
GDALClose((GDALDatasetH)poMSSDS);
GDALClose((GDALDatasetH)poOutDS);
return ;
} #endif
GPU核函数代码如下：
#pragma OPENCL EXTENSION cl_amd_printf:enable __kernel void ReSampleKernel(__global const float *mssData,
__global float *resampleData,
int mssWidth,
int mssHeight,
int mssBandCount,
int mssOffsetX,
int mssOffsetY,
int panWidth,
int panHeight,
float radioX,
float radioY,
float dfDstNoDataValue,
int MethodType)
{
int idx = get_global_id(); // 采样行
int idy = get_global_id(); // 采样列
float eps = 0.00001f;
if(idx < panHeight && idy < panWidth){
// 找到对应的MSS像素位置
float curX = (float)idx * radioX;
float curY = (float)idy * radioY;
int tmpP = (int)curX*mssWidth*mssBandCount + (int)curY;
if(mssData[tmpP] == dfDstNoDataValue)
{
int i = ;
while(i < mssBandCount){
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] = dfDstNoDataValue;
i++;
}
return;
}
if(MethodType == ){ // 最近邻
int nearX = (int)(curX + 0.5)>(int)curX?(int)(curX + ):(int)curX;
int nearY = (int)(curY + 0.5)>(int)curY?(int)(curY + ):(int)curY;
if(nearX >= mssHeight - ){
nearX = mssHeight - ;
}
if(nearY >= mssWidth - ){
nearY = mssWidth - ;
}
if(nearX < mssHeight && nearY < mssWidth){
int i = ;
while(i < mssBandCount){
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] =
mssData[nearX*mssWidth*mssBandCount + i*mssWidth + nearY];
i++;
}
}
}
if(MethodType == ){ // 双线性
float dataX = curX - (int)curX;
float dataY = curY - (int)curY;
if(dataX < eps){
dataX = 0.00001;
}
if(dataY < eps){
dataY = 0.00001;
}
int preX = (int)curX;
int preY = (int)curY;
int postX = (int)curX + ;
int postY = (int)curY + ;
if(postX >= mssHeight - ){
postX = mssHeight - ;
}
if(postY >= mssWidth - ){
postY = mssWidth - ;
} float Wx1 = - dataX;
float Wx2 = dataX;
float Wy1 = - dataY;
float Wy2 = dataY;
// 双线性差值核心代码
int i = ;
while(i < mssBandCount){
float pMssValue[] = {,,,};
pMssValue[] = mssData[preX*mssWidth*mssBandCount + i*mssWidth + preY];
pMssValue[] = mssData[preX*mssWidth*mssBandCount + i*mssWidth + postY];
pMssValue[] = mssData[postX*mssWidth*mssBandCount + i*mssWidth + preY];
pMssValue[] = mssData[postX*mssWidth*mssBandCount + i*mssWidth + postY];
resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] =
Wy1*(Wx1*pMssValue[] + Wx2*pMssValue[]) + Wy2*(Wx1*pMssValue[] + Wx2*pMssValue[]);
i++;
}
}
if(MethodType == ){ // 双三次卷积
float dataX = curX - (int)curX;
float dataY = curY - (int)curY;
//printf("dataX = %f dataY = %f\n",dataX,dataY);
int preX1 = (int)curX - ;
int preX2 = (int)curX;
int postX1 = (int)curX + ;
int postX2 = (int)curX + ;
int preY1 = (int)curY - ;
int preY2 = (int)curY;
int postY1 = (int)curY + ;
int postY2 = (int)curY + ;
if(preX1 < ) preX1 = ;
if(preY1 < ) preY1 = ;
if(postX1 > mssHeight - ) postX1 = mssHeight - ;
if(postX2 > mssHeight - ) postX2 = mssHeight - ;
if(postY1 > mssWidth - ) postY1 = mssWidth - ;
if(postY2 > mssWidth - ) postY2 = mssWidth - ; float Wx1 = -1.0f*dataX + *dataX*dataX - dataX*dataX*dataX;
float Wx2 = - *dataX*dataX + dataX*dataX*dataX;
float Wx3 = dataX + dataX*dataX - dataX*dataX*dataX;
float Wx4 = -1.0f*dataX*dataX + dataX*dataX*dataX;
float Wy1 = -1.0f*dataY + *dataY*dataY - dataY*dataY*dataY;
float Wy2 = - *dataY*dataY + dataY*dataY*dataY;
float Wy3 = dataY + dataY*dataY - dataY*dataY*dataY;
float Wy4 = -1.0f*dataY*dataY + dataY*dataY*dataY; //printf("preX1 = %d\n",preX1);
int i = ;
while(i < mssBandCount){
float pMssValue[] = {,,,,,,,,,,,,,,,};
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[preX1*mssWidth*mssBandCount + i*mssWidth + postY2]; pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[preX2*mssWidth*mssBandCount + i*mssWidth + postY2]; pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[postX1*mssWidth*mssBandCount + i*mssWidth + postY2]; pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + preY1];
pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + preY2];
pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + postY1];
pMssValue[] = mssData[postX2*mssWidth*mssBandCount + i*mssWidth + postY2]; resampleData[idx*panWidth*mssBandCount+i*panWidth + idy] =
Wy1*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[])+
Wy2*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[])+
Wy3*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[])+
Wy4*(Wx1*pMssValue[] + Wx2*pMssValue[] + Wx3*pMssValue[] + Wx4*pMssValue[]);
i++;
}
}
}
}
以上代码应该可以直接使用，欢迎大家一起交流探讨。
另外，我对GDAL、CPU和GPU版本的重采样算法效率进行了一下对比，GPU在三次卷积重采样算法上要明显的比CPU版本效率高很多。具体结果如下：
