簡述

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

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

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

相關知識參考：
遙感影像的“全色”與“多光譜”
高分一號（GF-1）衛星影像資料介紹
國內遙感衛星資源綜述
高分一號影像處理流程
遙感影像融合方法(英文)
ENVI 遙感影象融合

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.166;    // 藍
    opts.padfWeights[1] = 0.167;    // 綠
    opts.padfWeights[2] = 0.167;    // 紅
    opts.padfWeights[3] = 0.5;        // 近紅外

    // 設定全色波段(高解析度)
    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;
}
 

效果對比

GDAL融合效果和原始多光譜波段對比

GDAL融合效果和原始全色波段對比

ARCGIS融合效果與原始全色和多光譜對比

ArcGIS融合過程使用工具箱-->系統工具箱-->Data Management Tools-->柵格-->柵格處理-->建立全色銳化的柵格資料集。

左邊ArcGIS融合效果，右邊原始多光譜影像

GDAL融合效果與ArcGIS融合效果對比

左邊GDAL融合效果，右邊ArcGIS融合效果