在編寫重采樣圖像時，可以使用GDAL來讀寫圖像，然后自己編寫重采樣算法（最鄰近像元法，雙線性內(nèi)插法，三次立方卷積法等）【關(guān)于這采樣算法有時間我會單獨(dú)寫一篇文章來說明原理的】將計(jì)算的結(jié)果寫入圖像中來實(shí)現(xiàn)。

在GDAL的算法中，已經(jīng)提供了五種重采樣算法，其定義如下（位置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的源代碼發(fā)現(xiàn)，warp的功能非常強(qiáng)大，可以用來做投影轉(zhuǎn)換，重投影，投影定義，重采樣，鑲嵌，幾何精校正和影像配準(zhǔn)等。一句話，很好很強(qiáng)大。下面就看看其中的一點(diǎn)點(diǎn)皮毛，使用warp來編寫一個重采樣的接口，代碼如下：

/**
* 重采樣函數(shù)(GDAL)
* @param pszSrcFile        輸入文件的路徑
* @param pszOutFile        寫入的結(jié)果圖像的路徑
* @param fResX             X轉(zhuǎn)換采樣比，默認(rèn)大小為1.0，大于1圖像變大，小于1表示圖像縮小
* @param fResY             Y轉(zhuǎn)換采樣比，默認(rèn)大小為1.0
* @param nResampleMode     采樣模式，有五種，具體參見GDALResampleAlg定義，默認(rèn)為雙線性內(nèi)插
* @param pExtent           采樣范圍，為NULL表示計(jì)算全圖
* @param pBandIndex        指定的采樣波段序號，為NULL表示采樣全部波段
* @param pBandCount        采樣的波段個數(shù)，同pBandIndex一同使用，表示采樣波段的個數(shù)
* @param pszFormat         寫入的結(jié)果圖像的格式
* @param pProgress         進(jìn)度條指針
* @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("正?在?執(zhí)?行?重?采?樣?...");
    }

    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("不?能?創(chuàng)?建?該?格?式?的?文?件?！?");

        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("指?定?的?波?段?序?號?超?過?圖?像?的?波?段?數(shù)?，?請?檢?查?輸?入?參?數(shù)?！?");

            GDALClose((GDALDatasetH) pDSrc);
            return RE_PARAMERROR;
        }
        
        iNewBandCount = *pBandCount;
    }

    LT_Envelope enExtent;
    enExtent.setToNull();

    if (pExtent == NULL)    //全?圖?計(jì)?算?
    {
        double dPrj[4] = {0};    //x1,x2,y1,y2
        ImageRowCol2Projection(dGeoTrans, 0, 0, dPrj[0], dPrj[2]);
        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[1] = dGeoTrans[1] / fResX;
    dGeoTrans[5] = dGeoTrans[5] / fResY;

    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("創(chuàng)?建?輸?出?文?件?失?敗?！?");

        GDALClose((GDALDatasetH) pDSrc);
        return RE_CREATEFILE;
    }

    pDDst->SetProjection(strWkt.c_str());
    pDDst->SetGeoTransform(dGeoTrans);

    GDALResampleAlg eResample = (GDALResampleAlg) nResampleMode;

    if(pProgress != NULL)
    {
        pProgress->SetProgressTip("正?在?執(zhí)?行?重?采?樣?...");
        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("轉(zhuǎn)?換?參?數(shù)?錯?誤?！?");

        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("轉(zhuǎn)?換?參?數(shù)?錯?誤?！?");

        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;
}

PS：在使用Windows Live Writer來寫博客，使用VSPaste插件粘貼代碼的時候，發(fā)現(xiàn)會在漢字后面加一個“？”號，我懶得修改了，大家就湊合看看吧！

總結(jié)

以上是生活随笔為你收集整理的如何使用GDAL重采样图像的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。