簡述

最近在GDAL的代碼中看見了gdalpansharpen.cpp，於是簡單的嘗試了一下。

融合后的效果比較差，這應該是我對這個算法的使用還不熟悉，還有些地方沒有弄清楚。這個代碼比較新，是2.1版本才加上的，我在看的時候，代碼還有一些小問題，已經在github上提交了issuse了。

融合使用的數據是我在網上找到的高分一號的一景數據，先做了校正，形成全色波段TIFF(單波段)和多光譜波段TIFF(4波段)。

相關知識參考：

C++代碼

代碼基於當前https://github.com/OSGeo/gdal倉庫的master分支構建。

// g++ gdal_pansharpen.cpp -o gdal_pansharpen -I../include -L../lib -lgdal

#include <gdalpansharpen.h>
#include <cpl_auto_close.h>
#include <cpl_error.h>

int main()
{
    GDALAllRegister();
    // 打開全色波段(高分辨率)文件
    GDALDatasetH hPanDset = GDALOpen("/mnt/data/GF1_PMS2_E116.5_N39.6_20130501_L1A0000127213-PAN2_rpc.tiff",GA_ReadOnly);
    CPL_AUTO_CLOSE_WARP(hPanDset,GDALClose);
    VALIDATE_POINTER1(hPanDset,"Open Pansharpen file failed",1);
    // 打開多光譜(低分辨率)文件
    GDALDatasetH hMssDset = GDALOpen("/mnt/data/GF1_PMS2_E116.5_N39.6_20130501_L1A0000127213-MSS2_rpc.tiff",GA_ReadOnly);
    CPL_AUTO_CLOSE_WARP(hMssDset,GDALClose);
    VALIDATE_POINTER1(hPanDset,"Open Spectral Band file failed",1);
    int nSpectralBands = GDALGetRasterCount(hMssDset);

    GDALPansharpenOptions opts;
    opts.ePansharpenAlg = GDAL_PSH_WEIGHTED_BROVEY;     // 超分辨率貝葉斯法,當前僅支持brovery加權
    opts.eResampleAlg    = GRIORA_Cubic;                 // 重采樣至全色光譜波段分辨率的算法
    opts.nBitDepth     = 0;                            // 多光譜波段位深度，0表示默認
    opts.nWeightCount    = nSpectralBands;               // 權重系數數組元素個數(與輸入多光譜波段數一致)
    double* pWeightCount= static_cast<double*>(
        CPLMalloc(opts.nWeightCount * sizeof(double))); // 權重系數數組
    CPL_AUTO_CLOSE_WARP(pWeightCount,CPLFree);
    opts.padfWeights = pWeightCount;
    opts.padfWeights[0] = 0.334;    // 藍
    opts.padfWeights[1] = 0.333;    // 綠
    opts.padfWeights[2] = 0.333;    // 紅
    opts.padfWeights[3] = 0.0;        // 近紅外

    // 設置全色波段(高分辨率)
    opts.hPanchroBand = GDALGetRasterBand(hPanDset,1);
    // 設置多光譜波段
    opts.nInputSpectralBands = nSpectralBands;
    GDALRasterBandH* pInputSpectralBands = static_cast<GDALRasterBandH*>(
        CPLMalloc(sizeof(GDALRasterBandH) * nSpectralBands));
    CPL_AUTO_CLOSE_WARP(pInputSpectralBands,CPLFree);
    opts.pahInputSpectralBands = pInputSpectralBands;
    // std::generatr
    for(int i=0;i< nSpectralBands;++i) {
        opts.pahInputSpectralBands[i] = GDALGetRasterBand(hMssDset,i+1);
    }
    // 設置需要輸出到全色波段分辨率的波段
    opts.nOutPansharpenedBands = 4;
    // 這個數組里面存的是pahInputSpectralBands里波段的索引值
    int panOutPansharpenedBands[4] = { 2, 1, 0, 3}; // 紅、綠、藍、近紅外
    opts.panOutPansharpenedBands = panOutPansharpenedBands;

    opts.bHasNoData = FALSE;   // 全色和多光譜波段是否具有無效值(NoData值)
    opts.dfNoData   = 0.0;     // 全色和多光譜波段的無效值,也將作為輸出的NoData值
    opts.nThreads   = -1;      // 使用的線程數，-1表示使用CPU線程數
    // 設置多光譜波段與全色波段在像素上的移位(保證地理空間位置對齊)
    // 都是相對於全色波段的0,0像素的像素(全色波段像素大小)偏移
    // 也就是兩者的0,0像素的地理空間上的偏移量在全色波段分辨率相應的像素數
    double adfGTPan[6];
    GDALGetGeoTransform(hPanDset,adfGTPan);
    double adfGTMss[6];
    GDALGetGeoTransform(hPanDset,adfGTMss);
    opts.dfMSShiftX = (adfGTPan[0] - adfGTMss[0]) / adfGTPan[1];
    opts.dfMSShiftY = (adfGTPan[3] - adfGTMss[3]) / adfGTPan[5];

    GDALPansharpenOperation operation;
    CPLErr err = operation.Initialize(&opts);
    if(err != CE_None) {
        return -2;
    }
    // 創建輸出文件(和全色波段一樣大)
    int nOutXSize, nOutYSize;
    nOutXSize = GDALGetRasterBandXSize(opts.hPanchroBand);
    nOutYSize = GDALGetRasterBandYSize(opts.hPanchroBand);
    GDALDataType eBufDataType = GDALGetRasterDataType(opts.pahInputSpectralBands[0]);
    eBufDataType = GDT_Float64;
    GDALDriverH hDriver = GDALGetDriverByName("GTiff");
    CPLStringList CreateOption;
    CreateOption.AddNameValue("TILED", "YES");
    CreateOption.AddNameValue("BIGTIFF", "YES");
    CreateOption.AddNameValue("INTERLEAVE", "BAND");
    CreateOption.AddNameValue("COMPRESS", "LZW");   // 中度壓縮
    CreateOption.AddNameValue("ZLEVEL", "6");

    GDALDatasetH hOutDset = GDALCreate(hDriver,
                "/mnt/data/GF1_PMS2_E116.5_N39.6_20130501_L1A0000127213.tif",
                nOutXSize, nOutYSize, nSpectralBands, GDT_UInt16,
                CreateOption);
    CPL_AUTO_CLOSE_WARP(hOutDset,GDALClose);
    VALIDATE_POINTER1(hOutDset,"Create Output file error", -3);

    GDALSetGeoTransform(hOutDset, adfGTPan);
    GDALSetProjection(hOutDset,GDALGetProjectionRef(hPanDset));

     void* pData = CPLMalloc(256*256*GDALGetDataTypeSizeBytes(eBufDataType)*nSpectralBands);
     CPL_AUTO_CLOSE_WARP(pData,CPLFree);

    for(int nYOff = 0; nYOff < nOutYSize; nYOff += 256) {
        for(int nXOff = 0; nXOff < nOutXSize; nXOff += 256) {
            int nXSize = std::min(nOutXSize - nXOff,256);
            int nYSize = std::min(nOutYSize - nYOff,256);
            printf("Process[%6d,%6d,%6d,%6d]\r",nXOff,nYOff,nXOff+nXSize,nYOff+nYSize);

            err = operation.ProcessRegion(nXOff,nYOff,nXSize,nYSize,pData,eBufDataType);
            if(err == CE_Failure) {
                CPLError(err,CPLE_AppDefined,"operation.ProcessRegion");
                return -4;
            }
            int panBandMap[4] = { 1, 2, 3, 4};
            err = GDALDatasetRasterIO(hOutDset, GF_Write,
                        nXOff,nYOff,nXSize,nYSize,
                        pData,nXSize,nYSize,
                        eBufDataType,
                        4,panBandMap,
                        0,0,nXSize*nYSize*GDALGetDataTypeSizeBytes(eBufDataType));
        }
    }
    puts("\nPansharpen finished");

    return 0;
}