使用GDAL實現DEM的地貌暈渲圖(三)

charlee44發表於2019-07-27

1. 原理

之前在《使用GDAL實現DEM的地貌暈渲圖(一)》《使用GDAL實現DEM的地貌暈渲圖(二)》這兩篇文章中詳細介紹了DEM生成地貌暈渲圖的原理與實現。不過之前生成的都是暈渲強度值對應的灰度圖,而實際的應用過程中都會將DEM暈渲成彩色圖。

1) ArcMap生成彩色暈渲圖

可以通過ArcMap的做法來參考如何生成彩色暈渲圖(參考[1]),在ArcMap中生成彩色暈渲圖的步驟如下:

  1. 通過山體陰影工具生成灰度暈渲圖,這一點與前面文章介紹的相一致。
  2. 然後在原DEM圖的顯示中,選擇最大最小拉伸顯示,然後選擇一個合適的彩色色帶賦值。
  3. 最後,將步驟一的灰度暈渲圖設定一定的透明度,疊加到步驟二的彩色圖上,就生成了最終具有立體感的彩色暈渲圖。

ArcMap生成的彩色暈渲圖:
AraMap生成的彩色暈渲圖

2) 彩色色帶賦值

不難發現,生成彩色暈渲圖的關鍵是第二步:要選取合適的色帶,讓色帶根據對應的高程賦值。查閱了不少的資料,這個色帶應該沒有固定合適通用的模板,是需要自己根據具體的需要調整的。比如,海平面可以賦值成藍色;高山山頂賦值成白色;戈壁沙漠賦值成黃色;草原森林賦值成綠色,這些地貌特徵都具有一定的高程相關性,可以根據一定的絕對高程區間賦值。

我這裡採取的做法還是跟ArcMap一致,選取漸變色帶來賦值,將漸變色帶約束到DEM的最小最大高程。考慮到地貌的多變性,我這裡生成了藍-綠-黃-紅-紫的多段的漸變色帶。這樣DEM的暈渲效果就是越低越藍,越高越紫。

一般為了保證過渡效果會選擇漸變色帶,漸變色帶的生成也比較簡單,選擇頭尾兩個的顏色的RGB值和一定的漸變範圍,分別讓RGB的值勻速變換就行了。彩色色帶的生成演算法如下(可參考第二部分的具體實現來理解):

//顏色查詢表
vector<F_RGB> tableRGB(256);        

//生成漸變色
void Gradient(F_RGB &start, F_RGB &end, vector<F_RGB> &RGBList)
{
    float dr = (end.R - start.R) / RGBList.size();
    float dg = (end.G - start.G) / RGBList.size();
    float db = (end.B - start.B) / RGBList.size();
    for (size_t i = 0; i < RGBList.size(); i++)
    {
        RGBList[i].R = start.R + dr * i;
        RGBList[i].G = start.G + dg * i;
        RGBList[i].B = start.B + db * i;
    }
}

//初始化顏色查詢表
void InitColorTable()
{
    F_RGB blue(17, 60, 235);//藍色    
    F_RGB green(17, 235, 86);//綠色
    vector<F_RGB> RGBList(60);
    Gradient(blue, green, RGBList);
    for (int i = 0; i < 60; i++)
    {
        tableRGB[i] = RGBList[i];
    }
    
    F_RGB yellow(235, 173, 17);//黃色 
    RGBList.clear();
    RGBList.resize(60);
    Gradient(green, yellow, RGBList);
    for (int i = 0; i < 60; i++)
    {
        tableRGB[i+60] = RGBList[i];
    }
    
    F_RGB red(235, 60, 17);//紅色
    RGBList.clear();
    RGBList.resize(60);
    Gradient(yellow, red, RGBList);
    for (int i = 0; i < 60; i++)
    {
        tableRGB[i + 120] = RGBList[i];
    }
        
    F_RGB white(235, 17, 235);//紫色
    RGBList.clear();
    RGBList.resize(76);
    Gradient(red, white, RGBList);
    for (int i = 0; i < 76; i++)
    {
        tableRGB[i + 180] = RGBList[i];
    }
}

