基于GDAL的一个通用的3×3模板函数

在进行遥感图像处理时，经常会用到很多的模板算子，比如平滑锐化等，拉普拉斯算子，索伯尔算子等等。其实这些算法都一样，用一个模板窗口在图像上移动，然后把计算的结果写入图像中。

在查看GDAL源代码的时候，有个gdaldem的工具，里面有一个类似3×3的模板函数，我改造了一下，可以支持任意的3×3的模板运算。

/**
* @brief 3*3模板计算处理函数
* @param hSrcBand        输入图像波段
* @param hDstBand        输出图像波段
* @param pfnAlg            算法回调函数
* @param pData            算法回调函数参数
* @param pfnProgress    进度条回调函数
* @param pProgressData    进度条回调函数参数
*/
CPLErr GDALGeneric3x3Processing  ( GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand,
                                  GDALGeneric3x3ProcessingAlg pfnAlg, void* pData,
                                  GDALProgressFunc pfnProgress, void * pProgressData)
{
    CPLErr eErr;
    float *pafThreeLineWin;  /* 输入图像三行数据存储空间 */
    float *pafOutputBuf;     /* 输出图像一行数据存储空间 */
    int i, j;

    int bSrcHasNoData, bDstHasNoData;
    float fSrcNoDataValue = 0.0, fDstNoDataValue = 0.0;

    int nXSize = GDALGetRasterBandXSize(hSrcBand);
    int nYSize = GDALGetRasterBandYSize(hSrcBand);

    if (pfnProgress == NULL)
        pfnProgress = GDALDummyProgress;

    //初始化进度条计数器
    if( !pfnProgress( 0.0, NULL, pProgressData ) )
    {
        CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
        return CE_Failure;
    }

    //分配内存空间
    pafOutputBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
    pafThreeLineWin  = (float *) CPLMalloc(3*sizeof(float)*(nXSize+1));

    fSrcNoDataValue = (float) GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData);
    fDstNoDataValue = (float) GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData);
    if (!bDstHasNoData)
        fDstNoDataValue = 0.0;

    // 移动一个3x3的窗口pafWindow遍历每个象元
    // (中心象元编号为#4)
    // 
    //      0 1 2
    //      3 4 5
    //      6 7 8

    // 预先加载前两行数据
    for ( i = 0; i < 2 && i < nYSize; i++)
    {
        GDALRasterIO( hSrcBand,    GF_Read, 0, i, nXSize, 1,
            pafThreeLineWin + i * nXSize, nXSize, 1,
            GDT_Float32, 0, 0);
    }

    //不包括边界
    for (j = 0; j < nXSize; j++)
        pafOutputBuf[j] = fDstNoDataValue;

    GDALRasterIO(hDstBand, GF_Write, 0, 0, nXSize, 1, 
        pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);

    if (nYSize > 1)
    {
        GDALRasterIO(hDstBand, GF_Write, 0, nYSize - 1, nXSize, 1,
            pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);
    }

    int nLine1Off = 0*nXSize;
    int nLine2Off = 1*nXSize;
    int nLine3Off = 2*nXSize;

    for ( i = 1; i < nYSize-1; i++)
    {
        // 读取第三行数据
        eErr = GDALRasterIO( hSrcBand, GF_Read, 0, i+1, nXSize, 1,
            pafThreeLineWin + nLine3Off, nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)
            goto end;

        //不包括边界数据
        pafOutputBuf[0] = fDstNoDataValue;
        if (nXSize > 1)
            pafOutputBuf[nXSize - 1] = fDstNoDataValue;

// 这句使用OpenMP来加速
#pragma omp parallel for
        for (j = 1; j < nXSize - 1; j++)
        {
            float afWin[9];
            afWin[0] = pafThreeLineWin[nLine1Off + j-1];
            afWin[1] = pafThreeLineWin[nLine1Off + j];
            afWin[2] = pafThreeLineWin[nLine1Off + j+1];
            afWin[3] = pafThreeLineWin[nLine2Off + j-1];
            afWin[4] = pafThreeLineWin[nLine2Off + j];
            afWin[5] = pafThreeLineWin[nLine2Off + j+1];
            afWin[6] = pafThreeLineWin[nLine3Off + j-1];
            afWin[7] = pafThreeLineWin[nLine3Off + j];
            afWin[8] = pafThreeLineWin[nLine3Off + j+1];

            if (bSrcHasNoData && (
                afWin[0] == fSrcNoDataValue ||
                afWin[1] == fSrcNoDataValue ||
                afWin[2] == fSrcNoDataValue ||
                afWin[3] == fSrcNoDataValue ||
                afWin[4] == fSrcNoDataValue ||
                afWin[5] == fSrcNoDataValue ||
                afWin[6] == fSrcNoDataValue ||
                afWin[7] == fSrcNoDataValue ||
                afWin[8] == fSrcNoDataValue))
            {
                // 如果9个数据中有一个是NoData则将当前点设置为NoData
                pafOutputBuf[j] = fDstNoDataValue;
            }
            else
            {
                //一个合格的3*3窗口
                pafOutputBuf[j] = pfnAlg(afWin, fDstNoDataValue, pData);
            }
        }

        //写入一行数据
        eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1,
            pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0);

        if (eErr != CE_None)
            goto end;

        if( !pfnProgress( 1.0 * (i+1) / nYSize, NULL, pProgressData ) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
            eErr = CE_Failure;
            goto end;
        }

        int nTemp = nLine1Off;
        nLine1Off = nLine2Off;
        nLine2Off = nLine3Off;
        nLine3Off = nTemp;
    }

    pfnProgress( 1.0, NULL, pProgressData );
    eErr = CE_None;

