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

目錄

• 1. 原理
  • 1) 點法向量
  • 2) 日照方向
    • (1) 太陽高度角和太陽方位角
    • (2) 計算過程
  • 3) 暈渲強度
• 2. 實現
• 3. 參考

1. 原理

之前一直覺得對DEM的渲染就是簡單的根據DEM的高度不一樣賦予不一樣的顏色就能夠。後來實際這麼作的時候獲取的效果跟別的軟件相比，根本體現不出地形起伏的變化。若是要體現出地形的起伏變化，須要獲得地貌暈渲圖才行。暈渲法假設地形接受固定於某一位置光源的平行光線，隨坡面與光源方向的夾角不一樣，產生不一樣色調明暗效果。
根據文獻[1][2],能夠經過計算DEM格網點的法向量與日照方向的的夾角，來肯定該格網點的像素值。

1) 點法向量

咱們知道三點成面，面的法向量就是其三個頂點的法向量（三點成面計算其法向量可參看《已知三點求平面法向量》）。可是一個頂點可能會構成多個不一樣的面，那麼這種存在多個面的頂點的法向量怎麼求呢？其實很簡單，只須要把該點對應面的法向量相加就能夠了。能夠不用求平均，由於反正最後是要正規化的。
具體到DEM上來講，能夠將一個DEM的矩形網格分紅兩個一樣順序排列的三角形，每一個點涉及1到6個不等的面法向量。將這些面法向量相加並正則化，就獲得了每一個點的法向量。以下圖所示。

2) 日照方向

關於日照方向，我在《經過OSG實現對模型的日照模擬》這篇文章裏面有過詳細的表述，那麼這裏就直接搬運過來。

(1) 太陽高度角和太陽方位角

對於太陽光照來講，其方向並非隨便設置的。這裏須要引入太陽高度角和太陽方位角兩個概念，經過這兩個角度，能夠肯定日照的方向。
太陽高度角指的就是太陽光的入射方向和地平面之間的夾角；而太陽方位角略微複雜點，指的是太陽光線在地平面上的投影與當地子午線的夾角，可近似地看做是豎立在地面上的直線在陽光下的陰影與正南方向的夾角。其中方位角以正南方向爲0，由南向東向北爲負，有南向西向北爲正。例如太陽在正東方，則其方位角爲-90度；在正東北方時，方位角爲-135度；在正西方時，方位角是90度，在正西北方爲135度；固然在正北方時方位角能夠表示爲正負180度。
DEM渲染中通常將太陽高度角設置成45度，太陽方位角設置成315度，即西北照東南。

(2) 計算過程

根據上述定義，對於空間某一點的日照光線，能夠有以下示意圖。

令太陽光線長度L1=1，有以下推算過程：

α是太陽高度角，則日照方向Z長度L3=sin(α)；
L1在地平面（XY）平面的長度L2 = cos(α)；
β是太陽方位角，則日照方向X長度L5 = L2cos(β);
同時日照方向Y長度L4 = L2sin(β)。

所以，對於太陽高度角α和太陽方位角β，日照光線的單位向量n(x,y,z)爲：

X = cos(α)cos(β);
Y = cos(α)sin(β);
Z = sin(α);

3) 暈渲強度

在文獻[1][2]中提出由格網點法向量與光源方向的夾角，肯定當前格網點的暈渲強度值。其暈渲圖像素值i_cellvalue_hillshade計算公式以下所示（其中d_vectorvalue是法向量，a_rayvector是日照方向向量）：

這裏的夾角d_raytovector_angle的計算公式略微奇怪。其實夾角計算遠那麼複雜，若是法向量和日照方向向量都已經正則化，那麼其夾角能夠直接用向量點積公式：

d_vectorvalue * a_rayvector = |d_vectorvalue| * |a_rayvector|* cos(d_raytovector_angle) = cos(d_raytovector_angle)

即其夾角的餘弦值爲兩個正則化向量的點積。通過驗證，二者計算出來的夾角值是一致的。

2. 實現

根據上述原理，其具體實現以下。我這裏用到了GDAL來讀寫DEM和圖像，此外還有向量計算用到了osg庫裏面的內容，若是沒有osg，能夠本身簡單實現下，都是很簡單的數學知識。

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

using namespace std;
using namespace osg;

//計算三點成面的法向量
void Cal_Normal_3D(const Vec3d& v1, const Vec3d& v2, const Vec3d& v3, Vec3d &vn)
{
    //v1(n1,n2,n3);
    //平面方程: na * (x – n1) + nb * (y – n2) + nc * (z – n3) = 0 ;
    double na = (v2.y() - v1.y())*(v3.z() - v1.z()) - (v2.z() - v1.z())*(v3.y() - v1.y());
    double nb = (v2.z() - v1.z())*(v3.x() - v1.x()) - (v2.x() - v1.x())*(v3.z() - v1.z());
    double nc = (v2.x() - v1.x())*(v3.y() - v1.y()) - (v2.y() - v1.y())*(v3.x() - v1.x());

    //平面法向量
    vn.set(na, nb, nc);
}