3) 顏色疊加

第一步和第二步分別生成了暈渲強度圖和高程彩色色帶圖,第三步就是將兩者的顏色疊加,生成最終的效果圖。其實顏色疊加的原理特別簡單,對於暈渲強度圖的畫素值A,令其不透明度為α;對應的高程彩色色帶圖的畫素值B,那麼最後疊加的畫素值 C=αA+(1-α)B。可以這麼理解:光線到達A,由於透光性只產生了αA的效果,還有(1-α)強度的光線射到B,又產生了(1-α)B的效果,兩者疊加就是αA+(1-α)B。

2. 實現

繼續改造之前的程式碼,最終的實現過程如下:

#include <iostream>
#include <algorithm>
#include <gdal_priv.h>
#include <osg/Vec3d>
#include <fstream>

using namespace std;
using namespace osg;

//RGB顏色
struct F_RGB
{
    float R;
    float G;
    float B;

    F_RGB()
    {
    }

    F_RGB(float r, float g, float b)
    {
        R = r;
        G = g;
        B = b;
    }

    F_RGB(const F_RGB& rgb) 
    {
        R = rgb.R;
        G = rgb.G;
        B = rgb.B;
    }

    F_RGB& operator= (const F_RGB& rgb)
    {   
        R = rgb.R;
        G = rgb.G;
        B = rgb.B;
        return *this;
    }

};

//顏色查詢表
vector<F_RGB> tableRGB(256);        

//生成漸變色
void Gradient(F_RGB &start, F_RGB &end, vector<F_RGB> &RGBList)
{
    float dr = (end.R - start.R) / RGBList.size();
    float dg = (end.G - start.G) / RGBList.size();
    float db = (end.B - start.B) / RGBList.size();
    for (size_t i = 0; i < RGBList.size(); i++)
    {
        RGBList[i].R = start.R + dr * i;
        RGBList[i].G = start.G + dg * i;
        RGBList[i].B = start.B + db * i;
    }
}

//初始化顏色查詢表
void InitColorTable()
{
    F_RGB blue(17, 60, 235);//藍色    
    F_RGB green(17, 235, 86);//綠色
    vector<F_RGB> RGBList(60);
    Gradient(blue, green, RGBList);
    for (int i = 0; i < 60; i++)
    {
        tableRGB[i] = RGBList[i];
    }
    
    F_RGB yellow(235, 173, 17);//黃色 
    RGBList.clear();
    RGBList.resize(60);
    Gradient(green, yellow, RGBList);
    for (int i = 0; i < 60; i++)
    {
        tableRGB[i+60] = RGBList[i];
    }
    
    F_RGB red(235, 60, 17);//紅色
    RGBList.clear();
    RGBList.resize(60);
    Gradient(yellow, red, RGBList);
    for (int i = 0; i < 60; i++)
    {
        tableRGB[i + 120] = RGBList[i];
    }
        
    F_RGB white(235, 17, 235);//紫色
    RGBList.clear();
    RGBList.resize(76);
    Gradient(red, white, RGBList);
    for (int i = 0; i < 76; i++)
    {
        tableRGB[i + 180] = RGBList[i];
    }
}

//根據高程選顏色
inline int GetColorIndex(float z, float min_z, float max_z)
{
    int temp = floor((z - min_z) * 255 / (max_z - min_z) + 0.6);
    return temp;
}