end:
    CPLFree(pafOutputBuf);
    CPLFree(pafThreeLineWin);

    return eErr;
}

  为了使这个函数能够通用，所以前面定义了一个函数指针来设置算法，以及一个void类型的参数用来设置算法的一些参数。关于这个回调函数的原型是下面的代码：

/**
* @brief 3*3模板算法回调函数
* @param pafWindow            3*3窗口值数组
* @param fDstNoDataValue    结果NoData值
* @param pData                算法回调函数参数
* @return 3*3窗口计算结果值
*/
typedef float (*GDALGeneric3x3ProcessingAlg) (float* pafWindow, float fDstNoDataValue, void* pData);

关于函数的用法和参数说明都在代码里面了。下面来两个例子。第一个，做个拉普拉斯算子吧。首先先大概说明一下拉普拉斯算子。Laplacian算子的公式是：Delta^2f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)，写成模板的样子如图1所示：

图1 Laplacian算子

好了，下面我们需要实现一个模板运算的算法函数和一个算法参数结构体，这个算法函数的类型就是上面的回调函数的样子，算法参数结构体用来存储模板算子等。下面先是算法的参数结构体代码，很简单，就是一个double[9]数组。

/**
* @brief 空间滤波算法数据结构体
*/
typedef struct
{
    /*! 模板算子 */
    double dTemplate[9];
} SpaceFilterAlgData;

接下来是，模板算法的实现函数，我们知道模板算法是用模板算子的每一个值乘以图像的像元值然后求和，得到的结果就是结果图像的像元值。

float SpaceFilterAlg (float* afWin, float fDstNoDataValue, void* pData)
{
    SpaceFilterAlgData* psData = (SpaceFilterAlgData*)pData;

    double fResult = 0.0;
    for(int i=0; i<9; i++)
        fResult += (psData->dTemplate[i]*afWin[i]);
        
    return (float)fResult;
}

上面的算法函数也搞定了，下面就是开始调用的问题了。在调用之前，我们还需要写一个函数，就是创建算法的参数结构体的函数，废话不多说，上代码：

void* CreateSpaceFilterAlgData(double* adfTemplate)
{
    SpaceFilterAlgData* pData =
        (SpaceFilterAlgData*)CPLMalloc(sizeof(SpaceFilterAlgData));

    if(adfTemplate != NULL)
        memcpy(pData->dTemplate, adfTemplate, sizeof(double)*9);
    else
        memset(pData->dTemplate, 0, sizeof(double)*9);

    return pData;
}

有了上面的准备工作，那么下面就开始写一个空间滤波的函数了：

