如何使用GDAL重采样图像

来源：互联网发布：mac摔了下屏幕不亮编辑：程序博客网时间：2024/04/30 10:08

在编写重采样图像时，可以使用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]);        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("创?建?输?出?文?件?失?败?！?");        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;}

PS：在使用Windows Live Writer来写博客，使用VSPaste插件粘贴代码的时候，发现会在汉字后面加一个“？”号，我懒得修改了，大家就凑合看看吧！