// a b c
// d e f
// g h i
double CalHillshade(float *tmpBuf, double Zenith_rad, double Azimuth_rad, double dx, double dy, double z_factor)
{
    double dzdx = ((tmpBuf[2] + 2 * tmpBuf[5] + tmpBuf[8]) - (tmpBuf[0] + 2 * tmpBuf[3] + tmpBuf[6])) / (8 * dx);
    double dzdy = ((tmpBuf[6] + 2 * tmpBuf[7] + tmpBuf[8]) - (tmpBuf[0] + 2 * tmpBuf[1] + tmpBuf[2])) / (8 * dy);

    double Slope_rad = atan(z_factor * sqrt(dzdx*dzdx + dzdy*dzdy));
    double Aspect_rad = 0;
    if (abs(dzdx) > 1e-9)
    {
        Aspect_rad = atan2(dzdy, -dzdx);
        if (Aspect_rad < 0)
        {
            Aspect_rad = 2 * PI + Aspect_rad;
        }
    }
    else
    {
        if (dzdy > 0)
        {
            Aspect_rad = PI / 2;
        }
        else if (dzdy < 0)
        {
            Aspect_rad = 2 * PI - PI / 2;
        }
        else
        {
            Aspect_rad = Aspect_rad;
        }
    }

    double Hillshade = 255.0 * ((cos(Zenith_rad) * cos(Slope_rad)) + (sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad - Aspect_rad)));
    return Hillshade;
}