/**
* @brief 图像空间滤波，仅适用于3*3模板
* @param pszSrcFile        输入文件路径
* @param pszDstFile        输出文件路径
* @param pTemplate        3*3模板算子
* @param eDataType        输出数据类型，默认为GDT_Unknown，与原图像类型一致
* @param pszFormat        输出文件格式，详细参考GDAL支持数据类型
* @param pProcess        进度条指针
* @return 返回值，表示计算过程中出现的各种错误信息
*/
int ImageSpaceFilter(const char* pszSrcFile, const char* pszDstFile, double *pTemplate, 
                     GDALDataType eDataType, const char *pszFormat, CProcessBase* pProcess)
{
    if(pProcess != NULL)
    {
        pProcess->ReSetProcess();
        pProcess->SetMessage("正在执行空间滤波...");
    }

    if( pszSrcFile == NULL )
    {
        if(pProcess != NULL)
            pProcess->SetMessage("源文件路径为空！");

        return RE_FILENOTEXIST;
    }

    if( pszDstFile == NULL )
    {
        if(pProcess != NULL)
            pProcess->SetMessage("输出文件路径为空！");

        return RE_FILENOTEXIST;
    }

    GDALAllRegister();

    GDALDatasetH hSrcDataset = NULL;
    GDALDatasetH hDstDataset = NULL;
    GDALRasterBandH hSrcBand = NULL;
    GDALRasterBandH hDstBand = NULL;
    GDALDriverH hDriver = NULL;

    char **papszCreateOptions = NULL;
    GDALProgressFunc pfnProgress = ALGTermProgress;

    //打开图像
    hSrcDataset = GDALOpen( pszSrcFile, GA_ReadOnly );
    if( hSrcDataset == NULL )
    {
        if(pProcess != NULL)
            pProcess->SetMessage("输入文件不能打开！");

        return RE_FILENOTEXIST;
    }

    int nXSize = GDALGetRasterXSize(hSrcDataset);
    int nYSize = GDALGetRasterYSize(hSrcDataset);
    int nDstBands = GDALGetRasterCount(hSrcDataset);

    double  adfGeoTransform[6];
    GDALGetGeoTransform(hSrcDataset, adfGeoTransform);

    hDriver = GDALGetDriverByName(pszFormat);
    if( hDriver == NULL)
    {
        if(pProcess != NULL)
            pProcess->SetMessage("不能创建制定类型的文件，请检查该文件类型GDAL是否支持创建！");

        GDALClose( hSrcDataset );

        return RE_FILENOTSUPPORT;
    }

    hSrcBand = GDALGetRasterBand( hSrcDataset, 1 );
    GDALDataType eDstDataType = GDALGetRasterDataType(hSrcBand);

    if(eDataType != GDT_Unknown)
        eDstDataType = eDataType;

    hDstDataset = GDALCreate(hDriver,
        pszDstFile,
        nXSize,
        nYSize,
        nDstBands,
        eDstDataType,
        papszCreateOptions);

    if( hDstDataset == NULL )
    {
        if(pProcess != NULL)
            pProcess->SetMessage("创建输出文件失败，请检查文件路径是否正确！");

        GDALClose( hSrcDataset );
        return RE_CREATEFAILED;    //创建图像失败
    }

    GDALSetGeoTransform(hDstDataset, adfGeoTransform);
    GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset));

    double dfDstNoDataValue = 0;
    int bDstHasNoData = FALSE;
    void* pData = NULL;
    GDALGeneric3x3ProcessingAlg pfnAlg = NULL;

    pData = CreateSpaceFilterAlgData(pTemplate);
    pfnAlg = SpaceFilterAlg;

    for(int i=1; i<=nDstBands; i++)
    {
        hSrcBand = GDALGetRasterBand( hSrcDataset, i );
        hDstBand = GDALGetRasterBand( hDstDataset, i );

        if (bDstHasNoData)
            GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue);

        if (GDALGeneric3x3Processing(hSrcBand, hDstBand,
            pfnAlg, pData, pfnProgress, pProcess) != CE_None)
        {
            if(pProcess!=NULL)
            {
                if(!pProcess->m_bIsContinue)
                {
                    GDALClose(hSrcDataset);
                    GDALClose(hDstDataset);
                    CPLFree(pData);

                    CSLDestroy( papszCreateOptions );

                    return RE_USERCANCEL;
                }
            }

            return RE_FAILED;
        }
    }

    GDALClose(hSrcDataset);
    GDALClose(hDstDataset);
    CPLFree(pData);

    CSLDestroy( papszCreateOptions );

    if(pProcess != NULL)
        pProcess->SetMessage("计算完成！");

    return RE_SUCCESS;
}

关于这些函数的说明，代码里面都有注释的，基本全部依赖GDAL，没有其他的库，还有那些进度条啊，返回值啥的，参考我之前的博客里面的东西。到这里的话，我们就完成了空间滤波的所有函数，下面就差调用了。好吧，接上面，我们用那个拉普拉斯的算子。调用的代码如下：

int LaplacianFilter()
{
    CConsoleProcess *pProgress = new CConsoleProcess();  
    double dTemplate[9] = 
    {
        0, 1, 0,
        1, -4, 1,
        0, 1, 0
    };
    
    int iRev = ImageSpaceFilter("C:\\Test.img", "C:\\Lapiacian.tif", dTemplate, GDT_Unknown, "GTiff", pProgress);
    if(iRev == RE_SUCCESS)
        printf("Success!\n");
    else
        printf("Failed!\n");
        
    delete pProgress;
}

程序执行的原图和结果图，如图2和图3所示：

图2 原始图像

图3 Laplacian滤波后的图像

        明天接着这篇，写一个使用这个3×3的模板做DEM提取坡度坡向的示例。