int main()
{
    GDALAllRegister();          //GDAL全部操做都須要先註冊格式
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路徑

    const char* demPath = "D:/CloudSpace/個人技術文章/素材/DEM的渲染/dst.tif";
    //const char* demPath = "D:/Data/imgDemo/K51E001022/k51e001022dem/w001001.adf";
    
    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 = 1;
    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();

    vector<Vec3d> dotList;          //全部的頂點 
    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];
            dotList.push_back(Vec3d(x, y, z));

            if (abs(z - noValue) < 0.01 || z<-11034 || z>8844.43)
            {
                continue;
            }

            minZ = (std::min)(minZ, z);
            maxZ = (std::max)(maxZ, z);
        }
    }

    //計算每一個面的法向量
    multimap<size_t, size_t> dot_face;
    vector<Vec3d> faceNomalList;
    for (int yi = 0; yi < bufHeight - 1; yi++)
    {
        for (int xi = 0; xi < bufWidth - 1; xi++)
        {
            size_t y0x0 = (size_t)bufWidth * yi + xi;
            size_t y1x0 = (size_t)bufWidth *(yi + 1) + xi;
            size_t y0x1 = (size_t)bufWidth *yi + xi + 1;
            size_t y1x1 = (size_t)bufWidth *(yi + 1) + xi + 1;

            Vec3d vn;
            Cal_Normal_3D(dotList[y0x0], dotList[y1x0], dotList[y0x1], vn);
            dot_face.insert(make_pair(y0x0, faceNomalList.size()));
            dot_face.insert(make_pair(y1x0, faceNomalList.size()));
            dot_face.insert(make_pair(y0x1, faceNomalList.size()));
            faceNomalList.push_back(vn);

            Cal_Normal_3D(dotList[y1x0], dotList[y1x1], dotList[y0x1], vn);
            dot_face.insert(make_pair(y1x0, (int)faceNomalList.size()));
            dot_face.insert(make_pair(y1x1, (int)faceNomalList.size()));
            dot_face.insert(make_pair(y0x1, (int)faceNomalList.size()));
            faceNomalList.push_back(vn);
        }
    }


    //申請buf
    size_t dstBufNum = (size_t)bufWidth * bufHeight * dstBand * dstDepth;
    GByte *dstBuf = new GByte[dstBufNum];
    memset(dstBuf, 255, dstBufNum*sizeof(GByte));

    //設置方向：平行光
    double solarAltitude = 45.0;
    double solarAzimuth = 315.0;
    osg::Vec3d arrayvector(0.0f, 0.0f, -1.0f);
    double fAltitude = osg::DegreesToRadians(solarAltitude);                //光源高度角
    double fAzimuth = osg::DegreesToRadians(solarAzimuth);          //光源方位角
    arrayvector[0] = cos(fAltitude)*cos(fAzimuth);
    arrayvector[1] = cos(fAltitude)*sin(fAzimuth);
    arrayvector[2] = sin(fAltitude);
    
    vector<Vec3d> normalList;
    double alpha = 0.5;     //A不透明度 α*A+(1-α)*B
    for (int yi = 0; yi < bufHeight; yi++)
    {
        for (int xi = 0; xi < bufWidth; xi++)
        {
            size_t m = (size_t)bufWidth * yi + xi;

            auto beg = dot_face.lower_bound(m);
            auto end = dot_face.upper_bound(m);

            Vec3d n(0, 0, 0);
            int ci = 0;
            for (auto it = beg; it != end; ++it)
            {
                n = n + faceNomalList[it->second];
                ci++;
            }

            n.normalize();
            normalList.push_back(n);

            double angle = osg::RadiansToDegrees(acos(n * arrayvector));
            //double d_tmp = (n - arrayvector).length();
            //double angle = osg::RadiansToDegrees(asin(d_tmp / 2.0)) * 2;

            double value = (std::min)((std::max)(angle / 90 * 255, 0.0), 255.0);
            dstBuf[m] = (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裏面生成的暈渲效果比較以下，應該仍是比較接近的：

這裏只是獲得了暈渲的灰白強度圖，後續會繼續實現彩色暈渲圖的實現。

3. 參考

[1].地貌暈渲圖的生成原理與實現.丁宇萍,蔣球偉
[2].DEM-地貌暈渲圖的生成原理