int main()
{
    GDALAllRegister();          //GDAL所有操作都需要先註冊格式
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支援中文路徑
    CPLSetConfigOption("GDAL_DATA", "D:/Work/GDALBuild/gdal_build_result/data");  //支援中文路徑

    //const char* demPath = "D:/CloudSpace/我的技術文章/素材/DEM的渲染/dst.tif";
    const char* demPath = "D:/2.tif";
        
    GDALDataset* img = (GDALDataset *)GDALOpen(demPath, GA_ReadOnly);
    if (!img)
    {
        cout << "Can't Open Image!" << endl;
        return 1;
    }

    int imgWidth = img->GetRasterXSize();   //影象寬度
    int imgHeight = img->GetRasterYSize();  //影象高度
    int bandNum = img->GetRasterCount();    //波段數
    int depth = GDALGetDataTypeSize(img->GetRasterBand(1)->GetRasterDataType()) / 8;    //影象深度

    GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName("GTIFF"); //影象驅動
    char** ppszOptions = NULL;
    ppszOptions = CSLSetNameValue(ppszOptions, "BIGTIFF", "IF_NEEDED"); //配置影象資訊
    const char* dstPath = "D:\\dst.tif";
    int bufWidth = imgWidth;
    int bufHeight = imgHeight;
    int dstBand = 3;
    int dstDepth = 1;
    GDALDataset* dst = pDriver->Create(dstPath, bufWidth, bufHeight, dstBand, GDT_Byte, ppszOptions);
    if (!dst)
    {
        printf("Can't Write Image!");
        return false;
    }

    dst->SetProjection(img->GetProjectionRef());
    double padfTransform[6] = { 0 };
    if (CE_None == img->GetGeoTransform(padfTransform))
    {
        dst->SetGeoTransform(padfTransform);
    }

    //申請buf
    size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum * depth;
    float *imgBuf = new float[imgBufNum];
    //讀取
    img->RasterIO(GF_Read, 0, 0, bufWidth, bufHeight, imgBuf, bufWidth, bufHeight,
        GDT_Float32, bandNum, nullptr, bandNum*depth, bufWidth*bandNum*depth, depth);

    if (bandNum != 1)
    {
        return 1;
    }

    //
    double startX = padfTransform[0];           //左上角點座標X
    double dx = padfTransform[1];           //X方向的解析度
    double startY = padfTransform[3];           //左上角點座標Y
    double dy = padfTransform[5];           //Y方向的解析度

    //
    double minZ = DBL_MAX;
    double maxZ = -DBL_MAX;
    double noValue = img->GetRasterBand(1)->GetNoDataValue();
    
    //
    InitColorTable();
    for (int yi = 0; yi < bufHeight; yi++)
    {
        for (int xi = 0; xi < bufWidth; xi++)
        {
            size_t m = (size_t)bufWidth * yi + xi;
            double x = startX + xi * dx;
            double y = startY + yi * dy;
            double z = imgBuf[m];
        
            if (abs(z - noValue) < 0.01 || z<-11034 || z>8844.43)
            {
                continue;
            }

            minZ = (std::min)(minZ, z);
            maxZ = (std::max)(maxZ, z);
        }
    }
        
    //申請buf
    size_t dstBufNum = (size_t)bufWidth * bufHeight * dstBand;
    GByte *dstBuf = new GByte[dstBufNum];
    memset(dstBuf, 0, dstBufNum*sizeof(GByte));

    //設定方向:平行光
    double solarAltitude = 45.0;
    double solarAzimuth = 315.0;
    
    //
    double Zenith_rad = osg::DegreesToRadians(90 - solarAltitude);
    double Azimuth_math = 360.0 - solarAzimuth + 90;
    if (Azimuth_math >= 360.0)
    {
        Azimuth_math = Azimuth_math - 360.0;
    }   
    double Azimuth_rad = osg::DegreesToRadians(Azimuth_math);
        
    //a b c
    //d e f
    //g h i
    double z_factor = 2;
    double alpha = 0.3;     //A不透明度 α*A+(1-α)*B

    //
    for (int yi = 1; yi < bufHeight-1; yi++)
    {
        for (int xi = 1; xi < bufWidth-1; xi++)
        {           
            size_t e = (size_t)bufWidth * yi + xi;
            size_t f = e + 1;
            size_t d = e - 1;

            size_t b = e - bufWidth;
            size_t c = b + 1;
            size_t a = b - 1;

            size_t h = e + bufWidth;
            size_t i = h + 1;
            size_t g = h - 1;
            
            float tmpBuf[9] = { imgBuf[a], imgBuf[b], imgBuf[c], imgBuf[d], imgBuf[e], imgBuf[f], imgBuf[g],imgBuf[h], imgBuf[i] };
            double Hillshade = CalHillshade(tmpBuf, Zenith_rad, Azimuth_rad, dx, -dy, z_factor);
            GByte value = (GByte)(std::min(std::max(Hillshade, 0.0), 255.0));

            int index = GetColorIndex(imgBuf[e], minZ, maxZ);
            GByte rgb[3] = { (GByte)tableRGB[index].R, (GByte)tableRGB[index].G, (GByte)tableRGB[index].B };

            for (int ib = 0; ib < dstBand; ib++)
            {
                size_t n = (size_t)bufWidth * dstBand * yi + dstBand * xi + ib;
                double v = value * alpha + (1 - alpha) * rgb[ib];
                dstBuf[n] = (GByte)(std::min)((std::max)(v, 0.0), 255.0);
                //dstBuf[n] = (GByte)value;
            }   
        }
    }


    //寫入
    dst->RasterIO(GF_Write, 0, 0, bufWidth, bufHeight, dstBuf, bufWidth, bufHeight,
        GDT_Byte, dstBand, nullptr, dstBand*dstDepth, bufWidth*dstBand*dstDepth, dstDepth);
    
    //釋放
    delete[] imgBuf;
    imgBuf = nullptr;

    //釋放
    delete[] dstBuf;
    dstBuf = nullptr;

    //
    GDALClose(dst);
    dst = nullptr;

    GDALClose(img);
    img = nullptr;

    return 0;
}

最終實現的效果圖如下,可以看到確實實現了高程越低越藍,越高越紫的暈渲效果,同時具有深度感,能看清山脊地貌,效果與ArcMap基本一致,只是配色效果不同。
DEM的地貌暈渲圖

3. 結語

關於DEM的地貌暈渲圖的實現暫時告一段落了。應該來說還是有些模糊不清的地方,在查閱資料的時候就有所感覺,現在關於GIS的基礎原理資料總是不太清晰,原理概念一堆,但就是理解不了。但如果有新的問題或者發現,希望看到這幾篇文章的朋友能批評指正下。

4. 參考

[1]. ArcGIS製圖手冊(3-2)山體陰影和暈渲
[2]. RGB顏色插值漸變原理及演算法
[3]. 兩個RGBA四通道顏色的疊加計算方法與程式碼實現

相關文章