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

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
        depth = 4;
	size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum;
	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