1. 問題
之前我在《使用GDAL實現DEM的地貌暈渲圖(一)》這篇文章裡面講述了DEM暈渲圖的生成原理與實現,大體上來講是通過計算DEM格網點的法向量與日照方向的的夾角,來確定該格網點的暈渲強度值。但其實關於這一點我不是很理解,這樣做隨著坡面與光源方向的夾角不同,確實產生了不同色調明暗效果;但暈渲圖同時又有“陰坡面越陡越暗,陽坡面越陡越亮”的特性的,而陰陽坡面的劃分又是跟坡度和坡向相關,之前的生成方法能體現出這種特性嗎?
經過查閱資料,卻在ArcGIS的幫助文件《山體陰影工具的工作原理》(線上版本可檢視這篇文章《ArcGIS教程:山體陰影工作原理》)中查閱到了暈渲圖的另外一種生成演算法。利用直接利用坡度和坡向的關係,算出每個點的山體陰影值:
並且,在該文件中,還附帶了一個具體的計算示例:
具體的演算法過程說的很清楚了,可惜的就是不太明白具體的原理是什麼,在幫助中指向了一本1998的英文文獻[Burrough, P. A. and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), 190 pp]也實在是沒法深入查閱深究。而在查閱中文論文的時候,關於這一段的描述也是互相抄襲,摘錄如下:
這一段的論述反正我是沒看明白的,也就不多做論述了,希望看懂這個演算法的大神能指點我一下。
2. 實現
雖然更深入的原理沒弄明白,不過作為應用者卻足夠能夠實現其演算法了。我這裡通過GDAL實現了暈渲圖的生成:
#include <iostream>
#include <algorithm>
#include <gdal_priv.h>
#include <osg/Vec3d>
#include <fstream>
using namespace std;
using namespace osg;
// 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"); //支援中文路徑
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方向的解析度
//申請buf
size_t dstBufNum = (size_t)bufWidth * bufHeight * dstBand * dstDepth;
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 = 1;
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);
dstBuf[e] = (GByte)(std::min(std::max(Hillshade, 0.0),255.0));
}
}
//寫入
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]. ArcGIS幫助:山體陰影工具的工作原理。
[2]. 基於視覺表象的彩色暈渲地圖色彩設計.郭禮珍等.2004